2.1 Soil collection

Soils were collected from two long-term warming experiments at the Harvard Forest Long-Term Ecological Research (LTER) sites in Petersham MA USA (42°30′30″N, 72°12′28″W). The two experiments are located immediately adjacent to each other and had been warmed +5°C for 13 (Soil Warming and Nitrogen Addition Study [SWaN]; Contosta et al., 2011) or 28 years (Prospect Hill [PH]; Melillo et al., 2002). We sampled control and heated plots in SWaN, and disturbance control and heated plots in PH; these disturbance controls had heating cables buried as in the heated plots, but the power was never turned on. Experimental plots are 6 × 6 m at Prospect Hill, and 3 × 3 m at SWaN, and follow a randomized block design in the former and a completely randomized design in the latter.

The sites of the two warming experiments are co-located in the same forest stand, characterized by coarse-loamy inceptisols and the same dominant tree species: paper and black birch (Betula papyrifera and lenta), red maple (Acer rubrum), black and red oak (Quercus velutina and rubra), and American beech (Fagus grandifolia). SWaN plots are situated under a canopy gap relative to PH, resulting in slightly different understory vegetation between the two sites (Muth & Bazzaz, 2002). We pooled the control plot data from both sites for analysis since the control plots do not differ in the parameters measured here, with the exception of SOM lability index (I index). However, warmed treatment plots at SWaN and PH have not followed the same respiration trajectory over the course of the experiment, so we do not consider the field experiments as two stages in a chronosequence. Additional site details can be found in Table S1.

Soils were sampled on 15 July and 19 October 2019, timings which were chosen to reflect soil processes before and after autumn litter deposition (Munger & Wofsy, 2021). We observed that between July and October of 2019, most of the Acer trees had lost their leaves, resulting in an increase in forest floor litter material between the two time points.

Two 10 cm depth soil cores were collected from each of five plots per site, treatment, and time point using a 5.7 cm diameter tulip bulb planter. After removing and excluding the undecomposed forest floor litter material, these cores were separated into organic and mineral soil by color and texture. Soils were then sieved to 2 mm on site and transported back to the laboratory at ambient temperature, within 4 h of collection (Figure 1a). We ultimately processed and analyzed a total of 79 samples because the organic soil exceeded the length of the corer and prevented sampling of the mineral soil in one control SWaN plot in July. Immediately upon arrival at the laboratory, different subsamples of soil were dried to constant mass at 65°C overnight to determine soil moisture content or placed in the −80°C for metabolomic analysis. The remaining soil was left in plastic containers at room temperature (20°C) overnight. Samples which had greater than the target of 48% water holding capacity were left with the lids of the container ajar to allow them to dry slightly overnight (n = 2 samples). Soil was then weighed for the CUE, respiration, enzyme assays, microbiological assays, and SOM analyses described below (Figure 1).

2.2 Microbial community catabolic potential

We used BioLog Ecoplates to evaluate the intrinsic potential growth rate and substrate preferences of microbes in our soil samples. Within 24 h of collecting soils, we placed 0.25 g (organic) or 0.5 g (mineral) soil in a 125 ml Erlenmeyer flask with 50 ml of 0.9% sodium chloride and mixed on a shaker table at 180 rpm for 30 min. The flasks were then left for 5 min to allow the particles to settle, and then 100 μl of supernatant was transferred to each well of the EcoPlate. The initial absorbance of each plate was measured at 590 and 750 nm and after every 4–12 h afterwards to capture the full growth curves. The measurements at 750 nm were subtracted from the measurements at 590 nm to differentiate the change in dye reduction through time from the change in cell growth. Then the amount of dye color development in the water only well was subtracted out to quantify the metabolic activity due to the substrate itself. Negative control plates containing just buffer were used to verify the assays were not contaminated. We subsequently fit growth curves to both individual substrates and to the average color across the plate (mean well color development) using the growthcurver package (Sprouffske & Wagner, 2016). Growth rates not significantly different from zero (p > .05) were set to zero before all substrate-specific growth rates were used for calculating substrate class-specific growth rates and for principal coordinates analysis. Principal coordinates analysis was completed using the vegan package (Oksanen et al., 2019) with a Hellinger-transformed relative growth rate matrix rather than the raw growth rate data so that relative substrate preference could be separated from the overall differences in growth rate under the assay conditions. We refer to the first axis of this ordination as the “ecoplate substrate diversity.”

2.3 Extracellular enzyme activity

