A Genomic Reference Panel for Drosophila serrata.

Here we describe a collection of re-sequenced inbred lines of Drosophila serrata, sampled from a natural population situated deep within the species endemic distribution in Brisbane, Australia. D. serrata is a member of the speciose montium group whose members inhabit much of south east Asia and has been well studied for aspects of climatic adaptation, sexual selection, sexual dimorphism, and mate recognition. We sequenced 110 lines that were inbred via 17-20 generations of full-sib mating at an average coverage of 23.5x with paired-end Illumina reads. 15,228,692 biallelic SNPs passed quality control after being called using the Joint Genotyper for Inbred Lines (JGIL). Inbreeding was highly effective and the average levels of residual heterozygosity (0.86%) were well below theoretical expectations. As expected, linkage disequilibrium decayed rapidly, with r2 dropping below 0.1 within 100 base pairs. With the exception of four closely related pairs of lines which may have been due to technical errors, there was no statistical support for population substructure. Consistent with other endemic populations of other Drosophila species, preliminary population genetic analyses revealed high nucleotide diversity and, on average, negative Tajima's D values. A preliminary GWAS was performed on a cuticular hydrocarbon trait, 2-Me-C28 revealing 4 SNPs passing Bonferroni significance residing in or near genes. One gene Cht9 may be involved in the transport of CHCs from the site of production (oenocytes) to the cuticle. Our panel will facilitate broader population genomic and quantitative genetic studies of this species and serve as an important complement to existing D. melanogaster panels that can be used to test for the conservation of genetic architectures across the Drosophila genus.

The availability of whole genome sequence data for Drosophila species has greatly facilitated advances in the fields of genetics and evolutionary biology. For example, the sequencing of 12 Drosophila genomes (Clark ) was instrumental to new discoveries in comparative genomics (Stark ; Sturgill ; Zhang ). The advent of affordable genome sequencing has also allowed population geneticists to characterize genomic variation within and among natural populations, improving our understanding of the complex evolutionary histories of cosmopolitan species such as D. melanogaster and D. simulans (Begun ; Lack ; Langley ; Pool ). Most recently, multiple panels of re-sequenced inbred D. melanogaster lines have become available, facilitating the molecular dissection of complex trait variation (Grenier ; Huang ; King ; Mackay ). With these populations of reproducible genotypes, researchers have used genome-wide association analysis to identify genetic variants underlying variation in a broad range of traits including physiological traits (Burke ; Dembeck ; Gerken ; Unckless ; Weber ), behaviors (Shorter ), recombination rates (Hunter ), disease susceptibility (Magwire ), and traits related to human health (Harbison ; He ; King ; Kislukhin ; Marriage ).

Just as comparative genomic and population genetic studies of adaptation (e.g., Machado ; Zhao ) have been enhanced through the availability of multi-species genome resources, quantitative genetics may also benefit from the availability of multispecies genome panels . The development of panels of re-sequenced lines for Drosophila species beyond D. melanogaster may support broader lines of inquiry such as the conservation of genetic architectures among related taxa (Yassin ). To this end, we have developed a new genomic resource for D. serrata, a member of the montium group of species. The montium group has long been regarded as a subgroup within the melanogaster species group (Lemeunier ), but has more recently been considered as a species group of its own (Da Lage ; Yassin 2013). Although montium contains 98 species (Brake and Bachli 2008) and represents a significant fraction of all known Drosophila species, there have been very few genomic investigations of its members. Recently, genomic tools have been developed for D. serrata including an expressed sequence tag (EST) library (Frentiu ), a physical linkage map (Stocker ), and transcriptome-wide gene expression datasets (Allen ; Allen ; McGuigan ). Additionally, an assembled and annotated genome sequence (Allen ) make D. serrata only the second species in the montium group with a sequenced genome after D. kikkawai (NCBI Drosophila kikkawai Annotation Release 101). Coupled to this, D. serrata is one member of the montium group that has been extensively studied in the field of evolutionary genetics.

