Transcriptomic and functional genetic evidence for distinct ecophysiological responses across complex life cycle stages

Many organisms developing from juvenile life stages through adulthood are faced with changing environmental conditions that differ dramatically but predictably during development. These shifting conditions may include resource availability, predator/herbivore abundance and abiotic factors such as temperature (Krebs and Loeschcke, 1995; Ragland and Kingsolver, 2008; Woods, 2013). To survive these environmental changes, organisms may also dramatically change their morphology, behavior, and physiology across development. For example, juvenile stages often specialize for feeding and growth, while adults primarily (and sometimes exclusively) disperse and mate (Kingsolver et al., 2011; McGraw and Antonovics, 1983; Moran, 1994; Schluter et al., 1991). These developmentally variable environments and key fitness components (e.g. growth versus reproduction) lead to shifting natural selection, which may favor different trait combinations in different life stages (Haldane, 1932; Moran, 1994). This is perhaps most apparent in organisms that metamorphose such as amphibians and holometabolous insects. Their morphology has evolved independently in juvenile and adult stages that inhabit drastically different ecological niches. There are clear physiological differences across complex life cycle stages as well, in part because distinct developmental machinery underlies distinct morphologies and life history strategies across the life cycle (Arbeitman et al., 2002; Herrig et al., 2021; van Gestel et al., 2019). Such morphological and developmental decoupling supports the adaptive decoupling hypothesis, which posits that natural selection favors reduced genetic correlation across developmental stages to allow for stage-specific adaptation (Moran, 1994).

In addition to developmental differences in ‘baseline’ physiology, physiological responses to environmental perturbations may also vary across the life cycle. Many key studies have examined developmental variation in environmental responses by manipulating temperature, a nearly universal selective factor that often varies over the course of development (Bowler and Terblanche, 2008; Jensen et al., 2007; Klockmann et al., 2017). Most of these studies show that thermal responses (survival and various metrics of performance) have very low or absent genetic correlations between juvenile and adult stages of holometabolous insects (Dierks et al., 2012; Gilchrist et al., 1997; Loeschcke and Krebs, 1996; Tucić, 1979). Indeed, our recent studies show that the genetic correlation between juvenile and adult cold hardiness in the fly Drosophila melanogaster is not detectably higher than zero, with no evidence for pleiotropic effects of single nucleotide polymorphism (SNP) variation on thermal performance across metamorphosis (Freda et al., 2017, 2019).

We reason that there are two hypotheses that could explain such extreme genetic decoupling of thermal physiology across development: the ‘developmentally distinct physiology’ hypothesis and the ‘developmentally conserved physiology hypothesis’. The developmentally distinct physiology hypothesis posits that environmental responses may indeed be very different across life stages, mirroring the differences in developmental regulation. The developmentally conserved physiology hypothesis posits that thermal physiology could be largely conserved across development, with only a few stage-specific processes harboring segregating genetic variation. This explanation is less obvious, but still consistent with the observed lack of genetic correlations for environmental physiology across life stages, as we describe below.

Under the developmentally distinct physiology hypothesis, environmental responses would be distinct across life stages, mirroring the differences in developmental regulation across ontogeny. In this scenario, different genes would contribute to environmental responses across stages, with relatively low cross-stage pleiotropy. Though this hypothesis is consistent with the observed lack of genetic correlations across life stages, it would be somewhat at odds with predictions based on the conserved cellular stress response (CSR) (Kültz, 2005).

Conserved cellular and tissue-level responses would argue for the developmentally conserved physiology hypothesis, with many processes universally affecting thermal physiology across development. The CSR is an apparently conserved set of changes in cell physiology in response to a variety of environmental stressors (Kültz, 2005). For example, heat shock proteins (Hsps) and related chaperonins are upregulated in response to multiple stressors, including temperatures that are relatively hot or cold compared with an organism's optimal environmental temperature (Colinet et al., 2010b; Philip and Lee, 2010; Yocum, 2001). If these heat shock responses and other elements of the CSR have a substantial role in whole-organism level environmental responses, then many elements of environmental physiological responses should be very similar across the life cycle.

