ABSTRACT

Most studies exploring molecular and physiological responses to temperature have focused on constant temperature treatments. To gain a better understanding of the impact of fluctuating temperatures, we investigated the effects of increased temperature variation on Phanaeus vindex dung beetles across levels of biological organization. Specifically, we hypothesized that increased temperature variation is energetically demanding. We predicted that thermal sensitivity of metabolic rate and energetic reserves would be reduced with increasing fluctuation. To test this, we examined the responses of dung beetles to constant (20°C), low fluctuation (20±5°C), or high fluctuation (20±12°C) temperature treatments using respirometry, assessment of energetic reserves and HPLC-MS-based metabolomics. We found no significant differences in metabolic rate or energetic reserves, suggesting increased fluctuations were not energetically demanding. To understand why there was no effect of increased amplitude of temperature fluctuation on energetics, we assembled and annotated a de novo transcriptome, finding non-overlapping transcriptomic and metabolomic responses of beetles exposed to different fluctuations. We found that 58 metabolites increased in abundance in both fluctuation treatments, but 15 only did so in response to high-amplitude fluctuations. We found that 120 transcripts were significantly upregulated following acclimation to any fluctuation, but 174 were upregulated only in beetles from the high-amplitude fluctuation treatment. Several differentially expressed transcripts were associated with post-translational modifications to histones that support a more open chromatin structure. Our results demonstrate that acclimation to different temperature fluctuations is distinct and may be supported by increasing transcriptional plasticity. Our results indicate for the first time that histone modifications may underlie rapid acclimation to temperature variation.

INTRODUCTION

Temperature variation impacts the physiology, ecology and evolution of organisms (Angilletta, 2009) and is a key driver of species distributions and diversity, particularly in small ectotherms like insects (Janzen, 1967). A vast literature on biochemical and physiological responses to temperature has been produced over the past century (Somero et al., 2017). In particular, a theory has developed that the maintenance of optimal fluidity of enzymes, membranes and other important macromolecules underlies biochemical adaptation to changing temperatures (Somero, 2004; Ernst et al., 2016). However, this work has largely focused on the effects of changes in constant temperature. In real-world conditions, temperatures fluctuate, sometimes by 20°C or more daily (Marshall and Sinclair, 2012), and thus it is critical to understand biochemical adaptations to fluctuating temperatures. How organisms maintain optimal membrane and enzyme fluidity during fluctuating temperatures is unclear.

In theory, increased temperature variation should greatly affect physiological performance, increasing metabolic costs and potentially decreasing energetic reserves. In insects, metabolic rate scales non-linearly with temperature, first increasing exponentially from the critical thermal minimum (CTmin) to an inflection point below the thermal optimum (Topt), decelerating from the inflection point to Topt and, finally, declining rapidly to the critical thermal maximum (CTmax) (Angilletta, 2009). Because of this non-linearity, energy expenditure in fluctuating temperatures is different from what would be predicted by mean temperature even when temperatures are fluctuating symmetrically around the mean (Ruel and Ayres, 1999; Williams et al., 2012). An individual's energetic expenditure is influenced by its thermal sensitivity of metabolism (i.e. the slope of the relationship between temperature and metabolic rate) (Williams et al., 2012), the amplitude of the thermal fluctuations (Pásztor et al., 2000; Ruel and Ayres, 1999) and the average temperature of the environment (Vasseur et al., 2014). In response to increases in temperature variation, organisms may decrease their thermal sensitivity of metabolism and thus energetic costs across a broader range of temperatures.

Shifts in the thermal sensitivity of metabolism due to increased temperature variation should reduce energetic reserves – the macromolecules metabolized for energy – with consequences for survival and fitness. Insects utilize two main sources of energy for metabolism: lipids and carbohydrates, stored primarily as triacylglycerols and glycogen, respectively (Arrese and Soulages, 2010). Lipids are essential for growth, immunity, reproduction, survival during diapause and maintenance during prolonged bouts without feeding (Cheon et al., 2006; Dettloff et al., 2001; Hahn and Denlinger, 2007; Meyer-Fernandes et al., 2000; Mullen and Goldsworthy, 2003; Park et al., 2006). In addition, lipids are a key component of eggs and are the primary energy source for developing insect embryos (Arrese and Soulages, 2010; Park et al., 2006). Glycogen fuels prolonged flight and prevents cellular damage during thermal stress or drought (Beenakkers, 1969; Kaufmann and Brown, 2008; Storey, 1997; Watanabe et al., 2002). Thermally variable conditions that deplete lipids and carbohydrates thus reduce energetic reserves available for immunity, foraging and reproduction. Because egg development requires a substantial investment of maternal energetic reserves, females with limited energetic reserves may either reduce lipid and carbohydrate allocation to eggs or forgo breeding altogether (Elkin and Reid, 2005; van Noordwijk and de Jong, 1986).

Insects may also respond to temperature changes through molecular and biochemical shifts. Transcriptomic and metabolomic approaches have uncovered the molecular mechanisms used to cope with extreme cold or heat stress (MacMillan et al., 2016; Runcie et al., 2012; Schoville et al., 2012; Stillman and Tagmount, 2009). However, few studies have considered the genes and pathways used by species under increasingly variable temperature conditions. A notable exception is work using Drosophila melanogaster. One study examining the responses to repeated cold stress (−0.5°C) found that the number of cold exposures caused distinct transcriptomic responses relative to single cold exposures, with upregulation of genes associated with egg production and immune and stress responses (Zhang et al., 2011). A second study in D. melanogaster examined transcriptomic responses to relatively mild thermal fluctuations of 8°C, and found that fluctuation involved distinct responses centered on phototransduction (Sørensen et al., 2016). However, little is known about the effects of thermal variation in other insect species. Comprehensive characterization of transcriptomic and metabolomic changes is needed to uncover the genes and pathways associated with responses to increased temperature variation and help us understand an organism's ability to respond to the types of temperature shifts being faced as a result of climate change.

We quantified physiological and molecular responses to increased temperature variation using the dung beetle Phanaeus vindex Macleay 1819 (Coleoptera: Scarabaeinae), a diurnal species that occurs in North America primarily east of the Continental Divide and south of 43°N latitude. The species has a broad thermal physiology (Sheldon and Tewksbury, 2014), is able to tolerate a wide range of habitats and ecological conditions, and provides many ecosystem services, including nutrient cycling, soil aeration and secondary seed dispersal (Nichols et al., 2008). We hypothesized that greater temperature variation is energetically demanding for P. vindex and would therefore (1) reduce the thermal sensitivity of metabolism, (2) reduce carbohydrate and lipid reserves, and (3) lead to greater differential gene and metabolite expression compared with smaller temperature variation. In particular, we were interested in the molecular mechanisms that underlie metabolic adjustments to increased temperature variation. By identifying physiological and molecular responses to increased temperature variation, we can better understand how temperature stress associated with climate change affects organisms across biological levels of organization.

MATERIALS AND METHODS

Sampling sites and experimental design

We collected 51 P. vindex beetles (45 females, 6 males) outside Knoxville, TN, USA, in April 2017. We housed them for 2 weeks at 20°C in plastic containers with autoclaved soil and fed them ad libitum with autoclaved cow dung. We then transported the beetles to the University of Oklahoma, Norman, OK, USA, for laboratory tests. We followed all applicable national and/or institutional guidelines for the care and use of animals.