Populations of D. serrata have been recorded from as far north as Rabaul, Papua New Guinea (4.4°N) (Ayala 1965) to as far south as Woolongong, Australia (34.3°S) (Jenkins and Hoffmann 1999). This broad latitudinal range has made D. serrata an ideal model for population studies addressing the evolution of species borders (Blows and Hoffmann 1993; Hallas ; Magiafoglou ; van Heerwaarden ) and adaptation along latitudinal clines (Allen ; Frentiu and Chenoweth 2010; Kellermann ). D. serrata has also emerged as a powerful model for the application of quantitative genetic designs to investigate sexual selection (Gosden and Chenoweth 2014; Hine ; McGuigan ).

Here, we the report development of a panel of 110 re-sequenced inbred D. serrata lines that we have called the Drosophila serrata Genome Reference Panel (DsGRP). Similar to the DGRP (Mackay ), flies were sampled from a single large natural population with the exception that D. serrata was sampled from its endemic distribution. In this initial description, we estimate the degree of heterozygosity remaining in the lines after inbreeding, show the degree to which lines are genetically related to one another, estimate genome-wide levels of nucleotide diversity, and describe patterns of linkage disequilibrium. We also demonstrate how this panel of flies can be used to genetically dissect trait variation by performing a genome-wide association analysis on variation in a cuticular hydrocarbon (CHC) trait.

Methods

Collection and inbreeding

Drosophila serrata were collected from a wild population located at Bowman Park, Brisbane Australia (Latitude: -27.45922, Longitude: 152.97768) during October 2011. We established each line from a single, gravid female before applying 20 generations of inbreeding. Inbreeding was carried out each generation by pairing virgin brothers and sisters. 100 inbred lines, out of the initial 239 iso-female lines established, survived the full 20 generations inbreeding and a further 10 lines were established after 17 generations of inbreeding.

Sequencing

We sequenced the genomes of 110 inbred lines using 100 base-pair paired-end reads with a 500 base-pair insert on an Illumina Hisequation 2000 sequencing machine. Sequencing and library preparation were carried out by the Beijing Genomics Institute. DNA from each line was isolated from a pool of at least 30 virgin female flies using a standard phenol-chloroform extraction method.

Quality control and SNP calling

We received reads from the Beijing Genomics Institute for which approximately 95% of the bases from each line had a base quality score greater than or equal to 20 (Illumina GA Pipeline v1.5). Read quality was also assessed using FastQC v0.11.2 before being mapped to the Drosophila serrata reference genome (Allen ) using BWA-mem v0.7.10 (Li 2013) and were realigned around indels using the GATK IndelRealigner v3.2-2 (McKenna ). Genotypes for every line were inferred simultaneously using the Joint Genotyper for Inbred Lines (JGIL) v1.6 (Stone 2012). This probabilistic model was especially designed for genotyping large panels of inbred lines or strains and is considered to have high accuracy (Mackay ; Stone 2012). Genotype calls with a probability lower than 99% were treated as missing genotypes.

Residual heterozygosity

The residual heterozygosity per line was estimated as the genome-wide proportion of sites that remained heterozygous after 17-20 generations of inbreeding, more specifically, we summed all of the genotype calls that were heterozygous and expressed this statistic as a percentage of all genotyped sites. In addition, for each site in the genome that differed among the inbred lines, we calculated the percentage of lines that were heterozygous for that site. Site filtering based on minor allele frequency and coverage was not performed for this analysis.

Pairwise relatedness between lines

Pairwise relatedness between lines (j and k) was estimated using the–make-grm-inbred command of GCTA v1.24.2 (Yang ), which applies the expression:where is the number of copies of the reference allele for the ith SNP for individual j and is the population allele frequency. N is the total number SNPs. Only biallelic SNPs with a read depth between 5 and 60 and a minor allele frequency above 5% were used to estimate relatedness.

Estimating population substructure

To test whether the sample of lines exhibited any underlying substructure, we used the approach of Bryc which is founded on random matrix theory. Importantly, while the genomic relatedness matrix calculated in GCTA, A = WW’/N, is normalized with element Wij = , where is the number copies of the minor allele carried by individual j (Yang ) and is the population allele frequency, the genomic relatedness matrix used by the approach of Bryc is not. For this approach, the genomic relatedness matrix, X = CC’, where C is an M x N rectangular matrix with M corresponding to the number of individuals used to estimate X and N is the number of SNPs, here the element Cij = , the number of copies of the minor allele carried by individual j. X was scaled by equation 2.8 in Bryc using values of M= 110 and N = 3,709,328. Only SNPs without missing data were used for this analysis. The number of sub-populations was determined by the number of eigenvalues larger than that expected for a random relatedness matrix given the significance threshold t = (1 + F) /2 (Bryc ). Here, t corresponds to a value of 0.993 with our expected inbreeding coefficient (F) of 0.986 (Falconer and Mackay 1996) after 20 generations of full-sib mating. Upon seeing that a small number of lines were unusually highly related (see results), we repeated this analysis after removing four lines (line IDs: 29, 134, 159, 206) to verify that the significant results were driven only by these “outliers”.

