Terrestrial ecosystems play a vital role in the global carbon (C) cycle as soils are estimated to store twice as much C as the atmosphere and all vegetation combined (Schimel, 1995; Scharlemann et al., 2014). Soil C-dynamics is however tightly coupled to temperature changes making C sinks vulnerable to global warming. Studies have shown that warming-induced increases in soil respiration have the potential to convert soil from a C-sink to a C-source triggering a self-reinforcing C cycle-temperature feedback (Friedlingstein et al., 2006; Allison and Treseder, 2011; Tucker et al., 2013). Microbes play a pivotal role in soil C processes, acting to process incoming C to form stable soil organic matter (SOM) pools (Baldrian et al., 2012; Lladó et al., 2017) as well as acting to destabilize stored soil organic carbon (Bradford et al., 2008; Feng et al., 2008; Pisani et al., 2015). Increasing temperature as a result of climate change has significantly affected the microbial communities and microbial functions (Castro et al., 2010; Gray et al., 2011; Sheik et al., 2011; Zhou et al., 2012). However, we do not yet understand the microbial mechanisms driving warming-induced soil C losses well enough to predict these dynamics in a warming climate.

One challenge in identifying the mechanisms associated with microbial-driven soil C loss is the inherent variability in short- vs. long-term responses to warming. There is growing evidence that with long-term warming, there is a gradual decrease in warming effects on soil respiration (Luo et al., 2001, Melillo et al., 2002, 2017; Kirschbaum, 2004; Romero-Olivares et al., 2017). For example, in an ongoing soil warming experiment at the Harvard Forest in central Massachusetts, soils that are warmed 5°C above ambient temperatures fluctuated from being initially higher (Phase I; 1–10 years), then invariant (Phase II; 11–17 years), then higher (Phase III; 18–23 years) and currently approximately equal to (Phase IV) in respiration compared to control plots (Melillo et al., 2002, 2017). Such time-dependent variations can be due to acclimation of soil microorganisms to the increased temperature changes (Bradford et al., 2019; Dacal et al., 2019) and/or functional shifts in microbial community associated with decomposition of labile and recalcitrant C-pools based on their relative availability (Hartley et al., 2009; Schindlbacher et al., 2015; Pold et al., 2016).

Our recent work in Harvard Forest suggests that chronic warming has led bacteria to play an increasing role in SOM turnover (Frey et al., 2008; DeAngelis et al., 2015). Around 26 years of simulated warming was associated with sustained reduction in fungal biomass but not bacterial biomass (Frey et al., 2008; DeAngelis et al., 2015). We observed warming-accelerated respiration despite declining quality and quantity of SOM (Bradford et al., 2008; Pold et al., 2015, 2017). Periods of soil carbon decay were punctuated by periods of change in the microbial communities (Melillo et al., 2017), though the bacterial community composition did not appear to have changed substantially. Because two-thirds of soil respiration is induced by microbial activities (Melillo et al., 2002, 2011), we predicted that changes in microbial gene expression, as measured by meta-transcriptomics, could potentially detect the microbial mechanism of warming accelerated soil C loss observed in these forests. Further, studies have shown that short-term warming induces stress responses in bacteria triggering metabolic adjustments in protein productions (García-Descalzo et al., 2014) and growth vs. maintenance trade-offs by modifying their carbon use efficiency (CUE; Rodríguez-Verdugo et al., 2014, Pold et al., 2020). Taken together, these studies indicate the importance of studying the functional shifts in microbial communities experiencing long-term chronic warming. Inquiries as such can provide important insights on precise functional processes that dominate the energetic budgets of accelerated soil C loss in the face of climate change.

Along with long-term warming, annual temperature fluctuations due to seasonal variations in temperate forests can significantly alter the soil microbial community composition (Habekost et al., 2008; Cao et al., 2011) and the C allocation to soil microbes (Žifčáková et al., 2017). Few studies have looked into how seasonal variation in C allocation and microbial community composition can impact soil C dynamics (Lladó et al., 2017). Chronic soil warming was shown to have seasonal effects on soil moisture, respiration, and net nitrogen (N) mineralization (Contosta et al., 2011), suggesting that there is likely a seasonal component to chronic warming-accelerated soil carbohydrate degradation. Based on microarray analysis of soils from a grassland warming experiment, warming-associated gene expression was altered globally in warming treatments, with C degradation genes altered significantly with soil moisture (Xue et al., 2016). However, the relative importance of seasonality and its ultimate impact on microbially-accelerated soil C-dynamics in temperate forests is still largely unknown.

