Variant discovery

We generated a catalog of sequence variants using whole-genome Illumina next-generation sequencing data from 39 inbred lines from the Drosophila Genetic Reference Panel (DGRP) for which gene expression data are available for young adult flies [6]. First, we used PrInSeS-G [19], which uses de novo local assembly to generate a preliminary list of variant calls for each DGRP line; this tool works by first detecting regions with a fluctuation in sequencing coverage, then using a short fragment from the reference sequence as a seed to build a contig by extending it using the reads; the contig extending stops when another short fragment from the reference sequence is encountered. We then re-aligned each resulting genome against the reference genome to present variants in a coherent way among lines and to reduce variant fragmentation. Finally, we developed and used a genotyping algorithm, which uses the combined variant call set for all DGRP lines to improve variant discovery. This is particularly useful for genomes with low read coverage (Figure S1). This strategy allowed us to identify more than 3.6 million unique sequence variants across all 39 DGRP lines including 2.8 million SNPs, 0.6 million indels, and 0.2 million complex sequence variants (Figure 1A, Table 1). To validate our variant calls, we used five distinct approaches. First, we compared our variant catalogue to indel data from the FLYSNPdb database [20]. We validated 713 deletions and 923 insertions on the breakpoint and complete sequence level, thus covering in total 45% of the FLYSNPdb content. Second, we compared our SNP calls with those published by Mackay et al. [7], and found that 94% of calls made in that study match (Table S1). Third, we used whole genome Roche-454 reads available for the same lines to validate our SNP calls, resulting in confirmation of 98.7% of SNPs compared (Table S2). Fourth, we sequenced mRNA from three DGRP lines (see also below) and differentially aligned the reads to their respective transcriptomes and the reference transcriptome in order to validate variants in coding sequence. Thus, we examined, on average, 147,453 SNPs and 12,084 non-SNPs per line that were covered by mRNA sequencing reads, confirming on average 98.9% of all variant types within coding sequence for these lines (Table S3). The false discovery rate (FDR) is 0.5% for SNPs and 1.5% for indels and complex variants. Moreover, the FDR remains stable for indels up to 30 bp in size but increases to 8.9% for larger indels. Note that some variants could not be reliably assigned to the true or false positive categories due to low coverage. Finally, we validated 92 large indels (200 bp to 5.6 kb) in five DGRP lines using PCR (true positives 77.9%, false positives 7.2%, variant size different from predicted 14.9%, Table S4). Moreover, Sanger sequencing of 10 randomly selected large indels validated by PCR revealed a high breakpoint accuracy, whereby breakpoints of 9 out of 10 indels were perfectly reconstructed consistent with results based on query of smaller indels in the FLYSNPdb database. Together, these results indicate that our variant catalogue is of similar quality as variant data produced for other organisms using high-throughput sequencing since an overall false discovery rate of less than 10% and an overall accuracy of ∼95% or higher for all variant types (i.e., SNPs, indels and complex variants) is in line with numbers reported by these other studies [4], [21]–[23].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Overview of variants.

SNPs are shown in black/grey, insertions in red, deletions in blue, and complex variants in orange. (A) Number of base pairs affected by variants discovered per line, with lines ordered by depth of coverage (green dotted line). The line “Berkeley” is the reference line. (B) Number of unique variants by size (note that variants longer than 1,000 bp are grouped in a single x-coordinate). (C) Representation of variant density (0–10 SNPs/kb, 0–5 indels/kb, 0–5 complex variants/kb) across the euchromatic genome (concentric circles) in 50 kb bins. Large variants (>100 bp) mapping against a close homologous sequence (>90% sequence identify) are linked in the center with green lines representing intra-chromosomal- and black lines inter-chromosomal duplications. (D) Number of unique variants by number of lines.

https://doi.org/10.1371/journal.pgen.1003055.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Type classification of variants.

https://doi.org/10.1371/journal.pgen.1003136.t001

Indels and complex variants

