Soil Metatranscriptomes Under Long-Term Experimental Warming and Drying: Fungi Allocate Resources to Cell Metabolic Maintenance Rather Than Decay

Earth’s temperature is rising, and with this increase, fungal communities are responding and affecting soil carbon processes. At a long-term soil-warming experiment in a boreal forest in interior Alaska, warming and warming-associated drying alters the function of microbes, and thus, decomposition of carbon. But what genetic mechanisms and resource allocation strategies are behind these community shifts and soil carbon changes? Here, we evaluate fungal resource allocation efforts under long-term experimental warming (including associated drying) using soil metatranscriptomics. We profiled resource allocation efforts toward decomposition and cell metabolic maintenance, and we characterized community composition. We found that under the warming treatment, fungi allocate resources to cell metabolic maintenance at the expense of allocating resources to decomposition. In addition, we found that fungal orders that house taxa with stress-tolerant traits were more abundant under the warmed treatment compared to control conditions. Our results suggest that the warming treatment elicits an ecological tradeoff in resource allocation in the fungal communities, with potential to change ecosystem-scale carbon dynamics. Fungi preferentially invest in mechanisms that will ensure survival under warming and drying, such as cell metabolic maintenance, rather than in decomposition. Through metatranscriptomes, we provide mechanistic insight behind the response of fungi to climate change and consequences to soil carbon processes.


INTRODUCTION
Genetic mechanisms and strategies behind fungal resource allocation in response to warming remain unclear despite being critical to understand ecosystem-scale carbon (C) dynamics under global climate change. For example, warming-induced soil carbon losses in high-latitude ecosystems are predicted to provide a positive feedback to global warming (Allison and Treseder, 2011). However, litter mass loss and microbial CO 2 emissions seem to decrease after a decade or more exposed to warming (Melillo et al., 2017;Romero-Olivares et al., 2017a,b). Simultaneously, warming changes the community composition of microbes, and these changes have the potential to affect ecosystem-scale C losses (e.g., Pold et al., 2015;Treseder et al., 2016;Morrison et al., 2019). How do these community changes and ecosystem-scale C losses feedback to global warming? An intimate look into the energetic investment of fungi under warming can help elucidate priorities on resource allocation under global climate change and facilitate predictions on the fate of soil C.
Microbes are mediators of biogeochemical cycles in soil. In particular, fungi are decomposition power-houses (Berlemont, 2017). In high-latitude ecosystems, fungi are the major decomposers of forest litter (Baldrian, 2017). Previous studies in high-latitudes reported that warming treatments significantly alter the fungal community, favoring soil saprotrophs and ectomycorrhizal fungi, while disfavoring yeasts (Treseder et al., 2016;Solly et al., 2017). In this biome, warming can also cause soil drying, owing to higher evapotranspiration rates (Allison and Treseder, 2011). Soil drying can contribute to these shifts in fungal communities under warming treatments (Allison and Treseder, 2008;Treseder et al., 2016). But to understand how fungi respond to warming and drying, and to better predict ecosystem-scale C dynamics under global climate change, we must explore beyond community composition (Crowther et al., 2014).
Laboratory studies on the physiological responses of the model fungus Neurospora discreta show that under warming, N. discreta adapts its physiology and invests energy differently to increase its chances of survival (Romero-Olivares et al., 2015). More specifically, after 1500 generations exposed to warm temperature and a labile C substrate, N. discreta invested more energy toward reproduction rather than growth and decreased its carbon use efficiency (CUE) (Romero-Olivares et al., 2015). But when exposed to warm temperature, and either a labile or recalcitrant C substrate, N. discreta was less sensitive to temperature shifts and most sensitive to C substrates (Allison et al., 2018a). Moreover, N. discreta isolated from different parts of the world-and presumably adapted to different mean annual temperatures-did not differ in their response to temperature (Allison et al., 2018b). Altogether, these studies expose the complexity of evolutionary forces in soil and the different responses fungi may have. Temperature, as well as warming-induced biotic and abiotic environmental changes (e.g., warming-induced drying), may be strong evolutionary forces (Bennett and Lenski, 2007;Allison et al., 2018a).
Studies on the metabolic changes of fungi under environmental stress, warming included, have focused on short-term responses (i.e., 30-60 min) (e.g., Gasch et al., 2000;Albright et al., 2018). However, in ecosystem ecology, "stress" is considered a chronic condition that extends for long periods of time with physiological expenses for microbes (Schimel et al., 2007). These expenses are usually associated with reallocation of resources to ensure survival (i.e., tradeoffs) (Bennett and Lenski, 2007). The rate-yield tradeoff states that "extra energy devoted to resource acquisition speeds metabolic rate but reduces the net yield of energy" (Frank, 2010). Indeed, in microbes, CUE can decline with warming, because respiration maintenance processes are more temperature sensitive than growth processes (Farmer and Jones, 1976;Mainzer and Hempfling, 1976;Hall and Cotner, 2007;Allison, 2014;Rodríguez-Verdugo et al., 2014). On the long-term, species that consistently invest in cell metabolic maintenance mechanisms may be selected for under stressful conditions brought by global warming.
Certain fungal traits associated with stress-tolerance may be favorable under warming and drying. For example, the presence of melanin is known to be associated with drought-tolerance and is often found in fungi that inhabit ecosystems exposed to highradiation, such as deserts (Harutyunyan et al., 2008;Belozerskaya et al., 2010;Azua-Bustos et al., 2012;Fernandez and Koide, 2013;Casadevall et al., 2017) Similarly, trehalose protects fungal cell membranes from desiccation, freezing, and heat shock (Treseder and Lennon, 2015). This trait is commonly found in arctic fungi (e.g., Gunde-Cimerman et al., 2003) where water availability is low, as well as in lichenicolous fungi which are known to be xerotolerant (Mittermeier et al., 2015). Under warming and drying, fungal taxa with stress-tolerance traits may thrive over those without those traits.
Exploring the fungal community composition and its traits, paired with physiological and metabolic profiles, can provide an integrative overview of the response of the fungal community to warming and drying. In addition, it offers insight into fungal adaptation strategies in the wake of climate change, and facilitates predictions on the fate of soil C. In this study, we used metatranscriptomes to get a mechanistic understanding of the allocation of resources in fungal communities exposed to a long-term warming treatment (which includes drying). We tested the following hypotheses:

