The ghost of temperature past: interactive effects of previous and current thermal conditions on gene expression in Manduca sexta

Exposure to high temperatures can have significant negative consequences for the performance of ectotherms. However, an organism's response to high temperature stress depends not just on current conditions but also on previous thermal exposure of the individual and the evolutionary history of its population. Predicting the impact of thermal history on an organism's response is complex because it may depend on interactions between timing, duration and variation of past temperature exposures, i.e. time-dependent effects (Cavieres et al., 2018; Kingsolver et al., 2015; Kingsolver and Woods, 2016; Sinclair et al., 2016). Time-dependent temperature effects can have both positive and negative implications for an individual, depending on factors such as developmental stage or time scale over which thermal exposure occurs. For example, brief exposure to high but sub-lethal temperatures can be beneficial by increasing tolerance to subsequent heat stress, a form of short-term adaptive acclimation known as heat hardening (Huang et al., 2007; Kregel, 2002; Li et al., 2018). However, repeated exposure or a longer duration of exposure to the same sub-lethal temperatures can have a negative effect by increasing stress and leading to a reduction in other performance measures, such as decreased growth rate or survival (Folguera et al., 2011). Similarly, time-dependent effects have also been documented for cold temperature exposure (Colinet and Hoffmann, 2012; Marshall and Sinclair, 2012, 2015).

Many studies have assessed the impact of brief high temperature exposure over short time scales and constant high temperatures over longer periods (Angilletta et al., 2010; Huey et al., 2012); however, the role of diurnal temperature fluctuations during development for time-dependent effects has received less attention. Understanding the impact of fluctuations is important given that ectotherms in the field are likely to experience daily temperature variation. Studies comparing fluctuating versus constant acclimation temperatures have found substantial differences in phenotypic responses including thermotolerance (Bozinovic et al., 2011; Cavieres et al., 2016; Colinet et al., 2015; Davis et al., 2006; Ghaedi and Andrew, 2016; Salinas et al., 2019). Other recent experimental studies have shown that ignoring time-dependent effects can lead to poor predictions of growth and development rates for ectotherms in fluctuating temperatures (Carrington et al., 2013; Ketola and Saarinen, 2015; Kingsolver et al., 2015; Niehaus et al., 2012). In order to better understand the response of ectotherms to current variation in thermal conditions and predict the impact of future climate change, more work is needed to identify the specific mechanisms underlying time-dependent temperature effects using ecologically relevant temperature conditions.

Temperature and other environmental stresses cause major changes in gene expression and cell physiology; these changes are collectively known as the cellular stress response. Previous studies have identified specific changes in gene and protein expression in response to non-lethal heat shock for a variety of organisms. Most notably, the induction of heat shock proteins (HSPs) during thermal stress has been well documented in many cases (Feder and Hofmann, 1999; Hofmann, 1999; Lindquist and Craig, 1988; Sørensen et al., 2003; Tomanek, 2010; Tomanek and Somero, 2002). Expression levels for heat shock protein genes (hsp genes) tend to increase sharply initially and then decrease to pre-shock levels within several hours while elevated levels of the proteins (HSPs) can be detected in cells multiple days after heat exposure (Bahrndorff et al., 2009; Dahlgaard et al., 2002; Karl et al., 2012). While there is evidence that elevated levels of HSP are important for short-term heat hardening (Dong et al., 2010; Hu et al., 2014; Huang et al., 2007; Kregel, 2002), continued overexpression of HSPs over the long term can have deleterious consequences on performance (Chevin and Hoffmann, 2017; Kafri et al., 2016; Krebs and Feder, 1997). However, several studies have linked elevated baseline expression of HSPs to protection against heat stress over seasonal time scales (Hamdoun et al., 2003; Lund et al., 2006). Sørensen et al. (2017) found that basal HSP expression was variable among Drosophila melanogaster populations with differing selection regimes; therefore, constitutive HSP expression levels might be shaped by evolutionary pressure from specific environmental stresses. Given the conflicting evidence, however, it remains unclear to what extent HSPs may play a role in protection against longer-term heat exposure, such as during developmental acclimation.

