Temporal regulation of temperature tolerances and gene expression in an arctic insect

The environmental conditions in the polar regions pose significant challenges to many terrestrial arthropods. In some of these regions, the microhabitat temperature can be highly variable and frequently reaches both cold and warm extremes depending on season (Danks, 2004; Bahrndorff et al., 2021b). Coping with such extremes is critical for survival and reproduction of arthropods in these regions, and they have developed a suite of adaptations that allow coping with temporal changes in thermal conditions.

Many studies on arthropods from polar regions have focused on adaptations to cold temperatures that enable survival during long and cold winters (Young and Block, 1980; Sømme, 1981; Block and Sømme, 1982; Sømme and Block, 1982; Cannon and Schenker, 1985; Worland et al., 1998; Block and Convey, 2001; Convey et al., 2003; Danks, 2004; Purać et al., 2008; Teets et al., 2012a). However, with increasing mean temperatures and temperature variability (IPCC, 2014), changes in both cold and heat tolerance over short time scales might also be critical for behavior, survival, and reproduction during the short summer period in polar regions (Rinehart et al., 2000; Worland and Convey, 2001; Shreve et al., 2004). Accordingly, several studies have examined the potential for plasticity through rapid cold hardening (RCH) in polar terrestrial arthropods, and have shown how RCH can improve survival and reproduction after exposure to cold extremes (Lee et al., 1987; Worland and Convey, 2001; Sinclair et al., 2003a,b).

Less is known about rapid physiological adjustments that can alter heat tolerance at short temporal scales, such as heat hardening (HH), despite arctic and sub-arctic species being frequently exposed to high daily temperature peaks at a microhabitat scale (Mølgaard, 1982; Bahrndorff et al., 2021b; Noer et al., 2022b). Moreover, most knowledge on RCH and HH in arthropods, including species from the Arctic and sub-Arctic regions, is derived from studies conducted under highly controlled laboratory conditions (but see Koveos, 2001; Kelty, 2007; Loeschcke and Hoffmann, 2007; Overgaard and Sørensen, 2008; Schou et al., 2015; Jensen et al., 2019; Teets et al., 2020). However, acclimation and hardening responses to diurnal and seasonal variation in temperature and other environmental variables have been shown to differ between laboratory and field conditions, and more field studies are needed to fully understand thermal plasticity under field conditions (Kristensen et al., 2008; Vasseur et al., 2014; Colinet et al., 2015; Sørensen et al., 2016; Jensen et al., 2019; Noer et al., 2022b).

The Greenlandic seed bug Nysius groenlandicus inhabits Arctic and sub-Arctic regions, where it is exposed to extreme variable thermal conditions (Bahrndorff et al., 2021a). The species is univoltine and completes its life cycle within the short summer season, where local microhabitat temperatures can reach >30°C during the day and subzero temperatures at night. We know that N. groenlandicus has a broad thermal tolerance range (Böcher and Nachman, 2001; Sørensen et al., 2019; Bahrndorff et al., 2021b; Noer et al., 2022b) and that it can rapidly increase acute heat tolerance when exposed to high and stressful temperatures under laboratory (Sørensen et al., 2019) and field conditions (Noer et al., 2022b). However, we lack a comprehensive understanding of temporal changes in thermal tolerances of this species in the field and the molecular responses associated with diurnal environmental variation.

Here, we investigated temporal changes in thermal tolerances of N. groenlandicus collected in the field at different time points during days differing markedly in mean microhabitat temperatures and temperature variation. We further used an RNA-seq approach to: (1) identify transcripts that, under field conditions, differed in expression levels across days and time points within days, and (2) identify transcripts that were affected by gradual heating and cooling under controlled laboratory conditions to obtain gene profiles of laboratory hardening responses. With these data, we investigated the association between field temperatures, heat and cold tolerance, and gene expression patterns. Further, we contrasted gene expression patterns in the field and the laboratory to investigate whether similar or contrasting responses were observed.

To investigate how the environment affects thermal tolerances and potentially induces alterations in gene expression levels in Nysius groenlandicus (Zetterstedt 1840), we conducted two experiments. First, a field experiment (see Field collection, below) was carried out to examine how natural temperature dynamics in concert with omnipresent environmental variables affect heat tolerance, measured as heat knockdown time (HKDT), and cold tolerance, measured as temperature at chill coma recovery [T_recovery; procedure described in Noer et al. (2022b)]. In addition, the underlying gene expression levels were investigated in individuals collected simultaneously in the field using an RNA-seq approach. Second, a laboratory experiment (see Laboratory experiment, below) was conducted to examine transcriptional responses in individuals that were exposed to controlled heat and cold temperature treatments in the laboratory.