Among all lines, we identified a SNP every 43 bp and an indel or complex variant every 144 bp, together contributing a genetic marker every 33 bp on average. These findings illustrate a remarkably high density of molecular polymorphisms in Drosophila consistently greater than in humans [3], [24] and mice [4]. Deletions and complex variants affect in total 4.2% (5.0 Mb) of the reference euchromatic genome while insertions add 2.1 Mb (of which 0.5 Mb are in insertions ≥100 bp not found with similarity of 90% or above elsewhere in the reference genome). Non-SNP variants are thus a substantial source of genomic variation in Drosophila. Indel size ranges from 1 to 6,082 bp (Table 1, Figure 1B). Single base pair indels are abundant (32%) and small indels (2–10 bp) represent 50% of the indel repertoire. 6,419 structural variants (≥100 bp) represent 1% of all indels with a median size of 208 bp and encompass ∼2.5 Mb of sequence in total. 850 structural variants are homologous (≥90% sequence similarity) to another part of the assembled reference genome, representing either “young” variants or polymorphic forms of segmental duplications (i.e., copy number variants) [25]. Most (70%) of those structural variants correspond to a duplicon on the same chromosome, indicating that intra-chromosomal duplication events occur more frequently (Figure 1C). These observations support a recently proposed theoretical model governing the formation of segmental duplications in D. melanogaster, whereby, after a double strand break, a search for an ectopic homologous region which is preferentially located within the same chromosomal region triggers the repair mechanism [26].

The majority (73%) of complex variants are balanced multi-nucleotide substitutions and most (61%) constitute di-nucleotide substitutions. Balanced multi-nucleotide substitutions may arise as single nucleotide substitutions that happen to occur in adjacent positions or as multi-nucleotide mutational events [27]. Using simulation within all 1,000 bp genomic windows in all 39 lines, we find that there are on average 3.3 times more di- and 16.3 times more tri-nucleotide substitutions than the number of single nucleotide substitutions expected to be adjacent by chance, indicating that complex mutational events constribute substantially to balanced multi-nucleotide substitutions.

The allele frequency spectrum of indels and complex variants decays steeply, with 28.5% of all indels and complex variants present in only one line (Figure 1D). Comparison to neutral spectra (Figure S2) revealed an excess of low frequency alleles for deletions, consistent with previous results suggesting that more of these variants are under purifying selection [13]. We also observed a depletion of high frequency variants among indels and complex variants in protein-coding regions, splice junctions, and UTR sequences (Figure S3), again suggesting the action of pervasive purifying selection. In contrast, we found 452 deletions, 873 insertions, and 824 complex variants present among all 39 lines (Table S5). To evaluate whether these variants represent rare alleles in the reference genome, population-specific variants, or artifacts in our variant calling workflow, we focused on all 30 such insertions and deletions that were present in protein-coding regions. We found that 77% of these variants exactly recapitulate an ancestral allele (D. simulans, D. yakuba, or D. ananassae) and 90% when allowing at most one mismatch between the indel and the ancestral allele (Table S6). Therefore, these variants predominantly constitute rare alleles in the reference genome. The photoreceptor gene Rh6 is an interesting example: evolutionary conservation analysis around the two observed 17 bp and 2 bp insertions in its coding sequence revealed that both insertions perfectly match the ancestral allele of seven out of 11 Drosophila species (Figure S4). Moreover, the resulting gene model supports an Rh6 cDNA clone of the OregonR/white strain, thus revealing that the reference genome harbors a rare 19 bp deletion, which truncates the gene.

