Effects of larval crowding on quantitative variation for development time and viability in Drosophila melanogaster

Abstract Competition between individuals belonging to the same species is a universal feature of natural populations and is the process underpinning organismal adaptation. Despite its importance, still comparatively little is known about the genetic variation responsible for competitive traits. Here, we measured the phenotypic variation and quantitative genetics parameters for two fitness‐related traits—egg‐to‐adult viability and development time—across a panel of Drosophila strains under varying larval densities. Both traits exhibited substantial genetic variation at all larval densities, as well as significant genotype‐by‐environment interactions (GEIs). GEI was attributable to changes in the rank order of reaction norms for both traits, and additionally to differences in the between‐line variance for development time. The coefficient of genetic variation increased under stress conditions for development time, while it was higher at both high and low densities for viability. While development time also correlated negatively with fitness at high larval densities—meaning that fast developers have high fitness—there was no correlation with fitness at low density. This result suggests that GEI may be a common feature of fitness‐related genetic variation and, further, that trait values under noncompetitive conditions could be poor indicators of individual fitness. The latter point could have significant implications for animal and plant breeding programs, as well as for conservation genetics.


| 8461
HORVÁTH and KaLInKa nonadditive variation-particularly for competitive traits-is still not yet well understood, especially in the context of environmental changes.
Among temperature, salinity, and other abiotic factors, natural populations are also influenced by biotic variables, such as density. Yet, the way of how these biotic factors affect phenotypic traits has received far less attention. Larval density varies greatly in natural Drosophila populations, resulting in continuously changing density-dependent selection (Leips & Mackay, 2000), and affecting life span as well as other fitness-related traits, such as egg-to-adult development time (DT) and body size (Barker & Podger, 1970;Graves & Mueller, 1993;Miller & Thomas, 1958;Prout & McChesney, 1985). Studies found that heat-shock resistance, as a physiological trait, was insensitive of density (Bubliy, Imasheva, & Loeschcke, 1998), but phenotypic variation for morphological traits increased when exposed to high density (Imasheva & Bubliy, 2003). It has also been shown that larval density is a key factor in successfully selecting for altered life span, suggesting that high density can activate adaptive variation for a fitness-related trait (Clare & Luckinbill, 1985). Here, we studied whether larval density can similarly affect the expression of phenotypic variation for two further fitness-related traits, DT, and egg-to-adult viability (EAV).
DT is a complex trait, which correlates with many life-history traits, such as adult body weight-and size, age-specific fecundity or viability (Chippindale, Alipaz, Chen, & Rose, 1997;Chippindale, Alipaz, & Rose, 2004;Nunney, 1996;Prasad, Shakara, Anitha, Rajamani, & Joshi, 2001;Zwaan, Bijlsma, & Hoekstra, 1995). DT is also important ecologically, especially in Drosophila species in which the larvae typically live in ephemeral environments experiencing strong con-specific competition (Throckmorton, 1975), and hence is expected to have a major impact on fitness. As a result, DT is considered to be subject to strong directional selection, which would consequently reduce genetic variation (Nunney, 1996). Nevertheless, fly populations harbor considerable genetic variation for DT when exposed to selection (Cortese, Norry, Piccinali, & Hasson, 2002;Fanara et al., 2006;Neyfakh & Hartl, 1993;Nunney, 1996) or when developmental genes are disrupted (Mensch et al., 2008). In this latter study, Mensch et al. (2008) also found ample amount of GEI for DT when the temperature was manipulated, suggesting that the Drosophila genome has a large potential to respond to environmental factors. EAV has also been the subject of numerous studies, which provided valuable information on the genetics and environmental sensitivity of this trait. In Drosophila melanogaster, the genetic variation underpinning EAV can be shaped by both mutationselection balance and diversifying selection in different populations (Mukai & Nagano, 1983;Tachida, Matsuda, Kusakabe, & Mukai, 1983;Tachida & Mukai, 1985;Takano, Kusakabe, & Mukai, 1987). Moreover, it has also been suggested that loci underlying EAV are located primarily in noncoding regions of the genome (Takano et al., 1987). Studies of cactophilic Drosophila showed that polymorphic inversions also play an important role in the maintenance of variation for viability (Betrán, Santos, & Ruiz, 1998;Fernandez Iriarte & Hasson, 2000;Rodriguez, Fanara, & Hasson, 1999).
Here, we measured the phenotypic and quantitative genetic variation for DT and EAV among 31 wild-derived inbred Drosophila strains in several competitive conditions: under low, medium, and high larval density treatments. A central aim of the study was to explore how a dynamic component of the biotic environment affects the expression of major life-history and fitness-related traits. We demonstrated that GEI is a common characteristic of the two measured fitness traits, appearing as both rank-order change among the lines and increased phenotypic variance. We then looked at how the evolutionary potential of the traits-in the form of expressed genetic variation-changed along the environmental gradient. We show that DT and EAV respond differently to density stress, with DT showing increasing adaptive potential with the increased level of stress.
Fly stocks were ordered from the Bloomington Stock Center and are kept in standard molasses/soy-corn flour/agar media-containing vials.
They were maintained at 18°C in the dark and transferred to a fresh vial every 2 weeks. Each strain was propagated in 10 vials.

