Physiological responses to fluctuating temperatures are characterized by distinct transcriptional profiles in a solitary bee

In ectothermic animals such as insects, environmental temperatures play a crucial role in the animal's physiology. Exposure to stressful, low temperatures during development can have significant, detrimental effects on adult insects (Bale and Hayward, 2010; Bennett et al., 2013; O'Neill et al., 2011; Whitfield and Richards, 1992; Yocum et al., 1994, 2006). The degree to which insects respond to these temperature stressors is highly dependent upon the timing, duration and amplitude of exposure (Colinet et al., 2015). With global climate change threatening to increase temperature variance (Easterling et al., 2000; Hansen et al., 2012), insects will likely be faced with increased exposure to low, seasonally unpredictable temperatures (Vasseur et al., 2014). In chill-susceptible insects, prolonged exposure to low temperatures can lead to the accumulation of deleterious physiological effects called chill injury (Lee, 2010).

Chill injury causes significant increases in mortality (Colinet et al., 2011; Kostál, 2006; Lee, 2010; Renault et al., 2004). The downstream physiological consequences of chill injury are complex but appear to be caused in part by membrane phase transitions (Lee, 2010). Damage to cell membranes results in metabolic imbalances and perturbations in ion homeostasis (Kostál et al., 2004, 2006) that have been observed in chill-injured insects. Furthermore, temperature stress is associated with oxidative stress (Lalouette et al., 2011; Rojas and Leopold, 1996) and impairment of neuromuscular system function (Bennett et al., 2015; Yocum et al., 1994), both of which are implicated as downstream consequences of cell membrane damage (Lee, 2010).

The decreased survival and cellular damage associated with chill injury can be mitigated by exposure to transient, daily increases in temperature during chilling [termed fluctuating thermal regime (FTR)]. Exposure to these brief, daily temperature increases during chilling results in increases in survival in chill-injured insects (Chen and Denlinger, 1992; Colinet et al., 2006, 2007; Jing et al., 2005; Kostál et al., 2007, 2016; Leopold et al., 1998; Marshall and Sinclair, 2012; Renault et al., 2004, 2011; Yocum et al., 2012). Exposure to FTR during chilling allows for either the protection against or repair of damage caused by chill injury (Colinet et al., 2015; Lee, 2010). While it is possible that brief, daily exposure to warmer temperatures during chilling allows for a physiological preparation (i.e. protection) for chilling, the hypothesis that exposure to warmer temperatures allows reparative mechanisms to act has more empirical support (Colinet et al., 2006; Kostál et al., 2007; Nedvěd et al., 1998; Torson et al., 2015). FTR repairs chill injury damage in many taxa (Boardman et al., 2013; Jing et al., 2005; Kostál et al., 2007; Leopold et al., 1998; Nedvěd et al., 1998; Renault et al., 2004), but it is unclear whether the mechanisms governing the beneficial effects associated with FTR exposure are conserved.

Current management practices for the solitary alfalfa leafcutting bee, Megachile rotundata (Fabricius 1787), a valuable agricultural pollinator, can include stressful low temperature exposure at two life stages: prepupal post-diapause quiescence and pupal development. Megachile rotundata females emerge and provision offspring in the early summer. Larvae complete growth in about a week and enter a prepupal diapause for overwintering. Under commercial management, the overwintering period lasts 7–10 months (Pitts-Singer and Cane, 2011), but can be artificially extended to over a year (Rinehart et al., 2013). However, extended overwintering causes chill injury (Rinehart et al., 2013; Torson et al., 2015). In the spring, managers increased storage temperatures to initiate pupation and adult eclosion. If weather delays peak floral bloom, managers slow pupal development using low temperatures. Low temperature exposure during pupation also leads to chill injury (Rinehart et al., 2011). Diapausing prepupae and developing pupae are very different life stages both physiologically and molecularly (MacRae, 2010). Interestingly, exposure to FTR during low temperature exposure increases survival in both contexts (Rinehart et al., 2011, 2013), and reduces the sub-lethal effects of chill injury (Bennett et al., 2015).

