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

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.

Introduction

The potential for a population to cope with environmental change is contingent upon the maintenance of additive genetic variation for fitness. Although natural selection is expected to erode such variation, several studies have shown that there is abundant genetic variation for fitness‐related traits in natural populations (Fanara, Folguera, Iriarte, Mensch, & Hasson, 2006; Ledón‐Rettig, Pfennig, & Crespi, 2010; McGuigan, Nishimura, Currey, Hurwit, & Cresko, 2010; Paaby et al., 2015; Telonis‐Scott, McIntyre, & Wayne, 2005). This puzzling fact has been explained various ways (e.g., Barton & Turelli, 1989). A potential explanation lies in the highly polygenic nature of life‐history traits; it creates the opportunity for ample amount of nonadditive genetic variation—dominance, epistatic, and genotype‐by‐environment interactions (GEIs)—to exist. Recently, there have been studies quantifying dominance effects and epistasis (Armbruster, Bradshaw, & Holzapfel, 1997; Huang et al., 2012; Monnahan & Kelly, 2015), as well as the role of GEI in adaptation (e.g., Gutteling, Riksen, Bakker, & Kammenga, 2007; Jarosz & Lindquist, 2010; Rohner et al., 2013; for a review see Schlichting, 2008). Despite that, the pervasiveness of additive‐ and nonadditive variation—particularly for competitive traits—is still not yet well understood, especially in the context of environmental changes.

It has been proposed that environmental perturbations can increase both phenotypic and additive genetic variance for fitness traits (Parsons, 1987). To test this hypothesis and to study how changing environmental conditions can affect traits, empirical studies used several abiotic environmental factors to sensitize phenotypes. Among these, there is heat shock (Bateman, 1959; Berger, Bauerfeind, Blanckenhorn, & Schäfer, 2011; Rendel, 1959; Rutherford & Lindquist, 1998; Suzuki & Nijhout, 2006; Waddington, 1953), chemical agents, such as Hsp90 inhibitors (Jarosz & Lindquist, 2010; Queitsch, Sangster, & Lindquist, 2002; Rohner et al., 2013) or ether vapor (Gibson & Hogness, 1996; Waddington, 1956). Less extreme stresses have also been shown to increase phenotypic variation, such as altered salinity levels for the expression of body size variation in sticklebacks (McGuigan et al., 2010). Although mainly morphological traits were initially surveyed, some studies looked at the environmental sensitivity of fitness and fitness‐related traits too (Fowler, Semple, Barton, & Partridge, 1997; Gardner, Fowler, Barton, & Partridge, 2005; Paaby et al., 2014).

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 mutation–selection 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.

Materials and Methods

Drosophila melanogaster fly strains

We measured the phenotypic response of 31 DGRP strains (D. melanogaster Genetic Reference Panel, Mackay et al., 2012), which are inbred strains derived from a wild North American population with completely sequenced genomes. The DGRP strains are characterized by high genetic (Mackay et al., 2012; Massouras et al., 2012) and phenotypic variation (Ayroles et al., 2009; Durham, Magwire, Stone, & Leips, 2014; Ellis et al., 2014; Harbison, McCoy, & Mackay, 2013; Mackay et al., 2012; Unckless, Rottschaefer, & Lazzaro, 2015). 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 media‐containing 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).

Figure 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

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 (w1118) 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 with a logit link function. The following initial model was fitted to DT: Y = μ + L  + D  + S  + V  + L D  + L S  + D S  + L D S  + ε, where μ is the overall mean, L is the fixed effect of inbred strain, D is the fixed effect of density (low, medium, or high), S is the fixed effect of sex, V is the random effect of the replicate vial (to account for the nonindependence among data points) and ε 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  = μ + L  + D  + L D  + ε. 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 D ) or males and females (D S ) differ in their response to the changing environmental conditions, and indicate a genotype‐by‐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  = μ + L  + ε where μ is the mean, L tests the random effect of inbred strain (i = 1–31) and ε 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  = μ + L  + S  + L S  + V  + ε. Broad‐sense heritability (repeatability) was calculated as H 2  =  ) for DT (Ober et al., 2012), and as 4*  + π2/3) for EAV (Davies, Scarpino, Pongwarin, Scott, & Matz, 2015), where 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 CVP  = 100√(/X), where is the phenotypic variance (standard deviation of the mean, SD) and X is the population mean. Coefficient of genetic variation was calculated as CVG  = 100*√[(*H 2)/X], where H 2 is the broad‐sense heritability estimate.