Cell Metabolic Maintenance Hypothesis
Under the warming treatment, fungi allocate more resources into cell metabolic maintenance as a strategy to ensure survival compared to fungi under control conditions.

Decay Hypothesis
Under the warming treatment, fungi allocate less resources to decomposition compared to fungi under control conditions as a tradeoff for investing in cell metabolic maintenance.

Stress-Tolerance Traits Hypothesis
Under the warming treatment, fungal taxonomical orders that house taxa with stress-tolerant traits are more abundant compared to control.
Altogether, we hypothesized that under the warming treatment, there is an ecological tradeoff that favors the allocation of resources toward cell maintenance and stress-tolerance over the allocation of resources toward decomposition. This tradeoff is driven by the increase in abundance of stress-tolerant fungal taxa in response to warming and associate drying.

Field Site
The study area is located in a mature black spruce (Picea mariana) forest in Delta Junction, Alaska, United States (63 • 55 N, 145 • 44 W) on the Fort Greely military base. Here, the vegetation is dominated by black spruce and an understory of lichens, shrubs, and mosses. The annual precipitation is approximately 303 mm −1 and a mean annual temperature of −2 • C; the growing season starts around mid-May and ends in mid-September.
The warming experiment started in July 2005 as described in Allison and Treseder (2008). Greenhouses and neighboring control plots were established in pairs in a 1 km 2 area. Control plots were left under ambient conditions while greenhouses warmed passively using closed-top chambers. The top plastic panels of the greenhouses were removed at the end of the growing season and re-installed at the beginning of the growing season the following year to allow snow fall to enter the warmed plots. Gutters and tubing were installed to minimize drying by directing precipitation into the greenhouses during the growing season. On average, the air inside the greenhouses was 1.6 • C higher than in controls plots and the soil temperature at 5 cm depth was on average 0.5 • C higher inside the greenhouses compared to control plots. The warming treatment reduced soil moisture by 22%-on average-due to higher evapotranspiration (Allison and Treseder, 2008).

Sample Collection
We collected soil samples in the summer of 2015, after 10 years of the onset of warming. We used a 10 cm-tall corer to extract soil, homogenized 332 cm 3 of soil inside a plastic sterile Whirl-Pak R , collected approximately one gram, and immediately soaked it in 5 ml of LifeGuard TM Soil Preservation Solution (Qiagen, catalog 12868). Samples were kept in a cooler with ice for 24 h and then transferred to a −80 • C freezer until processed 1 week later.

