Genetic basis of natural variation in body pigmentation in Drosophila melanogaster

Abstract

Body pigmentation in insects and other organisms is typically variable within and between species and is often associated with fitness. Regulatory variants with large effects at bab1, t and e affect variation in abdominal pigmentation in several populations of Drosophila melanogaster. Recently, we performed a genome wide association (GWA) analysis of variation in abdominal pigmentation using the inbred, sequenced lines of the Drosophila Genetic Reference Panel (DGRP). We confirmed the large effects of regulatory variants in bab1, t and e; identified 81 additional candidate genes; and validated 17 candidate genes (out of 28 tested) using RNAi knockdown of gene expression and mutant alleles. However, these analyses are imperfect proxies for the effects of segregating variants. Here, we describe the results of an extreme quantitative trait locus (xQTL) GWA analysis of female body pigmentation in an outbred population derived from light and dark DGRP lines. We replicated the effects on pigmentation of 28 genes implicated by the DGRP GWA study, including bab1, t and e and 7 genes previously validated by RNAi and/or mutant analyses. We also identified many additional loci. The genetic architecture of Drosophila pigmentation is complex, with a few major genes and many other loci with smaller effects.

Keywords:

• Drosophila melanogaster Genetic Reference Panel
• extreme QTL mapping
• genome wide association study
• quantitative genetics
• replication study

Abbreviations

AIP	=	advanced intercross population
DGRP	=	Drosophila Genetic Reference Panel
GWA	=	genome wide association
LD	=	linkage disequilibrium
RNAi	=	RNA interference
xQTL	=	extreme quantitative trait locus

Introduction

Body pigmentation contributes to the spectacular biodiversity present in nature and mediates mate choice, mimicry, and physiological functions such as thermoregulation and UV resistance. Thus, pigmentation is a significant contributor to fitness.¹ In order to understand how complex traits such as pigmentation evolve, we must first identify the genetic variants underlying phenotypic variation.

The amount of pigmentation on the abdomen of Drosophila melanogaster is a relatively simple morphological trait which serves as a model for mapping the genetic basis of variation in complex traits. Previous studies identified a handful of genes with alleles of large effect that contribute to variation in female D. melanogaster abdominal pigmentation. Kopp et al.² found that a quantitative trait locus (QTL) corresponding to the bric-a-brac locus (bab1 and bab2) contributes to ˜60% of the color variation recorded among a panel of recombinant inbred lines. In addition to bab1, studies using flies collected from natural populations have identified natural variants in tan (t)³ and ebony (e).^3,4 Additional studies have confirmed the role of bab1 and bab2,^5,6 defined the function⁷ and cis-regulatory elements of e^8,9 and t,¹⁰ and characterized involvement of the regulatory genes double sex (dsx) and Abdominal-B (Abd-B) in patterning the D. melanogaster abdomen.^11,12

After tergite development and spatial patterning,^13,14 some of the enzymes and color-giving dopamine derivatives must be transported from the underlying epidermal cells to and from the cuticle. The epidermis modulates the production of black melanin by the exo- and endocytosis of Megalin-dependent vesicles containing Yellow, a protein necessary for dark melanization.^7,15,16 Cell adhesion proteins and chitin synthesis and binding proteins such as krotzkopf verkehrt (kkv), knickkopf (knk) and retroactive (rtv), are necessary for proper epidermis and cuticle development.¹⁷

In conjunction with classical mutant and biochemical analyses,^13,18-23 these studies have pieced together our current understanding of melanin (brown and black coloration) and sclerotin (yellow/ tan coloration) biosynthesis in D. melanogaster (Fig. 1).^1,24,25,17 However, this is only a small fraction of the 40+ genomic regions/ genes known to affect pigmentation from the catalog of mutant phenotypes previously reported.^12,17,18 Furthermore, genes outside of the pigment biosynthesis pathway may also contribute to natural variation in body coloration.

Figure 1. Construction of AIP for xQTL mapping. (A) Light and dark DGRP lines used for the diallel crosses. (B) Experimental protocol for maintaining the AIP.