At the University of Oklahoma, we placed all beetles in incubators (Panasonic MIR-154) at 20°C for a 3 day laboratory acclimation period prior to temperature treatments. After this period, 20 female P. vindex were randomly put in one of three temperature treatments for 6 days: constant (20°C), low fluctuation (20±5°C) or high fluctuation (20±12°C). The treatment temperatures fluctuated in a 24 h diel pattern with temperature set points changing every 2 h. Temperatures ramped from the previous temperature set point to the next point over the course of 120 min and passed the treatment's mean temperature twice in a 24 h period. We also placed 3 male P. vindex in both the constant and high-fluctuation treatments to see whether metabolomic responses varied by sex. The temperatures in the high-fluctuation treatment approach the highest and lowest temperatures the beetles would experience during their active time of year in Tennessee (April–September). We expected these temperatures to be stressful as the extremes are occurring on a diurnal cycle rather than across several months. In contrast, the low-fluctuation treatment is within the range of diurnal cycling during the active season. The constant temperature treatment represents the average temperature during the active season in Tennessee.

After 6 days in the temperature treatments, 5 females from each treatment were haphazardly selected and used for metabolic rate measurements, 5 were used for carbohydrate and lipid analyses, and the remaining 5 were used for both transcriptomics and metabolomics. The 3 males from the high-fluctuation and control treatments were used in metabolomic work for comparison with the females.

Measuring and analyzing metabolic rate

To examine the energetic demands of increased temperature variation, we measured the thermal sensitivity of metabolism (i.e. metabolic rate across a range of temperatures) for 5 females from each treatment. For metabolic trials, we used the Q-Box RP1LP Low Range Respiration Package from Qubit Systems and flow-through respirometry with a flow rate of 125–150 ml min−1. For each female, we measured the amount of CO2 produced at a series of increasing temperatures (at approximately 10, 15, 20 and 25°C). We used outdoor air as our baseline air source and scrubbed it of H2O using Drierite™ Drying Desiccants before experimental trials. For trials, we randomly selected a single female from one of the three temperature treatments, weighed the beetle, and placed it in a 15 ml glass chamber with tubing on either end for continuous air flow. The glass chamber was inside a darkened Styrofoam cooler with a viewing window, to minimize activity of the females but still allow us to see the beetle. We held the beetle in the chamber at 20°C for 10 min, then dropped the temperature by 1°C per minute to 10°C, and then held the beetles for another 10 min at 10°C prior to the start of the trial. At the start of a trial, we recorded a baseline reading of CO2 and then recorded the amount of CO2 produced by the beetle over approximately 5 min. We then took another baseline reading of CO2 and began raising the temperature in the chamber at 0.5°C per minute up to 15°C (i.e. the next test temperature). Once at 15°C, we again recorded baseline CO2, the CO2 produced by the beetle over approximately 5 min, and a final baseline CO2 reading. We continued this to obtain CO2 measurements at 20 and 25°C. We used the baseline readings to account for any drift. Because we held the beetles in a dark chamber and our test temperatures (10–25°C) were not stressful, the beetles were inactive during the metabolic trials. On the rare occasions when we saw activity, we waited for the beetle to stop moving before recording metabolic rate.

For data analysis of metabolic rate, we first converted CO2 measurements from parts per million concentration to μl min−1 using an equation suited to our experimental set-up:
formula
(1)
where FiCO2 is the incurrent fractional CO2 concentration and FeCO2 is the excurrent fractional CO2 concentration. For all calculations, we used an assumed respiratory quotient (RQ) of 0.85 as diet was not strictly controlled. We log-transformed final CO2 values to meet the assumptions of normality and homoscedasticity. We used a general linear mixed model (Proc Mixed, SAS v 9.4) that included temperature treatment, metabolic trial temperature and their interactions. We included beetle mass taken just prior to the respirometry trials as a covariate. We used a repeated statement specifying individual beetle as the subject to control for the non-independence of measurements from the same individual, and specified an unstructured covariance matrix. We fitted models using maximum likelihood estimation and used the Akaike information criterion (AIC) to choose the best-fit model.

Energetic reserves and analyses

For carbohydrate and lipid sampling, we weighed 5 female beetles from each treatment, then placed them individually in a 5 ml glass vial in liquid nitrogen to freeze thoroughly and stored them at −80°C until processing. For processing, we added 5 ml of homogenization mixture (5 parts Tween 20 to 995 parts RO water for a 0.05% Tween 20 solution) to each glass vial. We ground the beetles using a hand-held homogenizer, then made 2 ml aliquots of beetle homogenate and stored them at −80°C. We assayed the polar fraction for glycogen content using enzyme-linked spectrophotometric assay (Marshall and Sinclair, 2015). We measured protein content from the polar fraction using the BCA assay as a covariate (Tennessen et al., 2014). We measured the lipid content gravimetrically by first extracting lipids using Folch reagent (2:1 chloroform:methanol), then drying the lipid fraction in a Vacufuge speedvac, and then weighing the resulting lipid. We compared energetic content among experimental groups using ANOVA with fluctuation group as the only predictor. We also calculated effect size using Cohen's d, which reports the standardized difference between means.

Transcriptomics of beetles in constant and fluctuating treatments

For transcriptomic and metabolomic data, we took 5 females from each treatment and 3 males from the constant and high-fluctuation treatments and placed them individually in 5 ml glass vials in liquid nitrogen to freeze thoroughly, then stored them at −80°C until processing. We ground whole beetles in liquid nitrogen using a handheld homogenizer. From the 15 females, we extracted total RNA from half of each homogenate using a TRIzol protocol followed by cleanup from an RNAeasy Minikit (Qiagen). The other half of the homogenate from the 15 females and half of the homogenate from the 6 males was used for metabolomics (see below).

Sequencing libraries were prepared and sequenced by the Oklahoma Medical Research Foundation Genomics Facility (Oklahoma City, OK, USA). Libraries were prepared using the NEB mRNA stranded library prep kits and sequenced on a HiSeq 2000 producing 150 bp long paired-end reads. We used Trimmomatic quality control software (version 0.36.3; Bolger et al., 2014) to clip adapter sequences and remove low-quality reads (quality scores <26) and reads <100 bp. We assembled a transcriptome using Trinity (Haas et al., 2013), then clustered redundant sequences using CD-HIT-EST with a sequence identity cut-off of 0.95 (Fu et al., 2012). We aligned these sequences to the NCBI NR database using Diamond (Buchfink et al., 2015), then filtered out all reads that did not map to metazoans (viruses, algae, bacteria, etc.) using MEGAN6 (Huson et al., 2016) to produce a de novo transcriptome for P. vindex. Using Blast2Go software (BioBam, Valencia, Spain) we annotated transcripts with Interpro and GO classifications. We also annotated each sequence with the nearest KEGG (Kyoto Encyclopedia of Genes and Genomes) orthology term using the KEGG Automatic Annotation Server (Moriya et al., 2007).

To quantify gene expression, we mapped all quality-controlled reads back to the remaining transcripts using HISAT2 (Kim et al., 2015). We assembled transcripts using StringTie (Pertea et al., 2015), quantified expression using FeatureCounts (Liao et al., 2014), and then tested for differential expression using DESeq2 (Love et al., 2014). Using the list of significantly differentially expressed transcripts, we grouped transcripts by their upregulation or downregulation relative to beetles in the control treatment using Vennplex (Cai et al., 2013). Using Blast2Go, we examined enrichment of Gene Ontology (GO) categories in each group produced by Vennplex, then simplified the list of terms using ReviGO (Supek et al., 2011). We tested for KEGG pathway enrichment using the Gage and Pathview (Luo and Brouwer, 2013) packages in R. We also used MetaboAnalyst (Chong et al., 2019) to conduct sparse partial least squares discriminant analysis (sPLS-DA; Lê Cao et al., 2011) with 10 variables per component following mean-centering and dividing by the standard deviation of each variable on the most variable 5000 transcripts. To understand whether any function responded to amplitude in a dose-dependent way, we conducted a PatternHunter analysis in MetaboAnalyst to highlight transcripts that increased in relative expression with increasing amplitude of variation.

