Bioenergetic control of soil carbon dynamics across depth

Further details about radiocarbon and thermal analysis, isotopic partitioning procedures and quantification of their uncertainty, and statistical analyses can be found in Supplementary Methods.

Study soils, experimental design and soil sampling

We selected three soil types: eutric cambisol, chromic vertisol and silandic andosol⁷⁰. The three soil profiles studied were found in long-term semi-natural grasslands located relatively close to each other (<100 km) in Auvergne, France. They developed under a similar temperate semi-continental climate and mainly differed by their parent materials: granite, basalt, and trachyandesite for the cambisol, vertisol and andosol, respectively (Supplementary Table 7). The three soil profiles studied were well-aerated since none showed reductimorphic layers, and redoximorphic features were present only in the vertisol.

The experiment had a factorial design of two crossed factors: two soil layers including topsoil versus subsoil, and three soil types for a total of six treatments each including four replicates. We collected 20 cm high soil cores from the two layers for each soil profile. Intact soil columns of 8 cm diameter were extracted using a percussion core drill that can be opened from sideways. Topsoil cores were taken in the 5–25 cm depth part of the A horizon, which allowed us to remove both the native vegetation and a large proportion of their fresh litter. Subsoil cores were taken in the top 20 cm of the B horizon, that are at 40–60, 55–75 and 35–55 cm depth respectively for the cambisol, vertisol and andosol (Supplementary Fig. 5).

Soil biogeochemical analyses

Four soil cores of each treatment were used to characterise soil biogeochemical properties. These soil cores were first sieved at 2 mm and a portion was air-dried at room temperature. A portion of the air-dried soil was also ground (<250 μm) to homogeneity.

Soil C and N concentrations and δ¹³C was determined on ground soil containing ~1 mg C using an elemental analyser (EA, Carlo Erba, Rodana, Italy) coupled to an isotope-ratio mass spectrometer (IRMS; Elementar, Langenselbold, Hesse, Germany). None of the soil profiles contained carbonates, and total soil C content is interpreted as SOC only here. To determine the relative contribution of mineral-associated organic matter versus particulate organic matter to SOC, each soil sample was fractionated by particle size (50 μm) after full soil dispersion⁷¹. Briefly, 5 g of air-dried soil was shaken for 18 h in sodium hexametaphosphate (0.5%) with beads to completely disperse the soil. The dispersed soil was then rinsed onto a 50 μm sieve and the fraction passing through (<50 μm) was collected as MAOM, while the fraction remaining on the sieve was collected as POM. Each fraction was then analysed for C concentration using an EA.

The radiocarbon signature (Δ¹⁴C) of SOC was determined by Accelerator Mass Spectroscopy (AMS) on ground soil containing respectively ~1 and 0.14 mg C for topsoil and subsoil with a Mini Carbon Dating System (ECHoMICADAS) operated at LSCE (Climate and Environment Sciences Laboratory) in Gif-sur-Yvette, France. To assess SOC dynamics, we estimated SOC turnover time based on radiocarbon measurements using a modelling approach^72,73. The following time-dependent, homogeneous one-pool model was used:

where t is the time (in year), F^14C_SOC is the ¹⁴C content of SOC, F^14C_atm is the ¹⁴C content of CO₂ in the local atmosphere, τ is the mean SOC turnover time and λ is the radioactive decomposition constant for ¹⁴C (1.21 × 10⁻⁴ year⁻¹). We ran the model from 50 kyr BP until the year 2016 to calculate the predicted SOC Δ¹⁴C at the year of sampling for a range of τ values (1 to 30,000 years). We then derived τ values from our Δ¹⁴C measurements for each sample based on the relationship between τ and predicted Δ¹⁴C (Supplementary Fig. 6). Though the model assumption of SOC as a homogeneous one-pool is clearly an oversimplification, our approach remains useful to compare SOC dynamics across soil samples⁹.

