Central metabolic responses of microorganisms to years and 1 decades of soil warming

Microbial physiological responses to long-term warming are poorly understood. Here we applied 17 metatranscriptomics to investigate how microorganisms react to medium-term (8 years) and long- 18 term (>5 decades) subarctic grassland soil warming of +6 °C. 19 Decades, but not years, of warming induced changes in relative abundances of eukaryotic, 20 prokaryotic, and viral transcripts and reduced functional richness. However, irrespective of the 21 duration of warming, we observed a community-wide upregulation of central (carbon) metabolisms 22 and cell replication in the warmed soils, whereas essential energy metabolism and protein 23 biosynthesis complexes and pathways were downregulated. This coincided with a decrease of 24 microbial biomass and lower soil substrate concentrations (e.g. dissolved organic carbon and 25 phosphorus). 26 We conclude that permanently accelerated reaction rates at higher temperatures facilitate a 27 downregulation of energy metabolism and protein biosynthesis, potentially freeing energy and 28 matter for substrate acquisition and growth. This resource allocation seems to be a common 29 response in microorganisms and allows sustaining high metabolic activities and replication rates even 30 after decades of soil warming. the soil-climate feedbacks, in comprehensive study of the transcriptional response of microorganisms across all domains life, we show that downregulation of energy metabolism and protein biosynthesis is central in the microbial soil warming response that allow microorganisms to maintain high metabolic activities and cell division rates even after decades of warming. We suggest that accelerated biochemical reaction rates due to higher temperatures, have a positive feedback on central metabolisms via increased relative transcript abundances, if not constraint by reduced substrate concentrations, further accelerating microbial decomposition of SOM and the subsequent release of CO 2 to the atmosphere.


33
Global temperatures and atmospheric carbon dioxide (CO2) levels have increased steadily over the 34 last 100 years 1,2 . The terrestrial carbon (C) cycle feedback to the climate system represents a major 35 uncertainty in the prediction of future temperatures 3 . Soil microorganisms, including Bacteria, 36 Archaea, Fungi, and protists, are responsible for the turnover of soil organic matter (SOM) and the 37 subsequent release of CO2 from soils to the atmosphere 2 . Higher temperatures commonly lead to 38 higher microbial activities, so global warming should accelerate the decomposition of SOM to CO2 4 . 39 On the other hand, SOM consists largely of microbial necromass and warming may stimulate 40 microbial growth and thus necromass production 5,6 , promoting SOM formation. Whether soils will 41 ultimately act as C sinks or sources, thus depend on how microorganisms respond to long-term 42 warming. Nevertheless, microbial responses to global warming are currently poorly represented in 43 Earth system models 7 , which can, to some extent, be attributed to the challenges associated with 44 studying and quantifying microbial responses to environmental change in complex soil 45  Table 3). 111 112 DNA and RNA concentrations per unit of microbial C (CMO) had opposite trends, with higher DNA but 113 lower RNA concentrations in the warmed soils, irrespective of the warming duration (Fig. 1c). While 114 the RNA content per g soil positively correlated with CMO (rs(16) = 0.81, Pcorr < 0.001), the DNA content 115 did not (rs(16) = 0.12, Pcorr = 0.664) (Fig. 1c). CMO was likewise positively correlated with total soil C, N, 116 and P, dissolved organic C and P, microbial N and P, and water content (Supplementary Table 3).  117 To investigate what structural and functional changes in the soil microbiomes were associated with 118 these substantial differences in microbial biomass and soil chemistry we performed 119 metatranscriptomics. 120

