Different amplitudes of temperature fluctuation induce distinct transcriptomic and metabolomic responses in the dung beetle Phanaeus vindex

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 (CT_min) to an inflection point below the thermal optimum (T_opt), decelerating from the inflection point to T_opt and, finally, declining rapidly to the critical thermal maximum (CT_max) (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.

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.

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⁻¹. For each female, we measured the amount of CO₂ 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 H₂O 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 CO₂ and then recorded the amount of CO₂ produced by the beetle over approximately 5 min. We then took another baseline reading of CO₂ 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 CO₂, the CO₂ produced by the beetle over approximately 5 min, and a final baseline CO₂ reading. We continued this to obtain CO₂ 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 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.

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.

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⁻¹. 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).

As expected, metabolic rate increased as test temperature increased (F_1,11=8.75, P=0.01) but did not vary with female body size (F_1,11=3.82, P=0.08) (Fig. 1). Temperature fluctuation treatments did not affect the thermal sensitivity of metabolism (F_2,11=0.61, P=0.56) or overall metabolic rate (F_2,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 (F_2,12=1.17, P=0.344) or lipid (F_2,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).

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 N₅₀ 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).

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.

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.

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).

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).

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 1–3). 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’.

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.