Indels and SNPs are not uniformly distributed across the genome, with autosomal centromeric regions as well as the X chromosome containing fewer variants compared to other genomic regions (Figure 1C). Several models have been proposed to explain this pattern, ranging from purifying selection to recurrent hitchhiking to demography [7], [28], [29], although we cannot rule out that low coverage or read mapping quality issues in those regions has affected variant discovery. Moreover, we found that genome-wide SNP, indel, and complex variant densities are strongly correlated on a 1 Mb scale (Spearman's ρ_SNP-del = 0.927, ρ_SNP-ins = 0.90, ρ_SNP-complex = 0.946; Figure S5), suggesting that similar higher-order constraints such as local selection intensities, recombination, or mutation rates may act on all types of variants. To test the impact of recombination, we tested SNP, indel, and complex variant densities against recombination rates (1 cM/Mb), and found that recombination is strongly positively correlated with variant densities (Spearman's ρ_SNPs = 0.54, ρ_insertions = 0.55, ρ_deletions = 0.62, ρ_complex = 0.46), which is line with long-held observations regarding the actions of both pervasive positive selection and background selection in Drosophila [30], [31]. We further examined the SNP distribution around detected indels in more detail, as it has been proposed that indels may act as “mutators” of surrounding sequences [32], [33]. Consistent with this hypothesis, we observed that the SNP density is elevated in close proximity to indels independent of indel type, dropping quickly to background levels (Figure 2A, 2B, and 2C). On average 26% (40%) of all SNPs per line are within 40 bp (100 bp) of a non-SNP variant. An increase in SNP density around indels has previously been observed between species [32], [33]; here, we present genome-wide evidence of this effect within a single species. Thus, the accuracy of a traditional SNP-calling method based on simple read alignment may suffer materially since a considerable portion of the reads that would contribute to SNP calling may not be mapped correctly as they also contain indels; this problem is circumvented by the integrative variant calling approach employed here.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Variants in genomic context.

(A) Density of SNPs around variant breakpoints by variant type. The dashed lines show the SNP density at the same loci but in DGRP lines that do not have the variant. (B) and (C) Density of SNPs near indels with minor allele count 2 to 4 (B) and 11 to 19 (C). The dashed lines show the SNP density at the same locus for DGRP lines without the indel. If indels were mutagenic, one would expect enrichment for low allele count SNPs near the high allele count indels; instead, the allele count of the neighboring SNPs closely matches that of the indel. (D) Density of variants (reference bases affected per Mb) in selected genomic regions. (E) Number of indels in coding regions by indel size. Insertions are in red and deletions in blue. Bars representing indel sizes that are a multiple of three are coloured dark red and blue, respectively.

https://doi.org/10.1371/journal.pgen.1003055.g002

Interestingly, we observed a similar SNP density increase around other SNPs (Figure 2A), suggesting that either SNPs are also mutagenic or, more likely, that alternative explanations need to be considered for the increased SNP density around variants. We plotted the allele frequencies of SNPs around variants and found that high allele frequency SNPs cluster around other high allele frequency variants (including other SNPs) and, correspondingly, low frequency variants cluster together (Figure 2B and 2C). These findings provide compelling evidence that variants in general are likely not mutagenic, since otherwise we would expect to observe a greater number of rare, thus mostly recent variants, around high frequency alleles. Alternatively, this phenomenon could occur if some genomic regions are more susceptible to mutation than others [32], [34]. We therefore examined the SNP density in the same DNA regions, but in those of the 39 lines which do not contain the focal variants. We found no material enrichment in SNP density in these regions (Figure 2A, 2B, and 2C), providing little evidence of locally increased mutation rate.

Functional impact of indels and complex variants

Non-coding regions that are conserved between species contain significantly fewer variants than other non-coding regions (Figure 2D, Mann-Whitney U test P<2.2e-16), confirming that these regions experience intra-species purifying selection as well. We obtained similar results when integrating histone modification data from the modENCODE Project [9]. We observe fewer indels affecting regulatory genomic regions - especially promoters (marked by tri-methylation of histone H3 at lysine 4, H3K4me3) - than unmarked non-coding regions, (Mann-Whitney U test P = 3.1e-14) (Figure S6).

While the majority (>90%) of indels and complex variants fall into intronic or intergenic regions, 80% of all protein-coding genes are affected by indels or complex variants (Table 2). However, only 13% of non-SNP variants within coding regions result in exon disruptions or full gene deletions. More than 50% of indels leading to frame shifts are singletons versus 35% for non-frame shifting indels, indicating that purifying selection is acting against these likely deleterious mutational events. This is illustrated by our observation that the length of indels affecting coding regions is highly biased for sizes that are multiples of three, a pattern that is still visible even up to 51 bp (Figure 2E). Consistent with previous inter-specific comparisons [35], we found that genes affected by disruptive indels or complex variants, are significantly enriched for the functional categories of sensory perception of taste and smell, proteolysis, and innate immune response as well as pathways related to food and drug metabolism (Tables S7 and S8).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Potential functional indels and complex variants in genes.

