Downregulation of the microbial protein biosynthesis machinery in response to weeks, years, and decades of soil warming


 How soil microorganisms respond to global warming is key to infer future soil-climate feedbacks, yet poorly understood. Here we applied metatranscriptomics to investigate microbial physiological responses to medium- (8 years) and long-term (>50 years) subarctic grassland soil warming of +6 °C. Besides indications for a community-wide upregulation of central metabolisms and cell replication we observed a downregulation of the protein biosynthesis machinery in the warmed soils, coinciding with a lower microbial biomass, RNA, and soil substrate content. We conclude that permanently accelerated reaction rates at higher temperatures and reduced substrate concentrations results in a cellular reduction of ribosomes, the macromolecular complexes carrying out protein biosynthesis. Later efforts to test this, including a short-term warming experiment (6 weeks, +6 °C), further supported our conclusion. Downsizing the protein biosynthesis machinery facilitates liberation of energy and matter, allowing microorganisms to maintain high metabolic activities and cell division rates even after decades of warming.


Introduction
Global temperatures and atmospheric carbon dioxide (CO 2 ) levels have increased steadily over the last 100 years (1,2). Feedbacks from the terrestrial carbon (C) cycle to the climate system represent a major uncertainty in the prediction of future temperatures (3). Soil microorganisms, including Bacteria, Archaea, Fungi, and protists, play a key role in the terrestrial C cycle by being responsible for the turnover of soil organic matter (SOM) and the subsequent release of CO 2 from soils to the atmosphere (2). Higher temperatures commonly lead to higher microbial activities, so global warming should accelerate the decomposition of SOM to CO 2 (4). On the other hand, SOM consists largely of microbial necromass and warming may stimulate microbial growth and thus necromass production (5, 6), promoting SOM formation. Whether soils will ultimately act as C sinks or sources, thus depend on how microorganisms respond to long-term warming. Nevertheless, microbial responses to global warming are currently poorly represented in Earth system models (7), which can, to some extent, be attributed to the challenges associated with studying and quantifying microbial responses to environmental change in complex soil environments (4).
Microbial responses to long-term soil warming may include i) quantitative and compositional changes of microbial communities, ii) physiological adjustments of individual microorganisms, including changes in growth and resource use, via transcriptional and translational regulation, iii) shifts in microbial interactions and emergent properties of the community, and iv) microbial adaptations by genomic rearrangements and evolutionary changes of the genetic code.
In one of the few truly long-term warming studies (Harvard Forest Warming Experiment: 26 y, +5 °C, midlatitude forest soil) (8), the observed warming effects on the microbial community included a decrease in observations of higher biomass-speci c microbial growth and respiration rates are related to changes in microbial gene expression.

