A rapid return to normal: temporal gene expression patterns following cold exposure in the bumble bee Bombus impatiens

The ability to mount an appropriate physiological response to challenging abiotic conditions can determine the functional capacity and viability of organisms that encounter variable environments (West-Eberhard, 1989; Whitman and Agrawal, 2009). Under dynamic environmental conditions, individuals can benefit from flexible phenotypic responses that can match fluctuations in external stimuli (Lande, 2009; Schlichting and Pigliucci, 1993; Stomp et al., 2008). At the molecular level, individuals can exhibit such flexibility through differential gene expression, where the transcription of genes can be modified in response to environmental stressors, such as extreme temperature changes (Baniwal et al., 2004; DeSalvo et al., 2008; Schoville et al., 2012; Shearer et al., 2016; Sonna et al., 2002). Given ongoing temperature shifts associated with global climate change, there is substantial interest in characterizing the physiological and molecular responses of extreme thermal events (Buckley and Somero, 2009; DeSalvo et al., 2008; Harvey et al., 2023; Heikkila, 2017; Sonna et al., 2002).

For insects, thermal tolerance is a driving factor that shapes species distributions and life history strategies (Bale, 2002; Kellermann et al., 2012; Oyen et al., 2016; Stange and Ayres, 2010; Sunday et al., 2011), and gene expression modification is a key component of responses to temperature fluctuations (Des Marteaux et al., 2017; Kim et al., 2006; Lubawy et al., 2022; Malmendal et al., 2006; Sinclair et al., 2007; Teets and Denlinger, 2013; Xu et al., 2017; Zhao and Jones, 2012). Climate change has driven substantial research on tolerance of warming conditions; however, understanding the mechanisms underlying cold tolerance may be as important as studying upper thermal limits, as exposure to cold events (Teets et al., 2023) and the frequency and severity of extreme weather events (Ummenhofer and Meehl, 2017) will continue despite rising mean temperatures. Interestingly, while the upper range for heat tolerance (i.e. the maximum critical thermal limit, or CT_max) is relatively similar across insects (∼40°C; Neven, 2000) and generally stable within species (Sunday et al., 2012), species and populations can vary in cold tolerance (Bowler and Terblanche, 2008; Oyen et al., 2016; Pimsler et al., 2020; Sunday et al., 2014). For example, many insects exhibit large-scale gene expression alterations that occur during acclimation to gradual seasonal changes, such as the weeks or months leading to diapause (Costa et al., 2020; Denlinger et al., 2001; MacRae, 2010), although these drastic alterations can also occur in as little as 6 days (MacMillan et al., 2016). Sudden thermal stress and cold hardening tend to elicit more subtle and immediate transcriptional responses over minutes to hours (Teets et al., 2020, 2023). Such immediate responses can be crucial for species with life history strategies that expose individuals to more rapid shifts in temperature (Bale et al., 2002; Kingsolver et al., 2011; Sinclair, 2001), such as species active at different times of day or during seasonal transitions (Colinet et al., 2015). Several mechanisms that can influence cold tolerance are similar to those observed in response to heat exposure, such as the increase in metabolism (particularly glycolysis) and alterations in the expression of heat shock proteins (Teets and Denlinger, 2013). The production of cryoprotectants, shifts in cellular membrane composition and in ion transport, and cytoskeletal restructuring are some of the other mechanisms that may be controlled by rapid molecular changes to protect against rapid temperature reductions (MacMillan et al., 2014; Overgaard and MacMillan, 2017; Robertson et al., 2017; Teets and Denlinger, 2013). Such rapid responses of enhanced cold tolerance are protective but energetically expensive (Chen and Walker, 1994), and benefits can rapidly dissipate (Wang and Kang, 2003). Some work has evaluated shifts in gene expression for insects subjected to acute cold (Sinclair et al., 2007; Teets et al., 2012; Zhang et al., 2011), but for most species we still require a better characterization of transcriptomic profiles immediately after extreme thermal events to better understand both how organisms maintain or restore physiological functions and whether effects of thermal stress persist beyond the period of exposure.