Compositional changes of the functional microbial communities across all domains of life 121
Illumina paired-end sequencing produced an average of 6.69  10 6 mRNA reads per sample 122 (Supplementary Table 4). Bacteria dominated the mRNA transcript pools (93.36-98.52%), followed 123 by Eukaryota (1.28-6.45%), Archaea (0.5-1.6‰), and viruses (0.2-1.9‰). Overall, more than 1,000 7 different families were detected, and mean family richness was lower in the heated soils, albeit not 125 significant ( Supplementary Fig. 1a). An NMDS analysis of mRNA reads assigned at family level 126 indicated no clear overall separation between AT and ET samples and a greater variability between 127 MTW replicates than LTW replicates ( Supplementary Fig. 2a). Thus, we analysed the detailed 128 taxonomic assignments of LTW and MTW separately (Fig. 2). In the long-term warmed soils, but not 129 in the medium-term warmed soils, several taxa showed warming-induced differences in relative 130 transcript abundances (Fig. 2a, Supplementary Table 6). Those that may be affected by warming 131 include phylogenetic groups across all domains of life and included Fungi, protists, and Chloroflexi, 132 which showed lower mean relative abundances in LTW-ET, and Deltaproteobacteria, Planctomycetes, 133 Verrucomicrobia, as well as viruses, which showed higher mean relative abundances in LTW-ET 134 compared to LTW-AT (Fig. 2a, b). A taxonomic analysis of rRNA reads corroborated the bacterial 135 mRNA profiles, but the mean relative abundances of Eukaryotes differed between in the rRNA 136 datasets and the mRNA datasets ( Supplementary Fig. 3). 137 We further screened all taxonomic ranks down to family level (Fig. 2c). Five taxa showed lower 138 relative expression levels in LTW-ET than LTW-AT, and 17 taxa showed higher relative expression 139 levels in LTW-ET than LTW-AT (Fig. 2c). Differential gene expression analyses underpinned these 140 results, as many of the reported taxa were significantly different. In contrast, none of the six taxa 141 that displayed higher relative expression levels in MTW-ET than MTW-AT was significantly differential 142 expressed (Fig. 2c). We repeated the analysis on the bacterial fraction of the mRNAs to test whether 143 Eukaryota abundances influenced the bacterial patterns and found almost identical warming-induced 144 abundance patterns ( Supplementary Fig. 4). We also attempted to increase the taxonomic resolution 145 by assembling the metatranscriptome reads (Supplementary Methods). Our approach using 146 rnaSpades 22 , however, did not provide a sufficient number of long mRNA contigs (<10% of 147 functionally annotated mRNA contigs were longer than the unassembled reads).  Table 5). 174 To identify the nature of this functional response we explored in more detail the functional 175 assignments to KEGG1 categories Metabolism (65.4%), Genetic information processing (16.8%), 176 Cellular processes (9.0%), and Environmental information processing (8.9%) (Fig. 3). Transcripts for 177 Metabolism and major Metabolism subcategories such as Carbohydrate metabolism, Lipid 178 metabolism, and Nucleotide metabolism had higher relative abundances in the warmed soils, while 179 Genetic information processing subcategories such as Transcription and Translation had lower 180 relative abundances (Fig. 3, Supplementary Table 8). Transcript patterns for Cellular processes and 181 Environmental information processing were not consistent.  Table 8).