Dembeck et al.²⁶ assessed natural variation in female abdominal pigmentation in 175 sequenced inbred lines of the D. melanogaster Genetic Reference Panel (DGRP), derived from the Raleigh, NC population. The two most posterior abdominal segments, tergites 5 and 6 (T5, T6), were scored for the approximate percentage of the tergite colored with black or dark brown melanin. Consistent with other studies of D. melanogaster body pigmentation, there was significant genetic variation in the proportion of melanization and high broad-sense heritabilities for each tergite.²

The genome-wide association (GWA) studies identified 155 DNA variants in or near (± 1 kb) 84 genes associated with the proportion of melanization on T5 (84), T6 (34), and the difference between T5 and T6 (35).²⁶ Several of the top variants associated with variation in pigmentation were in tan (t), ebony (e), and bric-a-brac1 (bab1). Mutational analyses and targeted RNAi-knockdown showed that 17 out of 28 (61%) novel candidate genes implicated by the genome-wide association study affected abdominal pigmentation. Many of the novel genes affect other well-studied pathways and phenotypes, such as wing and bristle development, providing evidence for widespread pleiotropy. We named 4 of the novel, computationally predicted genes based on their mutant and/or RNAi phenotypes: pinstripe (pns, CG7852), triforce (tfc, CG9134), plush (ph, CG1887), and farmer (frm, CG10625). Variation in the novel candidate genes identified by GWA analysis may serve as targets for adaptive evolution and sexual selection in D. melanogaster.

The GWA findings suggested that naturally occurring genetic variation in pigmentation may affect multiple steps in different pathways involved in tergite development and melanization. However, we reported all GWA results based on a nominal P < 10⁻⁵ since quantile-quantile plots indicated a systematic departure from random expectation below this value. Very few variants had P-values less than a Bonferroni correction for multiple tests. Thus, although the top variants are enriched for true positives, we need to perform validation experiments to test all candidate genes to discern which genes are true positives. Two methods are often used to provide additional support for the GWA results, functional testing such as RNAi or mutant analysis or replication in other populations. Dembeck et al.²⁶ assessed changes in pigmentation using RNAi knockdown and mutant alleles of selected candidate genes. This approach is limited by the availability of these resources. Further, failure to confirm a candidate gene using RNAi and mutations does not necessarily mean that the candidate allele segregating in the DGRP does not affect the trait. Here, we present the results of a complementary study, an extreme quantitative trait locus (xQTL) GWA analysis²⁷ in an outbred advanced intercross population (AIP) derived from dark and light DGRP lines. These data provide additional evidence supporting the results of the DGRP GWA analysis and previous functional validation experiments.²⁶

Methodology

Drosophila stocks and phenotyping

From the 175 DGRP lines scored by Dembeck et al.²⁶ we selected 3 lightly and 3 darkly melanized lines that were scored on a scale from 0 (no dark melanin) to 4 (100% pigmented with dark melanin) for the proportion of dark melanin on T5 and T6. The light lines were DGRP_217 (T5, T6 line means: 0.8, 0.45); DGRP_808 (0.7, 0.45); and DGRP_181 (1.15, s1). The dark lines were DGRP_93 (2.5, 3.85); DGRP_153 (1.45, 3.85); and DGRP_233 (1.8, 4). Each DGRP line was reared at a controlled adult density (CAD) of 10 males and 10 females on cornmeal-molasses-agar medium at 25°C, 75% relative humidity, and a 12-h light-dark cycle. The flies were allowed to lay eggs for 3 d Virgin males and females from each line were then crossed using a full-diallel design (excluding homozygous parental lines) and reared as described above. Ten virgin female F₁ flies were aged to day 4 and phenotyped using the 0 – 4 proportion of melanization scale.

F₁ virgin males and females were then used to establish Generation (G) 0 of an outbred AIP. Thirty males and 30 females (one male and female of each F₁ genotype) were placed in each of 10 bottles (600 flies total). For G1, 30 male and 30 female virgins were collected from each G0 bottle, mixed according to sex, and then randomly distributed among 10 new rearing bottles. This cycle was then repeated until G14. At G10, G11, G13 and G14, ∼500 virgin females were collected and aged for 4 d The flies were phenotyped, and darkest ˜50 females and lightest ∼50 females were individually flash-frozen in a 0.2 mL microcentrifuge tube, and stored at −80°C until DNA extraction.