https://doi.org/10.1371/journal.pgen.1003055.t002

Although the density of variants around genes is lower than the overall background, we observed an increase in variant density almost to background levels directly upstream of the transcription start site (TSS), followed by a dramatic drop across the TSS (Figure 3A). This is consistent with strong selective constraint at the TSS but not the region immediately upstream, possibly because this region has high AT content and it is typically depleted of nucleosomes [36], [37]. Indeed, we observed a strong correlation between the local AT content and indel density around the TSS (Figure S7A), strengthening earlier observations that sequence context or chromatin structure may affect the indel rate [38], [39]. We observed a similar, albeit less striking effect at the transcription end site (TES) (Figure 3B and Figure S7B).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. Cis-associations of variants with gene expression.

(A) and (B) Variant density (blue) and significance of allele associations (red), in males around (A) the transcription start site (TSS) and (B) the transcription end site (TES) averaged out over all transcripts in a 10 kb window. The solid lines are cubic smoothing splines, fit to the data. Transcripts on both strands are orientated such that transcription takes place in the positive direction of the x-axis. The inlet in (A) corresponds to a 100 kb window length. (C) cis-eQTLs discovered in males, females, or both sexes (FDR<10%). (D) Breakdown of cis-eQTL-associated genes by sex. (E) Breakdown of cis-eQTL associated genes, discovered in males or females, by type of variant (i.e., SNP and non-SNP).

https://doi.org/10.1371/journal.pgen.1003055.g003

Identification and characterization of cis-eQTLs

We used published whole adult microarray gene expression data [6] to perform association analysis between the expression and variants within 10 kb of each gene. We initially identified 9,789 (9,434) genetically variable transcripts (FDR<0.001) in males (females) after removing probes from the microarray analysis that contain genomic variation. We then used the Kruskal-Wallis test followed by permutation-based multiple test correction to find significantly associated variants, termed cis-expression QTLs (cis-eQTLs), which may point to underlying functional regulatory elements. QTL studies have so far mostly associated bi-allelic variants with gene expression levels [40], [41]; given the availability of a high-resolution, comprehensive catalog of sequence variants, many of which (9.3%) are multi-allelic within the DGRP population, we grouped expression levels according to the corresponding allele as an input to each test. We conducted ∼3.8 million tests for each sex, and restricted our association analysis to variants present in at least three lines. We found 17,501 cis-eQTLs in 2,033 genes (26% of the 7,889 genes tested) at a false discovery rate of ∼10% (Figure 3C, 3D, and 3E; Tables S9 and S10; and results for 1% and 20% FDR thresholds are listed in Table S11), generating to our knowledge the first cis-eQTL map in Drosophila.

Surprisingly, the majority of cis-eQTLs were found to be sex-specific (58% specific to males and 17% to females; Figure 3C and 3D, and an example can be found in Figure 4A) with males having more than two-fold more cis-eQTL-associated genes compared to females. This result is not an artifact of the significance threshold, as a plot of the underlying association P-values in both sexes clearly reveals that the majority of cis-eQTLs are sex-specific (Figure S8). Further, the variance in gene expression among females is comparable to that among males as 56% of genetically variable transcripts without any cis-eQTL-associations display higher expression variance in females, making it unlikely that such bias affects our findings. In addition, only 38% of transcripts with male-specific cis-eQTLs exhibit greater gene expression variance in females than in males, refuting the hypothesis that the more than two-fold greater number of male- compared to female-specific cis-eQTL-associated genes may be due to a greater variability in gene expression in females compared to males. Intriguingly, we found that male-specific cis-eQTL-associated genes are significantly depleted from the X chromosome (P = 3.6e-11, χ² test), whereas female-specific cis-eQTL (P = 0.17) and unbiased (P = 0.02) genes show a more uniform chromosomal distribution, perhaps consistent with the observed depletion of male-biased genes on the X chromosome [42], [43]. Sex-unbiased cis-eQTL-associated genes were found to have more cis-eQTLs than sex-biased ones (Figure S9) and the effect of cis-eQTLs found in both sexes is larger (Mann-Whitney U test, P<2.2e-16; Figure S10). Sex-unbiased cis-eQTL-associated genes have also a greater tendency to be expressed only in somatic tissues (Figures S11 and S12), possibly explaining why the underlying variants affect gene expression in both sexes. Nevertheless, more than 20% of sex-specific cis-eQTL-associated genes are also only expressed in somatic tissues, indicating that sex-biased changes in gene expression are pervasive [44]. Furthermore, we found that 20% of male-specific cis-eQTL-associated genes are exclusively expressed in the testes or the accessory gland in males, whereas the expression of only few (5%) female-specific counterparts are restricted to the ovaries or spermathecae in females. Approximately 60% of male-specific cis-eQTL-associated genes are expressed in female reproductive tissues (whereby 28% having greater expression than in any somatic tissue) (Figure S13), yet these genes do not yield any cis-eQTLs in females. In other words, despite the fact that these genes are expressed both in males and females, genetic mutations only produce cis-eQTLs in males, indicating that these mutations do not significantly affect gene expression in females. Together with the lower number of female-specific cis-eQTL-associated genes in general, this suggests that the underlying regulatory architecture may be more constrained in females as compared to males.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Examples of cis-eQTLs and their associated genes.

