Introduction

Stressful environmental conditions not only act as a strong selective force, but also directly affect the genetic architecture of traits in several ways. Firstly, stressful conditions can influence the nature of environmental variation, VE (Hoffmann and Merilä 1999; Ørsted et al. 2018a). Secondly, these conditions can directly affect expression of additive genetic variation, VA (e.g., Czesak et al. 2006; Wilson et al. 2006). Both changes in VE and VA will influence heritability estimates (Gebhardt-Henrich and van Noordwijk 1991; Charmantier and Garant 2005; Wilson et al. 2006; Visscher et al. 2008; Dingemanse et al. 2009; Kristensen et al. 2015). Thirdly, environment-dependent gene action may exist, where the function of genes or gene networks, rather than overall genetic variation, change under stressful conditions (Telonis-Scott et al. 2009). This can lead to a situation where genes influencing a trait in one environment may not have much impact on that trait in a different environment (Sgrò and Hoffmann 2004), which in turn will alter genetic correlations between environments (de Jong 1990; Stearns et al. 1991; Vieira et al. 2000; Stinchcombe et al. 2010). The impact of these environmental effects on genetic variation can be important for evolutionary trajectories (Hoffmann and Merilä 1999; Hercus and Hoffmann 2000), but these mechanisms are rarely studied, as they are typically regarded as a nuisance in quantitative genetic models (but see Wilson et al. 2006; Husby et al. 2011; Wood and Brodie 2016).

Environmental conditions are likely to have particularly large effects on the genetic architecture of ecological traits critical for responding to variable environmental conditions. For this class of traits, the environment can influence VA, VE and gene action (Sgrò and Hoffmann 2004; Kristensen et al. 2015), affecting the extent to which organisms respond to selection and hence persist in the long run. Tolerance to stressful temperatures in ectotherms, especially in Drosophila, provides examples of traits whose genetic architecture is strongly affected by environmental variation, such as the impact of developmental temperature on the resistance to cold extremes (Schou et al. 2017; Ørsted et al. 2017).

The ways in which the environment influences trait expression across genotypes is often characterised through genotype-by-environment (GxE) interactions, where the direction and magnitude of GxE can depend on genotype (Rohde et al. 2017; Ørsted et al. 2018a) and/or as genetic correlations of a trait measured in several environments (Vieira et al. 2000; Sgrò and Hoffmann 2004; Agrawal and Stinchcombe 2009). Genotypic responses to environmental variation are usually expressed as differences in phenotypic plasticity, which is commonly quantified with linear reaction norms, i.e., the slope of trait values under continuous variation in an environmental variable (Valladares et al. 2006). However, while environmental effects on gene action have typically been investigated by comparing two environments (Telonis-Scott et al. 2009; Ellers and Driessen 2011; Levine et al. 2011; Ketola et al. 2012; Gerken et al. 2015), patterns across environmental gradients are not well established. If the slopes of different genotypes along a gradient vary, there will be GxE interactions and genetic correlations will be less than one. Under a continuous environmental gradient, a polygenic trait might be expected to show declining genetic correlations (increasing GxE interactions) as environments become more dissimilar, assuming that there is an increase in variance due to environment-specific gene expression (Via and Lande 1985; Sgrò and Blows 2004; Spichtig and Kawecki 2004). This could lead to a linear or non-linear decrease in genetic correlations with the environmental variable. Surprisingly, this rather simple prediction has rarely been studied experimentally. Knowledge of the effects of environmental similarity on genetic correlations of stress resistance traits can prove valuable when trying to predict adaptive responses in heterogeneous environments (Sheldon et al. 2003; Sgrò and Blows 2004), such as those arising from anthropogenic climate change (Etterson and Shaw 2001; Agrawal and Stinchcombe 2009).