| Experimental populations
Flies were kept on standard molasses/soy-corn flour/agar mediacontaining bottles and were maintained at 25°C, 80% humidity, and a 12-hr:12-hr light:dark schedule starting from 10 a.m. to 10 p.m. Prior to the start of the experiment, we reared flies under the above conditions for at least three generations to acclimatize them. To minimize variation associated with parental age and rearing conditions, the parents of flies used in this experiment were age-synchronized by collecting eggs within a 24-hr window and were grown at low density (10-15 flies/ml food) (Figure 1a). At 0-2 days posteclosion, 1,500-2,000 flies were used from each DGRP strain to set up large egg-laying cages (ø 100 mm, one cage/DGRP strain). We let the flies acclimatize for 2 days in the cages prior to the start of the experimental egg collections, while we fed them twice a day with fresh yeast paste, which was spread on 2.5% blackcurrant agar plates. We used blackcurrant juice for plate preparation instead of apple juice, because the dark pink color provides a better contrast for the subsequent egg counting (Figure 1b).

| Egg collection
On the first experimental day, caged flies were fed at 10 a.m. with a small amount of fresh yeast, and then were placed into carton boxes to allow egg laying in dark conditions, resulting in a better egg yield. Eggs collected on a fresh plate within the first hour (from 4 to 5 p.m.) were discarded because females tend to retain fertilized eggs in their reproductive tracts. To produce synchronized embryos (Neyfakh & Hartl, 1993), we used eggs from the subsequent 1-hr collection (from 5 to 6 p.m.). Plates were removed from the cages and embryos were carefully washed with water into egg collection chambers. Immediately after this, embryos were transferred back to the agar plates and were counted and sorted with a fine brush under a stereoscope. We repeated the assay three times on consecutive days, yielding three replicate measurements per strain.

| Density treatments
To get an estimate of the egg-to-adult competitive ability (DT and EAV) of these flies, we caged equal numbers of eggs from each strain and from a reference white-eye strain (w 1118 ) together under three different conditions-one in which the ratio of eggs to food was sufficient to support development of all the larvae (low-density treatment, 20-20 eggs/vial, 40 in total), one in which the majority of larvae could complete their development (medium-density treatment, 50-50 eggs/vial, 100 in total) and one in which the ratio of eggs to food supports a smaller fraction of adult eclosion (high-density treatment, 175-175 eggs/vial, 350 in total) (Figure 1). The white-eye strain served as a standardized intergenotype competitor in each experimental vial. These treatments, with the additional reference genotype, provide an ecologically relevant measurement of fitness (Joshi, 2001;Joshi & Mueller, 1996). All the experiments used 20-mm vials containing approximately 8 ml of the above-mentioned standard media. The number of white-and red-eyed, freshly emerged adults was recorded daily at 2 p.m. from the ninth day on, until no more flies eclosed. EAV was estimated as the proportion of emerging adults relative to the initial number of eggs seeded in each of the experimental vials (20, 50, or 175). DT was measured as the time elapsed from the day of egg collection until the emergence of the adult flies, and it was scored separately for males and females. We recorded DT for every fly until no more flies eclosed and used the mean DT across all the eclosed flies in a vial as the DT estimate for each replicate in each strain.