Metabolomics of beetles in constant and fluctuating treatments

We placed the other half of the homogenate generated from the 15 female and 6 male beetles in 3:1:1 chloroform:methanol:water (6.35 µl per mg tissue) and then prepared samples as in MacMillan et al. (2016). Briefly, homogenates were placed at −20°C for 60 min to precipitate protein, then centrifuged at 16000 RCF to remove debris. We then stored the polar supernatants at −80°C.

Samples were transferred to the University of Oklahoma Mass Spectrometry Facility where 5 µl of each sample was injected in technical duplicate on an Agilent 1290 Infinity UPLC system coupled to an Agilent 6545 quadrupole time-of-flight mass spectrometer (QTOF-MS; Agilent, Santa Clara, CA, USA) equipped with a dual-spray ESI source. An Acquity UPLC HSS C18 SB column was used for separation, with solvent gradient from 95% water (plus 0.1% formic acid) to 100% acetonitrile (plus 0.1% formic acid) at a flow rate of 0.4 ml min−1. Raw MS data were then processed as in Bonifay et al. (2017). XCMS Centwave (Tautenhahn et al., 2008) was used to detect peaks, and mzMatch (Scheltema et al., 2011) was used for peak alignment and filtering. Detected features were then matched against metabolites in the KEGG (Kanehisa and Goto, 2000), LIPID MAPS (O'Donnell et al., 2019) and HMDB (Wishart et al., 2007) databases with the criterion that the experimentally determined exact mass should be within 6 ppm of the theoretical mass of the metabolite.

We conducted all downstream analyses on the average of the two technical replicates. We filtered out any metabolites with a within-group relative standard deviation of >0.5. Using MetaboAnalyst, we then filtered non-informative (relatively constant) metabolites based on interquartile range. We scaled all metabolites by mean-centering and dividing by the standard deviation of each variable to follow assumptions of normality. Using MetaboAnalyst, we examined significant variation in metabolites using ANOVA followed by false discovery rate correction and Fisher's LSD post hoc tests. We also conducted sPLS-DA (Lê Cao et al., 2011) with 10 variables per component to examine the most informative metabolites. Because we found no differences in male and female responses based on metabolomic data, we focused our results on female responses (see Fig. S1 for comparison of males and females).

RESULTS

Metabolic rate and energetic reserves

As expected, metabolic rate increased as test temperature increased (F1,11=8.75, P=0.01) but did not vary with female body size (F1,11=3.82, P=0.08) (Fig. 1). Temperature fluctuation treatments did not affect the thermal sensitivity of metabolism (F2,11=0.61, P=0.56) or overall metabolic rate (F2,11=1.53, P=0.26).

Fig. 1.

Metabolic rate of Phanaeus vindex following acclimation to temperature fluctuation or control treatments. Metabolic rate was measured as CO2 production rate across a range of temperatures in P. vindex dung beetles (n=15 females) after acclimation for 6 days to the high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) or control (20°C; n=5) temperature treatments. Plotted are least squares means (±s.e.) of log-transformed data. Temperature treatments did not affect the thermal sensitivity of metabolism (F2,11=0.61, P=0.56) or overall metabolic rate (F2,11=1.53, P=0.26).

Fig. 1.

Metabolic rate of Phanaeus vindex following acclimation to temperature fluctuation or control treatments. Metabolic rate was measured as CO2 production rate across a range of temperatures in P. vindex dung beetles (n=15 females) after acclimation for 6 days to the high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) or control (20°C; n=5) temperature treatments. Plotted are least squares means (±s.e.) of log-transformed data. Temperature treatments did not affect the thermal sensitivity of metabolism (F2,11=0.61, P=0.56) or overall metabolic rate (F2,11=1.53, P=0.26).

For carbohydrate and lipid analyses, we found no effect of either total body mass or soluble protein content on glycogen or lipid content (P>0.3 in all cases), so these terms were dropped from further analyses. Among the three treatments, we found no significant differences in glycogen (F2,12=1.17, P=0.344) or lipid (F2,12=2.7, P=0.103) content of females. However, there was a large effect size of fluctuation amplitude on both glycogen (Cohen's d=0.76 for low amplitude; d=0.92 for high amplitude; Fig. 2A) and lipid content (Cohen's d=0.38 for low amplitude; d=1.82 for high amplitude; Fig. 2B).

Fig. 2.

Energetic reserves of P. vindex following acclimation to temperature fluctuation or control treatments. Total glycogen content (A) and total lipid content (B) of dung beetles (n=15 females) acclimated to high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) or control (20°C; n=5) temperature treatments for 6 days. The thick horizontal line indicates the group median and the limits of the solid box indicate the interquartile range.

Fig. 2.

Energetic reserves of P. vindex following acclimation to temperature fluctuation or control treatments. Total glycogen content (A) and total lipid content (B) of dung beetles (n=15 females) acclimated to high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) or control (20°C; n=5) temperature treatments for 6 days. The thick horizontal line indicates the group median and the limits of the solid box indicate the interquartile range.

Transcriptomics

Based on data from 15 females, a total of 332.5 million paired-end reads were produced from sequencing, with an average of 22.17 million paired-end reads per sample. Of these, an average of 13.15 million reads per sample passed our stringent quality test. These were then assembled into a transcriptome with 161,002 sequences, which was reduced to 106,200 sequences following CD-HIT-EST clustering. Filtering out all transcripts that were from non-metazoan sources reduced this number to 41,134 transcripts, with an N50 of 2052. Mapping rates to this transcriptome were generally high with an average of 72% of filtered reads mapping at least once to the de novo transcriptome.

Acclimation to both high- and low-amplitude temperature fluctuations caused significant differential transcript expression relative to the control temperature treatment. A total of 394 transcripts were upregulated following low-amplitude fluctuation acclimation, while 294 transcripts were upregulated following acclimation to high-amplitude fluctuations (Fig. 3). Out of these, a total of 120 were upregulated following any temperature fluctuation. Relatively fewer were downregulated: only 159 transcripts were downregulated following low-amplitude temperature fluctuations, and 112 following high-amplitude fluctuations, with only 11 of these downregulated in both groups. The transcripts most upregulated in both fluctuation groups included kelch 26 amino acid transporter, an AVT1A isoform and a transcript containing an ankyrin repeat domain (Table S1).

Fig. 3.

Differential expression of transcripts and metabolites in P. vindex following acclimation to temperature fluctuation treatments relative to a control treatment. Venn diagrams illustrating the number of differentially expressed transcripts (top) and metabolites (bottom) in dung beetles (n=15 females) following 6 days of acclimation to either low-amplitude (20±5°C; n=5) or high-amplitude (20±12°C; n=5) temperature fluctuations relative to control beetles maintained at 20°C (n=5). Up-arrows indicate significantly upregulated transcripts or metabolites, down-arrows indicate significantly downregulated transcripts or metabolites, and a tilde indicates contra-regulated metabolites which are upregulated in one treatment group but downregulated in another.

Fig. 3.

Differential expression of transcripts and metabolites in P. vindex following acclimation to temperature fluctuation treatments relative to a control treatment. Venn diagrams illustrating the number of differentially expressed transcripts (top) and metabolites (bottom) in dung beetles (n=15 females) following 6 days of acclimation to either low-amplitude (20±5°C; n=5) or high-amplitude (20±12°C; n=5) temperature fluctuations relative to control beetles maintained at 20°C (n=5). Up-arrows indicate significantly upregulated transcripts or metabolites, down-arrows indicate significantly downregulated transcripts or metabolites, and a tilde indicates contra-regulated metabolites which are upregulated in one treatment group but downregulated in another.

