Research ArticleGEOCHEMISTRY

Five years of whole-soil warming led to loss of subsoil carbon stocks and increased CO2 efflux

See allHide authors and affiliations

Science Advances  21 May 2021:
Vol. 7, no. 21, eabd1343
DOI: 10.1126/sciadv.abd1343


Subsoils below 20 cm are an important reservoir in the global carbon cycle, but little is known about their vulnerability under climate change. We measured a statistically significant loss of subsoil carbon (−33 ± 11%) in warmed plots of a conifer forest after 4.5 years of whole-soil warming (4°C). The loss of subsoil carbon was primarily from unprotected particulate organic matter. Warming also stimulated a sustained 30 ± 4% increase in soil CO2 efflux due to increased CO2 production through the whole-soil profile. The observed in situ decline in subsoil carbon stocks with warming is now definitive evidence of a positive soil carbon-climate feedback, which could not be concluded based on increases in CO2 effluxes alone. The high sensitivity of subsoil carbon and the different responses of soil organic matter pools suggest that models must represent these heterogeneous soil dynamics to accurately predict future feedbacks to warming.


Substantial uncertainty about Earth’s 21st century climate comes from the lack of evidence for whether warming-induced increases in soil respiration will be sustained over time and result in a net loss of soil carbon (1, 2), thereby creating a positive carbon-climate feedback. For example, forests are currently a net sink of carbon and their potential to absorb atmospheric CO2 over the coming decades has been touted as a powerful method of carbon sequestration (3). However, these projections ignore the impact of future climate change on forest soil carbon storage and on subsoil (>20 cm) carbon, in particular. Subsoils contain roughly half of global soil organic carbon (SOC) (4) and differ from surface soils in plant inputs, microbial communities, organic matter stabilization processes, and climate (5), but the dynamics of this large reservoir have not been explored in situ. The world’s soils are modeled to warm 4.5 ± 1.1°C down to 1 m by the end of the 21st century under Representative Concentration Pathway (RCP) 8.5 (6), yet we lack data from deep soil warming field experiments to test assumptions of how deep soils will respond. This lack of understanding has contributed to large uncertainties in model predictions of soil carbon responses to warming (1, 7). In initial results from the first whole-soil warming experiment in mineral soil, warming stimulated soil respiration (8). However, increased heterotrophic activity does not necessarily result in net soil carbon loss, for example, if it is offset by increased plant inputs (9). Evidence that subsoil carbon stocks can be lost due to warming would imply larger positive climate-carbon cycle feedbacks than currently estimated. A deeper understanding and quantification of subsoil carbon dynamics would therefore be needed to accurately predict and plan for long-term changes (7, 10).

Here, we show how 5 years of whole-soil warming has significantly changed subsoil carbon stocks, composition, and respiration rates in a replicated field experiment (8). The study site is a mixed-conifer temperate forest on well-developed Alfisol soil (11). Mean annual temperature is 12.5°C, mean annual precipitation is 166 cm, and the site is characterized by warm, dry summers with little rainfall between May and October. Our experimental design is described in detail by Pries et al. (8). Briefly, 2.4-m-long vertical heating cables were installed around three, 3-m-diameter plots to warm the soil evenly to 1 m depth. Shallow heater cables (5 cm deep) installed in concentric rings at 1 and 2 m in diameter within the plots help to maintain the temperature difference near the surface. Each warmed plot is paired with a control plot. This design resulted in a 4.0°C temperature difference from 0.2 to 1 m in depth and a 2.6°C temperature difference above 0.2 m (fig. S1).