Extreme QTL mapping

DNA was extracted from pools of the darkest ˜50 females and lightest ˜50 females at G10, G11, G13 and G14. Sequence libraries were prepared from each DNA sample, barcoded, multiplexed, and sequenced on the Illumina HiSeq 2500 platform using 100bp paired-end sequencing. Sequence reads were mapped to the Drosophila melanogaster reference genome (BDGP5) using BWA-MEM 0.7.10,²⁸ locally realigned, masked for PCR duplicates, and recalibrated for base quality using GATK 2.4.9.²⁹ Allele counts were obtained at polymorphic SNP sites where the 6 parental lines segregated according their DGRP genotype calls. Finally, the difference in allele frequencies in the dark and light pools was called on each site using a test statistic Z = d/s with a 2-sided P-value obtained from the standard normal distribution, where d is the allele frequency difference p₁ – p₂, ${\displaystyle s= \sqrt{p_0\left(1-p_0\right)\left(\frac{1}{n_1}+\frac{1}{n_2}+\frac{1}{h_1}+\frac{1}{h_2}\right)}}$ with p₀ = (p₁+p₂)/2; n₁, n₂ the sequ-ence coverages, and h₁, h₂ the numbers of haploid genomes in the pooled samples respectively. To combine results across multiple generations, a combi-ned Z score was obtained ${\displaystyle Z=\frac{C_1{d}_1+C_2{d}_2+C_3{d}_3+C_4{d}_4}{\sqrt{c_1^2{s}_1^2+c_2^2{s}_2^2+c_3^2{s}_3^2+c_4^2{s}_4^2}}}$ and the P-value drawn from standard normal, where the c₁, …, c₄ are the average coverages in the 4 generations. A total of 595,729 single nucleotide polymorphisms (SNPs) that remained polymorphic across all 4 generations were considered.

Results and Discussion

We created an outbred AIP from 3 darkly pigmented and 3 lightly pigmented DGRP lines. We initially crossed the lines in a full diallel design (excluding homozygous parental genotypes) to ensure equal representation of all lines, and maintained the populations at a census population size of 600 flies per generation, for a total of 14 generations (Fig. 1A, Fig. 1B, Fig. S1). At generations 10, 11, 13 and 14 we scored 500 females and selected the 10% darkest and lightest flies to use for xQTL mapping (Table S1). We sequenced the 8 DNA pools, and mapped QTLs affecting pigmentation by assessing the differences in allele frequency between the dark and light pools in each generation (Table S2). A qualitative visual assessment of the Manhattan plots showed that the QTL peaks are generally consistent across all analyses (Fig. S2). Therefore we combined data from all 4 generations and obtained a single P-value for each SNP (Fig. 2, Table S2). At a Bonferroni-corrected P-value of 8.39 × 10⁻⁸, we identified 7,069 variants in or near (± 1 kb) 1,884 genes from the xQTL analysis combined across generations. We mapped many genes and specific SNPs that overlap with the GWA results of Dembeck et al.,²⁶ providing independent validation for the results of this study. In addition, we also identified variants located within or near a new set of candidate genes that are biologically plausible with respect to variation in D. melanogaster pigmentation.

Figure 2. Manhattan plot for the xQTL mapping GWA analysis, combined across all generations. Gene symbols for loci implicated by the mapping are given. The Bonferroni corrected significance threshold (P = 8.39 × 10⁻⁸) is indicated by the horizontal dashed line.

Gene overlap between DGRP and xQTL GWA analyses