To examine functional enrichment of these differentially expressed transcripts, we started by conducting Fisher's exact tests followed by multiple testing correction on GO categories in Blast2Go. We found that there were no significantly enriched GO terms among the upregulated transcripts, but among the downregulated transcripts, there was a large but non-overlapping number of enriched GO terms in both acclimation groups following ReviGO simplification (Table 1). In beetles that received low-amplitude temperature fluctuations, we found that GO terms such as ‘glucosamine-containing compound metabolic process’, ‘amino sugar metabolic process’ and ‘protein dephosphorylation’ were significantly enriched among downregulated transcripts. In beetles that received high-amplitude temperature fluctuations, we found that GO terms such as ‘translation’, ‘cellular amide metabolic process’ and ‘gene expression’ were significantly enriched among downregulated transcripts. We observed a similar pattern in KEGG pathway enrichment (Table 2). In beetles that received low-amplitude fluctuations, we did not find any significantly upregulated pathways, but we found that ‘bacterial invasion of epithelial cells’ and ‘Fc gamma R-mediated phagocytosis’ were downregulated. In contrast, in beetles that received high-amplitude fluctuation, we discovered KEGG pathways such as ‘DNA replication’ and ‘cell cycle’ were upregulated, while ‘ribosome’ and ‘tyrosine metabolism’ were downregulated. None of the significantly enriched pathways were common to both amplitude fluctuation treatment groups.

Table 1.

Significantly enriched gene ontology (GO) terms in transcripts downregulated following acclimation for 6 days to fluctuating temperature treatments (see Fig. 3) relative to the control treatment in the dung beetle Phanaeus vindex

Significantly enriched gene ontology (GO) terms in transcripts downregulated following acclimation for 6 days to fluctuating temperature treatments (see Fig. 3) relative to the control treatment in the dung beetle Phanaeus vindex
Significantly enriched gene ontology (GO) terms in transcripts downregulated following acclimation for 6 days to fluctuating temperature treatments (see Fig. 3) relative to the control treatment in the dung beetle Phanaeus vindex
Table 2.

Significantly enriched KEGG terms in P. vindex following acclimation to fluctuating temperature treatments relative to the control

Significantly enriched KEGG terms in P. vindex following acclimation to fluctuating temperature treatments relative to the control
Significantly enriched KEGG terms in P. vindex following acclimation to fluctuating temperature treatments relative to the control

sPLS-DA was able to easily separate the samples by temperature treatment (Fig. 4) and so did hierarchical clustering (Fig. 5). Component 1 (18.2% of the variation explained) separated the control and low-amplitude fluctuation groups from the high-amplitude fluctuation group, while component 2 (12.2% of the variation explained) separated the control group from the low-amplitude fluctuation group, with the high-amplitude fluctuation group in the middle of the other two groups. Transcripts that loaded highly on component 1 (Table 3) included farnesyl pyrophosphate synthase-like (part of the juvenile hormone precursor pathway), a transcription factor and a ribosomal transcript. In contrast, transcripts that loaded highly on component 2 included a catenin delta-2 isoform, an integrin-like transcript and another currently unknown transcription factor.

Fig. 4.

Sparse partial least squares discriminant analysis (sPLS-DA) of transcripts (left) and metabolites (right). sPLS-DA easily differentiated among dung beetles (n=15 females) exposed to the high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) or control (20°C; n=5) temperature treatments.

Fig. 4.

Sparse partial least squares discriminant analysis (sPLS-DA) of transcripts (left) and metabolites (right). sPLS-DA easily differentiated among dung beetles (n=15 females) exposed to the high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) or control (20°C; n=5) temperature treatments.

Fig. 5.

Heat maps of transcript abundance (left) and metabolite abundance (right). The heat maps represent the top 50 differentially expressed transcripts and metabolites, each normalized by row, from dung beetles (n=15 females) in the high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) and control (20°C; n=5) temperature treatments. Clustering was performed with a Euclidean distance measure and by the Ward algorithm in MetaboAnalyst.

Fig. 5.

Heat maps of transcript abundance (left) and metabolite abundance (right). The heat maps represent the top 50 differentially expressed transcripts and metabolites, each normalized by row, from dung beetles (n=15 females) in the high-fluctuation (20±12°C; n=5), low-fluctuation (20±5°C; n=5) and control (20°C; n=5) temperature treatments. Clustering was performed with a Euclidean distance measure and by the Ward algorithm in MetaboAnalyst.

Table 3.

Top five component loadings from a sparse partial least squares analysis (sPLS-DA) of transcripts (Fig. 4) in P. vindex following acclimation to fluctuating temperature treatments relative to the control

Top five component loadings from a sparse partial least squares analysis (sPLS-DA) of transcripts (Fig. 4) in P. vindex following acclimation to fluctuating temperature treatments relative to the control
Top five component loadings from a sparse partial least squares analysis (sPLS-DA) of transcripts (Fig. 4) in P. vindex following acclimation to fluctuating temperature treatments relative to the control

Using PatternHunter to examine transcripts that were differentially regulated according to the amplitude of temperature fluctuation, we discovered many transcripts related to histone modifications, but these transcripts did not necessarily have informative GO terms (e.g. ‘protein binding’). After examining all significantly differentially expressed transcripts in either the low- or high-amplitude fluctuation groups that had the term ‘histone’ in the transcript name or the associated GO term, we found a dose dependency in the expression of these transcripts. Specifically, a greater amplitude of temperature fluctuation corresponded with greater differential expression of almost every single histone-related transcript (Fig. 6).

Fig. 6.

All histone-related terms in differentially expressed transcripts in P. vindex following acclimation to temperature fluctuation. Data are from dung beetles (n=15 females) in the high-fluctuation (20±12°C; n=5) and low-fluctuation (20±5°C; n=5) treatments relative to control (20°C; n=5) temperature treatment. The term ‘histone’ was found in either the transcript description or the assigned GO terms. Values represent doublings of each sample relative to the mean of the control. *Transcripts significantly differentially expressed relative to control in beetles exposed to high-amplitude temperature fluctuations (20±12°C); transcripts significantly differentially expressed relative to control in dung beetles exposed to low-amplitude temperature fluctuations (20±5°C).

Fig. 6.

All histone-related terms in differentially expressed transcripts in P. vindex following acclimation to temperature fluctuation. Data are from dung beetles (n=15 females) in the high-fluctuation (20±12°C; n=5) and low-fluctuation (20±5°C; n=5) treatments relative to control (20°C; n=5) temperature treatment. The term ‘histone’ was found in either the transcript description or the assigned GO terms. Values represent doublings of each sample relative to the mean of the control. *Transcripts significantly differentially expressed relative to control in beetles exposed to high-amplitude temperature fluctuations (20±12°C); transcripts significantly differentially expressed relative to control in dung beetles exposed to low-amplitude temperature fluctuations (20±5°C).

Metabolomics

HPLC-based metabolomics detected a total of 1752 peaks across our 15 female and 6 male samples. Of these, a total of 848 were assigned to a putative metabolite, and only 462 metabolites were retained following filtering for repeatability. Using ANOVA followed by post hoc tests then false discovery rate correction, we found a total of 93 metabolites whose abundance shifted compared with the control group following temperature fluctuation treatments, with the majority (58 metabolites) increasing in abundance regardless of amplitude (Fig. 3; Table S2). The patterns in male beetles strongly followed those in female beetles (see Fig. S1).