Soil RNA Extraction and Sequencing
Samples were thawed on ice and centrifuged at 2500 × g for 5 min to remove LifeGuard TM Soil Preservation Solution. We proceeded to extract RNA using RNA PowerSoil R following the manufacturer's instructions with modifications by Baldrian et al. (2012). Samples were cleaned using RNA clean & concentrator TM -25 kit (Zymo, catalog R1017) and DNAse treated with Turbo-DNA free TM kit (Life Technologies, catalog AM1907). RNA was checked for quality on an Agilent 2100 Bioanalyzer at the University of California Irvine genomics high-throughput facility. Good quality samples were prepared for sequencing and sequenced at the Joint Genome Institute (JGI) (Nordberg et al., 2014). The JGI prepared paired-end libraries using Illumina RNAseq stranded library preparation kit following low and ultra-low input RNAseq with rRNA depletion protocols. Shortly, rRNA was removed from either 100 or 10 ng of RNA using Ribo-zero TM rRNA removal kit (Epicentre). Stranded cDNA libraries were generated using Illumina TruSeq R RNAseq stranded kit. The rRNA depleted RNA was fragmented and reversed transcribed using random hexamers and SSII (Invitrogen) followed by a second strand synthesis. The fragmented cDNA was treated with end-pair, A-tailing, adapter ligation, and 10 or 15 cycles of PCR. Sequencing was carried out in a HiSeq 2500 system. Sequencing projects are deposited at the JGI with project ids: 1107-496, -499, -504, -507, -509, -514, -519, and -520.

COG Cell Metabolic Maintenance Genes
For an overview of the investment of fungi in cell metabolic maintenance in control versus warmed plots, we analyzed metatranscriptomes using the JGI Integrated Microbial Genomes and Microbiomes (IMG/M) platform (Markowitz et al., 2005;Chen et al., 2016), specifically, we used the functional categories of Clusters of Orthologous Groups (COG) at 90% similarity with an E value threshold of 10 −5 . We chose COG categories because this classification provides information on the functional characteristics of a community of microbes by using a group of proteins found to be orthologous across lineages (Tatusov et al., 2000). We used the COG functional category of "metabolism" because it is the only category that includes proteins involved exclusively in cell metabolic maintenance, such as the transport and metabolism of amino acids, nucleotides, carbohydrates, coenzymes, lipids, and inorganic ions; energy production and conversion; and secondary metabolite biosynthesis.

CAZy Genes
Metatranscriptomes were quality trimmed by removing adapters with Trimmomatic (v 0.35) using ILLUMINA TruSeq3-PE adapters with sliding window 4:15 and dropping reads below 25 bases long (Bolger et al., 2014). Quality control of trimmed samples was carried out with FastQC (v 0.11.2) (Andrews, 2010). Metatranscriptomes were filtered with sortmeRNA (v 2.1) to remove all rRNA sequences (Kopylova et al., 2012). Filtering was done against bacteria and archaea 16 s-and 23 s-, as well as eukaryote 18 s-and 28 s-, and 5 s-and 5.8 s-SILVA databases (Quast et al., 2013). The rRNA-free metatranscriptomes were used for downstream analyses.
Four metatranscriptomes from control plots and four metatranscriptomes from warmed plots were selected to assemble a de novo reference meta-transcriptome with Trinity (v 2.3.2) (Haas et al., 2013). We used bowtie2 (v 2.2.7) to map reads to our reference meta-transcriptome (Langmead and Salzberg, 2013) and Samtools (v 1.3) for sorting and indexing (Li et al., 2009). We used Trinity toolkit (i.e., Trinity transcript quantification) to quantify abundance of transcripts and ran differential expression analysis with edgeR (v 3.3.2) using a p-value cutoff of 0.01. We used the Carbohydrate Active Enzyme (CAZy) database (Lombard et al., 2014) and Transdecoder (v 2.0.1) 1 to find coding regions within transcripts. We created an annotated database with Trinotate (v 3.0.0) 2 and hmmer (v 3.0) (Bryant et al., 2017) with an E value threshold of 10 −5 . The CAZy database has six different classes based on structurally related functional domains of enzymes that degrade, create, or modify glycosidic bonds. We used classes involved in catabolism, and excluded classes involved in anabolism (i.e., glycosyl transferases). Briefly, auxiliary activities (AA) is a class that has ligninolytic enzymes and lytic polysaccharide monooxygenases involved in the breakdown of lignin. Carbohydrate-binding modules (CBM) has enzymes that assist in the decomposition of cellulose and other types of complex carbohydrates. Carbohydrate esterases (CE) includes enzymes involved in the degradation of hemicellulose and pectin. Glycoside hydrolases (GH) is a wide group of enzymes that hydrolyze glycosidic bonds between carbohydrates or other non-carbohydrate molecules. Finally, polysaccharide lyases (PL) participate in the decomposition of acidic polysaccharides such as starch and chitin (Lombard et al., 2014). The full pipeline and scripts were deposited at https://github.com/adriluromero/ warming_AK.

