Evolution of thermal tolerance in multifarious environments

Species extinction rates are many times greater than the direst predictions made two decades ago by environmentalists, largely because of human impact. Major concerns are associated with the predicted higher recurrence and severity of extreme events, such as heat waves. Although tolerance to these extreme events is instrumental to species survival, little is known whether and how it evolves in natural populations, and to what extent it is affected by other environmental stressors. Here, we study physiological and molecular mechanisms of thermal tolerance over evolutionary times in multifarious environments. Using the practice of “resurrection ecology” on the keystone grazer Daphnia magna, we quantified genetic and plastic differences in physiological and molecular traits linked to thermal tolerance in historical and modern genotypes of the same population. This population experienced an increase in average temperature and occurrence of heat waves, in addition to dramatic changes in water chemistry, over five decades. On genotypes resurrected across the five decades, we measured plastic and genetic differences in CTmax, body size, Hb content and differential expression of four heat shock proteins after exposure to temperature as single stress and in combination with food levels and insecticide loads. We observed evolution of the critical thermal maximum and plastic response in body size, HSP expression and Hb content over time in a warming only scenario. Molecular and physiological responses to extreme temperature in multifarious environments were not predictable from the response to warming alone. Underestimating the effect of multiple stressors on thermal tolerance can lead to wrong estimates of species evolvability and persistence.