Community-wide shifts in gene expression of central metabolisms and cellular functions 191
Since the gene expression patterns observed within the broad functional categories suggested 192 functional temperature-dependencies ( Fig. 3) we screened all KEGG3 categories to identify which 193 specific metabolic pathways and functional complexes were the basis for these overall patterns. 194 One third of all functionally annotated mRNAs were affiliated with KEGG3 categories that exhibited 195 warming-induced abundance patterns similar to those in Fig. 3 (Supplementary Fig. 5). Out of 55 196 KEGG3 categories with a visually distinct pattern, 13 showed lower relative expression levels in ET 197 than AT, while 42 KEGG3 categories showed higher relative expression levels in ET than AT. We further 198 investigated the most abundant of these 55 KEGG3 categories, each making up >1% of total 199 annotated mRNA reads. Eight of these abundant categories showed higher relative expression levels 200 11 in ET and encompassed central C metabolisms and metabolic pathways associated with amino acids 201 and nucleotides, while five showed higher relative expression levels in AT and included protein 202 biosynthesis related complexes and Oxidative phosphorylation (i.e. ATP formation) (Fig. 4a) . This 203 change in the transcript abundances for the most central metabolisms in organismal function was 204 observed for multiple taxa (Fig. 4b). In the long-term warmed soils, four phyla, eight classes, four 205 orders, and four families, representing 33.9%, 38.1%, 11.1%, and 1.8% of the LTW 206 metatranscriptomes, respectively, displayed this pattern, confirming that this is a taxonomically 207 broad type of warming response. This was supported by the analysis of the medium-term warmed 208 soils, where seven phyla, 14 classes, seven orders, and eleven families, representing 87.2%, 65.0%, 209 24.1%, and 6.0% of the MTW metatranscriptomes, respectively, displayed the same expression 210 pattern. It should be noted that the low percentages at family level is due to an average of only 211 22.1% of the annotated mRNA reads being assigned to a family-level taxon. Interestingly, the pattern 212 occurred in taxa ( Fig. 4b taxa in bold) regardless of whether these taxa became more active (higher 213 relative abundance of mRNA and rRNA) or less active with warming ( Fig. 2, Supplementary Fig. 3), 214 indicating that the physiology of these community members had been altered in a similar way. 215 Furthermore, our extensive screening revealed that on average this was true for nearly two thirds 216 (63.7%) of the taxa present in LTW and MTW (Fig. 4c, Supplementary Fig. 6). However, the 217 percentage of taxa expressing the pattern varied between 43 and 95% depending on the KEGG3 218 category and the warming duration and was highest for the KEGG3 categories related to RNA, 219 protein and energy metabolisms.  Table 9).

255
We observed trends of warming-induced transcriptional responses in a broad set of genes 256 fundamental to bacterial cell replication (Fig. 5, Supplementary Table 9). In LTW, we observed 257 14 increased relative transcript abundances for Peptidoglycan biosynthesis, DNA replication and FtsZ in 258 ET, albeit insignificant. However, in MTW-ET, relative transcript abundances for Peptidoglycan 259 biosynthesis, Glycerolipid metabolism, and DNA replication were significantly higher in MTW-Et than 260 MTW-AT (t-tests, n = 8, Pcorr < 0.05, Fig. 5), while the trend for FtsZ transcripts was similar but 261 insignificant. 262 While the relative abundances of cell replication-related transcripts increased with warming, 263 transcripts for energy conservation (Oxidative phosphorylation) showed the opposite trend. All four 264 enzyme complexes of the membrane-bound prokaryotic electron transport chain and ATP synthase 265 (producing ATP, the energy currency of cells) showed lower mean relative transcript abundances in ET 266 (Fig. 5a). In contrast, transcripts of enzymes providing inorganic phosphate (Pi) to the ATP synthase, 267 especially polyphosphate kinase, showed higher mean relative abundances in ET than AT. 268 Similar to Oxidative phosphorylation (i.e. ATP formation), transcripts for multiple complexes related 269 to protein biosynthesis showed a trend towards lower relative abundances in ET (Fig. 6b) were also observed for RNA polymerase, ribosomal proteins, ribosomal RNA, the main molecular 272 chaperones GroEL and DnaK (involved in protein folding), and the Sec-dependent pathway, the 273 dominant protein export pathway present in the metatranscriptomes. Contrary to all other 274 categories of the protein biosynthesis machinery, transcripts for RNA degrading complexes (RNA 275 degradosomes) showed higher relative abundances in ET than AT (Fig. 6b). 276 We also investigated in detail the gene expression profiles for central C metabolisms. While the 277 transcript abundances of KEGG3 categories Pyruvate metabolism, Glycolysis/Gluconeogenesis, 278 Methane metabolism, and C fixation pathways in prokaryotes pointed towards higher expression 279 levels in the warmed soils (Fig. 4), these differences were not significant. This suggests that while a 280 shift in community metabolism is indicated, this might not be reflected consistently in all pathway 281 steps and responsible bacteria that contribute to the above mentioned KEGG3 categories.  warming-induced abundance pattern towards lower transcript abundances in ET than AT as seen in (Fig. 4a).