DGRP lines in (A) and (B) are grouped by their allele. Male and female expression levels are depicted in blue and dark pink, respectively. (A) Sex-biased cis-eQTL. A SNP (3R:8,875,391) is associated with higher gene expression levels in females only. (B) Indel-based cis-eQTLs associated with gene expression. Two insertions (7 bp, 3L:332,512; 1 bp, 3L:332,594, r² = 0.20) are associated with markedly different expression levels in males and females. (C) cis-association overview. Plot illustrating the variant and association data for a single gene (mthl9) on a rolling window basis. The gene is shown on the top track, with UTRs in grey and coding regions in black. Significant cis-eQTLs are drawn below and color-coded by significance for each sex separately (red most significant). Linkage (r²>0.5) is shown by arcs, color-coded according to r², with higher values in red. Rows represent all 39 DGRP lines and the left column shows gene expression levels for each line and sex separately (red indicates the highest expression level and green the lowest). The grid contains a representation of variants in rolling 50 bp windows (successive windows overlapping by 45 bp) with net insertions in red, net deletions in blue, and variants not affecting the sequence size (mostly SNPs) in black. The height of each variant indicates the net size of variants with the window, up to 20 bp. The two shaded vertical bars mark the cis-eQTLs shown in (B).

https://doi.org/10.1371/journal.pgen.1003055.g004

We found that 5% of all cis-eQTL-associated genes do not exhibit detectable expression in any tissue of the corresponding sex (as represented in FlyAtlas [45]), indicating that the underlying genetic variation induces previously unreported gene expression in young adult flies. These cis-eQTLs show an enhanced effect compared to genes previously known to be expressed in adult fly tissues (Figure S14), thus strengthening the finding that these genes become strongly expressed given a specific genetic background. An example of this phenomenon involves Hsc70-2. This gene, ranking second among all female-specific associations, is known to be highly expressed in testes of adult males only [45], [46]. However, we find that females from several DGRP lines also exhibit high Hsc70-2 expression (Figure 4A), an observation which was validated using RNA-seq-based allele imbalance analyses (see below).

Of the cis-eQTL-associated genes, 53% are associated with SNPs only, 44% with SNPs and non-SNP variants, and 3% with non-SNP variants only. The latter is expected considering the increased density of SNPs near indels and complex variants (Figure 2A). In total, indels comprise 10% of all cis-eQTLs and complex variants 4%, and they have larger effect sizes than SNP cis-eQTLs (Mann-Whitney U test, P = 1.6e-10; Figure S15). One gene associated with indel cis-eQTLs is mthl9, where we identified two insertions with markedly different effects on gene expression in both males and females (Figure 4B).

Genomic distribution of cis-eQTLs