and North America due to atmospheric circulation patterns intensified by greenhouse gases (Meehl & Tebaldi, 2004), and to become more intense, longer lasting and/or more frequent (Karl & Trenberth, 2003).
Physiological and/or molecular plasticity has been shown to play an important role in organismal response to extreme temperatures.
Here, we study molecular and physiological mechanisms of thermal tolerance in an invertebrate ectotherm common to European lotic freshwater ecosystems, Daphnia magna (Miner, De Meester, Pfrender, Lampert, & Hairston, 2012). Daphnia spp are keystone species in lakes and ponds worldwide, and drivers of ecosystem dynamics (Altshuler et al., 2011;Miner et al., 2012). In the natural environment, D. magna is exposed to severe spatial and temporal environmental changes, including temperature, low oxygen and eutrophication, to which it responds via an ecoresponsive genome (Colbourne et al., 2011;Orsini, Gilbert et al., 2016Orsini et al., 2018. Daphnia magna has a parthenogenetic life cycle that allows the rearing of populations of genetically identical individuals (clones) from a single genotype. In unfavourable environmental conditions, D. magna switches from asexual to sexual reproduction, producing dormant embryos that arrest their development entering dormancy (Ebert, 2005). These dormant embryos sink into lake sediment and build up archives of living fossils that remain viable for decades or even centuries (Cousyn et al., 2001;Decaestecker et al., 2007;Frisch et al., 2014;Orsini, Spanier, & De Meester, 2012;Orsini, Marshall, et al., 2016). With the practice of "resurrection ecology" (Kerfoot & Weider, 2004), the dormant embryos can be resuscitated and historical and modern populations used in the same experimental settings.
These properties provide us with the unique opportunity to study mechanisms of response to environmental changes through evolutionary time and to disentangle the relative contribution of plastic and genetic response to environmental change.
Previous studies have shown adaptive response of D. magna to temperature changes via physiological and molecular mechanisms (Cambronero, Zeis, & Orsini, 2017;Geerts et al., 2015;Jansen et al., 2017). Geerts and co workers provided evidence of evolution of the temperature of maximum tolerance (CT max ) across few decades in response to warmer climates (Geerts et al., 2015). In a follow-up study, Jansen et al. showed that evolution in CT max was mediated by both plastic and evolutionary changes in gene expression at a number of candidate genes, including heat shock proteins . Cambronero et al. provided (Cambronero et al., 2017).
Three populations of D. magna separated in time were previously sampled from a well-characterized lake in Denmark, Lake Ring, which experienced an increase in average temperature and occurrence of heat waves across five decades  (Cuenca Cambronero & Orsini, 2018;Orsini, Marshall, et al., 2016). The lake also experienced changes in water chemistry transitioning from a severe event of eutrophication (1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970) to a partial recovery in modern times (>1999); the lake also suffered from increased agricultural land use between 1975 and 1985 (Cambronero et al., 2017;Orsini, Marshall, et al., 2016). In a previous study, Cambronero et al. (in review) used these three populations in common garden experiments to measure genetic and plastic response in life history traits to temperature as single stress and in combination with insecticide loads and food levels, mimicking major stressors that occurred in the lake (Cambronero et al., in review). Here, we assessed the impact of multifarious environments on thermal tolerance by measuring the critical thermal maximum (CT max ) and by quantifying the differential expression of four heat shock proteins (HSP20, HSP60, HSP70 and HSP90) at the upper thermal limit on the genotypes used in the common garden experiments. Furthermore, we reanalysed previous data generated on Hb content and body size collected on the same genotypes (Cambronero et al., 2017, in review).
We expected the physiological and molecular traits to evolve over time in response to an increase in average ambient temperature and occurrence of heat waves. Our results show that coping mechanisms to thermal stress over evolutionary time in the presence of multiple environmental stressors are not predictable from the response to temperature as a single stressor.

| Sampling site and Lake Ring paleolimnological profile
The study site is Lake Ring, a well-characterized lake in Jutland, Denmark (55°57′51.83″ N, 9°35′46.87″ E). In 2004, a sedimentary archive was sampled from the lake and stored in dark and cold (4°C) conditions. In 2015, a radiometric chronology of this archive was completed by ENSIS Ltd (UCL London) following standard protocols (Appleby, 2001). Environmental changes in the Lake were reconstructed from the paleolimnological analysis of sediment and historical records as described in Ref. (Cambronero et al., 2017;Cuenca Cambronero & Orsini, 2018). The paleolimnological analysis consisted of quantifying the organic matter in the sediment, measured as loss on ignition (LOI) (Heiri, Lotter, & Lemcke, 2001), and the invertebrate community (Cladocera) assemblage over time. Furthermore, total phosphorus, nitrogen and water transparency were measured in Lake Ring by the County of Vejle in the Jutland peninsula for the years 1971-1999 as part of a monitoring programme following standard protocols (Søndergaard et al., 1990). Temperature records were collected by the Danish Meteorological Institute at a weather station located 80 km from the Lake over the past century (http://www.dmi.dk/laer-om/ generelt/dmi-publikationer/2013/). Because air and water surface temperature have a positive correlation for shallow streams and lakes (Preudhomme & Stefan, 1992), especially for the summer months (e.g., (Livingstone & Lotter, 1998), the data from the weather station were previously used as estimates of the monthly water temperature in the lake (Cambronero et al., in review).
Historical records of pesticides used in Denmark were collected from the Danish national archives (www.middeldatabasen.dk).
According to these archives, carbamate pesticides (e.g., Carbaryl) were commonly used in the 1980s. The paleolimnological and the historical data revealed that temperature, high primary production and Carbaryl levels were key environmental factors impacting Lake Ring in the five decades analysed (Cambronero et al., in review;Michels, 2007). Based on the occurrence of these stressors over time, three main phases were identified as follows: (a) eutrophication phase (EP, 1960(EP, -1970 characterized by high primary production (high LOI) due to sewage inflow; (b) pesticides phase (PP, 1975(PP, -1985 driven by increase in agricultural land use. This phase is associated with heavy usage of carbamate insecticides among others; and (c) clear water or recovery phase (CWP, >1999), associated with a decrease in primary production, phosphorus and nitrogen, and increase in water transparency. This phase coincided with the diversion of sewage inflow and reduction in agricultural land use in modern times.
From each lake phase, Cambronero et al. (in review) resurrected D. magna populations separated in time following established protocols (Cuenca Cambronero & Orsini, 2018). Among the total hatched embryos from the sediment core (N = 262), 10 random genotypes from each lake phase (N = 30) were used in common garden experiments to assess genetic and plastic response in life history traits (size at maturity, age at maturity, fecundity and mortality) to temperature as single stress and in combination with loads of the insecticide Carbaryl or food levels. Prior to starting the CGEs, the genotypes were acclimated and synchronized for two generations in common garden conditions (16:8 light: dark regime, 16 ± 1°C and fed ad libitum 0.8 mg Carbon/L of C. vulgaris daily) to reduce interference from maternal effect. After two generations in these conditions, 24-48 hr old juveniles from the second or following broods were randomly assigned to the experimental conditions in which life history traits were measured in the time spanning an individual life cycle (until release of the second brood). In CGE1, the populations were exposed to 24 ± 1°C and 18 ± 1°C; the experimental animals were fed ad libitum with Chlorella vulgaris (0.8 mg C/L). In CGE2, the two experimental temperatures (18 ± 1°C and 24 ± 1°C) were crossed with two nutrient levels: 0.2 mg C/L and 2.4 mg C/L. In CGE3, the two experi-

| Physiological and molecular response to thermal stress
To assess the impact of multifarious environments on D. magna thermal tolerance, we studied physiological and molecular responses of the 30 genotypes resurrected from Lake Ring after exposure to temperature as single stress and in combination with either biotic stress (two food levels) or abiotic stress (two concentrations of the insecticide Carbaryl) in common garden experiments. At completion of the common garden experiments described above, we measured the critical thermal maximum (CT max ) on the 30 genotypes across all exposures. Furthermore, we quantified the expression of four heat shock proteins (HSP20, HSP60, HSP70 and HSP90) at the upper thermal limit. Differential expression was measured between genotypes at their upper thermal limit and the same genotypes reared in non stress conditions. HSP expression was measured on a subset of the genotypes from the common garden experiments exposed to temperature as a single stress and in combination with food limitation or one insecticide concentration. These specific conditions were previously identified as imposing severe sub lethal effects on life history traits (Cambronero et al., in review). Methods for measuring CT max and HSP expression are described below. In addition to data on CT max and HSP expression, we reanalysed data on body size from Cambronero et al. (in review) and data on haemoglobin induction available on the same genotypes (Cambronero et al., 2017).  (Table S1). These included exposure to temperature as single stress and in combination with biotic or abiotic stress.

| Critical thermal maximum
Following the protocol optimized by Geerts et al. (2015), we exposed individual adult Daphnia after they released their second brood to temperature ramping of 1°C increments every 20 ± 5 seconds until the animals stopped swimming (Geerts et al., 2015). Ten genotypes (biological replicates) per population and treatment, in a single clonal replicate, were used. The temperature ramping essay reflects a rapid increase in temperature associated with heat waves. In shallow lakes and ponds, air and water temperature have a positive correlation (Livingstone & Lotter, 1998;Schneider & Hook, 2010); hence, the simulated rapid increase in temperature is a realistic approximation of severe heat waves occurring in the natural environment. The individuals exposed to temperature ramping were placed in Eppendorf tubes (1.5 ml) within their own medium. The temperature ramping, starting at 16°C, was performed in a thermal block (Eppendorf ThermoMixer C) where the swimming activity was continuously monitored by an operator, who recorded the temperature of ceased swimming activity using a voice recording device. The temperature was monitored both on the display of the thermal block and with a thermometer. A second operator assisted by flash freezing animals in liquid nitrogen immediately after they stopped swimming. On these animals, we quantified the expression of four heat shock proteins using as reference the same genotypes not exposed to temperature ramping.

| Heat shock proteins
Differential expression (DE) of four heat shock proteins [HSP20, HSP60, HSP70 and HSP90, ] was quantified in a subset of genotypes from the common garden experiments. HSPs differential expression was measured in three genotypes (biological replicates) per population, in three technical replicates. These genotypes were from the common garden experiments assessing the impact of temperature as single stress and in combination with food limitation (0.2 mg C/L) and one concentration of Carbaryl (4μgL) (N = 1,440, Table S2).
Total RNA was extracted from single individuals using Agencourt ® RNAdvance™ Tissue kit, following the manufacturer's instructions. RNA was quantified using Nanodrop ® Technologies (USA) and cDNA synthesis was performed using the AffinityScript cDNA synthesis kit, following the manufacturer's instructions. qPCR was performed using a Roche LC96 qPCR machine, following the SYBR ® Premix Ex Taq (Tli, RNaseH Plus, Takara) protocol in 20 μL final volume. The cycling was as follows: initial denaturation for 1 min at 95°C, followed by 40 cycles consisting of 5″ at 95°C, 30″ at 60°C and 1″ minute at 72°C. A final extension of 5 min at 72°C was used. We calculated the mean CT (cycle threshold) value per sample averaging among replicates before rescaling the value. CT values per sample and per protein were rescaled using an interplate calibrator consisting of a pool of RNA extracted from different genotypes at different developmental stages and including genotypes from the three populations studied here following .

| Mechanisms of heat tolerance
We performed an analysis of variance on the physiological (CT max and body size) and molecular traits (HSP and Hb) to assess evolutionary mechanisms-plastic, genetic or a combination thereof-underpinning thermal tolerance. Haemoglobin content of each genotype from the three populations separated in time was quantified in a previous study at 20 ± 1°C and 30 ± 1°C (Cambronero et al., 2017).
The body size of adult Daphnia was quantified in the common garden experiments described above after exposure to temperature as single stressors and in combination with biotic and abiotic stress (Cambronero et al., in review) Results of the variance analysis are interpreted as follows: a) a significant population term indicates genetic differences among populations in the molecular or physiological trait; b) a significant response to treatment(s) indicates plasticity in the trait; c) a significant interaction between the population and the treatment(s) indicates evolution of plasticity in the trait. We analysed physiological traits via ANOVAs using linear mixed models (LMMs) in R v.3.3.3 (Team RC, 2017) and including genotypes as a random effect. We visualized the main effects of population and treatment(s) as well as population × treatment interactions on CT max and body size through univariate reaction norms, which describe the pattern of phenotypic expression of each genotype across treatments (Roff, 1997).
We analysed the expression of HSP proteins at the upper thermal limit as compared to non stressful conditions via multivariate statistics (MANOVA) followed by a univariate analysis per single protein (ANOVA). Both analyses were performed using linear mixed models (LMMs) in R v.3.3.3 (Team 2017). For these analyses, we included a random error structure in each model to account for genotypespecific differences in gene expression within each population. We fitted one model per gene. In this analysis, the term "population" represents the constitutive differences in gene expression among populations. The term "temperature ramping (TR)" assessed whether the single proteins (ANOVA) or the four proteins at once (MANOVA) displayed a significant change in expression in response to the treatment. This term measures the short-term plastic response of gene expression when genotypes are exposed to sudden and severe temperature increase. The interaction term "TR" × "population" quantifies the difference in gene expression among populations due to the treatment, and it reflects the evolution of plasticity in the expression of individual proteins or all HSPs considered at once.
We visualized the multivariate analysis results using phenotypic trajectories (PTA) (Collyer & Adams, 2007). For each population, we plotted (a) the magnitude of phenotypic change across the HSPs between TR treated and non treated populations; (b) and the direction of change in the phenotypic trajectory (Collyer & Adams, 2007;Dennis, Carter, Hentley, & Beckerman, 2011). Significant differences (p < 0.05) in magnitude and direction of phenotypic change between population pairs were derived from a residual randomization procedure following (Collyer & Adams, 2007). Scripts used for this analysis are in Ref. (Adams & Collyer, 2009). Individual HSPs differential expression was visualized through univariate reaction norms.
Population differential expression in Hb and supporting statistics are in Ref. (Cambronero et al., 2017). Here, we visualize these data as univariate reaction norms.

| CGE1 (temperature)
The effect of temperature treatment on CT max was significantly different among populations, whereas body size did not vary significantly among populations (Table 1-CGE1, Figure 1). Temperature treatment induced significant plastic response in CT max and body size in the three populations. Neither CT max nor body size showed significant evolution of plasticity over the five decades examined (Table 1-CGE1).
The univariate reaction norms, supporting the ANOVA analysis, showed that CT max was constitutively higher in the most recent populations (CWP) than in the other two populations; the three populations showed a higher CT max at 24°C than at 18°C (Figure 1-CGE1).
Body size was smaller at 24°C than at 18°C in all populations, and the plastic response was comparable across the three populations ( Figure 1-CGE1).

| CGE2 (temperature and food)
The effect of temperature combined with food levels varied significantly among populations for body size but not for CT max (Table 1-CGE2). We observed a significant plastic response of both traits to temperature and to food levels as single stressors (Table 1-CGE2, T and F). However, a significant interaction between temperature and food (T × F) was observed only for CT max , indicative of additive effect of these two factors (Table 1-CGE2). Evolution of plasticity was observed in CT max with significant interaction terms "temperature × population" and "food × population" (Table 2-CGE2; P × T and P × F).
The univariate reaction norms showed lower CT max in low food levels than in high food levels and a different plastic response of the populations within food levels (Figure 1-CGE2). In high food levels, TA B L E 1 Analysis of variance ANOVA for CT max and body size on the three populations of Daphnia magna after exposure in common garden experiments to temperature (CGE1), temperature and food levels (CGE2) and temperature and insecticide loads (CGE3). The population (P) term indicates genetic differences among populations in CT max and body size; the terms temperature (T), food (F) and insecticide (I) indicate plastic responses of CT max and body size to these treatment(s); the interaction terms indicate evolution of plasticity.

| CGE3 (temperature and pesticides)
The effect of temperature combined with the insecticide Carbaryl did not differ among populations in both physiological traits ( The reaction norms showed comparable plastic responses of CT max at both insecticide levels in the three populations (Figure 1-CGE3). The temperature of maximum tolerance was generally higher at 24°C than at 18°C, except for the PP population in low Carbaryl and the EP population in high Carbaryl (Figure 3-CGE3). Plastic responses in body size differed between high and low Carbaryl (Figure 1-CGE3). Specifically, body size was larger at 24°C than at 18°C in the presence of low Carbaryl in all three populations, whereas plastic response in body size differed among populations in high Carbaryl (Figure 1-CGE3, solid lines). In the presence of high Carbaryl, the PP population went extinct at 18°C. Therefore, we were unable to visualize the reaction norm for this population.

| Molecular response: HSPs and Hb
In the following, we present results on the HSPs multivariate analysis. In this analysis, the differential expression of the four HSPs is shown at 18°C and 24°C, and measured between genotypes exposed to the temperature ramping treatment (increase of 1°C every 20 ± 5 until swimming stops) and the same genotypes not exposed to this treatment (control).

| Warming
The effect of temperature ramping treatment (TR) on the HSPs did not vary significantly among populations both at 18°C and 24°C (   Table S3).
The univariate statistics (ANOVA) confirmed that HSPs differential expression did not vary significantly among populations at 18°C, except for HSP20, and that the expression of HSP60, HSP70 and HSP90 differed significantly among populations at 24°C F I G U R E 1 Reaction norms of physiological responses Population reaction norms, based on population means and SD (n = 10), showing CT max and body size changes after exposure in common garden experiments to temperature (CGE1), temperature combined with food levels (CGE2) and temperature combined with insecticide Carbaryl levels (CGE3). The populations are colour coded as follows: eutrophic population (EP, 1960(EP, -1970 in blue; pesticides population (PP, 1975(PP, -1985 (Figure 3-Warming, Table S4). Significant plastic response of the HSPs to TR treatment observed in the MANOVA was confirmed by the ANOVAs on individual proteins. We observed significant upregulation of HSP20, HSP60 and HSP70 at 18°C, and of HSP60, HSP70 and HSP90 at 24°C (Figure 3-Warming, Table S4). The univariate analysis also identified significant interaction between "TR" and "population" in HSP70 at 24°C (Figure 3-Warming, Table S4).

| Warming and food limitation
The effect of TR treatment on the HSPs varied significantly among populations at 18°C, but not at 24°C (   Table S3).
The univariate statistics confirmed a significant constitutive dif-  Table S4). A significant interaction term between "TR" and "population" was observed for all HSPs at 18°C and for HSP70 at 24°C (Figure 3-Warming and food limitation; Table S4).

| Warming and insecticide
The effect of TR treatment on the HSPs did not vary significantly among populations at both 18°C and 24°C (   Table S3). At 24°C, there was no significant difference in the magnitude of phenotypic trajectories TA B L E 2 Multivariate analysis of variance MANOVA showing multivariate response of four heat shock proteins to the temperature ramping treatment (TR). The TR treatment was performed after exposure of the genotypes to warming, temperature combined with food limitation (Warming and food), and temperature combined with one concentration of the insecticide Carbaryl (Warming and insecticide) in common garden experiments. The population (P) term indicates genetic differences in the overall HSP expression among populations; the TR term indicates a plastic response of the HSPs to temperature ramping treatment; the interaction term between P and TR indicates evolution of plasticity.  Table S3).
The univariate statistics confirmed that the constitutive expres-  Table   S4). Significant plastic responses identified by the MANOVA were confirmed by the univariate statistics (Table S4). We observed significant upregulation of HSP20, HSP60 and HSP70 at 18°C, and significant differential expression of HSP70 and HSP90 at 24°C (Figure 3-Warming and insecticide; Table S4). A significant interaction term, "TR" × "population," was observed in HSP20 and HSP60 at 18°C, and in HSP60 and HSP70 at 24°C (Figure 3, Table S4).

| Haemoglobin
Previously, Cambronero et al. (2017) showed that the constitutive Hb protein content did not significantly vary among populations at 20° C and 30°C, and that the synthesis of haemoglobin was higher at 30°C (Cambronero et al., 2017). These results are visualized here as population reaction norms (Figure 4).

| D ISCUSS I ON
Although tolerance to extreme temperatures is critical to survival and organisms' fitness, little progress has been made to elucidate whether it evolves in natural populations and by what mechanisms.
Many studies investigating heat tolerance tend to minimize confounding factors, searching for trends in relatively undisturbed systems (Anderson, Inouye, McKinney, Colautti, & Mitchell-Olds, 2012;Geerts et al., 2015;) but see (Brans et al., 2017). Limiting confounding factors has obvious advantages but can lead to wrong estimates of evolvability in the wild.
We studied the evolution of thermal tolerance in populations of the waterflea D. magna resurrected from different time points of the same biological archive. These populations had been exposed to an average increase in ambient temperature and occurrence of heat waves, as well as to changes in water chemistry, over five decades.
Capitalizing on common garden experiments previously conducted on the resurrected populations, we assessed the impact of temperature as single stressor and in combination with biotic (food) and abiotic (insecticide) stress on physiological and molecular mechanisms of thermal tolerance.
In experimental conditions mimicking warming as single stressor, the three populations separated in time showed a plastic response in CT max , with a positive correlation between upper thermal limit and plasticity in response to high temperature. Moreover, the most recent population in Lake Ring (CWP) showed a constitutive higher CT max than the two historical populations. Therefore, the most recent population showed both a constitutively higher thermal maximum and comparable plasticity to the other populations. This result suggests a positive correlation between constitutive and induced thermal tolerance, at least in the presence of warming as single stressor. This evolutionary response likely occurred in response to increase in average ambient temperature and occurrence of heath F I G U R E 2 Phenotypic trajectory analysis PTA on the three populations of Daphnia magna resurrected from Lake Ring, resulting from multivariate response of four heat shock proteins to temperature ramping treatment at 18°C and 24°C. Patterns for PC1 and PC2 are shown for three combinations of stressors to which the populations were exposed prior to the TR treatment: (a) warming; (b) warming and food limitation; (c) warming and insecticide. Full circles represent the control (animals not exposed to temperature ramping) and open circles represent the temperature ramping (TR) treatment. HSP expression is shown for each genotype and its replicates. Population centroids are connected by reaction norms (solid lines). Differences among populations, in terms of magnitude (M) and direction (θ) of plastic response, are shown for significant pairwise population comparisons. ns: non significant. Populations are colour coded as in Fig. 1: EP-blue; PP-green; CWP-red. The statistics supporting the PTA are in Table S3 - waves in the five decades studied here, recorded by a weather station adjacent to Lake Ring and supported by climate records in Europe (Committee on Climate Change 2017). In a previous study, the evolution of the critical thermal maximum over few decades was associated with the capacity of D. magna to tolerate higher temperatures (Geerts et al., 2015). Our results corroborate these previous findings and provide a further line of evidence supporting the evolution of CT max in D. magna, at least in a warming only scenario.
The evolution of CT max over time in a warming only scenario was not reflected in the evolution of higher Hb synthesis and HSP expression. Because of the direct link between temperature changes and oxygen solubility, changes in Hb content have been studied in association with thermal stress in ectotherms (Lamkemeyer et al., 2003;Paul et al., 2004;Pörtner, 2002;Verberk et al., 2016). Previous studies have shown that higher Hb content enables ectotherms to cope with hyper thermal stress after acclimation to high temperature (Cambronero et al., 2017;Lamkemeyer et al., 2003;Paul et al., 2004;Pörtner & Knust, 2007;Verberk et al., 2016), whereas the ability to cope with this stress is dampened in absence of temperature acclimation (Cambronero et al., 2017;Seidl, Pirow, & Paul, 2005;Williams, Dick, & Yampolsky, 2012). These findings, supported by the evidence presented here, suggest that haemoglobin expression and CT max evolve on different time scales and/or under different constraints.
F I G U R E 3 Differential expression of individual HSPs Population average (n = 3 genotypes × three technical replicates) and SD of the differential expression of four heat shock proteins after temperature ramping (TR) treatment at 18°C and 24°C. Expression of HSPs is measured on the same genotypes in absence of TR and after TR treatment. The TR treatment was imposed after exposure to (1) warming; (2) warming and food limitation; (3) warming and insecticide. Dotted lines indicate significant constitutive differences in gene expression among populations (population term in Table S4). Full circles indicate a significant plastic response in gene expression induced by TR treatment (TR term in Table S4). Asterisks (*) indicate significant interaction terms "population" × "TR" (P × TR in Table S4). Populations are colour coded as in Figure   F I G U R E 4 Haemoglobin differential expression population reaction norms based on population means (n = 10) and SD are shown for overall haemoglobin content measured at control temperature (20 ± 1°C) and after exposure to hyper thermal stress (30 ± 1°C). Populations are colour coded as in Figure 1: EP-blue; PP-green; CWP-red. Statistics supporting this plot and the Hb data used to generate these plots are from (Cambronero et al., 2017)  In natural ecosystems and especially in enclosed habitats (e.g., lakes and ponds), factors such as eutrophication or by-products of land use may play an important role in driving evolutionary responses by altering solubility of nutrient, conductivity and oxygen levels (Feuchtmayr et al., 2009). These factors may interact additively or synergistically with temperature. There is increasing evidence that multiple factors may alter trait-environment and genotype-environment interactions influencing responses to climate change (Chown et al., 2010;Karl & Trenberth, 2003); for example, CT max has been shown to change with diet (Bujan & Kaspari, 2017). Here, we show that in the presence of multiple stressors, a complex interplay among plastic and evolutionary responses both at physiological and molecular level underpinned population responses to extreme temperatures. Furthermore, we show that the evolutionary advantage of the most recent population, apparent in the constitutive higher CT max in the presence of warming as single stressor, was no longer evident when temperature co-occurred with food or insecticide. Moreover, the regulation of HSPs in the presence of multiple stressors showed unexpected patterns. Although exposure to extreme temperature stress, such as the temperature ramping treatment used here to mimic heat waves, may be expected to induce a shutdown of the regulatory machinery (Kristensen, Loeschcke, & Hoffmann, 2007), we observed both a downregulation and an upregulation of HSPs after the temperature ramping treatment in the presence of multiple stressors. These patterns have been previously observed in a temporal natural population and in an experimental population of D. magna exposed to temperature ramping treatment . Our results and the findings of Jansen et al. (2017) indicate that whereas the HSPs seem to play an important role in the response to extreme temperatures showing a high degree of plasticity, their direction of change may be both up-and downregulation. This may be the result of diverse stressors affecting HSP expression in different directions.
Overall, the direction and the magnitude of plastic changes at molecular and physiological levels in response to extreme temperatures in multi-stress environments were not predictable from populations' response to warming as single stressor. Our results contrasting physiological and molecular responses in single and multiple stress scenarios showed that the co-occurrence of other environmental stressors with temperature has the potential to affect the evolution of thermal tolerance in natural populations of D. magna, either physiologically or at molecular level or both, in ways not directly predictable from the response to warming alone.
Our study is pioneer in studying trade-offs between constitutive and induced thermal tolerance over evolutionary times and in investigating the effect of multiple stressors on thermal tolerance.
However, it suffers from a common limitation to resurrection ecology studies. Whereas these studies provide important insights into the evolutionary processes underlying adaptation in the wild, they require large efforts, and, for this reason, suffer from low replication. These limitations require caution when extrapolating results to other species.
Our results suggest that underestimating the effect of multiple stressors on thermal tolerance can lead to wrong estimates of species evolvability and persistence to future global change. To develop more realistic predictions about the biological impacts of climate change on species persistence, interactions between the mean and variance of environmental temperature, as well as the impact of biotic and abiotic stressors on thermal tolerance, should be considered (Bozinovic, Medina, Alruiz, Cavieres, & Sabat, 2016;Nadeau, Urban, & Bridle, 2017).

DATA ACCE SS I B I LIT Y
Data associated with this study are deposited in the DRYAD databank at the following entry: https://doi.org/10.5061/dryad.3fb854n.