We compared the top candidate genes identified in the DGRP GWA analysis of Dembeck et al.²⁶ and the xQTL mapping GWA analysis for all generations combined, and identified 28 overlapping genes. These included genes previously known to affect Drosophila pigmentation (bab1, bab2, e, osa, t, and the 2 genes in the cis-regulatory region of t, Gustatory receptor 8a (Gr8a) and CG15370); and 11 of the 28 novel candidate genes identified and tested in Dembeck et al.²⁶: roughoid (ru), klarsicht (klar), multiple wing hairs (mwh), pns, ph, tfc, Vesicular monoamine transporter (Vmat), Glucose transporter 1 (Glut1), krotzkopf verkehrt (kkv), CG15803 and CG42340. All but 4 (Vmat, CG42340, CG15803, mwh) of the latter genes caused changes to T5 or T6 pigmentation when analyzed with RNAi and/or mutants. The additional genes that overlapped but were not previously tested were CG5888, CG30485, CG9194, CG42663, empty spiracles (ems), CG13330, Protein tyrosine phosphatase 61F (Ptp61F), CG2493, mir-288, RNA-binding protein 6 (Rbp6) and E5.

Although we had not previously tested Ptp61F, we note that it overlaps with ru, which led to large reductions in T5 and T6 pigmentation when knocked down. The variants in Ptp61F associated with pigmentation could well correspond to a cis-regulatory element for ru. We replicated effects on pigmentation in the xQTL mapping GWA analysis for which mutations and RNAi knockdown did not yield significant effects on pigmentation. We speculate that the naturally segregating variant(s) in these genes affect pigmentation and that the effects are not recapitulated by the available mutations or by reducing expression ubiquitously or with a pannier-GAL4 driver.

SNP overlap between DGRP and xQTL GWA analyses and other studies

Three SNPs located in the cis-regulatory region of t were mapped in the xQTL and DGRP GWA analyses and the xQTL study of variation in tergite 7 pigmentation of wild caught D. melanogaster females by Bastide et al.³ One of these SNPS, X_9121177_SNP, as highly significant in the combined generations analysis (P = 1.34 × 10⁻²⁹). Two other SNPs were significant following Bonferroni correction in G10 and G11 (X_9121129_SNP) or in G13 and G14 (X_9121094_SNP). These SNPs are promising causative variants given their replication across the 3 studies.

Four additional SNPs were shared by the DGRP GWAS and the combined generations xQTL analysis. Two SNPs, 3L_1085642_SNP (P = 1.56 × 10⁻¹⁴) and 3L_1084129_SNP (P = 1.08 × 10⁻¹¹), are located within the intron of bab1. These SNPs fall on each side of the dimorphic element described by Williams et al.¹¹, which contains 3 other bab1 variants, 3L_1084980_SNP, 3L_1085137_SNP, and 3L_1085230_SNP mapped by Dembeck et al. and Rogers et al.³⁰ (Fig. 3). 3L_1084129_SNP is also within a dl transcription factor binding site; again highlighting the association of variants within dl binding sites as was seen in the DGRP GWAS. This overlap and the proximity of these associated SNPs to the cis-regulatory dimorphic element suggest that these SNPs may also functionally affect the regulatory region. 3L_1381936_SNP (P = 1.94 × 10⁻¹⁴) is located within the introns of the overlapping genes ru and Ptp61F. 2R_9404826_SNP is located downstream of Vmat (P = 9.42 × 10⁻⁹). These results provide additional support for variants with small effect sizes contributing to natural variation in female D. melanogaster pigmentation.

Figure 3. Schematic showing the location of the 2 SNPs in bric-a-brac1 that were associated with variation in body pigmentation in independent studies.

Even more candidate genes

One feature of the Manhattan plot from the combined xQTL analysis is the bifurcation of the large QTL peak on the X chromosome into 2 very close peaks. Upon inspection of the genes in these regions, it is evident that the first peak results from the association of genetic variants in t or its cis-regulatory regions. The second peak, somewhat to our surprise, is centered upon variants within or near megalin (mgl), a gene known to be involved in wing cuticle pigmentation in D. melanogaster. mgl mediates the endocytic removal of the Yellow protein, which is required to produce black melanin, from the developing cuticle. When mgl expression is reduced, Yellow accumulates, leading to the ectopic formation of abnormally large melanin granules.¹⁵ To our knowledge, this is the first study to associate natural genetic variation in mgl with variation in abdominal pigmentation in D. melanogaster.