Taxonomical Composition Analysis
We assessed fungal community composition in control and warmed plots using metatranscriptomes and M5NR database within MG-RAST (Meyer et al., 2008). Metatranscriptomes were deposited in the MG-RAST metagenomics analysis server with the following ID: mgm4761695.3, mgm4762031.3, mgm4762032.3, mgm4762033.3, mgm4762030.3, mgm4762034.3, mgm4762035.3, mgm4762036.3. We filtered results using an E value threshold of 10 −5 and kept only hits for fungi. We used data at the order level because previous work by Treseder et al. (2016) showed that variance in warming treatment response was greatest at this level. In addition (Randle-Boggis et al., 2016) recommends focusing on order when using M5NR in MG-RAST, because this taxonomic level yields that highest numbers of correctly annotated sequences.

Statistical Test
For all statistical tests we used R (R Core Team, 2018) with ranked data, since the unranked data did not meet assumptions of normality. We performed a linear mixed-effects model using the function lmer from the R package lme4 (Bates et al., 2015). In our model, plot was random effect, treatment was fixed effect, and transcription (i.e., gene counts for COG and CAZy) or relative abundance (for taxonomical order) was the dependent variable. COG metabolic maintenance subcategories, CAZy categories, or taxonomical order were categorical variables. We performed a Type III ANOVA in our model, to test for interactions; we used Satterthwaite's method within lmerTest package (Kuznetsova et al., 2017) in R. This method provides p-values for the F values produced with the lmer function. We used post hoc t-test to estimate significant differences in gene counts of each categorical variable between control and warmed plots. Since we ranked the data, all statistical tests were nonparametric. All statistical scripts and test results are deposited at github.com/adriluromero/warming_AK.

RESULTS
Under the warming treatment (which includes drying), fungi invested more in the transcription of cell metabolic maintenance genes compared to controls (Figure 1A). We accepted our first hypothesis, because the overall transcription of COG metabolic genes was significantly higher in the warmed treatment compared to control plots (Supplementary Table S1, P < 0.01). Moreover, there was no significant interaction between COG metabolic subcategory and warming treatment, indicating that COG categories responded similarly to the warming treatment (P = 0.97). Regardless of treatment, certain COG subcategories (e.g., energy production and conversion, amino acid transport and metabolism) were transcribed significantly more than others, as suggested by a significant main effect of COG subcategory (P < 0.01).
Conversely, under the warming treatment, fungi invested less in the transcription of genes for enzymes catabolizing glycosidic bonds (CAZy genes) compared to control temperatures ( Figure 1B). Accordingly, our second hypothesis was supported (P < 0.01). The various CAZy classes responded similarly to the warming treatment, based on a non-significant treatment by CAZy class interaction (Supplementary Table S2, P = 0.36).
FIGURE 1 | Cumulative average of transcripts in controls compared to the warming treatment for (A) metabolic clusters of orthologous groups (COG) and (B) carbohydrate-active enzyme classes (CAZy). Fungi in the warming treatment had on average significantly more transcripts of COG genes and significantly fewer transcripts of CAZy genes compared to controls. Bars are cumulative averages of four plots (n = 4) ± 1SE, P < 0.01 for treatment. Transcript counts of each sample and each COG subcategory and CAZy category are reported in Supplementary Tables S1, S2, respectively. Frontiers in Microbiology | www.frontiersin.org Across treatments, glycoside hydrolases was the most transcribed CAZy class, and PL was the least (main effect of CAZy class: P = 0.02).