Linkage disequilibrium

We estimated linkage disequilibrium as the square of the inter-variant allele count correlation (r2) using PLINK v1.9 (Chang ). We estimated r2 in non-overlapping 500 base-pair windows across the entire genome. The analysis was performed on biallelic SNPs that had an average read depth of between 5 and 60 and a minimum minor allele frequency of 5%. The four highly related lines (r >= 0.1) were removed prior to the analysis.

Nucleotide diversity and neutrality

We estimated nucleotide diversity (π) and Tajima’s D (Tajima 1989) statistic in 50 kilobase non-overlapping sliding windows and took the mean for each of the major chromosome arms (2L, 2R, 3L, 3R, and X) using vcftools v0.1.15 (Danecek ). The SNP data used for this analysis had an average read depth across all lines of between 5 and 60 but no threshold on minor allele frequency was applied. Again, the four highly related lines were removed prior to this analysis.

Genome-wide association of female CHC expression

As proof-of-concept we performed a GWAS on a single cuticular hydrocarbon, 2-methyloctacosane, (2-Me-C28). CHCs are waxy substances that are secreted on the cuticle with 2-Me-C28 being one of a suite of CHCs that have been extensively studied in this species due to their role in species recognition, mate choice, and desiccation resistance. For each of 94 lines, we extracted CHCs from four virgin females, across two replicate vials using individual whole-body washes in 100μl of the solvent hexane. We used a standard gas chromatography method to quantify the amount of 2-Me-C28 (Blows and Allan 1998). To maintain the trait scale used in previous studies, we transformed the amount of 2-Me-C28 into a log-contrast value following Aitchison (1986), using an additional trait, the CHC 9-hexacosane, 9-C26:1, as the divisor. This transformation turns the expression of CHCs into a proportional measure and provides an internal control for other sources of variation including body size and condition.

Our GWAS contrasts with other analyses performed on the DGRP in that we model trait variation at the individual, rather than line mean level. We applied a single marker mixed effects association analysis, where the following model was fit for every biallelic SNP that had mean coverage between 5 and 60 and sample MAF of 5%:Here, CHC expression (y) of replicate individual j on genotype is modeled as a function of the mean term (), the additive fixed effect of the candidate SNP, the polygenic random effect that is captured by the genomic relatedness matrix () where has a ∼N(0,) distribution, and the residual error (). This model was specifically designed for populations where identical genotypes can be measured independently in multiple organisms such as inbred lines (Kruijer ) and was fit using AsReml-R v3 (VSN International) for a total of 3,318,503 SNPs. There are a couple of differences between the approach outlined above and other mixed modeling approaches to GWAS implemented in programs such as GEMMA (Zhou and Stephens 2012), FaST-LMM (Lippert ), and GCTA (Yang ). First, the use of individual-level opposed to line mean level observations, allows for estimation of the genomic heritability (Kruijer ). Although mapping power is unlikely to be significantly boosted through the use of individual level data per se when working with inbred lines (Kruijer ), a second aspect to our approach does potentially increase power. Our, albeit slower, approach re-estimates the polygenic variance component when each SNP tested and results in an exact calculation of Wald’s test statistic (Zhou and Stephens 2012). This contrasts with the approach used by many mixed model GWAS programs where, to save computation time, this variance component is estimated once in a null model with no fixed effect of SNP and then held constant for each SNP tested. Such an approach produces an approximate value of the test statistic which can result in power loss under some circumstances (Zhou and Stephens 2012). As these circumstances are difficult to predict beforehand, we chose to re-estimate the polygenic random effect despite the computational cost.