Developmentally conserved physiology could result in low genetic correlations across life stages if the genetic loci that regulate thermal physiology (e.g. CSR) are highly conserved, and thus not genetically variable. Genetic correlations only assess whether variants at loci affect two traits (e.g. juvenile and adult performance), not whether a given locus itself affects the traits. Thus, these conserved loci would not influence measures of genetic correlations. Rather, a subset of an environmental response may be stage-specific and mediated by genetically variable loci. This scenario could also generate low genetic correlations across life stages.

To test these two hypotheses, we examined physiological responses to cold across the life cycle in D. melanogaster, using two approaches to compare larvae (juveniles) and adults separated by a major metamorphic transition. First, we tested whether whole transcriptome responses to low temperature exposure differ in identity of responding transcripts and/or their temporal patterns of differential expression. Transcriptome sequencing provides a broad snapshot of organismal physiology, and allowed us to assess the similarity of the environmental (temperature) response across the two life cycle stages. Second, we tested whether knocking down a set of nine candidate genes affected response to low temperature in larvae, adults or both. We selected these candidates based on a previous study that found evidence for knockdown effects on cold performance in adult D. melanogaster (Teets and Hahn, 2018). Though the sample of nine genes is relatively small, it provides a first functional test for the presence of stage-specific (consistent with the developmentally distinct physiology hypothesis) or cross-stage (consistent with the developmentally conserved hypothesis) genetic effects on environmental physiology regardless of genetic variability.

We obtained all D. melanogaster (Meigen 1830) lines (Table 1) from the Bloomington Drosophila Stock Center (BDSC; Bloomington, IN, USA) at Indiana University – specific lines used in this study are described below. We reared flies at 25°C, 12 h:12 h light:dark in narrow vials on medium containing agar, cornmeal, molasses, yeast and the antimicrobial agents propionic acid and Tegosept (Genesee Scientific, Morrisville, NC, USA), as described previously (Freda et al., 2017, 2019). We sorted parental flies from appropriate lines (details below) under light CO₂ anesthesia and transferred them into fresh vials containing medium sprinkled with dry, active yeast to facilitate oviposition. We then transferred the parents each day for four consecutive days into fresh vials to produce offspring for use in experiments. The vials from the first egg-laying day were discarded to remove any residual effect of anesthesia on oviposition. We collected third instar larvae and 5-day-old adults for use in both experiments described below. We extracted experimental third instar feeding larvae from cultures 5 days post-oviposition using a 20% w/v sucrose solution and following the protocol of Freda et al. (2017). Experimental adults were collected and sorted into fresh vials under light CO₂ anesthesia 10–12 days post-oviposition (within 1–2 days of eclosion). These flies were held at 25°C, 12 h:12 h light:dark until they were 5 days old to limit any carryover effects of CO₂ exposure (Nilson et al., 2006).

To obtain a transcriptomic metric for how physiology changes during cold exposure and subsequent recovery under benign conditions, we sampled whole-body transcriptomes of third instar larvae and 5-day-old adult D. melanogaster prior to, during and after exposure to a cold temperature (Fig. 1A). We include example (but not a comprehensive list of) predictions of transcriptional patterns that illustrate conserved versus stage-specific responses to thermal stress (Fig. 1A).

We crossed five male and five virgin female flies from each of six Drosophila Genetic Reference Panel (DGRP; Mackay et al., 2012; Huang et al., 2014) isogenic lines (Table 1) to produce offspring for use in this experiment. We initially chose these six lines in order to compare three lines exhibiting high cold tolerance in adults but not larvae, and three lines exhibiting high cold tolerance in larvae but not adults (Table 1; Freda et al., 2017). However, initial analyses revealed little evidence for transcriptomic variation tied to differences between these two classes (high and low cold tolerance) of fly lines, and any weak evidence was highly influenced by outlier lines (Fig. S1). Thus, we treated line (six levels) as a fixed effect (random effects cannot be modeled using the methods that we applied), providing replication across genetic backgrounds, and did not model phenotypic effects in any of our subsequent analyses.