DISCUSSION
In our study, we found that experimental warming treatments (which include drying) elicits an ecological tradeoff in fungal resource allocation which favors cell metabolic maintenance over decay (Figure 1). We also found that warming and drying causes an increase in abundance of fungal orders with stresstolerant traits (Figure 2), which may be drivers of the observed ecological tradeoff.
Previous studies have shown that cell maintenance costs increase with temperature (Joergensen et al., 1990), causing higher energetic demand to maintain cellular function (Steinweg et al., 2013). This increase in energetic demand results in increased death rates of microbes (Alvarez et al., 1995;Curtin et al., 2012). Consequently, increased temperature that leads to microbial death or resource allocation tradeoffs can potentially shape fungal communities and ecosystem-scale C processes. In our study, fungi under warming and drying invested resources in cell metabolic maintenance potentially to ensure their survival under global climate change.
Other studies have also documented microbial ecological tradeoffs. Treseder et al. (2011) found a tradeoff between recalcitrant C users and nitrogen availability; recalcitrant C users were more prevalent under environments high in nitrates and low in organic nitrogen. When two processes are energetically expensive, an ecological tradeoff will favor acquisition of one trait that favors survival at the expense of losing another one. Similarly, Romero-Olivares et al. (2015) found that the model fungus N. discreta adapted to elevated temperatures by producing more spores at the expense of biomass. Because producing spores is energetically expensive, N. discreta could not afford investing in both, spore and biomass production. Thus, investing in resources which will secure reproduction was an effective survival strategy. Treseder et al. (2016) studied community shifts in our field experiment by high-throughput sequencing of soil DNA. In the current study, we used a highly sensitive sequencing method based on soil RNA and found similar results; the Eurotiales, Pezizales, and Sordariales increased under warming and drying. Most of the orders we found at higher abundance under warming were taxa known to house lichenized, free-living, melanized, stress-tolerant, endophytic, and pathogenic fungi. Although some taxa displayed decreases in abundance in response to warming and drying, none of these were significantly different (e.g., Morteriellales, Cystofilobasidiales) ( Supplementary Table  S3). Interestingly, in our dataset, the Mortierellales, Mucorales, and other high-sporulating taxa were the most abundant taxa in the warming treatment and controls alike ( Supplementary  Table S3). Indeed, in culture isolation efforts carried out in these same soils by AL Romero-Olivares, Umbelopsis spp. Ambomucor spp., and Penicillium spp.-all high-sporulating taxa-were the most commonly isolated species in each treatment (data not published). This could indicate that the orders that were equally abundant in control and warming treatments were the result of spore banks in soils, which can be metabolically active (Novodvorska et al., 2016). Surprisingly, Basidiomycete fungi, known to play a key role in decomposition, were not significantly affected by the warming treatment. Specifically, no Basidiomycete orders decreased significantly under the warming treatment. In fact, the only Basidiomycete order that responded significantly to the warming treatment was the Malasseziales, which are known to house pathogenic yeasts (Figure 2 and Supplementary Table S3). Alternatively, Ascomycete orders were significantly affected by the warming treatment. Similar trends have been reported under long-term chronic nitrogen fertilization (Hesse et al., 2015).
Our findings are in line with mass loss litterbag experiments done in our field site. Here, we reported that resource allocation toward decomposition was lower in the warming treatment compared to controls. Romero-Olivares et al. (2017a) reported that mass loss and extracellular enzyme activity was lower in the warming treatment compared to controls. In contrast, they noted that the proportional loss of recalcitrant C in the warming treatment was greater than in controls. Differential expression analysis showed that AA genes, which are involved in the decomposition of recalcitrant C, tended to be transcribed at (nonsignificantly) higher rates in the warming treatment compared to controls (Supplementary Table S2 and Supplementary Figure S1). In contrast, CE and GH genes tended to be transcribed (significantly) more in controls than in the warming treatment. While the area of soil metatranscriptomes continues to develop, we suggest accompanying metatranscriptomes with field observations to draw ecologically relevant conclusions. FIGURE 2 | Fungal orders that displayed significantly different warming treatment response as percentage increase in relative abundance in response to warming. Bars are means ± 1SE of n = 4. * P < 0.05, * * P < 0.01, * * * P < 0.001. Transcript counts and relative abundance of all orders in this study are reported in Supplementary Table S3.
As the area of soil metatranscriptomics develops, sampling efforts should likewise increase. This goal is important, since fungi communities residing in different ecosystems might vary in their responses to warming and drying (Wallenstein and Hall, 2012). For example, deserts, savannas, tropical forests, and other particularly warm ecosystems may already be dominated by stress tolerant taxa. In those cases, shifts between cell maintenance and decay might not be as strongly evident. At the moment, sampling efforts are hampered by insufficient random access memory (i.e., RAM), even for high-performance computer clusters. The amount of sequencing data in metatranscriptomes is massive (i.e., ∼0.5 TB of sequencing data per sample), and analysis platforms require high RAM in order to successfully process the data. As such, we limited our sample collections to one sampling time in the middle of the summer. Thus, our results represent a snapshot of the time and conditions of that particular day and time in this specific ecosystem and generalization is impossible. In the future, transcript abundance should be analyzed throughout a day, week, months, and years. This is especially important because similar studies examining temporal changes of microbial activity have shown that fungal transcription in soil differs among seasons (Žifčáková et al., 2016). In fact, CAZy GH and AA genes are transcribed at higher rates in summer compared to winter (Žifčáková et al., 2017).
Moreover, preservation methods for RNA in soil samples are also critical. Flash freezing in liquid nitrogen may be the best preservation method for RNA in soil samples (e.g., Žifčáková et al., 2016, 2017). However, this approach was not feasible in our remote field site. Instead, we used LifeGuard TM Soil Preservation Solution. We caution that the preservation liquid may not reach every single cell in the soil fast enough to capture all RNA. In fact, the effect of preservation methods in forest soils has not been thoroughly compared yet (Tatangelo et al., 2014); when flash freezing in liquid nitrogen is not a feasible option, LifeGuard TM Soil Preservation Solution seems to be the most widely used option (e.g., Carvalhais et al., 2012;Che et al., 2017;Jurburg et al., 2017;Meredith et al., 2018;Rampadarath et al., 2018;Thompson et al., 2018). In addition, pipelines for processing and analyzing metatranscriptomes in the context of fungi and soil C processes are still developing (Baldrian and López-Mondéjar, 2014). Vast efforts are underway to develop pipelines for de novo metatranscriptomes (e.g., Haas et al., 2013) and improve fungal annotations (Grigoriev et al., 2014). Our work is the result of the existing information available in databases and pipelines for processing and analyzing transcriptomes without reference genomes.
In conclusion, our data suggest that under long-term warming and drying, fungal communities allocate resources toward cell metabolic maintenance at the expense of decomposition. In addition, fungi with stress-tolerance traits, such as melanized, xerotolerant, and lichenized fungi, seem to prevail under warming and drying. The metatranscriptome approach allowed us to pin-point, at the gene level, the mechanisms behind ecosystem-scale observations on soil C processes. We present a pipeline that paired with field observations, provides insight into the genetic mechanisms behind resource allocation efforts of fungi under long-term experimental warming and drying. It can help in predicting the fate of soil C under global climate change. In our study, soil metatranscriptomes provide a mechanistic snapshot of C-related ecological processes happening under climate change and offer insight into microbial resource allocation which are in-line with field observations.

AUTHOR CONTRIBUTIONS
AR-O conceived and designed the study, performed the research, analyzed the data, and wrote the manuscript. GM-C provided support with the data processing and helped to revise the manuscript. AL-L helped in the data processing and revising the manuscript. KT conceived and designed the study, and helped to revise the manuscript.

FUNDING
Financial support for this study was provided by the Consejo Nacional de Ciencia y Tecnología de México (CONACyT), UCMEXUS (scholarship no. 216015), and UCMEXUS-CICESE short-term research training fellowship to AR-O and the National Science Foundation grants DEB-1256896, DEB-1457160, and EAR-1411942.

ACKNOWLEDGMENTS
We thank Harry Magalam from the Office of Information Technology at UC Irvine and R. O. Romero-Olivares for scripting advice, J. Catalán-Dibene for the field assistance, C. I. Looby for the valuable feedback on the manuscript, the reviewers for the valuable comments that improved the quality of our manuscript, and Fort Greely and the U.S. Army for providing access to field sites.
TABLE S1 | Total transcripts of predicted metabolic clusters of orthologous groups (COG) genes in control and warmed samples (n=4).