| Statistical analysis
Generalized linear mixed-effects models (GLMM, R package "lme4," Bates, Maechler, Bolker, & Walker, 2015) were used to determine which variables explain variation in DT and EAV. These models allow for analyzing non-normal, heterogeneous nested data with repeated measures (Zuur, Ieno, Walker, Saveliev, & Smith, 2009). We modeled DT as a Poisson-distributed variable and used log as the link function in the models. For EAV, we utilized the binomial distribution F I G U R E 1 Experimental scheme and setup. (a) To avoid age effects of the parental assay population, we set up 10 synchronously developing bottles for every DGRP strain, individually containing 50 ml standard agar media and ~500 eggs in each bottle. On the 11th day post-egg transfer, 1,500-2,000 age-synchronized flies were used to set up the assay populations. (b) We used embryos from an hour-long egg collection on each day to set up replicates of the three treatments with a logit link function. The following initial model was fitted to DT: Y ijklm = μ + L i + D j + S k + V l + L i D j + L i S k + D j S k + L i D j S k + ε ijklm , where μ is the overall mean, L i is the fixed effect of inbred strain, D j is the fixed effect of density (low, medium, or high), S k is the fixed effect of sex, V l is the random effect of the replicate vial (to account for the nonindependence among data points) and ε ijklm represents the error term. We fitted this model to a dataset that contained the DT of individual flies instead of replicate means, and thus, the residual variance corresponded to individual measurement errors. To test the effects of the above-mentioned parameters, we proceeded with dropping terms starting with the three-way interaction, to find the minimal model that explains the data best. For comparing models with and without the predictors, we used chi-square tests on the model AIC values (Akaike's information criterion; Crawley, 2007;Zuur et al., 2009;Field, Miles, & Field, 2012). For EAV, our trait of interest was the proportion of eggs that survived to adulthood. To analyze this data, we added a matrix to our data frame containing the number of survivors in each replicate vial in the first column (successes) and the number of nonsurvivors in the second column (failures). The fitted binomial GLMM for EAV did not include the sex effect because the sex of the eggs was unknown: Y ijk = μ + L i + D j + L i D j + ε ijk . A significant strain effect indicates genetic differences between the inbred strains for the trait of interest, while a significant effect of density indicates phenotypic plasticity. The interaction terms test whether inbred strains (L i D j ) or males and females (D j S k ) differ in their response to the changing environmental conditions, and indicate a genotypeby-environment interaction (GEI) if significant (Fanara et al., 2006).
We also used reduced linear mixed models to partition sources of variance and to estimate heritability for each trait and for each density treatment separately. The following model was fitted to the EAV data for all the three densities: Y ij = μ + L i + ε ij, where μ is the mean, L i tests the random effect of inbred strain (i = 1-31) and ε ij is the error term. For the log-transformed DT data, the model included the fix effect of sex and the random effect of replicate vial as well: Y ijkl = μ + L i + S j + L i S j + V k + ε ijkl . Broad-sense heritability (repeatability) was calculated as H 2 = σ 2 L ∕(σ 2 L +σ 2 LS +σ 2 V +σ 2 residual ) for DT (Ober et al., 2012), and as 4*σ 2 L ∕(σ 2 L + π 2 /3) for EAV (Davies, Scarpino, Pongwarin, Scott, & Matz, 2015), where σ 2 L is the variance of the random line effect and π 2 /3 is the variance of the logistic distribution. To eliminate the scale effects and therefore have a comparable measure of variability of the traits across the densities, we used coefficients of variation (CVs; Houle, 1992). We calculated the phenotypic coefficient of variation as CV P = 100√(σ 2 P /X), where σ 2 P is the phenotypic variance (standard deviation of the mean, SD) and X is the population mean. Coefficient of genetic variation was calculated as CV G = 100*√[(σ 2 P *H 2 )/X], where H 2 is the broad-sense heritability estimate.
A significant GEI can be the result of either changing betweenline variances or change in the rank order of genotypes (phenotypic plasticity; Fanara et al., 2006). To get a measure for the former one, and to test for inhomogeneity of variances across low, medium, and high densities, we used Levene's test, applying the robust Brown-Forsythe variant of the test (Brown & Forsythe, 1974), which uses deviations from median values instead of means. The test was applied to untransformed DT and EAV values using the following formulae (Dworkin, 2005): l iD = |y iD − y D |, where l i refers to the Levene statistic for the i th DGRP strain, y i is the median phenotype of the i th strain (across replicates), y is the population median, and subscript D indicates the density treatment (low, medium, or high density). Levene's statistic tests the null hypothesis that there is no heterogeneity among the variances at high, medium, and low densities. To test for rank-order changes across the three densities, we subtracted line-mean ordered genotype IDs for both DT and EAV, ranked them by using low-density data as a reference, and calculated the Kendall's coefficient of concordance (Kendall's W; Friedman, 1937). We also examined correlations between DT, EAV, and other traits that have been measured in the DGRP by other research groups (e.g., Ayroles et al., 2009;Durham et al., 2014;Ellis et al., 2014;Harbison et al., 2013;Unckless et al., 2015). All correlation coefficients reported in the results are Spearman's ρ. Statistical analyses were carried out in R (version 3.1.2, R Development Core Team, 2016).