Surprisingly, the most significantly associated SNP on the X chromosome in the combined xQTL analysis was not in the cis-regulatory region of t or in mgl. Instead, the SNP (P = 6.4 × 10⁻³⁴) was located in an intergenic region upstream of suppressor of rudimentary (su(r)). su(r) is located downstream from t and it is possible that the SNP may be located in another t cis-regulatory region. However, this intergenic region is a strong candidate as a cis-regulatory element for su(r). Insects use β-alanine for the synthesis of N-β-alanyl dopamine, which contributes to the yellowish-brown color of the cuticle. There are 2 primary pathways used to synthesize β-alanine. It can be formed from aspartate by aspartate decarboxylase (encoded by black) or from uracil by dihydropyrimidine dehydrogenase (Pyd1, encoded by su(r)), dihydropyrimidase (Pyd2), and β-ureidopropionase (Pyd3).³¹ This supports our previous findings that genetic variation at a number of steps in biosynthetic, regulatory, and developmental pathways may contribute to natural variation in abdominal pigmentation.

The xQTL mapping analyses reveals 3 QTL peaks on chromosome 2R that are consistent in G10, G11 and G14 (Fig. S2). The first of these 2R QTL peaks contains a SNP within an intron of laccase2. laccase2 is a multicopper oxidase that is essential for cuticular sclerotinization.^20,32,33 As with mgl, to the best of our knowledge, this is the first study to associate natural genetic variation in laccase2 with variation in female abdominal pigmentation in D. melanogaster. Similarly, on chromosome 2L, there was a SNP (P = 6.30 × 10⁻¹⁰) 67-bp downstream of Dopa decarboxylase (Ddc), which could affect regulation of Ddc expression.

We identified several other notable and plausible candidate genes based on the xQTL mapping GWA analysis. One of these variants was in Tripartite motif containing 9 (Trim9) (P = 5.8 × 10⁻⁸). Trim proteins are found in all metazoans and have many intracellular functions. Trim9 binds to microtubules and regulates synaptic vesicle exocytosis.³⁴ The most strongly associated variant on chromosome 2L was in the intron of skywalker (sky) (P = 1.7 × 10⁻¹²). sky has a role in the negative regulation of synaptic vesicle recycling.³⁵ Trim9 and sky may be involved in the vesicle dynamics along with mgl and Vmat that transport pigment biosynthesis enzymes or dopamine-derivatives to and from the cuticle. Another gene implicated by a significant SNP (P = 1.2 × 10⁻¹¹), CG10211, is located between Trim9 and sky. It is predicted to have peroxidase activity and reported to be lethal when gene expression is reduced using RNAi.^36,37 Peroxidases have been suggested to facilitate sclerotinization of the cuticle by catalyzing the oxidation of NADA and NBAD to ortho-quinones.³³ There are also several variants on 2L in and around CG9336. RNAi knockdown of this gene results in body color differences on the central notum and the surrounding area.³⁶

There were a handful of other candidate genes each containing several associated SNPs from the combined xQTL mapping analyses implicated by smaller QTL peaks on chromosomes 2R, 3L, and 3R. These include CG30015, fruitless (fru), and obstructor-I (obst-I), yellow-g and yellow-g2. While no functional annotation has been reported for CG30015, RNAi resulted in cuticle color differences on the notum.³⁶ fru has a vital role in sex determination and development and its variants may ultimately affect the expression of sex determination genes known to affect female pigmentation such as doublesex (dsx) and Abdominal-B (Abd-B). Abd-B was also implicated to affect pigmentation in the xQTL analysis. The three closely linked genes obst-I, yellow-g and yellow-g2 are all reported to be associated with cuticle development. obst-I is a chitin binding protein expressed primarily in the larval digestive system.³⁷ yellow-g and yellow-g2 are expressed in the female follicular cells during egg development and are likely to contribute to hardening the eggshell.³⁷ Given that these genes are not known be functional in adult cuticle pigmentation, we hypothesize that the SNPs in the region may indicate the presence of a cis-regulatory element for ecdysoneless, which is nearby and known to affect cuticle development.³⁷