A sPLS-DA found that component 1 (22% of total variation explained) separated experimental groups along a gradient of increasing amplitude of variation (Fig. 4). Metabolites that loaded strongly on this component included 5S-HETE diendoperoxide (a byproduct of lipid oxidation), a few tripeptides, and 3-methoxy-4-hydroxyphenylglycol, which is a byproduct of noradrenaline (norepinephrine) degradation (Table 4).

Table 4.

Top five component loadings from a sPLS-DA of metabolites (Fig. 4) in P. vindex following acclimation to fluctuating temperature treatments relative to the control

Top five component loadings from a sPLS-DA of metabolites (Fig. 4) in P. vindex following acclimation to fluctuating temperature treatments relative to the control
Top five component loadings from a sPLS-DA of metabolites (Fig. 4) in P. vindex following acclimation to fluctuating temperature treatments relative to the control

DISCUSSION

Here, we demonstrate that temperature fluctuation elicits distinct responses at the transcriptomic and metabolomic levels compared with constant temperature. Interestingly, we did not observe the energetic costs we would have predicted on the whole-organism level, but it is possible that feeding rates, which we did not measure, may have compensated for increases in energetic demands in the fluctuating treatments. We did, however, find that the nature of the underlying transcriptomic and metabolomic mechanisms depend significantly on the amplitude of temperature fluctuation. These results suggest that simple models using laboratory acclimation to constant temperatures are unlikely to accurately predict organismal responses to more variable thermal environments.

Experimental studies of the impact of fluctuating temperatures on whole-organism metabolic rates remain relatively few given the complexity of potential experimental designs (Colinet et al., 2015). One study that used a naturalistic temperature regime found that a greater amplitude of temperature variation caused a significant depression in the thermal sensitivity of metabolism in caterpillars of the hesperid Erynnis propertius (Williams et al., 2012). By contrast, a recent review of the effects of fluctuating temperature on metabolic rate in ectothermic vertebrates found that, generally, resting metabolic rate increases with increased temperature fluctuation (Morash et al., 2018). Given that P. vindex is a temperate beetle that experiences high seasonality in temperature, it has a relatively broad thermal physiology (Sheldon and Tewksbury, 2014) and is expected to have a high degree of thermal plasticity to respond to increased temperature variation (Sheldon et al., 2018). We thus expected that acclimation to a greater amplitude of fluctuation would cause P. vindex beetles to decrease their thermal sensitivity of metabolism relative to acclimation to lower amplitude of fluctuation (Williams et al., 2012), but this was not the case. It is possible, however, that our sample size (5 beetles per treatment) and the large variation inherent in metabolic rate data hindered our ability to detect differences. In addition, we were unable to account for the age of the beetles, which has been shown to influence metabolic rate in other insect species (Hack, 1997; Piiroinen et al., 2010). However, we trapped dung beetles at the same time and randomly assigned them to treatments, thus reducing the potential for age to confound results. Importantly, we do not know the thermal history of the beetles (e.g. developmental temperatures), which could potentially impact metabolic rate. However, we acclimated adult beetles to the same conditions prior to experimental treatments. In addition, previous work with a different dung beetle species, Onthophagus taurus, showed that the temperatures that individuals experienced during the adult stage, but not egg to pupal stages, affected adult thermal plasticity of metabolism (Carter and Sheldon, 2020). Similarly, Gray (2013) found that thermal plasticity of metabolism in mosquitoes was affected by adult but not larval treatment temperatures. This suggests that although we could not account for thermal history of the beetles, the metabolic rate of adults may be driven largely by the temperatures they experienced in experimental trials.

The two different amplitudes of temperature fluctuation caused distinct transcriptomic responses both in individual transcripts and in GO enrichment (Tables 13). We found that sPLS-DA was able to separate beetles in both temperature fluctuation treatments from those in the control treatment. In addition, there was no component in which beetles from the temperature fluctuation treatments were rank ordered by amplitude of temperature fluctuation (Fig. 4), suggesting that increasing amplitude of fluctuation is not an additive effect at the transcriptomic level. Instead, component 1 describes the separation of beetles from the high-fluctuation treatment from beetles in the control and low-fluctuation treatments, and component 2 describes the separation of beetles in the low-fluctuation treatment from those in the control treatment (Fig. 4).

The transcripts that loaded highly on component 1 include several transcription factors as well as the enzyme farnesyl pyrophosphate synthase-like (Table 3). This enzyme catalyzes the reaction leading to the production of farnesyl pyrophosphate, which is the compound that straddles the boundary between the melavonate and juvenile hormone pathways (Noriega, 2014). Indeed, several isoforms of this enzyme appear to be upregulated in the high-amplitude fluctuation group, but not the low-amplitude fluctuation group (Tables S3 and S4). Farnesoic acid is the immediate precursor to the active juvenile hormone III, which is the key juvenile hormone in beetles (Shinoda and Itoyama, 2003). This suggests that juvenile hormone is involved in the response to high-amplitude temperature fluctuations but not in the response to low-amplitude temperature fluctuations.

Based on the GO and KEGG term functional enrichment, it is clear that high-amplitude fluctuations also significantly downregulated transcripts associated with translation (Tables 1 and 2). In the low-amplitude group, there were no significantly upregulated pathways using either GO terms or KEGG pathways, although we observed several pathways that were downregulated, including cell-to-cell adhesions, carbohydrate metabolic processes and protein dephosphorylation (Tables 1 and 2).

Though we found significant differentiation between our three experimental groups at the transcriptomics level, there were a total of 131 transcripts differentially expressed relative to controls in the same direction in response to any temperature fluctuation (Fig. 3). These included several amino acid transporters (Table S1), which may be part of a general stress response that has been identified in plants (Wan et al., 2017). Taken together, our data suggest that any amplitude of temperature fluctuation induces responses at the transcriptomic level that are unique from responses to constant temperatures.

Much of the literature exploring the effects of temperature fluctuation on insects have used a ‘fluctuating thermal regime’ (FTR), in which a long-term cold stress is interrupted by a brief rewarming to a benign temperature (reviewed by Colinet et al., 2015). These experiments have generally found that transcriptomic responses to FTR include multiple repair mechanisms, such as scavenging for reactive oxygen species, rebuilding ion homeostasis and mobilization of compatible solutes that function as osmoprotectants and cryoprotectants. We did not see any evidence for these processes with our design, which is more focused on understanding the effects of temperature fluctuation rather than on the mechanisms of repair following cold exposure.

Given that high- and low-amplitude fluctuations seemed to induce separate transcriptional responses, we were interested in exploring whether any functions responded to amplitude of fluctuation in a dose-dependent way. We found that several genes for histone protein modifications were upregulated in both fluctuation treatments but this tended to be higher in females from the high-fluctuation treatment (Fig. 6). Histone modifications can increase or decrease the packing of chromatin, thereby increasing or decreasing the availability of RNA polymerase II binding sites, ultimately modulating gene expression (Bannister and Kouzarides, 2011). These upregulated histone modification enzymes all appear to support greater chromatin accessibility (particularly the histone acetyltransferases). Increased expression of enzymes that modify histones to increase chromatin packing, and therefore decrease gene expression, has been observed in response to anoxia in the turtle Trachemys scripta elegans (Wijenayake et al., 2018) and appears to be a driver of diapause in the flesh fly Sarcophaga crassipalpis (Reynolds et al., 2016). Transcriptional studies of fluctuating temperatures are few, but in the annual killifish Austrofundulus limneaus, acclimation to fluctuating temperatures appears to cyclically regulate the high mobility group b1 protein, which binds to histones and thereby modulates gene expression (Podrabsky and Somero, 2004). Taken together, it appears that histone modification may be a relatively common way to acclimate to increased temperature variation. We thus speculate that modulating gene expression ‘on the fly’ through thermal cycles may be an important mechanism for survival in the naturally occurring thermal cycles in the wild.