In this study, we applied techniques in meta-transcriptomics to identify changes in bacterial structure and function when exposed to both seasonal fluctuations and long-term simulated warming in Harvard Forest, MA. We first compared the magnitude of seasonal vs. long-term effects of warming on the active bacterial communities to understand long-term changes in C-energetics in this forest. Second, we identified metabolic pathways that showed a prominent “warming effect” to understand processes that facilitated the observed warming responses in these communities over time. Specifically, we tested the hypothesis that accelerated mineralization of recalcitrant C pools is due to a shift in the microbial community structure and functional potential. Finally, we used quantitative gene expression data to generate new hypotheses about the mechanistic nature of microbial adaptation of warming that extends beyond the expected increased utilization of soil carbon pools. Meta-transcriptomic data is limited to providing a relative account of the active transcripts, which can lead the libraries to be dominated by metabolic pathways shared by most cells instead of representing unique expression pathways characteristic of the given environmental conditions (Gifford et al., 2011). To overcome this, we employed the use of an internal standard to estimate the absolute number of transcripts in our samples compared to the number of transcripts sequenced (Gifford et al., 2011; Satinsky et al., 2013). Because of the previously observed changes in soil bacterial and fungal abundances, absolute measures of transcripts can help elucidate how gene expression might be linked to chronic warming.

Forests are significant global C sinks (Pan et al., 2011; Zhu et al., 2019), and understanding their ecology helps to manage and predict C cycling processes and their associated feedback to global climate. Microbial decomposition of soil organic matter is expected to increase in response to warming. The rates and mechanisms of these microbial driven soil C losses are, however, not fully understood. In this study, we used transcriptome (functional) data to identify changes in the structure and function of soil microbial communities in response to 30 years of simulated warming. Previously, we observed increased total soil carbon loss (Melillo et al., 2017; Pold et al., 2017), and evidence of increased microbial activity and an increasing role of bacteria in soil C degradation (Pold et al., 2015, 2017). However, meta-transcriptomes revealed an overall small effect of chronic warming on bacterial gene transcript abundances (Figure 2; Table 3). Although, we did not see a significant warming effect in major KEGG and CAZy classes (Figures 2, 3), several transcripts of enzymes related to carbohydrate, fatty acid, and lipid metabolism exhibited higher abundances in heated plots, primarily in the organic horizon (Figures 2C, 3C). These results are consistent with previous studies that showed an increased concentration of lipids (Pold et al., 2017) and overall decline in the quality and quantity of SOM with 3 decades of warming in these forests (Bradford et al., 2008; Pold et al., 2015, 2017). Pold et al. (2016) reported “warming effect” on similar carbohydrate-degrading enzyme genes relative abundances as seen in this study, but in mineral horizons in soils collected in Phase III when respiration was slightly higher in the heated plots (Melillo et al., 2017). They suggested an incomplete degradation and subsequent translocation of SOM from organic into the mineral horizon. Earlier depletion of labile C in the organic layer might have led to functional adjustments in these communities that are now able to access more recalcitrant SOM. This can subsequently result in complete or near-complete degradation of organic matter within the top layer, accounting for the “warming effect” in organic horizon observed here on soils collected in Phase IV (respiration similar or equal between plots; Melillo et al., 2017). Previously, it has been observed that compared to mineral horizon, organic soil exhibits greater shifts in structural and functional potential in response to environmental pressures (Cardenas et al., 2015; DeAngelis et al., 2015). Further, changes in SOM quality and quantity impacts the functional diversity in soil microbes to a greater extent compared to taxonomic diversity (Baldrian et al., 2012; Uroz et al., 2013). Taken together, the observed functional shifts in organic layer communities toward increased carbohydrate and lipid degrading potential not only supports our initial hypotheses of increased degradation of organic matter in the heated plots, but also provides evidence of a more responsive organic soil type to long-term warming.

Microbial C-dynamics appear to be dominated by different metabolic processes during different seasons and soil layers (See Supplementary Data: Appendices 5, 7), as evidenced by the fact that across both time points and horizons, the shifts in these transcript abundances are mostly unique. These differences in microbial processes can be due to variations in soil microclimatic parameters leading to transient changes in resource availability driven by diurnal, seasonal, and annual changes in these forests (Luo and Zhou, 2006; Edwards et al., 2007; Blankinship et al., 2011; Contosta et al., 2011; Zhang et al., 2011). Peng et al. (2018) found that changes in microclimatic parameters were the primary drivers of loss in soil organic carbon in a grassland ecosystem in the Qinghai-Tibet Plateau. Additionally, WGCNA analyses indicated a complex relationship between gene abundances and community function, where although several transcripts showed strong association to soil enzyme activities, no direct correlation was observed (Figure 5). This decoupling can be the result of microclimatic and temporal variability between samples, although, the broad substrate specificity of the KEGG enzymes cannot be ignored as a possible explanation for the same. Therefore, additional annual and seasonal time points should be sampled to further test whether the same processes dominate the microbial C-cycling consistently over time. Nevertheless, such spatial and temporal heterogeneity in microbial C-cycling genes can affect the functional capacities of soil overtime and is crucial to understand and predict warming-induced C losses in a warming climate.