Environment-dependent gene action also raises the question of how plasticity is linked to stress tolerance in a particular environment. Phenotypic plasticity can be costly, because it requires energy and flexibility of the organism on a number of biological levels (Auld et al. 2010; Tonsor et al. 2013; Murren et al. 2015), and might become maladaptive if environments change rapidly and in unpredictable ways (Kristensen et al. 2008). Several studies have found an inverse association between stress resistance in one environment (here defined as “basal” tolerance) and the capacity for plasticity (Stillman 2003; Nyamukondiwa et al. 2011; Gerken et al. 2015; Comte and Olden 2017; Noh et al. 2017), suggesting a trade-off between the ability to cope with current conditions and performance in an altered environment in the future (Murren et al. 2015). However, other studies have not found evidence for such trade-offs (Kellett et al. 2005; Calosi et al. 2008; Franke et al. 2012; Gunderson et al. 2015; Ørsted et al. 2018a). Thus, it remains unclear if plasticity constrains basal thermal tolerance and vice versa. Adding to the ambiguity, some studies have found that trade-offs are species-specific even across related taxa (Nyamukondiwa et al. 2011; Strachan et al. 2011; Comte and Olden 2017), while others have found seasonal variation in trade-offs (Noh et al. 2017), highlighting the influence of environmental conditions on these patterns.

Here, we investigate the effect of developmental temperature on the genetic architecture of cold tolerance in Drosophila melanogaster. The aims were to explore the environmental dependency of quantitative genetic parameters of cold tolerance, and potential trade-offs between cold tolerance in each environment and plasticity in this trait. We used the Drosophila melanogaster Genetic Reference Panel (DGRP; Mackay et al. 2012; Huang et al. 2014), which was established from a natural population of D. melanogaster from Raleigh, USA. The DGRP consists of a collection of lines that has been inbred to an expected inbreeding coefficient of F ~1, resulting in extremely high homozygosity throughout the genome within each line. Thus, by assessing cold tolerance of the same DGRP genotype reared in different thermal environments, knowledge about the genetic aspects of cold tolerance, and the association between cold tolerance at individual temperatures and plasticity, can be obtained. Low temperatures limit the geographical distribution of many species more than warm temperatures (Sunday et al. 2011; Williams et al. 2014), and cold tolerance is therefore a good predictor of present and future distributions of species (Kimura 2004; Overgaard et al. 2011; Araújo et al. 2013). The temperature regimes employed in this study were well within the range of what the cosmopolitan D. melanogaster will experience in its natural habitats (Hoffmann et al. 2003), thus providing ecological relevance to our experimental setup. However, it can be difficult to ascertain which of the thermal environments are optimal and which imposes a level of stress, as this is highly dependent on the biological context and traits investigated (David et al. 1997; Hoffmann 2010).

We reared 166 DGRP lines at five different developmental temperatures and quantified cold tolerance in adults from each developmental temperature. We performed quantitative genetic analyses to estimate genetic and environmental variance components, heritabilities, evolvabilities and phenotypic/genetic correlations of cold tolerance in each developmental environment and for (slope-based) plasticity. This allowed three main questions to be answered. (1) To what extent is cold tolerance phenotypically and genetically correlated between developmental environments, and are these correlations dependent on environmental similarity? (2) Do the developmental temperatures impact on variance components, heritability and evolvability of cold tolerance within developmental temperatures, and on plastic cold tolerance? (3) Is there a trade-off between cold tolerance at individual developmental temperatures and the plasticity of cold tolerance, and is this trade-off dependent on the environment?

Materials and methods

Experimental setup