To increase computational speed, we nested this linear model within an R loop that allows the access of multiple cores using the “foreach” and “doMC” packages (Revolution Analytics and Steve Weston 2015). Significant SNPs were identified as those with p-values that passed Bonferroni multiple test correction –log10(p) > 7.8, we also report SNPs with –log10(p) > 5 for comparison to other Drosophila GWAS where this arbitrary threshold value is used (Mackay ). We took statistically significant SNPs and annotated them to the current version of the D. serrata genome (Allen ). If a significant SNP was located within a gene, we blasted the D. serrata gene sequence to the D. melanogaster genome to determine gene orthology using Flybase (Attrill ).

Data availability

Raw reads for all sequenced lines are available from the NCBI short read archive under Bioproject ID: PRJNA419238. The genomic relatedness matrices used in the population structure analysis are provided in File S1 (grm_full_Bryc.txt, grm_reduced_Bryc.txt, grm_full_gcta.txt, and grm_reduced_gcta.txt). We have provided the R code, Bryc.R, that implements the test for large eigenvalues (Bryc ). The SNP list used to analyze the data in this study is available from Dryad https://doi.org/10.5061/dryad.vh0h80g and also from www.chenowethlab.org/resources.html. The CHC phenotype file is provided as the supplementary file pheno.txt. The linear model used to fit the GWAS in ASREML/R model is provided in the file asreml_gwas.R.

Results and Discussion

Identification of SNPs

We established a panel of 110 inbred lines of Drosophila serrata from wild females caught from a single population in Brisbane Australia and sequenced their genomes. 100 base-pair paired-end reads were mapped to the Drosophila serrata reference genome (Allen ) with a mean coverage of 23.5 ± 0.5 reads per line. Using the Joint Genotyper for Inbred Lines (Stone 2012), we identified 15,228,692 biallelic single nucleotide polymorphisms (SNPs) applying a 99% probability threshold. 13,959,239 of these SNPs had a median coverage between 5 and 60 in which over 80% of the lines were genotyped for that variant. Most SNPs segregate at low frequencies (Figure 1) with 6,090,058 instances of singletons, where a SNP was present in only one line. The majority (62%) of the SNPs were annotated to intergenic regions of the genome while 12% of SNPs annotated to exonic regions and 26% were found to be in introns. A total of 3,709,329 SNPs met the minimum allele frequency threshold (MAF) of 5% to be used for genome-wide association analysis and the estimation of relatedness between lines.

Figure 1

The allele frequency spectrum of SNPs annotated as intergenic (red), intronic (blue), or exonic (green). Singletons (MAF = 0.009) are not shown. We identified 3,748,429 singletons from intergenic regions, 1,535,651 from intronic regions, and 754,075 from exonic regions of the genome.

Despite the application of inbreeding for many generations, inbred Drosophila lines often contain regions of residual heterozygosity (King ; Lack ; Mackay ; Nuzhdin ). After 17-20 generations of inbreeding, residual heterozygosity in our lines was very low and we observed only a small proportion of segregating sites within lines, suggesting that inbreeding had successfully fixed variation across these genomes. Of the 110 inbred lines, 104 had fewer than 2% segregating SNPs and 82 lines had less than 1% segregating SNPs (Figure 2). Across lines, the average proportion of segregating SNPs was 0.86% ± 0.11%, less than the theoretical expectation of 1.4% for lines that have experienced 20 generations of full-sib mating which corresponds to an expected inbreeding coefficient of F = 0.986 (Falconer and Mackay 1996). This slightly lower than expected level of residual heterozygosity may simply reflect sampling variation, be due to SNPs on the X chromosome, or may indicate purging of partially deleterious alleles during the inbreeding process (García-Dorado 2008; Garcia-Dorado 2012). There was no detectable difference in the fraction of heterozygous SNPs between the lines inbred for 17 and 20 generations (ANOVA: F1,108 = 0.0944, P = 0.76). This result suggests that 17 generations of inbreeding may be sufficient for future line development.

Figure 2

The distribution of residual heterozygosity as measured by the percentage of total genotyped biallelic SNPs (loci) that were called as heterozygous within each inbred line by the Joint Genotyper of Inbred Lines (JGIL).