Each experimental replicate consisted of 10 offspring (10 larvae or five male+five female adults). To minimize stochastic environmental effect, each replicate group of 10 offspring was homogenized together to create pools for RNA sequencing. As larvae are prone to desiccation outside of food media, they were extracted from vials in which they were reared and placed in vials containing only a vial flug (Genesee Scientific, catalog no. 49-102) moistened with ddH₂O to inhibit desiccation while adults were tested in a completely empty vial. In this way, we eliminated any effect of thermal inertia that may have been introduced by food while facilitating immediate snap-freezing of both larvae and adults (Freda, 2018). We took an initial sample at 25°C prior to cold exposure (time zero, t₀), then exposed all remaining replicates to −5°C by immediately immersing fly vials in a temperature-controlled recirculating bath (ECO RE 2025, Lauda Corporation, Lauda-Königshofen, Germany). Vials were almost fully immersed (except for the portion containing the flug), ensuring all individuals in a vial were exposed to similar temperatures, i.e. there was no warm air space. We then took samples at 30 and 60 min during the cold exposure (t₃₀ and t₆₀, respectively). At 60 min, all remaining vials were removed from the bath and placed back at 25°C, and one final sample was taken 30 min after this transfer (30 min of recovery, or 90 min total, t₉₀). Flies were transferred to Eppendorf tubes and snap-frozen in liquid N₂ within 2 min after exposure and stored at −80°C for later RNA extraction. The overall experimental design included six lines by two stages by four time points by three replicates, yielding 144 total samples.

To extract RNA from DGRP lines for RNASeq, we homogenized each sample (pool of 10 individuals) with a plastic micropestle in Tri-reagent (Zymo) and used the Zymo Direct-zol total RNA extraction kit according to the manufacturer's instructions. We prepared our cDNA libraries using an RNA-tag sequencing approach, as described previously (Lohman et al., 2016). The resulting libraries were sequenced on five lanes as 100 bp single-end reads on an Illumina HiSeq 2500 at Kansas University's Genome Sequencing Core Laboratory, resulting in an average of 6 million reads per sample. The library size for each sample is available in Table S1. We used STAR (Dobin et al., 2013) to map reads to the D. melanogaster reference genome (version 6.06) obtained from FlyBase (Gramates et al., 2017), with >95% total mapped reads across all samples. Read counts per gene and per isoform were generated using RSEM (Li and Dewey, 2011). After filtering out all gene models not covered by at least one read in 50% of samples, we retained 13,242 genes. Following normalization of read counts across libraries using the weighted trimmed mean of M-values (TMM) method (Robinson and Oshlack, 2010), we examined variation among libraries using a multi-dimensional scaling (MDS) plot generated using the 500 genes with the highest root mean square log₂ fold-change among samples (Robinson et al., 2010). We removed 10 samples that were very clear outliers on the MDS plot (Fig. S2) and exhibited low read counts (less than 200,000 reads) compared with the median read count of 4,808,878 before outliers were removed (Table S1). After removing outliers, all stage×time×line combinations were represented by at least two replicates (Table S2).

Our main goal in experiment 1 was to quantify whether and how the transcriptional response to low temperatures varied between larval and adult life history stages. We expected that many transcripts would be differentially expressed (DE) between life stages because they have very distinct tissue compositions (Arbeitman et al., 2002). Thus, although we estimated gross life stage differences and other contrasts, the parameter of primary interest was a stage×time interaction, indicating stage-specific thermal response during and/or after low temperature exposure (see Fig. 1A for example predicted gene expression trajectories). Below, we detail nested, ad hoc model selection to best estimate that parameter and characterize thermal response trajectories for transcripts with stage-specific expression patterns. The code for this analysis is also publicly available (see Data availability section). We recognize that transcripts/effects removed from these models may also be of interest, but they were not the focus of this study.