Thermal analyses were used to characterise the bioenergetic signature of SOM. The activation energy (E_a) of SOC thermal decomposition was measured by evolved gas analysis during ramped combustion using Rock-Eval® thermal analysis. Sequential ramping by pyrolysis and oxidation was performed on ~60 mg of ground soil using a Rock-Eval® 6 Turbo device (Vinci Technologies, France). Hydrocarbon effluents were quantified by flame ionisation detection during the pyrolysis phase, while CO and CO₂ were quantified by infra-red detection during both ramping phases. Based on evolved SOC kinetics, a regularised inverse method was used to determine the continuous distribution of E_a that best predicts the profile of SOC decay measured during ramped combustion⁷⁴ (Supplementary Fig. 1a, b). The distribution of E_a was then integrated to calculate the mean (µE_a) and standard deviation (σE_a) of activation energy. Additionally, the temperature at which 90% of hydrocarbons was evolved during pyrolysis (T₉₀-H_xC_y-pyrolysis, °C), and the temperatures at which 50% of CO₂ was evolved during pyrolysis and oxidation (T₅₀-CO₂-pyrolysis and T₅₀-CO₂-oxidation, °C) were used as indices of SOC thermal stability⁷⁵. The hydrogen index (HI) and oxygen index (OI) providing information about elemental H:C and O:C ratios of SOM were calculated respectively as the amount of hydrocarbons and the amount of CO₂ and CO formed during pyrolysis divided by SOC concentration.

Energy density of SOC (ΔE) was measured by differential scanning calorimetry (DSC)^28,76. Oxidation ramping was performed on ~60 mg of ground soil using a DSC thermal analyser (TGA-DSC 3+ model, Mettler-Toledo, Greifensee, Switzerland) to measure the net energy released by SOM combustion (enthalpy of combustion), knowing that some of the energy applied to the sample is consumed by the breakdown of the organo-mineral associations. The ΔE was calculated as the net energy released determined by integration of heat flux over the exothermic region associated to SOC combustion (185–600 °C, Supplementary Fig. 1c), divided by SOC concentration and multiplied by SOC molar mass (Supplementary Table 8, Supplementary Methods). Calorimetry also allows to estimate the degree of reduction of SOC (γ_SOC)⁷⁷, which was calculated here as ΔE divided by Q₀, the oxycaloric quotient representing the ratio between the enthalpy of combustion and the degree of reduction⁷⁶. We used a Q₀ value of 109.04 kJ mol⁻¹ SOC degree of reduction⁻¹, obtained from the average of the heat of combustion of a large set of organic compounds⁷⁸. Additionally, a return-on-energy-investment (ROI) parameter was calculated as ΔE divided by µE_a²⁶ and was used here as an index of SOC energy quality^27,28.

Soil microbial biomass was measured on fresh sieved soil by the chloroform-fumigation-extraction method⁷⁹. We performed an extraction with 50 mL 0.5 M K₂SO₄ on 10 g of fresh soil. A second set of samples was placed in a vacuum desiccator and fumigated with chloroform for 24 h prior to K₂SO₄ extraction as above. Extractable C was analysed using an automated analyser (TOC-L analyser, Shimadzu, Milton Keynes, UK). Microbial biomass C was calculated as the difference between the fumigated and unfumigated C extracts and a correction for extraction efficiency was applied by dividing with a coefficient of 0.45⁸⁰. Initial soil mineral N content was measured from 25 g of fresh soil after extraction in 2 M KCl using a continuous-flow analyser (AA3, Bran + Luebbe, Norderstedt, Germany). All roots retained by sieving of soil cores and visible roots in sieved fresh soil were handpicked and washed with tap water. After drying at 60 °C for 48 h, roots were weighted and root density was calculated as root dry mass divided by the soil core volume.