Results And Discussion
Warming effects on soil physicochemical and biological properties We analyzed subarctic grassland soil samples from eight stable natural soil temperature gradients located in the same landscape at two adjacent sites being subject to sustained soil warming for years (8 y; MTW, medium-term warming site; 4 gradients) and decades (>50 y; LTW, long-term warming site; 4 gradients) (Fig. 1A), respectively. The analyzed samples comprised two distinct soil temperature regimes (A T , ambient soil temperatures and E T , +6 °C above ambient soil temperatures). To evaluate sample selection and characterize the individual sample groups we analyzed soil physicochemical and biological properties (Table S1). On average, total C, N and P (C tot , N tot , P tot ), dissolved organic C and P (DOC, DOP), and microbial C, N and P (C MO , N MO , P MO ), as well as RNA and water content were higher in A T than in E T ( Fig. 1B). Dissolved organic N (DON) and DNA content varied considerably between MTW and LTW, resulting in a lack of consistent differences between A T and E T for these variables. The pH was slightly higher in E T . These properties were individually tested for signi cant differences between de ned groups ( Fig. 1B, right panel, Table S2). Signi cant differences were mainly observed between the two temperature groups A T and E T . While C tot , N tot , and P tot , did not differ signi cantly between A T and E T (P adj < 0.097), pH, DOC and DOP, as well as most biological properties (C MO , N MO , P MO , and RNA content, respectively) did. Besides DON, no signi cant differences between the two non-warmed control groups (MTW-A T and LTW-A T ) were observed (Fig. 1B). These results largely mirror previous studies conducted at the same sites that also report lower substrate concentrations and microbial biomass contents in the warmed soils (15)(16)(17)(18), indicating high reproducibility and suitability of the experimental sites for a comparative analysis of medium-and long-term warming effects.
In addition to the lower RNA contents in the warmed plots (Fig. 1B) and a strong linear correlation between microbial biomass (using microbial C as proxy) and RNA per g dry weight soil (Fig. 1C, scatterplot, Table S3), we observed a trend towards lower RNA content per unit of microbial biomass in the warmed soils (P adj = 0.083). This was contrasted by a signi cantly higher DNA content per unit of microbial biomass (Fig. 1C, boxplots).
Taken together, our observations are in line with previous studies which characterized these warmed plots as signi cantly divergent from their non-warmed counterparts (15)(16)(17)(18)(19). The new observations also support the proposition that an initial acceleration of biotic activity due to warming leads to decreases in C, N, and P pools within the rst years after the onset of warming, followed by a decrease of microbial biomass (15). Warming also affected the ratio between nucleic acids and microbial biomass, an observation that was not made previously.
Warming effects on microbial gene expression Illumina paired-end sequencing of total RNA from 16 soil samples produced an average of 6.69 ´ 10 6 mRNA reads per sample (Table S4). Bacteria dominated the mRNA read pools (93.36-98.52%), followed by Eukaryota (1.28-6.45%), Archaea (0.5-1.6‰), and viruses (0.2-1.9‰). An average of 2.66 ´ 10 6 mRNA reads per sample was assigned to a molecular function (KO number) de ned in the KEGG Orthology database (21) (Table S4). These reads were also dominated by Bacteria, accounting for 90.82-98.47% of the functionally annotated mRNA reads (see Table S5 for a detailed comparison). We further analysed all functionally annotated mRNA reads assigned to a KEGG metabolic pathway or functional complex. Rare ed abundance counts (Table S6) indicate no signi cant effect of warming on functional richness (Fig. S1, Table S7) while a PERMANOVA analysis revealed a signi cant effect of warming on the composition of expressed KOs assigned to a KEGG metabolic pathway or functional complex (Fig. S2, Table S8). To identify the nature of this functional response we explored in more detail the functional assignments to KEGG categories and performed a differential gene expression analysis. KEGG offers a hierarchical structure with four layers, KEGG1 (broad functional categories), KEGG2 (functional subcategories), KEGG3 (functional pathways and complexes), KO (molecular functions, i.e. orthologous genes encoding proteins, enzymes, and enzyme subunits). In total, 3164 unique KOs assigned to a KEGG metabolic pathway or functional complex were detected in the metatranscriptomes. Of all functionally annotated mRNA reads, 78.9%, represented by 2732 KOs, were assigned to only one KEGG1 category (Fig.  S3), while the remaining mRNA reads were assigned to two or more KEGG1 categories. To avoid redundancy, we next explored the reads assigned to a single KEGG1 category and corresponding KEGG2 sub-categories. Their distribution across the soil temperature groups led to four major observations ( Fig.  S3): (i) Higher relative gene expression levels for major sub-categories such as Carbohydrate metabolism, Amino acid metabolism and Lipid metabolism, in the warmed soils (KEGG1: Metabolism). (ii) Lower relative gene expression levels for Energy metabolism in the warmed soils (KEGG1: Metabolism). (iii) Lower relative gene expression levels for sub-categories associated with protein biosynthesis, i.e. Translation, Folding, sorting & degradation, and Transcription in the warmed soils (KEGG1: Genetic information processing). (iv) Higher relative gene expression levels for Replication & repair in the warmed soils, albeit most pronounced in MTW-E T (KEGG1: Genetic information processing). Only within the KEGG1 category, Metabolism, a substantial fraction of the reads was assigned to two or more KEGG2 categories. Within the other KEGG1 categories, the fractions of ambiguous assignments were <1% of all KEGG2 assignments (Table S9). The very same observations (i.e. i-iv) were made if all ambiguous assignments were included (Fig. S4).
These initial ndings pointed towards substantial overall changes in the most basic cellular functions such as protein biosynthesis, replication, energy conservation, and central metabolisms. Next, we employed an abundance-pattern lter (see Materials and Methods) to identify KEGG3 pathways and functional complexes exhibiting putative warming-induced differential expression. More than 40% of all mRNA reads assigned to a KEGG metabolic pathway or functional complex were associated with KEGG3 categories that exhibited warming-induced differential abundance patterns (Fig. S6, Table S10). Out of 56 KEGG3 categories with a visually distinct pattern, 13 showed lower relative gene expression levels in E T than A T , while 43 KEGG3 categories showed higher relative gene expression levels in E T than A T . The comparatively high relative abundances of eukaryotic mRNA reads in the LTW-A T samples (4.7% ± 1.8 SD) compared to MTW-A T , MTW-E T , and LTW-E T (1.5% ± 0.6 SD) may skew or mask biologically and ecologically relevant changes. Thus, we complemented the abundance-pattern analysis with a differential gene expression (DGE) analysis on the KO level (orthologous genes), excluding eukaryotic assignments. A total of 787 KOs were included in the DGE analysis (see Materials and Methods), out of which 66 KOs (i.e. 8.4%) were differentially expressed (false discovery rate (FDR) < 0.05) when comparing all A T with all E T soils (Table S11). Comparing gene expression pro les of MTW-A T and LTW-A T revealed no signi cant differences between the two non-warmed control groups (Table S11). We then assigned the differentially expressed KOs to the hierarchical KEGG structure (Table S12). The majority of the differentially expressed KOs fell within KEGG3 categories associated with protein biosynthesis, i.e. Ribosome, Protein export, RNA polymerase, and RNA degradation. In addition, several KOs assigned to these functional categories with an FDR between 0.05 and 0.1 shared the same expression pattern. The summarized expression of genes for all of these functions was lower in the warmed soils (Fig. S7) indicating a down-regulation of the protein biosynthesis machinery in warmed soils. Furthermore, several differentially expressed KOs were assigned to the Energy metabolism KEGG3 category Oxidative phosphorylation, the lower summarized expression of these genes in the warmed soils (  (Table S13-S15).

Down-regulated enzymes and enzyme complexes
The protein biosynthesis machinery. Due to the high number of differentially expressed KOs encoding enzymes involved in protein biosynthesis (Fig. S7) we decided to analyze this aspect of cellular metabolism in more detail. Our analysis included (i) the bacterial DNA-directed RNA polymerase, responsible for the transcription of DNA into rRNA (ribosomal RNA), mRNA (messenger RNA), and small non-coding RNAs, (ii) the bacterial ribosome, consisting of ribosomal proteins and rRNA and the site of protein synthesis (translation of mRNA into peptides), (iii) bacterial RNA degradosomes, responsible for the degradation and recycling of RNAs, (iv) bacterial protein folding complexes, responsible for the proper folding of peptides into mature proteins, and (v) the bacterial Sec dependent pathway of posttranslational translocation of proteins across membranes ( Fig. 2A, drawings). The relative abundances of mRNA reads encoding RNA polymerase subunits, ribosomal proteins, GroEL and DnaK (protein folding), and Sec translocase subunits (protein export) were consistently lower in the warmed soils ( Fig. 2A,  boxplots). For all these enzyme complexes we found a substantial fraction of KOs with signi cantly lower relative expression levels in the warmed soils, the pattern being largely consistent between MTW and LTW ( Fig. 2A, table, Table S16). KOs encoding ribosomal proteins accounted for most of the differentially expressed KOs, with signi cantly lower relative expression levels in the warmed soils (overall 19 KOs). Within MTW, 22 ribosomal proteins showed signi cantly lower relative KO expression levels in the warmed soils, supported by eight ribosomal proteins with lower but not signi cant relative KO expression levels in the warmed soils (FDR 0.05-0.1). In LTW, we observed signi cantly lower relative KO expression levels in the warmed soils for seven ribosomal proteins. Supporting this observation, were 10 ribosomal protein encoding genes with lower, but not signi cantly different, relative KO expression levels in the warmed soils (FDR 0.05-0.1) (Table S16). Two ribosomal proteins showed signi cantly higher relative expression levels in the warmed soils of LTW, indicating that not all ribosomal proteins are regulated by the same feedback mechanisms (22). The expression of Sec subunits (protein export) showed highly consistent temperature responses in MTW (Table S16). While less clear in LTW, the summarized mean relative mRNA abundances being overall lower in the warmed soils ( Fig. 2A, Protein export-boxplot). A contrasting temperature response was shown for RNA degradosome sub-units, with overall higher, but not signi cantly different, expression in the warmed soils relative to the ambient soils ( Fig. 2A). Taken together our results suggest that the entire protein biosynthesis machinery, including enzyme complexes involved in transcription, translation, and protein processing, may be downregulated in the microbiomes subject to soil warming of +6 °C above ambient temperatures, irrespectively of the duration of warming (8 y vs. >50 y). Notably, the lower relative abundances of mRNA reads encoding ribosomal proteins suggests a downregulation of ribosome synthesis in the warmed soils. This was supported by the lower total RNA content per unit of soil ( Fig. 2A, RNA-boxplots) and total RNA content per unit of microbial C (Fig. 1C) in both medium-and long-term warmed soils.
Energy metabolism. We observed lower expression of KOs associated with the Energy metabolism subcategory Oxidative phosphorylation in the warmed soils (Fig. S7). Oxidative phosphorylation is the nal step in aerobic respiration (Fig. 2B, drawings). Electrons harvested from the oxidation of carbohydrates are transferred via a membrane-bound electron transport chain, O 2 acts as terminal electron acceptor, and the generated proton-motive force is used by the ATP synthase to phosphorylate ADP and produce ATP, the major energy currency of all cells. Even though respiratory electron transport chains vary among Bacteria (23), we observed lower mean relative abundances of multiple common complexes including ATP synthase in the warmed soils ( Fig. 2B, boxplots). Of all investigated membrane bound complexes, only succinate dehydrogenase, and only in MTW, did not exhibit a warming-dependent expression pattern. However, succinate dehydrogenase represents a non-pathway-speci c enzyme (21), offering a potential explanation. Despite few genes being signi cantly differentially expressed (Fig. 2B, table), the large number of genes encoding metabolically connected membrane-bound enzyme complexes displaying a lower mean relative expression level in the warmed soils suggests that downregulation of the oxidative phosphorylation complexes may be a microbial response to soil warming.

Up-regulated enzymes and enzyme complexes
Central (carbohydrate) metabolisms. KEGG2 and KEGG3 categories with overall higher relative expression levels in the warmed soils were to a large extent observed within the KEGG1 category Metabolism (Fig.  S3, S6, and Fig. S7), prompting us to analyze community transcript pro les focusing on central (carbohydrate) metabolisms (Fig. 3A). We observed only few differentially expressed KOs but all had higher relative expression levels in the warmed soils. Individual KEGG metabolic pathway maps (21) re ecting KEGG3 categories (e.g. Pyruvate metabolism, Glycolysis / Gluconeogenesis, and Methane metabolism) overlap and the few differentially expressed KOs mainly represent KOs assigned to several KEGG3 categories. Nevertheless, a few noteworthy observations were made (Fig. 3A).
Enzymatic reaction steps leading from pyruvate to Fatty acid biosynthesis via acetyl-CoA showed higher relative expression levels in the warmed soils while steps leading from pyruvate to Oxidative phosphorylation showed the opposite trend, irrespectively of the duration of warming (8 y vs. >50 y) (Fig.  3A). This matches our detailed analysis on Fatty acid biosynthesis (Fig. S8) and Oxidative phosphorylation (Fig. 2B), as an upregulation of fatty acid biosynthesis would require an increased ow of C from the central intermediates, acetyl-CoA and pyruvate, while a downregulation of Oxidative phosphorylation might be linked to a downregulation of succinate supply. Some integral pathways exhibited inconsistent temperature responses. For example, the rst steps in glycolysis (from α-D-and β-D-glucose to β-D-fructose-6P) showed higher relative expression levels in the warmed soils and include differentially expressed KOs (Fig. 3A). However, most subsequent steps did not reveal a clear temperature trend. A previous in situ study on samples from the same LTW plots, taken one summer earlier, reported signi cantly higher biomass-speci c organic C uptake rates (mg C g -1 C MO h -1 ) and higher biomassspeci c respiration rates (mg C g -1 C MO h -1 ) in the warmed soils (E T , +6 °C) compared to their ambient counterparts (16). However, unlike the very consistent patterns within protein biosynthesis described above, our more detailed investigation of central (carbohydrate) metabolisms did not reveal su ciently convincing patterns to conclude that a common microbial upregulation of central (carbohydrate) metabolisms is occurring due to soil warming and ambiguous KOs complicate interpretation.
Cell replication. We observed higher relative expression levels for amino acid and nucleotide metabolisms (Fig. S7), Glycerolipid metabolism ( Fig. S7), Fatty acid biosynthesis (Fig. S8), and Peptidoglycan biosynthesis (Fig. S6), important pathways in the production of major biomolecule building-blocks and cell membrane and cell wall biosynthesis. Taken together these results prompted us to have a closer look on cell replication (Fig. 3BC). Bacterial cell replication can be divided into two closely interacting cycles, (i) the DNA cycle, i.e. DNA replication and chromosome segregation, and (ii) the division cycle, i.e. cytokinesis and cell separation (24). FtsZ is a key enzyme and regulator of bacterial cell division and an upregulation of its synthesis has been shown to correlate with the formation of the Z-ring, a ring-like structure at the division site constricting during cell division (24,25). We saw indications for higher average expression levels of DNA polymerase III subunits and DNA replication initiation factors (Fig. 3B) as well as higher average expression levels of FtsZ in the warmed soils (Fig. 3C), but these patterns were not signi cant. Cellular DNA content depends on the stage of growth, usually being higher in actively replicating cells (26). The total DNA content per unit of microbial C was signi cantly higher in the warmed soils (Fig. 1C), further indicating a higher relative abundance of cells being replicated in the warmed soils. This interpretation is strongly supported by a previous study reporting signi cantly higher biomassspeci c microbial growth rates (mg C g -1 C MO h -1 ; measured via DNA replication rates) in samples from the same LTW plots (E T , +6 °C), taken one summer earlier, compared to their ambient counterparts (16).
However, extracellular or relic DNA stemming from dead microorganisms, which can represent up to 40% of prokaryotic and fungal DNA in soil (27), could also have contributed to the higher DNA:biomass ratios in the warmed soils.
Decoupling of microbial community structure and functional temperature response A simple explanation for changes in gene expression between ambient and warmed microbiomes could be a compositional change in the active microbial communities as a response to warming. Therefore, we analyzed the taxonomic composition of mRNA read pools in detail, combining diversity and variance analysis with abundance-pattern and DGE analyses. An average of 6.02 ´ 10 6 mRNA reads per sample was taxonomically classi ed (Table S4). As described above, mRNA read pools were largely dominated by Bacteria (93.36-98.52%; Table S5, Fig. S9). Overall, more than 1,000 different prokaryotic, eukaryotic, and viral families were detected in the metatranscriptomes (Table S6). Mean family richness did not differ signi cantly between warmed soils and their ambient counterparts (Table S7), albeit showing a trend towards a reduced richness in the warmed soils (Fig. S1B). Furthermore, a PERMANOVA analysis revealed no signi cant effect of warming on the taxonomic composition of mRNA reads assigned on family level (Fig. S2B, Table S8), contrasting the observation of a signi cant effect of warming on the composition of expressed KOs (Fig. S2A). Correspondingly, a differential expression analysis of all mRNA reads assigned on family level did not reveal any overall signi cant differences in the taxonomic composition between ambient and warmed soils (A T vs. E T , Table S19). Repeating the differential expression analysis on subsets representing bacterial families and/or considering only functionally and taxonomically classi ed reads identi ed up to three families with differential abundances between A T and E T (Table S19). However, these taxa accounted for less than 0.05% of the whole communities, a rming that these minor changes in the taxonomic composition were not responsible for the gene expression changes. Furthermore, no signi cant differences in the community composition of the mRNA read pools between the ambient sites (MTW-A T vs. LTW-A T ) that could potentially explain the lack of taxonomic changes were observed (Table S20). In contrast, signi cant differences between the warmed sites were observed (MTW-E T vs. LTW-E T , Table S20) prompting us to investigate taxonomic composition of mRNA read pools, abundance pattern, and changes in relative abundances of both sites individually (Fig. S9). While the taxonomic composition of mRNA reads assigned on family level did not change in response to medium-term warming (MTW-A T vs. MTW-E T ), it did in response to long-term warming (LTW- (Table S21). Up to 44 families, 5 orders, 3 classes, and 6 phyla showed a differential abundance between LTW-A T and LTW-E T . These taxa accounted for less than 2% of the whole communities (Table S21) and up to 75% of the families displaying a differential abundance between LTW-A T and LTW-E T represented fungal families with lower relative mRNA read abundance in the warmed soils (Table S21). Thus, the earlier-mentioned functional changes observed within the bacterial community in response to warming (Fig. 2, Fig. 3) were not related to a shift in the taxonomic composition (of mRNA reads assigned on family level). We attempted to increase the taxonomic resolution by assembling the metatranscriptome reads (see Supplementary Materials). Our approach using rnaSPAdes (28), however, did not provide a su cient number of long mRNA contigs (<10% of functionally annotated mRNA contigs were longer than the unassembled reads). Also, previous studies employing amplicon sequencing of rRNA genes did not indicate warming-induced changes of the microbial community composition at genus and operational taxonomic unit (OTU) level (15,16,29). Taken together these results provide evidence for a decoupling of microbial community structure and functions, as recently indicated in studies of different ecosystems including short-term warmed peat soil (30)(31)(32)(33)(34). In line with this, an extensive screening of all phyla, classes, orders, and families with a mean relative abundance ≥1‰ in MTW-A T , MTW-E T , LTW-A T , or LTW-E T revealed that on average nearly two thirds of all taxa (63.7%) showed taxon-speci c expression patterns re ecting the overall temperaturedependent expression patterns of central carbohydrate metabolisms, energy metabolism, and protein biosynthesis (Fig. 3D, Fig. S11). However, the percentage of taxa exhibiting the pattern varied between 43 and 95% depending on the KEGG3 category and the warming duration and was clearly highest for the KEGG3 categories related to protein biosynthesis and energy metabolism.
Downregulation of the protein biosynthesis machinery as common physiological response of microorganisms to soil warming? Downregulation of the protein biosynthesis machinery across the bacterial community (Fig. 3D), especially of ribosomal proteins, was the most pronounced response of the microbial communities to soil warming of +6 °C observed in this study ( Fig. 2A). Combined with the reduced total RNA content per unit of microbial C in the warmed soils (Fig. 1C), these data suggested that the ribosome content per bacterial cell is lower in the warmed soils than in their ambient counterparts, as rRNA can account for >90% of the total cellular RNA content (35). Besides higher temperatures, the warmed ForHot soils are also characterised by a lower C, N, and P content (16-18). Starving Escherichia coli and Salmonella spp. cells reduce their ribosomal content (36, 37), suggesting that a downregulation of the translation machinery, is metabolically favorable in nutrient-limited ecosystems (38). Since the translation machinery accounts for up to 40% of total cellular proteins (39) and protein biosynthesis is the costliest type of macromolecular synthesis (40-42) a reduced number of ribosomes would furthermore reduce the energy costs (ATP) related to cell maintenance and replication. The indicated warming-induced downregulation of membrane-bound energy harvesting complexes including ATP synthase (Fig. 2B) might therefore be directly linked to the downregulation of the protein biosynthesis machinery. The described downregulation of the protein biosynthesis machinery and energy harvesting complexes do not necessarily re ect a downregulation of the processes itself (i.e. protein biosynthesis and energy generation). Protein and ATP synthesis rates (per ribosome and per ATP synthase complex, respectively) are also highly affected by temperature, with higher temperatures resulting in higher rates. This ts to the higher relative expression levels of genes for nucleotide, amino acid, and inorganic phosphate (Pi) supply, to the protein and ATP synthesis machineries (Fig. S7 and Fig. 2B). Consequently, the warming-induced downregulation of the protein biosynthesis machinery and energy harvesting complexes is not contradictory to a warming-induced upregulation of cell replication and increased biomass-speci c growth rates (16). It is long known that increased temperatures accelerate mRNA synthesis and protein synthesis rate per ribosome (peptide chain elongation rate) in the model organism E. coli (43). Conversely, it was indicated more recently that lower growth rates induced by lower protein synthesis rates are associated with an increased content of ribosomal proteins (39,44).
Protein biosynthesis is a central cellular process and related complexes, especially ribosomes, are conserved in all organisms. Our observation of a warming induced downregulation of the protein biosynthesis machinery across the soil bacterial community led to the following hypothesis: Reduction of cellular ribosome content is a physiological key response of soil microorganisms exposed to warming (Fig. 4A).
To test this hypothesis, we conducted two supplementary experiments. We extracted RNA from soil samples collected from the same A T and E T plots of the geothermal soil temperature gradients in a different season (autumn) four years later (October 2020) and conducted a short-term warming experiment using LTW-A T soil sampled at the same time (Fig. 4BCD). Our rationale for the rst test was that an observation of the same pattern of RNA content relative to soil dry weight, during a different timepoint and season, would be a rst indication that ribosome reduction is a long-term warming response.
The rationale for the experiment was that if this response is also triggered by short-term temperature change, it is directly in uenced by temperature, not just indirectly due to altered nutrient concentrations (Fig. 1B).
The in situ RNA content per unit of soil extracted in autumn 2020 re ected the pattern observed in summer 2016, with a lower RNA content in the warmed soils (Fig. 4C). This indicates ribosome reduction throughout the seasons as a response of microorganisms to medium-and long-term soil warming (i.e. years and decades of warming, respectively). Likewise, we observed a massive reduction to approximately half the RNA content after one week of incubation of ambient soils (A T ) at higher temperature (+6°C), remaining at the lower level throughout the experiment (Fig. 4D, upper bar plots). The decrease of RNA content per unit of microbial C (a proxy for microbial biomass) matched the decrease per gram dry weight soil (Fig. 4D lower bar plots). Thus, cellular ribosome reduction appears to be a fast but stable response to soil warming. We think that higher reaction rates, including higher protein synthesis rates per ribosome (43), caused by warming might allow microorganisms to reduce their cellular ribosome content and increase serial protein synthesis. A resulting ratio change between ribosomal and metabolic proteins was supported by the observed transcriptional patterns, although patterns for most central (carbohydrate) metabolisms were not signi cant. Such re-allocation of energy and matter could potentially accelerate overall metabolism and growth. This presumption is in line with theoretical and experimental work pointing at the importance of covariation of microbial protein composition and growth rate and suggesting that a careful tuning of highly abundant proteins, such as ribosomal proteins, liberates most resources for acceleration of other reactions ((39) and references therein). Very recently, optimal cellular resource allocation was described to be rather constant across non-extreme temperatures, assuming the same temperature-dependency of metabolic and translational rates, but resource allocation and ribosome content will vary with the nutrient status (predictive modelling approach calibrated with E. coil quantitative proteome data) (45). As such, cellular ribosome reduction in medium-and long-term warmed ForHot soils could be mainly caused by the lower C, N, and P content in the warmed soils. We did not observe a reduction in total C and N content during the short-term warming experiment (Table S24), indicating that the increased temperature is the main driver of cellular ribosome reduction on the short-term, but the quality of substrate may have changed. However, studying and quantifying substrate availability and accessibility in heterogeneous soil ecosystems and the nutrient status of individual members of complex soil microbial communities is very challenging. Thus, determining the relative effect of environmental variables affected by warming, such as soil temperature and substrate availability (and quality), on metabolic and physiological changes of microbial cells remains. Furthermore, cultivation-independent in situ studies need to be complemented by pure culture studies on ecologically relevant soil microorganisms, to assess how widespread the observed changes are among individual community members and ecologically relevant functional guilds. Additionally, microbial cells reducing their cellular ribosome content might subsequently reduce their cell sizes. It has been shown that E. coli cell sizes decrease slightly at higher temperatures (46). Thus, the carbon turnover effects of changing cellular density, cell composition, and cell size, potential key microbial responses to global warming, should be targeted in future studies.

Conclusion
How warming affects the microbial control of the global C cycle is a key question to better understand soil-climate feedbacks, an answer to which is urgently needed (4,47). We propose that a downregulation of the protein biosynthesis machinery is a central physiological response during acclimation of microorganisms to short-, medium-, and long-term soil warming. This is suggested by metatranscriptomic data indicating regulatory changes in several connected pathways in two different sites and durations of warming, nucleic acid extractions from two different seasons and years, and the results from a short-term warming experiment. The regulation of protein biosynthesis and the protein biosynthesis machinery allows microorganisms to maintain high metabolic activities and cell division rates even after decades of warming.

Materials And Methods
Grassland sites and soil sampling Soil samples were collected from a natural grassland near Hveragerði (64°0′N, 21°11′W), Iceland, in July 2016. The grassland is part of the ForHot experiment (14) and features two sites, each consisting of replicated soil temperature gradients (Fig. 1A) that were formed by natural geothermal activity. One site has experienced warming for >50 y, possibly since before 1708 (14). Geothermal activities may have varied over time, but warming has been stable in the area since at least 1963, and the warming intensities of the temperature gradients have not varied since detailed measurements began (14)(15)(16). The second site recently developed similar temperature gradients after an earthquake in 2008. We collected 16 samples that were later analysed in detail, i.e. four samples of long-term warmed soils (LTW) exposed to elevated temperatures (+6 °C above ambient, LTW-E T ) and four corresponding controls at ambient temperatures (LTW-A T ) as well as four samples of medium-term (8 y) warmed soils (MTW) exposed to +6 °C above ambient temperatures (MTW-E T ) and four corresponding controls (MTW-A T ). The grassland soils are classi ed as Silandic Andosols, and both sites are dominated by Agrostis capillaris with varying undergrowth and vascular pant and moss cover (14). We took soil samples (0-10 cm depth) from ambient (A T ) and +6 °C (E T ) plots of four replicated gradients per site (i.e. 8 gradients in total; 4 x LTW and 4 x MTW) at one time point (Fig. 1A and Table S1). A T and E T plots of each gradient were consecutively sampled. Samples were immediately frozen in liquid N 2 for subsequent nucleic acid extraction and metatranscriptomics after sieving to 2 mm.

Physicochemical soil properties
Fresh aliquots of each soil sample were extracted with either 1 M KCl or 0.5 M NaHCO 3 solution for 30 min at room temperature or fumigated with chloroform for 48 h and subsequently extracted using the same extractants. Various soil C and N compounds and soil P compounds (Table S1) were quanti ed in the KCl and NaHCO 3 extracts, respectively, using standard procedures (48). C MO , N MO , and P MO were calculated as the differences between DOC, TDN, and TDP contents in the fumigated and non-fumigated extracts. Total C and N contents were analysed in dried soil aliquots using an elemental analyser coupled to an isotope ratio mass spectrometer (EA-IRMS; EA1110 coupled via a ConFlo III interface to a DeltaPLUS IRMS, Thermo Fisher Scienti c). Soil pH was determined from sieved soil samples in 0.05 M CaCl 2 solution and gravimetric water content was determined.

Nucleic acid extractions and sequencing
We extracted total nucleic acids from 16 sieved soil samples (Table S4), i.e. four replicates from each of the four sampled soil groups (LTW-A T , LTW-E T , MTW-A T , MTW-E T ), representing two warming intensities (A T , +0 °C, and E T , +6 °C) of >50 y warmed soils (LTW) and 8 y warmed soils (MTW). These samples were selected because 6 °C above ambient is within the predicted and most severe range of (soil) warming over the next 60-100 y (49,50). Each sample was extracted in triplicate (technical replicates), as previously described (51). Brie y, a phosphate buffer, a detergent solution containing CTAB, and TE saturated phenol (pH 8) were added to 0.3 g of soil in a lysis matrix E tube (MP Biomedicals) containing silica beads and shaken vigorously for 30 s (6.5 m s -1 ) in a FastPrep machine (MP Biomedicals). After centrifugation, the aqueous supernatant was retained. This procedure was repeated two more times using fresh buffer, detergent, and phenol. The three supernatants were then pooled, followed by phenol:chloroform:isoamylalcohol (25:24:1) and chloroform:isoamylalcohol (24:1) extraction and precipitation of the nucleic acids using PEG 8000. The nucleic acids were treated with DNase (RQ1, Promega) before the metatranscriptomes were generated. DNA (pre-treated samples) and RNA quantity and quality were assessed by automated agarose gel electrophoresis (Bioanalyzer 2100, Agilent), a NanoDrop TM spectrophotometer (ND-1000, Peqlab), and the PicoGreen TM and RiboGreen TM assay kits (Thermo Fisher Scienti c). Total RNA was ampli ed using the MessageAmp II-Bacteria Kit (Ambion Life Technologies) with an input of 100 ng RNA, following the kit protocol, except the RNA was linearly ampli ed for 14 h. The kit uses E. coli poly(A)-polymerase; it preferentially polyadenylates mRNA and by that facilitates and unbiased ampli cation an enrichment mRNA (52). The three technical replicates of each of the 16 samples were pooled prior to the ampli cation. Overlapping paired-end 125-bp reads were sequenced using the HighSeq2500 sequencer (Illumina; library generation: NEBNext Ultra TM RNA Library Prep Kit for Illumina) at the Vienna Biocenter Core Facilities, Vienna, Austria.
Sequence data pre-processing We performed the following pre-processing steps and the majority of the subsequent analyses using the Life Science Compute Cluster (LiSC) run by CUBE (Division of Computational Systems Biology), Department of Microbiology and Ecosystem Science, University of Vienna, Austria. PEAR v.0.9.10 (Paired-End reAd mergeR, default settings) was used to merge the raw paired-end reads (53). We subsequently used SortMeRNA v.2.1 to lter out non-rRNA reads from the total RNA reads (54). The non-rRNA reads were quality ltered (-min_len 180 -max_len 240 -min_qual_mean 30 -ns_max_n 5 -trim_tail_right 15 -trim_tail_left 15) using prinseq-lite v.9.20.4 (55). A second ltering step was performed to obtain mRNA reads and reduce the size of the dataset for later analyses. All non-rRNA reads that aligned to any sequence in the NCBI nr (56) database (as of September 2018) using a non-conservative DIAMOND blastx (57) search (v.0.9.18, -k 1 -e 0.001) were kept. See Table S4 for more details.

Functional and taxonomic annotation
We functionally annotated the total mRNA reads (Table S4) by aligning them to the KEGG database (21) (as of April 2019) using a DIAMOND blastx (57) search (v.0.9.18, -k 1 --min-score 52.5). A prior analysis indicated that at a bit score of 52.5 equalled an e-value of ≤0.0001 (DIAMOND blastx of the same dataset against the KEGG database (21) as of February 2019). We, therefore, used a bit score of 52.5 rather than an e-value of 0.0001 for all following analyses to obtain a better comparability between different database searches, independent of increasing database sizes. We performed a DIAMOND blastx (57) search (v.0.9.25, -k 25 --min-score 52.5) against the NCBI nr (56) database (as of October 2019) to obtain taxonomic information for the mRNA reads. We used the program blast2lca, the standalone implementation of the MEGAN (v.6.11.1) LCA (lowest common ancestor)-assignment algorithm, to obtain one taxonomic assignment for each read based on the 25 nr database hits (58). The following parameters were used to obtain the LCAs: -ms 50 -me 0.01 -top 5 -mid 0.0 (mapping le: prot_acc2tax-Jul2019X1.abin).
Filtering and normalisation. First, we normalised the data to library size and transferred them to permille (‰). Second, only functions (KOs, i.e. KEGG database entries, i.e. orthologues genes) present in all four replicates of at least one group (LTW-A T , LTW-E T , MTW-A T , or MTW-E T ) were kept (removing 0.6‰ of the data). Third, low abundant functions were removed by ltering out KOs with a total relative abundance (sum of all 16 datasets) <0.1‰ (removing 4.6‰ of the data remaining after the previous step), which equalled a mean relative abundance of a speci c function per sample of <0.00625‰ (mean relative abundance [‰] = 0.1‰/16). A second dataset containing all taxonomically classi ed mRNA reads (without functional annotation) was likewise normalised and transferred. The taxonomy-only dataset (tax.data) was ltered by removing all families (taxonomic strings from domain to family) that were not present in all replicates of at least one group (LTW-A T , LTW-E T , MTW-A T , or MTW-E T ) (removing 0.07‰ of the data); no second lter was applied. Viral reads were summarised prior to ltering depending on the subsequent analyses.
Analysis of contextual data. Total, dissolved, and microbial C, N and, P concentrations, nucleic acid concentrations, soil pH, and water content were individually tested for signi cant differences between de ned groups (groups = A T and E T , MTW and LTW, MTW-A T , MTW-E T , LTW-A T , and LTW-E T ) using twoway ANOVAs (basic R function aov) and computing Tukey Honest Signi cant Differences (basic R function TukeyHSD) to evaluate sample selection and grouping. X-fold changes were calculated in R: (mean E T /mean A T ).
Correlations. The basic R function cor.test was used to identify associations between microbial biomass (using microbial C as proxy) and nucleic acid contents by applying Spearman's rank correlation (twosided; P values were corrected for multiple testing using Benjamini-Hochberg procedure, basic R function p.adjust). The ggplot function geom_smooth was used to indicate correlations (method = lm).
Richness estimates. Taxonomic (family) and functional (KO) richness was estimated from raw read counts (families and KOs not present in all four replicates of at least one group (LTW-A T , LTW-E T , MTW-A T , or MTW-E T ) were considered as noise and excluded) using the vegan function rarefy. The basic R function t.test was used to perform two-sided Student's t-Tests to identify signi cant differences between A T and E T (of LTW and MTW, respectively, or across LTW and MTW). Obtained P values were corrected for multiple testing (Benjamini-Hochberg procedure, basic R function p.adjust).
Non-metric multidimensional scaling (NMDS) and permutational multivariate analysis of variance (PERMANOVA). NMDS was used to obtain ordination plots depicting (dis)similarities between expressed microbial functions (KOs) and microbial community structures of the samples (Fig. S2). We used the metaMDS function implemented in the R package vegan, two dimensions, and a maximum of 10,000 random starts in search of a stable solution. The sequencing data were normalised and ltered as described above prior to the NMDS analyses. GUSTA ME (GUide to STatistical Analysis in Microbial Ecology) (60) was consulted for selecting the appropriate dissimilarity index (i.e. "canberra"). PERMANOVA (vegan function adonis) was used to identify the effect of warming (A T and E T ) and warming duration (LTW and MTW), respectively, on the distribution of samples by gene expression (9,999 permutations, dissimilarity index "canberra").
Alluvial plots. Alluvial plots (Sankey diagrams) were created using the R package ggalluvial. Individual Sankey diagrams were manually merged if more than two levels were shown.
Heatmaps. Heatmaps were generated using the geom_tile function of the ggplot2 R package. Before plotting, normalised data were transformed by z-scoring, either over all 16 samples or separately for LTW and MTW. Two explorative abundance-pattern lters were subsequently applied for selecting patterns of interest. Abundance-pattern lter applied on KEGG3 categories (Fig. S6, Fig. S7): The relative abundances of KEGG3 categories between A T and E T across all samples were compared. Patterns were retained only if: i) at least six samples of one temperature group (A T or E T ) were higher than the third highest sample of the other temperature group (E T or A T ) and at least seven samples of one temperature group (E T or A T ) were lower than the third lowest sample of the other temperature group (A T or E T ) or ii) at least seven samples of one temperature group (A T or E T ) were higher than the third highest sample of the other temperature group (E T or A T ) and at least six samples of one temperature group (E T or A T ) were lower than the third lowest sample of the other temperature group (A T or E T ). Therefore, the critical threshold to pass the abundance pattern lter lies at 80% consensus with the most stringent warming-associated distribution (i.e. the eight highest relative transcript abundances are found in one temperature group and the eight lowest in the opposite temperature group). Abundance-pattern lter applied on individual subsets (MTW and LTW): A more stringent lter (referred to as 4/4-lter) was used to analyse the MTW and LTW individually; (Fig. S9) only taxa with higher or lower relative abundances in all four replicates of one group relative to their counterparts passed the lter threshold (e.g. a taxon passes the lter if the four highest values are found in LTW-A T and the four lowest in LTW-E T ). Log2-fold changes in mRNA abundances of selected taxonomic groups were calculated in R: log2(mean E T /mean A T ).
Differential gene expression (DGE). To identify differences in gene expression (mRNA abundances) between A T and E T (of LTW and MTW, respectively, or across LTW and MTW) DGE analyses were performed using edgeR (function glmQLFTest), recommended for a small number of replicates (<12) (61). Raw read counts were ltered as described in the "Filtering and normalisation" section but normalized using the default trimmed mean of M-values normalization (TMM) method implemented in edgeR. DGE analysis on functionally annotated mRNA reads (KO level, i.e. orthologous genes) was performed as follows: Eukaryotic mRNA was excluded from the analysis (to eliminate potential biases related to taxonomic composition) and a minimum of 500 physical read counts per sample was used as a threshold for a KO to be considered in the DGE analysis (to exclude infrequent transcripts and allow a taxon-speci c breakdown). KOs with a false discovery rate (FDR) between 0.05 and 0.1 were considering as transcripts of interest when being part of pathways or complexes containing signi cantly differentially expressed transcripts (FDR < 0.05). DGE analysis on taxonomically annotated mRNA reads (domain/kingdom, phylum, class, order, and family level): DGE analyses were performed on taxonomically annotated datasets (func.tax.data and tax.data). The threshold for a taxon to be considered in the DGE analysis varied depending on the taxonomic rank and ranged from 5-500 physical mRNA reads.
Supplementary experiments RNA extractions. The same A T and E T plots that were analyzed in detail in our metatranscriptomics study were re-sampled in October 2020, four years after the metatranscriptomics sampling campaign (Fig. 4B, Table S22). After removing the plant cover, soil from the upper 5 cm was sampled and kept on dry ice or at -80 °C until processing. The samples were homogenized and grinded on liquid N. Afterwards total nucleic acids were extracted as described above. Total RNA content was quanti ed using Qubit RNA HS Assay Kit (Thermo Fisher Scienti c). The basic R function t.test was used to perform a one-sided Student's t-Test to test the hypothesis that warming leads to a reduction in the RNA content.
Short-term-warming experiment. A short-term warming experiment (Fig. 4D) was performed using LTW-A T soil sampled in October 2020. The sample (~ 270 g) was shipped (on cool packs), sieved (2 mm), homogenized, and pre-incubated (3 weeks) in a glass beaker (in the dark, tightly closed with tin foil) at the mean ambient October soil temperature (7 °C). After pre-incubation the homogenized soil was sampled (t0 sampling; starting point of the experiment; n=5). Subsequently, the remaining soil was distributed to 10 serum bottles (closed with a rubber stopper). The regularly aerated bottles were incubated for 6 weeks at 7 °C (control, n=5) and 13 °C (warming treatment, +6 °C above ambient, n=5), respectively, and sampled after one, three and six weeks of incubation. Total nucleic acids were extracted using the phenol:chloroform triple bead-beating protocol described above. Total RNA content was quanti ed using Qubit RNA HS Assay Kit (Thermo Fisher Scienti c). Microbial C (C MO ) and total soil C and N were determined as described above. The basic R function t.test was used to perform two-sided Student's t-Tests to identify signi cant differences between control and warming treatment at the different sampling timepoints.
Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The raw sequence data are available at the NCBI Sequence Read Archive (SRA); BioProject ID: PRJNA663238, accession numbers SAMN16124403-SAMN16124422 (release date: 1 st of October 2021). Other underlying data (if not available as Supplementary Table) and scripts are available upon request (and will be published on GitHub when the paper is published). with signi cant differences in the respectively aligned physicochemical and biological properties (twoway ANOVAs on each individual property, TukeyHSD test to correct for multiple testing, see Table S2 for the full statistical analysis and exact P values). (C) Nucleic acid (DNA and RNA) content per unit of microbial C (boxplots) and correlation between microbial C and RNA content per g dry weight soil (scatterplot). See Table S2 and S3 for the full statistical analysis and exact P values (*, Padj < 0.05).  Table  S16 for relative abundances of all individual KOs summarized in the boxplots and details on the DGE analysis.

Figure 3
Enzymes and enzyme complexes involved in central carbohydrate metabolisms and cell replication and taxonomic distribution of overall expression patterns. (A) A metabolic pathway map focusing on central (carbohydrate) metabolisms, associated reactions, and relative mRNA read abundances assigned to these reactions. Reactions are represented by their EC numbers (EC-no.) and derived from the KEGG pathway drawings. Color code indicates mean relative expression levels (red, higher relative expression levels in ET; blue, higher relative expression levels in AT); intensities additionally re ect the presence of differentially expressed KOs and KOs with an FDR between 0.05 and 0.1. Bold fonts indicate Metabolism sub-categories (KEGG3); dashed arrows indicate connections to sub-categories for which no details are given; the color code (same as above) of the dashed frames around these categories indicate if the whole sub-category shows indications of a temperature-dependent response. See Table S17 for