Between the collection and completion of extracellular enzyme assays (within 4 days), we stored soils at 15°C, which was near the middle of field temperature range. We assayed the cellulose-degrading enzyme beta-glucosidase (BG), the chitin- and peptidoglycan-degrading enzyme N-acetylglucosaminidase (NAG), and the oxidative enzyme pool (phenol oxidase + peroxidase; OX) at 4, 10, 15, 20, 25, and 30°C, covering most of the range of growing season temperatures experienced by heated and control plot soil communities at these sites. These enzymes were selected for their propensity to target chemically similar (BG) or chemically diverse (OX) bonds, as well as microbial necromass (NAG) (Burns & Dick, 2002).

Slurries were prepared with 1.25 g fresh wet weight soil in 175 ml 50 mm pH 4.7 sodium acetate buffer, using a Waring blender set to high for 1 min. The slurry was stirred at 300 rpm during pipetting to ensure an even distribution of soil. Fourteen technical replicates of 200 μl slurry were pipetted into black plates for BG and NAG, and 500 μl into deep well 2 ml plates for total oxidative enzyme pool. Plates were placed at the assay temperature for at least 25 min before substrate addition to allow them to reach the specific target temperature. Fifty microliter of 4000 μm 4-methylumbelliferyl β-d-glucopyranoside or 2000 μm 4-methylumbelliferyl N-acetyl-β-d-glucosaminide, or 500 μl 25 mm L-DOPA +0.03% H₂O₂ were added to each well, which we previously deemed to be sufficient to attain the maximum activity rate. Each hydrolytic plate contained a slurry-only control (200 μl slurry, 50 μl buffer) and a standard curve (50 μl 0–500 μm in twofold dilution, in 200 μl buffer). A separate plate for mineral and organic soil quenches (200 μl slurry, 50 μl standard) and substrate-only controls (200 μl buffer, 50 μl substrate) was prepared for each temperature.

Hydrolytic assay plates were measured at an excitation/emission wavelength pair of 360/450 nm immediately following substrate addition and after an additional 2, 4, and 8 h. Greatest activity was observed between 0 and 2 h, so we used data from this period in our analyses. A 100 μl aliquot of the supernatant from oxidative assay plates was removed to a clear polystyrene plate after 4 h of incubation and read at 460 nm. We used a SpectraMax M2 plate reader controlled by the SoftMaxPro v. 5.4 software for all measurements. Extracellular enzyme rates were calculated as previously (German et al., 2011) and normalized by the sample-specific microbial biomass carbon (MBC) measurement (described below).

2.4 Respiration, growth, and CUE

We used the substrate-independent (H₂¹⁸O-CUE) method (Spohn et al., 2016) to evaluate soil microbial CUE. Weighing for growth and respiration measurements began the day following soil collection. Three 0.3 g (mineral) or 0.15 g (organic) replicate aliquots of soil were weighed into small vials, placed in a larger tube, and sealed with parafilm to prevent additional moisture loss. Once all samples were weighed (2 days after soil collection), water was added to bring them to 60% water holding capacity. Two replicates received water so that 20% of the total water was present as ¹⁸O-water, while the remaining replicate received all water as ¹⁶O-water to account for natural abundance ¹⁸O. Tubes were stoppered with a neoprene bung and immediately placed in the appropriate incubator (4, 10, 15, 20, 25, or 30°C). Empty tubes were also sealed after every 5–10 tubes in order to measure the starting CO₂ levels in tubes. After 24 (15–30°C) or 48 (4 and 10°C) hours, the CO₂ was measured in the tube using a 30 ml headspace sample injected into a Quantek instruments infrared gas analyzer with 10 ppm sensitivity. We used a longer time for the lower temperatures because preliminary incubations indicated that it was necessary to reliably detect respiration. The soil samples were placed at −80°C until DNA extraction. We also incubated larger quantities of soil in order to validate respiration measurements from the CUE incubations. One gram (organic) or 2 g (mineral) soil were brought to 60% water holding capacity with ¹⁶O-water, and the tubes were incubated alongside the CUE samples.

CUE measurements using the (H₂¹⁸O-CUE) method estimate the new microbial biomass produced during the incubation period based on ¹⁸O-DNA at the end of the incubation period. DNA was extracted from all soils incubated with ¹⁸O-water and a subset of soils incubated with ¹⁶Owater using the Qiagen Powersoil kit. Technical duplicates were pooled before quantification using Qubit. The ¹⁸O enrichment of the DNA was measured using TC/EA-IRMS (Delta V Advantage, Thermo Fisher) at the UC Davis Stable Isotope Facility. CUE was calculated as per Spohn et al. (2016) but using a sample-specific conversion factor rather than the overall average because large differences in MBC:DNA ratios across community types can bias CUE measurements (Pold, Domeignoz-Horta & DeAngelis, 2020).

2.5 Temperature response curve fitting