Although we considered only cis-associations, we identified 1,162 cis-eQTLs (6.6%) associated with two genes, and 58 with three genes. We found on average seven (median three) significant cis-associations per transcript (males and females combined), with the physical distance between two consecutive cis-eQTLs being 689 bp (median 121 bp), thus narrowing putatively causal variants down to a few hundred base pairs. However, cis-eQTLs associated with the expression of the same gene were found to be often in strong linkage disequilibrium (average r² = 0.88).

We examined the genomic location of significant associations relative to TSSs and TESs (Figure 3A and 3B). Both the TSS and TES show the most significant associations independent of the variant type. We observed a quasi-symmetric distribution around TSSs in marked contrast to TESs, where upstream associations are, on average, more significant than downstream associations. These findings are consistent with previous studies in other metazoans including humans [40], [47], [48]. The high density of genetic markers in this study affords however greater resolution, clearly positioning the most significant associations within a window of less than 500 bp around the TSS and 1 kb upstream of the TES, consistent with findings in yeast [49]. The results are also consistent with enrichment for significant associations in H3K4me3 regions (2.8-fold for females and 2.6-fold for males) as compared to the rest of the genome (Table S12). Non-coding conserved regions are not enriched for cis-eQTLs, in agreement with findings in humans, where eQTLs are underrepresented in such regions [50]. One explanation could be that variants in conserved non-coding regions that dramatically affected gene expression were purged from the genome. Regions marked by tri-methylated histone H3 at lysine 27 (H3K27me3) are underrepresented for cis-eQTLs, consistent with the modification being a heterochromatin mark [51]. However, while we observe an enrichment of cis-eQTLs in functionally annotated regions, the majority of cis-eQTLs are outside of previously annotated regulatory regions. High-resolution maps of these cis-eQTL-derived putative regulatory elements (Figure 4C) may in this regard constitute a powerful resource for dissecting the regulatory architecture of gene expression in further detail.

Inheritance of gene expression

To validate cis-eQTLs genome-wide, we performed mRNA sequencing (RNA-seq) of reciprocal F₁ female hybrids from two crosses involving three DGRP lines. Specifically, we aimed to identify whether transcripts predicted to be regulated by cis-eQTLs exhibit a significant allele-specific bias in gene expression. Since both alleles act in the cross in the same trans environment, differential expression in the F₁ is a direct measure of cis-regulatory activity [52], [53]. As a quality control step, we first analyzed the correlation of the expression data between the reciprocal crosses, obtaining a Spearman correlation coefficient of >99% for both F1_{362–765-initial}/F1_{362–765-reciprocal} and F1_{517–765-initial}/F1_{517–765-reciprocal}, respectively. We then tested on average ∼7,100 transcripts for allelic imbalance in the two crosses. Using a 10% FDR, we found that 6% (443) of the tested transcripts show significant allele-specific gene expression (Figure 5A) with a median allelic expression difference of 1.5-fold. Of the cis-eQTLs for which the parental lines have different alleles, and whose population mean effect on microarray expression was two-, three- and four-fold among all 39 lines, we found that on average 39%, 60% and 75% of the associated transcripts respectively exhibited a significant allelic imbalance in the crosses (Figure 5B). These results are in line with similar, albeit smaller-scale analyses in yeast [49] and mouse [48], and thus provide support for our cis-eQTL map, revealing that the majority of strong cis-eQTLs can be validated in heterozygous individuals.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Validation of cis-associations in F1.

(A) Allelic imbalance measured for ∼7,100 transcripts in F1 (362/765 and reciprocal) with RNA-seq. Dots represent the fold-change (log2) between allele-specific reads counts. Red dots indicate transcripts with significant allelic imbalance in both crosses at a false discovery rate of 10%. Circles mark transcripts that demonstrate significant allelic imbalance and that were found to be associated with cis-eQTL in females (note that only cis-eQTLs were considered when the allele between both parental lines was not the same). (B) The proportion of cis-eQTL-associated transcripts that show allelic expression imbalance in F1s scales with the strength of the cis-eQTL (P<0.001, permutation-based, see Methods for details).

https://doi.org/10.1371/journal.pgen.1003055.g005

Variant prediction

We used Release 5 of the Drosophila melanogaster reference genome for all analyses.

Gene expression analysis