293
Finally, we examined, within all KEGG3 categories, the differential expression of single functions (KO). 294 Notably, the observed patterns, including those not significant at the KEGG3 level, were confirmed by 295 multiple functions (KO) that were differentially expressed and significant after multiple testing 296 16 correction (Supplementary data S2). For functions related to protein biosynthesis and oxidative 297 phosphorylation (ATP formation) the patterns were particularly strong and supportive of the above 298 described trends. Less, although substantial, support was found for central carbon and amino acid 299 metabolism patterns. 300 301

302
Here we used metatranscriptomics to elucidate how soil microorganisms change their gene 303 expression of central metabolic functions and cellular processes in response to warming. We showed 304 that long-term (>50 y), but not medium-term (8 y), soil warming (+6 °C) resulted in a significant, 305 albeit small, decrease of unique molecular functions encoded by soil Bacteria, Archaea, Fungi, and 306 microbial Eukaryotes. However, irrespectively of the duration of warming, physiological responses to 307 warming as shown by the transcriptional tuning of central metabolisms for C, energy, protein 308 biosynthesis, and growth, were consistent and common across many community members. 309 The apparent overexpression of genes involved in Pyruvate metabolism, Glycolysis/Gluconeogenesis, 310 Methane metabolism, and C fixation pathways in prokaryotes indicate that long-term as well as 311 medium-term warming affect the central cellular C metabolism of soil microorganisms. However, 312 inconsistencies in the relative transcript abundances of subjacent KOs indicated that the perceived 313 upregulation of these metabolisms does not involve all associated enzymes, reflecting either how 314 metabolic pathway efficiency can be controlled by only a few rate-limiting steps 24 or that the 315 patterns created by the distinct up-regulation of pathways in some organisms are diluted by other 316 organisms that share pathway steps but not the response pattern. 317 In contrast, the reduction in relative transcript abundances of oxidative phosphorylation complexes, 318 RNA polymerases, protein processing enzymes, and ribosomal proteins were consistent across 319 multiple enzyme subunits and pathway steps. A downregulation of the protein biosynthesis 320 machinery is not immediately recognizable as a strategy to counteract substrate limitation or 321 maximize growth rates. However, a reduction of ribosomes may be biochemically favourable under 322 the warming conditions. Ribosomes are macromolecular complexes that are the sites of protein 323 synthesis (translation). They consist of ribosomal proteins and rRNAs and are present in tens of 324 thousands of copies per bacterial cell 25 . Thus, ribosomes can comprise over one third of the dry cell 325 mass 26 and rRNAs can account for >90% of the total cellular RNA content 27 . Starving Escherichia coli 326 and Salmonella spp. cells reduce their ribosomal content 28,29 , suggesting that a downregulation of the 327 translation machinery, which accounts for up to 40% of total cellular proteins 30 , is metabolically 328 favorable in nutrient-limited ecosystems 31 . Furthermore, carefully tuning concentrations of abundant 329 proteins such as ribosomal proteins can also free resources for accelerating other reactions 30 . We not 330 only demonstrated lower relative transcript abundances of ribosomal proteins and protein 331 biosynthesis-related enzymes in the warmed soils, but our results also showed lower microbial 332 biomass per gram of soil dry weight, and a decrease in RNA per gram dry soil that correlated 333 significantly with the biomass. A significant decrease in microbial biomass per gram of soil dry weight 334 in the warmed soils has been observed previously at the same site 16 . Furthermore, the transcriptional 335 patterns for protein biosynthesis were consistent for many community members, regardless of their 336 relative abundances of total mRNAs and rRNAs in the long-and medium-term warmed soils and the 337 controls. Thus, a downregulation of protein biosynthesis-related enzymes and a reduction of 338 ribosomes might be a common response of microorganisms, especially bacteria, to ecosystem 339 warming and partly responsible for the lower microbial biomass in the warmed soils. It seems 340 contradictory that we observed, at the same time, higher relative transcript abundances for enzymes 341 related to cell division and biosynthesis of amino acid, the building-blocks of proteins. However, a 342 reduced ribosomal content is not necessarily linked to lower protein synthesis rates or lower growth 343 rates. It is long known that increased temperatures accelerate mRNA synthesis and the protein 344 synthesis rate per ribosome (peptide chain elongation rate) in the model organism Escherichia coli 32 . 345 Conversely, it was indicated more recently that lower growth rates induced by slower reaction rates 346 are associated with an increased content of ribosomal proteins 30,33 . Likewise, an acceleration of 347 reaction rates in the membrane-bound electron transport chain and ATP synthesis (Oxidative 348 phosphorylation) despite lower expression might also have been driven by temperature and might be 349 to metabolically more active microbial cells driven by temperature and a different partitioning of 359 cellular resources. Previous observations of higher biomass-specific growth, respiration, organic C 360 uptake, and turnover rates in the same long-term warmed soils 16 corroborates this interpretation. 361 These results are clearly not consistent with the suggested thermal acclimation of soil 362 microorganisms used to explain a return to pre-warming SOM degradation rates after a few years of 363 warming 11,12 (and references therein). Rather, the apparent thermal acclimation in soil respiration 364 rates might be driven by a reduction in microbial biomass caused by reduced carbon and nutrient 365 concentrations, as previously suggested 16 . However, our results point at a different type of microbial 366 acclimation to warming where the physiological adjustments allow the microbial community 367 members to maintain a high activity even after decades of warming, despite more limiting 368 conditions. 369 The microbial functional responses were more pronounced and widespread in the microbial 370 communities in the medium-term warmed soils than in the long-term warmed soils reflecting the 371 19 recently reported systemic overreaction to years versus decades of warming in the same soils 15 . The 372 authors proposed that an initial acceleration of biotic activity due to warming led to rapid decreases 373 in C, N, and P pools within the first years after the onset of warming, followed by a decrease of 374 microbial and fungal biomass and system stabilisation within decades. Our taxonomic analyses of the 375 mRNA transcript pools extended our insight into this proposed soil warming chronology; long-term, 376 but not medium-term, soil warming resulted in differential abundances. Relative mRNA transcript 377 abundances of several bacterial taxa, fungi, protists, and viruses were affected by long-term 378 warming. This indicated shifts in trophic interactions (e.g. from top-down to bottom-up control of the 379 prokaryotic community) and possibly a reduced importance of fungi, with repercussion on organic 380 matter decomposition. Reduced Fungi:Bacteria ratios have been described previously in long-and 381 short-term warming experiments of forest and grassland soils 8,37 and Fungi appear to be more 382 abundant and active at lower temperatures 4,38 and in soils with low pH 39,40 . Besides the higher 383 temperatures and slightly higher pH in the warmed ForHot grassland soils, the lower concentrations 384 of various soil C, N and P compounds and decreased soil aggregate sizes 19 may explain the lower 385 relative abundances and diversity of fungal mRNAs. The few and non-significant taxonomic 386 differences between ambient and medium-term warmed soils suggest that while the soil 387 microorganisms respond functionally, there is little effect on microbial community structure. 388 However, we cannot exclude that the apparent lack of taxonomic response after eight years resulted 389 from the high variation between the biological replicates. Our results somewhat contrasts previous 390 studies on the ForHot grassland sites that reported generally little taxonomic shifts with warming 16,20 . 391 However, these studies applied DNA-based methods (i.e., 16S gene and ITS amplicon sequencing) 392 targeting the potential microbial community, including active cells, dormant cells, and relic DNA. 393 These methods are also prone to biases from PCR and primers. In contrast, metatranscriptomics is 394 PCR-free and deploys random hexamer-primers targeting prokaryotic, eukaryotic, and viral 395 transcripts in parallel, thus allowing a comprehensive analysis of the active soil microbiome. 396 Furthermore, our analyses did demonstrate a stronger effect of warming on microbial functions 397 20 (KOs) than on taxonomy. This functional response was similar in medium-term and long-term 398 warmed soils. Thus, our study provides evidence for a decoupling of microbial community structure 399 and functions, as recently suggested from several ecosystems including soil [41][42][43][44] . 400 How warming affects the microbial control of the global C cycle is a key question to better 401 understand the soil-climate feedbacks, an answer to which is urgently needed 4,45   Heatmaps. Heatmaps were generated using the geom_tile function of the ggplot2 R package. Before 515 plotting, normalised data were transformed by z-scoring, either over all 16 samples or separately for 516 TLW and MTW. Two explorative filters were subsequently applied for selecting patterns of interest. 517 Abundance patterns of taxa present in LTW and MTW, respectively, were subject to a stringent filter 518 (referred to as 4/4-filter); only taxa with higher or lower relative abundances in all four replicates of 519 one group relative to their counterparts passed the filter threshold (e.g. a taxon passes the filter if 520 the four highest values are found in LTW-AT and the four lowest in LTW-ET). A less stringent filter 521 (referred to as 13/16-filter) was applied when the relative abundances of KEGG functions between AT 522 and ET across all samples were compared (e.g. Figure 4). Patterns were retained only if: i) at least six 523 samples of one temperature group (AT or ET) were higher than the third highest sample of the other 524 temperature group (ET or AT) and at least seven samples of one temperature group (ET or AT) were 525 lower than the third lowest sample of the other temperature group (AT or ET) or ii) at least seven 526 samples of one temperature group (AT or ET) were higher than the third highest sample of the other 527 temperature group (ET or AT) and at least six samples of one temperature group (ET or AT) were lower 528 than the third lowest sample of the other temperature group (AT or ET). Therefore, the critical 529 threshold to pass the 13/16-filter lies at 80% consensus with the most stringent warming-associated 530 distribution (i.e. the eight highest relative transcript abundances are found in one temperature group 531 and the eight lowest in the opposite temperature group). 532 Alluvial plots. Alluvial plots (Sankey diagrams) were created using the R package ggalluvial. Individual 533 Sankey diagrams were manually merged if more than two levels were shown. 534 Boxplots. Boxplots were generated using geom_boxplot (R package ggplot2). 535 Non-metric multidimensional scaling (NMDS). NMDS was used to obtain ordination plots depicting 536 (dis)similarities between the microbial functions and microbial community structures of the samples. 537 We used the metaMDS function implemented in the R package vegan, two dimensions, and a 538 26 maximum of 10,000 random starts in search of a stable solution. The sequencing data were 539 normalised and filtered as described above prior to the NMDS analyses. GUSTA ME (GUide to 540 STatistical Analysis in Microbial Ecology) 57 was consulted for selecting the appropriate dissimilarity 541 index (i.e. "canberra"). 542 Statistics and post hoc analyses. The basic R function t.test was used to perform two-sided Student's 543 t-Tests to identify significant differences between AT and ET (of LTW and MTW, respectively, or across 544 LTW and MTW). Obtained p-values were corrected for multiple testing (Benjamini-Hochberg 545 procedure, basic R function p.adjust). Corrected p-values (Pcorr, q-values) < 0.05 were considered to 546 indicate significant differences. We deliberately chose a parametric test (decreasing the chance of 547 making a type II error) combined with multiple-testing adjustment and considered corrected p-values 548 < 0.1 to indicate a temperature-dependent trend, reflecting the explorative nature 58 of our study. 549 The basic R function cor.test was used to identify associations between microbial biomass and RNA 550 content, DNA content, and various C, N, and P concentrations by applying Spearman's rank 551 correlation (two-sided). The ggplot function geom_smooth was used to indicate correlations (method 552 = lm). Taxonomic