In a field study, we show that sustained whole-soil warming led to a loss of subsoil SOC stocks (Fig. 1A), along with similar losses of soil nitrogen (fig. S2 and table S1). The subsoil contained 45% of the total 18.6 ± 1.8 kg C m−2 (mean ± SE, here and throughout text) in the top 90 cm at this site (fig. S2 and table S1). After 4.5 years of treatment, the warmed plots had lost 3.21 ± 1.0 kg C m−2 in the subsoil compared to their paired control plots, a deficit of 33 ± 11% in SOC stored between 20 and 90 cm (P1:35 = 0.0026; Fig. 1A) (12). In contrast, there was a statistically nonsignificant increase in SOC and soil nitrogen stocks in the top 20 cm, although at 0 to 10 cm, the mean SOC stocks were higher in warmed (8.24 ± 1.7 kg C m−2) than control plots (5.00 kg ± 0.6 kg C m−2; P1:5 = 0.14; (Fig. 1A and fig S2). Because of the decline in SOC stocks with depth (fig. S2), the greatest absolute difference in SOC stocks was found at 20 to 40 cm, while the greatest relative difference was at 60 to 70 cm (Fig. 1A). Thus, warming had different effects on SOC stocks in the subsoil (net loss) and topsoil (nonstatistically significant gain).

Fig. 1 Soil profile bulk soil and respiration responses to warming.

(A) Mean and SE of the difference in SOC stocks between the warmed and control plots expressed in % [(Warmed − Control) / Control × 100], negative values highlighted in red, and positive values highlighted in blue. (B) Mean and SE of δ13C (‰) of the bulk soil. Soils collected in April 2018, 4.5 years after the onset of soil warming (n = 6). (C) Average of 5 years of mean annual apparent Q10 (Qa10) of CO2 production from the soil profile, positive warming effect highlighted in blue. Error bars are SE.

In addition to the impact on SOC stocks, warming also led to changes in SOC chemical and physical composition. The δ13C of SOC is an indicator of the degree of organic matter degradation, since microorganisms preferentially respire 12C during decomposition (13, 14). In the control plots, δ13C increased with depth down to about 50 cm, below which it remained relatively constant (depth, P1:51 < 0.0001; Fig. 1B). In the warmed plots, the δ13C profile was similar to the control plots down to 60 cm but became more enriched than the control plots below that (warming, P1:23 < 0.0001; depth × warming, P1:51 = 0.0004). The enrichment of 13C in SOC with soil depth is typical for soil profiles (5), but the greater 13C enrichment in warmed soils relative to control soils below 60 cm indicates a shift toward a more degraded SOC (15).

One major unknown regarding the vulnerability of soil organic matter is how molecular composition (16) and physical stabilization mechanisms (17) are affected by warming. Soil organic matter is heterogeneous, spanning partially decomposed plant material to microbial residues, and can be physically protected from microbial decomposition inside soil aggregates or via adsorption to reactive mineral surfaces (18). To determine how warming affected soil organic matter pools of differing physical composition, we fractionated the soil into free particulate organic matter (fPOM), which resembles plant material and is generally accessible to microbial enzymes; aggregate-occluded particulate organic matter (oPOM), which is protected within soil aggregates; and mineral-associated organic matter (MAOM), which contains more microbially processed organic matter that is protected within microaggregates or via association with mineral surfaces (19). As is typical, fPOM was mainly plant-derived (fPOM C:N ratio, 61.5 ± 2.2; plant litter layer, CN 65.3 ± 5.4), with a modern radiocarbon signature indicating decadal turnover times (8), while MAOM had a radiocarbon signature indicative of older organic matter (8), with C:N ratios signifying microbial origin (MAOM CN ratio, 16.1 ± 0.69) (20).

The unprotected, plant-derived soil organic matter (fPOM) was the most sensitive to warming. In control soils, the distribution of bulk soil carbon among the three soil organic matter fractions did not change with depth (Fig. 2) (19). In contrast, after 3.5 years of warming, there was more fPOM in shallow soils and a large loss of fPOM at depth (Fig. 2A and fig. S3A) such that the proportion of total carbon in the fPOM was slightly higher at the surface (P1:2 = 0.074) and was significantly lower at depth (fPOM, P1:2 = 0.027; Fig. 2A), resulting from more shallow-fPOM carbon in the warmed (727 ± 103 g C m−2) than in the control (317 ± 86 g C m−2; fig. S3A) plots. Thus, warming had a different effect on decomposition of different types of organic matter across the profile (warming × depth interaction for fPOM g C g−1 soil C, P2:10 = 0.021). The role of fPOM decomposition in driving the net loss of SOC stocks in the subsoil was confirmed by deep soil 13C-enrichment (Fig. 1C), since fPOM was more 13C-depleted (δ13C = −25.9 ± 0.46%) than MAOM (δ13C = −22.3 ± 0.34%) below 60 cm (fig. S4).

Fig. 2 Distribution of soil carbon among soil organic matter fractions.

(A) fPOM, (B) oPOM, and (C) MAOM, (mean ± SE, n = 3) sampled in June 2017, 3.5 years after the onset of warming. P are for any significant effects of treatment, depth, and their interaction (Trt:D). The asterisk indicates a significant treatment effect (P1:2 = 0.02) at 85 to 89 cm in the fPOM fraction after within-depth, post hoc analysis. Stocks of carbon in each fraction are shown in fig. S3.

Warming-induced losses of SOC stocks were accompanied by 30 ± 4% greater CO2 effluxes over all 5 years, with no significant attenuation over time (warming, P1:130 < 0.0001; year, not significant; Fig. 3). Surface CO2 effluxes were measured in seven locations per plot monthly (12). This increased CO2 efflux rate equates to a mean apparent Q10 [Q10a (12)] of 2.3 ± 0.2, which emerges from many underlying microbial and abiotic (e.g., mineral surface interactions) processes. Sustained enhancement of soil CO2 efflux due to warming over 5 years, also reported in shallow soil warming experiments (21, 22), is evidence that the loss in subsoil carbon stocks due to warming was associated with an increase in net belowground CO2 emissions to the atmosphere. If these elevated emissions continue with climate change, then they would create a positive feedback to climate change.

Fig. 3 Mean annual soil surface efflux.

Warmed and control plots (circles), and the effect of warming on mean annual soil respiration rates per year in % difference between warmed and control plots (squares). n = 24 to 36 per year (10). Error bars are SE.

Seasonal trends in soil efflux were not obvious in this Mediterranean climate zone (fig. S5) due to the positive relationship with soil temperature (P1:167 < 0.0001 at 15 cm) and nonlinear, unimodal relationship with soil moisture (P1:167 = 0.0296 at 10 cm), with an apparent optimum volumetric water content of 0.22 cm3 cm−3. In contrast, the warming effect on soil efflux did vary seasonally (fig. S6). The warming effect was amplified by soil moisture (P1:116 = 0.0011, slope = 5.0 × 10−4, r2 = 0.08) and diminished at higher temperatures (P1:116 < 0.0001, slope = −0.05, r2 = 0.16), resulting in the largest warming effects during the spring and lowest during the dry summer and fall seasons (fig. s6). Thus, while warming had a consistent effect on mean annual soil efflux (Fig. 3), the distinct dry season response is indicative of moisture limitation to the warming impact. This moisture limitation could help to explain the greater CO2 efflux response to warming observed in more mesic forests than in shrublands and grasslands (8, 9).

Soil respiration at all depths was responsive to warming (Fig. 1C), although absolute CO2 production (m−2) was much lower in deeper than shallower soils (depth, P1:1163 < 0.0001; Fig. 4). The mean warming effect on CO2 production from the whole profile and all sampling dates (n = 536) was 42.6 ± 4.5%, more than was estimated by surface efflux. To quantify the effect of warming at each sampling depth, we calculated Q10a using the observed temperature difference between each paired plot at each depth. The average Q10a was 3.47 ± 0.23, higher than is prescribed by most models for decomposition, which varied with depth, and was highest (4.5 ± 0.67) below 70 cm (Fig. 1C). Thus, the responsiveness documented after 2 years of warming, in which all depths had Q10a > 2 (8), was sustained over 5 years (Fig. 1C). There was no significant change in Q10a over time (P > 0.05) for either soil profile CO2 (fig. S7) production or surface efflux (Fig. 3), demonstrating that the impact of warming on soil respiration has been sustained throughout the whole-soil profile for 5 years.

Fig. 4 CO2 production throughout the soil profile.

Average of 5 years of mean annual production inset shows deeper depths on finer scale x axis. Error bars are SE.


There are at least four possible explanations for the large loss of SOC and nitrogen stocks in the subsoil. First, microbial decomposition of soil organic matter was increased by warming in the subsoil as evidenced by increased CO2 production at depth. Second, warming led to increased concentrations of dissolved organic carbon (DOC) (P1:165 < 0.0001) in soil pore water after rain and snowmelt events at our site (fig. S8). The transport of this DOC to deeper depths could have alleviated carbon limitation (23) and accelerated decomposition of deep soil organic matter stocks, since it is hypothesized that decomposition of deep soil carbon, including fPOM, is inhibited by a lack of labile compounds to prime microbial activity and enzyme production (5, 24). Third, the proportion of microbial biomass to SOC (Cmic:Corg) was 2.5 times greater in the subsoil (2.5 ± 0.55%) than the surface soil (1.1 ± 0.15%) at our site (fig. S9). Since soil microbial activity is intrinsically temperature sensitive (25, 26), the greater Cmic:Corg in deep soil may have made it more likely for the stimulation of microbial activity by warming to lead to a measurable effect on SOC stocks. For example, a 30% increase in microbial activity in the subsoil could have 2.5 times the impact on SOC in the subsoil than in the topsoil based on the Cmic:Corg ratios. Last, the decomposition of more chemically complex (higher activation energy) subsoil carbon has been hypothesized to be more temperature sensitive (16), which, if true, would explain the losses of older subsoil SOC; however, we did not observe higher temperature sensitivity in respiration [this study and (8)].

Although warming led to higher DOC concentrations in soil pore-water samples, we did not find a significant warming effect on dissolved nitrogen concentrations (P1:177 > 0.05; fig. S10). Nearly all of the dissolved nitrogen was organic, rather than inorganic, indicating a high demand for inorganic nitrogen in this system. At this site, 80% of fine root biomass are concentrated in the top 30 cm of the soil (8). Dissolved nitrogen concentrations increased between 20 and 70 cm (P1:16 = 0.0049), attributed to high nitrogen demand by fine roots. The lack of an effect of warming on dissolved nitrogen, despite the effect on DOC concentrations, suggests that excess mineralized nitrogen due to warming was rapidly assimilated by plants in this forest (22).

Surface SOC and nitrogen stocks had a different response to warming for several possible reasons. First, we achieved a lower level of warming near the surface, so expected less impact (fig. S1A). Second, warming could have increased surface inputs of fine roots (27), which are concentrated toward the soil surface (19). Third, the warmed surface soils were slightly drier during summer (fig. S1B), which could have inhibited the decomposition of plant litter and fPOM and decreased transport from the surface to deeper horizons. Surface drying is expected to coincide with warming, except where altered precipitation occurs, and their interactive effects merit further research. Last, changes in the microbial community composition with depth could imply different temperature sensitivities in surface and subsoils. The average increase in SOC stock in the top 10 cm, although not statistically significant, is on the higher end of the large variation in soil warming responses found at sites with a similar initial SOC stock in the top 10 cm (2).

A critical uncertainty in climate-carbon feedbacks is whether deep SOC is as vulnerable as faster cycling shallow SOC to perturbations such as altered climate (18). Recent incubation experiments have reported both higher (8) and lower (27) temperature sensitivity of subsoils compared to surface soils. However, incubations introduce many artifacts (28) including disruption of soil structure, loss of plant-soil interactions, and loss of diurnal climate variability, which differentiate shallow and deep soil environments (29). Thus, laboratory incubations alone are not sufficient for quantifying warming effects. In a field experiment, we showed that warming can lead to a net loss of SOC stocks below 20 cm, mainly due to the decomposition of unprotected particulate organic matter, and to comparable increases in CO2 effluxes from the whole-soil profile. Subsoil CO2 losses were stimulated by warming without a compensating increase in inputs. Although the results from our in situ experiment were limited to three replicated paired plots and two soil cores from each plot, we found statistically significant effects of warming on soil carbon and nitrogen stocks. Apparently, divergent responses of subsoil and surface-soil carbon to warming highlight the need for more experimental data from different ecosystems and explicit vertical resolution in models, which only 3 of >20 Coupled Model Intercomparison Project phase 6 (CMIP6) models currently have (28). Without the representation of heterogeneous SOM pools and subsoil environments, models will likely fail to predict the long-term response of soil carbon to warming.

Soils contain an estimated 1400 Pg C in the top 1 m, over half of which is stored in the subsoil (4). Temperate forest soils in the same soil order as our study site (Alfisols) store 8% of global soil carbon (29) and are predicted to warm by 4.1° ± 1.0°C by the end of the 21st century under RCP 8.5 (7). To explore the order-of-magnitude importance, we scaled our findings (33% of subsoil carbon lost in 5 years) to all Alfisols and estimate that the global loss would be 18 Pg C from Alfisols alone, which is not currently accounted for in predictions of future forest carbon sinks (3). If a similar pattern holds in other ecosystems, then roughly 231 Pg C could be lost from subsoils, with a strong amplification effect on climate change that should not be ignored in predictions of future climate scenarios. Results from whole-soil warming experiments in different ecosystems are needed to provide more confidence in these extrapolative exercises. These extrapolations do not take into account predicted responses to other changing conditions, such as rising atmospheric CO2 concentrations, which were not tested in this experiment but will likely affect forest carbon and nutrient cycling over the 21st century, along with warming (30).

Understanding subsoil carbon dynamics is important not only for climate predictions but also for terrestrial carbon sequestration. Many land-based sequestration strategies count on the stability and longevity of deeper soil carbon (31). We have shown here that subsoil carbon, including recent plant inputs, is vulnerable to warming, implying that the longevity of sequestered carbon depends on future climate and that terrestrial sequestration will be more effective if it is accompanied by early action to mitigate anthropogenic emissions. Thus, understanding carbon dynamics, including subsoil vulnerability, is critical for intentional sequestration, the terrestrial carbon sink, and climate prediction over the coming decades to centuries.


Study site

The University of California Blodgett Experimental Forest is located in the Sierra Nevada foothills of California at 1370 m above sea level (120°39′40′′W; 38°54′43′′N). It is a mixed coniferous forest of mainly ponderosa pine (Pinus ponderosa), sugar pine (Pinus lamberiana), incense cedar (Calodefrus decurrens), white fir (Abies concolor), and douglas fir (Pseudotsuga menziesii), located below the permanent winter snowline. Mean annual precipitation is 166 cm, with most of it occurring from November to April. The mean annual temperature is approximately 12.5°C. The soils are classified generally as Alfisols and are described as Holland series, fine-loamy, mixed, superactive, mesic Ultic Haploxeralfs of granitic origin with thick, >5-cm O horizons (8).

Experimental design

The soil warming experiment, described in detail by Pries et al. (8), began in October 2013 when heaters were turned on to warm the soil 4°C above ambient to 1 m depth while maintaining the natural temperature gradient following the design of Hanson et al. (32). Briefly, the experiment consists of three blocks of replicated, warmed-control paired plots, and each plot is 3 m in diameter. Warming is maintained throughout the plot and down to 1 m depth with 22, 2.4-m-long resistance heating cables (BriskHeat, Ohio, USA) installed vertically in 1″ diameter metal conduit spaced 50 cm apart and 0.25 m outside the edge of the plots. Two concentric rings of heater cable at 1 and 2 m in diameter and 5 cm in depth were installed within the plots to compensate for surface heat loss. Similarly, we installed vertical conduit around and two concentric rings of unheated cables in the control plots.

We used temperature sensors (Omega 44005) to monitor soil temperature every 10 s at 5, 15, 20, 30, 50, 70, 75, and 100 cm in all plots. We used enviroSCAN (Sentek, Australia) probes fitted with capacitance sensors at 10, 30, 50, and 90 cm in each plot, which were calibrated five times in both dry and wet periods to monitor soil volumetric water content. We used dataloggers (CR1000, Campbell Scientific, Utah, USA) to record soil temperature and moisture values continuously and created a program that adjusted the power to the heaters every 10 min based on soil temperatures in the control plots to maintain a targeted 4°C warming in the paired heated plots. Temperature sensors at 15 and 20 cm were queried to control the power supply to the shallow heaters, while temperature sensors at 70 and 75 cm were queried to control power supply to deep heaters.

Each plot contains five, ¼′′ diameter stainless steel, septa-topped, gas sampling tubes permanently installed at a 45° angle. These are used to sample soil gas from 15, 30, 50, 70, and 90 cm and to model CO2 production from between these depth increments within the soil profile. Each plot is also equipped with two porous quartz tipped sampling tubes with an average pore size of 2 μm (Prenart Super Quartz Sampler, Prenart equipment, Frederiksberg, Denmark) used to collect soil pore water samples to measure dissolved organic matter concentrations. The pore water samplers were initially coated with silica flower and installed at a 30° angle, one at 20 cm and one at 70 cm, to collect shallow and deep soil pore water samples via vacuum suction. We collected soil pore water samples into acid washed vials when possible after rain and snow melt events in the winter and spring (five times in winter 2013–2014, three times in winter 2014–2015, eight times in 2015–2016, zero times in 2016–2017, and five times in 2017–2018).

Bulk soil analysis

We collected two soil cores (4.8 cm in diameter) divided into approximately 10-cm-depth increments down to 90 cm from each plot on 4 April 2018. The first core was stored at 4°C until it could be sieved, which was completed within 1 week of sampling. We sieved the samples through a 2-mm mesh to separate soil (<2 mm) from roots and stones (>2 mm). Roots longer than 1 cm that passed through the 2-mm sieve were also picked from the soil with forceps. We weighed a subsample of the moist soil and after drying at 105°C to calculate gravimetric water content. A subsample of the soil was ground with stainless steel balls on a roller mill for 48 hours and then analyzed by Elemental Analyzer-Isotope Ratio Mass Spectrometry (EA-IRMS; Isoprime 100 IRMS in line with a Vario micro cube EA, Isoprime, UK) for carbon concentration and δ13C. The second core was freeze-dried and then sieved to 2 mm. A subsample from each depth increment of the second core was ground in a ball mill (MM400, Retsch, Haan, Germany) and analyzed by EA-IRMS (Flash, 2000-HT Plus, linked by Conflo IV to Delta V Plus IRMS) for carbon concentration and δ13C. Results from both cores are presented as mean values for each plot.

Bulk density measurements were taken from three soil pits adjacent to the experimental plots in 2014. Each pit was 1 m3, and horizontal soil cores with a diameter of 6 cm were carefully excavated to measure soil bulk density. Total bulk density was used for calculation of diffusivity for CO2 production measurements, but bulk density of the soil <2 mm after sieving was used for calculations of SOC stocks. We used the mean bulk density from the three pit measurements as one estimate for bulk density for all plots in our experiment. We calculated SOC stocks by multiplying this single bulk density value (g soil cm−3) by carbon concentration (g C g soil−1) from each soil core. For calculation of the SE bars shown in Fig. 1, we first calculated the difference in SOC stocks between warmed and control cores for each core pair within each plot pair. An average value for each plot pair was then used to calculate SD and SE (n = 3) from each plot. We could not take accurate bulk density measurements from our long-term warming experiment plots themselves, because digging a pit large enough to take horizontal cores would be too much of a disturbance to these long-term plots. Bulk density is typically inversely related to soil organic matter content. Therefore, it is possible that the loss of SOC stocks due to warming reported here is slightly underestimated if warming also led to increased bulk density. Total SOC stocks below 20 cm were calculated by summing the SOC stocks measured in 10-cm increments from 20 to 90 cm for each plot.

Soil organic matter fractionation

In June 2017, one soil core (4.4 cm in diameter) from each field plot was collected for bulk soil analysis and soil organic matter physical fractionation. Bulk soil samples from three depths—10 to 14, 44 to 49, and 85 to 89 cm—were sieved to 2 mm to separate bulk soil from roots and stones, following which bulk soil and roots were freeze dried. A 20-g sample of freeze dried soil was then separated into fPOM (free light fraction), aggregate oPOM (occluded light fraction), and MAOM (dense fraction) using a sodium polytungstate solution (1.7 g ml−1) and sonication (19). Mass recovery after fractionation was ±2%. A bulk soil sample and each soil organic matter fraction were ground and analyzed for carbon concentration, δ13C, and δ15N by EA-IRMS (Isoprime 100 IRMS in line with a Vario micro cube EA, Isoprime, UK).

Soil CO2 efflux

We measured surface soil respiration (efflux) approximately monthly from March 2014 to December 2018 using an infrared gas analyzer (Licor LI-6400, Licor, Lincoln, Nebraska, USA) or a cross-calibrated cavity ring down spectrometer (Picarro, Santa Clara, California, USA), fitted with a 10-cm-diameter opaque chamber on seven permanently installed PVC collars in each plot. There were 8 sampling dates in 2014, 11 in 2015, 8 in 2016, 9 in 2017, and 10 in 2018.

Soil profile CO2 production

On the same days as the surface soil respiration measurements, we collected soil gas samples from the gas sampling tubes with a syringe and stored them in evacuated glass vials. We sampled soil gas from 6.35-mm-diameter stainless steel tubes inserted to 15, 30, 50, 70, and 90 cm within the soil profile and sealed at the surface with swage pipefittings (Swagelok Ohio, USA) and septa. We measured CO2 concentration on the soil gas samples by syringe injection via gas chromatography (Shimadzu 2014, Shimadzu Scientific Instruments Inc., USA).

We modeled depth-resolved CO2 production using the flux gradient approach (8, 33, 34) from soil CO2 concentrations collected approximately monthly. The CO2 production for the volume of soil between the shallowest soil gas sample (15 cm) and the soil surface was calculated using the flux estimate at the shallowest soil gas sample depth and the surface flux measurement, which was measured on the same day (33, 34). The mean annual CO2 production warming response from each plot was calculated as the mean warming effect of all sampling dates for each plot per year.

DOC and dissolved nitrogen

Soil pore water samples were stored frozen until analysis. They were diluted with Nanopure water and analyzed for DOC and total dissolved nitrogen concentrations (Shimadzu TOC-L, Shimadzu Scientific Instruments Inc., USA). Dissolved inorganic carbon was measured on the samples collected in 2014 (Shimadzu TOC-VCPH analyzer, Shimadzu Corporation, Japan) and was found to comprise <1% of the total dissolved carbon, so we report the total dissolved carbon as DOC concentrations. Inorganic nitrogen (NO3, NO2, and NH4+) were measured on an Alpkem Flow Solution IV Automated wet chemistry system (O.I. Analytical, College Station, TX), and dissolved organic nitrogen was calculated as the difference between total and inorganic nitrogen. Inorganic nitrogen concentrations were often below detection limit (0.5 mg N/liter).

Microbial biomass

We took two 10-g subsamples of the field-moist soil from 0- to 10-cm-, 10-to 20-cm-, 30- to 40-cm-, 50- to 60-cm-, and 80- to 90-cm-depth increments from the first core of the April 2018 sampling to measure microbial biomass carbon using chloroform fumigation extraction [modified from (35)]. All extractions were completed within 1 week of field sampling, and samples were stored in air tight bags at 4°C until extraction. Briefly, one subsample was extracted by shaking in 0.05 M K2SO4 for 1 hour and filtering over a Whatman #1, ash-free filter. The second subsample was first fumigated with chloroform under vacuum for 5 days and then extracted in the same way. The extracts were analyzed for total organic carbon on a Shimadzu TOC-L (Shimadzu Scientific Instruments Inc., USA). Microbial biomass carbon was calculated as the difference between the chloroform-fumigated and nonfumigated concentrations on a dry-soil basis. We used a 0.45 extraction efficiency correction factor (36) to report microbial biomass values at comparable scales to those reported by Xu et al. (37).

Data analysis

We expressed the warming response as the percentage difference between warmed and control, normalized by the control value for soil carbon stocks and surface fluxesWarming response (%)=WarmControlControl×100

The warming response was always calculated for each plot pair or for each depth within each plot pair when relevant. We also compared warmed and control soil surface efflux and CO2 production from the soil profile using an apparent Q10 (Q10a)Q10a=(RWRC)10(TWTC)

Where RW and RC are the soil efflux or CO2 production from the warmed or control plot, respectively, and TW and TC are the temperatures from the warmed and control plot, respectively. Q10a was calculated for each plot pair at each sampling date and depth (when applicable). For surface efflux, soil temperature at 10 cm was used, since this represents the soil depth with the greatest contribution to overall CO2 production (Fig. 4). This equation accounts for the exact TW and TC at any given depth and time point. This is an “apparent Q10 value, because it is not associated with a specific biological response to temperature but instead is a comparison of the emergent CO2 efflux or production, which derives from many processes, which may be changing over time due to the treatment.

Soil surface efflux and CO2 production varied seasonally. To test whether the warming response and Q10a were sustained over 5 years, we first calculated the mean annual warming response or Q10a based on the mean of the soil efflux values for each plot at each sampling date. We then calculated the mean annual Q10a for each of the 5 years of the experiment. We then compared the mean annual efflux or production value for each of the three blocks over 5 years (n = 15).

In general, we ran statistical analyses by building mixed effects models using the nlme package (38) in RStudio (39) using restricted maximum likelihood with block (n = 3) as a random effect to test the effect of warming treatment on various measured responses. For all models, we visually checked for normality of distribution and homoscedasticity using Q-Q plots and histograms of residuals and applied a log transformation to meet model assumptions when needed, such as for SOC stocks and dissolved organic matter concentrations.

The bulk soil SOC stock, bulk soil δ13C, soil organic matter fractions, surface soil efflux rates, and soil profile CO2 production datasets each required slightly different statistical model structures due to their different data structures. For SOC stocks and δ13C, our model included treatment, depth and their interaction as fixed effects, and soil core nested within depth and block as random effects. The model was first run for all 10-cm-depth increments from 0 to 90 cm. Since a significant depth × treatment interaction was found for SOC stocks (table S2), we ran a post hoc analysis separating the profile into the topsoil (0 to 20 cm) and subsoil (20 to 90 cm). Total stocks in the topsoil and subsoil were calculated by summing up the SOC stocks over the entire depth interval. For each soil organic matter fraction, our model included treatment, depth and their interaction as fixed effects, and block as a random effect. When a significant treatment × depth interaction was found, as for the fPOM, we ran post hoc analyses within each depth increment.

For soil surface CO2 efflux rates, we took the mean of the seven CO2 flux measurements from each plot as the soil flux value for each plot at each sampling date. To test the treatment and sampling date effects on soil surface CO2 flux rates, our model included treatment, sampling date and their interaction as fixed effects, and block nested within date as random effects to account for repeated measures. To test the effect of soil moisture (10 cm) and temperature (15 cm) on control CO2 fluxes (not from the warmed plots), our model included volumetric water content, temperature and their interaction as fixed effects, and block nested within sampling date as random effects. To test whether warming affected the mean annual flux response, we first calculated the warming response for each block at each sampling date. We then modeled the fixed effect of year on the mean annual warming response, with block nested within sampling date as random effects to account for repeated measures.

We tested the fixed effects of treatment, depth, and their interaction on soil CO2 production from within the soil profile using block nested within sampling date as random effects. As we did with surface fluxes, we tested whether warming affected the mean annual CO2 production from the profile over 5 years by first calculating the warming response for each depth, block, and sampling date and then applying a mixed effects model with year, depth and their interaction as fixed effects, and block nested within sampling date as random effects to account for repeated measures. For DOC concentrations in soil pore water, our model included treatment, depth, year and their interactions as fixed effects, and depth nested within block nested within sampling date as random effects.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We wish to acknowledge the contributions of the Lawrence Berkeley National Laboratory’s Geosciences Measurement Facility and its support of custom instrumentation, as well as the University of California, Berkeley Blodgett Forest Field Research Station staff for their excellent support, especially A. Thompson, A. Mason, and J. York. Funding: This work is supported by the U.S. Department of Energy Office of Science, Office of Biological and Environmental Research Terrestrial Ecosystem Science Program under award DE-SC-0001234, Swiss National Science Foundation project 200021_172744, and Colorado State University. Author contributions: Conceptualization and methodology: J.S., M.S.T., and C.E.H.P. Investigation: J.S., C.E.H.P, C.C., R.C.P., N.O., and M.W.I.S. Formal analysis: J.S. and W.J.R. Writing (original draft): J.S. Writing (review and editing): All authors. Funding acquisition: M.S.T., W.J.R., and M.W.I.S. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Published data will be available on the ESS-DIVE repository Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article