We anticipated that respiration and extracellular enzyme activity would increase non-linearly with incubation temperature. Therefore, we compared fits of the Taylor exponential model (Lloyd & Taylor, 1994; ln[rate] = ln[a] + b × Temperature), and the Ratkowsky square root model (Ratkowsky et al., 1982; sqrt[rate] = b[T−T₀]) to these data. Models for respiration were fit to individual soil samples using the lm() function in R where the response variable (i.e., respiration or extracellular enzyme activity) was either natural log- or square root transformed. We extracted the slope (i.e., temperature sensitivity; b in equations above) estimates from the model for each sample to compare among experimental factors.

Visually, the Ratkowsky model consistently underestimated the respiration at higher temperatures compared to the Taylor model, and the R² statistic was also consistently lower (median R² = .91 [range 0.65–0.98] versus 0.96 [range 0.66–0.99], median difference = 0.052; Table S3). Therefore, we proceeded to extract the slope parameter from the Taylor model and used this as a metric of temperature responsiveness for downstream statistical analyses for both respiration and extracellular enzyme activity. For BG, the R² for the Taylor model varied from 0.42 to 0.99, with a median of 0.85. For NAG, R² varied from 0.04 to 0.99 with a median of 0.85 ( Table S4). Since all R² below .4 were restricted to October, we excluded NAG activity from this time for our analysis. The Taylor model similarly showed a poor fit to the oxidative enzyme activity (median R² = .88, range .00–.99) with instances of poor fit seen in data from both time points. Therefore, we did not further analyze oxidative enzyme activity temperature sensitivity (Table S4). The growth and CUE data were very variable compared to the respiration data and individual samples did not follow a consistent pattern. Therefore, we could not fit a curve to the data and instead used ranked CUE and growth rates at different temperatures for each soil sample. We considered the CUE temperature optimum to be the temperature where we measured its highest value (Bölscher et al., 2020).

2.6 Quantitative PCR

The abundance of total fungi and total bacteria was assessed using quantitative PCR (qPCR) with ITS primers (ITS1: 5′-TCCGTAGGTGAACCTGCGG-3′ and 5.8 S: 5′-CGCTGCGTTCTTCATCG-3′)and 16 S ribosomal RNA primers (Eub338: 5′-ACTCCTACGGGAGGCAGCAG-3′ and Eub518: 5′ATTACCGCGGCTGCTGG-3′) (Fierer et al., 2005), respectively. The abundance in each soil sample was based on increasing fluorescence intensity of the SYBR Green dye during amplification. The qPCR assay was carried out in a 15 μl reaction volume containing 2 ng of DNA, 7.5 μl of SYBR Green PCR master mix (Qiagen quantifast SYBR kit), and each primer at a concentration of 1 μm. Inhibition tests were performed by running serial dilutions of DNA extractions and did not indicate inhibition of amplification. For each sample at least two independent qPCR assays were performed for each gene with technical duplicates within each assay. The qPCR efficiencies for both genes ranged between 78% and 110%. Values are reported as gene copy number g⁻¹ dry soil.

2.7 Microbial biomass carbon

Quadruplicate 0.5 g (organic) or 1 g (mineral) soil samples for MBC were weighed out with the CUE samples. We kept the weighed out soil samples at 15°C (the middle of the temperature range used for CUE and enzyme measurements) until the end of CUE incubations so growth and biomass measurements could be taken on similarly treated soils. Two replicates were exposed to chloroform vapor fumigation under vacuum for 24 h, while the other two replicates were placed at 4°C for the duration. Soils with and without chloroform were subsequently extracted in 15 ml of 0.05 m K₂SO₄, which a preliminary trial indicated led to similar levels of extractable organic carbon and MBC detected as the more standard 0.5 m K₂SO₄ for our soils. DOC was measured on a Shimadzu TOC analyzer.

2.8 SOM quality

We used ramped thermal rock-eval® pyrolysis (RE) to evaluate SOM quality (Soucémarianadin et al., 2018). During RE, carbon oxides are quantified as they come off a soil sample subject to increasing temperatures, thereby providing a metric of SOM intrinsic thermal stability. Compounds with high thermal stability include aromatic and phenolic non-lignin compounds, while lipids and polysaccharides tend to have lower thermal stability (Sanderman & Grandy, 2020). Mineral and organic soils were dried at 65°C and crushed to a fine powder in a mortar and pestle. Between 50and 70 mg soil was pyrolyzed over a temperature ramp from 200 to 650°C, followed by combustion to 850°C using a rock-eval 6 pyrolyzer (Vinci technologies) at the Institute of Earth Sciences of the University of Lausanne (Switzerland). Hydrocarbons released during this process were measured by a flame ionization detector. The resultant thermogram was used to calculate the I index (“labile carbon fraction”) and R index (“recalcitrant carbon fraction”) as previously (Sebag et al., 2016). We also used the thermogram of each sample to calculate a Bray–Curtis distance matrix of all samples, taking the peak height at each 1°C increment as an input value to create the distance matrix between all samples (Domeignoz-Horta et al., 2021). This approach allowed us to evaluate SOM composition in addition to C quantities in the different soils. The first axis of the NMDS was used as a proxy for SOM quality as it was strongly correlated with the I index (R² = .89, p < .0001).