In the metabolome, we observed similar patterns with a clear separation between the beetles that experienced high- and low-amplitude fluctuations (Figs 4 and 5). However, we did observe greater rank ordering of metabolites compared than for the transcriptome, with component 1 of the sPLS-DA (22% of the variation explained) separating the groups by amplitude and component 2 (13.4% of variation explained), clearly separating the high- and low-amplitude groups. The top weighted metabolite on component 1, 5S-HETE di-endoperoxide, is the product of degradation of 5S-HETE by COX-2 (Griesser et al., 2009). This metabolite is notable because 5S-HETE is a byproduct of the oxidative degradation of arachidonic acid, a common fatty acid in insects. This suggests that the observed increase in 5S-HETE di-endoperoxide with greater temperature fluctuation may be related to greater oxidative damage. There was also a large series of compounds (58) that increased in abundance in response to both high and low fluctuations relative to the control treatment (Fig. 3). The compound that increased the most relative to the control was a cardiolipin (Table S5). These lipids make up about 20% of the lipids in the inner mitochondrial membrane and are essential for abiotic stress responses as well as the formation of reactive oxygen species (de Paepe et al., 2014). While we observed greater overlap between high- and low-amplitude treatment groups in the metabolome than in the transcriptome, there was still clear functional separation between the groups, with greater abundance increases of small phospholipid compounds in the high-amplitude treatment group (Table S6).

Taken together, our data demonstrate that relatively short acclimation to even low-amplitude fluctuating temperatures induces a distinct physiological response in dung beetles relative to a constant temperature treatment. This suggests that models built using static laboratory acclimation temperatures do not represent insect responses in the wild. In addition, the amplitude of temperature fluctuation may induce acclimation through differing physiological pathways, with only higher amplitudes of fluctuation leading to suppression of protein translational activity. Finally, our data on transcripts associated with histone modifications suggest a novel mechanism for acclimation to temperature fluctuations: in the face of a highly variable environment, increased transcriptional plasticity may be beneficial for adjusting transcription ‘on the fly’.

Acknowledgements

We are grateful to the late Dr Rosemary Knapp who was an inspiring mentor for K.S.S. and K.E.M. and who provided equipment for our metabolic work. We thank Brittney Williams for laboratory assistance and Vincent Bonifay for helpful discussions.

Footnotes

Author contributions

Conceptualization: K.S.S., K.E.M.; Methodology: K.S.S., M.P., K.E.M.; Validation: K.S.S., K.E.M.; Formal analysis: K.S.S., A.W.C., K.E.M.; Investigation: K.S.S., M.P., K.E.M.; Resources: K.S.S., K.E.M.; Data curation: K.S.S., K.E.M.; Writing - original draft: K.S.S., K.E.M., A.W.C.; Writing - review & editing: K.S.S., M.P., A.W.C., K.E.M.; Visualization: K.S.S., A.W.C., K.E.M.; Supervision: K.S.S., K.E.M.; Project administration: K.S.S., K.E.M.; Funding acquisition: K.S.S., K.E.M.

Funding

This research was supported by the US National Science Foundation Division of Integrative Organismal Systems [IOS-1930829 to K.S.S. and K.E.M.].

Data availability