Although there are several mechanisms that can inhibit the fixation of an allele within an inbred line, when short-read re-sequencing technology is used for genotyping, loci can falsely appear to be segregating due to the presence of paralogous genes or other repetitive DNA sequences (Treangen and Salzberg 2012). If paralogous genes are not represented in the reference genome, DNA sequences from the original gene and a divergent duplicate gene are mapped to the same region of the genome, causing the appearance of segregating loci in the population. When this phenomenon occurs, it is expected that regions of the genome with high “apparent heterozygosity” will associate with a higher read depth than the genome-wide average. Across the genome we found a weak, positive correlation between the level of residual heterozygosity and read depth (Spearman’s ρ = 0.036, P = 2.2 × 10−16). Plots of these two factors however, clearly show an alignment of regions with both high levels of residual heterozygosity and read depth, suggesting that the D. serrata reference genome could be missing some duplications (Figure 3). Alternatively, this result could be due to copy number variation among the re-sequenced lines and/or between the reference genome and the 110 lines. We hope that further work and improvement of our genome for this species will elucidate these small regions of residual heterozygosity.

Figure 3

The proportion of lines that are heterozygous at every site along the largest scaffold (scf7180000003208) of the reference genome (black) overlayed with the average read depth at each site shown in red. Peaks in red lines represent potential gene duplications that have collapsed to the same region during assembly.

Relatedness between lines

Population structure and cryptic relatedness are well known to confound genetic association studies, potentially generating false positive genotype-phenotype associations (Kittles ; Knowler ). Conceptually, these confounding factors can be described as the unobserved pedigree of the sampled individuals caused by distant relationships (Astle and Balding 2009). The sources of these relationships are varied but include population admixture, inadvertent sampling of close relatives, and the presence of shared chromosomal inversions. Fortunately, the unobserved pedigree can be estimated using marker based approaches. Here we used such an approach to estimate the pairwise relatedness between all lines in the form of a genomic relatedness matrix.

The structure of the genomic relatedness matrix shows that the majority of the DsGRP lines are unrelated as would be expected of a sample from a large, randomly mating population (Figure 4). We found a pair of lines that were 100% related to one another, most likely due to contamination either during the inbreeding process or DNA extraction and subsequent library preparation. Generally however, this panel of flies exhibits a lower level of relatedness compared to the DGRP, where 2.7% of the 20,910 possible pairs of lines had estimates of pairwise relatedness over 0.05 (Huang ) compared to 0.08% of 5,995 pairs of lines reported here. The discrepancy between the levels of relatedness between the DsGRP and DGRP is potentially due to the different demographic histories of the founding populations that generate population structure. North American populations of D. melanogaster have relatively complex demographic histories with admixture of African and European ancestors and instances of secondary contact (Pool ) compared to the endemic population of D. serrata.

Figure 4

Heat-map of the genomic relatedness matrix for each pair of the 110 inbred DsGRP lines. Each colored squared is an estimate of genomic relatedness between a pair of inbred lines estimated by GCTA (Yang ). The shade of color represents the degree of relatedness, with light blue showing low levels of relatedness and red high levels of relatedness.

Another likely explanation for the increased levels of relatedness in the DGRP is the presence of common segregating inversions. While chromosomal inversions are known to segregate in D. serrata, their frequency and number tends to increase in populations approaching the equator (Stocker ). Therefore, it may be the case that founding the DsGRP from the higher latitude of Brisbane has resulted in sampling relatively few inversions. As of yet, these lines have not been karyotyped; however the low levels of relatedness and the lack of any bimodal distribution for residual heterozyosity, such as the one found in the DGRP (Huang ), where a portion of the lines had high levels of segregating SNP loci (15–20%), suggests that segregating inversions are negligible in this population.

We performed an eigendecomposition of the genomic relatedness matrix to test for the presence of population structure using the approach outlined in Bryc . This analysis revealed two large eigenvalues ( and ) that were greater than that expected for a random relatedness matrix of equal size (Threshold = 0.993) (Figure 5). There is therefore evidence that the full set of 110 of lines contain substructure in the form of two subpopulations. We reasoned that the second large eigenvalue was likely caused by the four pairs of lines that were highly related to each other ( = 0.29, 0.38, 0.39, and 1.04; Figure 4). To test this, we repeated the analysis after randomly removing one line from each of the four pairs of closely related lines. Confirming the prediction, there was only one significantly large eigenvalue in this second analysis (). Such a result is expected when the data includes only a single population. To summarize, after the four highly-related lines have been removed, there is no clear evidence for population structure in the DsGRP that would lead to spurious genotype-phenotype associations in genome-wide association analysis.