2.9 Metabolomics

Water-extractable polar metabolites (here, “water-extractable organic matter,” or WEOM) were obtained by shaking 1 g dry soil equivalent with 5 ml (mineral soil) or 10 ml (organic soil) of LC/MS-grade water at 200 rpm and 4°C for 1 h. Every extraction was accompanied with a water-only control. The soil extracts and the water controls were filtered through a 45 μm PTFE filter. The filtrate was lyophilized and sent for analysis at the DOE Joint Genome User Facility at Lawrence Berkeley National Laboratory. Polar metabolites were resuspended in 170 μl of 100% methanol containing 13C-15N labeled amino acids (30 m, 767964, Sigma). UHPLC normal phase chromatography was performed using an Agilent 1290 LC stack, with MS and MS/MS data collected using a Q Exactive HF Orbitrap MS (Thermo Scientific). Full MS spectra were collected from m/z 70–1050 at 60 k resolution in both positive and negative ionization mode, with MS/MS fragmentation data acquired using stepped 10, 20, and 40 eV collision energies at 17,500 resolution. Mass spectrometer source settings included a sheath gas flow rate of 55 (au), auxiliary gas flow of 20 (au), sweep gas flow of 2 (au), spray voltage of 3 kV, and capillary temperature of 400°C. Normal phase chromatography was performed using a HILIC column (InfinityLab Poroshell 120 HILIC-Z, 2.1150 mm, 2.7 m, Agilent, 683775-924) at a flow rate of 0.45 ml min⁻¹ with a 2 uL injection volume. To detect metabolites, samples were run on the column at 40°C equilibrated with 100% buffer B (99.8% 95:5 v/v ACN:H₂O and 0.2% acetic acid, w/ 5 mm ammonium acetate) for 1 min, diluting buffer B down to 89% with buffer A (99.8% H₂O and 0.2% acetic acid, w/ 5 mm ammonium acetate and 5 m methylene-di-phosphonic acid) over 10 min, down to 70% over 4.75 min, down to 20% over 0.5 min, and isocratic elution for 2.25 min, followed by column re-equilibration by returning to 100% B over 0.1 min and isocratic elution for 3.9 min. The injection order of the 79 experimental samples and 5 extraction controls was randomized and an injection blank (2 uL of 100% MeOH) run between each sample. Metabolomics data were analyzed using Metabolite Atlas (Yao et al., 2015) and with in-house Python scripts to obtain extracted ion chromatograms and peak heights for each metabolite.

Metabolite identifications were verified with authentic chemical standards and validated based on three metrics comparing (1) detected versus theoretical m/z (<5 ppm), (2) retention time (≤0.5 min), and (3) fragmentation spectra similarity to a chemical standard run using the same chromatography and LC–MS/MS method. Data from internal standards and quality control samples (included throughout the run) were analyzed to ensure consistent peak heights and retention times. Peak height of different samples was calculated out of four replicates and we excluded metabolites that showed a quality score under 1 and were not measured in every sample. Where multiple ions existed for the same compound within a sample analyses we summed the peak heights. Where compounds appeared in both the positive and the negatively charged compound summary table, we used the greatest value. After quality control we had 227 compounds. The WEOM matrix was normalized by sum and pareto scaled before calculating Bray–Curtis distance and completing principal coordinates analysis using the MetabolAnalyze, vegan and ape packages, respectively (Nyamundanda et al., 2010; Oksanen et al., 2019; Paradis & Schliep, 2019). Pareto scaling divides each compound by the square root of the variance. We selected this scaling approach to moderately reduce the importance of the extremely dominant compounds and remove the large differences in total peak height between mineral and organic soil samples without completely abolishing all relative abundance information.

Compound complexity was determined based on the Bertz/Hendrickson/Ihlenfeldt complexity reported in PubChem. The complexity score is higher in polymers compared to monomers of the same compound, and in asymmetric and/or heteropolymeric compared to symmetric and/or homopolymeric compounds. Therefore, the complexity score can be considered an imprecise proxy for how difficult a compound is to break down. These values were downloaded from PubChem. Mean nominal oxidation state was calculated as per (LaRowe & Van Cappellen, 2011) (NOSC = -((−Z + 4C + H−3 N−2O + 5P -2 S)/ C) + 4). NOSC was selected as a proxy for the quality of DOC because oxidation of compounds with higher NOSCs is associated with greater electron transfer (LaRowe & Van Cappellen, 2011), while compounds with lower NOSCs have greater potential energy but require a larger investment to acquire it. As such, the NOSC of SOM generally decreases as decomposition proceeds (Gunina & Kuzyakov, 2022). The NOSC and complexity indices for each compound are reported in Table S5.