While changes in gene expression based on chronic warming were not large, long-term warming caused a substantial decrease in microbial biomass (Frey et al., 2008; DeAngelis et al., 2015) and changes in abundance of certain microbial populations (DeAngelis et al., 2015; Pold et al., 2016). In this study, we used both mRNA and rRNA to identify changes in bacterial communities in response to warming (Figure 4). Our results showed that mRNA and rRNA annotations provided similar estimates of the dominant phyla, except that Actinobacteria were overrepresented in the mRNA sequences, possibly suggesting their increased metabolic activity relative to other groups. mRNA annotations were not able to classify the broad diversity that the rRNA was able to capture especially for rarer phyla such as Crenarcheota, Firmicutus, and Deltaproteobacteria. Cohen’s d values for mRNA and rRNA varied between soil types, which might be due to differences in annotation techniques used in either case (See Supplementary Data: Appendices 7, 8). Our data therefore suggests that compared to rRNA classification, taxonomic annotation of protein coding reads can only provide a cursory estimate of the abundant soil bacterial communities.

By separately analyzing rRNA vs. mRNA, we were able to detect changes in function that are separate from changes in community structure. Consistent with previous studies (DeAngelis et al., 2015; Pold et al., 2016), we observed relatively small changes in community structure with warming (Figure 1C), though both methods are biased toward over-representing bacteria. Actinobacteria increased, while the phyla Crenarchaeota, Planctomycetes, and Bacteroidetes decreased in relative abundance in the heated plots compared to controls (Figure 4). Changes in phyla abundance were only observed in the organic horizon, with a significant decrease in diversity in the heated plots. The three phyla depleted with long-term warming were all Gram-negative, and the enriched Actinobacteria were Gram-positive, supporting the previous observation of community shifts toward Gram-positive in the heated plots (Frey et al., 2008). The increased abundances of Actinobacteria might indicate a greater role of this taxa in driving the carbohydrate and lipid degradation in the heated plots. Although, we did not observe any differences in ligninolytic enzyme activities in our functional data, several members of Actinobacteria are associated with production of bacterial laccases that are important enzymes in lignin degradation (de Gonzalo et al., 2016). As such, their increase in the heated plots might imply a shift toward communities that are able to access more recalcitrant SOM, such as lignin and lignin-derived SOM, as labile C pools continue to deplete in response to warming. Members of Actinobacteria are generally over-represented in genomic databases; but, the consistent enrichment of Actinobacteria in heated plots in the present and other warming studies (Castro et al., 2010; Sheik et al., 2011; Shade et al., 2012) might indicate a greater tolerance to warming and its associated changes in soil physical parameters (Barnard et al., 2013).

Long-term warming in these forests has resulted in a decrease of fungal biomass with no substantial change in bacterial biomass in response to warming (Frey et al., 2008; DeAngelis et al., 2015). Enzyme activities were observed to be the same per gram of soil but decreased per unit biomass because of warming-induced decline in microbial biomass (Pold et al., 2017). Further, long-term warming also decreases soil moisture by lowering the soil water holding capacity (Rustad et al., 2001; Werner et al., 2020), which in turn can drive changes in microbial biomass and enzyme productions. While gene expression on a per gram soil basis may show small changes, as we present here, these changes are probably occurring against a backdrop of a smaller or less abundant microbial community in the heated plots. The use of internal standards was an attempt to distinguish between increased functional capacity vs. functional shifts in microbial communities by measuring absolute transcript copy numbers (Conant et al., 2011). However, we were unable to estimate the absolute transcript copy numbers as four of the 32 samples showed no recovery of internal standards (Table 1). This could be due to differences in sequencing artifacts, which can be avoided in the future by either increasing technical replicates or using a combination of multiple internal standards to increase robustness (van de Peppel et al., 2003). Nonetheless, measuring absolute transcript copy number in future functional studies will be valuable to estimate not only the role of under-represented but crucial metabolic pathways (Gifford et al., 2011), but also to evaluate the role of resilience and/or adaptation in functional potential of these communities in response to chronic long-term warming. To our knowledge, this is the first report on the application of internal standards to measure absolute transcript abundances in soil systems. Further refinement would be needed to determine if this method is applicable to soil systems.