The fieldwork was conducted in July–August 2018 in Narsarsuaq, Greenland (61.160°N, 45.424°W). Insect collections and thermal tolerance assays were performed as described in Noer et al. (2022b) with minor modifications. Briefly, adult individuals of N. groenlandicus were collected at four time points during the day for 5 days. The time of collection was dependent on weather conditions and the abundance of the species; for exact collection dates and times, see Table S1. Specimens were placed individually in 4 ml screw-cap glass vials (45×14.7 mm), the sex of individuals was then assessed by eye, and the vials were transferred as quickly as possible to the laboratory (all within 30–45 min of collection). During this time period, bugs were kept at ambient field temperatures as we wanted to test them immediately after they were taken from whatever field temperature they experienced at the time of capture. Immediately after returning to the laboratory, HKDT and T_recovery were scored using 20 females and 20 males for each assay, collection date and time. HKDT was scored by submerging the vials containing the bugs into a temperature-controlled water bath (PolyScience MX Immersion Circulator: MX-CA12E) pre-set at 48°C. The individuals were observed and stimulated with flashes of light and gentle taps on the vial caps with a metal rod. The time at which all movement had stopped was recorded as the HKDT. T_recovery was scored by submerging the vials into another temperature-controlled water bath (LAUDA Proline Edition X RP 1845-C) with a glycol–water solution maintained at −3°C, at which the bugs quickly went into a chill coma. The glycol–water was then gradually heated at a rate of 0.2°C min⁻¹. The temperature at which the bugs recovered from the chill-induced coma, i.e. first movement observed, was recorded as T_recovery. Abrupt exposure to high or low stressful temperatures as in the HKDT and T_recovery assays was used to reduce potential acclimation effects during thermal exposure. Further, T_recovery was used because N. groenlandicus, like most insects, is chill-susceptible, and entered a state of reversible coma at temperatures above those causing freezing and recovered gradually when returned to benign temperatures in a time-dependent manner (Bahrndorff et al., 2021a; Andersen et al., 2017).

For RNA-seq, female individuals were collected at the same four sampling times (08:00, 12:00, 16:00 and 20:00 h) on 2 days (i.e. 22 and 27 August 2018; see Table S1), representing a day with high and low temperature variation, respectively. At each collection time, six samples of 10 females were collected and transferred directly into ice-cold RNAlater-ICE Frozen Tissue Transition Solution that had been stored at −20°C (Invitrogen, cat. no. AM7030) to prevent degradation of RNA during subsequent storage at −20°C in Narsarsuaq for approximately 1 week. The samples were then transferred to our laboratory in Denmark, where they were stored at −80°C and later used for RNA sequencing.

The air temperature at the collection site was continuously recorded at 5-min intervals in the shade using Easylog USB data loggers (LASCAR Electronics, EL-USB-2⁺). The loggers were placed 15 cm above the soil surface to reflect the thermal environment at the top of the grasses where N. groenlandicus was most abundant and caught with a sweeping net.

The laboratory experiment was conducted using N. groenlandicus collected in 2019 from the same field site as used for the field experiment (see previous section). The specimens were collected from the field and transported to the laboratory in Aalborg, Denmark, within 2 days of collection. The bugs were placed randomly in large Petri dishes (145×20 mm) at a density of approximately 40 individuals per Petri dish, and fed with sunflower seeds, grasses from the field site, and a 10% sucrose solution in Eppendorf tubes plugged with small cotton balls. The Petri dishes with animals were kept under constant rearing temperature at 23°C and a 12 h:12 h light:dark cycle for 3 days before initiating the experiment.

Three temperature treatments were initiated simultaneously: gradual cold exposure, gradual heat exposure, and a control group kept at the constant rearing temperature of 23°C. The heat and cold treatments were initiated simultaneously with the sampling of the control group for RNA-seq for this group to serve as a control for the stress treatments. Females were placed in screw-cap glass vials and submerged into two different water baths at 25°C. The temperatures in the baths were then slowly increased or decreased at a rate of 0.2°C min⁻¹. Four replicates of five females were then sampled at two different temperatures during heat ramping (33 and 43°C) and two temperatures during cold ramping (13 and 3°C) and snap-frozen in liquid nitrogen. All samples from the laboratory experiment were stored at −80°C until subsequent RNA-seq analysis.