2.10 Metagenome sequencing and annotation

Metagenomes were obtained from organic soil collected in July from both sites, totaling 20 samples. DNA extracted from soils incubated at 20°C was sequenced at the Joint Genome Institute (Methods S1 and Table S2). Reads were assembled and PFAMs were annotated in the assembled fraction of reads using the DOE-JGI Metagenome Annotation Pipeline (Huntemann et al., 2016). This led to an average of 1.72 Gbp assembled and 2,807,085 genes predicted per metagenome. We then mapped PFAMs to CAZys (Pold et al., 2016) and used mean mapping depth for CAZys relative to the mapping depth to the single copy rpoB gene (PF04563, PF04561, PF04565, PF10385, PF00562, PF04560) as our proxy of CAZy abundance across samples. We elected to focus on CAZymes as they are responsible for the polymerization and depolymerization of many of the dominant compounds found in leaf litter and microbial cell walls, including cellulose, lignin, and chitin. We subsequently categorized CAZymes based on the primary substrates they decompose (Berlemont et al., 2015).

2.11 Data analysis

Analysis was completed in R version 3.6.3 using packages reshape2 (Wickham, 2007), vegan (Oksanen et al., 2019), lme4 (Pinheiro, et al., R Core Team, 2019), agricolae (Mendiburu, 2019), grofit (Kahm et al., 2010), dplyr (Wickham et al., 2022), lattice (Sarkar, 2008), MASS (Venables & Ripley, 2002), and permute (Simpson, 2022). Variables were log-transformed where necessary to fit model assumptions. We calculated Hedge's G (Cohen, 1988) using the effect size package (Torchiano, 2020) to compare the effects of warming treatment and season on biotic and abiotic soil properties; we considered drivers to have a significant effect when the 90% confidence intervals surrounding the effect size estimate did not overlap zero (Figure 2). We elected to use an alpha of 0.1 (and 90% confidence intervals) rather than 0.05 (or 95% confidence intervals) because we were limited by the number of field replicates we could collect and a power analysis using means for microbial biomass, soil C stocks, and extracellular enzyme data collected from the site in 2014 (Pold et al., 2017) indicated this would be needed to detect changes in mean values of 25% with a power of 80%. Figures were made using the ggplot2 package (Wickham, 2016).

2.12 Structural equation modeling

We used structural equation modeling (SEM) to test the hypothesis that long-term warming affects microbial CUE and temperature sensitivity of respiration via changes in C quantity and quality. The hypothesized path structure for the temperature sensitivity of respiration was based on the premise that long-term warming, season and soil depth had direct effects on C quantity, quality and microbial communities (fungal: bacterial ratio). These variables are then hypothesized to affect the size and temperature sensitivity of the extracellular enzyme pool and the composition of the water-extractable SOM. Finally, we hypothesized that the temperature sensitivity of respiration would be controlled by the enzyme pool, the enzyme temperature sensitivity, and the quantity and quality of soil C (Figure S8). We also completed SEM analysis on microbial CUE at 20°C. This model had a similar overall layout of the temperature sensitivity of respiration model, except that we additionally hypothesized that diversity of EcoPlate substrate consumption would be affected by season and carbon quality, and that respiration and growth would be impacted by carbon quantity and quality (Figure S9).

The SEM model path fit was performed using the piecewiseSEM package (Lefcheck, 2016). Cumulative enzyme activity and temperature sensitivity are represented in the model as composite variables (Lefcheck, 2016), using the mean of BG, NAG, and oxidative enzymes for the former and BG alone for the latter. We kept the model that explained the most variation in temperature sensitivity of respiration and CUE, and had a non-significant Chi-squared test (p > .05), low Akaike information criterion (AIC) and high Comparative Fit Index (CFI >.9). If the test of direct separation (Lefcheck, 2016) identified missing paths in our hypothesized model we added these paths into the model.

3.1 Long-term warming seasonally affects soil C substrates

Warming reduced SOM quality and quantity by most measures, but this effect was principally visible during the July sampling (Figure 2). We evaluated the chemistry of soil DOM by analyzing the polar water-extractable fraction of soil organic matter (WEOM) using metabolomics, and the thermostability of SOM using rock-eval®.