Gene expression comparisons could reveal whether the response to FTR is conserved across different developmental stages of M. rotundata. An RNA-seq study conducted in post-diapause quiescent prepupal M. rotundata provided support for the physiological evidence in other chill-injured insects and offered insight into mechanisms driving the reparative effect of FTR versus exposure to a constant 6°C (Torson et al., 2015). Individuals exposed to FTR had increased expression of transcripts functioning in the maintenance of ion homeostasis, response to oxidative stress and an increase in metabolic function; all congruent with the hypothesis that FTR allows for repair of damage during the warm phase. These transcriptional profiles could represent a conserved mechanism. Given that FTR has similar protective effects on survival in both overwintering prepupae and pupae exposed to cold during development, we hypothesized that a conserved suite of genes would underlie the response. The goals of our study were to (1) identify differential gene expression in individuals exposed to FTR protocol during pupal development and (2) determine whether expression profiles observed in pupal FTR support a conserved reparative mechanism of FTR across life stages in M. rotundata. Developing pupae were exposed to either a constant, low-temperature stress (static thermal regime, STR) or the same low temperature with a daily, 1-h warm-temperature pulse for 7 days (FTR). RNA was harvested from bees exposed to either treatment after the warm pulse on day seven. We hypothesized that, as a result of a conserved mechanism driving the reparative effect of FTR across taxa, the exposure to fluctuating temperatures during low-temperature stress during this life stage would elicit transcriptional responses similar to prepupal responses observed previously (Torson et al., 2015).

Megachile rotundata were purchased from JWM Leafcutter (Nampa, ID, USA), and were of Canadian origin. Diapausing prepupae arrived as loose cells and were housed in Percival model I-30BLL reach-in incubators at 6±0.5°C under darkness until initiation of the experiment. To end overwinting and initiate development, prepupae were transferred to a Percival model I-30BLL incubator at 29°C. To assess the progression of pupal development prior to treatment, bees not destined for experimental treatment were dissected out of their brood cells and placed in 24-well plates to function as guide plates. The guide plates were stored in plastic containers with breathable lids and a saturated NaCl solution was used to maintain humidity at 75% (as described in Winston and Bates, 1960). Experimental bees were maintained at 29°C until ∼50% of the bees in the guide plates reached a developmental stage characterized by melanization of the eyes, called the ‘red-eye’ stage. At the red-eye stage, bees were grouped into one of two treatments: FTR or STR (Fig. 1). Bees reared under the FTR protocol were exposed to 6°C with a daily warm pulse of 20°C consisting of a 1 h ramp up to 20°C (0.23°C min⁻¹), a 1 h incubation at 20°C, and a 1 h ramp down back to 6°C for each of the 7 days of treatment. Those reared under the STR treatment were maintained at 6°C for the duration of the 7-day treatment.

Bees used in the experiments were collected over two field seasons. The RNA-seq and validation was conducted on bees derived from the 2013 summer field season, which arrived in the laboratory in late March 2014. The qPCR time series was conducted on bees derived from the 2014 summer field season, which arrived in the laboratory in late March 2015. Before arriving at the laboratory, prepupal bees were overwintered at a constant temperature (4–6°C) by the supplier (JWM Leafcutter).

Total RNA was collected from whole M. rotundata pupae exposed to either the FTR or STR treatment (three biological replicates per treatment; one biological replicate is one individual). To ensure the assessment of response to low temperature exposure and not a warm-temperature exposure, individuals from both treatments were harvested for RNA during the cold phase (6°C; Fig. 1) and maintained on ice until the animals were harvested. Animals were harvested during the cold phase (after FTR warm pulse) to ensure a measurement of low temperature response in both treatments. Samples were harvested shortly after (∼1 h) the warm pulse in the event that the effect was transient. Pupae were transferred directly from ice and homogenized in 500 μl of TRIzol using a polypropylene pestle in a nuclease-free 1.5 ml microcentrifuge tube. The Invitrogen TRIzol protocol (Carlsbad, CA, USA) was followed for RNA extraction. Final product was stored as an ethanol precipitate at −80°C until needed.