Both field and laboratory collected samples were prepared for RNA-seq using the following protocol. Total RNA was extracted from pools of 10 female N. groenlandicus using TRIzol^TM (Invitrogen, cat. no. 15596018) and PureLink™ RNA mini kit (Invitrogen, cat. no. 12183018A) following the manufacturer's protocol with few modifications. Pools of 10 females (15–20 mg tissue from field-collected samples and 8–12 mg tissue from laboratory-collected samples) were homogenized in 1 ml of TRIzol using a FP-120 FastPrep bead beater for 2×30 s at 6500 rpm and 5 s rest in between. The samples were incubated for 5 min on ice. Subsequently, 0.2 ml of chloroform was added, and the samples were shaken and incubated for 2–3 min. The samples were centrifuged for 15 min at 12,000 g at 4°C. The upper aqueous phase (∼400–600 μl) was transferred to new tubes and an equal volume (400–600 µl) of 70% ethanol was added and mixed. Total RNA was obtained following the protocol from this step but discarding the use of on-column purelink DNase treatment. RNA concentrations were examined with a NanoDrop 2000 (Thermo Fisher Scientific) and RNA quality was spot-checked on random samples with an Agilent 2200 TapeStation system.

From the RNA extract, mRNA was poly-A selected and reverse transcribed to cDNA libraries using a TruSeq Stranded mRNA Library Prep kit (Illumina, cat. no. RS-122-2101). Each library was given an individual 8-base barcode adapter to allow multiplexing and subsequently all libraries were pooled. The library pool was checked for the correct insert size (∼150 bp) on an Agilent 2200 TapeStation system and paired-end sequenced on an Illumina HiSeq 3000/4000 platform for 150 cycles.

Pre-processing and assembly were performed with CLC Genomics Workbench v. 20.0 (Qiagen, RRID:SCR_011853). RNA-seq reads were trimmed for quality (Phred score limit of 0.05) and sequences with ambiguous nucleotides (N) larger than 2 were removed. Homopolymers (Ts) were removed from the 3′ ends and poly-As were allowed to remain at the 5′ ends. Sequences shorter than 40 bp were removed to avoid sequences generated from technical artifacts and improve the de novo assembly.