The nominal oxidation state of compounds (NOSCs) identified in the WEOM varied from −1.5 to 3, although the mean value was less than zero for all samples (Figure S7). The effect of warming on NOSC was site and time point dependent, only increasing with warming in SWaN plots in July (Figure 2, Figure S7). Warming and sampling time point did not affect the molecular complexity of the WEOM, but it was greater in mineral compared to organic soil (Figure 2a, Figure S7, mean of 143 vs 135). Overall dissimilarity in WEOM composition was primarily explained by horizon (Adonis R² = .46; ADONIS F(1,78) = 87.50, p < .001), followed by warming treatment (Adonis R² = .083; F(2,78) = 7.82, p < .001) and time point (Adonis R² = .045; F(1,78) = 8.46, p < .001). There were no significant interactions between these variables, and no evidence of heterogeneous group dispersions (i.e., betadisper p > .05).

Long-term warming reduced the total dissolved organic carbon in the mineral soil at the 28-year-old warming experiment at both time points, and reduced the lability index of the total SOM pool in July in both field experiments (Figure 2). Warming explained just 2.2% of the variation in SOM chemical composition based on ramped thermal rock-eval® (F(1,76) = 7.30, p = .006), with much of the remaining variation explained by soil horizon (Adonis R² = .71, F(1,76) = 231.22, p < .001) and, to a lesser degree, season (R² = .020, F(1,76) = 6.76, p = .008) (Figure 2). However, the effect of warming depended on horizon (warming × horizon interaction R² = .011, F(1,76) = 3.76, p = .04) and sampling time point (sampling × warming interaction R² = .008, F(1,76) = 2.70, p = .09). The first axis of this ordination was strongly positively correlated with the thermal labile SOM (I index) (cor: 0.89 p < .001).

3.2 Long-term warming affects microbial communities and activity

Effects of warming on microbial communities and activity were most apparent in July (Figure 2). MBC was reduced by up to 45% in July (Figure 2a and Figure S1). Warming also reduced the July respiration at 20°C in mineral soils at both field experiments, and in organic soil samples taken from the 13-year old field experiment (Figure 2). Oxidative enzyme activity per gram of soil was also reduced by warming in the mineral soil of both sites in July, but oxidative and BG enzyme activities per unit microbial biomass were enhanced by warming in organic soils of the older field experiment (Figure 2b). Based on cultivation in EcoPlates, potential growth rate of microbial communities was unaffected by warming, except for a reduction in the organic horizon of the 13-year-old site in July (Figure 2b).

To understand how microbial communities might adapt to long-term warming in these soils we measured the temperature response of select microbial processes between 4 and 30°C. We only observed the expected increases in process rates with incubation temperature for respiration and extracellular enzyme activity (Figure 3a,b; Figures S4–S6). The temperature sensitivity of respiration was reduced by long-term warming in July in both long-term warming experiments and soil horizons (Figures 2, 3a and Figure S5).

The large variation in the ¹⁸O signal detected in DNA meant that growth and CUE data did not fit the expected unimodal or monotonic decline pattern, and we were unable to fit quadratic curves to the data to assess subtle shifts in the temperature optimum of growth or CUE (Figure S4). However, median CUE was greatest at 15°C across July samples and 15 or 20°C in October samples (Figure S4).

3.3 Microbial activity reflects SOM chemistry

We used correlations between microbial activity per unit microbial biomass and soil properties to evaluate changes in microbial physiology. This allowed us to evaluate how processes rate and temperature sensitivity of microbial activity are associated with SOM chemistry (Figure 4). The temperature sensitivity of microbial activity was generally positively correlated with the concentration of total and dissolved SOM, and the thermal lability of SOM (I index), and negatively correlated with the NOSC and chemical complexity of WEOM (Figure 5a). Microbial biomass-specific respiration rate at 20°C showed similar correlations to SOM indices as the temperature sensitivity, but the biomass-specific rates of BG activity showed contrasting patterns. Additionally, mass-specific oxidative enzyme activity was greatest in soils characterized by chemically complex and high NOSC WEOM and with low SOM lability indices (Figure 5a). We also observed a significant positive correlation between the temperature sensitivity of respiration and the temperature sensitivity of BG activity (Figure 5b), indicating that soils with a pool of cellulolytic enzymes capable of responding to increasing temperatures also respired more in response to changing temperatures.

3.4 Effects of warming on microbial physiology play out through its effects on soil C

We used path analysis to integrate the relationships between the abiotic and biotic measurements performed here and test the hypothesis that long-term warming drives microbial activity through its effect on soil C. The path analysis confirms this hypothesis and shows that long-term warming treatment, season and soil horizons influenced the carbon quantity and quality in the soil while also affecting respiration temperature sensitivity (Figure 6). However, long-term warming duration has no direct effect on microbial physiology or WEOM and SOM quality. Unexpectedly, the nominal oxidation state of water-extractable organic matter was not a key integrator of shifts in SOM chemistry and did not drive respiration temperature sensitivity (Figure 6). Rather, temperature sensitivity of respiration was directly correlated with quantity of carbon and the size of the extracellular enzyme pool and temperature sensitivity of BG.