Bumble bees (Hymenoptera: Apidae: Bombus) provide an interesting model organism with which to evaluate responses to cold stress. Bombus evolution and biogeography has been shaped by global cooling periods (Hines, 2008; Martinet et al., 2018; Williams et al., 2018), and species diversity is high in temperate and cold (e.g. montane, arctic) environments. Although bumble bees establish thermally controlled colonies in cavities during the growing season, individuals of different castes can face varying levels of cold exposure risk while outside the nest. For example, queens emerging in the spring may spend weeks unprotected above ground, risking exposure to early spring snowstorms before settling underground in their established colonies (Makinson et al., 2019). Workers can forage >1 km from the colony and occasionally remain outside the colony overnight, increasing their risk of experiencing extreme cold (Free, 1955; Müller and Schmid-Hempel, 1993; Walther-Hellwig and Frankl, 2000). To promote heat regulation and retention in the face of thermal challenges, bumble bees are facultatively endothermic insects (Heinrich, 2004), with morphological features such as large body size and dense pile that aid in temperature regulation (Heinrich, 2004; Maebe et al., 2021b; Peters et al., 2016). Bumble bees can thus maintain elevated body temperatures in conditions that are often too cold for many ectothermic insects to function and are capable of activities earlier in the day or season compared with other pollinators (Goulson, 2010). The metabolic powerhouse of heat production in bumble bees lies in their thoracic muscles (Heinrich, 2004), which not only power flight but also allow bumble bees to warm up prior to flight and other activities (Esch et al., 1991; Heinrich, 1993; Heinrich and Esch, 1994). However, when cooled sufficiently, bees and other insects reach their critical thermal minimum (CT_min), when loss of central nervous system function and impairment of peripheral neuromuscular junctions leads to a state of arrested movement known as chill coma (Andersen et al., 2015; MacMillan and Sinclair, 2011). Previous studies have demonstrated regional differences in bumble bee CT_min among both species (Maebe et al., 2021a; Oyen et al., 2016) and populations (Pimsler et al., 2020). For example, Bombus vosnesenskii workers reared from queens sampled across an elevational and latitudinal gradient had CT_min strongly correlated with temperatures at queen collection sites, whereas CT_max values were relatively invariable (Pimsler et al., 2020). Bees from different regions displayed variations in gene expression in response to cold stress, while responses to heat stress were strong but conserved among regions (Pimsler et al., 2020). Because bumble bees are intolerant of freezing (Keaveny et al., 2022), and must navigate variable conditions without approaching their relatively high freezing point to prevent cold injury and death, molecular mechanisms for tolerating periods of cold stress may be an important component of bumble bee environmental adaptation.

In this study, we further investigated differential gene expression responses associated with cold stress in the common eastern bumble bee, Bombus impatiens Cresson 1863, a wide-ranging generalist pollinator found across eastern and central North America (Williams et al., 2014). We aimed to test whether acute cold exposure has lasting effects on gene expression or whether the high metabolic rate (and, therefore, fast metabolic turnover) of facultatively endothermic bumble bees leads to a rapid and short-lived molecular response. We induced chill coma and assessed differences in partial mRNA transcripts using 3′ tag RNA sequencing (TagSeq; Lohman et al., 2016; Meyer et al., 2011) at pre-determined time points between cold-exposed and control bees to capture potential effects immediately after cold exposure and up to several hours after recovery from chill coma. If differentially expressed genes are detectable earlier in a sampling time series, we hypothesize that such genes indicate immediate and short-term B. impatiens chill coma recovery processes. Alternatively, if differentially expressed genes are important for long-term recovery and are potentially beneficial for protection from subsequent cold exposures, expression of such markers should persist and be detectable in later sampling periods.