We also characterized the mineral composition of soil samples. Soil pH_water was measured in a 1:5 soil:solution ratio after 1-h end-over-end shaking. Soil clay content was measured using the pipette method⁸¹. Phyllosilicate mineralogy was determined by X-ray diffraction⁸². Random oriented powders with a Philips PW 3710 X-ray diffractometer with Cu-Kα radiation at 40 kV and 40 mA were used to obtained spectra (Supplementary Fig. 7). A counting time of 13 s per 0.02° step was used for 2θ in the range 3.5–80°. Cation exchange capacity (CEC) as well as major exchangeable cations (Ca, Mg, Na, K, Al, Mn, Fe) were determined determined using the cobalt hexamine exchange method⁸³. The concentrations of divalent cations involved in organo-mineral cation bridging was calculated as the sum of exchangeable calcium and magnesium (Ca_ex+Mg_ex). Pedogenic reactive metals were quantify by selective dissolution procedures of Fe, Al and Si using standard methods of citrate–dithionite (d), acid ammonium oxalate (o), and sodium pyrophosphate (p) extractions⁸⁴. Crystalline Fe minerals was then calculated as the difference between Fe_d and Fe_o (Fe_d–o). We used Fe_o to quantify short-range-ordered Fe-oxyhydroxides such as ferrihydrite, while Al_o+Si_o were used to quantify organo-metal complexes and short-range-ordered aluminosilicates such as allophane. Al_p allows to quantify organo-metal complexes, but pyrophosphate is not completely selective and can also extract Al from silicates, which could be indicated by Si_p⁸⁴. Organo-metal complexes was thus calculated as Al_p–xSi_p, where x is the Al:Si ratio of the dominant clay mineral (x = 0.5 for subsoil of the vertisol dominated by montmorillonite, and x = 1 for all other layers dominated by vermiculite, kaolinite or halloysite).

Plant isotopic labelling and soil incubation experiment

Soil cores collected for the experiment were immediately proceeded in the field following sampling to establish microcosms with a new soil column made of intact soil cores derived exclusively from either topsoil or subsoil (Supplementary Fig. 2). For each microcosm, three soil cores of the same layer were gently stacked vertically and tightly sealed together within a polyethylene sheath before transfer into a 60 cm high PVC pot (diameter 10 cm, height 60 cm, with a permeable bottom-cap to allow drainage). For each the six treatments including the three soil types and the two layers, we included four planted replicates and four unplanted controls for a total of 48 microcosms. The experiment started within 2 months after sampling and microcosm were stored at 4 °C until then.

The experiment was performed for 279 days, from late August 2016 until early July 2017. Two weeks before starting the incubation experiment, the microcosms were transferred at ambient temperature, irrigated until water saturation and weighed after 48 h of water percolation to measure the soil water-holding capacity (WHC). Planted microcosms were sown (1400 seeds m⁻²) with Dactylis glomerata, a fast-growing grass species with a dense and deep root system commonly found in temperate grassland⁸⁵. After germination, microcosms were transferred to a greenhouse exposed to natural light and temperature conditions (Clermont‐Ferrand, temperate semi‐continental climate). The greenhouse was coupled to a continuous dual-labelling (¹³C/¹⁴C) system^86,87. Labelled air depleted in both ¹³C and ¹⁴C was produced by injecting fossil fuel-derived CO₂ (δ¹³C = −35.23 ± 0.02 ‰, Δ¹⁴C~0‰) in CO₂-free air ([CO₂] < 20 ppm) up to reach ambient CO₂ concentration (400 ppm), and the greenhouse was continuously supplied with labelled air during daytime, with a flow renewing the greenhouse volume once every 2 min to maintain constant CO₂ concentration, δ¹³C and Δ¹⁴C⁸⁶. Soil water content was monitored daily using soil moisture sensors (ECH2O EC-5, Decagon®, USA) inserted at 5 cm and drip irrigation was adjusted individually for each treatment as to maintain moisture around 85 ± 5% of WHC. In order to compensate for the low nutrient availability and plant growth potential expected in subsoil relative to topsoil, we fertilised the planted subsoil microcosms for each soil type on day 51 after planting. Unplanted subsoil microcosms were kept unfertilised to avoid excessive concentrations of mineral nutrients, which already tend to accumulate in soils in absence of rhizodeposition and nutrient uptake by plant roots^41,88. The fertilisation solution was composed of inorganic N (NH₄NO₃), P, S, K and Mg, with a dose of 11.5, 0.6, 1.0, 1.5 and 0.7 g m⁻², respectively. The fertilization was adjusted to compensate for initial differences in measured soil mineral N concentrations compared with planted topsoil microcosms, with a single dose added for each soil type and a second dose including only N applied to the andosol (Supplementary Table 9). This allowed to reach similar soil mineral N concentrations, plant N concentration and biomass between planted topsoil and subsoil microcosms following fertilization (Supplementary Tables 9, 10).