Whole-adult gene expression microarray raw data files (i.e., .CEL) for 39 inbred lines were downloaded from the EBI ArrayExpress Archive (accession number E-MEXP-1594). Briefly, Ayroles et al. [6] derived inbred lines from the Raleigh (USA) population by 20 generation full-sib mating and hybridized RNA from two independent pools (25 flies/line/sex) on Affymetrix Drosophila 2.0 microarrays. The raw data set consists of gene expression measurements for males and female in two replicates. We used a four-step pipeline to analyse the microarray data:

Tissue-specific gene expression

Tissue-specific gene expression data for young adult flies was obtained from FlyBase (ftp://flybase.org/flybase/associated_files/Gelbart.2010.10.13.tar.gz). A gene was considered to be expressed in a tissue if the expression level was above 50.

Cis-associations

We grouped overlapping variants for all 39 lines as alternative alleles. For each such group of variants, we calculated associations with any transcript whose either end is within 10 kb, and for which microarray expression data was available. We employed two statistical methods: the first one, Alignment Score Association, is linear regression between the rank of alignment scores of all alleles and the rank of expression; in this context we defined alignment score to be the maximum of the number of bases inserted and that removed by the variant. Note that, since we used a non-parametric model, there was little advantage of using a more complex alignment score. This method thus associates the size of the variants rather than the exact genotype, which is important for the cases where more than two alleles have variants of different lengths. For the second method, Allele Association, we grouped the rank of expression of each line by its allele at each variant locus; we then performed a Kruskal-Wallis test. For both methods we corrected for multiple tests by repeating each test for 10,000 random permutations of gene expression, in a similar fashion to [40]. We then took for each permutation the lowest P-value for all tests in the same transcript, thus obtaining a vector of 10,000 such P-values. This vector was then sorted, and we obtained the adjusted P-value threshold from the value indexed by our unadjusted threshold. For example, the adjusted (multiple-test corrected) P-value threshold for a nominal threshold of 0.05 is the 500'th element of the sorted vector (10,000×0.05 = 500). The FDR is the ratio of the number of transcripts expected to pass the adjusted P-value threshold by chance over the number that actually passes. This, of course, implies that a different adjusted threshold is used for each transcript for the same nominal threshold. Further, since the permutations are random, each run of the above workflow will produce a slightly different list of eQTLs, since the P-value thresholds will be slightly different each time. We note that the overlap between the two association methods was more than 95% for all metrics, as both give approximately the same significance measure for bi-allelic variants, hence we only presented the results of the Allele Association method in the main manuscript.

Variant validation

We used five complementary methods to validate our variant calls. First, we compared our variant catalogue against publicly available indel data from the FLYSNPdb database, which features indels (1 bp to 360 bp) from five previously established D. melanogaster lines (but none of the DGRP lines) [20]. Second, we compared our calls with those published by Mackay et al. [7]. We used two methods (comparing all variants and comparing SNPs not within 30 bps of non-SNPs, Table S1). In many cases, SNP calls by the other study overlapped or were close to indel calls, which indicates a possible false positive since methods based on direct read alignment are confounded by indels. Third, we aligned Roche-454 whole genome sequencing reads from Mackay et al. available for all 39 lines [7] to validate our SNP calls (Massouras et al., in preparation), including those that (a) were originally supported by two thirds or more of the Illumina reads and (b) were covered by at least two Roche-454 reads (Table S2). Fourth, we sequenced mRNA from three DGRP lines at the young adult stage (see below and main manuscript for results) (Table S3). Finally, we examined 92 large indels (200 bp to 5.6 kb) in five DGRP lines using PCR after which we randomly selected 10 for subsequent Sanger sequencing to validate the predicted variant breakpoints (Table S4).

Variant validation using RNA–seq and allelic expression imbalance in F1

All fly stocks were grown at 25°C and a 12 h light-dark cycle on corn-meal fly medium. We collected virgin flies from three of the DGRP lines (lines 362, 517, and 765) during three days and we set up the following crosses: ♀ 362×♂ 517 (and reciprocal) and ♀ 517×♂ 765 (and reciprocal). We performed RNA-seq to assess transcript expression profiles of 3–5 day old F1 females. 25 flies per sample were frozen between 1–3 pm. Total RNA was extracted using the combined Trizol/RNeasy (Qiagen, http://www.qiagen.com/) protocol. RNA quality was measured using RNA Labchips and Bioanalyzer from Agilent Technologies (http://www.chem.agilent.com/). RNA-seq libraries were prepared using the RNA-True seq kit (Illumina). Prepared libraries were sequenced with an Illumina Genome Analyzer 2 DNA Sequencing Platform (GTF, Lausanne). We used the UCSC transcript annotation to derive the sequences of the annotated transcriptome from the reference genome. We mapped the reads to both the genomes and the annotated transcriptomes; any reads aligning to more than one genomic mapping coordinate were discarded. The tool we built specifically for this purpose takes as input the sequencing reads, two genome haplotypes, and the UCSC-annotated transcript coordinates in order to derive two transcriptomes. For every read in turn it attempts to align it to both transcriptome haplotypes. Reads that align to the same position in a transcript of the two haplotypes are then checked to see if the alignment overlaps with any variants. If it does, and the variants are different between the two haplotypes, the read is deemed to “support” the allele that produces the lowest number of mismatches. If the read aligns to only one transcriptome haplotype, it is also deemed to support the corresponding allele. In this way the tool measures allele-specific differential gene expression in the same line.

For the parental lines, we used the aforementioned tool, supplying it with the Stage 3 genome (see Variant prediction paragraph for more details) and the reference genome. For the crosses we gave it the two parental Stage 3 genomes. In each case, for each position with a variant we counted the number of reads best aligning to each allele. As a result, for the parental lines, this method provides a measurement of the true positive and false positive rate with regard to our variant discovery, since it measures the differential alignment of reads to the reference sequence and the target haplotype of the line. For the F1s, this method provides a measurement for differential expression for each transcript.

For variant validation, we considered a variant true positive (TP) when either all aligned reads support this variant, as opposed to the reference or when the number of reads that support it compared to that supporting the reference results in P<0.05 in a two-tailed binomial test. We considered a variant a false positive (FP) when either all aligned reads supported the reference at this position, or the number supporting the reference was significant using the same binomial test.

Allelic expression imbalance analysis

We tested for significant allelic expression imbalance between alleles in F1 using binomial exact tests (two-sided). About half of all tested transcripts were covered by at least 200 discriminative reads, thus providing a reasonable power to detect even small allelic imbalances (≥1.5-fold allelic expression differences). P-values were corrected for multiple hypothesis testing with the Benjamini & Hochberg procedure (implemented in the function p.adjust, R). Only transcripts that were tested for cis-associations in females and that were covered by at least 20 informative reads (i.e., allele1+allele2≥20) were tested for significant allelic imbalance. Transcripts that passed with 10% FDR in both reciprocal crosses were considered to have a significant allelic imbalance in gene expression.

The increase in percent-validated (Figure 5B) was tested by permutation analysis. For each “cis-eQTL strength cutoff” we obtained a random set of transcripts (from all tested cis-eQTL transcripts) and calculated how many of those where validated in F1. We repeated this sampling procedure 1,000 times and obtained an empirical P-value for each “cis-eQTL strength cutoff” by counting how many random transcript sets scored higher than the real set of transcripts and divided by the number of permutations (i.e., 1,000).

Recombination estimates

The Fiston-Lavier [59] recombination rate calculator was used to estimate the recombination rate (cM/Mb) in 1 Mb windows along each chromosome (only chromosome 4 excluded). The rate at the center of each interval was used.

Genomic distribution of cis-eQTL–associated genes

Genes were grouped into four different categories (i.e., genes with no-, male-specific-, female-specific-, and sex-unbiased cis-eQTLs) and within each category the percentage of genes located on the X chromosome was calculated. For cis-eQTL associated genes, we calculated whether the fraction of X-linked genes deviates significantly from non-cis-eQTL associated genes using the Chi-square test with Yates' correction.

Data availability

The full catalog of variants can be obtained from dgrp.epfl.ch/downloads. The RNA-seq data are available in the ArrayExpress Database using accession number E-MTAB-1266.