Expression patterns of other genes including those involved in metabolism, antioxidant production, transport activity, immune response and transcription also appear to be strongly altered by short-term heat stress (Colinet et al., 2013; Cui et al., 2019; DeBiasse and Kelly, 2016; Franke et al., 2019; Gleason and Burton, 2015; López-Maury et al., 2008; Schoville et al., 2012; Sørensen et al., 2005; Vihervaara et al., 2018). Similar pathways may also be involved in longer-term acclimation under constant temperatures (Metzger and Schulte, 2018). However, most transcriptomic and proteomic studies have only considered the effects of single heat shocks or response to constant thermal stress. Only a handful of studies have assessed how fluctuations in thermal conditions during development impact the cellular stress response, especially for genes other than hsp genes (Bellantuono et al., 2012; Barshis et al., 2013; Kenkel and Matz, 2017; Podrabsky and Somero, 2004; Lima and Willett, 2017). Even fewer studies have evaluated gene expression in fluctuating thermal conditions for terrestrial ectotherms (but see Sørensen et al., 2016b). Further investigation is warranted given that Sørensen et al. (2016b) found that thermal acclimation induced by differences in mean temperature versus fluctuations in temperature appears to be achieved via independent developmental pathways in Drosophila. Additionally, many of the previous studies evaluated gene expression using a mix of tissue types, rather than a tissue-specific approach. While mixed-tissue sequencing is useful for determining the overall transcriptomic response of an organism (and may be necessitated by a particular study organism's body size), such an approach could potentially miss important, fine-scale expression patterns if different tissues vary in their responses to stress.

To better understand how diurnal fluctuations throughout development impact gene expression patterns in ectotherms, we used the tobacco hornworm, Manduca sexta (Linnaeus 1763) (Lepidoptera: Sphindidae). Manduca sexta larvae have been used extensively to study development, physiology and thermal biology in ectotherms (Casey, 1977; Heinrich, 1970; Kingsolver and Woods, 1997; Nijout, 1999; Reynolds and Nottingham, 1985; Riddiford et al., 2003). These experimental studies have primarily focused on measures of growth, development rate and survival to quantify thermal responses, and have documented time-dependent effects of constant and fluctuating rearing temperatures on thermal tolerance and physiological acclimation (Kingsolver et al., 2015; Kingsolver and Woods, 2016). Additionally, M. sexta larvae exhibit heat-hardening responses (increased critical thermal maximum, CT_max) to short-term (24 h) heat shocks of 35–42°C, and show upregulation of HSPs at 38°C with maximal HSP synthesis at 42°C (Fittinghoff and Riddiford, 1990); higher mean or daily maximum temperatures during development can also increase CT_max (Kingsolver et al., 2016). The recently published M. sexta genome (Kanost et al., 2016) opens the door for next-generation sequencing-based studies of the genetic mechanisms involved in complex temperature responses. Here, we considered two populations of M. sexta that show differences in their thermal sensitivity: a field population collected from central North Carolina, and a lab colony originally derived from a genetically similar North Carolina field population but maintained at constant temperature (25–26°C) in the lab for more than 300 generations. Under field conditions for M. sexta in the southeastern and southwestern USA, larvae can regularly experience fluctuations in body temperature of 20°C or more over the course of a day and can reach body temperatures of up to 38–40°C during the summer (Casey, 1976; Kingsolver and Nagle, 2007). The lab population experiences little to no temperature variation and shows reduced ability to deal with high temperature stress, likely due to relaxed selection to temperature (Kingsolver and Nagle, 2007; Kingsolver, 2007). There is also evidence that the two populations show differences in their thermal acclimation responses (Kingsolver et al., 2015; Kingsolver lab, unpublished).

In this study, we used RNA sequencing (RNA-Seq) of fat body tissue to evaluate the effect of diurnal temperature fluctuations during development on M. sexta gene expression both independently and in conjunction with subsequent temperature stress. We also compared these results between the lab and field populations to determine the impact of evolutionary history on these organisms' thermal stress responses. We used fat body tissue specifically because it shows a robust HSP response to heat stress in M. sexta larvae (Fittinghoff and Riddiford, 1990) and plays a central role in many insect metabolic functions that might be affected by temperature (Arrese and Soulages, 2010). Based on previous thermal studies of these two M. sexta populations, we proposed several hypotheses for our current study. First, given that lab larvae show increased sensitivity to thermal stress compared with field larvae, we expected to see more differential expression of known stress pathways for lab versus field caterpillars during thermal stress treatments within relevant ecological ranges. Similarly, we expected that repeated high temperature fluctuations should have greater negative consequences for lab caterpillars, potentially resulting in long-term perturbations of gene expression patterns extending beyond immediate recovery from a thermal exposure. Finally, if acclimation ability is reflective of selective pressure for the environmental conditions experienced (Masel et al., 2007), we predicted that the field population should show greater capacity to acclimate. This would be demonstrated by fewer stress-related genes differentially expressed for field caterpillars at high temperatures after repeated exposure compared with novel exposure. However, some past studies have failed to find support for the prediction that greater environmental variability necessarily results in increased plasticity (Fragata et al., 2016; Manenti et al., 2015; Manenti et al., 2017; Sørensen et al., 2016a). In a broader context, this study will help to build an improved framework to predict temperature effects outside of controlled laboratory conditions. Transcriptomic studies using variable, ecologically relevant temperature conditions are still underrepresented in the literature, but are important for understanding the responses of ectotherms to current variation in thermal conditions and predicting the impact of future climate change.

Eggs were collected from two M. sexta populations: a field population [Upper Coastal Plain Research Station in Rocky Mount, NC, USA (35.89°N, 77.68°W) during July–August of 2016] and a laboratory colony (maintained at the University of North Carolina – Chapel Hill since the 1980s and originally established from field populations near Raleigh, NC, USA, in the 1960s). Eggs collected from the field were kept in coolers during transport back to the lab to prevent overheating and eggs from both populations were maintained at 25°C until hatching. Sex of larvae was not determined, but sex ratios are assumed to be relative equal for all treatments.

In order to assess how previous exposure to diurnal temperature fluctuations affects gene expression under thermal stress, and also to determine the effects of rearing temperature fluctuations and thermal stress independently, our experimental design included eight different temperature treatments (Fig. 1). Various combinations of these treatments were designed to test for specific comparisons of interest. These comparisons were: rearing temperature, test temperature, interaction of rearing and test temperature, and sampling time (as a control for samples collected at different time points). By also comparing lab and field populations for each of these comparisons, we were able to evaluate the impact of evolutionary history on gene expression due to the recent relaxation of thermal functional constraint.

Upon hatching, M. sexta larvae from both lab and field populations were placed on Lepidopteran artificial diet containing a small amount of dried tobacco leaf (Kingsolver and Nagle, 2007) and assigned to one of three rearing temperature regimes carried out in programmable environmental chambers (Percival I-36VL): constant 30°C (control), fluctuating 30±5°C or fluctuating 30±10°C. 30°C was chosen as the control/mean temperature because it is near M. sexta’s thermal optimum (Kingsolver and Nagle, 2007; Reynolds and Nottingham, 1985); 35 and 40°C were chosen as daily maximal temperatures because 30±10°C during development strongly reduces larval growth and developmental rates, especially for the laboratory population (Kingsolver et al., 2015). Prolonged exposure to constant 35°C conditions results in significant mortality of lab larvae compared with field larvae, but short-term exposure has no detectable negative effects for either population (Kingsolver and Nagle, 2007; Reynolds and Nottingham, 1985). By contrast, short-term exposure to 40°C reduces growth and development rates for both populations, and prolonged exposure leads to death within a few days (Kingsolver lab, unpublished). Caterpillars in fluctuating rearing conditions recurrently experienced 2 h at maximum and minimum temperatures over a 24 h period (maximum at 13:00–15:00 h and minimum at 01:00–03:00 h) with linear ramping in between. A 14 h light:10 h dark photocycle was used for all rearing treatments. Caterpillars were reared under these conditions until reaching the second day of the 5th instar stage. At that time (at 08:00 h when all rearing treatments were at 30°C), a subset of caterpillars were frozen in liquid nitrogen (time point 1 – hereafter T1) and the remaining caterpillars were assigned to a second temperature treatment phase (Fig. 1).

The remaining caterpillars were then exposed to a constant 30°C (control), or a linear ramp up to either 35 or 40°C (total time: 7 h, 2 h at maximum temperature) before being collected and frozen (time point 2 – hereafter T2). For some of the larvae, 35 or 40°C was a novel temperature exposure (heat shock), while other larvae had previously been exposed to these temperatures during the rearing phase. Although there were a total of three rearing treatment temperatures and three test treatment temperatures, not all possible combinations of these were used (Fig. 1). Nine field and nine lab caterpillars were subjected to each of the eight treatment conditions (72 caterpillars from each population, 144 caterpillars in total). The eight temperature treatments are represented using a three-part naming system in which the first part indicates the magnitude of temperature fluctuations that were experienced throughout larval development, the second part denotes the time of sampling, and the third part refers to the temperature at sampling (Fig. 1C,D). For example, (±0)T1.30 indicates that the larval sample experienced no diurnal fluctuations during the rearing phase and was frozen at time point 1 at a temperature of 30°C.

RNA was obtained from fat body tissue of each of the frozen caterpillar samples. A tissue-specific sequencing approach was chosen to ensure that gene expression patterns were not obscured if different tissues varied in the magnitude or direction of their response to thermal stress. Fat bodies were dissected from the last two or three pro-leg segments and stored in 1.5 ml Eppendorf tubes in 100 µl RNA-later ICE at −20°C. RNA extractions were performed using the Direct-zol RNA MiniPrep Kit (Zymo Research) according to the manufacturer's protocol. RNA samples were quantified using a Qubit 1.0 fluorometer (Invitrogen/Thermo Fisher Scientific). Equal amounts of RNA from three caterpillars from the same treatment conditions were pooled into final samples with a total RNA concentration of 2 µg. This resulted in three biological replicates, each made up of three pooled individuals, for each of the different treatment combinations (8 temperature treatments×2 populations×3 biological replicates=48 samples in total). Final quality and quantity of RNA samples were checked with a Bioanalyzer 2100 (Agilent Technologies). Library preparation and 125 bp single-end Illumina Hiseq 2500 sequencing were conducted at the North Carolina State University Genomic Sciences Laboratory (Raleigh, NC, USA). Stranded mRNA libraries were generated using the NEBNext Ultra Directional RNA Library Prep Kit (New England BioLabs) according to the manufacturer's instructions. The 48 samples were randomly split into two groups to prevent sequencing bias and each was barcoded with one of 24 unique sequence indexes (TruSeq LT adapters, Illumina). These pools of 24 samples were each run on three different lanes of a single flow cell (six lanes total), resulting in three technical replicates for each unique biological sample.

Quality of the sequences was checked using FastQC v0.11.5 (Andrews, 2010). Adapter sequences and low quality bases (quality score cutoff of 20 for bases at either end of read and average score over 5 base sliding windows) were trimmed using Trimmomatic v0.36 (Bolger et al., 2014). Trimmed reads shorter than 50 bp were discarded. Reads were then mapped to the M. sexta reference genome (Kanost et al., 2016) using HISAT2 v2.1.0 (Kim et al., 2015). A list of known splice sites generated from M. sexta transcriptome assembly was used to aid alignment program, but otherwise default alignment and scoring options were used. Counts (number of reads that mapped to each gene) were generated using HTSeq v0.6.1 (Anders et al., 2015) using the M. sexta transcriptome assembly. Counts for alternatively spliced transcripts were collapsed to the gene level. Reads with multiple alignments or that map to regions that overlap more than one gene were excluded. Because the reference genome is likely more similar to the lab population, we verified that there were no differences in the proportion of secondary mappings between reads of lab and field samples. This ruled out the possibility that differences in mapping efficiency introduced bias to the differential expression analysis.

Identification of differentially expressed (DE) genes was performed using DESeq2 v1.16.1 (Love et al., 2014) in R v3.4.1 (http://www.R-project.org/). Read count data was normalized across samples to account for variation in sequencing depth and the three technical replicates for each unique sample were collapsed together after confirming that there was no significant effect of sequencing lane. Differential gene expression was tested using a negative binomial model which accounted for both population and temperature treatment effects (i.e. 16 different treatment variables, each with three biological replicates). The number of significantly DE genes for a subset of the possible pairwise comparisons between the 16 different treatment combinations was assessed using a stringent false discovery rate (FDR) cutoff of <0.01 and a log₂ fold-change (LFC) in expression >1. This cutoff was selected because a more permissive cutoff resulted in an unacceptably high number of predicted false positives; however, the specific cutoff used did not substantially affect interpretation of the results. Analysis of DE genes shared between different pairwise comparisons was done using an online Venn diagram web tool, Venny v2.1 (https://bioinfogp.cnb.csic.es/tools/venny/index.html).

Functional annotation of M. sexta genes was performed by blasting (BLASTX) to NCBI's non-redundant (nr) database (e-value threshold of 1.0E−3 and a maximum of 10 target sequences). The highest blast hit was used to assign the gene description. Blast2GO v5.1.1 (Conesa et al., 2005) was then used to assign gene ontology (GO) terms to blast results (e-value threshold of 1.0E−6, annotation cut-off >55 and GO weight >5). Annotation of GO terms was increased by running an InterPro search against all EBI databases to retrieve GO terms associated with domain/motif information. GO enrichment analysis of significant DE genes was assessed with Fisher's exact test feature in Blast2GO. (A significance cutoff of FDR <0.05 often excluded most results, so GO terms with FDR <0.1 were also considered.)

Illumina sequencing resulted in an average final read count of 7,043,708 single-end reads per sample that uniquely mapped to the M. sexta reference transcriptome after filtering. Once technical replicates were collapsed, the average number of reads per sample was 21,131,123. Approximately 89% of trimmed reads uniquely mapped to one place in the genome, and ∼67% of trimmed reads were included in finalized counts – this excludes reads that mapped to the genome but not the transcriptome and ambiguously mapped reads. The final dataset used for all subsequent differential expression analysis contained 14,756 genes with non-zero read counts out of 15,542 total annotated genes in the M. sexta reference transcriptome; 68 of these genes were identified as outliers by DESeq2 using Cook's distance method and excluded from differential expression analysis. Of the total annotated genes, 14,727 (94.8%) had BLAST hits and 11,600 (74.6%) were assigned GO terms.

Overall gene expression patterns were markedly divergent between the M. sexta field and laboratory populations, as indicated by a principal component analysis conducted using the top 500 most variable genes (Fig. 2). Pairwise comparisons between field and lab samples identify large numbers of significant differentially expressed (DE) genes across each of the eight treatment conditions (Table S1), although some treatments showed more similarity in response than others (average number of significant DE genes: 912.5, range: 213–2475). More specifically, treatments involving 35°C seemed to result in the largest differences in gene expression between field and lab populations. In most cases, more genes tended to be significantly upregulated in the field compared with the lab samples. Upregulated DE genes in field versus lab comparisons were enriched in GO terms associated with oxidoreductase and monooxygenase catalytic activity, iron/heme binding, DNA integration and chitin metabolism. Downregulated DE genes were not found to be enriched for any annotated GO category, but were enriched for genes with unknown functions (Table S2).

In order to account for any effects of sampling time, we tested for changes in gene expression between T1 and T2 at identical rearing and test temperatures (Fig. 1). Relatively few genes (62 and 67 genes for lab and field, respectively) were found to be significantly DE at T2 compared with T1 when sampled under control temperature conditions [(±0)T2.30 versus (±0)T1.30; Table 1]. Of these genes, only 6 were shared in both lab and field comparisons and neither gene set was enriched for any GO categories. Because the effect of sampling time on gene expression appeared to be minor, we ignored this effect for the remaining analyses.

The effects of diurnal temperature fluctuations during development on gene expression were assessed by comparing rearing treatments with fluctuating versus constant temperatures. These treatments were all sampled at the same time (T1) and temperature (30°C) (Fig. 1). Temperature fluctuations alone did not dramatically alter gene expression on a longer scale (17 h after high temperature exposure) although there was some variation in the number of DE found across different treatment comparisons [(±5)T1.30 and (±10)T1.30 versus (±0)T1.30; Table 1]. The 30±10°C rearing treatment for the lab population showed the largest response, with 186 genes significantly DE compared with the constant control (Fig. S1). Notably, almost no hsp genes were found to be DE in any of these rearing comparisons. Significant GO term enrichment was found only for the lab population exposed to a 30±10°C rearing treatment. Significantly upregulated genes in this comparison were enriched for amino acid and nucleic acid biosynthesis, proteolysis and extracellular matrix components (Table S2).

To determine the effect of test temperature alone, without any effect of rearing temperature, we compared caterpillars that were all reared under constant 30°C conditions and then exposed to a higher novel test temperature (Fig. 1). Large numbers of genes (333–778) were found to be significantly DE in response to a heat shock of 35 or 40°C for both populations [(±0)T2.35 and (±0)T2.40 versus (±0)T2.30; Table 1]. Of these genes, those encoding HSPs tended to show the largest LFC in expression (Fig. 3). While all four comparisons show significant changes in gene expression due to heat shock, surprisingly only 23 genes were shared across the four gene sets (Fig. 3E; Table S4). Of those genes involved in the shared heat shock response, more than half (12) were hsp genes. The two lab population comparisons had the most similarity in gene expression response (204 shared genes), followed by the 40°C comparisons (143 shared genes). The response of the field population at 35°C was notably unique with the fewest DE genes shared with other novel heat shock comparisons. Of the significantly enriched GO terms found, none were shared among all sets of upregulated genes other than the hsp response. All four sets of downregulated genes were enriched for general metabolic processes, but no specific processes were shared by all sets. Some functional responses were shared by two of the sets (Table S3). The unique response of the field population to a 35°C heat shock was largely due to the upregulation of a large number of genes (many of these genes had LFC of a magnitude similar to or greater than that of the hsp genes); however, this gene set was only found to be enriched for chitin metabolism (which makes up a very small proportion of total significantly upregulated genes) and genes with unknown functions (Table S2).

Thermal history (exposure to temperature fluctuations during rearing) had little effect on gene expression in response to subsequent exposure to 35°C for the lab population and to 40°C for both lab and field populations. This is indicated by small numbers of DE genes (7–33) for comparisons between novel heat shock and thermal history [(±0)T2.35 versus (±5)T2.35 and (±0)T2.40 versus (±10)T2.40; Table 1 and Fig. 4]. However, the gene expression response of the field population to 35°C did depend on thermal history (i.e. ±5°C fluctuations during rearing reduced the heat shock response at 35°C). For this comparison, 258 genes were found to be significantly DE among samples that differed only in rearing treatment conditions. Of these genes, 198 (∼77%) were the same as those DE for the field population in the (±0)T2.35 versus (±0)T2.30 comparison. There were no significant GO enrichment results for any novel heat shock versus thermal history comparison except for the field population at 35°C. Upregulated genes in this comparison were enriched for chitin metabolism and unknown function, and downregulated genes were enriched for defense response (Table S2).

We predicted that M. sexta lab and field individuals would show differences in gene expression during heat stress as a result of greater expression changes to stress-related pathways in the lab population. However, an alternative prediction could be argued – that the field population might instead show a stronger stress response, indicating a greater capacity to accommodate stress. While these populations did have substantial divergence in their gene expression response at high temperatures, the differences cannot be fully explained by DE of stress-induced genes alone; non-stressful (30°C) temperatures also showed surprisingly large gene expression changes between populations. The greatest difference in response between populations was seen after exposure to 35°C, and there were relatively fewer differences at 40°C (Table S1). This pattern could be explained by observed thermal phenotypes for these populations which indicate that 35°C causes more stress to lab compared with field larvae, but 40°C seems more uniformly stressful to both (Kingsolver and Nagle, 2007; Kingsolver lab, unpublished). Such differences in thermal sensitivity at 35°C compared with 40°C between populations may therefore correspond to divergence and convergence of expression in thermal stress response pathways, respectively. It is unclear, however, which population exhibited a stronger overall stress-related expression response to 35°C because many of the DE genes do not have clear functional annotations.

Even under non-stressful control conditions, hundreds of genes are significantly DE between the two populations. This is consistent with findings that in addition to differences in thermal tolerance, lab and field populations have diverged in multiple other phenotypic traits as a result of domestication. Lab larvae have evolved traits including increased size and growth rate via active selection during domestication, while traits such as immunity and host chemical response have been reduced, likely as a result of relaxed selection in the laboratory environment (D'Amico et al., 2001; Diamond and Kingsolver, 2011; Kingsolver and Nagle, 2007; Kingsolver, 2007). Interestingly, a larger proportion of DE genes showed higher expression in the field than in the lab population (Table S1). Divergent gene expression in domesticated versus wild species has also been observed in rainbow trout (Tymchuk et al., 2009); however, the number of upregulated and downregulated genes in domesticated compared with wild fish was tissue specific, with some tissues showing more genes upregulated in domesticated trout. Another study found an overall pattern of decreased gene expression diversity for multiple domesticated versus wild plant and animal species, which is possibly associated with decreased genetic diversity in domesticated organisms (Liu et al., 2019). Because we used only fat body tissue in our study, we cannot confirm that the pattern of greater expression in field compared with lab populations holds for other tissue types. Some of the divergence in gene expression patterns between field and lab populations could also be explained by effects of genetic background, such as higher levels of drift and inbreeding in the laboratory colony (Sarup et al., 2011), and other unpredictable changes to thermal phenotypes associated with lab adaptation (Maclean et al., 2018). Although we attempted to limit environmental differences between the lab and the field by rearing both populations from eggs in controlled laboratory conditions, it is possible that pathogens, early egg environment differences or transgenerational effects could also contribute to differences in gene expression between the lab and field populations. Some immune genes did show interesting expression differences between populations. Antimicrobial peptide activation pathways (Toll and IMD) were upregulated in the field compared with the lab population under control temperatures, potentially indicating stronger innate antimicrobial defenses in the field population. Additionally, the field population showed downregulation of many immune genes compared with the lab population at 35°C and to some extent at 40°C. This appears to be the combined result of the field downregulating these immune genes specifically at 35°C, and the lab upregulating immune genes possibly in response to temperature stress at 35°C. Alternatively, expression of immune genes could be highly sensitive to other environmental perturbations and may not necessarily reflect a functional response to temperature. While not the focus of this study, future work could look at additional genes and pathways underlying phenotypic differences in lab and field M. sexta, as well as patterns of genetic divergence between populations.

We hypothesized that if repeated fluctuations altered gene expression beyond immediate heat exposure, these effects should be greater for the lab population. Our results do provide some support for this hypothesis but, overall, fluctuating rearing conditions alone had relatively little effect on gene expression patterns over the long term. While we saw large-scale changes in gene expression when caterpillars were sampled at 35 or 40°C, these changes were short lived and most DE genes returned to control levels within hours after thermal exposure. In particular, the upregulation of hsp genes was not maintained over the 17 h period after a heat shock, mostly consistent with the findings of previous studies that show hsp expression returning to control levels 23–32 h after heat exposure (Bahrndorff et al., 2009; Dahlgaard et al., 1998; Karl et al., 2012). However, because hsp gene expression is not necessarily directly related to HSP protein levels in cells (Bahrndorff et al., 2009), future work could quantify differences in HSP protein levels between populations in fluctuating conditions.

While three out of our four rearing treatment comparisons had low numbers of significantly DE genes, the lab population reared at 30±10°C did show greater differences compared with the constant-temperature control. Consistent with our prediction that fluctuations will increase the stress response for the lab compared with the field population, a relatively large proportion of these genes (87 out of 186) were the same as those DE in at least one novel heat shock response comparison. This could indicate that the lab population has a harder time recovering from repeated heat shocks of 40°C, resulting in an elevated thermal stress response maintained for longer periods of time. Although only nine (out of 48) DE genes for the field 30±10°C rearing comparison were also shared with other novel heat shock comparisons, we cannot rule out the possibility that ±10°C fluctuations for the field population also result in somewhat increased stress. Even if there is an indication of a slightly elevated stress response in the field population, it is still much less than that seen in the lab population for these conditions. Sørensen et al. (2016b) found that diurnal fluctuating rearing temperatures caused relatively small transcriptional effects in Drosophila, while Lima and Willett (2017) saw larger effects of fluctuations on gene expression in copepods. Our results are somewhat mixed, but suggest that small transcriptional changes due to fluctuations can still result in substantial differences in gene expression response to later thermal exposure for conditions in which developmental acclimation is possible.

Both lab and field populations showed a thermal stress response during novel exposure to 35 and 40°C, with significant upregulation of many of the same hsp genes in all cases. Previous studies have found that HSPs can be induced at 38°C in M. sexta (Fittinghoff and Riddiford, 1990). Here, we saw induction of hsp transcripts at 35°C, even in field individuals. The number and fold-change of upregulated hsp genes was similar for the two populations at 40°C, but was slightly elevated (in number and to some extent fold-change) for field compared with lab individuals at 35°C. This pattern is contrary to our original prediction, but offers some support for the alternative interpretation that a larger stress-related response, for hsp genes at least, correlates with a greater capacity to accommodate stress. Increased hsp gene expression has been associated with increased heat tolerance in Drosophila larvae and intertidal copepods (Krebs and Feder, 1997; Lima and Willett, 2017; Schoville et al., 2012). The elevated hsp gene expression levels at 35°C for field individuals in our study are therefore consistent with their increased thermal tolerance compared with lab individuals over long-term exposure to this temperature (Kingsolver and Nagle, 2007).

An unexpected finding was that novel heat shock at 35°C for field individuals caused changes in gene expression patterns that were strikingly different from the other three comparisons (Fig. 3C). Many of these significantly upregulated genes showed fold-changes that were large in magnitude, similar to those of hsp genes, and the majority of these significantly DE genes were not shared by any other heat shock comparisons and were not enriched for other known thermal stress pathways. Several cuticular proteins (CPs) were among the list of DE genes with notably high expression. Other studies have also found some CP genes upregulated in response to thermal stress, especially in populations with higher thermal tolerance limits (Chen et al., 2018; Schoville et al., 2012). CPs could play a role in protection against stress by thickening the cuticle to prevent water loss. While CPs are known to be expressed in fat body tissue as well, it is unclear how they might contribute to stress tolerance in this context. Of the 427 genes significantly DE, 143 (approximately 1/3) were genes without GO terms; this likely explains the lack of significant GO enrichment results for this comparison, but this prevented us from determining which pathways are involved in the unique gene expression response of field larvae to a novel exposure to 35°C.

Of the significantly DE genes other than hsps, few were shared across all four novel heat shock comparisons (Fig. 2E). This was in part due to the notably unique gene expression response of field individuals to 35°C, but large numbers of genes (154–326) were unique to other comparisons as well, and only 69 genes were shared between any three sets. The absence of a conserved thermal stress response (other than hsps) could be caused by a lack of power to detect DE genes in field population comparisons as a result of more variation between biological replicates. However, this alone does not seem like an adequate explanation given the large numbers of genes we did detect as significantly DE in field individuals. One potential confounding effect of our experimental design is that ramp rates differed between 35 and 40°C treatments (as a result of keeping ramp duration consistent). This could affect gene expression patterns and potentially explains some of the differences we found between the novel heat shock treatments. Despite the small numbers of shared genes, GO enrichment results for individual comparisons were generally consistent with the findings of other studies about pathways involved in the thermal stress response in other organisms (López-Maury et al., 2008; Niskanen and Palvimo, 2017; Richter et al., 2010; Verger et al., 2003; Vihervaara et al., 2018). These include upregulation of signaling and transport proteins and pathways involved in sumoylation, ubiquitination and oxidative stress response, as well as downregulation of DNA replication, cell growth regulators and various metabolic processes. Therefore, there is significantly more conservation in heat shock responses on the level of pathways than in the response of specific genes. Because we took a tissue-specific approach for this study, the interpretations of our results are limited to genes expressed in fat body tissue. Fat body tissue is known to show robust expression of hsp genes as well many genes involved in a variety of important metabolic functions that may be altered by thermal stress (Arrese and Soulages, 2010; Fittinghoff and Riddiford, 1990), therefore we have no reason to think that overall stress-related gene expression patterns would differ substantially from those of other tissues. However, some GO terms from our analysis, such as lipid storage and transport, could be related to fat body-specific functions.

We found partial support for our prediction that the field population would show changes to gene expression reflective of a greater capacity to acclimate to temperature fluctuations compared with the lab population. Fluctuating rearing conditions of ±5°C affected the gene expression response to 35°C test temperatures for the field but not lab population of M. sexta. However, this trend did not hold for acclimation to 40°C and we have only very limited support for the more general prediction that acclimation capability is maintained as a direct consequence of selection due to greater environmental temperature variation. Gene expression patterns from this study parallel different physiological responses of field versus lab populations to diurnal fluctuations at 30°C mean temperature (Kingsolver et al., 2015; Kingsolver lab, unpublished). Compared with constant (30°C) temperatures, lab larvae reared at 30±5°C have longer development times and slower short-term growth rates at intermediate (20–35°C) temperatures. In contrast, for field larvae, these same fluctuating rearing conditions had no effect on development time and increased short-term growth rates at 20–35°C, indicative of an acclimation response, further supporting our conclusion that the difference in gene expression we saw in the field population at 35°C between fluctuating and constant temperatures may be an adaptive acclimation response. For lab caterpillars, the elevated gene expression stress response we found after ±10°C fluctuations during rearing could coincide with reduced short-term growth rates as there are likely metabolic costs associated with maintaining an elevated stress response over longer periods of time. However, given that acclimation (increase in short-term growth rate) has been observed for lab caterpillars in 25±10°C rearing conditions (Kingsolver et al., 2015), greater acclimation in the field population may not be a generalizable pattern for other fluctuating thermal conditions. Instead, the specific thermal history conditions required for inducing an acclimation response may be shifted for these two populations. Future work could evaluate plasticity for other fluctuation temperature regimes and determine whether acclimation at other temperatures results in similar gene expression patterns to the ones seen in this study.

Past temperature exposure can lead to complex interactions with current temperature in ectotherms; however, our current knowledge about these time-dependent effects for realistic temperature conditions is lacking. In this study, we used M. sexta larvae from laboratory and field populations to determine how diurnal temperature fluctuations affect patterns of gene expression both under control (non-stressful) conditions and for later thermal stress. We found evidence that repeated high temperature exposure can have diverse consequences depending on the specific thermal conditions and the populations tested. We predicted that the lab population should show greater DE of stress-related genes both during thermal stress and after repeated exposure compared with the more thermotolerant field population. Diurnal fluctuations of ±5°C did not alter long-term gene expression patterns for either population, but ±10°C fluctuations did appear to cause more stress for the lab relative to the field population. Novel temperature exposure resulted in upregulation of many of the same hsp genes in all cases, but the expression patterns of other genes were more variable across temperatures and populations. Notably, exposure to a 35°C novel test temperature resulted in strikingly unique gene expression profiles of field individuals, which could be related to the higher thermal tolerance observed in the field compared with the lab population. We also predicted that diurnal fluctuations could affect the gene expression response to later heat stress, with an acclimation-like response more likely for field larvae. Our results somewhat supported this expectation – changes to gene expression consistent with acclimation were observed for the field population at 35°C. However, greater acclimation capacity in the field is not necessarily a consistent pattern across other temperature conditions. These results highlight the complexity of making generalizable predictions about the physiological consequences of temperature fluctuations and the importance of continued exploration into the underlying mechanisms involved in time-dependent effects.

We thank current and former lab members J. Higgins, K. Augustine, M. Nielson, G. Legault, A. Parker, A. Deconinck and K. Malinski, and many others including J. Umbanhowar, L. Sligar, J. McGirr, E. Richards, M. St John and C. Chen, for helpful feedback at various stages of this research. Thanks also to instructors and participants at MDIBL Environmental Genomics Workshop 2017 for help with RNA-Seq analysis, and to two anonymous reviewers for constructive comments that improved the manuscript.