Seasonal variation in abiotic factors, such as temperature and moisture content, can both directly and indirectly impact soil C dynamics (Benbi et al., 2014; Lladó et al., 2017; Žifčáková et al., 2017) and changes in soil microbial community (Habekost et al., 2008; Cao et al., 2011). We therefore predicted a seasonal impact on warming-accelerated microbial decomposition. We chose two time points, T2 (June) and T6 (October), based on Pold et al. (2017) that demonstrated a strong seasonal trend in extracellular enzyme activities, with most enzymes peaking during T2 and showing a decline or a 2nd peak during T6. Contrary to our expectations, we detected little effect of season on community response and community structure. In fact, the seasonal effect was the smallest factor driving global gene expression, carbohydrate-associated gene expression, and taxonomic marker transcripts (Table 2). Our results are consistent with another meta-transcriptome study from the same site that found no seasonal variation in community structure across six different time points (Rodríguez-Reillo, 2019). Previous studies have reported increased respiration and enzyme activities in summer (highest temperature), with lower microbial decomposition rates in spring and winter (lowest temperature; Contosta et al., 2011; Žifčáková et al., 2017). Winter communities decomposed more recalcitrant substrates such as lignin and cellulose, as compared to summer communities (Koranda et al., 2013; Žifčáková et al., 2017). To our knowledge, studies on seasonal variations on microbial gene expression (meta-transcriptome) are still lacking. The absence of a seasonal effect in this study is possibly due to the collection of soils at time points that were experiencing similar temperatures at the time of sampling (Pold et al., 2017). Further, seasonal variation in microbial activity is related to availability of different C-sources (Contosta et al., 2011), which can decrease under long-term warming due to depletion of labile C and reduction in microbial biomass (Frey et al., 2008). This, however, cannot explain the lack of seasonal variation in control plots. Furthermore, Žifčáková et al. (2017) indicated that there is temporal variation in bacterial vs. fungal gene expression related to decomposition of recalcitrant C-pools. Comparing the bacterial gene expression from this study with fungal gene expression (Romero-Olivares et al., 2019) across these time points might reveal seasonal patterns in microbial functions that are otherwise obscured. Additionally, extrapolating seasonal trends from gene expression measured at a single time point is perhaps not feasible given the observed diurnal variations in abiotic and biotic factors in temperate ecosystems (Zhang et al., 2011). In the future, gene expression should be measured and compared recurrently over extended periods to understand and capture changes in microbial functional potential in response to long-term chronic warming. Although, analyses of such vast sequencing data are technologically challenging, understanding, and capturing the nuances in microbial processes will help us better predict microbial driven C cycle-temperature feedbacks in the future.

In this study, we did a comparative meta-transcriptome analysis to understand the consequences of long-term warming on terrestrial soil bacterial communities and the related effects on soil C dynamics. Although, the warming effect on overall gene expressions was small, we did see increased expression of several enzymes involved in carbohydrate and lipid metabolism in heated plots, providing support to our main hypothesis of warming-induced acceleration of recalcitrant C-pool. Our results also show a strong effect of soil type on carbohydrate-associated gene expression, with greater warming effect on microbial transcription in the top organic horizon. We did not find any substantial changes in microbial community structure in response to warming. Both community structure and function showed no seasonal variation in our data. The meta-transcriptome approach allowed us to identify unique metabolic processes that controlled microbial decomposition and in turn soil C fluxes in different soil types in response to long-term warming. The lack of an overall seasonal gene expression pattern in this study highlights the need for comparing microbial functional profiles over an extended period of time to understand the contribution of these mechanisms to soil C-cycling in a chronically warmed forest. Information as such will substantially increase our predictive power to model correlated changes in microbial functions and temperature fluctuations in a warming climate.

The Next-seq RNA sequence data from this study have been uploaded to the DDBJ Sequence Read Archive (SRA) under the accession number PRJNA624774.

KD and JB conceived and designed the study. JM provided access to the study site and assisted with soil sampling. AB conducted RNA isolation from samples. PR and LA conducted library construction. PR performed Illumina sequencing and data analyses. SG performed the taxonomic analyses. JS performed part of the CAZy analyses. PR and KD wrote the paper. All authors contributed to the article and approved the submitted version.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.