This study included 166 lines from the D. melanogaster Genetic Reference Panel (DGRP) obtained from Bloomington Drosophila Stock Centre (NIH P40OD018537), maintained at 23 °C, 50% relative humidity and a 12:12 h photoperiod on a standard oatmeal–sugar–agar–yeast medium for two generations before starting the experiment. The experimental procedure follows the setup previously described in Ørsted et al. (2018a). In summary, we reared each DGRP line at five different temperatures: 17, 20, 23, 26 and 29 °C, from eggs to adults in vials containing 7 mL standard medium (for details on medium composition see Ørsted et al. 2018b). Eggs were laid by approximately 20 adult flies in three 12 h periods. Vials were checked daily and eclosed flies were sexed under CO2 anaesthesia. Males were transferred to vials containing fresh medium, and kept at the respective developmental temperature for 48 h. At age 60 ± 12 h, the cold tolerance of around 10 males (for exact numbers see Table S1) per DGRP line per temperature were assessed using the dynamic measure Critical Thermal Minimum (CTmin; Overgaard et al. 2012). In this standardised procedure, individual flies were placed in glass vials in a water bath pre-set to 23 °C. Water in the bath was gradually cooled at a rate of 0.1 °C/min. The temperature at which a fly enters chill coma, i.e., when all movement ceases, was recorded and defines CTmin. In total, 7690 flies were scored for CTmin by the same person to reduce potential scoring bias. We chose males because virginity is not likely to impact on males to the same extent as it does to females (Goenaga et al. 2012), and we were not able to assure that all flies were virgins when assessed.

Phenotypic plasticity of cold tolerance was estimated as the regression coefficient of a linear regression, i.e., slope of CTmin on rearing temperature. Linear regression analysis provided the best fit compared to polynomial regression, as assessed by the coefficient of determination (results not shown). In addition, linear norms-of-reaction are the most commonly used in studies of phenotypic plasticity (Valladares et al. 2006), and especially for the relationship between CTmin and developmental temperature, which is typically linear in Drosophila (Schou et al. 2017; Ørsted et al. 2018a).

Quantitative genetic analyses

Genetic analyses was performed with the average information restricted maximum-likelihood (AI-REML) procedure (Madsen et al. 1994) utilising the R package QGG (http://psoerensen.github.io/qgg/) that links to DMU (http://dmu.agrsci.dk/DMU/). We estimated variance components, genetic correlations and broad and narrow sense heritabilities (H2 and h2, respectively) based on multi-trait linear mixed models. For CTmin we fitted a five-trait multi-trait model including all developmental temperatures, and to estimate genetic correlations between CTmin at individual temperatures and plasticity we fitted a bivariate model. The multi-trait model including i-traits was generalised as:

$$\left[ {\begin{array}{*{20}{c}} {{\mathbf{y}}_1} \cr \vdots \cr {{\mathbf{y}}_i} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathbf{X}}_1} \cr \vdots \cr {{\mathbf{X}}_i} \end{array}\begin{array}{*{20}{c}} {{\mathbf{b}}_1} \cr \vdots \cr {{\mathbf{b}}_i} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathbf{Z}}_1} \cr \vdots \cr {{\mathbf{Z}}_i} \end{array}\begin{array}{*{20}{c}} {{\mathbf{g}}_1} \cr \vdots \cr {{\mathbf{g}}_i} \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} {{\mathbf{e}}_1} \cr \vdots \cr {{\mathbf{e}}_i} \end{array}} \right]$$

