Research article | Open | Published:
A novel method for quantifying the rate of embryogenesis uncovers considerable genetic variation for the duration of embryonic development in Drosophila melanogaster
BMC Evolutionary Biologyvolume 16, Article number: 200 (2016)
Embryogenesis is a highly conserved, canalized process, and variation in the duration of embryogenesis (DOE), i.e., time from egg lay to hatching, has a potentially profound effect on the outcome of within- and between-species competition. There is both intra- and inter-specific variation in this trait, which may provide important fuel for evolutionary processes, particularly adaptation. However, while genetic variation underlying simpler morphological traits, or with large phenotypic effects is well described in the literature, less is known about the underlying genetics of traits, such as DOE, partly due to a lack of tools with which to study them.
Here, we establish a novel microscope-based assay to survey genetic variation for the duration of embryogenesis (DOE). First, to establish the potential importance of DOE in competitive fitness, we performed a set of experiments where we experimentally manipulated the time until hatching, and show that short hatching times result in priority effect in the form of improved larval competitive ability. We then use our assay to measure DOE for 43 strains from the Drosophila Genetic Reference Panel (DGRP). Our assay greatly simplifies the measurement of DOE, making it possible to precisely quantify this trait for 59,295 individual embryos (mean ± S.D. of 1103 ± 293 per DGRP strain, and 1002 ± 203 per control). We find extensive genetic variation in DOE, with a 15 % difference in rate between the slowest and fastest strains measured, and 89 % of phenotypic variation due to DGRP strain. Using sequence information from the DGRP, we perform a genome-wide association study, which suggests that some well-known developmental genes affect the speed of embryonic development.
We showed that the duration of embryogenesis (DOE) can be efficiently and precisely measured in Drosophila, and that the DGRP strains show remarkable variation in DOE. A genome-wide analysis suggests that some well-known developmental genes are potentially associated with DOE. Further functional assays, or transcriptomic analysis of embryos from the DGRP, can validate the role of our candidates in early developmental processes.
The diversity of traits directly connected to an organism’s life cycle– life-history traits– is endless, and so are the adaptive strategies that organisms can use to increase their fitness. Because adaptation depends on the genetic variation available for selection, the genetics of life-history traits has long been a major concern of evolutionary biology, and subject to serious controversy and debate. Advances in theoretical [1, 2] and empirical quantitative genetics [3–7] have changed our view of life-history traits, enabling more realistic conclusions about the components of phenotypic variation and the forces maintaining this variation. The low heritability of complex quantitative traits is not necessarily the consequence of depleted additive genetic variation, as it had been thought previously [8–10], but often the result of extensive non-additive and environmental variance ([2, 11–13]; for a review see ). While many studies have contributed to our understanding of fitness and its phenotypic trade-offs, gaps in our knowledge about the complex network of fitness characters remain.
One life-history trait that has received much attention is development time (DT) in holometabolous, ectotherm organisms, such as Drosophila and other insects [15–24]. DT is one of the best-studied and documented life-history traits, particularly regarding the effects of environmental cues, such as temperature [15–17, 23] and crowding [18–20]. Genetic studies showed that mutations in certain genes can also influence DT . Much of the interest in this trait is because it is believed to be closely related to fitness, as has been shown for Drosophila [9, 20, 25], particularly when larval competition is high . Furthermore, because Drosophila is a genetic model system, it provides an especially good opportunity for understanding the genetic basis of developmental rate, with developmental control genes and their downstream counterparts being well-documented , as are those affecting developmental rate when disrupted . Mutations resulting in hormone signaling defects in Drosophila [27–29], for example, can affect the larval-pupal transition, delaying or disrupting normal development. Moreover, Drosophila is also a useful model for understanding natural variation in developmental rate, with selection experiments revealing considerable genetic variation for DT in Drosophila [30–32].
Here, we are interested in the initial developmental stage of insects, taking place from the fertilization of an egg until hatching of the larva, or the duration of embryogenesis (DOE). DOE has often been treated as a fixed, species-specific developmental event, and variation in the rate of embryonic development has been considered as a passive, unavoidable consequence of abiotic factors such as temperature. In spite of that, there is a great variety of hatching patterns, and embryos as a responsive, evolving stage in the life-cycle have gained attention recently . There is significant between- [17, 34, 35] and within-species variation in this trait [36, 37]. Importantly, DOE varies greatly among insects, ranging from months (e.g., 4.5 months postdiapause for Aulocara at 15 °C)  to weeks (2 weeks for Grillus at 29 °C) , days (3 days at 25 °C, Heliconius ; 7 days at 25 °C, Parasteatoda ) or hours (~22 h at 25 °C for Drosophila ).
As for DT generally, Drosophila is a particularly good model for understanding genetic variation in DOE, as embryogenesis itself has been the subject of intense study [17, 26, 30, 35, 36]. But, while egg-to-adult and larval development have both been well-characterized [15, 16, 18–24], variation in DOE has been relatively neglected, in part because it is not an easy trait to measure. Previous studies quantified DOE using a technique in which eggshells of the hatching larvae were counted manually [30, 35]. As this method is very labor-intensive and prone to observational errors, it also limits work to only one or two strains/populations at a time, and experimenters have only one chance to record hatching. Given these limitations, a more reliable and repeatable method adaptable to high-throughput is clearly needed to investigate this trait, and has been developed in the present study. That is, to determine whether there is within-species genetic variation for DOE in Drosophila, and to precisely quantify it, we established a new microscope-based phenotypic assay and measured the trait in 43 DGRP strains. We found extensive genetic variation for DOE among the strains, and showed that the majority of the identified polymorphisms retard the time to hatching.
Establishing the ecological importance of embryonic development time
To establish the importance of DOE for larval competition, we experimentally manipulated differences between strains in DOE by modifying the age of the eggs in one strain (RAL555: DGRP line 555) while keeping the age fixed in the other (white-eye reference strain, w1118). We chose this DGRP line because previous experiments showed that it is similar to the reference strain with respect to egg-to-adult viability (EAV) and DT (data not shown). We tested the competitive advantage of a head-start in DOE under both high density (HD; 175 reference/175 DGRP eggs; 350 eggs in total on 8 ml standard media) and low density (LD; 50/50 eggs; 100 eggs in total). We gave the DGRP strain different levels of advantage in DOE by collecting eggs from this strain earlier than those of the control strain. We tested three levels of advantage: (a) no advantage (“0 h”), (b) 4 h of advantage (“4 h”) and (c) 7 h of advantage (“7 h”). We assessed the fitness of the two competing strains by measuring EAV and DT under these conditions, as both traits are closely related to overall fitness [9, 20, 42]. An increase in EAV and/or decrease in DT in the strain with the advantage given, or the opposite effects for the reference strain, would confirm that DOE contributes to overall fitness in D. melanogaster.
D. melanogaster fly strains for DOE measurements
The phenotypic response of 43 DGRP strains (D. melanogaster Genetic Reference Panel)  was measured, including all 40 strains of the initially sequenced and phenotyped core-set. These inbred Drosophila strains allowed us to phenotype a large number of individuals with the same genetic background, making it possible to measure a precise and robust trait mean for each strain. All flies were reared on standard molasses/soy-corn flour/agar media containing vials under low density conditions, in a single 25 °C incubator at 70–80 % relative humidity with a 12/12 h light/dark cycle.
Experimental populations and egg collection
To measure DOE, we first age-synchronized the parental population. To do this, we set up large egg-laying population cages at 25 °C (100 mm, one cage/DGRP strain), and collected eggs on 100 mm Petri dishes filled with blackcurrant juice and agar media, sprinkled with fresh yeast. Eggs were collected over a 12 to 24 h period, and set up in bottles on standard Drosophila media at low density (10–15 flies/ml food). For measuring DOE, 1 to 2 day-old flies were transferred to population cages and were acclimatized at 25 °C for 2 days prior to the experiments. On days when measurements were taken, synchronized embryos were collected for an hour (between 6 and 7 pm) on blackcurrant juice plates. Two 1-h long pre-lays preceded this egg collection in order to encourage females to lay eggs they may have incubated for varying time periods . Eggs were carefully washed off the plates into egg collection chambers, and transferred with a fine brush onto a 24-well cell-culture plate (see below) for imaging. Each well of the plate can accommodate ~100 embryos, however, numbers varied depending on the number of embryos collected. Each round of phenotyping was conducted with randomly chosen blocks of strains, with 2–4 of them assayed simultaneously on each array. We also phenotyped a reference strain (w1118, D. melanogaster) in parallel with the DGRP strains. The w1118 strain acted as a control to account for experimental noise.
Imaging plate and imaging
The transfer and subsequent imaging with a stereoscope camera provides far more accurate and repeatable measures than those obtained manually. It reduces experimental and observational error, as well as increases the throughput of the experiment, since in this way more individuals can be measured simultaneously. Five randomly distributed wells were used for each DGRP strain, and four for the reference strain. We imaged embryos on a 24-well cell culture plate in which all wells were previously filled with 1 ml 2.5 % blackcurrant agar. Using a precisely measured amount of agar in each well facilitates subsequent imaging, because one focus setting can be used throughout the experiment, and the blackcurrant juice provides a contrasting background for the images of eggs.
Sequences of images were taken at sixteen predefined time points throughout the day, between 10 am and 8 pm (Additional file 1: Table S1). Imaging started at 10 am, equivalent to 15.5 h of DT (taking the midpoint of the 1-h long egg collection as starting point). A strong distinction exists between non-retained and retained embryos, with hatching times following a bimodal distribution . Embryos that hatched before 10 am are presumably from eggs that began embryonic development while still retained in the female, and therefore were discarded from the analysis . Excluding these embryos from the analysis allowed us to avoid confounding effects due to heterogeneous egg retention behavior of the strains.
Scoring was done with image analysis software (Fiji, ) by counting empty egg shells of the hatched larvae on consecutive images (Additional file 2: Figure S1). First, a time-lapse image stack was created from the individual pictures of each batch of eggs, then this image stack was aligned using the eggs themselves as landmarks. For counting the hatched eggs, we used the “cell-counter” plugin (http://rsbweb.nih.gov/ij/plugins/cell-counter.html).
We partitioned hatched embryos into 17 developmental categories (Additional file 1: Table S1). After excluding early-hatching embryos (<15.5 h) from the analysis, and combining late-hatching ones (>25.5 h) with the last measured time point (25.5 h), we had measurements of DOE for every strain and every replicate separately. We calculated weighted mean DOE by multiplying the number of eggs hatched at a given time point with the corresponding DOE. We calculated relative DOE by subtracting the mean reference value from the strain means. Negative values represent faster DOE compared to the reference strain, while positive values mean slower development.
Quantitative genetic analysis
All statistical analyses were conducted using R version 3.1.2 (R Core Team 2014). To analyze variation in DOE among DGRP strains, we performed ANOVA on the untransformed data, since the normality test did not indicate non-normal residual distribution (Shapiro-test, W = 0.98, p = 0.07). Specifically, we fit the following linear model: Y ij = μ + L i + ε ij , where Y ij is the relative mean DOE measure for the j th replicate of the i th line, L i is the effect of line i (i = 1–43), ε ij is the error term for the i th line and the j th replicate and j is the replicate (j = 1–3). The σL 2 provides an estimate of genetic variance (σG 2), and the error (ε) provides an estimate of the environmental variance (σE 2). Phenotypic variance (σP 2) is calculated according to the formula σP 2 = σG 2+ σE 2, and broad-sense heritability (H 2) was estimated as H 2 = σG 2/ σP 2. The phenotypic and genetic coefficients of variation were calculated using raw data, as CVP = 100*√ (σP 2/X) and CVG = 100*√ [(σP 2*H 2)/X], where X and σP 2 are the mean and variance of raw DOE estimates across all strains. To perform post-hoc analysis of the data, we used the Tukey-Kramer test implemented in the HSD.test test in R (package “agricolae”).
To assess the ecological importance of DOE, we first fitted the maximal, four-way ANOVA model on the Box-Cox transformed DT. The Box-Cox transformation was used to satisfy the normality and homogeneity of variance assumptions for the residuals of the fitted model, and lambda was calculated with the powerTransform function (package “car”). The initial, full model included four main effects [a line effect, L (i = 1–2); density treatment, D (j = 2); sex, S (k = 2) and the hours of advantage given to the DGRP strain, T (l = 3)], all of their interaction terms, and an error term (ε). After sequentially dropping non-significant terms, including several interaction terms and one main effect (sex), we obtained the simplified model below, which was used for subsequent analysis: Y ijkl = μ + L i + D j + T k + L i D j + D j T k + ε ijkl , where Y ijkl is the transformed DT (Table 2A). As an independent confirmation for our model reduction, we used the step function in R, which uses AIC as the model selection criteria, and which produced the same reduced model as above. We also analyzed the data separately by density and advantage given to the DGRP strain (Additional file 3: Table S2).
To analyze the EAV data, we performed the same steps as for DT, but using the arcsine square-root transformed EAV ratios. In this case, however, sex was not included in the full model. The minimal adequate model obtained for EAV was Y ijk = μ + L i + D j + L i D j + ε ijk (Table 2B). The linear models for DT and EAV fit the transformed data well, as indicated by the non-significant Shapiro-test results performed on the model residuals (DT: Shapiro-test, W = 0.99, p = 0.84; EAV: Shapiro-test, W = 0.97, p = 0.38).
We used correlations to gain insight into the underlying basis for the genetic variation.
If the variation is maintained due to antagonistic pleiotropy—with genes affecting one trait favorably and another trait negatively—we expect negative correlations between fitness components. If it is the result of deleterious mutations exposed by inbreeding, we instead expect positive genetic correlations between fitness components . To examine such scenarios, we used a Spearman rank test for quantifying correlations between DOE and other DGRP phenotypes.
Genome wide association study
We tested for associations between SNPs and the relative DOE measures (median of the three replicates of the relative mean DOE). Single trait measures were uploaded in February 2016 to the DGRP webtool (http://dgrp2.gnets.ncsu.edu) , which performed the analysis between the trait measures and the corresponding DGRP strain sequence variants. Briefly, we tested SNPs and indel markers for an effect on the Wolbachia—and major inversion polymorphism–adjusted single trait measures, restricting the analysis to markers with a minor allele count ≥ 4 and a per-strain coverage between 2 and 30. For every marker, we fit the model Y ij = μ + M i + ε ij , where M was the effect of the polymorphic marker (SNP or indel). Multiple regression and partial least squares were used to correct for biased allelic effect estimates.
We performed calculations to assess the power of detecting large effect alleles at the level of p < 10−5 influencing DOE, following the procedure of Ivanov et al. . For this, we used the pwr.t2n.test function in R (package “pwr”), which needs minor (n 1 ) and major allele counts (n 2 ) as well as a Cohen’s effect size estimate (d). d is calculated by dividing the mean phenotypic difference between the two genotypes (μ1–μ2) by the standard deviation, σp, which is the square root of the sample variance, calculated by √[(n 1 -1)σ2 1 + (n 2 -1)σ2 2 /(n 1 + n 2 -2)] .
We also performed GO analysis with GOWINDA  to test for any enrichment in functional gene groups with SNPs affecting DOE. With a relaxed p-value of 10−4 we included 268 variants as candidates in the analysis, for which we used the program’s gene mode, included SNPs 2000 bp up/downstream of a gene, used minimum gene number of 5 and 1,000,000 simulations. To compare rates of molecular evolution of candidate DOE genes relative to others in the Drosophila genome, we used estimates of d N , d S , and d N /d S (ω) along the lineage leading to D. melanogaster from FlyBase (ftp://ftp.flybase.net/genomes/12_species_analysis/clark_eisen/paml/).
Ecological importance of embryonic development time
To assess the fitness consequences of shortened/prolonged embryogenesis, we manipulated the DOE for one of the strains, and measured DT and EAV for two competing genotypes under two larval densities (Table 1). Prior to data analysis, we corrected for the given advantages to the DGRP line in DOE (4 and 7 h), and subtracted the corresponding values from the overall DT of the experimentally delayed reference strain. The experiments revealed significant effects of density and genotype, as well as significant interactions between density and genotype and density and embryo age (Table 2). We found no significant effect of sex on either of the measured traits. The strain effect for DT was much smaller than that for EAV, measured as the fraction of phenotypic variance explained by DGRP line. This appears to be mainly due to the strong effect of density on variation in DT; when restricting the analysis to single density treatments, the strain effects were more similar between the two traits (Additional file 3: Table S2; phenotypic variance explained by strain for DT, low density: 56.8 %, DT, high density: 83.1 %; EAV, low density: 90.8 %, EAV high density: 91.6 %). To have a simple estimate of the fitness effect of embryo age, we calculated relative DT and EAV values by taking the difference between the two competing strains (DTw–DT555; EAVw–EAV555; Fig. 1). Higher relative DT values imply faster DT in the DGRP strain, while a smaller difference between the EAV of the two strains indicates increased viability for the DGRP strain (since viability of the reference strain was higher in all comparisons). Under benign, low-density conditions, DT was marginally faster in the reference strain, regardless of the advantage given to the DGRP strain [mean relative DT ± S.D. (in hours): −3.12 ± 11.74 (“0 h”); −3.43 ± 5.68 (“4 h”); −4.81 ± 2.25 (“7 h”); F1,42 = 4.300, p = 0.0443; Fig. 1, Additional file 3: Table S2A.1], but the embryo age did not appear to have a significant effect. However, under competitive, high-density conditions, ANOVA revealed significant main effects of both modified embryo age (F1,42 = 6.29, p = 0.0041) and genotype (F1,42 = 38.74, p = 1.89E-07; Fig. 1, Additional file 3: Table S2A.2). The difference between the two strains increased gradually from “0 h” to “7 h”, in keeping with the idea that DOE per se has an impact on fitness, measured as DT (Additional file 3: Table S2A.3–5). That is, increasing the advantage to the DGRP strain from “0 h” to “4 h” increased the DT difference from 34.07 ± 33.9 h (“0 h”) to 47.45 ± 26.16 h (“4 h”) (ANOVA, F1,20 = 3.13, p = 0.092). With 7 h of desynchronization, this gap became even more pronounced, the relative DT difference between the two strains increased to 56.27 ± 27.45 h (ANOVA, F1,20 = 12.72, p = 0.0019).
Development time and viability of the two competing strains (RAL555: DGRP line 555, w: w1118). Adv is the advantage given to the DGRP strain (0, 4 or 7 h), n 555 and n w represents the number of individual flies measured for each genotype. Development time is calculated as the time elapsed from the day of egg collection until the emergence of the adult flies (scored daily and for males and females separately), and viability was the ratio of emerged adult flies to the initial egg numbers (50 or 175).
On the other hand, we found no evidence for a significant embryo age effect for EAV at either low or high density (ANOVA, LD: F1,18 = 0.033, p = 0.967; HD: F1,18 = 1.633, p = 0.223), and consequently no considerable increase/decrease in viability in the competing genotypes due to the given advantage/delay (Fig. 1, Additional file 3: Table S2B). Despite the fact that viability of the DGRP strain showed an increase with the increasing levels of advantage at both low and high densities (LD: 0.56–0.58–0.6; HD: 0.4–0.41–0.45), the relative viability measures showed no trend due to greater, undirected changes in viabilities of the reference strain (Tables 1 and 2).
Genetic variation for the duration of embryogenesis in the DGRP
To characterize natural genetic variation in DOE, we measured time to hatching relative to a reference (D. melanogaster, w1118 strain) for 43 DGRP strains. In all, we were able to measure DOE for 59,295 individual embryos, with 47,435 belonging to the DGRP strains and 11,860 for the reference strain. A mean ± S.D. of 1103 ± 293 embryos was measured per DGRP strain, ranging from 209 (strain 303) to 1491 (strain 321). Mean DOE across all strains was 21.69 ± 0.65 h (Additional file 4: Figure S2). The vast majority of the strains showed similarity in their DOE, with 29 of 43 DGRP strains hatching within 1 h [±0.5 h of the reference (Fig. 2, Table 3; Additional file 4: Figure S2)]. However, we found one strain (732) that was considerably faster than the reference, showing a difference of 1.471 h, and three that were up to 1.38 h slower (427, 304 and 714; Fig. 2, Table 3).
We found significant differences between strains in their DOE (ANOVA, F1,42 = 8.311, p = 1.29E-16), with a high broad-sense heritability, H2, of 0.89. The genetic coefficient of variance (CVG) was 2.782 while the phenotypic coefficient of variance (CVP) was 2.975. The between-strain variation that underlies variation in DOE appears to be potentially complex; there were 19 different trait means among the 43 lines according to a post-hoc Tukey-Kramer test.
Using the publicly available webtool (http://dgrp2.gnets.ncsu.edu) , we conducted a GWAS analysis on 1,526,387 genetic variants found among the 43 strains, for which the minor allele was represented in four or more strains, using single-locus analyses. At P < 10−5, we found 46 variants (45 SNPs + an indel) associated with the rate of embryogenesis (Fig. 3). The small number of strains makes our study underpowered to draw many conclusions about the genetic architecture of this trait. However, our precise and robust phenotype measures based on hundreds of individual embryos per strain can identify genes and regions putatively affecting DOE and provide a basis for subsequent follow-up studies .
To assess whether there is enough power to detect alleles with a large phenotypic effect on DOE, we performed power calculations, similar to Ivanov et al. . We used mean phenotype values and variances of the two genotypes from our highest ranked SNP (SNP: 3L_2,989,334 (7/43), p = 1.26 × 10−7; Additional file 5: Figure S4, first panel; μ1 = 0.89, σ2 1 = 0.138, μ2 = −0.086, σ2 2 = 0.18). Despite using only 43 strains, we showed that we had some power to detect variants with large effects, depending on the frequency of the minor allele (MAF) –which ranged from 9 to 43 % for the 46 candidate SNPs uncovered in this analysis. For example, using a P-value cut-off of 10−5, we had 31 % power to detect a variant with an effect size equal to 1 h of difference in DOE and a MAF of 10 %. For alleles with higher MAFs, we had more power: we had 86 % power to detect variants of this effect size with a MAF of 20 %, up to 98 % power with a MAF of 40 % (Additional file 6: Figure S5).
SNPs significantly associated with variation in DOE are presented in Additional file 7: Table S3. Half of the significant SNPs (25 of 46) occurred in close proximity to each other, and are likely non-independent (or in linkage disequilibrium), resulting in a total of 27 candidate genes (Fig. 3; Additional file 7: Table S3). Among these, 11 have no known function, though the biological functions of 16 genes have been described previously. Most of the annotated polymorphisms (52 %) are intronic; if causal, the effect of these intronic SNPs may be due to their effects on gene expression, as introns can harbor regulatory elements .
We also examined the data for a relationship between DOE and the presence of the intracellular bacteria Wolbachia and common large chromosomal inversions, both of which occur in these lines. We found no effect of Wolbachia infection status on the rate of embryonic DT (ANOVA, F1,2 = 3.35, p > 0.05; Additional file 8: Figure S3). The DGRP lines have been typed for several chromosomal inversions : among the 16 inversions that have been identified in the DGRP , 11 are monomorphic in our sample of 43 lines. The remaining 5 inversions (In.2 L.t, In.2R.NS, In.3R.P, In.3R.K and In.3R.Mo) show some extent of polymorphism, and two of the 5 were associated with DOE (In.2R.NS, ANOVA, F1,3 = 3.434, p = 0.04; In.3R.K, F1,3 = 12.39, p = 0.001; Additional file 9: Figure S6). We evaluated if these inversions correspond to the genomic location of any candidate genes. Two genes, Mur89F (in In.3R.K) and Sema-2b (in In.2R.NS) occur within the significant inversions, and in both cases the identified variant is prevalent in the sample (MAF = 46.5 % and MAF = 27.5 %, respectively). However, for both inversions, very few of the lines carry the inverted karyotypes– only 4 strains harbor In.2R.NS in either homozygous or heterozygous form, and only one line carries In.3R.K. Thus, the significant association between DOE and inversion karyotype should be treated with caution.
The majority of the identified candidate genes, 15 of 27 are embryonic genes, having the highest expression in early, intermediate or late embryos among all the developmental stages (Additional file 7: Table S3; data from ). Most SNPs affected DOE negatively (Additional file 5: Figure S4). However, in two cases (one variant in the Mur89F gene, and one associated with genes CG7341 and CG32195), the derived allele was associated with a significant decrease in DOE in the lines carrying the allele (Additional file 5: Figure S4). While Mur89F is only expressed in late-stage embryos (after stage 13 ) and in later developmental stages (in larvae and pre-pupae ), both CG7341 and CG32195 are embryonic genes, with CG7341 having the highest life-time expression in 4–6 h old embryos (Additional file 7: Table S3). GO analysis (p < 10−4: 268 variants) does not reveal any significantly enriched categories after correcting for multiple testing. However, three-quarters of the 268 mutations occurred in genes that were related to development, regulation and morphogenesis (data not shown). We also studied the molecular evolution of our candidate DOE genes identified by means of the ratio of non-synonymous (d N ) to synonymous (d S ) substitutions rates (ω), and compared it to other genes in the Drosophila genome. Consistent with purifying selection acting on these genes, embryonic genes had ω < < 1, similar to that of other protein-coding sequences in the genome (median rate for embryonic gene = 0.075, other genes = 0.063, W = 67,946, p = 0.675).
Among the candidate genes with defined biological functions, there was an association with a non-synonymous mutation in the cana gene (chr2L, pos: 11,276,472, GWAS p < 2.8 × 10−6). This embryonic gene might play a crucial role in the timing of development as it is involved in biological processes such as metaphase/anaphase transition of mitotic cell cycle and microtubule-based movements. A non-synonymous SNP was found in the Ranbp9 gene (chr3R, pos: 7,238,756, GWAS p < 4.6 × 10−6), a gene involved in nuclear protein transport activity. Six of 43 strains carried four intronic SNPs in the wnt6 gene, which affects processes such as neuron differentiation and cell fate commitment. An intronic SNP was found in the hid(W) gene in 18 of 43 strains. This gene is a well-described apoptosis activator (proapoptotic) gene. Moreover, this gene is also associated with head involution, the mechanism of the rearrangement of head segments at the anterior tip of the fly embryo .
Of the significant GWAS hits, eight of 46 fall into the same, 1.2 kB long region of the 3 L chromosome arm, including the five SNPs with the lowest p-values (chr3L, pos: 2,989,334–2,990,517). This region lacks annotation, and no known genes are in the immediate proximity of these SNPs (±5000 bp). The eight SNPs form a distinct haplotype (with the same seven DGRP lines carrying all of the mutations). To investigate the possible impact of these SNPs on gene expression, we used PROMO [54, 55] to map putative transcription factor binding sites (TFBS) in our region. We compared two haplotypes: sequence of the DGRP lines with the major alleles (in all cases, the reference D. mel sequence in Flybase) and the sequence of the 7 lines carrying the 8 minor alleles. The minor haplotype shows loss of two doublesex binding sites, and the gain of several diverse TFBSs: deformed (2), tailless (2), paired (1), ftz (1) and dorsal (1), proteins with crucial functions in normal development (The Interactive Fly database). Further, of the two closest flanking genes, one (cpt2, carnitine palmitoyltransferase 2) has its peak expression in early, 0–6 h embryos. (The downstream gene, CG2113 has no detectable expression in embryos ). Though not conclusive, these results suggest a possible mechanism for the effect of these SNPs on DOE, via regulatory changes affecting cpt2 expression.
Correlations with other traits
To determine the relation between DOE and other known characteristics of these strains, we correlated DOE (single, median strain values) with other traits measured in the DGRP. We found a significant negative relationship between egg viability (measured as the ratio of hatched vs. total egg number in our phenotype measures) and DOE (ρ = −0.46, p = 0.002, Fig. 4a), implying that slow embryonic development is coupled with low egg viability. Startle response  is also negatively correlated with DOE (ρ = −0.39, p = 0.03, Fig. 4b). We found no correlation between DOE and DT , but we found that strains with the longest DOE are also the ones with the longest DT at both high and low densities (data not shown). However, short DOE did not result in shorter overall DT. Similarly, we tested for a relationship between EAV and egg viability measured for the same subset of lines . We found a strong positive correlation between the EAV and egg viability at both low (ρ = 0.49, p = 0.005, Fig. 4c), medium (ρ = 0.87, p = 1.38E-10, Fig. 4d) and high density (ρ = 0.66, p = 7.55E-05, Fig. 4e). These results suggest that egg viability may be the determining factor in overall pre-adult viability in these lines.
It has been shown previously that fly species with longer DTs tend to have larger genomes . As TE content is a major determinant of genome size, we hypothesized that overall TE content might be an important factor altering DOE. However, DOE did not correlate with either total TE content, or the number of novel or unique insertions  in these lines. We tested if there is a direct relationship between genome size  and DOE, but we did not find one (ρ = −0.012, p = 0.94), even though the subset of 43 lines was equally variable in genome size as the entire set of 211 strains (W = 4247, p = 0.51), providing us with some power to detect a possible correlation.
Egg retention– the ability of female flies to withhold fertilized eggs in their reproductive tracts –can have a major impact on DOE. Therefore, we also developed a way to measure this in the DGRP (in preparation), and tested whether the traits are correlated. We found no relationship between egg retention and DOE, when excluding early hatching embryos (<15.5 h) from the analysis (ρ = −0.15, p = 0.331). However, the relation between DOE and egg retention is significant when including early hatching embryos (ρ = −0.64, p = 5.01E-06). The lack of correlation between DOE and egg retention confirms that we chose an appropriate cut-off for excluding retained embryos, therefore our DOE measures are not confounded by the retention phenotype differences between strains.
We developed a novel assay to measure the duration of embryogenesis (DOE) and found significant genetic variation for this trait in the DGRP. Compared with the previously used, time- and labor-intensive method , our assay makes it possible to quantify DOE for large numbers of embryos and strains. While previous measurements were based on a small number of observations (44 – 113 per species in ; 81–470 per population in ), we were able to obtain robust phenotype measures with our method (mean ± S.D. of 1103 ± 293 per DGRP strain). Another advantage of imaging over the manual counting is that images are taken in rapid succession, minimizing noise and error in the measurements. Moreover, since the images can be stored for a long time, data can be subtracted and analyzed at a later time point, which is not possible with the manual counting.
We used our newly developed assay to study genetic variation in the DGRP for DOE, and to find possible candidate genes and proteins contributing to variation in the rate of embryogenesis. The DGRP resource  is valuable for studying how variation in the genome maps to phenotypic differences. The vast amount of phenotypic information available for this panel of strains, combined with genetic and transcriptomic data, makes this a powerful resource. Here, we found extensive variation in DOE in a subset of DGRP strains, translating to a 15 % difference in developmental rate between the slowest and fastest strains, comparable to the reduction of egg-to-adult DT in multigenerational selection experiments (in  17 %; in  24–32 %). Moreover, the variation was largely genetic in this experiment, with strain accounting for 89 % of the total phenotypic variance. This strong strain effect on DT has also been observed in lines containing P-element induced mutations (84 % in ). Such a strong strain effect is not unexpected for broad-sense heritabilites, measured on inbred lines under strictly controlled conditions.
Because complete genome sequences are available for all DGRP strains, we were able to conduct genome-wide association mapping to identify potential candidate genes that may influence the rate of Drosophila embryogenesis. Developmental timing is based on molecular mechanisms and molecules such as cell cycle components, cell-signaling factors and hormones . Accordingly, we can find genes with previously described biological functions among our candidates that play a part in cell-cycle (cana), receptor activity (Octbeta2R, Ranbp9, tho2), cell signaling (Wnt6, dlp) and transcription (Wnt6, Zpr1). Naturally, genes with known functions in development and morphogenesis would be strong candidates for controlling DOE. Seven of the 16 genes with known functions (Jon25Bii, cana, Wnt6, pad, hid(W), dlp and Lar) are involved in developmental processes, indicating that besides their previously described roles they may be important in modulating DOE. Moreover, the available developmental transcriptome data  suggest that the majority of the genes, 15 of the 27 have their highest expression in embryos among all the developmental stages.
We found that the majority of significant sequence variants caused retardation in DOE. However, two alleles caused decrease in DOE (Mur89F intronic variant; CG7341 intronic SNP). Mur89F is a mucin-related protein, known to be expressed in salivary glands of Drosophila embryos . As many mucins are expressed in various tissues throughout embryogenesis, their potential role in organ development has been proposed . While the function of the CG7341 gene is yet to be described, its predominant embryonic expression suggests its role lying in early developmental processes. Mur89F is also one of the two candidate genes that can be found within segregating inversions. Many fitness traits show associations with inversion genotypes [50, 60], and in concordance with this, we found association between DOE and two inversion polymorphisms. Inversions span large regions of the genome, and can have a much higher amount of substitutions then the genome average, therefore taking into account inversions is necessary when assessing the genetic basis of traits.
Genes that are involved in highly conserved functions, such as developmental patterning and morphogenesis, are expected to be highly constrained . Instead, altering the expression of crucial developmental genes through chromatin modifications or changes in transcription factor binding activity is a safer, yet effective way to control phenotypes. Consistent with this idea, most candidate SNPs, with only two exceptions, are found in non-coding gene regions (introns, up- and downstream from genes) and are potentially affecting embryogenesis via gene expression. Half of the SNPs in our sample can be found in intronic gene regions, and thus, any effect they have are likely to be mediated by gene expression . Neyfakh and Hartl  proposed, that the plastic relationship between developmental rate and temperature suggests that variation affecting embryonic development time must be common in natural populations. In agreement with this, the average minor allele frequency in our sample was 23 % (9–46.5 %).
Our results are consistent with a general trend in the literature suggesting that rapid development, at least in Drosophila, is associated with high fitness. Zwaan et al.  found that in selection experiments, slow-developing flies showed reduced viability compared to the control and fast-developing strains. Similarly, we found a negative correlation between egg viability and relative DOE, which suggests that the relationship between high viability and rapid development also holds for embryos. Further, we mimicked the advantages of rapid embryonic development and showed that this resulted in faster overall development time when larval competition was high. The advantage experienced by a strain established early in a culture can be thought of as a kind of intraspecific priority effect [63, 64]. Priority effects describe the effect of a community initiating species on the later arriving ones, and are well known and studied in ecology literature (for a review see ), in particular for plant [65–67] and aquatic communities [68, 69]. While most studies focus on interspecific interactions, the importance of intraspecific effects has been described as well [68, 70]. The observed delayed/de-synchronized development is indicative of a strong competitive interaction, and can be seen as a form of intraspecific inhibitive priority effect [63, 64]. More importantly, we also found that this effect is context- and trait-dependent: the head-start caused no difference in DT of the competing lines under benign conditions (LD), but had pronounced effects when competition was high (HD). On the other hand, we observed a small but gradual increase in viability of the line with the head-start, regardless of the experimental conditions. Such context-dependency has also been shown for Daphnia species  and planktonic protists , and can have important implications for species relying on rapid growth, such as Drosophila. Priority effects increase the competitive strength of organisms and can facilitate the establishment and adaptation of otherwise less fit genotypes . Finally, the strong positive correlation between egg viability and egg-to-adult viability in the studied lines suggests that embryonic viability is an important determinant of overall viability. Our results show that embryogenesis is important for fitness, and altered DOE can have consequences for other life-history traits.
In this work, we studied genetic variation in a trait that we argue is ecologically relevant. For this, we developed a new, microscope- and computer-based method for quantifying genetic variation for the duration of embryogenesis (DOE) in Drosophila melanogaster, using the DGRP mapping panel. By characterizing DOE for more than 59,000 embryos, we found extensive variation for this trait in a subset of 43 lines. The obtained 3-h difference between the slowest and fastest strain translates to a 15 % difference in developmental rate within a single Drosophila species, comparable to the reduction of egg-to-adult development time in previous multigenerational selection experiments. We also show that DOE itself is ecologically important, and altered DOE can cause changes in other fitness components. With our GWAS analysis, we identified genes that may influence DOE, with the majority of these being embryonic genes and playing part in cell-cycle, receptor activity and cell signaling. We also suggest the role of being involved in developmental processes for genes with unknown biological functions. Further functional assays, or transcriptomic analysis of embryos from the DGRP can validate the role of our candidates in early developmental processes.
Analysis of variance
Drosophila Genetic Reference Panel
Duration of embryogenesis
Egg-to-adult development time
Minor allele frequency
Price T, Schluter D. On the low heritability of life history traits. Evolution. 1991;45:853–61.
Houle D. Comparing evolvability and variability of quantitative traits. Genetics. 1992;130:195–204.
Armbruster P, Bradshaw WE, Holzapfel CM. Evolution of the genetic architecture underlying fitness in the pitcher-plant mosquito, Wyeomyia smithii. Evolution. 1997;51:451–8.
Fowler K, Semple C, Barton NH, Partridge L. Genetic variation for total fitness in Drosophila melanogaster. Proc R Soc B. 1997;264:191–9.
Telonis-Scott M, McIntyre LM, Wayne ML. Genetic architecture of two fitness-related traits in Drosophila melanogaster: ovariole number and thorax length. Genetica. 2005;125:211–22.
Yang J, Benyamin B, McEvoy BP, Gordon S, Henders AK, et al. Common SNPs explain a large proportion of the heritability for human height. Nat Genet. 2010;42:565–9.
Ellis LL, Huang W, Quinn AM, Ahuja A, Alfrejd B, Gomez FE, et al. Intrapopulation genome size variation in D melanogaster reflects life history variation and plasticity. PLoS Genet. 2014;10:e1004522.
Gustafsson L. Lifetime reproductive success and heritability: empirical support for Fisher's fundamental theorem. Am Nat. 1986;128:761–4.
Roff DA, Mousseau TA. Quantitative genetics and fitness: lessons from Drosophila. Heredity. 1987;59:103–18.
Mousseau TA, Roff DA. Natural selection and heritability of fitness components. Heredity. 1987;59:181–97.
Crnokrak P, Roff DA. Dominance variance: associations with selection and fitness. Heredity. 1995;75:530–40.
Falconer DS, Mackay TFC. Introduction to Quantitative Genetics. 4th ed. Longmans Green, Harlow, Essex, UK; 1996.
Houle D. The character problem in life history evolution. In: Wagner GP, editor. The Character Concept in Evolutionary Biology. Academic Press; 2001. p. 109–140.
Merilä J, Sheldon BC. Genetic architecture of fitness and nonfitness traits: empirical patterns and development of ideas. Heredity. 1999;83:103–9.
Montchamp-Moreau C. Interspecific competition between Drosophila melanogaster and Drosophila simulans: temperature effect on competitive ability and fitness components. Genet Sel Evol. 1983;15:367–78.
Gibert P, De Jong G. Temperature dependence of development rate and adult size in Drosophila species: biophysical parameters. J Evol Biol. 2001;14:267–76.
Kuntz SG, Eisen MB. Drosophila embryogenesis scales uniformly across temperature in developmentally diverse species. PLoS Genet. 2014;10:e1004293.
Miller RS. Larval competition in Drosophila melanogaster and D. simulans. Ecology. 1964;45:132–48.
Barker JSF, Podger RN. Interspecific competition between Drosophila melanogaster and Drosophila simulans: Effects of larval density on viability, developmental period and adult body weight. Ecology. 1970;51:170–89.
Horváth B, Kalinka TA. Effects of larval crowding on quantitative variation for development time and viability in Drosophila melanogaster. Ecol Evol. 2016; accepted manuscript.
Mensch J, Lavagnino N, Carreira VP, Massaldi A, Hasson E, Fanara JJ. Identifying candidate genes affecting developmental time in Drosophila melanogaster: pervasive pleiotropy and gene-by-environment interaction. BMC Dev Biol. 2008;8:78.
Fanara JJ, Folguera G, Iriarte PF, Mensch J, Hasson E. Genotype by environment interactions in viability and developmental time in populations of cactophilic Drosophila. J Evol Biol. 2006;19:900–8.
Barton M, Sunnucks P, Norgate M, Murray N, Kearney M. Co-gradient variation in growth rate and development time of a broadly distributed butterfly. PLoS One. 2014;9:e95258.
Morales-Ramos JA, Rojas MG, Shapiro-Ilan DI, Tedders WL. Developmental plasticity in Tenebrio molitor (Coleoptera: Tenebrionidae): Analysis of instar variation in number and development time under different diets. J Entomol Sci. 2010;45:75–90.
Throckmorton LH. The phylogeny, ecology and geography of Drosophila. Handb Genet. 1975;3:421–69.
Campos-Ortega J, Hartenstein V. The embryonic development of Drosophila melanogaster. New York: Springer; 1985.
Sliter TJ, Gilbert LI. Developmental arrest and ecdysteroid deficiency resulting from mutations at the dre4 locus of Drosophila. Genetics. 1992;130:555–68.
Venkatesh K, Hasan G. Disruption of the IP3 receptor gene of Drosophila affects larval metamorphosis and ecdysone release. Curr Biol. 1997;7:500–9.
Hall BL, Thummel CS. The RXR homolog Ultraspiracle is an essential component of the Drosophila ecdysone receptor. Development. 1998;125:4709–17.
Neyfakh AA, Hartl DA. Genetic control of the rate of embryonic development: selection for faster development at elevated temperatures. Evolution. 1993;47:1625–31.
Nunney L. The response to selection for fast larval development in Drosophila melanogaster and its effect on adult weight: an example of a fitness trade-off. Evolution. 1996;50:1193–204.
Cortese MD, Norry FM, Piccinali R, Hasson E. Direct and correlated responses to artificial selection on developmental time and wing length in Drosophila buzzatii. Evolution. 2002;56:2541–7.
Warkentin KM. Environmentally cued hatching across taxa: embryos respond to risk and opportunity. Integr Comp Biol. 2011;51:14–25.
Kim J, Kerr JQ, Min GS. Molecular heterochrony in the early development of Drosophila. Proc Natl Acad Sci U S A. 2000;97:212–6.
Markow TA, Beall S, Matzkin LM. Egg size, embryonic development time and ovoviviparity in Drosophila species. J Evol Biol. 2009;22:430–4.
King JC. Differences between populations in embryonic developmental rates. Am Nat. 1959;93:171–80.
Zwaan B, Bijlsma R, Hoekstra RF. Artificial selection for developmental time in Drosophila melanogaster in relation to the evolution of aging: direct and correlated responses. Evolution. 1995;49:635–48.
Fisher JR, Kemp WP, Pierson FB, Wight JR. Grasshopper Egg Development: the Role of temperature in predicting Egg Hatch. Cap. IV.2. In: Grasshoppers Integrated Pest Manager User Handbook United States Departament of Agrigulture, Animal, and Plant Health Inspection Service. 1996;1809:1-7.
Donoughe S, Extavour CG. Embryonic development of the cricket Gryllus bimaculatus. Dev Biol. 2016;411:140–56.
Aymone AC, Lothhammer N, Valente VL, Araujo AM. Embryogenesis of Heliconius erato (Lepidoptera, Nymphalidae): a contribution to the anatomical development of an evo-devo model organism. Dev Growth Differ. 2014;56:448–59.
Mittmann B, Wolff C. Embryonic development and staging of the cobweb spider Parasteatoda tepidariorum C. L. Koch, 1841 (syn.: Achaearanea tepidariorum; Araneomorphae; Theridiidae). Dev Genes Evol. 2012;222:189–216.
Stearns SC. The evolution of life histories. Oxford: Oxford University Press; 1992.
MacKay TFC, Richards S, Stone EA, Barbadilla A, Ayroles JF, et al. The Drosophila melanogaster Genetic Reference Panel. Nature. 2012;482:173–8.
Stanfield SW, Helinski DR. An improved method for collecting highly synchronous D. melanogaster eggs. Drosoph Inf Serv. 1977;52:172–3.
Schindelin J, Arganda-Carreras I, Frise E, Kaynig V, Longair M, et al. Fiji: an open-source platform for biological-image analysis. Nat Methods. 2012;9:676–82.
Ivanov DK, Escott-Price V, Ziehm M, Magwire MM, MacKay TFC, et al. Longevity GWAS using the Drosophila Genetic Reference Panel. J Gerontol A Biol Sci Med Sci. 2015;70:1470–8.
Kofler R, Schlötterer C. Gowinda: unbiased analysis of gene set enrichment for genome-wide association studies. Bioinformatics. 2012;28:2084–5.
Chow CY, Wolfner MF, Clark AG. Large neurological component to genetic differences underlying biased sperm use in Drosophila. Genetics. 2013;193:177–85.
Hershberg R, Petrov DA. Selection on codon bias. Annu Rev Genet. 2008;42:87–99.
Huang W, Massouras A, Inoue Y, Peiffer J, Ràmia M, et al. Natural variation in genome architecture among 205 Drosophila melanogaster Genetic Reference Panel lines. Gen Res. 2014;24:1193–208.
Graveley BR, Brooks AN, Carlson JW, Duff MO, Landolin JM, et al. The developmental transcriptome of Drosophila melanogaster. Nature. 2011;471:473–9.
Syed ZA, Härd T, Uv A, van Dijk-Härd IF. A potential role for Drosophila mucins in development and physiology. PLoS One. 2008;3:e3041.
Vanhook A, Letsou A. Head involution in Drosophila: Genetic and morphogenetic connections to dorsal closure. Dev Dyn. 2008;237:28–38.
Messeguer X, Escudero R, Farré D, Núñez O, Martínez J, Albà MM. PROMO: detection of known transcription regulatory elements using species-tailored searches. Bioinformatics. 2002;18:333–4.
Farré D, Roset R, Huerta M, Adsuara JE, Roselló L, Albà MM, et al. Identification of patterns in biological sequences at the ALGGEN server: PROMO and MALGEN. Nucleic Acids Res. 2003;31:3651–3.
Gregory TR, Johnston JS. Genome size diversity in the family Drosophilidae. Heredity. 2008;101:228–38.
Chippindale AK, Alipaz JA, Chen H-W, Rose MR. Experimental evolution of accelerated development in Drosophila. 1. Developmental speed and larval survival. Evolution. 1997;51:1536–51.
Prasad NG, Shakarad M, Gohil VM, Sheeba V, Rajamani M, Joshi A. Evolution of reduced pre-adult viability and larval growth rate in laboratory populations of Drosophila melanogaster selected for shorter development time. Genet Res. 2000;76:249–59.
Moss EG. Heterochronic genes and the nature of developmental time. Curr Biol. 2007;17:R425–34.
Hoffmann AA, Rieseberg LH. Revisiting the impact of inversions in evolution: from population genetic markers to drivers ofadaptive shifts and speciation? Annu Rev Ecol Evol Syst. 2008;39:21–42.
Castellano D, Coronado-Zamora M, Campos JL, Barbadilla A, Eyre-Walker A. Adaptive evolution is substantially impeded by Hill–Robertson interference in Drosophila. Mol Biol Evol. 2016;33:442–55.
Millar DZ, Horan M, Chuzhanova NA, Cooper DN. Characterisation of a functional intronic polymorphism in the human growth hormone (GH1) gene. Hum Gen. 2010;4:289–301.
Fukami T. Historical contingency in community assembly: integrating niches, species pools, and priority effects. Ecol Evol Syst. 2015;46:1–23.
Devevey G, Dang T, Graves CJ, Murray S, Brisson D. First arrived takes all: inhibitory priority effects dominate competition between co-infecting Borrelia burgdorferi strains. BMC Microbiol. 2015;15:61.
Fukami T, Bezemer TM, Mortimer SR, van der Putten WH. Species divergence and trait convergence in experimental plant community assembly. Ecol Lett. 2005;8:1283–90.
Helsen K, Hermy M, Honnay O. A test of priority effect persistence in semi-natural grasslands through the removal of plant functional groups during community assembly. BMC Ecol. 2016;16:22.
von Gillhaussen P, Rascher U, Jablonowski ND, Plückers C, Beierkuhnlein C, et al. Priority effects of time of arrival of plant functional groups override sowing interval or density effects: a grassland experiment. PLoS One. 2014;1:e86906.
Sefbom J, Sassenhagen I, Rengefors K, Godhe A. Priority effects in a planktonic bloom-forming marine diatom. Biol Lett. 2015;11:20150184.
Almany GR. Priority effects in coral reef fish communities of the Great Barrier Reef. Ecology. 2004;85:2872–80.
Dibble CJ, Hall SR, Rudolf VHW. Intraspecific priority effects and disease interact to alter population growth. Ecology. 2014;95:3354–63.
We are grateful to Brian Charlesworth, Claus Vogl and Tom Ellis for their valuable comments on the manuscript. We thank Christian Schlötterer for providing us with materials and analysis tools. We thank Francois Mallard for providing us with the Java code of the Fiji implemented plugin to create image stacks.
AJB is funded by the FWF grant P27048. BH is supported by intramural funding of the Vetmeduni Vienna awarded to C. Schlötterer.
Availability of data and materials
The datasets supporting the conclusions of this article are included within the article and its additional files.
BH and ATK conceived of the experiments. BH measured DOE, DT and EAV in the DGRP, and BH and AJB analyzed the data. BH wrote the manuscript with support from AJB and ATK. All authors read and approved the final manuscript.
No competing interests declared.
Consent for publication
Ethics approval and consent to participate
Assay of measuring the duration of embryogenesis—imaging time points and corresponding embryonic development times (PDF 51 kb)
Time-lapse images of a sub-replicate measuring the embryonic development time of a DGRP strain (335) (PDF 719 kb)
Results of ANOVA on DT and EAV (PDF 118 kb)
Representation of the raw phenotypic variation for embryonic development time in 43 DGRP strains (PDF 119 kb)
Box plots of the significant SNP associations in the GWA analysis (PDF 94 kb)
Results of the power calculations to detect associated SNPs with effect size of 1 h of difference in the duration of embryogenesis (PDF 110 kb)
Top GWAS results for embryogenesis length and the corresponding DGRP phenotypes (PDF 124 kb)
The effects of Wolbachia pipientis infection on the measured phenotype (PDF 111 kb)
Inversion polymorphisms and their effect on DOE (PDF 170 kb)