A significant GEI can be the result of either changing between‐line 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  = |y  − y |, where l refers to the Levene statistic for the i DGRP strain, y is the median phenotype of the i 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).

Results

Larval density affects mean DT and EAV

We scored a total of 24,028 flies derived from 48,050 eggs for EAV and DT under low (40 larvae, LD), medium (100 larvae, MD), and high (350 larvae, HD) larval densities. Overall mean EAV was 49.5%, which declined at high density (EAVLD  = 56.2 ± 23% (±SD), EAVMD  = 56.1 ± 18%, and EAVHD  = 36.6 ± 16%). Mean DT increased with increasing larval densities (DTLD  = 10.13 ± 0.58 days, DTMD  = 10.8 ± 0.94 days, and DTHD  = 14.38 ± 2.85 days [Tables 1 and S1]).

Table 1

Summary statistics for development time (DT) and viability (EAV)

Trait	Sex	Density:	LD (20–20)	MD (50–50)	HD (175–175)
Development time	M	Mean	10.185	10.832	14.358
		SD	0.537	0.979	2.822
		CVP (%)	22.957	30.068	44.333
		CVG (%)	10.670	15.529	16.833
		H 2	0.216	0.267	0.144
		N	804	1,256	2,847
	F	Mean	10.078	10.722	14.631
		SD	0.621	0.905	2.876
		CVP (%)	24.827	29.056	44.332
		CVG (%)	11.710	16.245	17.492
		H 2	0.222	0.313	0.156
		N	881	1,381	3,110
Viability	Merged	Mean	0.562	0.561	0.366
		SD	0.23	0.18	0.16
		CVP (%)	63.380	55.676	65.062
		CVG (%)	56.884	41.131	53.478
		H 2	0.806	0.546	0.645
		N	1,685	2,637	5,957

SD refers to the phenotypic standard deviation; CVP is the phenotypic coefficient of variation, calculated as CVP = 100√(/X), CVG is the genetic coefficient of variation, calculated as CVG = 100√(*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.

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, p = .027), whereby females developed faster under low‐density conditions and slower at high‐density [mean±SE: HD: 14.36 ± 0.11 (M) vs. 14.63 ± 0.17 days (F); MD: 10.83 ± 0.13 (M) vs. 10.72 ± 0.18 days (F); LD: 10.19 ± 0.12 (M) vs. 10.08 ± 0.15 days (F)]. We performed model fitting on the mean‐centered DT data as well, but the results remained unchanged, and thus, we report results on the original data. By visually inspecting it, the model residuals did not deviate from a normal distribution. Effect sizes were calculated on the untransformed DT values. Density expectedly had a major effect on DT (ΔAIC = 730, p = 2.95E‐40; Table 2), increasing it from low density to medium by 0.66 days (±0.18, SE), and by 4.48 days (±0.17) when density was high. Similarly, the strain effect was also significant itself (ΔAIC = 100, p = 1.66E‐22) (Table 2, Figures 2 and 3).

Table 2

GLMM analysis of DT and EAV to test for significant genotype‐by‐environment interactions

Trait	Variable	Δ AIC	p Value
Development time	Sex × density × genotype	−106	1
	Sex × genotype	−23	.2124
	Sex × density	3	.0274*
	Density × genotype	7	1.25E‐06***
	Sex	2	.0542
	Genotype	100	1.66E‐22***
	Density	730	2.95E‐40***
	Vial	−2	1
Viability	Density × genotype	157.8	1.01E‐33***
	Genotype	59.1	1.53E‐12***
	Density	62.1	4.52E‐15***

To test for significant GEI, we step by step removed interaction terms from a saturated model and retained variables if they led to a significant increase in AIC value. Δ AIC indicates the change in AIC value when a particular interaction term is removed from the model. To get significance for the main effects, we compared a model with all the main effects and the model lacking the particular variable. We retained the random experimental vial effect in every model in order to utilize the same mixed‐effects models for the comparisons.