where, yi was a vector of phenotypic observations (i.e., either temperature specific CTmin or plasticity), Xi and Zi are design matrices linking fixed and random genetic effects to the phenotypes. The bi is a vector of fixed effects (experimental day, water bath, time on day (for CTmin), Wolbachia infection status, and five major chromosomal inversions (only when estimating the additive genetic variance), and ei is a vector of random residuals (\(\boldsymbol {e}_i\sim N\left( {0, \boldsymbol {I}\sigma _{e_i}^2} \right)\)). Estimates of the genetic variance components (and therefore H2) were obtained by assuming the DGRP lines to be independent, which we modelled by an identity block matrix (IL) as (co)variance structure of the genetic values, \(g{\mathrm{\sim }}N\left( {0,{{\boldsymbol {I}}}_L\sigma _g^2} \right)\). To estimate the additive genetic variance components (including h2) we estimated the realised relationship among DGRP lines from publicly available single nucleotide polymorphisms (SNPs). Thus h2 was the proportion of phenotypic variance captured by common segregating SNPs within the DGRP. We used standard filtering for the DGRP, i.e., only considering biallelic SNPs with minor allele frequencies ≥0.05, and for which the Phred scaled variant quality was greater than 500, and the genotype call ≥0.8 (Mackay et al. 2012; Huang et al. 2014). This yielded 1,725,755 genetic markers distributed on the six chromosome arms (2L, 2R, 3L, 3R, 4 and X). The notation is similar as above, however, here we included five major polymorphic inversions present in the DGRP system; I2Lt, I2Rns, I3Rp, I3Rk and I3RMo), gi was a vector of random genetic effects, \({\boldsymbol{g}}_i{\mathrm{\sim }}N\left( {0,{\boldsymbol{G}}\sigma _{g_i}^2} \right)\), where G was the additive genomic relationship matrix, and ei was a random residual, \({\boldsymbol{e}}_i{\mathrm{\sim }}N\left( {0,{\boldsymbol{I}}\sigma _{e_i}^2} \right)\). The G matrix was computed as G = WW′/m (VanRaden 2008), where m was the number of genetic markers, and W was a centred and scaled genotype matrix. Each column vector of W, \({\boldsymbol{w}}_i = \frac{{{\boldsymbol{a}}_i - 2p_i}}{{\sqrt {2p_i\left( {1 - p_i} \right)} }}\), pi was the frequency of the minor allele at locus i, and ai was the ith column vector of the allele count matrix, A, with the genotypes encoded as 0, 2, referring to the number of the minor allele. The SNP data, including karyotype inversions and Wolbachia infection status, can be obtained from the DGRP2 website (http://dgrp2.gnets.ncsu.edu). To achieve a replicated structure of the plasticity measures, which was needed to obtain accurate estimates of the genetic parameters, we grouped the CTmin data within each DGRP line into three groups of individuals according to day of assay, time of day, and water bath. Based on these replicate groups, we obtained three replicates measures of slope per DGRP line. The slope estimates were robust against re-sampling of random individuals into groups of the same size (for details see Ørsted et al. 2018a).

Phenotypic correlations were computed as Pearson’s correlation using within DGRP line means, and standard errors were estimated as: \({\mathrm{SE}}_{r} = \sqrt {\frac{{1 - {r}^2}}{{n - 2}}}\), where n was the sample size. We also computed correlations between the rank order of the DGRP lines (rr), in which lines are ranked within rearing regimes from 1 to 166.

All variance components were tested under the null hypothesis of being zero using Wald’s test. In addition, we tested if correlations differed from zero, i.e., the estimate deviated more than 1.645 × SE from 0 (p < 0.05). We performed pairwise comparisons of genetic correlations, rg among CTmin in the five rearing temperatures (e.g., rg (17, 20 °C) vs. rg (17, 23 °C)), as well as between correlations between CTmin in each environment and plasticity across environments. These comparisons were based on Fisher’s z transformation. Evolvability (IA) was computed at each developmental temperature as: \({I}_{\mathrm{A}} = \widehat {V}_{\mathrm{A}}/\left( {\overline {{\mathrm{CT}}} _{{\mathrm{min}}}} \right)^2 \times 100\), where \(\overline {{\mathrm{CT}}} _{{\mathrm{min}}}\) is the mean CTmin across all DGRP lines in a given environment (Houle 1992; Hoffmann et al. 2016). All analyses were done within the R programming environment (R Core Team 2017).

Results

Heritable variation for cold tolerance at the individual rearing temperatures and plasticity

We observed a significant increase in line mean CTmin, i.e., a decrease in cold tolerance, with increasing rearing temperature (Fig. 1, Table 1, F4,804 = 3956; p < 0.001). Average CTmin (±SE) across DGRP lines spanned from 2.81 ± 0.03 °C for flies reared at 17 to 7.80 ± 0.04 °C for flies reared at 29 °C, with considerable variation among DGRP lines within rearing environment (Fig. 2, Table 1 and Table S1). The estimates of H2 and h2 for CTmin were consistent across the developmental temperatures, with estimates in the range of 0.41–0.48 and 0.27–0.34 for H2 and h2, respectively (Table 2). Plasticity of cold tolerance was estimated as the slope of a linear regression of CTmin as a function of developmental temperature. The slope was significantly different from zero in all DGRP lines (t(3–21) > 2.13; p < 0.001). A significant interaction between line and developmental temperature (two-way ANOVA; F(165,477) = 2.22, p < 0.001) indicated a GxE interaction, reflecting differences among DGRP lines in slope ranging from 0.32 to 0.57 °C in CTmin per 1 °C increase in developmental temperature. The heritability estimates for plasticity of CTmin were 0.65 and 0.51 for H2 and h2, respectively. We observed a higher total genetic variance \(\left( {\widehat {V}_{\mathrm{G}}} \right)\), additive genetic variance \(\left( {\widehat {V}_{\mathrm{A}}} \right)\) and environmental variance \(\left( {\widehat {V}_{\mathrm{E}}} \right)\) of CTmin at 29 °C than at 26 °C (Table 2). At 17, 20 and 23 °C, we observed a decrease in all variance components as compared to 29 °C (Table 2). In contrast, for evolvability (IA; mean corrected \(\widehat {V}_{\mathrm{A}}\)) of CTmin, we observed higher values at lower temperatures, i.e., from 0.67 at 17 to 0.18 at 29 °C (Table 2).

Fig. 1
figure 1

Linear regressions of critical thermal minimum (CTmin) across the five developmental temperatures. Each grey line represents one DGRP line (n = 166). The two full black lines represent the DGRP line with the steepest slope and the DGRP line with the flattest slope, with the difference in CTmin from developing at 17–29 °C indicated as ΔCTmin. The two dashed black lines represent the DGRP lines with the highest and lowest average CTmin across the temperatures indicated as \(\overline {{\mathrm{CT}}} _{\min }\). These latter two DGRP lines represent a genotype that performs poorly across environments and a genotype that performs well across environments

Full size image
Table 1 Mean ± SE, minimum, and maximum CTmin at the five developmental temperatures (17, 20, 23, 26 and 29 °C) and slope of CTmin across the thermal environments for all DGRP lines (n = 166)
Full size table
Fig. 2
figure 2

Correlations of CTmin across developmental thermal environments. Diagonal elements show the DGRP line mean CTmin ( ± SE) sorted by increasing values within rearing regime. Plots above the diagonal show the phenotypic correlations, rp (SE in parenthesis) along with the printed genetic correlations, rg (SE in parenthesis) between CTmin for pairs of environments. Plots below the diagonal show the correlations between the rank order of DGRP lines, rr (SE in parenthesis). In these plots, the DGRP lines are ranked within rearing regimes from 1 to 166. All correlations are Pearson’s product moment correlations. In all plots above and below the diagonal, the grey dots represent one DGRP line. Black lines are the regression lines visualising the correlations. All correlations are significantly different from zero, i.e., the estimate deviate more than 1.645 × SE from 0 (p < 0.05). In the lower right corners of the above diagonal plots, the genetic correlations, rg are compared pairwise within rows, and significantly different correlations are denoted with separate letters, e.g., rg (17, 20 °C) vs. rg (17, 23 °C). In the below diagonal plots, the letters represent significant differences in pairwise comparisons of rank correlations, rr within columns. Some correlations have more letters as they are similar to some but different from other correlations. Comparisons of these correlations are based on Fisher’s z transformation (n = 166)

Full size image
Table 2 Estimated variance components (genetic variance,\(\widehat {V}_{\mathrm{G}}\), environmental variance, \(\widehat {V}_{\mathrm{E}}\), and additive genetic variance, \(\widehat {V}_{\mathrm{A}}\)) and broad and narrow sense heritability estimates (H2 and h2, respectively) ± standard errors (SE), as well as evolvability (IA = \(\widehat {V}_{\mathrm{A}}\)/(mean CTmin)2 × 100) for CTmin at the five developmental temperatures and for the slope of CTmin across temperatures
Full size table

Correlations between CTmin from different developmental temperatures increased with environmental similarity

We found strong positive genetic correlations (rg) between CTmin for DGRP lines at the different rearing temperatures (Fig. 2). The rg values were inversely proportional to the difference between developmental temperatures, i.e., correlations between similar developmental temperatures were larger than between dissimilar developmental temperatures. This pattern was confirmed by statistical pairwise comparisons of rg. The phenotypic correlations (rp) followed the same pattern, decreasing with increasing environmental dissimilarity (Fig. 2). Similarly, the rank order of the DGRP lines was more positively correlated across thermal environments that were more similar to each other (Fig. 2).

Trade-offs between CTmin at individual developmental temperatures and plastic cold tolerance depended on the environment

The phenotypic (rp) correlations between CTmin at the individual temperatures and slope of CTmin across temperatures depended on the environment. At 17 °C, the correlation was significantly negative (Fig. 3), i.e., in this environment genotypes with a high cold tolerance (low CTmin) tended to have a high plasticity of cold tolerance. However, we did not find this association at 20 °C. At higher temperatures, correlations were significantly positive (Fig. 3), thus genotypes with a high cold tolerance (low CTmin) tended to have low plasticity of cold tolerance, suggesting a trade-off between inherent cold tolerance at these developmental temperatures and plasticity. The genetic correlations (rg) followed the same pattern, particularly when involving CTmin scores following development at 17, 26 and 29 °C (Fig. 3).

Fig. 3
figure 3

Correlations between CTmin at the individual temperatures and slope of CTmin across temperatures: a 17 °C, b 20 °C, c 23 °C, d 26 °C and e 29 °C. Each dot represents one DGRP line, and the full black lines are the regression lines visualising the phenotypic correlations. In each environment Pearson’s product moment correlation coefficients are shown for the phenotypic correlations, rp (SE in parenthesis) along with the printed genetic correlations, rg (SE in parenthesis) between CTmin. Asterisks designate correlations that are significantly different from 0, i.e., the estimate deviate more than 1.645 × SE from 0 (p < 0.05). Correlations that are not different from 0 are designated “ns”. In pairwise comparisons of the correlations, significantly different correlations are denoted with separate letters after the correlation coefficients, e.g., rg (CTmin 17 °C, slope) vs. rg (CTmin 20 °C, slope). Phenotypic and genetic correlations are compared separately. Comparisons of these correlations are based on Fisher’s z transformation (n = 166). The two arrows in panel a shows the direction of increased fitness (here opposite of increased CTmin) and increased plasticity of CTmin, respectively, for all panels

Full size image

Discussion

In the present study, we investigated the effects of developmental temperature on quantitative genetic parameters of cold tolerance, which is a key stress tolerance trait. We showed that environmental conditions can directly influence the quantitative genetic basis of cold tolerance by affecting variance components, especially \(\widehat {V}_{\mathrm{A}}\). It is expected that if standing genetic variation is environmentally dependent, selection responses will also be affected (Roff 2002; Bijlsma and Loeschcke 2005; Wilson et al. 2006; Husby et al. 2011; Wood and Brodie 2016). This has been demonstrated in recent studies showing that genetic variation and evolutionary potential of stress resistance, morphological traits and life-history traits depend on the environmental conditions (van Heerwaarden and Sgrò 2011; Bubliy et al. 2012; Kristensen et al. 2015; Bastide et al. 2016; Zhao et al. 2016; Feldman et al. 2017). Such findings can be explained by environment-specific gene expression, i.e., the impact of genes on a trait can depend on prior environmental exposures. This highlights the need to perform experiments and field studies under environmental conditions that are relevant to organisms at the time traits are under selection (Hoffmann and Parsons 1991; Hoffmann and Merilä 1999). This is important not only from an ecological and evolutionary perspective (Bay et al. 2017) but also for animal and plant breeding, in which effects of the environment on the genetic basis of production traits pose a challenge, such as for genomic selection in heterogeneous environments (Kadarmideen et al. 2006; Zhao et al. 2016).

It has long been debated whether some types of environmental conditions, and particularly stressful or benign conditions, affect heritable variation in distinct ways (Hoffmann and Parsons 1991; Møller and Swaddle 1998). Here we showed an increase in \(\widehat {V}_{\mathrm{A}}\) of cold resistance with increasing developmental temperature, and heritability estimates being similar across developmental temperatures (Table 2). However, for mean corrected \(\widehat {V}_{\mathrm{A}}\), (i.e., evolvability (IA)), we found higher values at 17 and 20 °C than in the range 23–29 °C, suggesting that the evolutionary capacity for increasing cold tolerance might be more restricted in warm than in cold environments. This is supported by recent Drosophila studies showing that the potential for evolution is lower in relatively warmer conditions (Schou et al. 2014; Kristensen et al. 2015). In nature, individuals reared under cold conditions are more likely to encounter subsequent extreme cold conditions, therefore stronger selection for cold tolerance might be expected in natural populations experiencing low developmental temperatures compared to in hotter environments (Ayrinhac et al. 2004). Thus, in the case of cold tolerance, developmental environments with low temperatures might be considered more optimal than those with high temperatures. However, because we did not have a direct measure of stress in each thermal environment, we cannot confidently determine whether the patterns in \(\widehat {V}_{\mathrm{A}}\) and IA are associated with unfavourable environmental conditions (Hoffmann and Merilä 1999; Swindell and Bouzat 2006; Talloen et al. 2009).

The environment can also affect genetic variation via environment-dependent gene action, where genes influencing a certain trait in one environment may not be important for that trait in a different environment. In such cases selection responses can be slowed or even become maladaptive (Steinger et al. 2003). This will often be expressed as altered genetic correlations either between different traits in one environment or between the same trait in a range of environments (Vieira et al. 2000; Sgrò and Hoffmann 2004; Agrawal and Stinchcombe 2009). Here we found positive genetic and phenotypic correlations between cold tolerance across environments, i.e., genotypes with high cold tolerance when developed at one temperature also tended to have high cold tolerance when developed at other temperatures, implying an overlap in the sets of alleles influencing cold tolerance across developmental thermal environments (Via and Lande 1985; Falconer and Mackay 1996; Agrawal and Stinchcombe 2009). However, we also showed that the magnitude of correlations between thermal environments decreased as developmental temperatures diverged. This suggests a shared genetic basis for cold tolerance between rearing environments, which contributes less to genetic variation as environments become more dissimilar. Although evidence for such patterns has rarely been investigated, a few studies have found similar results; Sgrò and Blows (2004) found declining genetic correlations with environmental dissimilarity in three life-history traits in Drosophila serrata, while Stinchcombe et al. (2010) found declining genetic correlations with environmental distance for growth rate in jewelweed.

Our findings are in line with other studies suggesting that the genetic architecture of fitness components is specific to an environment (Bourret and Garant 2015; Rellstab et al. 2017). This represents a challenge when testing for adaptation at specific loci which often seem to be environment and population specific (Zhao et al. 2011; Manenti et al. 2016). Knowledge about the genetic architecture of many complex traits has been rapidly accumulating recently, aided by developments in genomics and quantitative genetics. Typically, such studies identify genes, variance components and heritabilities for a given trait in one environment to predict evolutionary trajectories. In animal and plant breeding, such knowledge may be used to predict selection responses and devise genomic selection procedures. However, our results show that the architecture of an important stress resistance trait is highly dependent on environmental conditions, challenging the notion that there are candidate genes and a particular architecture applicable across a gradient of environmental conditions.

An important issue for the evolution of thermal resistance is whether an organism’s thermotolerance is constrained by plasticity, particularly in response to anthropogenic climate change (Stillman 2003; Chown et al. 2010; Levine et al. 2011; Gunderson et al. 2015; Comte and Olden 2017). Some studies point to a negative association between an organism’s thermal tolerance and plasticity in tolerance (Hoffmann et al. 2003; Kellett et al. 2005; Nyamukondiwa et al. 2011; Gerken et al. 2015; Noh et al. 2017). Here we showed that any association between cold tolerance and plasticity depend on the environment; both the phenotypic and genetic correlations shifted in sign as developmental temperatures shifted from 17 to 29 °C (Fig. 3). Given that the lowest CTmin occurred at the lowest developmental temperature, genotypes with a flat slope would likely have a fitness advantage when exposed to low temperatures, as they would remain the most cold tolerant (i.e., maintain relatively lower CTmin at higher developmental temperatures). However, as indicated by the negative correlation between slope and CTmin at 17 °C, DGRP lines with low CTmin reared at this temperature, typically had a high CTmin (low cold tolerance) when they developed at higher temperatures, indicative of a trade-off (Fig. 3a). At low developmental temperatures, lines with low CTmin and a steep slope could be considered “specialists”; i.e., they performed relatively better in one environment, but relatively poorly across environments. At higher developmental temperatures (23, 26 and 29 °C), we found a positive correlation between CTmin and the plastic response, suggesting that high cold tolerance (low CTmin) at high developmental temperatures was associated with low plasticity (Fig. 3c–e). Thus, genotypes with low CTmin at higher developmental temperatures could also be considered as the “specialists” because they performed relatively better in one environment. This shift in sign of the correlations reflects a reversal in the ordering of genotypes for cold resistance across the thermal gradient, with a tipping point at 20 °C, when there was no correlation between cold tolerance and plasticity.

These genotypic differences in constraints may reflect differential environmental sensitivity and suggests that the fitness effect of plasticity depends heavily on the environment where genotypes are reared. Thus, the fitness landscape of cold tolerance needs to be considered in the context of environmental frequencies, i.e., how likely are individuals reared under one set of thermal conditions to encounter subsequent extreme cold conditions (Ayrinhac et al. 2004)? Our findings support those of Levine et al. (2011) who suggested a “directionality” of plasticity in gene expression depending on the thermal environment commonly encountered in the native range of populations of D. melanogaster. This has implications for understanding developmental processes in variable environments (Stearns et al. 1991), and also for understanding local adaptation in thermal tolerances in populations with a broad geographic distribution along an environmental cline (Ayrinhac et al. 2004; Levine et al. 2011; Cooper et al. 2012; Kellermann et al. 2012; Kristensen et al. 2015). Our results indicate that the quantitative genetic basis of basal cold tolerance and plastic cold tolerance is highly environment specific, which is important for predicting selection responses in natural and variable environments (Sgrò and Hoffmann 2004; Agrawal and Stinchcombe 2009).

Most studies investigating phenotypic and genetic correlations and/or candidate genes between thermal environments focus on only two temperatures (Guerra et al. 1997; Ellers and Driessen 2011; Ketola et al. 2012) or investigate two settings such as field vs. lab conditions (Kristensen et al. 2008; Overgaard et al. 2010), and thus do not provide a comprehensive picture of changing genetic correlations and architectures across environments. A meta-analysis by Charmantier and Garant (2005) specifically called for approaches that consider genetic variability under multiple conditions, and ideally along a continuous environmental gradient such as in the present study. Because correlations within and between traits can change between environments or even shift in sign (Sgrò and Hoffmann 2004; Agrawal and Stinchcombe 2009) as evident in this study, we reiterate this call and caution about genetic constraints and trade-offs deduced from studies that consider only one or two environments.

In conclusion, we showed that the expression of additive genetic variation for cold tolerance was environment dependent, and that evolvability of cold resistance decreased with increasing developmental temperatures. We also showed that the relationship between basal and plastic cold tolerance was affected by the environment, and even switched sign across a gradient of developmental temperatures. We further provided evidence for an evolutionary trade-off between cold tolerance in one environment and cold tolerance across environment. It will be interesting to test for these patterns in relevant ecological settings, such as in populations exposed to strong diurnal and seasonal temperature fluctuations.