For the first incubation series, we measured CO₂ fluxes of each microcosm 76, 139, 174, 201, 242 and 272 days after planting, from late fall until early summer (n = 288 incubations). Microcosms were sealed in opaque airtight chambers and incubated for 24 h at 21.5 °C. Chamber gas was sampled at the end of incubation, and its CO₂ concentration and δ¹³C were measured using a gas chromatograph (Clarus 480, Perkin Elmer, Waltham, MA, USA) and an isotopic analyser (G2201-i, Picarro, Santa Clara, CA, USA). The amount and δ¹³C of CO₂ derived from plant-soil respiration were corrected for background atmospheric CO₂.

For the second incubation series performed 279 days after planting, the soil column was gently extracted from each microcosm. Shoots were cut at the soil surface for planted microcosms, and the soil column was vertically sliced to recover the three original soil cores. Each soil core was transferred into a 3 L flask as gently as possible to maintain the soil core structure intact (n = 144 incubations). After a preincubation period of 24 h at 21.5 °C, flasks were airtight‐sealed and incubated for 24 h at 21.5 °C. The evolved CO₂ was trapped in NaOH and its concentration was measured using an automated analyser (TOC-L analyser, Shimadzu, Milton Keynes, UK). After carbonate precipitation with BaCl₂ and filtration, the δ¹³C of evolved CO₂ was analysed by an EA-IRMS. We also measured the Δ¹⁴C of evolved CO₂ by AMS analyses as described above for planted subsoil cores of 0–20 and 40–60 cm depth, and unplanted subsoil cores of 40–60 cm depth. Given the high cost of AMS analyses, Δ¹⁴C–CO₂ measurements were restricted to the cambisol and andosol that a priori featured the most contrasting mineral reactivity (n = 24 incubations).

After the incubation, each planted soil core was sieved at 2 mm. Roots retained by sieve and all visible roots in sieved soil were handpicked and washed with tap water. Shoot, root and soil materials were dried at 60 °C, weighed, grounded and analysed separately for C and N concentrations and δ¹³C using an EA-IRMS as described above. For the cambisol and andosol, Δ¹⁴C of root biomass in subsoil cores was measured as described above.

Isotopic partitioning and calculations

For both incubation series, the continuous isotopic labelling of plants with ¹³C-depleted air allowed us to partition total respiration into its soil and plant/root sources (~25‰ average difference in δ¹³C). It was calculated using the following equations:

where R_total and δ¹³C_total are respectively the total CO₂ flux and its δ¹³C from plant/root-soil respiration at the microcosm/core scale; R_soil and δ¹³C_soil are respectively the CO₂ flux and its δ¹³C from microbial respiration of unlabelled native (pre-existing) SOC and root litter; and R_plant and δ¹³C_plant are respectively the CO₂ flux and its δ¹³C from respiration of labelled (recently fixed) plant/root-derived organic carbon (OC). For the first incubation series, we used as δ¹³C_plant the mass-weighted δ¹³C of the mesocosm shoot and living root biomass, assuming negligible ¹³C fractionation during whole-plant respiration⁸⁹. A parallel experiment running during the same period and using a common labelling system found that the δ¹³C of plant biomass was constant through time⁸⁶, ensuring that the labelling remained homogeneous throughout the experiment. For the second incubation series, we used as δ¹³C_plant the δ¹³C of living root biomass corrected by a δ¹³C fractionation factor of root respiration, which was assumed to be −0.61‰ for grass species based a previous study⁸⁷. For both series of incubations, we used as δ¹³C_soil the δ¹³C of CO₂ derived from the respiration of native SOC in unplanted controls.

In the second incubation series, the continuous dual-labelling (¹³C/¹⁴C) of plants allowed us to quantify the radiocarbon signature of native SOC respired by decomposers (Δ¹⁴C_soil) by partitioning the Δ¹⁴C signature of total respiration into its soil and root sources^17,20. It was calculated using the following equation:

where R_total and Δ¹⁴C_total are respectively the total CO₂ flux and its Δ¹⁴C of root-soil respiration at the core scale, Δ¹⁴C_plant is the Δ¹⁴C of root-derived OC respiration, and R_plant and R_soil are the CO₂ fluxes of respectively root-derived and soil-derived OC respiration as calculated in Eqs. 2, 3. Assuming no fractionation of ¹⁴C during root respiration, we used the Δ¹⁴C of living root biomass as Δ¹⁴C_plant.