Season had a direct effect on the dissolved fraction of SOM and C quality while increasing soil depth decreased the SOM quality (e.g., labile fraction). Respiration temperature sensitivity was driven directly by the enzyme pool in soils, carbon quantity and the enzyme temperature sensitivity. The soil potential extracellular enzyme activity was driven by C quantity and season while the enzyme temperature sensitivity was driven by both season and carbon quantity. Although microbial communities can modulate thermal acclimation, the fungal to bacterial ratio was only driven by soil horizon and did not influence other parameters in our model. Our SEM model explained 27% of the variance in the temperature sensitivity of respiration, although many of the drivers of temperature sensitivity were explained to a much greater degree.

Path analysis shows that CUE at 20°C is directly affected by respiration and growth as expected, but also by other factors such as the dissolved SOM fraction and season, with higher CUE measured in October compared to July (Figure 7, Figure S4). The mass-specific microbial respiration component of CUE was much better explained than the growth component (R² = .77 vs. .14), with 86% of the variation in CUE explained overall (Figure 7). This was expected given the large variation in our measurements of growth rate. Distinct factors also drove respiration and growth differently. The mass-specific enzyme pool only drives respiration while carbon quality correlates positively with mass-specific respiration and negatively with mass-specific growth.

In order to avoid the potential problems associated with predicting CUE from its constituent components (i.e., growth and respiration), we also built a SEM without these variables to evaluate how much the other variables we measured directly explained CUE (Figure S10). This model captures 39% of CUE variance and the direct drivers were season, mass-specific enzyme pool, and carbon quality. High-quality carbon is associated with lower CUE in both models because it accelerates mass-specific respiration while slowing mass-specific growth (Figure 7 and Figure S10).

4.1 Long-term warming impact on SOM chemistry

Here, we evaluated the hypothesis that the effect of warming on microbial physiology plays out through its effects on the quantity and quality of SOM. We found that while soil chemistry had predictable correlations with microbial activity, the effects of warming on SOM chemistry—and, in turn, microbial physiology—were mostly limited to summer. Therefore, the effect of warming on microbial physiology depends on its effects on SOM chemistry.

The effect of warming on SOM quality was most apparent in the mineral soil, where soil lability index was reduced (Figure 2b). This is consistent with previous studies at these sites which show that warming-associated changes in soil organic matter quality are especially apparent in the mineral soil and typical plant-derived compounds such as lignin and polysaccharides were more depleted after 4 and 23 years of warming (Pisani et al., 2015; Pold et al., 2017). Greater decomposition with warming has also increased overall relative abundance of lipids (Pec et al., 2021; Pisani et al., 2015; Pold et al., 2017; VandenEnden et al., 2021) despite a reduction in plant lipids (Pisani et al., 2015), indicating further microbial processing. Lignin and lipids are both among the least oxidized (LaRowe & Van Cappellen, 2011) and most readily pyrolyzed (Sanderman & Grandy, 2020) compounds found in soil, such that a potential shift toward being more lipid-rich and lignin-poor SOM under warming in the organic soil would not be detected. However, it is more probable that any signal of increased microbial processing was drowned out in the high plant input: microbial processing ratio in the organic soil. This is consistent with the observation that warming does not consistently alter SOM quality throughout seasons, even in the more microbial necromass-dominated mineral soil (Figure 2b and Figure S1). Rather, differences were only apparent in our summer sampling when belowground inputs likely dominated carbon inputs (Abramoff & Finzi, 2015; Cheeke et al., 2021), and not following the autumn leaf litter influx. This is consistent with a wealth of studies indicating that season strongly regulates the input of plant material and nutrients in temperate forest soils and that this seasonality is associated with variation in microbial activity and biomass (Cheeke et al., 2021; Kaiser et al., 2010; Pold et al., 2017; Rasche et al., 2011; Wardle et al., 2004).

4.2 SOM chemistry impact on microbial thermal acclimation, biomass, and activity