*p < .05, ***p < .001.

Figure 2

Between‐strain variation for DT (a, c) and EAV (b, d) at low, medium, and high larval densities. (a) and (b) show the overall variation with medians taken across the replicates for each strain (five replicates for low, three replicates for medium and high densities). (c) and (d) illustrate changes in the reaction norms of the 31 DGRP strains using median strain values (sexes pooled for EAV). The nonparallel reaction norms indicate genotype‐environment interactions. LD, low density; MD, medium density; HD, high density; M, males; F, females

Figure 3

Phenotypic variation for DT and EAV across 31 DGRP strains showing values for five or three replicates in each strain and ranking strains by median phenotype values. (a) EAV, low‐density treatment. (b) EAV, medium‐density treatment. (c) EAV, high‐density treatment. (d) DT at low density, (e) DT at medium density, (f) DT at high density. LD, low density; MD, medium density; HD, high density

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 + w1118 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 interaction (ΔAIC = −7, p > .1), despite the observation that above the density D3 (100–100 flies) females tended to slow down in two of the three assayed lines (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.

Table 3

Summary statistics for the density experiments

Trait	Sex	Density	D1	D2	D3	D4	D5	D6
Development time	M	Mean	9.594	9.956	10.688	11.803	14.382	16.795
		SD	0.665	0.633	0.788	1.265	2.379	3.331
	F	Mean	9.457	9.730	10.679	11.925	14.433	17.527
		SD	0.561	0.574	0.771	1.319	2.389	3.397
Viability	Merged	Mean	0.744	0.728	0.653	0.597	0.523	0.366
		SD	0.224	0.177	0.131	0.086	0.095	0.162

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: 10–10 (D1), 20–20 (D2), 50–50 (D3), 100–100 (D4), 175–175 (D5), and 250–250 flies (D6) [DGRP line + w1118 reference competitor].

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 genotype‐by‐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). To test for inhomogeneity of variances across the three treatments, we applied Levene's test to both EAV and DT. DT showed significantly higher variance at both high (F 1,59 = 20.09, p = 3.46E‐05) and medium (F 1,59 = 5.84, p = .0186) densities when compared to low density. The contrast remained for DT when comparing high density with medium density (F 1,60 = 8.58, p = .0048; Figure 2a,c). For EAV, the effect of larval density was not as apparent; EAV showed no difference in variance between low‐ and medium‐density (F 1,59 = 2.59, p = .113; Figure 2b,d) or medium‐ and high‐density treatments (F 1,60 = 1.12, p = .295). However, the variance was significantly elevated in the low‐density treatment when compared to the high‐density treatment (F 1,59 = 7.68, p = .007).

Table 4

Reduced mixed‐effects models for DT, providing estimates of variance components and heritability

Trait	Density	Variable	Δ AIC	p Value	σ2 (%)
Development time	LD	Sex × genotype	−24.9	2.089E‐07***	7.6
		Sex	−11.5	2.458E‐04***	Fixed
		Genotype	−68.0	6.173E‐17***	14.7
		Vial	−6.7	.0032**	3.7
		Residuals			74.0
	MD	Sex × genotype	−2.9	.027*	1.2
		Sex	−13.0	1.096E‐04***	Fixed
		Genotype	−81.0	8.42E‐20***	28.4
		Vial	−5.9	.0048**	1.8
		Residuals			68.6
	HD	Sex × genotype	−26.2	1.106E‐07***	1.5
		Sex	−20.3	2.426E‐06***	Fixed
		Genotype	−62.1	1.171E‐15***	11.4
		Vial	−30.5	1.211E‐08***	1.6
		Residuals			85.5

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.