Raw metabolomics data, metabolic rate data, and carbohydrate and lipid data are available from Open Science Framework (doi:10.17605/OSF.IO/H34J2, https://osf.io/h34j2/). Raw transcriptomic read data have been archived in the NCBI database under Bioproject accession PRJNA667000.

References

Angilletta
,
M. J.
(
2009
).
Thermal Adaptation: A Theoretical and Empirical Synthesis
.
OUP
.
Arrese
,
E. L.
and
Soulages
,
J. L.
(
2010
).
Insect fat body: energy, metabolism, and regulation
.
Ann. Rev. Entomol.
55
,
207
-
225
.
Bannister
,
A. J.
and
Kouzarides
,
T.
(
2011
).
Regulation of chromatin by histone modifications
.
Cell Res.
,
21
,
381
-
395
.
Beenakkers
,
A. M. T.
(
1969
).
Carbohydrate and fat as a fuel for insect flight. A comparative study
.
J. Insect Physiol.
15
,
353
-
361
.
Bolger
,
A. M.
,
Lohse
,
M.
and
Usadel
,
B.
(
2014
).
Trimmomatic: a flexible trimmer for Illumina sequence data
.
Bioinformatics (Oxford, England)
30
,
2114
-
2120
.
Bonifay
,
V.
,
Wawrik
,
B.
,
Sunner
,
J.
,
Snodgrass
,
E. C.
,
Aydin
,
E.
,
Duncan
,
K. E.
,
Callaghan
,
A. V.
,
Oldham
,
A.
,
Liengen
,
T.
and
Beech
,
I.
(
2017
).
Metabolomic and metagenomic analysis of two crude oil production pipelines experiencing differential rates of corrosion
.
Front. Microbiol.
8
,
99
.
Buchfink
,
B.
,
Xie
,
C.
and
Huson
,
D. H.
(
2015
).
Fast and sensitive protein alignment using DIAMOND
.
Nat. Methods
12
,
59
-
60
.
Cai
,
H.
,
Chen
,
H.
,
Yi
,
T.
,
Daimon
,
C. M.
,
Boyle
,
J. P.
,
Peers
,
C.
,
Maudsley
,
S.
and
Martin
,
B.
(
2013
).
VennPlex--a novel Venn diagram program for comparing and visualizing datasets with differentially regulated datapoints
.
PLoS ONE
8
,
e53388
.
Carter
,
A. W.
and
Sheldon
,
K. S.
(
2020
).
Life stages of the dung beetle Onthophagus taurus differ in plasticity to temperature fluctuations and uniquely contribute to adult phenotype
.
J. Exp. Biol.
223
,
jeb227884
.
Cheon
,
H.-M.
,
Shin
,
S. W.
,
Bian
,
G.
,
Park
,
J.-H.
and
Raikhel
,
A. S.
(
2006
).
Regulation of lipid metabolism genes, lipid carrier protein lipophorin, and its receptor during immune challenge in the mosquito Aedes aegypti
.
J. Biol. Chem.
281
,
8426
-
8435
.
Chong
,
J.
,
Yamamoto
,
M.
and
Xia
,
J.
(
2019
).
MetaboAnalystR 2.0: from raw spectra to biological insights
.
Metabolites
9
,
57
.
Colinet
,
H.
,
Sinclair
,
B. J.
,
Vernon
,
P.
and
Renault
,
D.
(
2015
).
Insects in fluctuating thermal environments
.
Ann. Rev. Entomol.
60
,
123
-
140
.
de Paepe
,
R.
,
Lemaire
,
S. D.
and
Danon
,
A.
(
2014
).
Cardiolipin at the heart of stress response across kingdoms
.
Plant Signal. Behav
.
9
,
e29228
.
Dettloff
,
M.
,
Wittwer
,
D.
,
Weise
,
C.
and
Wiesner
,
A.
(
2001
).
Lipophorin of lower density is formed during immune responses in the lepidopteran insect Galleria mellonella
.
Cell Tissue Res.
306
,
449
-
458
.
Elkin
,
C. M.
and
Reid
,
M. L.
(
2005
).
Low energy reserves and energy allocation decisions affect reproduction by Mountain Pine Beetles, Dendroctonus ponderosae
.
Funct. Ecol.
19
,
102
-
109
.
Ernst
,
R.
,
Ejsing
,
C. S.
and
Antonny
,
B.
(
2016
).
Homeoviscous adaptation and the regulation of membrane lipids
.
J. Mol. Biol.
428
,
4776
-
4791
.
Fu
,
L.
,
Niu
,
B.
,
Zhu
,
Z.
,
Wu
,
S.
and
Li
,
W.
(
2012
).
CD-HIT: accelerated for clustering the next-generation sequencing data
.
Bioinformatics (Oxford, England)
28
,
3150
-
3152
.
Gray
,
E. M.
(
2013
).
Thermal acclimation in a complex life cycle: The effects of larvaland adult thermal conditions on metabolic rate and heat resistance in Culex pipiens (Diptera: Culicidae)
.
J. Insect Physiol.
59
,
1001
-
1007
.
Griesser
,
M.
,
Boeglin
,
W. E.
,
Suzuki
,
T.
and
Schneider
,
C.
(
2009
).
Convergence of the 5-LOX and COX-2 pathways: heme-catalyzed cleavage of the 5S-HETE-derived di-endoperoxide into aldehyde fragments
.
J. Lipid Res.
50
,
2455
.
Haas
,
B. J.
,
Papanicolaou
,
A.
,
Yassour
,
M.
,
Grabherr
,
M.
,
Blood
,
P. D.
,
Bowden
,
J.
,
Couger
,
M. B.
,
Eccles
,
D.
,
Li
,
B.
,
Lieber
,
M.
, et al. 
(
2013
).
De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis
.
Nat. Protoc.
8
,
1494
-
1512
.
Hack
,
M. A.
(
1997
).
The effects of mass and age on standard metabolic rate in house crickets
.
Physiol. Entomol.
22
,
325
-
331
.
Hahn
,
D. A.
and
Denlinger
,
D. L.
(
2007
).
Meeting the energetic demands of insect diapause: nutrient storage and utilization
.
J. Insect Physiol.
53
,
760
-
773
.
Huson
,
D. H.
,
Beier
,
S.
,
Flade
,
I.
,
Górska
,
A.
,
El-Hadidi
,
M.
,
Mitra
,
S.
,
Ruscheweyh
,
H.-J.
and
Tappu
,
R.
(
2016
).
MEGAN community edition - interactive exploration and analysis of large-scale microbiome sequencing data
.
PLoS Comput. Biol.
12
,
e1004957
.
Janzen
,
D.
(
1967
).
Why mountain passes are higher in tropics
.
Am. Nat.
101
,
233
-
249
.
Kanehisa
,
M.
and
Goto
,
S.
(
2000
).
KEGG: Kyoto encyclopedia of genes and genomes
.
Nucleic Acids Res.
28
,
27
-
30
.
Kaufmann
,
C.
and
Brown
,
M. R.
(
2008
).
Regulation of carbohydrate metabolism and flight performance by a hypertrehalosaemic hormone in the mosquito Anopheles gambiae
.
J. Insect Physiol.
54
,
367
-
377
.
Kim
,
D.
,
Langmead
,
B.
and
Salzberg
,
S. L.
(
2015
).
HISAT: a fast spliced aligner with low memory requirements
.
Nat. Methods
12
,
357
-
360
.
Lê Cao
,
K.-A.
,
Boitard
,
S.
and
Besse
,
P.
(
2011
).
Sparse PLS discriminant analysis: biologically relevant feature selection and graphical displays for multiclass problems
.
BMC Bioinformatics
12
,
253
.
Liao
,
Y.
,
Smyth
,
G. K.
and
Shi
,
W.
(
2014
).
featureCounts: an efficient general purpose program for assigning sequence reads to genomic features
.
Bioinformatics
30
,
923
-
930
.
Love
,
M. I.
,
Huber
,
W.
and
Anders
,
S.
(
2014
).
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2
.
Genome Biol.
15
,
550
.
Luo
,
W.
and
Brouwer
,
C.
(
2013
).
Pathview: an R/Bioconductor package for pathway-based data integration and visualization
.
Bioinformatics (Oxford, England)
29
,
1830
-
1831
.
MacMillan
,
H. A.
,
Knee
,
J. M.
,
Dennis
,
A. B.
,
Udaka
,
H.
,
Marshall
,
K. E.
,
Merritt
,
T. J.
and
Sinclair
,
B. J.
(
2016
).
Cold acclimation wholly reorganizes the Drosophila melanogaster transcriptome and metabolome
.
Sci. Rep.
6
,
28999
.
Marshall
,
K. E.
and
Sinclair
,
B. J.
(
2012
).
The impacts of repeated cold exposure on insects
.
J. Exp. Biol.
215
,
1607
-
1613
.
Marshall
,
K. E.
and
Sinclair
,
B. J.
(
2015
).
The relative importance of number, duration and intensity of cold stress events in determining survival and energetics of an overwintering insect
.
Funct. Ecol.
29
,
357
-
366
.
Meyer-Fernandes
,
J. R.
,
Gondim
,
K. C.
and
Wells
,
M. A.
(
2000
).
Developmental changes in the response of larval Manduca sexta fat body glycogen phosphorylase to starvation, stress and octopamine
.
Insect Biochem. Mol. Biol.
30
,
415
-
422
.
Morash
,
A. J.
,
Neufeld
,
C.
,
MacCormack
,
T. J.
and
Currie
,
S.
(
2018
).
The importance of incorporating natural thermal variation when evaluating physiological performance in wild species
.
J. Exp. Biol.
,
221
,
jeb164673
.
Moriya
,
Y.
,
Itoh
,
M.
,
Okuda
,
S.
,
Yoshizawa
,
A. C.
and
Kanehisa
,
M.
(
2007
).
KAAS: an automatic genome annotation and pathway reconstruction server
.
Nucleic Acids Res.
35
,
W182
-
W185
.
Mullen
,
L.
and
Goldsworthy
,
G.
(
2003
).
Changes in lipophorins are related to the activation of phenoloxidase in the haemolymph of Locusta migratoria in response to injection of immunogens
.
Insect Biochem. Mol. Biol.
33
,
661
-
670
.
Nichols
,
E.
,
Spector
,
S.
,
Louzada
,
J.
,
Larsen
,
T.
,
Amezquita
,
S.
and
Favila
,
M. E.
(
2008
).
Ecological functions and ecosystem services provided by Scarabaeinae dung beetles
.
Biol. Conserv.
141
,
1461
-
1474
.
Noriega
,
F. G.
(
2014
).
Juvenile hormone biosynthesis in insects: what is new, what do we know, and what questions remain?
Research Article
2014
,
967361
.
O'Donnell
,
V. B.
,
Dennis
,
E. A.
,
Wakelam
,
M. J. O.
and
Subramaniam
,
S.
(
2019
).
LIPID MAPS: Serving the next generation of lipid researchers with tools, resources, data, and training
.
Sci. Signal.
12
,
eaaw2964
.
Park
,
J.-H.
,
Attardo
,
G. M.
,
Hansen
,
I. A.
and
Raikhel
,
A. S.
(
2006
).
GATA factor translation is the final downstream step in the amino acid/target-of-rapamycin-mediated vitellogenin gene expression in the anautogenous mosquito Aedes aegypti
.
J. Biol. Chem.
281
,
11167
-
11176
.
Pásztor
,
L.
,
Kisdi
,
É.
and
Meszéna
,
G.
(
2000
).
Jensen's inequality and optimal life history strategies in stochastic environments
.
Trends Ecol. Evol.
15
,
117
-
118
.
Pertea
,
M.
,
Pertea
,
G. M.
,
Antonescu
,
C. M.
,
Chang
,
T.-C.
,
Mendell
,
J. T.
and
Salzberg
,
S. L.
(
2015
).
StringTie enables improved reconstruction of a transcriptome from RNA-seq reads
.
Nat. Biotechnol.
33
,
290
-
295
.
Piiroinen
,
S.
,
Lindström
,
L.
and
Lyytinen
,
A.
(
2010
).
Resting metabolic rate can vary with age independently from body mass changes in the Colorado potato beetle, Leptinotarsa decemlineata
.
J. Insect Physiol.
56
,
277
-
282
.
Podrabsky
,
J. E.
and
Somero
,
G. N.
(
2004
).
Changes in gene expression associated with acclimation to constant temperatures and fluctuating daily temperatures in an annual killifish Austrofundulus limnaeus
.
J. Exp. Biol.
207
,
2237
-
2254
.
Reynolds
,
J. A.
,
Bautista-Jimenez
,
R.
and
Denlinger
,
D. L.
(
2016
).
Changes in histone acetylation as potential mediators of pupal diapause in the flesh fly, Sarcophaga bullata
.
Insect Biochem. Mol. Biol.
76
,
29
-
37
.
Ruel
,
J. J.
and
Ayres
,
M. P.
(
1999
).
Jensen's inequality predicts effects of environmental variation
.
Trends Ecol. Evol.
14
,
361
-
366
.
Runcie
,
D. E.
,
Garfield
,
D. A.
,
Babbitt
,
C. C.
,
Wygoda
,
J. A.
,
Mukherjee
,
S.
and
Wray
,
G. A.
(
2012
).
Genetics of gene expression responses to temperature stress in a sea urchin gene network
.
Mol. Ecol.
21
,
4547
-
4562
.
Scheltema
,
R. A.
,
Jankevics
,
A.
,
Jansen
,
R. C.
,
Swertz
,
M. A.
and
Breitling
,
R.
(
2011
).
PeakML/mzMatch: a file format, Java library, R library, and tool-chain for mass spectrometry data analysis
.
Anal. Chem.
83
,
2786
-
2793
.
Schoville
,
S. D.
,
Barreto
,
F. S.
,
Moy
,
G. W.
,
Wolff
,
A.
and
Burton
,
R. S.
(
2012
).
Investigating the molecular basis of local adaptation to thermal stress: population differences in gene expression across the transcriptome of the copepod Tigriopus californicus
.
BMC Evol. Biol.
12
,
170
.
Sheldon
,
K. S.
and
Tewksbury
,
J. J.
(
2014
).
The impact of seasonality in temperature on thermal tolerance and elevational range size
.
Ecology
95
,
2134
-
2143
.
Sheldon
,
K. S.
,
Huey
,
R. B.
,
Kaspari
,
M.
and
Sanders
,
N. J.
(
2018
).
Fifty years of mountain passes: a perspective on Dan Janzen's classic article
.
Am. Nat.
191
,
553
-
565
.
Shinoda
,
T.
and
Itoyama
,
K.
(
2003
).
Juvenile hormone acid methyltransferase: a key regulatory enzyme for insect metamorphosis
.
Proc. Natl Acad. Sci. USA
100
,
11986
-
11991
.
Somero
,
G. N.
(
2004
).
Adaptation of enzymes to temperature: searching for basic “strategies”
.
Comp. Biochem. Physiol. B Biochem. Mol. Biol.
139
,
321
-
333
.
Somero
,
G. N.
,
Lockwood
,
B. L.
and
Tomanek
,
L.
(
2017
).
Biochemical Adaptation: Response to Environmental Challenges from Life's Origins to the Anthropocene
, p.
572
.
Sunderland, MA
:
Sinauer Associates
.
Sørensen
,
J. G.
,
Schou
,
M. F.
,
Kristensen
,
T. N.
and
Loeschcke
,
V.
(
2016
).
Thermal fluctuations affect the transcriptome through mechanisms independent of average temperature
.
Sci. Rep.
6
,
30975
.
Stillman
,
J. H.
and
Tagmount
,
A.
(
2009
).
Seasonal and latitudinal acclimatization of cardiac transcriptome responses to thermal stress in porcelain crabs, Petrolisthes cinctipes: acclimatization of thermal stress responses
.
Mol. Ecol.
,
18
,
4206
-
4226
.
Storey
,
K. B.
(
1997
).
Organic solutes in freezing tolerance
.
Comp. Biochem. Physiol. A Physiol.
117
,
319
-
326
.
Supek
,
F.
,
Bošnjak
,
M.
,
Škunca
,
N.
and
Šmuc
,
T.
(
2011
).
REVIGO summarizes and visualizes long lists of gene ontology terms
.
PLoS ONE
6
,
e21800
.
Tautenhahn
,
R.
,
Böttcher
,
C.
and
Neumann
,
S.
(
2008
).
Highly sensitive feature detection for high resolution LC/MS
.
BMC Bioinformatics
9
,
504
.
Tennessen
,
J. M.
,
Barry
,
W. E.
,
Cox
,
J.
and
Thummel
,
C. S.
(
2014
).
Methods for studying metabolism in Drosophila
.
Methods (San Diego, Calif.)
68
,
105
-
115
.
van Noordwijk
,
A. J.
and
de Jong
,
G.
(
1986
).
Acquisition and allocation of resources: their influence on variation in life history tactics
.
Am. Nat.
128
,
137
-
142
.
Vasseur
,
D. A.
,
DeLong
,
J. P.
,
Gilbert
,
B.
,
Greig
,
H. S.
,
Harley
,
C. D. G.
,
McCann
,
K. S.
,
Savage
,
V.
,
Tunney
,
T. D.
and
O'Connor
,
M. I.
(
2014
).
Increased temperature variation poses a greater risk to species than climate warming
.
Proc. R. Soc. B Biol. Sci.
281
,
20132612-20132612
.
Wan
,
Y.
,
King
,
R.
,
Mitchell
,
R. A. C.
,
Hassani-Pak
,
K.
and
Hawkesford
,
M. J.
(
2017
).
Spatiotemporal expression patterns of wheat amino acid transporters reveal their putative roles in nitrogen transport and responses to abiotic stress
.
Sci. Rep.
7
,
1
-
13
.
Watanabe
,
M.
,
Kikawada
,
T.
,
Minagawa
,
N.
,
Yukuhiro
,
F.
and
Okuda
,
T.
(
2002
).
Mechanism allowing an insect to survive complete dehydration and extreme temperatures
.
J. Exp. Biol.
205
,
2799
-
2802
.
Wijenayake
,
S.
,
Hawkins
,
L. J.
and
Storey
,
K. B.
(
2018
).
Dynamic regulation of six histone H3 lysine (K) methyltransferases in response to prolonged anoxia exposure in a freshwater turtle
.
Gene
649
,
50
-
57
.
Williams
,
C. M.
,
Marshall
,
K. E.
,
MacMillan
,
H. A.
,
Dzurisin
,
J. D. K.
,
Hellmann
,
J. J.
and
Sinclair
,
B. J.
(
2012
).
Thermal variability increases the impact of autumnal warming and drives metabolic depression in an overwintering butterfly
.
PLoS ONE
7
,
e34470-e34470
.
Wishart
,
D. S.
,
Tzur
,
D.
,
Knox
,
C.
,
Eisner
,
R.
,
Guo
,
A. C.
,
Young
,
N.
,
Cheng
,
D.
,
Jewell
,
K.
,
Arndt
,
D.
,
Sawhney
,
S.
, et al. 
(
2007
).
HMDB: the human metabolome database
.
Nucleic Acids Res.
35
,
D521
-
D526
.
Zhang
,
J.
,
Marshall
,
K. E.
,
Westwood
,
J. T.
,
Clark
,
M. S.
and
Sinclair
,
B. J.
(
2011
).
Divergent transcriptomic responses to repeated and single cold exposures in Drosophila melanogaster
.
J. Exp. Biol.
214
,
4021
-
4029
.

Competing interests

The authors declare no competing or financial interests.

Supplementary information