As the root material harvested for topsoil was composed of both pre-existing root litter (unlabelled) and living roots (labelled) that could not be clearly visually sorted, we used an isotopic partitioning method to estimate living root biomass (Root_living) for each planted topsoil core using the following equation:

where Root_total and δ¹³C_total are respectively the biomass and δ¹³C of both dead and living roots; and δ¹³C_living and δ¹³C_dead are the δ¹³C of respectively living and dead roots.

We also quantified the net rhizodeposition corresponding to root-derived OC remaining in the soil after microbial utilisation. It was estimated for each planted soil core using the following equation:

where SOC_total and δ¹³C_SOC-final are respectively the concentration and δ¹³C of SOC in planted soil core at the end of the experiment, δ¹³C_SOC-initial is the average δ¹³C of SOC from initial soil cores, and δ¹³C_root is the δ¹³C of living root biomass in planted microcosm.

In order to evaluate the uncertainty associated with our isotopic mixing model assumptions, we performed sensitivity analyses where we quantified the error in source proportion related to a 1‰ variation in the δ¹³C of both sources for each isotopic partitioning (Supplementary Methods). Low levels of uncertainty were found for every isotopic partitionings, providing evidence that our results were robust (Supplementary Tables 10–14, see Supplementary Methods for further details).

Assuming first-order decomposition kinetics, native SOC decomposition rate (k_SOC) was calculated as native SOC respiration (R_soil) divided by initial SOC concentration. The corresponding rhizosphere priming effect (RPE, in % of change in k_SOC for planted relative to unplanted soils) was calculated using the following equation:

The respiration of root-derived OC (R_plant) in the second incubation series includes both the autotrophic respiration of roots and the heterotrophic respiration of soil microbial decomposers derived from rhizodeposits. We thus used it here as a proxy of plant C allocation belowground to root activity, including root maintenance, growth and uptake of water and nutrients, as well as rhizodeposition of fresh OC to soil which is then quickly taken up and metabolised by rhizosphere microbes³³. This allowed us to assess how root effects on SOC decomposition could depend on plant inputs to soil of energy that could take different forms, including food substrates for microbes (metabolic energy), ligand exudates desorbing SOC from minerals (chemical energy) and root uptake of water and nutrients breaking up soil aggregates by physical disturbances (physical energy). After converting R_plant from an amount of C in CO₂ (mg) into an amount of glucose (mol), we calculated the plant energy investment into root activity (J dm⁻³ day⁻¹) as R_plant multiplied by the heat of combustion of glucose, that is 2807 kJ per mol⁷⁸. This estimation relies on the assumption that glucose composed most of the fresh photosynthate-C metabolised by both roots and soil microbes³².

Statistical analyses

All analyses were performed using R v3.4.3. We used a rotated principal component analysis (rPCA) to explore soil properties covariance and divergence between treatments. Analyses of variance were used to partition the variance explained by factors ‘soil layer’, ‘soil type’ and their interaction for the two first axis scores and soil properties.

For each incubation series, the responses of k_SOC, RPE and Δ¹⁴C_soil to predictors were assessed by regression for each treatment. We tested linear (Y = a + bX), polynomial (Y = a + bX + cX²) and power (Y = aX^b) regression functions, where Y is the response variable and X is the predictor. The k_SOC and RPE values were standardised for each treatment to a common high value of the following predictors: ‘respiration of plant-derived OC’ for the first incubation series, ‘living root density’ and ‘plant-derived OC respiration’ for the second incubation series. Additionally, we used analyses of covariance including the quantitative explanatory variables, the factors ‘soil layers’ and ‘soil type’, and their interactions as fixed factors to test their effects and quantify the proportion of variance they explain. To deal with the repeated measures design in both incubation series, we used linear mixed-effect models including ‘microcosm’ as random factor in regression and analyses of covariance.

We explored bivariate relationships between variation in SOC dynamics and SOM properties across depth using partial regression and correlation analyses controlling for soil types. Additionally, we performed an ordination of SOC dynamics variables constrained by soil biogeochemical predictors using a redundancy analysis.

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

Source: Ecology - nature.com