Figure 5

Distribution of eigenvalues from an eigendecomposition of the genomic relatedness matrix for all 110 lines excluding one large eigenvalue where . Eigendecomposition of the genomic relatedness matrix, X, which was scaled by equation 2.8 in Bryc using values of N = 3,709,328 and M = 110. Here, N corresponds to the number of SNPs used to estimate X and M is the number of lines. Only SNPs with without missing data were used for this analysis. The red vertical line corresponds to the significance threshold, t, for declaring an eigenvalue larger than that expected for a random relatedness matrix. t = (1 + F) /2 and corresponds to a value of 0.993 with our expected inbreeding coefficient (F) of 0.986 after 20 generations of full-sib mating. The two largest eigenvalues were significant, (not plotted) and .

The rapid decay of linkage disequilibrium with genomic distance is a common feature of Drosophila species with r2 dropping below 0.1 within 100 base pairs (Long ; Mackay ). This allows for higher resolution mapping compared to other species such as maize and humans in which the equivalent decay does not occur until after approximately 2000 base pairs (Remington ) and 50,000 base pairs (Koch ), respectively. In the DsGRP, linkage disequilibrium decays rapidly with r2, on average, dropping below 0.1 after 75 base pairs. Surprisingly, we observe faster decay on the X chromosome compared to the autosomes (Figure 6), contrary to Mackay , despite the fact that the X chromosome has a smaller effective population size than the autosomes.

Figure 6

Decay of linkage disequilibrium (r2) between SNPs with genomic distance (bp) in the DsGRP. Values are averaged across each chromosome.

For the two cosmopolitan species of Drosophila that have been studied extensively, D. melanogaster and D. simulans, the ancestral populations from Africa consistently exhibit higher levels of polymorphism compared to the derived populations from America and Europe (Andolfatto 2001; Baudry ; Begun and Aquadro 1993; Grenier ; Lack ). Presumably, nucleotide diversity is reduced during bottleneck events associated with the colonization of new habitat. The relatively high estimates of nucleotide diversity for D. mauritiana, an endemic species from Mauritius, bolster this trend (Garrigan ). We therefore expected that our population of D. serrata, founded from the species’ ancestral range, would exhibit relatively high levels of nucleotide diversity. We estimated nucleotide diversity (π) along the major chromosome arms 2L, 2R, 3L, 3R, and X using a 50 kilobase non-overlapping sliding window approach (Figure 7). Averaged across the genome, we estimated that π = 0.0079, which is consistent with the pattern seen in other species of relatively increased levels of nucleotide diversity for populations from ancestral ranges compared to more recently established population outside of the ancestral range (Andolfatto 2001; Baudry ; Begun and Aquadro 1993; Grenier ; Lack ).

Figure 7

Boxplots of nucleotide diversity and Tajima’s D by chromosome arm. Shown are the estimates from 50 kilobase non-overlapping sliding windows with the breakdown of the number of windows per chromosome as follows: 2L = 691, 2R = 664, 3L = 641, 3R = 741, and X = 628.

Using the same sliding window approach, we tested for departures from neutrality using Tajima’s D (Figure 7) (Tajima 1989). As with other populations of Drosophila, the DGRP (Mackay ), the Zimbabwean population in the Global Diversity Lines (Grenier ), and D. mauritiana (Garrigan ), Tajima’s D was negative across the entire genome (Tajima’s D = -1.27). This is consistent with an abundance of rare alleles and could be indicative of population expansion or the occurrence of selective sweeps, however this statistic cannot distinguish the effects of demography from selection. In the DsGRP, chromosome arm 2L has higher estimates of Tajima’s D compared to the other chromosome arms (Figure 7). The causes for this pattern will hopefully be resolved through a more in-depth population genomic analysis, as chromosome 2L may potentially harbor more genomic regions that experience balancing selection, which would increase the average value of Tajima’s D.

Genome-wide association analysis of female CHC expression