We started with a full generalized linear model with binomial error fitted using the edgeR package (Robinson et al., 2010) in R (https://www.r-project.org/) to predict the mean read count for each transcript, then removed effects and transcripts to estimate stage×time (interaction) effects that did not depend on DGRP line. The full model included regression coefficients modelling the effects of stage, time, line, and all two-way interactions and the three-way interaction of these variables. Statistical inferences from this model identified 130 transcripts with a significant (FDR<0.05) three-way interaction term. We then removed all transcripts with significant three-way interactions, then fit a reduced model omitting the three-way interaction, which identified 19 transcripts with a significant two-way interaction between time and line. We removed these transcripts, then fit our final, reduced model including all main effects plus the stage×time and stage×line two-way interactions.

The transcripts of primary interest in our final, reduced model were those that had either (1) a significant main effect of time but no stage×time interaction, or (2) a significant stage×time interaction. The former are transcripts that respond to low temperature in similar ways in both life history stages, while the latter are transcripts that exhibit distinct responses to cold in larvae versus adult flies. We used linear contrasts to estimate the trajectories of differential expression over time for all transcripts in both of these categories by estimating the log₂ fold-change between each time point relative to time zero (t₀). This model also allowed us to identify transcripts that had a significant main effect of stage or a stage×line interaction, but no stage×time interaction. These were not of primary interest, but allowed us to estimate how much of the transcriptome was differentially expressed between life history stages but not responsive to cold.

We used the DAVID functional annotation tool (Huang et al., 2009a,b) to identify functional categories enriched in the set of transcripts illustrating stage-specific responses to cold temperatures (i.e. those transcripts with a stage and a stage×time interaction from the modelling above). This included 763 genes that were differentially expressed in larvae and 121 genes differentially expressed in adults (Table S3). Functional categories included Uniprot keyword searches (UPK), Gene Ontology groups (GO), Interpro protein domains (INTERPRO) and Kyoto Encyclopedia of Genes and Genomes pathways (KEGG). A functional category was considered ‘enriched’ when differentially expressed genes were overrepresented in that category (FDR<0.05).

Transcriptomics from whole bodies are coarse measurements that ignore tissue-specificity of gene expression, and in this case differential expression in response to changing temperatures. However, they provide a comprehensive snapshot of whole-organism physiological responses. We could not directly assess how differences in tissue composition across stages might contribute to different transcriptomic responses without tissue-specific RNA libraries. Rather, we tested whether genes that exhibit high levels of tissue-specific expression in D. melanogaster were overrepresented in sets of transcripts that we identified as differentially expressed between life stages, or exhibiting stage×time interactions.

In order to functionally test whether genes can have stage-specific effects on cold tolerance, we compared the effect of knocking down target gene expression on survival of third instar larvae and 5-day-old adult females and males following a cold stress (Fig. 1B). We include example (but not a comprehensive list of) predictions of how gene knockdown can affect survival if that gene is important for conserved versus stage-specific responses to cold shock (Fig. 1B). Gene knockdown was achieved using TRiP (Transgenic RNAi Project) lines (Table 1) as described in Teets and Hahn (2018). Briefly, five virgin females from each TRiP line carrying dsRNA under the control of a UAS promoter were crossed to five males of a driver line carrying the GAL4 gene under the control of an actin promoter to produce F₁ offspring for experiments. The GAL4 driver promotes expression of dsRNA in all tissues to knock down expression of the target gene in the F₁ generation. We measured survival of groups of 20 larvae or 20 adults (10 adult females and 10 adult males) kept in single fly vials after a 60 min exposure to −5°C, with at least three replicate vials of each stage per line (Fig. 1B). The cold treatment was chosen because 40–60% of control flies (no RNAi) survived this temperature, allowing us to detect effects of RNAi that either increased or decreased survival relative to the control.

We exposed flies to a −5°C cold stress by immersing vials of larvae and adults in a temperature-controlled Arctic A40 recirculating bath (Thermo Fisher Scientific, Denver, CO, USA) containing 50% (v/v) propylene glycol in water. The fly vials for larvae contained fresh medium, and larvae were allowed to burrow into the food prior to cold treatment via holes poked in the medium; the fly vials for adults were empty (Freda et al., 2017). We verified the temperature in vials using a 36-AWG type-T copper-constantan thermocouple (Omega Engineering, Norwalk, CT, USA) interfaced with Picolog v6 software (Pico Technology, Cambridge, UK) via a Pico Technology TC-08 unit. After a 60 min exposure to −5°C, we returned groups of larvae or adult flies to 25°C, 12 h:12 h light:dark to recover. Larvae recovered from cold exposure in the same vials and were monitored for adult eclosion over the next 10 days. We classified larval survivors as those that completed development and eclosed as adults (Freda et al., 2017). Adults recovered in small Petri dishes containing an approximately 1 cm³ piece of fly food medium. We classified adult survivors as those that were motile (could walk/fly independently) 24 h post-cold stress (Jakobs et al., 2015).

Experiment 2 included 11 TRiP lines (Table 1) whose cold tolerance in adult females has previously been characterized: two control (non-RNAi) lines and nine lines that each knocked down expression of a target gene (Teets and Hahn, 2018). None of the genes from the TRiP lines that we chose were differentially expressed in the RNASeq experiment above (Table S3). We reasoned that these genes previously had observable effects on adult cold responses, and thus would provide an appropriate test for whether those responses carry over to other life history stages. Two control lines were required because the dsRNA insertion site (and therefore the genetic background) differed among RNAi lines: four lines (+ one control) had an attP2 insertion site, while five lines (+ one control) had an attP40 insertion site (Table 1).

For each of the nine target genes in experiment 2 (Table 1), we compared the survival post-cold stress of RNAi and control flies with the same genetic background (attP2 or attP40 insertion sites). We used the nlme function in the lme4 package in R (Bates et al., 2014) to fit generalized linear mixed models with binomial error and a logit link function. We modeled survival as a function of the fixed effects of line (control/RNAi), stage (larvae/adult female/adult male) and their interaction, and random (subject level) effects of vial. Example predictions for the effect of RNAi on survival for genes with stage-specific function are given in Fig. 1B.

A large number of transcripts were significantly (FDR<0.05) differentially expressed between larval and adult life stages regardless of time sampled during cold treatment (n=10,966; Fig. 2). A smaller, but still sizeable, number of transcripts changed in abundance over time. However, only 21 transcripts changed in a similar pattern in both life stages (significant main effect of time, no stage×time interaction), while the bulk of the temperature-sensitive transcripts changed over time in a stage-specific manner (n=880 with significant stage×time interaction).

Patterns of change over time (i.e. those induced by cold stress) were also distinct between life stages. Using linear contrasts, we identified many more transcripts that were significantly DE across at least one time point in larvae (n=763) compared with adults (n=121). Subdividing these into transcripts upregulated or downregulated over time on average revealed that most cold-sensitive transcripts in larvae were upregulated or downregulated during the cold exposure and remained at similar levels during recovery (Fig. 3A). In contrast, far fewer transcripts were cold-sensitive in adults, and these were mainly upregulated only during recovery (Fig. 3B). Transcripts that changed significantly over time in larvae did not change over time in adults (Fig. 3A; adult trajectories remain flat). Transcripts significantly upregulated over time in adults did tend to be upregulated in larvae as well (Fig. 3B). However, those larval expression trajectories did not demonstrate the same pronounced upregulation during recovery observed in adults. The small number of transcripts (n=21) with significant main effects of time but no stage×time interaction were upregulated over time in various patterns during cold exposure and recovery (Fig. 4).

Closer examination of genes in several functional categories identify candidate mechanisms underlying the cold response that are also stage-specific. Many functional categories were overrepresented (FDR<0.05) when we performed functional enrichment analysis of the genes that were significantly differentially expressed over time in larvae or adults, and had significant stage×time interactions (Table S3). Below, we focus on members of select categories enriched in either larvae (GO Autophagy, INTERPRO Basic leucine zipper, KEGG Fatty Acid metabolism), adults (UPK Stress response) or both (GO Response to bacterium).

Transcripts participating in autophagy, often involved in clearance of cellular damage and nutrient recycling during energy stress (Kroemer et al., 2010), were mainly upregulated during and after cold exposure in larvae (Fig. 5A). Transcripts with basic leucine zipper domains, largely transcription factors playing roles in developmental regulation, exhibited similar patterns in larvae (Fig. 5B). Transcripts participating in fatty acid metabolism, potentially influencing lipid metabolism or temperature-induced changes in membrane fluidity (Clark and Worland, 2008; Koštál, 2010), were mainly downregulated during and after cold exposure in larvae (Fig. 5C). All of the transcripts in these three functional categories demonstrated little change during and after cold exposure in adult flies. In contrast, transcripts associated with stress response, mainly chaperonins, exhibited the most pronounced changes only during recovery in adults (Fig. 5D). Though some of these transcripts also changed over time in larvae, many were downregulated, including multiple copies of the well-known, temperature-inducible stress response gene Hsp70.

Like transcripts in the stress response category, some transcripts associated with the immune response (within the Response to Bacterium GO group) responded to cold in larvae and adults, though again demonstrating stage-specific patterns (Fig. 5E). The immune response has previously been implicated in responses to thermal extremes in insects (Ferguson et al., 2018; Salehipour-shirazi et al., 2017; Sinclair et al., 2013). All but one transcript in the category were substantially upregulated during cold exposure in larvae, but tended to decrease in relative abundance during recovery. In adults, transcripts for Attacin-C (AttC), two Cecropins (CecA1 and CecA2), Diptericin A (DptA) and Metchnikowin (Mtk) were upregulated at 30 min during cold exposure, downregulated by 60 min, then upregulated again during recovery. One additional Cecropin (CecC) was upregulated over time in a similar pattern for larvae versus adults (main effect of time but no stage×time interaction; Fig. 4).

We found no evidence that transcripts differentially expressed between whole-body extracts from different stages tended to be more tissue specific. Rather, we found a slight tendency for that set of transcripts to exhibit less tissue specificity than chance expectations. The median τ for the set of 10,931 transcripts with stage or stage×line effects (FDR<0.05) was 0.88, and we did not observe a value this small in 10,000 random samples of 10,931 transcripts (median τ of random samples=0.90; P<0.0001). For reference, the distribution of τ is heavily left skewed in D. melanogaster, with most transcripts having high tissue specificity (Fig. S3).

Similarly, we found no evidence that transcripts with stage×time interactions (different responses to the temperature treatments across stages) tended to be more tissue specific. The median τ for the set of 849 transcripts with stage×time effects (FDR<0.05) was 0.82, and we did not observe a value this small in 10,000 random samples of 849 transcripts (median τ of random samples=0.90; P<0.0001).

We observed that the effect of RNAi of target genes on cold tolerance could be stage-specific, although this effect was not universal and was complicated by sex. Three genes of the nine genes tested in this study exhibited clear stage-specific effects of RNAi on cold hardiness (Fig. 6, Table 2). Knockdown of CG10505 or Clk decreased adult, but not larval, survival relative to control flies after a cold stress, suggesting these two genes are important for adult cold tolerance only. Conversely, knockdown of Ir85a decreased larval, but not adult, survival, suggesting this gene is important for larval cold tolerance only. RNAi of three other genes had both stage- and sex-specific effects on cold hardiness (Fig. 6, Table 2). klu knockdown only increased female adult survival, but had no effect on larvae or male adults. CG32533 or mthl15 knockdown had opposite effects on larvae (low survival) and adults (high survival) of one sex only – either female (CG32533) or male (mthl15). The expression of these two genes therefore seems important for larval cold tolerance but detrimental to either female or male adult cold tolerance. We observed no significant effect of RNAi on cold hardiness for the remaining three genes (NtR, pigs or  psh), although pigs and psh trended toward stage-specific effects (Fig. 6, Table 2).

Our results generally support the developmentally distinct physiology hypothesis, showing that both the expression and function of genes pertinent to cold hardiness differ dramatically across development in D. melanogaster. Transcriptional responses to cold in larvae and adults differed in timing (during versus after cold stress), magnitude (many more DE transcripts in larvae) and constituent genes. In addition, of the nine genes whose expression we knocked down via RNAi, most of them (six) affected adult and larval cold hardiness differently. Though differences in tissue composition across life stages probably have some influence on transcriptional responses to cold, they do not appear to account for the majority of whole-organism transcriptional differences in the thermal response across stages. Other studies have demonstrated transcriptional differences across stages in a complex life cycle (Arbeitman et al., 2002; Chevalier et al., 2006; Sanil et al., 2014; Strode et al., 2006), but this is the first study to our knowledge that demonstrates distinct transcriptome-wide environmental responses across life stages, with additional support from functional genetics experiments.

Although classic CSR genes (e.g. Hsps) were not similarly regulated in response to cold in both adults and larvae, we had minor support for the developmentally conserved physiology hypothesis based on transcription of immune response genes. Immunity-related genes have been identified as cold-responsive in a number of other studies of adult drosophilid flies (MacMillan et al., 2016; Sinclair et al., 2013; Vermeulen et al., 2013). However, to our knowledge, this is the first study to find similar results in adults and larvae, both of which upregulated antimicrobial genes. The function of immunity genes in cold-mediated responses remains unknown, though Vermeulen et al. (2013) suggest that some constituent genes may play a role in repair of cellular damage through their known effects on wound healing. The consistency with which these genes are observed in cold responses across species (Cheng et al., 2017; Salehipour-shirazi et al., 2017; Su et al., 2019; Sun et al., 2019), and here across stages, suggests that they play a specific role in cold physiology, and are not just a general stress response in the style of the CSR.

Though some changes in transcription in response to environmental stress undoubtedly have important, adaptive benefits (Chen et al., 2018; Feder, 1999; Feder and Hofmann, 1999; Feder and Krebs, 1998), differences in baseline (unperturbed) physiology may be equally important. In particular, organisms may have higher fitness when exposed to stress because they are physiologically better prepared prior to stress exposure (Hercus et al., 2003; Lecheta et al., 2020; Krebs and Loeschcke, 1994). In other words, genetic variation in genes that directly affect morphogenesis and development may prepare an organism for stress without being differentially expressed during stress. To be sure, we have shown that many transcripts differ in expression between stages in benign (baseline) conditions, but this largely reflects the massive developmental differences between the stages. These data do not allow us to identify which of these differences might contribute to differences in expression during and after stress, or to whole organism performance in response to stress, for that matter. However, to the extent that baseline transcriptomes heavily influence transcriptomic responses to a stressor, this still implies that (baseline) physiology affecting cold performance is distinct between life stages.

We propose that differences between life stages in transcriptomic responses to cold stress are partially driven by differences in cold stress resistance between life stages. Our reasoning is based on an emerging pattern from studies that compare transcriptomic responses to stress among species or populations: species or populations that are the most stress resistant are also the least transcriptomically responsive to environmental stressors. Or, more generally, species or populations that more frequently encounter a given environment tend to have more muted transcriptomic responses to that environment. This has been observed in a Drosophila DGRP population responding to extreme cold (Garcia et al., 2020); Trinidadian guppies responding to predator cues (Ghalambor et al., 2015); fruit-feeding flies responding to different host fruits (Ragland et al., 2015); and in marine invertebrates (Lockwood et al., 2010; Schoville et al., 2012), rice plants (Zhang et al., 2012) and other drosophilid flies (Königer and Grath, 2018; Parker et al., 2015) responding to thermal stressors.

In our study, the more cold-tolerant life stage (adults) is the least transcriptionally responsive to cold, similar to stress-tolerant populations and species listed above. Adult D. melanogaster survive cold stressors better than larvae (Freda et al., 2017, 2019; Jensen et al., 2007), and we observed relatively few cold-sensitive transcripts in adults in the present study. Moreover, the identity of transcripts involved in the larval transcriptomic response suggest more severe cold-induced damage in larvae compared with adults. Larvae differentially expressed autophagy genes during cold stress, suggesting that larvae need to mitigate cellular damage (i.e. degrade damaged cellular components; Kroemer et al., 2010) or to redistribute macromolecules and energy needed for cell differentiation or growth (Neufeld, 2012; Wang and Levine, 2010). We note that Drosophila larvae are quite sensitive to starvation (Juhasz and Neufeld, 2008), so upregulation of autophagy could partially reflect a response to food limitation during the cold exposure. However, larvae rapidly entered chill coma during the cold exposure assays, and immediate nutritional stress was likely minimal. In contrast, adults did not upregulate autophagy-related transcripts and mainly upregulated chaperonins during recovery to preserve cellular function rather than clearing highly damaged cells (Colinet et al., 2010a; Frydenberg et al., 2003; Koštál and Tollarová-Borovanská, 2009). Finally, we note that D. melanogaster larvae are not susceptible to all stressors; they are more heat-tolerant than adults (Freda et al., 2019), likely because they feed in fruits that can become substantially hotter than air temperatures experienced by adults (Feder et al., 1997). Therefore, we do not think our results reflect general stress susceptibility in larvae, e.g. simply because they are undergoing rapid cellular growth, division and differentiation compared with adults. Indeed, if future studies compare transcriptomic responses of D. melanogaster to heat stress, we predict that larvae would have a muted response relative to the less heat-tolerant adults.

Larvae rapidly differentially regulated a relatively large number of transcripts both during and following cold exposure. These changes likely include active regulation in response to cellular damage, as evidenced by the aforementioned autophagy response. We also observed differential regulation of lipid metabolism in larvae. Fatty acids are important in energy storage (as part of triacylglycerides) and membrane fluidity (as part of phospholipids) (Denlinger et al., 2010). Larvae downregulated several desaturases (e.g. Desat1, CG8630, CG9743), suggesting that they are not increasing the abundance of unsaturated fatty acids in phospholipids to maintain membrane fluidity at low temperatures (Ohtsu et al., 1998; Overgaard et al., 2005). However, the downregulation of several enzymes associated with fatty acid catabolism (e.g. ACOX1) and fatty acid synthesis (e.g. ACC, acsl, bgm) is consistent with restructuring of lipid metabolism. The downregulation of fatty acid synthesis genes suggests that catabolic (rather than anabolic) metabolism may be favored during and after cold exposure to potentially support growth or recovery from stress (Sinclair and Marshall, 2018).

In contrast, adults had relatively muted transcriptomic responses during cold exposure, with a limited (in number of transcripts) but robust (in the degree of differential expression) response during recovery. Upregulation of genes following (rather than during) both cold and heat exposure has been documented in other studies of D. melanogaster adults (Colinet et al., 2010a; Sinclair et al., 2007; Sørensen et al., 2005). The best-characterized gene expression response to temperature stress, hot or cold, is upregulation of Hsps and other chaperonins during recovery after exposure to a stressor (Colinet et al., 2010b; Philip and Lee, 2010; Yocum, 2001). This was the most prominent adult response in our study as well, with no detectable changes in Hsp expression during cold exposure. As mentioned above, it is likely that the relative stability of gene expression during stress in adults reflects less severe perturbations from homeostasis and more restricted cellular damage.

Though the transcriptome is only one metric of physiology, the scale of the differences across stages in this study suggests that allelic variants in many genes could strongly affect environmental sensitivity of one stage, while having little effect on other stages. Our results are entirely consistent with empirical studies that repeatedly show little to no genetic correlation in environmental (thermal) sensitivity across metamorphosis in insects (Dierks et al., 2012; Freda et al., 2017, 2019; Gilchrist et al., 1997; Loeschcke and Krebs, 1996; Tucić, 1979). In combination, these results support the developmentally distinct physiology hypothesis and suggest that strong genetic decoupling of environmental sensitivity is relatively common for organisms with complex life cycles, likely facilitating adaptation/acclimation of different life stages to different thermal environments.

Being so widespread, differences in stage-specific thermal tolerance might not appear so surprising. However, temperature is fundamental to limiting species' spatial distributions (Bale, 2002; Bale et al., 2002), and thus thermal performance must be constrained in some ways. The results of our RNAi knockdown experiments suggest that cross-stage pleiotropy for environmental sensitivity is not widespread, which (similar to the transcriptomic results) is also consistent with the developmentally distinct physiology hypothesis. However, we expect that some cross-stage pleiotropy exists, and will constrain the limits of thermal flexibility across life stages. For example, genetic modifications to increase Hsp70 copy number (and subsequently expression) affected both larval and adult thermal tolerance in D. melanogaster (Krebs and Bettencourt, 1999). In that instance, the genetic differences between modified and non-modified lines were relatively extreme (12 extra gene copies). However, there is some evidence for cross-stage effects of naturally segregating genetic variants in plants. Quantitative trait locus (QTL) studies in rice have identified QTL associated with cold tolerance at multiple developmental stages, though most QTL only affect a single developmental stage (Yang et al., 2020).

Given the polygenic architecture of environmental tolerance in general (Healy et al., 2018) and thermal tolerance specifically (Barghi et al., 2019; Freda et al., 2017; Sanghera et al., 2011), it is unlikely that further, detailed analysis of single-locus pleiotropy will fully address questions about the limits of stage-independent adaptations to environmental stressors. Rather, comparative studies leveraging existing variation in stage-specific adaptation or selection studies generating relevant phenotypic variation would seem to be the most promising avenues for further research.

The authors would like to thank Maizey Funk and Lahari Gadey for help with D. melanogaster rearing, and the ant-fly-beetle group at CU Denver for comments on earlier versions of the manuscript.