| Substantial genetic variation and G × E interactions for DT
To test for possible genotype-by-environment interactions, we fitted GLMMs for both DT and EAV data. Prior to model fitting, we pruned the data and retained only complete observations (i.e., where data were available for a strain at all three densities, and at least one observation was recorded for both males and females). p Values were obtained by chi-square tests on the model AIC values, where models with and without the effect in question were compared. For DT, we found that the minimal adequate model included two of the three possible two-way interaction terms, indicating that GEI is a common, important feature of DT. The three-way interaction term including sex, genotype, and density effects was dropped as being nonsignificant (p > .05), and so was the genotype-by-sex interaction (p > .05). The genotype-by-density interaction was highly significant, indicated by both the chi-square test (ΔAIC = 7, p = 1.25E-06; Table 2) and the crossover reaction norms in Figure 2, with 7.6% of the genetic variance being due to this interaction effect. Moreover, we obtained a significant sex-by-density interaction for DT (ΔAIC = 3,

| Sex-by-density interaction for DT
We obtained a significant sex-by-density interaction for DT, whereby females completed their development faster in low-density conditions, but were slower compared with males when density was high.
To further test this interaction, we set up experiments with additional lower and higher densities for three chosen DGRP lines, to have a higher resolution of where the shift in DT occurs. Among the three lines, RAL 301 and RAL 307 followed the observed pattern of females being faster at low density in the original DT dataset, while the pattern for line RAL 555 was not so clear. We used densities of 10-10 (D1), 20-20 (D2), 50-50 (D3), 100-100 (D4), 175-175 (D5), and 250-250 flies (D6) [DGRP line + w 1118 reference competitor] and again measured DT the same way as above. To test for significant GEI, we fitted an identical GLMM model for these data as described above for DT.
We found that females were again faster in low experimental densities [e.g., D1: 9.59 ± 0.67 (M) vs. 9.46 ± 0.56 days (F); Table 3, Figure   S1). However, the models did not reveal a significant sex-by-density T A B L E 2 GLMM analysis of DT and EAV to test for significant genotype-by-environment interactions SD refers to the phenotypic standard deviation; CV P is the phenotypic coefficient of variation, calculated as CV P = 100√(σ 2 P /X), CV G is the genetic coefficient of variation, calculated as CV G = 100√(σ 2 P *H 2 /X). N is the number of flies. DT is in days, and viability is given as proportion of survivals. Individual fly values were used for DT, while the replicate averages for EAV.  (Table 3, Figure S1). The difference between mean DT of males and females at the highest density for line 301 was as substantial as 2.12 days (~50 hr; Figure S4). Additionally, we found a significant genotype-by-sex interaction (ΔAIC = 2, p = .0395), and a highly significant genotype-by-density effect as well (ΔAIC = 22, p = 7E-06). The reason for not finding the significant sex-by-density effect might be the combination of high interindividual variation in DT tested for only a few lines, as shown in Figure S1.

| Strong density-by-strain interaction for EAV
Chi-square statistics on deviances of the binomial GLMMs showed that including both density (ΔAIC = 62.1, p = 4.52E-15) and genotype effect (ΔAIC = 59.1, p = 1.53E-12) produced a significant improvement in the model fit for EAV. Fitting the genotype-by-density interaction term substantially decreased the residual variance, suggesting that the GEI term was also significant (ΔAIC = 157.8, p = 1.01E-33). In addition, we obtained the lowest AIC value for this full model. The genetic variance explained by the genotype-by-density interaction was 18.9% (Table 2).

| Test for rank-order changes and inhomogeneous variances as signs of GEI
We showed above that both DT and EAV expressed strong genotypeby-environment interactions (Tables 2 and 4). A significant GEI can be the result of either change in the rank order of lines in the different environments, or by change in the among-line variances (Fanara et al., 2006 These results are consistent with an increase in genetic variation in stressful conditions, but there is an alternative explanation. Greater developmental noise within strains could contribute to random, nongenetic differences between strains for DT, possibly producing the same signature. This does not appear to be the case here, however. At high density, the between-strain variance is significantly greater than interindividual variance (Kruskal-Wallis test; χ 2 = 926 (30 df), p = 2E-175; Figure S2), suggesting that between-strain, genetic variance is the major contributor to the signature of GEI. Furthermore, if the significant GEI is due to noise, then we expect DT measures at high density to show more variance, but be uncorrelated with those at lower densities. Instead, there is a significant positive correlation between the strain means at medium-and high-density treatments for DT (ρ = 0.60, p = 3.1E-04), indicating that DT measures at high density are genetic properties of the strain, rather than due to increased developmental noise. On the other hand, we found no correlation between low and medium densities for DT (ρ = −0.05, p = .78), and similarly, no correlation between low and high densities (ρ = −0.22, p = .25) indicating that the low-density treatment is not stressful. A correlation between medium-and high-density treatments suggests that the medium-density treatment is also stressful (100 larvae in total), although not as stressful as the high-density  Table 3 shows averages for the three DGRP lines studied (RAL 301,307,and 555). SD refers to the phenotypic standard deviation, and the studied densities were as follows: EAV data. The correspondence in rank orders was generally low for both traits in every comparison-we did not obtain any significant correlation, suggesting that substantial changes occurred in the rankings of the studied lines (see also Figure 2). For DT, we found the least agreement in rank orders between low and high densities (W = 0.167, p > .1), followed by the low-density-medium-density comparison (W = 0.197, p > .1).
The highest agreement was found for the medium and high densities (W = 0.256, p > .1). For EAV, the correspondence was generally higher

| Coefficients of variation and heritability for DT and EAV
We utilized reduced models to infer quantitative genetics parameters for both DT and EAV. These models revealed significant genotype, sex, and genotype-by-sex effects for DT (Table 4), with broad-sense heritability being similar for low and high densities (H 2 LD = 14.7%; H 2 HD = 11.4%), but increased when larval density was medium (H 2 MD = 28.4%). The highest genotype-by-sex contribution was observed for low density (7.6%), followed by high (1.5%) and medium densities (1.2%). Broad-sense heritability estimates were generally higher for EAV than for DT, but the estimates showed a different pattern: Highest heritability was estimated for low density (H 2 LD = 0.806), while both medium and high densities showed intermediate heritability values (H 2 MD = 0.546, H 2 HD = 0.645). Phenotypic (CV P ) and genetic (CV G ) coefficients of variation were also calculated for DT and EAV by using variances, trait mean values, and estimated heritability across the three environments in order to predict whether the evolutionary potential of the traits change with the stress level. CVs provided us with unbiased variance estimates correcting for the link between trait means and variances (scaling effect). For DT, we found that CV P was lowest for low density, and gradually increased to high density (CV P-LD = 24.03%; CV P-MD = 29.58%; CV P-HD = 44.36%; Table 1). We observed the same trend regardless if we used interindividual trait measure variation in the equation (as reported above), or included only single median trait values in the coefficient calculation (CV P-LD = 16.12%; CV P-MD = 22.84%; CV P-HD = 28.76%), or even if we used replicate information (CV P-LD = 20.06%; CV P-MD = 24.49%; CV P-HD = 30.99%). Similarly, we calculated the genetic coefficients of variation (CV G ) using the heritability estimates from the reduced models. We showed that CV G was lowest for low density and increased to high density (CV G-LD = 9.18%; CV G-MD = 15.68%; CV G-HD = 16.42%).
This result is an independent confirmation of the above-described increase in genetic variation. Moreover, we found this despite the increase in environmental variance along with the genetic variances ( Figure 4), suggesting that density as a stress factor can activate adaptive variation for a fitness-related trait and have a pronounced effect on evolvability.
For EAV, we observed a different pattern: CV P was the lowest under medium-density conditions and highest at high density (CV P-LD = 60.29%; CV P-MD = 55.78%; CV P-HD = 65.06%). Similar to DT, including not only median measures, but the used three or five replicates per line did not change this trend (CV P-LD = 63.38%; CV P-MD = 55.68%; CV P-HD = 66.58%), suggesting that the quality of measurements for both traits were adequate. CV G values showed high similarity for low and high densities (CV G-LD = 56.88%; CV G-HD = 53.48%), and a substantial reduction for medium density (CV G-MD = 41.13%, Table 1). Models were run on log-transformed DT data. To estimate the effects of variables, we compared model AIC values with and without the variable in question. Variance components (σ 2 ) were calculated from the full models fitting genotype, genotype-by-sex and replicate vial as random effects. Broad-sense heritabilities can be calculated as variance due to genotype divided by the full variance (Ober et al., 2012) and are marked in the table with bold numbers. *p < .05, **p < .01, ***p < .001.
T A B L E 4 Reduced mixed-effects models for DT, providing estimates of variance components and heritability

| Sex ratios
Although we did not include sex in the EAV model, we tested whether there was a deviation from the expected 50:50 sex ratio in any of the ). Not surprisingly, the per-line replicates were more consistent for medium and high densities ( Figure S4). The strongest density effect was observed for line RAL 820, where under high-density conditions there were four times more females emerging than males ( Figure S4).

| Correlation with adult-to-adult competitive fitness
DT exhibited a significant negative correlation with adult-to-adult competitive, reproductive fitness (measured using the competitive index technique in Ayroles et al., 2009) at both high (ρ = −0.56, p = .0016) and medium densities (ρ = −0.54, p = .0026), but not at low density (ρ = 0.19, p = .336), in keeping with the notion that fast egg-to-adult development in Drosophila is closely related to overall fitness (Throckmorton, 1975), especially under competitive conditions. This result also demonstrates that fitness-related genetic variation is uncovered by the stressful medium-and high-density treatments. In contrast, there was no relationship between EAV and fitness at high (ρ = 0.26, p = .16), medium (ρ = 0.17, p = .36) or low density (ρ = 0.16, p = .44). In addition, there was no correlation between EAV and DT.

| Negative correlation between larval viability and adult body size
Previous studies have found a decrease in larval viability when larger adult body size is artificially selected (Partridge & Fowler, 1993).
Similarly, our results revealed a significant negative correlation between EAV and adult body size (data obtained from Durham et al., 2014) under high-density conditions (ρ = −0.41, p = .029). To determine whether egg or larval viability was predominantly related to body size, we also measured egg viability for all 31 strains (data not shown). Observations suggested that pupal viabilities do not vary substantially between strains; therefore, we did not account for differences in pupal survival. After correcting for egg viability and contrasting only larval viability with body size, we found that the correlation is predominantly driven by viability of the larvae (larval viability: ρ = −0.49, p = .0074; egg viability: ρ = −0.09, p = .63; Figure S5). A relationship between larval viability and adult body size suggests that small larval body size can be advantageous under competitive conditions, perhaps because the larvae of small-bodied adults need to feed less in food containing metabolic waste products.

| DISCUSSION
Effects of competition and crowding have been widely studied on various aspects of Drosophila fitness. Early experiments have demonstrated the importance of frequency-dependent selection in both intraspecific (Levene, Pavlovsky, & Dobzhansky, 1954Lewontin & Matsuo, 1963) and interspecific competition (Barker & Podger, 1970), as well as how density affects the expression of life-history traits (Barker & Podger, 1970;Graves & Mueller, 1993;Miller & Thomas, 1958;Prout & McChesney, 1985). However, very few studies have assayed natu- To test these questions, we performed thorough phenotyping and found for both DT and EAV that the studied genotypes expressed F I G U R E 4 Line and residual variances from the reduced models for DT. Estimated variance components were subtracted from the reduced DT and EAV models to compare changes in genetic and environmental variances across densities. Genetic/residual variance proportions for DT were 0.198 for low density, 0.414 for medium density, and 0.129 for high density alternative phenotypes as a function of the environment, resulting in significant GEI. GEI is a key factor in evolvability, as it can maintain variation in natural populations; with GEI being present, the strength and direction of selection might be dramatically different between environments (Zajitschek & Bonduriansky, 2014). Two major forms of GEI are distinguished: significant crossover GEI that result in changes in the rank order of genotypes from one environment to another, and GEI manifested in increased overall expressed variation. The former one is the consequence of some genotypes being superior in one environment but inferior in another (Haldane, 1947), while variance increase is likely due to uncovered adaptive genetic variation. While we found evidence for the latter for DT, only crossover reaction norms were observed for EAV. It is possible that variants in genes underlying essential pathways exhibit deleterious effects that cannot be fully ameliorated by benign environmental conditions. EAV, and viability in general, may be predominantly underpinned by processes that are essential in development and physiology, and this might explain why we observed only rank-order changes for this trait. The absence of marked changes in variance with increasing densities has been also observed for body weight in D. melanogaster and D. simulans (Barker & Podger, 1970). This suggests that there is trait-to-trait variation in the propensity for expressing GEI with traits involving essential processes being more likely to show phenotypic plasticity without changes in the overall phenotypic variance.
Among the significant GEI, we observed a significant interaction between sex and larval density for DT. At low larval density, females tended to develop faster, while at high larval densities development was faster for males. Miller (1964) showed that D. melanogaster females indeed develop somewhat faster when larval density is low (<60), but found no evidence for the opposite when density was high.
Further studies have found evidence for faster development of females (Paranjpe, Anitha, Chandrashekaran, Joshi, & Sharma, 2005) and males (Mensch et al., 2008) under differing circumstances, supporting the notion that the sex that develops fastest depends on the environmental context, and presumably also the genetic architecture of the flies. The acquisition of sufficient resources to support egg production is the constraint that most likely influences the DT of females, possibly explaining slower development of females at higher densities. The adult body size of males is known to influence female mate choice with larger males receiving more matings (Friberg & Arnqvist, 2003;Pitnick, 1991), and it is possible that faster male development at higher densities reflects a tension between sexual and natural selection-the opportunity for sexual selection may be stronger in low-density populations (Coltman et al., 1999). Further studies of development at different densities together with measures of adult body size, mating success, and life-time fecundity will be needed to determine the underlying causes.
Theoretical studies indicate that temporally fluctuating environments can support the maintenance of genetic variation (Kawecki, 2007;Sasaki & Ellner, 1997;Svardal, Rueffler, & Hermisson, 2011), an idea that was first proposed in a verbal model by Fisher (1930), arguing that the expected continual degradation of a population's environment would ensure that additive genetic variance for fitness would never be exhausted. Furthermore, while the expression of hidden, potentially deleterious, genetic variation when simultaneously encountering a novel environment might act to push a small population to extinction, the expression of GEI in a large population that is experiencing a fluctuating degeneration of their environment is much less likely to do so (Reed, Lowe, Briscoes, & Frankham, 2003). In this sense, we envisage a role for GEI in the long-term stability of a population when experiencing recurring stresses, such as high density, and not just a role limited to an enhanced potential for adaptation to entirely novel environments. Several studies have found that Drosophila populations can readily respond to artificial selection for shorter DT with decreases of up to 40 hr from egg to adult (Burke et al., 2010;Chippindale et al., 1997;Nunney, 1996;Prasad et al., 2001;Zwaan et al., 1995), indicating that abundant additive genetic variance for DT is available in natural populations.
Coefficients of variation (CVs) allow us to compare the evolvability of traits in an unbiased way, correcting for potential interdependence between phenotype means and variances (Houle, 1992). By comparing CVs across different environments, we can predict selection outcomes and therefore choose optimal selection regimes, which may be of a great importance for not only breeders but also in conservation biology. We observed that exposure to high larval density increased both CV P and CV G substantially for DT. This is concordant with the observations of Gebhardt and Stearns (1992); that is, decreasing yeast concentration of fly media increases variation for DT. Developing on poorer food resource potentially stresses flies similar to high density, posing stronger competition on conspecifics for limited resources. On the other hand, EAV showed a different pattern: While both low and high densities were characterized by high CV values, the mildly stressful, medium-density conditions showed comparably little evolutionary potential.
Although it is clear that the genetic variance expressed under crowded larval conditions has adaptive potential, it is necessary to resolve the presence of genetic variance for DT at low density. Relative to morphological traits, DT is complex and highly polygenic, increasing the number of routes by which the trait could be influenced, and thereby increasing the difficulty confronting any potential canalizing mechanism. Recessive mutations will also be exposed in the inbred strains and are likely to lengthen DT to varying degrees in the different strains (Hollingsworth & Maynard Smith, 1955). In addition, the pooling of sexes could potentially inflate the variance. Although there may be reasons why genetic variance is expressed for DT at low density, the relationship between DT and competitive fitness suggests an alternative approach to the problem. If the variance that is of evolutionary importance is not variance in DT per se, but the variance in DT that contributes to genetic variance in competitive fitness, then this trait-competitive fitness-has no relevant variance expressed in DT at low density. From this perspective, the variance in competitive fitness that is expressed via DT at medium and high densities is cryptic (Dworkin, Palsson, Birdsall, & Gibson, 2003).
An increased variance in DT under high larval density was also reported by Miller (1964) and Barker and Podger (1970) (Willmore & Hallgrimsson, 2005). In support of this interpretation, we found evidence for both variance increase and developmental noise in our DT measures. Miller (1964) proposed that the increased variation in DT is a consequence of some larvae pupating prematurely under conditions of starvation with other larvae feeding for extended periods of time in nutrient-poor conditions. Such a scenario predicts a positive relationship between the length of development of a larva and its adult body size under conditions of high larval density, a relationship that has been observed when artificially selecting for decreased DT (Nunney, 1996). If true, our results would suggest that the plastic response of premature pupation has a genetic basis with variation for this trait segregating in natural populations. It is also possible, however, that some strains have a higher tolerance of larval waste products perhaps allowing them to more efficiently extract nutrients under conditions of crowding. Experimental evolution of Drosophila populations under high-density larval conditions has resulted in the evolution of larvae with an increased tolerance of urea and ammonia (Betrán et al., 1998;Shakarad et al., 2005;Shiotsugu, Leroi, Yashiro, Rose, & Mueller, 1997), possibly suggesting that strains with greater tolerance of waste products might also have a reduced DT. However, the ability to tolerate metabolic waste might instead incur a cost of decreased efficiency of nutrient extraction (Joshi & Mueller, 1996), thereby leading to either an extended DT or a reduced adult body size, or both. Hence, a relationship between tolerance of crowding and tolerance of waste products does not necessarily predict any relationship between DT and adult body size under high density, although tolerance of waste products and premature pupation response need not be mutually exclusive traits.
Although we found that DT in the context of larval crowding correlates with adult-to-adult competitive fitness, there was no such correlation for EAV. This is a counterintuitive result because the viability of any given genotype will determine the chance that this genotype survives to reproduce in the next generation, and, hence, would be expected to impact the fitness of the genotype. However, if there is a trade-off between EAV and fecundity, this would obscure the relationship between EAV and fitness. Such a trade-off might arise if EAV is strongly dependent on the mother's investment into the yolk of the egg, and if this investment limits the total number of eggs that she can lay (Einum & Fleming, 2000;Mappes & Koskela, 2007). The existence of such a relationship would imply that a significant fraction of the variation in EAV is a maternal effect (Heath, Fox, & Heath, 1999;Kern, Ackermann, Stearns, & Kawecki, 2007).
After controlling for egg viability, we found a significant negative correlation between larval viability and adult body size. Again this appears to be a counterintuitive result as it might be expected that larvae with the best chance of survival are those that have larger body sizes.
However, when competing with other larvae for food and space, those that feed for shorter periods of time, and, thus, have smaller adult body sizes, will also be less affected by the negative consequences of competition. Hence, our result hints at the existence of a tradeoff between viability and body size that is exerted in the preadult stage in Drosophila under competitive conditions. A potential tradeoff between viability and body size has been proposed to explain the persistence of small-bodied species (Blanckenhorn, 2000), and our findings suggest that such a trade-off might be limited to individual life-cycle stages.

ACKNOWLEDGMENTS
We are grateful to Tom Ellis for statistical help and for his valuable comments on the manuscript, and to Gökce Aköz for her comments on the manuscript. We thank the two anonymous reviewers whose comments/suggestions helped to greatly improve and clarify this manuscript. We also thank Iva Kelava and Tom Hill for contributing to