The transcript assembly model was made using the de novo assembly module in CLC Genomics Workbench, using an implementation of the de Bruijn algorithm and 373,541,899 preprocessed RNA-seq reads as input. Bubble and word size were optimized automatically and contigs shorter than 200 bp were omitted. This resulted in 215,802 transcript models. However, many of these contigs were supported by only a few reads and are unlikely to represent true transcripts. Therefore, all RNA-seq reads (373,541,899) were mapped to the draft transcript model using a similarity parameter of 0.8 and a length parameter of 0.5. Overall, 346,731,148 (92.8%) were mapped to the draft model. Transcript models with an average read coverage of >1000 were included in the final transcript set that contained 7005 transcript models (FASTA files are available from figshare: https://doi.org/10.6084/m9.figshare.20310711). Although the final set contained only 3.2% of the naïve transcripts, 77.7% of the reads mapping to the draft transcript model mapped to the final transcript model, consistent with the final set being an inclusive representation of the true transcriptome. The final transcript model set was searched against the NCBI protein database specified for all organisms using BLASTx. Significant positive matches (e-value <10⁻⁴) were identified and gene ontology (GO) terms assessed with the CLC workbench plugin Blast2GO (Conesa et al., 2005). BUSCO v5 analysis was conducted on the 7005 transcripts in the 7kTranscript set using eukaroyta_odb10 as database and eukarota as lineage.

The temporal scale at which adjustments of thermal tolerances are affected by field temperature was examined using the following procedure. First, the mean temperature and temperature variability (coefficient of variation, CV) were extracted from the field-site microhabitat temperature recordings. These metrics were extracted in a range of time intervals, spanning 1 to 12 h, preceding the time points at which N. groenlandicus specimens were collected from the field (exact collection time points in Table S1) and tested for thermal tolerances. In detail, each metric was extracted in time ‘windows’ by moving the past time boundary back in time in 1-h intervals, thereby accumulating the data input for each measure. Next, we corrected the phenotypic data for long-term acclimation effects to separate the short-term (minutes to hours) effects of temperature on physiological adjustments from long-term acclimation effects (days to months). This was done by regressing mean HKDT and T_recovery on ‘test day’, whereafter the residuals from the regressions were extracted and used for further analysis. Thus, only hourly adjustments in thermal tolerances were considered further. Last, the residuals were correlated with the different temperature metrics at each of the extracted time intervals using Pearson's correlations based on the linear relationship between the two continuous variables (for further description of the method, see Noer et al., 2022a). The correlation between the raw mean HKDT and T_recovery were also examined by Pearson's correlation and are included in Fig. S1. Mean HKDT and T_recovery were normally distributed for both sexes, except for female mean HKDT, where data were slightly right-skewed. All microclimate data and phenotype data were analyzed using the software R (https://www.r-project.org/).

Differentially expressed genes (DEGs) were determined using the CLC Genomics Workbench. The CLC algorithm is an implementation of EdgeR and models each transcript using a generalized linear model assuming read counts follow a negative binomial distribution. Pairwise comparisons on TMM normalized read counts were performed on laboratory groups with ‘temperature treatment’ as the independent variable. In addition, pairwise comparisons were performed to test the effect of ‘sampling time’ within each sampling day, and finally pairwise comparison of the effect of ‘day’ within each sampling time. A false discovery rate (FDR)-corrected P-value of <0.05 and absolute log₂ fold change of >1 were used as significance criteria of DEGs. Enriched GO terms were determined using Fisher's exact test in the R package TopGO (https://bioconductor.org/packages/topGO/) with restricted output to GO terms with ‘Biological Process’ ontology only.

Finally, we compiled a list of specific genes associated with temperature stress responses in insects based on literature (Table S2) and used this to search the lists of DEGs for indications of stress responses.

The correlations between residuals of both HKDT and T_recovery and the mean temperature in the field varied over time for both sexes but differed in the strength and temporal scale of adjustments (Fig. 1). Generally, in females we observed stronger associations between HKDT and field temperatures than for males, i.e. stronger adjustments in tolerances with changes in environmental temperatures. Heat tolerance of males showed no association with short-term mean field temperature, though there was a weak positive relationship with variation in temperature (CV; Fig. 1A). In contrast, the correlation between T_recovery and field temperature was stronger for males (Fig. 1B). Hence, T_recovery for males was very sensitive to changes in field temperature and thus a high mean field temperature up to 8–10 h prior to sampling resulted in longer recovery times from chill coma. The association between HKDT, T_recovery and the mean temperature in the field ahead of collecting specimens was strongest for females in the time windows closest to the collection time. Thus, females were more heat tolerant (higher HKDT) when the field temperatures experienced just prior to testing were high, and less heat tolerant when field temperatures experienced before testing were low (Fig. 1A). However, these correlations were not significant for HKDT. The pattern of association between T_recovery and the field temperature was similarly highest at the time point leading up to collection of specimens for females with a gradual decline with time (Fig. 1B). Hence, chill coma recovery temperature was lowest (fast recovery) at cold environmental temperatures for females and vice versa at warm temperatures (Fig. 1B). These associations were significant for males, but for females only so when the data were not adjusted for ‘effects of day on T_recovery’ (Fig. S1). The association between HKDT and temperature variation (CV) peaked at 10–12 h prior to test (Fig. 1A), and similarly for T_recovery, but with slightly weaker associations (Fig. 1B).

To pinpoint transcriptomic responses to gradual changes in temperature, we used RNA-seq to first construct de novo transcript models (see Materials and Methods for details) and second, quantify transcript expression levels in female individuals that had been exposed to increasingly warm or cold temperatures in the laboratory within the range of temperatures occurring in the natural habitat of N. groenlandicus during summer. As expected, the RNAseq based transcript set is not complete because a cut-off was employed in order to avoid false transcript models at the expense of excluding true transcripts that are not highly expressed and therefore weakly supported. Overall, BUSCO analysis showed 99 complete BUSCOs (38.8%), of which 98 were single copy as expected, 12 (4.7%) partial BUSCOs and 144 (56.5%) missing BUSCOs.

The treatment causing the most substantial changes in expressed transcripts was warm temperature ramping from 23 to 43°C, which resulted in 292 differentially expressed (DE) transcripts of the 7005 transcript models relative to the 23°C control group (Fig. 2). A temperature increase from 23 to 33°C or a decrease from 23 to 13°C resulted in 122 and 106 DE transcripts, respectively, and the lowest response of 59 DE transcripts was observed at the low temperature ramping from 23 to 3°C. In the groups exposed to the most extreme temperature treatments, i.e. temperature ramping from 23 to 3°C and from 23 to 43°C, 24 DE transcripts were shared, of which five obtained a significant match during the BLASTx search of the protein database: dolichyl-phosphate beta-glucosyltransferase, hypothetical protein acetyl-L-carnitine (ALC)60_11283, cytochrome oxidase subunit I, and probable ATP-dependent rRNA helicase RRP3 and Cc8L18.2-like protein. The first three of these were found in all treatments together with nine other shared transcripts, but with no annotation (Fig. S2).

The GO enrichment analysis detected 28 enriched terms with a P-value <0.01 for the heat treatments (Table S3) and 10 enriched terms for the cold treatments. The top 10 GO enrichment terms for individuals exposed to heat ramping (combined 33 and 43°C treatments) were dominated by processes involved in protein function (protein folding, translation, peptide metabolic process, organonitrogen compound biosynthetic process, cellular amide metabolic process, amide biosynthetic process) as well as glycolytic energy processes (glycolytic process, pyruvate metabolic process and ATP generation from ADP). Together, this suggests that the main challenge for individuals exposed to heat stress is related to maintaining protein function and stability, at least in part by synthesizing novel proteins and providing energy to this effort. In contrast, the combined cold response (13 and 3°C) showed an abundance of GO terms related to DNA structure and conformational change (DNA conformation change, DNA geometric change, DNA duplex unwinding and chromosome organization), as well as providing energy to these processes during cellular respiration and energy production (cellular respiration, energy derivation by oxidation of organic compounds, aerobic respiration, ATP metabolic process and oxidative phosphorylation). This suggests that a challenge to be met during cold treatment is maintaining DNA dynamics, such as the ability to transcribe.

Targeted search for genes involved in the temperature stress responses (heat and cold stress, oxidative stress, desiccation, immune response, etc.) was performed using the list of genes compiled from the literature in Table S2. None of these candidate genes were found or significantly differentially expressed at 3 and 13°C compared with 23°C controls, whereas expression of 16 candidate stress-related genes was significantly increased at 43°C compared with 23°C, and one candidate gene was significantly increased at 33°C (Table 1). The majority of these were molecular chaperones and vitellogenin (vg) 1 and 2. The expression level at the high temperature of some of these genes was notable, including a 37-fold increase in the inducible Hsp70, followed by a 24-fold increase in vg2, a 15-fold increase in Hsp83 and a 9-fold increase in one Hsp90 variant.

Despite not finding stress-related genes in the other temperature contrasts, expression of some transcripts changed substantially. The ‘reverse transcriptase rna-directed dna polymerase from mobile element jockey-like protein’ and the ‘membrane component hypothetical protein ALC60_11283’ increased expression 256-fold and 393-fold in the individuals ramped to 13°C compared with the 23°C-exposed individuals, respectively. Likewise, three transcripts involved in oxidative phosphorylation in mitochondria, cytochrome b and two cytochrome oxidase subunit 1, responded strongly to ramping to 3°C, with 11-, 182- and 380-fold decreases in expression compared with the 23°C-exposed individuals, suggesting decreased ATP production at lower temperature.

The microhabitat temperatures varied markedly between the two experimental days on which samples were collected for RNA-seq (Fig. 3). The temperature on the day with low (day 1) temperature variation was characterized by the lowest recorded daily variation, whereas temperature on the day with high (day 2) temperature variation was characterized by a high temperature variation with an variation of 25°C; the average temperatures at the four collection times differed by 4.9°C on day 1 and 21.2°C on day 2. To understand the effect of natural field temperatures on the transcriptional response, we: (1) compared the transcripts and expression levels at each collection time point between the two days (Fig. 3), (2) quantified the number of transcripts differentially expressed between the collection time points within each day (Fig. 4) and (3) examined whether any of the transcripts expressed differentially in the laboratory ramped groups were also differentially expressed during the field sampling days (Fig. 5).

Differential expression examined across days at each time point showed that at all collection time points, except for the morning group (08:00 h), more transcripts were upregulated on day 2 (high variation) relative to day 1 (Fig. 3). Notably, the number of transcripts upregulated in the evening group (20:00 h) was much higher on day 2 relative to day 1 despite the temperature being almost identical for both days at this time point. In the morning, more transcripts were downregulated on day 2 relative to day 1, though not many transcripts were differentially expressed at this time. Downstream GO enrichment showed that transcripts that were differentially expressed between the two days in the morning and evening (08:00 and 20:00 h) were associated with cellular respiration and proton transmembrane transport related to the electron transport chain (Table S4). At the warmest time point (12:00 h), lipid storage/localization and ATP-synthesis dominated the GO terms, whereas gene expression/translation and cellular respiration transcripts were overrepresented in the afternoon (16:00 h).

The number of transcripts differentially expressed between time points in field-collected individuals varied little on the thermally stable day, ranging between 34 and 39 transcripts (Fig. 4A, day 1). Conversely, the number of DE transcripts changed substantially, for example, 198 upregulated transcripts (Fig. 4B), in individuals collected between the two time points by which the temperature changed the most (Fig. 3). A large proportion of these upregulated transcripts overlapped with the transcripts expressed in the laboratory groups exposed to gradual cooling and heating, namely individuals ramped to 43, 33 and 13°C (Fig. 5). Notably, between the morning and noon, 50 DE transcripts overlapped with the DE transcripts in the females that were heat ramped from 23 to 43°C in the laboratory, and 10 transcripts overlapped with DE transcripts in females exposed to 33°C heat ramping. The extensive overlap between field-induced transcription responses and transcriptional response to the 43°C heat ramping vanished between noon and the afternoon (16:00 h) (Fig. 5). This indicates that many of the pathways regulated at high temperatures in the laboratory are shared with responses to temperature fluctuations in the field, but at time points associated with lower temperatures in the field.

The ability of animals to cope with environmental stress partly depends on past experiences, yet studies that examine the effect of environmental factors influencing thermal tolerance and underlying physiological responses under natural conditions remain sparse (Koveos, 2001; Kelty, 2007; Loeschcke and Hoffmann, 2007; Kristensen et al., 2008; Overgaard and Sørensen, 2008; Schou et al., 2015; Gleason et al., 2017; Jensen et al., 2019). Here, we show that recent (within hours) thermal exposures drive changes in the transcriptome and thermal tolerances of N. groenlandicus. That is, past thermal experience, leading to hardening or acclimation responses, affects thermal tolerances in a trait- and sex-specific manner (Fig. 1). These results add to earlier work on this species showing that HKDT and T_recovery of field-collected individuals varied markedly across days and time of day (Noer et al., 2022b) (Fig. S2, Table S2). Together, these results suggest that female N. groenlandicus can change their cold and heat tolerance rapidly through plastic responses within a few hours, and that thermal performance is strongly linked to temperatures experienced in the period prior to testing. These results are in accordance with earlier work suggesting that heat hardening responses in this species can take place within 45 min post temperature hardening (Sørensen et al., 2019). Conversely, male heat tolerance did not respond strongly to short-term changes in field temperature (Fig. 1). Sexual dimorphism in thermal tolerances has been described by several studies, revealing that for many arthropod species, males have lower tolerances and plastic capacity than their female counterparts (Bradley, 1978; Willett, 2010; Schou et al., 2017; Ørsted et al., 2018; Noer et al., 2022a). Such dimorphism can lead to, and cause, sex-specific evolution of life-history traits including age of sexual maturity, body size, and longevity (De Block and Stoks, 2003; Stillwell and Fox, 2007; Stillwell et al., 2009; Rogell et al., 2014). Thus, increasing temperatures can pose a risk of disparate evolutionary trajectories for male and female bugs. Alternatively, this discrepancy between plasticity in heat tolerance for male and female N. groenlandicus might reflect behavioral differences in the natural environment, for instance females basking in the sun to complete their reproductive life cycle during the short arctic summer.

On a longer time scale, the correlation between mean field temperature and HKDT was negative for both males and females; thus, colder temperatures >12 h prior to test resulted in higher HKDTs. Studies using temperature fluctuations show that interruption from high temperature stress can have positive effects on heat tolerance (Renault et al., 2004; Colinet et al., 2006; Koštál et al., 2007; Ørsted et al., 2022). This is because alleviation from high temperatures allows repair of accumulated cell injuries and could explain the observed patterns in our study, suggesting a low thermal limit for injury accumulation for this species.

The associations between T_recovery and mean field temperature suggest that males, like females, are highly plastic for this trait. This is in accordance with hypotheses suggesting that cold tolerance is under stronger selective pressures than heat tolerance (Huey et al., 2003; Bahrndorff et al., 2007; Sunday et al., 2011; 2014). Species living in cold and highly seasonal environments in particular are expected to be exposed to more variation in lower critical temperatures compared with warm extremes and have evolved stronger phenotypic responses to cold conditions (Gaston and Chown, 1999; Sinclair et al., 2003a; Ghalambor et al., 2006; Sunday et al., 2011; Bahrndorff et al., 2021b).

With climate change, temperatures are expected to become more variable, and the frequency with which small organisms are exposed to both cold and warm extremes within short time frames is likely to increase (IPCC, 2014). Studies suggest that species at high latitudes have a broader thermal tolerance and are living in climates that are currently cooler than their physiological limits (e.g. Deutsch et al., 2008). However, polar arthropods are understudied in this context, and we therefore lack knowledge on this aspect. Further, living in temporally heterogeneous environments, e.g. for species such as N. groenlandicus, where frequent induction of energetically costly heat and cold responses can reduce energy allocated to other thermal defenses, reproduction and longevity, can potentially have negative impacts on fitness (Feder and Hofmann, 1999; Sørensen et al., 2003).

The gradual laboratory heat and cold ramping treatments led to mainly unique but also shared transcriptional responses (Fig. 2, Table 2). This indicates that the mechanisms operating during cold and warm temperature stress are partly different, which is also reflected in the trade-off between heat and cold tolerances found for female individuals. The most substantial increase in DE transcripts occurred in females heat ramped to 43°C, which induced three to six times more DE transcripts than the other treatments, and the high temperature treatment initiated part of the heat shock response (Fig. 2, Tables 1 and 2). This supports the observation that a major challenge of increased heat is maintaining protein functioning. The heat shock response is the cellular protective response induced by exposure to stressful high temperatures (and other environmental variables) and is found to be a conserved mechanism across animal and plant kingdoms (Feder et al., 1992; Feder and Hofmann, 1999; Sørensen et al., 2001; Bahrndorff et al., 2009; González-Tokman et al., 2020). During heat stress, levels of molecular chaperones increase in the cells, where they stabilize and refold unstable or denatured proteins and macromolecules. We found an increase in 12 different molecular chaperones, including four different transcripts of Hsp90, with heat ramping ending at 43°C. The largest fold change occurred in inducible Hsp70, which is one of the major heat shock proteins that has been studied extensively in relation to heat stress in arthropods (e.g. Feder et al., 1992; Dahlgaard et al., 1998; Krebs and Feder, 1998; Sørensen et al., 2001; Bettencourt et al., 2008; Bahrndorff et al., 2009; King and Macrae, 2015), followed by Hsp90 and Hsp83. No chaperones were found when ramping the temperature to 33°C, so apparently short-term exposure to this temperature in the laboratory is not within the thermal range inducing this stress response. Further, at 43°C, the expression of vg1 and vg2 was upregulated almost 27-fold. Vgs are apolipoproteins and egg-yolk precursors that transport lipids and proteins to the oocyte and are crucial for ovarian development (Wu et al., 2021). Thus, our data support that maintaining proper egg development, by maintaining a suitable amount of functional vg, even under increased temperature promoting protein denaturation, is a fundamentally important biological process, for which a protective response has evolved. This rapid increase in vg1 and vg2 might also aid females to complete their life cycle rapidly during the short, warm summer. Additionally, in bees, vg aids in building fat bodies for energy storage (Corona et al., 2007) and this could be of vital importance for insects in cold environments, where the energy reservoirs are critical for synthesizing cryoprotectants that protect against low temperatures during the freezing arctic nights (Storey, 1997; Arrese and Soulages, 2010). Finally, vg has antioxidant properties and can neutralize free radicals, which are byproducts of increased metabolism at higher temperature (Seehuus et al., 2006; Corona et al., 2007; Havukainen et al., 2013; Salmela et al., 2016). We cannot deduce whether vg has an impact on the lifecycle or heat tolerance from our results, but the different proposed functions of vg are not mutually exclusive, and vg might have several critical functions in this species under high temperature exposures.

We did not detect DE transcripts previously associated with cold hardening/acclimation in our study where we ramped temperatures from 23 to 3°C. The critical lethal temperature (CT_min) of N. groenlandicus ranges from −3.2 to 3.4°C (Bahrndorff et al., 2021b); thus, we would expect that the cold ramping would induce a cold shock response as the temperature approaches CT_min. It is not uncommon that the transcriptional cold response is relatively weak despite a marked functional phenotypic response to cold hardening/acclimation (Sinclair et al., 2007; Sørensen et al., 2007; Von Heckel et al., 2016; Königer and Grath, 2018). For instance, Teets et al. (2012b) found no upregulated transcripts during RCH in the flesh fly (Sarcophaga bullata). Despite this, the study found that metabolic pathways increasing gluconeogenesis (glucose synthesis), and synthesis of amino acids and polyols were affected by the cold shock treatment. Thus, our results support the notion that RCH does not require synthesis of new gene products (Sinclair et al., 2007; Teets et al., 2012b). Instead, posttranslational processes and intracellular signaling are sufficient to induce RCH, hence the response should be detectable at higher levels of biological organization such as metabolite level, as shown in Noer et al. (2022b; see also Teets et al., 2012b). In this context, our results indicate that DNA dynamics (conformational changes) rather than gene expression might be critical for adjustments of cold tolerance. Chromatin structure, histones and other DNA modifications participate actively in regulation of gene activity, and these might be crucial for maintaining gene regulation at low temperatures, e.g. chromatin and histone modifications help access of regulatory proteins necessary for regulation of cold responsive genes (Zeng et al., 2019).

In the field, expression patterns changed in individuals collected across time and days. We found a strong effect of environmental temperature on expression patterns (Figs 3 and 4) with more transcripts expressed differentially across all time points, except the morning, on the day with large temperature fluctuations compared with the day with less temperature variation (Fig. 4). Several abiotic and biotic factors co-vary with changes in temperature across time points, e.g. humidity, and might affect expression. However, our observation that many of the same transcripts are differentially expressed after exposure to high temperatures both when assessed in the laboratory and the field provide evidence that the transcriptional responses are indeed linked to temperature variation.

Even though a change in transcription at high temperatures occurs across many transcripts, it might not have a physiological impact. mRNA transcripts are labile and break down faster at higher body temperature. Hence, increased/constant transcription is required to maintain the same levels of transcripts as at lower temperature. Likewise, higher transcription levels might be a product of greater turnover of proteins. Thus, the protein levels might not change if the transcription and protein turnover are tightly coupled (Podrabsky and Somero, 2004). These observations are supported by the GO terms enriched during the laboratory heat ramping from 23 to 33 and 43°C combined, which suggests that protein function and stability were challenged and dominated by synthesis of novel proteins and energy. Downstream GO enrichment also showed that transcripts that were differentially expressed between the two days in the morning and evening (08:00 and 20:00 h) were associated with cellular respiration and proton transmembrane transport, indicating a high energy demand, which is in accordance with previous findings showing higher levels of various carbohydrates including glucose on the thermally variable day, especially in the morning and evening (Noer et al., 2022b). Between days, there was also a significant increase in transcripts related to lipid storage and localization, including increased expression of vg1 and vg2 on the warm day. Thus, vg also seems to have an ecological relevance in natural environments, but at time points where temperatures were lower than the exposure temperatures used in the laboratory.

Within days, by far the greatest response was observed when comparing morning and noon of the temperature-variable day, which are also the time points with the largest change in temperature. Changes in thermal responses can take place within hours but are also affected by prior exposure. It is thus currently unclear to what degree long-term exposure to temperature extremes also affects expression of transcripts. Between the morning and noon time points, several transcripts found to be differentially expressed in the laboratory were also differentially expressed in the field. For instance, between the earliest sampling point and noon, we found a small increase in expression of the chaperones Hsp83 and Hsp90. Although the expression levels were magnitudes smaller than observed in the laboratory, it is notable that the induction temperature of Hsp expression was much lower than in the laboratory. There could be several reasons for this. First, the temperature of induction of the heat shock response varies according to evolutionary adaptation (Gehring and Wehner, 1995; Boher et al., 2016) or acclimatization to the environment (Buckley et al., 2001; Karl et al., 2009). The individuals used for the laboratory experiment had been acclimated to 23°C for 3 days prior to the heat ramping. Based on the sensitive plastic acclimation responses of N. groenlandicus, this acclimation period might have been sufficient to raise the Hsp induction temperature. Second, previous studies have shown that during recovery from cold shock, the level of Hsps, including Hsp22, Hsp23 and Hsp70, in the cells increases as a consequence of cellular cold damage (Colinet et al., 2010; Teets et al., 2012b). This might be the case here as the night-time/early morning temperature is around the freezing point (Fig. 3). Third, multiple stressors act in concert in the natural environment, and a range of abiotic (e.g. temperature and desiccation) and biotic (e.g. predation, competition, starvation) stressors can affect the expression of Hsps in the field.

In summary, our results show marked thermal plasticity in N. groenlandicus, a widespread arctic arthropod, adapted to extreme thermal conditions. We show that the investigated population of this species continuously tracks the thermal environment and adjusts its thermal tolerance and transcript levels to match the environmental temperature. At the transcriptional level, individuals collected in the field on a day with high temperature variation showed a distinct expression profile compared with individuals collected on a day with low temperature variation; thus, increasing temperature variability leads to strong plastic responses which are likely adaptive, but which might also have energetic costs, especially in increasingly heterogeneous environments. Finally, we show that transcriptional pathways activated by high temperatures in the laboratory are shared but induced at time points with lower temperatures in the field, highlighting the importance of performing ecologically relevant field studies to understand how species cope with thermal variation in their natural habitat.

We thank Toke Thomas Høye, Department of Ecoscience, Aarhus University, for collecting Nysius groenlandicus from the field site in Narsarsuaq, 2019, for use in controlled laboratory experiments, Finn Thomsen, Mittarfik, Narsarsuaq and Storch Lund Polar-tut, Narsarsuaq, for access to field sites and laboratory facilities in Narsarsuaq, Greenland. The Results and Discussion in this paper are reproduced from the PhD thesis of Natasja Krog Noer (Noer, 2022).