Three commercial B. impatiens colonies (Koppert Biological Systems, Howell, MI, USA) were maintained at 28°C with 60% relative humidity (RH) on a 12 h:12 h light:dark cycle in an incubator (Model #I36VL, Percival Scientific, Perry, IA, USA) in Laramie, WY, USA. They were fed nectar ad libitum via a wick and reservoir provided by the supplier in addition to 1 tbsp (15 ml) ground fresh frozen pollen (Prairie River Honey Farm, Grand Island, NE, USA) every other day. Seventy-two hours prior to experiments, 20 workers from each of the three commercial colonies (i.e. six microcolonies) were randomly selected and placed in colony-specific control and treatment chambers (Fig. 1A). The bees were fed a pollen ball (mix of ground pollen and nectar) and nectar [1:1:2 sucrose, invert sugar and distilled water, with <0.1% sorbic acid, Honey B Healthy (Honey B Healthy, Cumberland, MD, USA), Amino-B Booster (Honey B Healthy) and citric acid added] ad libitum through a wick inserted into a nectar-filled vial attached to the chamber wall; a small piece of comb from the respective natal commercial colony was added to each microcolony chamber. While still in the chambers, the six microcolonies were returned to the incubator after changing the light:dark cycle to 16 h:8 h, where they remained for a 72 h settlement period.

After the 72 h settlement period, we placed the three treatment microcolonies (one from each commercial colony) into a small plastic container (to ensure that the microcolonies stayed dry) that was submerged in an ice bath inside a cooler (ambient air measured 0°C) for 75 min (Fig. 1B). Conditions for the remaining three microcolonies were held constant in the 28°C incubator to serve as the control. At the end of the 75 min cold exposure, three bees from each of the six microcolonies (three control and three treatment microcolonies) were randomly selected and flash-frozen (time 0 after cold exposure, n=18, referred to as 0 min), before returning the cold-treated microcolonies with all bees still ventrally curled from chill coma to the incubator. Once in the incubator, we started a timer in the view of a video camera (Go-Pro, San Mateo, CA, USA) mounted above the cold-treated microcolonies, which allowed measurement of recovery time by monitoring movements of individual bees. All bees sampled at 0 min were in chill coma whereas all other bees sampled were fully recovered from chill coma. Three bees from each of the six microcolonies were flash-frozen as above, at 10, 30, 120 and 720 min after the cold exposure ended (n=9 bees per sampling time, per treatment type; Fig. 1). All bees were stored at −80°C until preparation for RNA extraction. To split bees, a metal stage was positioned on top of a block of dry ice on a sterilized surface. One bee was removed from the cryovial, weighed, placed on the stage, and cut along the sagittal plane using a sterile razor and forceps. There were no significant differences in mass among treatments or time points (ANOVA, all P>0.148). Each half was placed into a separate cryovial and returned to −80°C, with one half for RNA isolation and the other reserved for archiving and possible follow-up experiments.

Because of the significance of flight muscle in thermoregulation (see Introduction), RNA was isolated from thoracic tissue (Invitrogen RNAqueous-Micro Total RNA Isolation Kit, Thermo Fisher Scientific) following the manufacturer's instructions, with the extraction order of samples randomly selected across colonies and treatments to minimize potential batch effects. Aliquots of RNA at 10 ng µl⁻¹ were sent on dry ice to The University of Texas at Austin Genomic Sequencing and Analysis Facility (Austin, TX, USA) for quality control and 3′ Tag-based RNA sequencing (TagSeq) library preparation and sequencing (Lohman et al., 2016; Meyer et al., 2011). Details are provided in Lohman et al. (2016) but, briefly, TagSeq is a reduced representation approach for RNA sequencing of the 3′ ends of transcripts. cDNA libraries are prepared by isolating mRNA and subjecting the transcripts to a heat-based fragmentation to reduce length bias. A unique oligo-dT primer targeting the 3′ end of fragmented transcripts is used for first strand synthesis, and the libraries are PCR amplified for sequencing (Lohman et al., 2016; Meyer et al., 2011). Sequencing (100 bp, single-end) was conducted across two lanes on a NovaSeq 6000 SR100 (Illumina) flowcell. Samples used in analyses were sequenced to a mean±s.d. of 12.8±2.3 million reads per sample; a small number of samples that did not produce successful libraries or meet our sequence data threshold (<8 million reads) were excluded (see Table S1 for final sample sizes per colony per treatment). The final number of samples used in the analyses per time point, per treatment was as follows: 0 min control bees n=4, 0 min treatment bees n=8; 10 min control bees n=8, 10 min treatment bees n=7; 30 min control bees n=8, 30 min treatment bees n=8; 120 min control bees n=8, 120 min treatment bees n=9; 720 min control bees n=7, 720 min control bees n=7).

Quality assessment of reads was conducted using FastQC v0.11.9 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc). The two lanes were concatenated into a single FASTQ file per sample. As recommended for TagSeq data, the tagseq_clipper.pl script from the ‘tag-based_RNAseq’ pipeline by Mikhail Matz (https://github.com/z0on/tag-based_RNAseq, pipeline version Dec 29, 2017) with cutadapt v1.18 (Martin, 2011) was used to clean and deduplicate reads with the following command for each sample: tagseq_clipper.pl {SAMPLE_ID}.fq | cutadapt - -a AAAAAAAA -a AGATCGG -q 15 -m 25 -o {SAMPLE_ID}.trim. We used STAR v2.7.6a (Dobin et al., 2013) to align the cleaned reads to the B. impatiens genome (BIMP_2.2, GCF_000188095.3) (Sadd et al., 2015) and generate transcript counts per locus (sense strand). The R v4.0.3 (http://www.R-project.org/) package DESeq2 v1.30.1 (Love et al., 2014) was used for differential expression analysis from the transcript counts matrix. The DESeq() function from DESeq2 estimates size factors while accounting for differences in sequencing depth, calculates dispersion estimates for how much the counts of each sample will vary around expected values, and fits dispersion estimates to a generalized liner model, which provides the estimates of the log fold-change between contrasts. Because variance was substantially different among time points, we elected to analyze each time point separately for the DESeq2 analyses. To pinpoint the time at which gene expression differences between treatment and control bees become evident and/or dissipate, we analyzed models to test for a treatment effect at each time point, with source colony specified as a covariate. DESeq2 uses Wald testing to determine P-values in the gene expression analysis, and the Benjamini–Hochberg method (Benjamini and Hochberg, 1995) to account for multiple testing and provides an adjusted P-value. Our significance threshold for significant differentially expressed genes (DEGs) was an adjusted P-value ≤0.05. Plots and figures were created using DESeq2, the R package VennDiagram v1.7.3 (Chen and Boutros, 2011), and ggplot2 v3.3.6 (Wickham, 2011). The R package GeneOverlap v1.34.0 (https://github.com/shenlab-sinai/GeneOverlap) was used to determine whether DEG overlaps between time series were greater than expected using Fisher's exact tests.

As a secondary approach to identify genes that respond to cold treatment across time points, we identified co-expression patterns among genes using a weighted gene co-expression network analysis with the R package WGCNA v1.70-3 (Langfelder and Horvath, 2008). We used Pearson correlation coefficients and a scale-free network to create a dissimilarity matrix across the entire gene dataset that identified ‘modules’ of co-expressed genes with relative expression patterns that were summarized by module eigengene values. Modules were identified using a dynamic cut tree approach (with a minimum cluster size of 30) on a dendrogram created through hierarchical clustering of the dissimilarity data. The WGCNA analysis is complementary to the DESeq2 analysis above and has the benefit of reducing the number of statistical tests and risk of false positives or negatives associated with differential expression analysis of individual transcripts. To test for associations between module eigengene values and experimental factors, we used linear mixed effects models in the R package lme4 (Bates et al., 2015). The full model tested for each module included treatment (cold or control), time point, and a treatment×time interaction, with colony specified as a random effect: [Eigengene∼Treatment+Time+Treatment*Time+(1|Colony)]. We also tested alternative models without the interaction (Treatment+Time), or with only time (Time Only) or treatment (Treatment Only) as fixed effects. We identified the best-fit model for each module using the Akaike information criterion (AIC) and only retained modules that had a significant fixed effect of interest determined with the R package lmerTest v3.13 (Kuznetsova et al., 2017) (Satterthwaite approximation, P≤0.05).

We performed gene ontology (GO) enrichment analyses with the R package GOfuncR v1.10.0 (https://github.com/sgrote/GOfuncR). Hypergeometric tests were conducted for each time period by comparing a list of DEGs with a background list created using B. impatiens-specific GO annotations downloaded from the Hymenoptera Genome Database (Walsh et al., 2022). GO terms generated for each gene set of interest were assigned biological process, cellular component and molecular function ontology categories. Terms that were significant (P≤0.05) and had a family-wise error rate (FWER)≤0.1 were passed through Revigo (Supek et al., 2011; web-based interface at http://revigo.irb.hr/) to summarize and simplify the GO term list (allowed similarity=small, or 0.5). The R package GOplot v1.0.2 (Walter et al., 2015) was used to visualize the results. For the WGCNA module analysis, the gene sets within each module were evaluated for GO enrichment utilizing the GOfuncR-based approach while maintaining the same criteria and thresholds mentioned in the previous section. Results for the WGCNA modules were further visualized using the R package simplifyEnrichment v1.0.0 (Gu and Hübschmann, 2022), which clusters functional terms in a similarity matrix using the binary cut method. As gene ontology analyses can provide long lists of functional terms with high levels of redundancy, programs such as simplifyEnrichment attempt to use semantic measures to cluster terms into groups that have high within-group similarity and mutual exclusivity between groups. This can greatly simplify the results of functional analyses and aid biological interpretations. The simplifyGO() function of this package performs the clustering and produces a heatmap for each defined cluster, along with word cloud annotations composed of key functional terms associated with each cluster.

All workers recovered from chill coma within 10 min of returning to room temperature (Fig. S1). Only workers collected at the 10 and 30 min sampling points after the cold exposure showed significant transcriptional responses compared with controls (Fig. 2A; Table S4A,B). In the following, we refer to ‘upregulated’ or ‘downregulated’ transcripts as those with relatively high or low expression, respectively, in cold treatment samples relative to controls. At the 10 min time point, 371 significant DEGs were detected (231 upregulated, 140 downregulated; Fig. 2B; Table S4A), and at the 30 min time point, 51 DEGs were detected (36 upregulated, 15 downregulated) (adjusted P<0.05). The full set of significant genes is presented in Table S4A, but shared genes upregulated at both the 10 min and 30 min time points were glycogen phosphorylase (LOC100746613; Fig. 2C), serine/arginine repetitive matrix protein 1 (LOC100744870), facilitated trehalose transporter Tret1 (LOC100742930), and two uncharacterized loci (LOC100741559, LOC105680363), while shared downregulated genes were phospholipase A1 (LOC100747885; Fig. 2C) and monocarboxylate transporter 13-like (LOC100748182).

Genes involved in metabolic processes provided the clearest response following recovery from cold exposure. At the 10 min sampling point, eight significant biological process GO terms were identified (Table S4C), primarily associated with carbohydrate synthesis and metabolism (e.g. terms included ‘hexose metabolic process’, ‘carbohydrate metabolic process’ and ‘generation of precursor metabolites and energy’; Fig. 3). Significantly upregulated DEGs (Table 1, Fig. 3) associated with carbohydrate-related GO terms included phosphoenolpyruvate carboxykinase 1 (LOC105680266), glycogen phosphorylase (LOC100746613), pyruvate kinase (LOC100743581), enolase (LOC100746870), pyruvate carboxylase (LOC100749336) and facilitated trehalose transporter Tret1 (LOC100742930). The broader ‘catabolic process’ GO term was also overrepresented and included multiple upregulated genes associated with autophagy such as syntaxin-17 (LOC100748805), autophagy-related protein 9A (LOC100745297) and calcium-dependent protein kinase 16-like (LOC100741566).

Conversely, we observed a downregulation of genes from the 10 min sampling period that function in alternative metabolic pathways. Phospholipase A1 (LOC100747885), a gene critical in the breakdown of phospholipids into fatty acids, had a strong relative downregulation (log₂ fold-change of −4.28; Table 2), along with other genes relevant in lipid metabolism, including pancreatic triacylglycerol lipase-like (LOC100741829), short-chain specific acyl-CoA dehydrogenase, mitochondrial (LOC100741796) and probable enoyl-CoA hydratase, mitochondrial (LOC100746130). Biological process-related GO terms such as ‘purine nucleobase metabolic process’ and ‘carbohydrate derivative metabolic process’ were identified from the downregulated gene list and included DEGs such as phosphoribosylformylglycinamidine synthase (LOC100749333), GART trifunctional enzyme (LOC100745134) and xanthine dehydrogenase (LOC100745073; Fig. 3, Table 2). Notably, ‘small molecule metabolic process’ was also statistically overrepresented in the downregulated gene set’ however, this GO category is broad and also contains upregulated DEGs (Fig. 3).

Outside the scope of metabolic pathway regulation, there were other significant DEGs of interest identified from the upregulated 10 min gene list, including those relating to muscular function [uncharacterized locus (LOC100741447, D. melanogaster homolog myofililn), troponin C, isoallergen Bla g 6.0201 (LOC10075020) and troponin T (LOC100744757)], cytoskeletal restructuring [PDZ and LIM domain protein 4-like (LOC100748245), ankyrin-2 (LOC100744281) and uncharacterized locus (LOC100749880, D. melanogaster homolog supervillin)] and ion homeostasis [sodium/potassium-transporting ATPase subunit beta-1-interacting protein (LOC100744623)]. There was also a mix of upregulated and downregulated DEGs that play a role in the regulation of misfolded proteins during environmental stress [i.e. heat shock proteins (HSPs)] such as protein lethal (2) essential for life (LOC100742443; downregulated), heat shock protein 60A (LOC100748180; downregulated), dnaJ homolog subfamily C member 1 (LOC100748155; downregulated), 97 kDa heat shock protein (LOC100742698; downregulated) and Chaperone protein DnaJ-like (LOC105682145; upregulated).

Because of the smaller DEG set at 30 min, only a single significantly overrepresented term was recovered, which was related to lipase activity (molecular function: 1-acyl-2-lysophosphatidylserine acylhydrolase activity; Table S4C).

WGCNA clustered expression data into 20 co-varying eigengene expression modules (Table S2). Five of these modules had a significant relationship with a main effect (Table 3; Table S3A–E). In contrast to the clear temporal expression differences observed in the DESeq2 analysis, AIC supported ‘Treatment Only’ as the best-fitting model for all five modules. Nonetheless, this analysis allowed us to identify broader sets of genes and associated terms responding in parallel to the cold treatment without the need for gene-specific statistical tests. These modules contained many GO terms that overlapped with the DEG analyses discussed above, as well as several additional interesting terms (Table S4F). For example, differential gene expression of HSPs was observed in the 10 min sampling period, but GO terms relating to misfolding proteins were much more strongly represented in the WGCNA analyses for two modules; Module 19 was the smallest module (Table 3, Fig. 4), with 57 genes that were associated with ‘endoplasmic reticulum stress’, ‘protein folding’, ‘cellular response to topologically incorrect protein’, ‘chaperone binding’, ‘misfolded protein binding’ and ‘unfolded protein binding’. Module 12 (n=507 genes) was also associated with protein folding regulation (e.g. ‘protein glycosylation’, ‘protein localization to the endoplasmic reticulum’; Fig. 4) and included other GO terms represented in the 10 min downregulated DEGs, such as ‘carbohydrate derivative metabolic process’ and ‘catalytic activity’. Module 8 (n=361 genes) was strongly associated with energy production and mitochondrial function, including GO terms such as ‘quinone biosynthetic process’, ‘ATP metabolic process’ and ‘mitochondrial respiratory chain complex assembly’ (Fig. 4). The remaining significant modules (modules 1 and 15) represented a wider array of GO terms that included a mix of gene expression regulation and cell signaling terms (Fig. 4), including one notable observation for Module 1 (n=1283): GO terms important for neuromuscular function and synaptic maintenance that are expected to be crucial for chill coma recovery (Kelty et al., 1996). Out of all the modules, Module 1 also contained the most DEGs (n=104), which included top DESeq2 DEGs such as glycogen phosphorylase (LOC100746613), facilitated trehalose transporter Tret1 (LOC100742930) and sodium/potassium-transporting ATPase subunit beta-1-interacting protein (LOC100744623).

We investigated temporal patterns of gene expression during recovery from an acute chill coma-inducing event in the bumble bee B. impatiens to test the degree to which the recovery response is immediate and short-lived or might also comprise longer-lived transcriptional modifications that could prime individuals for future cold exposure (Teets et al., 2023). Our results indicate that B. impatiens workers undergo a rapid gene regulatory response during recovery from chill coma and after only 30 min exhibit no detectable DEGs from bees maintained at ambient temperatures, with most transcriptional changes observed at 10 min after removal from cold exposure, shortly after recovery from chill coma (Fig. S1). Interestingly, no prolonged differential gene expression was detectable at the 120 and 720 min time points. Overall, this suggests that B. impatiens workers mount a rapid and short-lived gene-expression response after cold exposure. These results imply that bumble bees may not cold harden, at least in relation to resisting and recovering from chill coma, though physiological experiments with additional timing and temperature treatments in conjunction with gene expression and proteomic analyses will be necessary to confirm this. That said, only one other study, to our knowledge, has reported potential cold-hardening responses in bumble bees: bees were more likely to survive exposure to −5°C if they had previously experienced a 1 h exposure to 0°C or if they were gradually ramped to −5°C (Owen et al., 2013). Mechanisms providing protection from lethal freezing (Colinet et al., 2018) may differ from those determining CT_min and recovery from chill coma (MacMillan et al., 2009; MacMillan and Sinclair, 2011; Overgaard et al., 2005; Teets and Denlinger, 2013; Teets et al., 2019). Further, the short-lived transcriptomic responses documented here may lead to longer lived alterations, for example in membrane properties or metabolite abundances that confer protection against subsequent cold exposure (Colinet et al., 2012; Hazel, 1995). Nevertheless, transcriptomic responses to cold exposure in other (ectothermic) insects often persist for hours to days with clear links to subsequent cold hardening (Overgaard et al., 2007; Teets et al., 2019); heterothermic bumble bees with high metabolic turnover may rely on alternative strategies to persist in highly variable thermal environments.

The identification of GO terms involved in carbohydrate synthesis and metabolism from the upregulated gene set indicates a shift in the metabolic activities of cold exposure treatment individuals as an integral part of the cold recovery process, consistent with other transcriptomic studies in Hymenoptera (Liu et al., 2020; Xu et al., 2017). One key feature in bumble bee cold tolerance is a capacity to contract thoracic muscles to generate heat (Beenakkers, 1969; Candy et al., 1997; Heinrich, 2004), which is powered by glucose metabolism (Surholt et al., 1991). The gene product of glycogen phosphorylase catalyzes glycogenolysis to produce glucose-1-phosphate within muscular tissue (Zhang et al., 2017), while facilitated trehalose transporter Tret1 is responsible for the mobilization of trehalose, the major disaccharide in insect hemolymph that is cleaved to provide glucose for flight and recovery from abiotic stress (Shukla et al., 2015). Trehalose is also a key cryoprotectant important for insect cold hardiness (Arrese and Soulages, 2010; Zeng et al., 2020). In our DEG dataset, upregulation of both genes persisted from 10 min to 30 min, suggesting the importance of carbohydrate utilization in bees recovering from cold exposure. Interestingly, protein phosphatase 1 regulatory subunit 3C-B, which stimulates glycogen synthesis and inhibits glycogen phosphorylase activity (Agius, 2015), was also upregulated in bees from 10 min but not 30 min. This might suggest that cold-exposed bees transition from glycogen synthesis to glycogenolysis within 10 min post-chill coma recovery, with the catabolism of glycogen for energy persisting for at least an additional 20 min. While other energy reserves such as lipids offer more metabolic energy per molecule compared with carbohydrates (Arrese and Soulages, 2010), glucose in thoracic muscles and circulating in the hemolymph is more easily accessible and faster to metabolize for energy when rapidly recovering from cold exposure. Thus, the presence of downregulated DEGs associated with lipid metabolism [e.g. phospholipase A1 (Richmond and Smith, 2011), pancreatic triacylglycerol lipase-like (Lowe, 1997)] was not surprising.

Corroborating the inference that energy derived from glucose plays a major role in the initial muscle activation, phosphoenolpyruvate carboxykinase 1 (Pepck1), an essential rate-limiting enzyme for gluconeogenesis that helps convert oxaloacetate into phosphoenolpyruvate (PEP; Spacht et al., 2018), was also strongly upregulated at 10 min. The upregulation of Pepck1 after acute cold exposure is consistent with findings in other insects responding to cold stress and preparing for diapause (Spacht et al., 2018). Other key upregulated genes include enolase, the gene product of which can process PEP into 2-phosphoglycerate in the gluconeogenesis pathway (Melkonian et al., 2019) or conversely serve in the glycolysis pathway by reversing the conversion of 2-phosphoglycerate back to PEP (Avilán et al., 2011), pyruvate kinase, which serves as a catalyst for the conversion of PEP into pyruvate (Hsiao et al., 2002; Wiskich, 1980), and pyruvate carboxylase, which carboxylates pyruvate to form oxaloacetate in the first step of gluconeogenesis (Melkonian et al., 2019). Module 8 in the WGCNA analysis comprised multiple GO terms associated with mitochondrial activity and ATP synthesis (Fig. 4), in which pyruvate plays an important role in the citric acid cycle (Akram, 2014). This trend is also observed in other insects that experience both seasonal and rapid cold hardening (Teets and Denlinger, 2013), suggesting that anaerobic respiration is involved in cold tolerance and subsequent responses across varying time scales (from acute to chronic exposure). Although the above genes do not represent all DEGs observed relating to gluconeogenesis and glycolysis (Table 1; Table S4A), they emphasize the importance of energy derived from glucose for fueling bumble bee thoracic muscles as they recover from chill coma.

Several genes directly associated with muscle and neuromuscular function were detected in our DESeq2 and WGCNA approach, including LOC100741447 (D. melanogaster homolog myofilin), which is required for striated muscle formation in insects (Qiu et al., 2005), and troponin C, isoallergen Bla g 6.0201 and troponin T, which facilitate muscle contractions via the binding of calcium (Ca²⁺) ions, similar to observations in cold-acclimated Drosophila (Bayley et al., 2020, 2018; MacMillan et al., 2016; Overgaard and MacMillan, 2017). Intracellular Ca²⁺/K⁺ ion homeostasis is vital for the induction of and recovery from chill coma (Overgaard and MacMillan, 2017) and genes associated with ion balance have been detected in transcriptomic and genomic scans for selection signatures across bumble bee populations exposed to cold temperatures (Jackson et al., 2018; Pimsler et al., 2020). The upregulation of one ion transport gene, sodium/potassium-transporting ATPase subunit beta-1-interacting protein (Hilbers et al., 2016), could be indicative of one mechanism of chill coma recovery by repolarizing the membrane potential and restoring neuromuscular function and excitability (Overgaard and MacMillan, 2017). A single gene does not confirm this phenomenon is a major component of the chill coma recovery in our experiment, and differential expression of additional ion transport genes would be required to achieve full muscular repolarization, but the detection is at least consistent with the hypothesis that repolarization may be involved in bees recovering from chill coma. Additionally, expression of the cytoskeletal binding protein encoding ankryn-2, suggested to play a role in the stability and modification of Drosophila neuromuscular junction synaptic connections (Cunha and Mohler, 2009; Koch et al., 2008), was also significantly increased in our bees recovering from chill coma. Other upregulated genes linked to cytoskeletal restructuring included PDZ and LIM domain protein 4-like and supervillin (a well characterized Drosophila homolog), which encode for proteins that can alter the actin cytoskeleton and influence cell motility and adhesion (Guryanova et al., 2011; Krcmery et al., 2010; Oh et al., 2003; Smith et al., 2010). Similarly, molecular indicators of cold-induced cytoskeletal reorganization have also been identified in other insect species (Des Marteaux et al., 2017, 2018; Kim et al., 2006; MacMillan et al., 2016; Teets et al., 2012), implying an important and possibly conserved nature of this response. From the WGCNA analysis and GO clustering approaches, Module 1 consisted of GO terms such as ‘cellular signaling’, ‘synaptic development and regulation’, ‘locomotion’ and other biological regulators and processes that respond to external stressors (Fig. 4; Table S3). An increase of eigengene values for this module (treatment relative to control; Table S4D) suggests that these functions and processes are actively being restored during the time of rapid recovery, coinciding with the differential expression of associated genes observed in muscular and cytoskeletal genes.

The upregulation of several genes involved in autophagy highlights the role of this process in cold recovery. Autophagy can break down glycogen to provide glucose for glycolysis, and also plays an important role in mitigating cellular damage sustained during cold shock (Kaur and Debnath, 2015). Autophagy has been observed as a crucial process in other insects exposed to cold acclimation treatments (Des Marteaux et al., 2017). We also detected differential expression of HSPs, although GO analysis only showed significant enrichment in the WGCNA modules. HSPs aid in protein stabilization and refolding after cold exposure (Joplin et al., 1990; King and MacRae, 2015; Lubawy et al., 2022; Rinehart et al., 2007; Vierling, 1991), and have been implicated in winter diapause (King and MacRae, 2015; Koštál and Tollarová-Borovanska, 2009; Rinehart et al., 2007), cold acclimation (Des Marteaux et al., 2017; MacMillan et al., 2016; Štětina et al., 2015) and cold shock (Colinet et al., 2010; Lu et al., 2016; Michaud and Denlinger, 2004) for a variety of insect species. Chaperone protein DnaJ-like and dnaJ homolog subfamily C member 1, both mediators of the Hsp70 family of proteins (Cyr et al., 1994), were differentially expressed in our cold-exposed B. impatiens, with DnaJ-like also identified as a DEG in B. vosnesenskii cooled to CT_min (Pimsler et al., 2020). Conversely, expression was reduced for protein lethal (2) essential for life, Heat shock protein 60A and 97 kDa heat shock protein, indicating that different HSPs are likely to be important in different ways during the cold exposure recovery process. Coupled with the decreased eigengene values (treatment relative to control; Table S4D) of Module 12, a module composed of ontologies and terms such as ‘unfolded protein response’, ‘glycosylation’ and ‘glycoprotein’, the differential expression patterns of HSPs may indicate issues with mechanisms regulating protein folding in cold-exposed bees. Module 19, which also had reduced eigengene values (Table S4D), was overrepresented by GO terms that also relate to protein folding regulation and, more importantly, ‘response to endoplasmic reticulum stress’ (ER stress). As an abundance of unfolded proteins aggregate in the endoplasmic reticulum, ER stress can elicit the unfolded protein response (Ron and Walter, 2007), which can be reversed through programmed cell death to purge afflicted cells. It is possible that our treatment bees experienced ER stress early in their recovery from chill coma, coinciding with the observed HSP expression patterns, the identification of Module 12 and the upregulation of autophagy-related genes.

We have attempted to highlight some of the most important DEGs and gene ontologies that have clear relevance to the chill coma recovery response. However, numerous other processes were detected. For example, we detected the metabolism of purines, which serve as key energy carrier molecules (e.g. ATP and GTP) and play an important role in facilitating proper cell signaling (e.g. cAMP and GMP) (Holland et al., 2011; Ramsey et al., 2010; Tang et al., 2019), and which have also been observed in cold-acclimated Drosophila suzukii (Enriquez and Colinet, 2019). The differential expression of xanthine dehydrogenase, a catalyst for oxidizing hypoxanthine and xanthine to NADH2 and uric acid, is particularly intriguing, as similar genes [xanthine dehydrogenase/oxidase-like (LOC100741462 and LOC117161180)] have repeatedly been detected in selection scans of bumble bees from different environments, although the function of this gene in Bombus adaptation remains elusive (Ghisbain et al., 2020; Heraghty et al., 2022; Pimsler et al., 2017).

This study provides insight into transcriptional changes that alter particular biological processes and molecular functions in B. impatiens recovering from cold exposure. Notably, differentially expressed genes identified in this study may help elucidate drivers of population variation in thermal tolerance across the geographic ranges of bumble bee species. Many bumble bee species are active early in spring and can fly and forage at colder temperatures than many other insects. Emerging spring gynes (reproductive females), early season workers, and males that have left the colony are especially prone to acute cold events as they venture away from the temperature-regulated colony. When bumble bees experience intermittent thermal fluctuations, the ability to thermoregulate, endure cold snaps and rapidly return to biological homeostasis is essential for individual survival. Our study reveals a rapid molecular response following acute cold exposure. High metabolic rates may be sufficient for bumble bees to quickly regain function without the need to maintain long-term transcriptional shifts in anticipation of subsequent cold events. That said, we recognize the possibility that long-term markers of cold hardening may in fact be present, but may exist at more subtle log fold-changes requiring larger sample sizes per sampling period to enhance detection probability.

For future research, capturing the intermittent transcriptomic profiles of bees recovering from chill coma, particularly between the 0 min and 10 min interval, may be beneficial in characterizing even finer temporal scale recovery responses, prior to, during and immediately after regaining central nervous system and neuromuscular function. Gene expression differences of alternative tissue types, such as brain and fat body, might also yield interesting results, given their roles in bumble bee sensory and motor function (Mares et al., 2005; Paulk et al., 2008), and energy storage/production (Alford, 1969; Arrese and Soulages, 2010; Costa et al., 2020; Röseler and Röseler, 1986), respectively. An experimental design incorporating individuals across castes (i.e. queens, workers and males) might also be an interesting research avenue as caste-specific responses may inhibit or enhance function given their respective life history strategies (Goulson, 2010). Ultimately, incorporating an integrative multi-omics approach, such as supplementing gene expression data with metabolomic and proteomic datasets and utilization of network-based integration approaches (Zhou et al., 2020), would be ideal to establish connections across multiple levels of biological organization and underlying molecular processes and pathways implicated in chill coma recovery.

We thank The University of Alabama College of Arts and Sciences and the University of Wyoming College of Agriculture, Life Sciences and Natural Resources. We are also grateful to The University of Texas at Austin Genomic Sequencing and Analysis Facility for preparing our TagSeq libraries, sequencing the libraries on their NovaSeq 6000 SR100, and performing quality control assessments.