Prior to sequencing, RNA pellets were dissolved in DEPC-treated H₂O and shipped at a concentration of no less than 20 ng µl⁻¹ (∼10 µg total RNA) on dry ice overnight to the University of Georgia. Upon arrival to the University of Georgia Genomics Facility (Athens, GA, USA), RNA samples from both FTR and STR treatments were checked for quality using an Aligent Bioanalyzer (Santa Clara, CA, USA). A total of 4 μg of total RNA was used for input in the Illumina Truseq Stranded mRNA kit. Samples were pooled and run on a NextSeq500 using a Version 1 Mid-Output 300 cycle kit with PE150 settings. Demultiplexing and adapter trimming were done by default when data were processed in BaseSpace (Illumina, San Diego, CA, USA).

Raw sequence data (BioProject accession: PRJNA352846) were quality checked using FastQC (version 0.10.1 http://www.bioinformatics.babraham.ac.uk/projects/fastqc/; Babraham Bioinformatics). TopHat (version 2.0.9; Trapnell et al., 2012) was used to align the raw reads against the M. rotundata reference genome (BioProject accession: PRJNA66515). Assembly of mapped reads was conducted using Cufflinks (version 2.0.2; Trapnell et al., 2012) and differential expression analysis was carried out using Cuffdiff (version 2.0.2; Trapnell et al., 2012) via the iPlant Collaborative Discovery Environment (Stanzione, 2011). Transcripts were considered significantly differentially expressed when the q-value (Benjamini–Hochberg-corrected P-value) was below the false discovery rate (FDR) of 0.05. The R package cummeRbund (Trapnell et al., 2012) and the Java-based functional annotation program Blast2GO (Conesa et al., 2005) were used for downstream analysis and the generation of differential expression figures.

To assess whether significant overlap between gene expression profiles exist during low-temperature stress between post-diapause quiescent prepupae (Torson et al., 2015) and developing pupae, Venn diagrams were constructed using the R package VennDiagram (Chen and Boutros, 2011). Statistical significance of the overlap of differentially expressed transcripts between data sets was determined via hypergeometric distribution using the R function phyper. The total number of transcripts for comparison was 11,195, which is the total number of transcripts in the reference transcriptome for M. rotundata.

Residual RNA from each RNA-seq sample (2013 field season; described above) was aliquoted and stored at −80°C under ethanol and used for the qPCR validation of the computational methods described above. The 10 most upregulated transcripts in response to FTR were selected for validation.

RNA samples were diluted to a concentration of 0.333 μg μl⁻¹ and subsequently treated with Invitrogen DNAase I to remove genomic DNA contaminants. A total of 1 μg RNA for each of the three biological replicates per treatment was used as template for first-strand cDNA synthesis using the Invitrogen Super Script III first strand synthesis system for RT-PCR.

Real-time quantitative PCR reactions were conducted on a Roche LightCycler 480 II (Indianapolis, ID, USA) using Roche FastStart Universal SYBR Green I Master Mix with ROX. Primers for all targets (10 most upregulated transcripts from RNA-seq experiment; Table S1) and eight reference gene candidates were designed using the Integrated DNA Technologies IDR program (Coralville, IA, USA). Candidate reference genes were selected using this RNA-seq data set and libraries generated in Torson et al. (2015) (18 libraries in total) to identify stably expressed transcripts. Expression of differentially expressed transcripts in qPCR runs were normalized against two most stably expressed reference targets identified using the algorithm geNorm (Vandesompele et al., 2002).

To ensure quality control between and within plates, three controls were used on each plate: a ‘no template control’ consisting of all enzymatic components without cDNA template, a ‘negative RT control’ lacking reverse transcriptase, and a ‘positive RT control’ using a control cDNA template and control primers for each plate. The positive RT control also served as calibrator reactions in the event that comparisons spanned multiple plates.

Biogazelle qbase+ (Ghent, Belgium) and Prism GraphPad (La Jolla, CA, USA) were used for statistical analysis and graphical representations of the data, respectively. Statistical significance for each of the 10 upregulated and eight non-differentially expressed transcripts was determined by t-test of the calibrated normalized relative quantity (CNRQ) of each transcript from qbase+.

Megachile rotundata from the 2014 field season were reared at 29°C until eye pigmentation stage and then split into either the FTR or STR protocol described above. Total RNA for the qPCR time series was collected from whole M. rotundata pupae exposed to either the FTR or STR treatment (six biological replicates; one individual per replicate). Before exposure to treatment, pupae were harvested to serve as a baseline for gene expression (T₀; Fig. 1). After placement into treatment, individuals were harvested during the cold phase (6°C, shortly after the warm pulse in FTR) for each of the 7 days of treatment (T₁–T₇; Fig. 1). Animals were placed on ice during dissection from brood cells until they were homogenized for RNA extraction to prevent an additional warm-temperature exposure. The Invitrogen TRIzol protocol was followed for RNA extraction. Final product was stored as an ethanol precipitate at −80°C until needed. Sample preparation and qPCR protocols were similar to RNA-seq validation (see previous section). Two-way ANOVA with Bonferroni's correction for multiple comparisons determined statistical significance of differential expression, and relationships across time, within treatments were determined using linear regression.

Three biological replicates were harvested from each treatment, totaling six individual transcriptomic libraries (BioProject accession: PRJNA352846). Each library averaged a total of 47,279,330 reads (Table 1). Raw reads were mapped to the M. rotundata reference genome (BioProject accession: PRJNA66515) with an average alignment rate of 69.11%. FastQC (version 0.10.1) was used to determine baseline quality of raw reads, and percentage of reads mapped to the reference genome and mate pair concordance (Table 1) were used as a metric of assembly quality.

Using the Tuxedo protocol (Trapnell et al., 2012), expression analysis revealed a total of 434 differentially expressed transcripts between STR and FTR treatments after 7 days of exposure to treatment (307 downregulated and 127 upregulated in FTR). Real-time quantitative PCR confirmed the expression profiles of eight transcripts upregulated in FTR, relative to STR, and six non-differentially expressed transcripts that we used as reference gene candidates (Table S2). Differential expression was confirmed in all cases with at least one primer pair. We were able to confirm expression of eight of the top 10 FTR-upregulated transcripts (Tables S2 and S3) and six of eight non-differentially expressed transcripts (Table S3).

A Benjamini–Hochberg-corrected Fisher's exact test (one-tailed, FDR=0.05) of upregulated and downregulated transcripts in response to FTR revealed no overrepresented or underrepresented functional classes relative to the gene ontology (GO) distribution in the reference transcriptome. Owing to the lack of enrichment for GO terms, relative to the reference transcriptome, we used direct counts of GO terms of the differentially expressed transcripts in this study for downstream analyses and discussions. Direct counts of the transcripts differentially expressed between FTR and STR treatments on T₇ of treatment revealed classes of transcripts functioning in biological processes including regulation of transcription, oxidation–reduction, signal transduction and various metabolic processes (Fig. S1). Also, a direct count of biological process GO terms (Fig. 2) revealed the upregulation in FTR of transcripts functioning in signal transduction and trans-membrane transport. Furthermore, GO terms associated with the construction of cellular components, including plasma membrane and various membrane components, were upregulated in response to FTR (Fig. 2).

When the identities of differentially expressed transcripts observed in this study are compared with those presented during the two time points assessed during post-diapause quiescence (Torson et al., 2015), 30 (6.9%) differentially expressed transcripts in the pupal response (434 transcripts) were also differentially expressed during the early time point (215 transcripts), before mortality diverged in prepupal STR/FTR comparison. This overlap is more than we would expect by chance, based on hypergeometric distribution (P=1.27E-10). Seventeen of the 434 differentially expressed transcripts (3.9%) overlapped with the late time point (after mortality diverged, 64 transcripts) in prepupae, also more than expected by chance (P=1.68E–11). Six (1.4%) transcripts differentially expressed in pupae were differentially expressed in both prepupal time points (Fig. 3A).

While differentially expressed transcripts overlap more than expected by chance, when direction of regulation (whether the transcript is upregulated or downregulated in FTR-exposed individuals) is considered, all overlap among the three differential expression profiles disappears (Fig. 3B,C). In fact, the only differential expression patterns conserved between the two life stages are regulated in opposite directions (Fig. 3D).

Exposure to daily fluctuations in temperature (FTR) during chilling can mitigate the deleterious physiological effects associated with chill injury (Colinet et al., 2015). FTR appears to provide a protective benefit across all taxa studied so far, but we lack evidence for a conserved mechanism driving this protective effect. In this study, we exposed developing M. rotundata pupae to an STR, known to lead to chill injury, or to an FTR, which diminishes the negative effects of chilling (Bennett et al., 2015). An RNA-seq analysis after 7 days of treatment revealed 434 transcripts that were differentially expressed between the two low-temperature treatments. These differentially expressed transcripts provide evidence for cellular damage, oxidative stress and repair of damage (Table 2; Table S4). Additionally, we observed a lack of conservation in differential gene expression in response to FTR exposure across life stages in M. rotundata, suggesting a robust physiological response driven by variable transcriptional mechanisms.

Plants, microbes, and other animals change the physical composition of cell membranes in response to changes in ambient temperature, namely the proportion of unsaturated fatty acids in the phospholipid bilayer (Cossins, 1983). In insects, direct evidence of shifts in membrane phase transitions have yet to be observed (Hayward et al., 2014; Lee, 2010), but strong correlative evidence has been provided (Boardman et al., 2011; Kostál et al., 2006, 2007; Macmillan et al., 2012). Chill injury in insects has been associated with the collapse of ion gradients presumably caused by damage at the cellular level and possibly resulting in decreased activity of membrane-bound enzymes (Hazel, 1995; MacMillan and Sinclair, 2011). Extracellular increases of potassium and diminished levels of sodium and magnesium have been observed in chill-injured tropical cockroaches (Kostál et al., 2004, 2006, 2007). These changes in hemolymph (extracellular) ion concentrations could explain the loss of neuromuscular coordination observed in chill-injured insects (Andersen et al., 2016; MacMillan and Sinclair, 2011) and may be a direct effect of damage to cell membranes. In this study, transcripts such as sodium potassium calcium exchanger 6, a mitochondrial potassium-dependent sodium–calcium antiporter (Table 2), are upregulated in individuals exposed to the STR protocol. Mitochondrial function decreases in Drosophila melanogaster exposed to constant low temperature (Colinet et al., 2017). Proper mitochondrial function would be heavily dependent on membrane composition, so the differential expression of this sodium–potassium–calcium exchanger suggests that STR-exposed individuals may be attempting to restore proper ion gradients in response to cellular damage. Additionally, the upregulation of transcripts with GO terms associated with the construction of various cellular components, including cell membranes and integral membrane components (Fig. 2, Table 2), in response to FTR exposure indicates that one of the FTR's mechanisms of action is to repair the membrane damage that has accumulated over the duration of low-temperature stress.

The upregulation of aquaporins, water-selective transmembrane channels, has been observed in freeze-tolerant insects (Philip et al., 2011), and the blockage of those channels has been shown to decrease freeze tolerance owing to an inability of the cells to move water out of the cell; subsequently causing damage as a result of osmotic pressure (Philip et al., 2008). We observed the upregulation of aquaporins during constant low-temperature stress (STR), suggesting that these individuals may be attempting to return to proper osmotic conditions. Furthermore, STR-exposed M. rotundata showed increased expression of two isoforms coding for the facilitated trehalose transporter TRET1. This transmembrane protein mediates the movement of trehalose synthesized in the fat body bidirectionally across the cell membrane (Kanamori et al., 2010). The presence of two isoforms for Tret-1 is consistent with orthologs in D. melanogaster (Zhao et al., 2000). The enrichment of Tret-1 in STR-reared individuals not only indicates changes in membrane composition, but suggests that these individuals may be transporting trehalose into the hemolymph because of its known role as a cryoprotectant (Bale and Hayward, 2010). Interestingly, transcripts coding for a trehalase precursor, an enzyme needed to convert trehalose into glucose (Shukla et al., 2015) and a facilitated glucose transporter (Table 2) were upregulated in FTR-reared individuals, indicating a shift in hemolymph composition between the two treatments.

In addition to changes in structural components within membranes, such as aquaporins and trehalose transporters, an increase in unsaturated fatty acids has been observed in response to decreased temperatures across taxa (Cossins, 1983). In our dataset, three different transcripts coding for acyl-delta desaturases, enzymes functioning in the transition between saturated and unsaturated fatty acids, were upregulated in STR-reared individuals (Table 2). The upregulation of transcripts coding for desaturase enzymes has been associated with cold hardiness in the onion maggot, Delia antiqua (Kayukawa et al., 2007), indicating this response is conserved. Increased expression of these transcripts suggests that either bees exposed to STR are trying to compensate for a decrease in temperature or FTR-reared individuals may have less of a strain on their cell membranes because of the temporary increases in temperature.

Oxidative stress has been implicated in the occurrence of chill injury in many insect species (Joanisse and Storey, 1996, 1998; Lalouette et al., 2011; Rojas and Leopold, 1996; Torson et al., 2015) and has been hypothesized to be a downstream consequence of membrane phase transitions (Lee, 2010). Quiescent post-wandering D. melanogaster larvae exposed to extended low-temperature regimes also show enrichment in transcripts that suggest redox imbalance and possibly oxidative stress (Koštál et al., 2016).

Gene expression indicates that STR-exposed individuals are repairing damage, possibly in response to oxidative stress. The upregulation of a gene coding for a myb domain-containing protein in STR-reared individuals (Table S4) implies that M. rotundata is indeed experiencing a stress response. In plant species, this class of transcription factors is known to respond to abiotic stressors, including decreased temperature. The functional unit of the translated proteins is highly conserved, implying a similar function in insects (Ambawat et al., 2013). Furthermore, increased expression of oxidase-peroxidase (Table S4) is indicative of STR-reared individuals experiencing increased levels of peroxide, a consequence of malfunction in mitochondrial respiration (Prasad et al., 1994). If antioxidant systems cannot combat the accumulation of free radicals, such as peroxide, damage to cellular structures would be likely.

The upregulation of transcripts in STR-reared bees such as myofilin isoform b, sestrin-like and DNA damage-binding protein 1 (Table S4) suggests that STR-reared individuals experience DNA damage, potentially caused by increased levels of reactive oxygen species in the bees during chilling. In Drosophila, expression levels of sestrins increase in response to target of rapamycin (TOR) under increased levels of reactive oxygen species (Lee et al., 2010a). The increased expression of these sestrins has also been linked to genotoxic stress including DNA damage (Hay, 2008; Lee et al., 2010b). While sestrins, a class of highly conserved cytoplasmic proteins, have not been shown to repair DNA damage directly, they may be part of a larger signaling cascade. Direct assessments of DNA damage will be required to assess whether oxidative stress may be playing a significant role in the deleterious effects of chill injury or whether high levels of antioxidant expression in FTR-reared individuals is just an artifact of increased metabolic rate owing to transient increases in ambient temperature.

In addition to the benefits of FTR observed during pupal development in M. rotundata, significant increases in survival have been observed when prepupal bees in post-diapause quiescence are exposed to a similar treatment (Rinehart et al., 2013). For the duration of diapause and post-diapause quiescence in both natural and managed populations, M. rotundata are exposed to low temperatures and have mechanisms to cope with these temperatures. Conversely, other life stages, such as developing pupae, which we have assessed in this study, may not be adapted for low-temperature exposure and may have different mechanisms driving low-temperature stress responses, and subsequently the beneficial effects of FTR. Adult bees observed after STR treatment in both life stages exhibit diminished neuromuscular function, behavioral abnormalities and diminished flight capacity (Bennett et al., 2013, 2015). These phenotypes suggest chill injury accumulation in both contexts. Initially, we predicted that the temperature stress during the pupal stage would elicit transcriptional responses similar to the prepupal responses observed in Torson et al. (2015). This hypothesis was not supported. Markedly different differential expression profiles (Fig. 3) between the two life stages indicates that either (1) the mechanism responsible for a conserved FTR benefit vary throughout development or (2) the benefit of FTR exposure changes throughout development.

Even within the pupal life stage, temperature exposure can have variable effects on gene expression depending on context. Survival benefits of FTR, relative to STR, are robust to field season (Bennett et al., 2015; Rinehart et al., 2016), but gene expression changes at a finer scale. Expression profiles of the most differentially expressed transcripts identified in this RNA-seq study (Table S2) were not maintained in M. rotundata pupae from the following field season (2014; Fig. S2). The direction of regulation, relative to STR, was maintained, but expression was highly variable, even with a doubling of the sample size. High variability between individuals suggests that slight differences in developmental timing are introducing variation or season-specific environmental history is influencing gene expression, or both. However, it should be noted that the individuals used for qPCR were exposed to one additional month of overwintering prior to pupal development. This short extension of the overwintering period, compared with the RNA-seq individuals, could have also influenced gene expression during subsequent pupal development. Gene expression in response to FTR exposure now appears to be highly context dependent (i.e. environmental history and developmental stage) and provides further evidence that the conservation of the benefit of FTR benefit may not be at the level of gene expression, but at a higher level of organization. This result provides the basis for our hypothesis that chill injury repair in response to FTR is a robust physiological response that is driven by life-stage-specific mechanisms.

In biological systems, an animal's ability to maintain a certain characteristic response under perturbations in environmental conditions is known as a robust response (Stelling et al., 2004). With similar physiological responses having been observed across taxa, and, thus far, a lack of evidence of a conserved mechanism at the level of gene expression, we hypothesize that the beneficial effects of chill injury repair during FTR exposure represents a robust physiological response driven by variable gene expression. Redundancy is one basic mechanism of a robust system. Redundancy provides alternative molecular strategies to carry out an action (Hartman et al., 2001). For example, glycosylphosphatidylinositol (GPI)-linked glycoproteins in yeast, termed fungal adhesions, have functionally distinct roles separated temporally and spatially (Guo et al., 2000), but members of this protein family can replace one another when inappropriately expressed in the organism. Similar to our results in M. rotundata, work in yeast has shown little overlap in gene expression profiles in response to stressors (Berry et al., 2011), suggesting that variable downstream effectors (i.e. gene expression response) may enact similar physiological responses (i.e. redundancy or genetic buffering).

Between the larval and pupal life stages in M. rotundata, a complete restructuring of the animal takes place, resulting in entirely different physiologies. Redundancy in transcript function, or ‘buffering’, may explain why variable gene expression patterns are observed in what appear to be similar physiological responses to FTR exposure. However, changes in gene or protein expression and metabolite accumulation might be secondary effects rather than direct or adaptive responses to chilling and/or FTR exposure. As a result of this notion, an alternative hypothesis is that the variable gene expression profiles between life stages in M. rotundata is not evidence of different life-stage-specific mechanisms driving the same response at the physiological level, but actually different responses driven by stimuli at a higher level of organization. Additional functional genomic work to identify potentially redundant genes and physiological assays to verify similar physiological responses to FTR exposure will be necessary to provide support for the hypothesis that redundancy allows for conservation of the FTR response across life stages in M. rotundata.

Expression of a subset of the most significantly upregulated and downregulated transcripts in the RNA-seq study after 7 days of treatment (Table S2) was assessed across the duration of low-temperature stress (Fig. 1) using qPCR. Generally, direction of regulation between treatments at T₇ was maintained between this time series (different field season) and the RNA-seq time point, but the qPCR time series had greater variation in expression between biological replicates within treatment and time point, resulting in a lack of statistical significance even after sample sizes were increased from three to six. Expression of the transcript coding for a zinc finger SCAN domain-containing protein had a significant time effect (F_7,70=3.752, P=0.0012), showing an increase in expression in individuals exposed to STR (F_1,46=17.04, P=0.0002) over time and a significant difference between treatments at T₆ (P=0.0183). Expression profiles of Cabut for STR and FTR both decreased over time (STR, F_1,46=6.026, P=0.0179; FTR, F_1,46=8.779, P=0.0048), with time (F_7,70=3.609, P=0.0022) and treatment (F_7,70=7.180, P=0.0231) explaining a significant amount of variation in gene expression, but only T₂ was significantly different between treatments (P=0.0085; Fig. S2). The serine protease Stubble was upregulated in FTR-reared individuals at T₂ (P=0.0196).

The authors thank Marnie Larson of the USDA-ARS, Fargo, ND, USA, for her technical assistance and advice. We also thank Kendra J. Greenlee of North Dakota State University for her comments on the manuscript.