A major motivation for developing the DsGRP was to begin connecting molecular variation with standing variation for some of the well-studied quantitative traits of D. serrata. Here, we applied genome-wide association analysis to the expression of the CHC 2-Me-C28 in females and identify new candidate genes that might influence trait variation. Using our mixed model approach, we found 4 SNPs that passed the 0.05 significance threshold after Bonferroni multiple test correction (Figure 8). Two of the SNPs are situated within genes, while the other two lie within 3kb of genes. We found a further 189 SNPs with p-values lower than a suggestive threshold of 10−5.

Figure 8

Genome-wide association for 2-Me-C28 expression in female Drosophila serrata. Red line indicates Bonferroni threshold corresponding to P = 0.05 and the blue indicates an arbitrary significance threshold of P = 10−5. A total of 3,318,503 biallelic SNPs were analyzed with a minimum minor allele frequency of 0.5. Alternating color shading within chromosomes indicates different contigs in the D. serrata genome assembly. Contig order has not yet been established for the genome.

We note that a GWAS performed on line mean data using the GCTA program (Yang ), which like many other mixed model GWAS applications, estimates the polygenic variance only once, detected no SNPs above Bonferroni threshold and only 34 under the arbitrary threshold of 10−5. We also applied our ASREML approach to line means rather than individuals and found that it detected exactly the same 4 SNPs above Bonferroni as our approach in equation 2 (190 SNPs had p-values lower than 10−5). It therefore appears that the increase in detection rate is in this case mainly due to the ASREML model re-estimating the polygenic variance for each SNP tested which results in an exact, rather than approximate, calculation of the test statistic (Zhou and Stephens 2012).

The majority of the literature regarding the expression of CHCs has identified genes that are related to their production within specialized cells, oenocytes. These genes constitute the major biosynthetic pathway known for CHC production and are involved with fatty-acid synthesis, elongation, desaturation, and reduction (Chertemps ; Chertemps ; Chung ; Fang ; Labeur ; Marcillac ; Wicker-Thomas ). Although none of the genes associated with the statistically significant SNPs found in this study are involved in this biosynthetic pathway, there are other biological processes involved with CHC expression, as measured by hexane washes from the cuticle. How CHCs are transported from the oenocytes to the cuticle is unknown, this study provides a potential candidate gene involved in this process. One of the significant SNPs resides in the gene Cht9, a chitinase found on chromosome 2R. Cht9, along with a number of other chitinases and imaginal-disc-growth-factors are important for the development of epithelial apical extracellular matrix, which controls the development and maintenance of wound healing, cell signaling, and organ morphogenesis in Drosophila (Galko and Krasnow 2004; Turner 2009). Knocking out the expression of Cht9 with RNAi leads to deformed cuticles, inability to heal wounds, and defects in larval and adult molting (Pesch ), and here, we provide evidence that variation in this gene may also influence other cuticular traits such as CHC abundance. Notwithstanding the small number of lines, the genome-wide association analysis presented here, combined with a previous study that identified the major role of the transcription factor POU domain motif 3 (pdm3) for polymorphic female-limited abdominal pigmentation (Yassin ), illustrate the potential of the DsGRP to discover novel regions of the genome that underpin the genetic architecture of traits.

Conclusion

We have assembled a new resource for the study of quantitative traits and population genomic variation in a non-model Drosophila species within its endemic distribution. These reproducible genotypes sampled from a single population not only provide a rich genomic dataset suitable for population genomic studies, but also provide a critical resource for the discovery of genetic variants underlying ecologically important quantitative traits. We hope that the DsGRP will provide a useful complement to other Drosophila resources such as the DGRP (Mackay ), the DSPR (King ), and the Drosophila Genome Nexus (Lack ).

In this first characterization of the DsGRP at the genomic level, we have shown that the inbreeding process has been successful in homogenizing the majority of the genome of each of the lines. Through the estimation of the genomic relatedness matrix we have shown that the DsGRP represents a random sample from a large population that contains very low levels of cryptic relatedness. These characteristics, along with rapid decay of linkage disequilibrium, make the DsGRP an ideal resource for the application of genome-wide association analysis and for generating new multifounder QTL mapping populations that will boost mapping power.

Supplementary Material

Supplemental Material is available online at www.g3journal.org/lookup/suppl/doi:10.1534/g3.117.300487/-/DC1

Click here for additional data file.