We used the Levene's test to further test whether the within‐strain interindividual variance was greater at high density than at low or medium density for DT. After applying a Bonferroni correction, 30 of 31 strains showed significantly greater variance between individuals at high vs. low, as well as high density vs. medium density for DT (Figures S2 and S3, Table S2; no equivalent test can be applied to EAV as each estimate is an aggregate over all of the individuals that eclose). We also found significant differences in variance for 10/30 lines between the two lower densities (Figures S2 and S3, Table S2). Greater interindividual variance at high density for DT indicates that high larval density elicits developmental noise in the rate of development of an individual larva (Willmore & Hallgrimsson, 2005).

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 treatment: The standard deviation of hatching times at medium density (σ = 0.56 days) is only half that at high density (σ = 1.19 days). In contrast, there is a nearly fivefold difference between low and high densities (low density: σ = 0.26 days). Such a dramatic increase in hatching delay is likely to have significant consequences in natural populations.

To test whether the observed GEI results from rank‐order changes among the strains, we calculated Kendall's W for all the three possible comparisons (LD vs. MD, LD vs. HD, and MD vs. HD), for both DT and 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 than for DT, and all comparisons resulted in similar values (LD‐MD: W = 0.337, p > .1; LD‐HD: W = 0.346, p > .1; MD‐HD: W = 0.337, p > .1).

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 ( = 14.7%;  = 11.4%), but increased when larval density was medium ( = 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 ( = 0.806), while both medium and high densities showed intermediate heritability values (  = 0.546,  = 0.645).

Phenotypic (CVP) and genetic (CVG) 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 CVP was lowest for low density, and gradually increased to high density (CVP‐LD = 24.03%; CVP‐MD = 29.58%; CVP‐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 (CVP‐LD = 16.12%; CVP‐MD = 22.84%; CVP‐HD = 28.76%), or even if we used replicate information (CVP‐LD = 20.06%; CVP‐MD = 24.49%; CVP‐HD = 30.99%). Similarly, we calculated the genetic coefficients of variation (CVG) using the heritability estimates from the reduced models. We showed that CVG was lowest for low density and increased to high density (CVG‐LD = 9.18%; CVG‐MD = 15.68%; CVG‐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.

Figure 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

For EAV, we observed a different pattern: CVP was the lowest under medium‐density conditions and highest at high density (CVP‐LD = 60.29%; CVP‐MD = 55.78%; CVP‐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 (CVP‐LD = 63.38%; CVP‐MD = 55.68%; CVP‐HD = 66.58%), suggesting that the quality of measurements for both traits were adequate. CVG values showed high similarity for low and high densities (CVG‐LD = 56.88%; CVG‐HD = 53.48%), and a substantial reduction for medium density (CVG‐MD = 41.13%, Table 1).

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 three environments. Overall, there were more females eclosing from the experimental vials than males (52.26% [F], 47.74% [M]). This ratio is unchanged when the data are split into low density (52.28% [F], 47.72% [M]), medium density (52.37% [F], 47.63% [M]), and high density (52.21% [F], 47.79% [M]). 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, 1954, 1958; Lewontin & 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 natural fly populations—mainly various mutant and inversion‐baring strains have been studied, or experiments were performed with a single isofemale line. By utilizing the DGRP, we could test the effects of competition on naturally occurring genetic variation (Mackay et al., 2012) and draw ecologically relevant conclusions about variation for the studied traits. Moreover, studying a panel of genotypes provided us with the opportunity to test for important evolutionary questions, such as how does evolvability of traits change with the increasing level of density stress? Or how pervasive GEIs are for fitness‐related traits?

To test these questions, we performed thorough phenotyping and found for both DT and EAV that the studied genotypes expressed 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 CVP and CVG 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) in studies of interspecies competition between D. melanogaster and D. simulans larvae. However, both results were based on one laboratory strain of D. melanogaster derived from a single wild‐caught female without subsequent generations of dedicated sib–sib mating, as was used to produce the DGRP strains used in the present study. Hence, it is possible that the strain used by Miller harbored substantially more genetic variation than individual DGRP strains. For this reason, Miller's result could partially reflect the increased variance case of GEI, despite the use of a single strain. Nonetheless, it is also possible that some of the variation measured by Miller arose from stochastic differences between genetically near‐identical individuals, a phenomenon known as developmental noise (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 trade‐off between viability and body size that is exerted in the preadult stage in Drosophila under competitive conditions. A potential trade‐off 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.

Conflict of Interest

None declared.

Click here for additional data file.