Juvenile hormone (JH) has profound effects during insect development. It regulates embryogenesis, development and eclosion, and stimulates sexual maturation in adults.³⁸ Two SNPs in the haemolymph JH binding protein CG15497 were strongly associated (P = 1.1 × 10⁻¹⁸ and P = 1.3 × 10⁻¹⁸) with pigmentation in the combined xQTL mapping analysis. These proteins are responsible for transporting JH from the corpora allata through the haemolymph to specific target tissues.³⁹ As these variants are located over 400,000 bp away from the known pigmentation gene e, we feel confident that their association with pigmentation is not due to LD with e variants.

Genetic architecture of D. melanogaster abdominal pigmentation

The GWA analyses in the DGRP and the AIP derived from extreme DGRP lines have complementary properties that can be leveraged to replicate true positive candidate genes and variants in the case of additive gene action. The DGRP has the favorable properties of minimal relatedness and low local LD, but with ∼200 lines, only variants with large effects can be detected with high statistical confidence. xQTL mapping in an AIP can be done with very large numbers of individuals, such that top variants can achieve significance following a strict Bonferroni correction for multiple tests; but unless many generations of recombination are allowed, local LD is much higher than in the DGRP, decreasing the mapping precision. In addition, association mapping in the DGRP can only be performed using common variants (minor allele frequencies > 0.05),⁴⁰ whereas the effects of private variants in any of the 6 parental lines used to construct the AIP can be assessed. However, this replication approach will only work if gene action of the causal variants is additive; in the presence of epistatic gene action, differences in allele frequency between the DGRP and the AIP will lead to failure of replication.^6-8

Since we were able to replicate many of the GWA results of Dembeck et al.²⁶ at the levels of both genes and SNPs, we infer the genetic architecture of abdominal pigmentation is largely additive. Prior to the study of Dembeck et al.²⁶ and this report, the genetic basis of natural variation in D. melanogaster abdominal pigmentation was thought to be oligogenic, controlled by variants with large effects at a few loci. Our studies indicated that reveal a more complex genetic basis than this, with many additional loci with smaller effects. Replicated associations of genetic variants in the genes that were tested previously²⁶ using RNAi knockdown or mutant analyses but that did not affect pigmentation suggest that these genes may indeed contribute to variation in pigmentation, and that naturally occurring variants can cause subtle phenotypic effects that cannot be reproduced by mutagenesis or RNAi knockdown. Since most of the associated SNPs are in intronic and intergenic regions, we hypothesize that they affect tissue- or stage-specific gene regulatory motifs. These results highlight the known and potential contributions of cis-regulatory evolution, which likely limits negative pleiotropic effects, to phenotypic variation within the DGRP population.

Previously,²⁶ we hypothesized that Efa6, Klp61F, klar, and pns facilitate vesicle transport of enzymes and/or dopamine-derivatives to and from the cuticle. The replication of the association of variants in these genes in addition to those in Vmat, mgl, sky, and Trim9 provide more support for this hypothesis. Together these genes may represent components of the transport mechanism for cuticle sclerotization and melanization in D. melanogaster.

In summary, we confirm our previous results and provide additional evidence that genetic variation at a number of steps in biosynthetic, regulatory, developmental, and transport may pathways contribute to natural variation in abdominal pigmentation. Our results hold open the door for new hypotheses to be tested on the epithelial cell vesicle dynamics, the locations of cis-regulatory elements, and how metabolic and hormonal regulation can lead to variation in pigmentation. Future studies will provide much needed mechanistic inferences about these novel candidate genes and SNPs and their contributions to the evolution of pigmentation.

Disclosure of Potential Conflicts of Interest

No potential conflicts of interest were disclosed.

Supplemental Material

Supplemental data for this article can be accessed on the publisher's website.

Funding

This work was funded by National Institutes of Health grant R01 GM45146 to TFCM.