The seasonal dependence of warming effects on SOM chemistry was further associated with temporally inconsistent effects of warming on microbial physiology. Limited substrate availability has been cited as necessary condition for microbial community physiological acclimation under warming (Alster et al., 2022; Bradford, 2013; Fierer et al., 2006; Kirschbaum, 2004; Luo et al., 2001; Moinet et al., 2021). We only observed apparent thermal acclimation of C cycling processes when the input of fresh plant material was limited (i.e., in summer), and not once a pulse of fresh aboveground plant input reaches the soil (i.e., in autumn; Figures 2b and 3). Temperature is an important driver controlling how soil microbes degrade organic matter (Bradford, 2013), and to a certain degree elevated temperatures may help microbes overcome substrate limitation by accelerating the rate at which extracellular enzymes process SOM. However, as warmer temperatures sustain accelerated microbial activity through time, substrates ultimately become depleted and the temporary warming-induced increase in respiration ceases (Luo et al., 2001; Romero-Olivares et al., 2017; Walker et al., 2018). This substrate depletion leads to a reduction in microbial biomass followed by a decrease in microbial activity (Walker et al., 2018). However, as seen here, seasonal flushes of fresh aboveground litter inputs may temporarily mitigate this warming-induced substrate depletion, allowing microbial biomass and therefore process rates to somewhat recover to the levels observed in the absence of chronic warming. This concept of substrate limitation induced reductions in microbial activity via reductions in microbial biomass is supported by the observation that biomass-specific C processes rates are less affected by warming compared to the activities measured on a per gram of soil basis (Figure 2b). Accordingly, two metrics strongly correlated with microbial biomass—total organic carbon and size of the extracellular enzymatic pool—along with temperature sensitivity of extracellular enzyme activity, were the main drivers of respiration temperature sensitivity (Figure 6).

Part of the greater respiration temperature sensitivity with larger enzyme pools may be linked to greater temperature sensitivity of extracellular enzyme activity, enabling greater supply rate of products for microbial metabolism without greater investment in those enzymes. However, the connection between the temperature sensitivity of respiration and extracellular enzyme activity is not equally linked to all SOM pools equally; soil samples with a bigger pool of cellulolytic enzymes were more capable of responding to increasing temperatures and consequently respired more in response to changing temperatures (Figure 6). Therefore, our results implicate the relevance of SOM quantity in driving microbial community activity (Kirschbaum, 2004; Moinet et al., 2021) and functioning (Domeignoz-Horta et al., 2021).

4.3 Relationships between SOM quantity and quality to microbial CUE

Carbon use efficiency is an important integrator of C quality with microbial metabolism, both because of the costs associated with extracellular enzyme production necessary to access those substrates (Domeignoz-Horta et al., 2020; Malik et al., 2019), and because of differences in the communities able to use and the energetic yield associated with the use of those different substrates (Pold, Domeignoz-Horta, Morrison, et al., 2020). The respiration component of CUE was more responsive to changes in SOM and WEOM pools compared to growth, indicating that growth has a relatively fixed cost compared to respiration. This is in line with the observation that microbes may downregulate the protein production machinery in response to warming-induced substrate depletion, allowing them to sustain growth (Sollinger et al., 2022) while the respiration component could then absorb changes in substrate quantity and quality. This downregulation of the protein machinery does not require a change in the active community (Sollinger et al., 2022). Changes in the microbial community have historically been variable and limited at our study sites (DeAngelis et al., 2015; Frey et al., 2008; Pec et al., 2021; Roy Chowdhury et al., 2021), and we did not observe a change in the F:B ratio under warming here (Figure S3). Nor did we observe a correlation between F:B ratio and CUE. This reinforces the idea that substrates quality and quantity rather than microbial identity have the dominant effect on CUE and its response to long-term warming in these soils (Figures 5, 7 and Figure S10).

Microbial activity measurements are more strongly correlated with the thermal lability (I index) of the total organic matter measured by rock-eval® than the water-extractable pool measured by polar metabolomics. It is uncertain whether this is due to the latter method missing the compounds most important for microbial activity over the course of the assays period, or whether the chemistry of bulk SOM rather than the water-extractable component is truly a better predictor of microbial physiology. Nonetheless, our results show that SOM thermal-lability index based on rock-eval® captures organic matter pools which vary between soils characterized by divergent microbial community metabolisms. The warming impact on resource availability is inconsistent due to seasonality in inputs (Figures 2 and 3), which might result in a variable selection regime (McDonald, 2019) limiting the degree of intrinsic thermal acclimation of the microbial community to long-term warming in these soils. As such, our results are consistent with previous findings that the response of microbial communities to global warming and its potential impact on soil carbon stocks will depend on concurrent changes in soil organic matter chemistry and availability to microorganisms (Bradford et al., 2008; Melillo et al., 2017). While we still lack a clear understanding of the controls that drive variation in the response of soil carbon to warmer temperatures globally (Van Gestel et al., 2018), our results identify the direct drivers of apparent thermal acclimation of C processes drawing attention to the importance to monitor soil C chemistry to better understand possible soil C feedback mechanisms to the atmosphere as the world warms.