Boreal-forest soil chemistry drives soil organic carbon bioreactivity along a 314-year fire chronosequence

Following a wildfire, organic carbon (C) accumulates in boreal-forest soils. The long-term patterns of accumulation as well as the mechanisms responsible for continuous soil C stabilization or sequestration are poorly known. We evaluated post-fire C stock changes in functional reservoirs (bioreactive and recalcitrant) using the proportion of C mineralized in CO₂ by microbes in a long-term lab incubation, as well as the proportion of C resistant to acid hydrolysis. We found that all soil C pools increased linearly with the time since fire. The bioreactive and acid-insoluble soil C pools increased at a rate of 0.02 and 0.12 MgC ha⁻¹ yr⁻¹, respectively, and their proportions relative to total soil C stock remained constant with the time since fire (8 % and 46 %, respectively). We quantified direct and indirect causal relationships among variables and C bioreactivity to disentangle the relative contribution of climate, moss dominance, soil particle size distribution and soil chemical properties (pH, exchangeable manganese and aluminum, and metal oxides) to the variation structure of in vitro soil C bioreactivity. Our analyses showed that the chemical properties of podzolic soils that characterize the study area were the best predictors of soil C bioreactivity. For the O layer, pH and exchangeable manganese were the most important (model-averaged estimator for both of 0.34) factors directly related to soil organic C bioreactivity, followed by the time since fire (0.24), moss dominance (0.08), and climate and texture (0 for both). For the mineral soil, exchangeable aluminum was the most important factor (model-averaged estimator of −0.32), followed by metal oxide (−0.27), pH (−0.25), the time since fire (0.05), climate and texture (∼0 for both). Of the four climate factors examined in this study (i.e., mean annual temperature, growing degree-days above 5 ^∘C, mean annual precipitation and water balance) only those related to water availability – and not to temperature – had an indirect effect (O layer) or a marginal indirect effect (mineral soil) on soil C bioreactivity. Given that predictions of the impact of climate change on soil C balance are strongly linked to the size and the bioreactivity of soil C pools, our study stresses the need to include the direct effects of soil chemistry and the indirect effects of climate and soil texture on soil organic matter decomposition in Earth system models to forecast the response of boreal soils to global warming.

1 Introduction

Soil is the largest terrestrial carbon (C) reservoir (Scharlemann et al., 2014) and a major source of uncertainty in ecosystem C predictions (Shaw et al., 2014). Therefore, an advanced mechanistic understanding of soil C processes needs to be investigated and integrated into forecast models to reduce uncertainties in global C cycle feedback projections and to better predict the effects of global change on the soil C reservoir (Bradford et al., 2016; Schmidt et al., 2011). The maintenance of the vast soil C reservoir is partly under microbial control (Cotrufo et al., 2013) and could respond to variations in environmental conditions (Davidson and Janssens, 2006). Hence, the C quality-temperature hypothesis states that more “recalcitrant” soil organic matter should have a higher temperature sensitivity (Craine et al., 2010; Fierer et al., 2005). According to this hypothesis, it is important to distinguish the recalcitrant portion of the soil organic matter from the active portion in order to predict the impact of a rise in temperature on soil heterotrophic respiration. Furthermore, the bioreactive C fraction of soil organic matter is cycled on timescales relevant to global warming. Thus, quantifying the size of the bioreactive soil C reservoir and understanding the controlling factors of soil C bioreactivity (C_BioR) are key to inform models of the C cycle in the face of and to better anticipate global warming.

Wildfire is a major natural disturbance in boreal forests that drives the ecosystem C balance (Bond-Lamberty et al., 2007; Kurz et al., 2013) and is known to impact several soil properties, including organic matter quantity and quality (Certini, 2005; Knicker, 2007). Key soil properties, some evolving following a fire (e.g., soil acidity) and some not (e.g., particle size distribution), interact with climate and vegetation composition in complex causal direct and indirect relationships to regulate post-fire soil C accumulation (Andrieux et al., 2018). A saturation of soil C accumulation, especially for its recalcitrant portion, is often observed in soils when the rates of organic matter input to the soil are increased (Stewart et al., 2007; Hassink, 1996). Saturation of recalcitrant C is believed to come from the finite capacity of stabilization mechanisms in soils, such as chemical protection from the decomposition of soil organic matter by clay surfaces (Six et al., 2002). At a steady state and when the accumulation of physicochemically stabilized soil C occurs, soil C can only accumulate as nonprotected C (Castellano et al., 2015) that is more prone to decompose quickly. However, the long-term patterns of change in soil C quality and the accumulation pattern of recalcitrant and bioreactive C pools as well as the mechanisms responsible for continuous accumulation or the stabilization of soil C reservoirs are poorly known and have not been explicitly integrated into soil biogeochemistry (Luo et al., 2016). Most models of soil C dynamics divide soil organic matter into several conceptual pools and simulate decomposition as a first-order decay process (Luo et al., 2016). As part of this study, we characterized the acid-insoluble and bioreactive soil organic C pools (C_slow and C_fast, respectively; expressed as stocks) that accumulate following a wildfire. The acid-insoluble soil C fraction is assumed to be recalcitrant or resistant to biological degradation (Paul et al., 2006; Xu et al., 1997). Hereafter, we define C_BioR as the proportion of C mineralized in CO₂ by microbes at a constant temperature and constant water content over a long period of time as a relative measure of soil C lability (Laganière et al., 2015; Xu et al., 1997). Assessing the sizes of resistant and bioreactive soil C fractions through a fire chronosequence would help modelers to enhance the current and future C balance of landscapes prone to wildfires.

Besides its direct role in C cycling, climate has been shown to be an indirect predictor of soil C storage (quantity and accessibility for microbial decomposition) through its effects on geochemistry (Doetterl et al., 2015). In addition, vegetation types determine the quantity, quality and vertical distribution of soil litter inputs and so lead to differential mechanisms of soil C protection and stabilization (Jobbagy and Jackson, 2000; Laganière et al., 2017), with the moss stratum being a major source of soil C inputs in boreal ecosystems (Preston et al., 2006). Although many of these processes have been investigated separately, we are not aware of any empirical study so far that has quantified all these processes simultaneously and assessed the relative contribution of climate, time since fire (TSF), vegetation attributes and soil physicochemistry to soil C_BioR.

The objectives of this study are to fill these knowledge gaps by quantifying changes in boreal-forest soil bioreactive C stock with TSF (from 2 to 314 years) and disentangling the direct and indirect relative contributions of climate, moss dominance, soil particle size distribution and soil chemical properties (pH; exchangeable manganese, Mn, and aluminum, Al; and metal oxides) to soil C_BioR across the spruce–feather moss bioclimatic domain in eastern North America. Focusing on the complex interplay between climatic and non-climatic factors and their direct or indirect influence on soil C_BioR, we addressed the following questions. (i) Does soil C_BioR reservoir change with the soil organic C accumulation observed with TSF? (ii) To what extent do direct and indirect relationships among TSF, climate, physicochemical soil properties and bryophyte dominance influence soil C_BioR? We framed our study within the state factor model of ecosystems (Amundson and Jenny, 1997), which emphasizes soil physicochemical properties understood to be important to the pedogenesis of podzolic soils (Schaetzl and Anderson, 2005) that occur on the sampled sites. From there, our general hypotheses are that once site factors such as overstory composition, surficial deposits and soil drainage are accounted for, as they were in the present study, (1) soil C_BioR increases as forest stands get older, leading to a buildup of soil bioreactive C stock under the cold conditions of the boreal forest; (2) alternatively, if the bioreactive soil C stock reaches a new equilibrium because of rapid turnover, the proportion of bioreactive C stock should decline as total soil C stocks increases with TSF; and (3) soil C_BioR is primarily controlled by TSF and moss dominance in the O layer and by soil physicochemistry in the mineral soil (see Sect. 2.3 for detailed hypotheses).

Table 1 General characteristics of the sampling sites.

^∗ According to IUSS Working Group WRB (2015).

Figure 1 Map of the study area showing the location of the sample plots. Mean annual temperature (a) and mean annual precipitation (b) are interpolations from the 1981–2010 Canadian climate normals on a 10×10 km pixel grid (Chaste et al., 2018). The lower panel presents the location of the sample plots in relation to the time since fire (year). NAD: North American Datum.

2 Materials and methods

2.1 Site selection, sampling design and fieldwork

To account for the effects of TSF and climate on soil C pools, we established sample plots across both a chronosequence and a climosequence (Fig. 1; see Andrieux et al., 2018, for a description of the study area). Using numeric forest inventory maps compiled by the Ministère des Forêts, de la Faune et des Parcs du Québec (MFFPQ), we selected stands with as many similarities as possible in terms of canopy composition (black spruce – Picea mariana – stands), surficial deposits (thick till) and mesic drainage conditions. The soils that develop under these cool and humid climate conditions with acidic litter inputs typically belong to the podzolic order (Table 1). In boreal forests, podzolic soils often have a thick organic surface horizon (O layer), an eluviated A horizon (Ae) from where leached materials accumulate in the illuviated B horizon (Fig. S1). Within mesic drainage conditions, soil texture was quite variable (Table 1). These stands were overlaid with fire maps produced by the MFFPQ and other published dendrochronological surveys (Belisle et al., 2011; Bouchard et al., 2008; Cyr et al., 2012; Frégeau et al., 2015; Le Goff et al., 2007, 2008; Portier et al., 2016) to establish the chronosequence. We assumed that the canopy composition of black spruce did not change significantly with time and that the forest cyclically returned to the dominance of black spruce after a fire, in so-called recurrent dynamics, as previously described for these forests in a paleological survey (Frégeau et al., 2015). Then, while studying this ecosystem with a single and cyclic successional trajectory and low vegetation diversity and by carefully selecting stable permanent site conditions, we guarded against pitfall conclusions associated with space–time substitutions when using a chronosequence approach (Walker et al., 2010; Kenkel et al., 1997).

For field inventory and soil sampling, we followed Canada's National Forest Inventory ground plot guidelines (Canadian Forest Service, 2008). Stand biophysical description and soil sampling took place in a single 314 m² circular plot (10 m radius) in each stand. Slope inclination and orientation were recorded from the center of each plot with a clinometer and a compass, respectively. Every 2 m along two orthogonal transects oriented following the main cardinal directions and for a total of 20 records per plot, the thickness of the O layer was measured on a sample taken with a soil auger, and the dominant moss types (Sphagnum spp. or feather mosses) were identified using 400 cm² microplots (Fig. S1 in the Supplement). After litter and green living mosses were removed, the O layer was sampled at the edge of the plot in three 400 cm² microplots that were spaced 15 m from each other, from which we extracted volumetric mineral soil samples (top 15 cm) with a metallic cylinder (diameter of 4.7 cm; height of 15 cm). One soil pit was dug at the plot edge and at the same location where we sampled one of the three O-layer samples, down to the bottom of the podzolic B horizon or to the bedrock when possible, for a soil description and to collect the mineral soil from 15 to 35 cm under the forest-floor–mineral-soil boundary, as well as in the top 15 cm B horizon with a metallic cylinder (diameter of 4.7 cm). The significant stone content at one site prevented us from sampling the mineral soil from 15 to 35 cm; thus, no analyses could be provided for this layer of soil. Samples were kept in the dark and brought to the lab within 15 d for each region. They were kept in the dark at 2 ^∘C until processing during the fall.

The fieldwork was conducted in 2015, from 15 June to 8 September. The sampling effort covered 72 sites in black-spruce-dominated stands where a fire had burned 2 to 314 years ago. Climate data were interpolated at the plot level using BioSim v10.3.2 (Régnière et al., 2013) together with 1981–2010 climate normal series (http://climat.meteo.gc.ca/, last access: 12 November 2018) from surrounding weather stations and considering local slope attributes measured in the field as correcting factors (Régnière, 1996). Soil characteristics are summarized below (Table 1).

2.2 Laboratory analyses

2.2.1 Soil preparation

First, we prepared a composite of soil materials obtained from every sampled microplot (N=3), by plot and soil layer (O layer or the top 15 cm of mineral soil), to create representative samples for each layer in each of the 72 sample plots. O-layer samples were sieved through a 6 mm mesh before being oven-dried (60 ^∘C), whereas mineral soil samples (either the top 15 cm of the mineral soil, mineral soil from 15 to 35 cm or the top 15 cm of the B horizon) were dried by air and passed through a 2 mm sieve. Bulk density was determined after weighing the dried samples, assuming there were no coarse fragments in the O layer, and corrected for fragments >2 mm for the mineral soil. A part of each sample was retained for soil incubation. We used the <2 mm fraction to determine pH, exchangeable cation and texture (mineral soil only for the latter). Finely ground subsamples (<0.5 mm) were used for C concentration, pyrophosphate-extractable Fe and Al (top 15 cm of the B horizon only), and acid hydrolysis analyses.

2.2.2 Soil physicochemistry

C concentration of each sample was analyzed by dry combustion (Skjemstad and Baldock, 2007) using a LECO TruMac (LECO Corp, St. Joseph, MI, USA). Exchangeable cations were extracted using a Mehlich 3 solution and were analyzed by inductively coupled plasma atomic emission spectroscopy (Ziadi and Sen Tran, 2007), using an Optima 7300 DV (Perkin Elmer Inc., Waltham, MA, USA). Pyrophosphate-extractable Fe and Al (i.e., organically complexed metals; Mpy, hereafter defined as metal oxides) were extracted with a 0.1 N Na₄P₂O₇ solution before analysis with the Optima 7300 DV (Courchesne and Turmel, 2007). O-layer and mineral soil pH were determined in a soil : water solution by weights of 1:10 and 1:2 (Hendershot and Lalande, 2007), respectively, using a pH meter (Orion 2 Star). Particle size distribution of the mineral soil was assessed using a standard hydrometer method (Kroetsch and Wang, 2007).

2.2.3 Incubation settings

Soil incubation followed the method described in Paré et al. (2011). We prepared a total of 215 microcosms (72 sites for three soil layers minus one sample in the 15–35 cm mineral soil). We used 9 g of oven-dried O-layer and 50 g of air-dried mineral soil. Dried soil was used to ensure that the initial incubation moisture conditions were similar. Soil samples were placed in 100 mL bottom-perforated plastic containers. The containers were previously filled with glass wool (to avoid material losses during moisture adjustments) and prewashed with HCl (0.1 M) followed by deionized water. The microcosms were saturated with deionized water, drained for 24 h at 2 ^∘C and weighed to determine their water-holding capacity. Over approximately 50 weeks of experiments, microcosms were placed under constant air temperature (26 ^∘C) and humidity (100 %) in a growth chamber, and, when necessary, deionized water was periodically added to adjust the soil moisture to 85 % of the water-holding capacity. Except during CO₂ production measurements, each microcosm was stored in a 500 mL Mason jar kept open to maintain aerobic conditions and to prevent CO₂ accumulation to toxic levels. A rubber septum was installed on the metal lid for gas sampling when measuring CO₂ production.

2.2.4 CO₂ production measurements

CO₂ produced by each microcosm was measured periodically (at days 20, 26, 48, 64, 108, 126, 154, 227, 264 and 340 for the O layer and at days 8, 14, 21, 29, 36, 43, 51, 57, 72, 79, 86, 101, 113, 140, 203, 238 and 358 for mineral soil layers; Fig. S2), using a LI-6400 portable photosynthesis system (LI-COR, Lincoln, NE, USA) connected to an N₂ carrier gas (LI-COR application note 134). The flow rate of the carrier gas was set to 100 mL min⁻¹ using the gas flow meter FMA1812A (Omega Engineering, Inc., Norwalk, CO, USA). Initial CO₂ measurements were taken directly after hermetically sealing a jar with a metal lid, and final CO₂ measurements were taken after 4 to 24 h, depending on the soil layer and the progress of the experiment; these measurements were carried out using 2.5 or 10 mL (for O-layer or mineral soil, respectively) of air volume, extracted from the jar headspace with a syringe through the rubber septum. This gas sample was injected through the carrier gas into the LI-6400 infrared analyzer. The CO₂ concentration (µmol mol⁻¹) was predicted using the linear regression of a sample's measured CO₂ peak against calibration curves obtained from a benchmark gas (CO₂ at 800 and 3000 ppm). The first measurement (initial CO₂ concentration) accounted for the CO₂ concentration of the ambient air when closing the jars. This value was subtracted from the final CO₂ concentration to account solely for the CO₂ produced by the microcosm.

All data were subsequently standardized to a 24 h period to provide a daily respiration rate and to calculate cumulative C mineralization (Fig. S2; Paré et al., 2006). In short, we applied the gas law to convert CO₂ concentration (µmol mol⁻¹) to a C mass basis, using a constant pressure at 101.3 kPa and the specific head space volume of each sample (total volume of the jar minus soil volume and container). Cumulative respiration was calculated according to the following Eq. (1):

$$\begin{array}{}\text{(1)}& M_t=M_{t-\mathrm{1}}+{\displaystyle \frac{\left(R_t+R_{t-\mathrm{1}}\right)}{\mathrm{2}}}\times \left(t-t_{-\mathrm{1}}\right),\end{array}$$

where M_t (mg CO₂-C) is the cumulative mass of C mineralized at time t, R_t (mg CO₂-C d⁻¹) is the daily respiration rate at time t and t is the Julian day (d). M_t was divided by the initial C mass (g) of each sample to compute the specific respiration rate (Rs; mg CO₂-C g⁻¹ C_org). Then, dividing Rs by 10 gave the percentage of the initial mass of soil C lost through microbial respiration or C_BioR.

2.2.5 Acid hydrolysis

We used acid hydrolysis as an index of biologically recalcitrant soil C (Xu et al., 1997), which has been proposed as an indicator of a slowly cycling soil C pool (Paul et al., 2006). Hydrolysis was carried out by refluxing 2 g of soil with 50 mL HCl (6 M) brought to the boiling point on a hot plate. We used a 2 h reaction time because the majority of soil organic matter is hydrolyzed during the first 2 h, and longer reaction times do not significantly change C release (Silveira et al., 2008; Xu et al., 1997). Acid-insoluble residues were separated from hydrolysates by filtering the solution on inert paper filters, rinsed three times with 50 mL of deionized water to remove any chlorine residues, oven-dried at 60 ^∘C overnight and weighed before C concentration analysis by dry combustion (see Sect. 2.3.2). Based on the total C concentration of the acid-insoluble residues and mass loss of the samples during the treatments, the hydrolysability (Plante et al., 2006) of a sample was calculated based on the following Eq. (2):

$$\begin{array}{}\text{(2)}& {\mathrm{C}}_{\mathrm{AI}}=\left({\displaystyle \frac{\left[{\mathrm{C}}_{\mathrm{AI}}\right]\times M_{\mathrm{AI}}}{\left[{\mathrm{C}}_{\mathrm{i}}\right]\times M_{\mathrm{i}}}}\right)\times \mathrm{100},\end{array}$$

where C_AI is the percentage of the acid-insoluble C (%), [C_AI] and [C_i] are the C concentration of the acid-insoluble residues and of the initial soil (%), and M_AI and M_i are the masses of the acid-insoluble residues and of the initial soil sample (g), respectively.

2.3 Ecological a priori hypotheses

According to our third general hypothesis and to address the complex interplay among climatic and non-climatic factors, we first selected the following environmental variables documented in the literature as being important drivers of soil C_BioR and the pedogenesis of podzolic soils: climate (temperature and water supply), soil texture, TSF, dominance of the moss functional type, soil pH, and the concentration of metal oxides and of exchangeable elements (Mn and Al; Wiesmeier et al., 2019; Schaetzl and Anderson, 2005). As in other ecosystems (Fierer et al., 2003; Salomé et al., 2010), the boreal-forest soil microbial community as well as the chemical and physical environment change with soil depth (Clemmensen et al., 2013; Hynes and Germida, 2013), suggesting different drivers of the decomposition process in the O layer and in the mineral soil layer. Hence, for each of these layers, we built two separate sets of a priori ecological hypotheses expressed as direct acyclic graphs (DAGs) representing different causal relationships among environmental variables and soil C_BioR (Fig. 2). Therein, we tested the validity of four competing a priori ecological hypotheses represented by a DAG. This hypothetico-deductive approach, in which each a priori hypothesis was supported by ecological knowledge, allows for testing an alternative causal explanation in a falsifiable form regarding the underlying mechanisms of soil C_BioR in the two soil layers. For each soil layer, the first hypothesis assumed only direct relationships between environmental variables and soil C_BioR (baseline hypotheses for the O layer and the mineral layer, O1 and MIN1 in Fig. 2) such that they mirrored the widespread assumptions used in soil C prediction models based on multiple regression or ANOVA analyses. In addition, framed within Jenny's factor model of soil formation (Jenny, 1994), these baseline hypotheses assumed independence among environmental variables. Alternatively, we formulated a priori competing hypotheses in which both direct and indirect effects among variables and soil C_BioR were explicit (alternative hypotheses for the O layer and the mineral layer, O2 and MIN2 in Fig. 2). Justifications for each a priori ecological hypothesis are listed below.

2.3.1 Baseline hypothesis for the O layer (O1)

This hypothesis assumes that soils that have developed under cooler conditions limiting microbial activity should have greater C_BioR (Laganière et al., 2015) once temperature constraints have been removed. Rainfall under good drainage conditions (such as in this study) should promote greater decomposition rates and hence lower C_BioR. Because wildfire induces polymerization and polycondensation of organic compounds, resulting in residues that are more resistant to biological degradation (Certini, 2005; Gonzalez-Perez et al., 2004; Knicker, 2007), TSF is expected to have a direct and positive effect on C_BioR which was anticipated to increase with TSF. Soil pH also had a direct effect on soil bioreactivity because it regulates the microbial community (Fierer and Jackson, 2006) and is a key determinant of the decomposition process (Prescott et al., 2000; Zhang et al., 2008). We expected decreasing C_BioR with decreasing pH because acidic soil conditions limit the activity of soil decomposers. Compared to Sphagnum spp., feather mosses are more palatable to microbes (Fenton et al., 2010; Lang et al., 2009), so we expected lower C_BioR with more Sphagnum. Also, Mn availability has been shown to be a good predictor of boreal-soil C stocks (Stendahl et al., 2017); hence, with Mn being a co-metabolic compound of lignin degradation, we assumed that Mn has a direct positive effect on C_BioR.

2.3.2 Alternative hypothesis for the O layer (O2)

As in hypothesis O1, TSF, pH, moss functional type and Mn had direct effects on C bioreactivity. However, this hypothesis differed from O1 in that TSF and moss dominance also had indirect effects on C_BioR through changes in pH conditions. We expected decreasing pH with increasing Sphagnum spp. dominance because some physiological characteristics of these organisms lead to environment acidification (Andrus, 1986). Also, in the short term, fire modifies pH through the liming effect (Gonzalez-Perez et al., 2004; Knicker, 2007). In the long term, soils acidify with TSF as a result of vegetation regrowth, which involves the exchange of protons against cations to maintain the physiological electroneutrality of the vegetation (Driscoll and Likens, 1982). Contrary to hypothesis O1, which postulated that climate and soil texture had direct effects on C_BioR, hypothesis O2 assumed that these drivers had only indirect effects on C_BioR through their influence on moss dominance. Indeed, we expected that Sphagnum spp. would dominate over feather mosses under wetter conditions induced by greater precipitation, fined-texture soils holding more water or both because of their greater dependence to high soil water content.

2.3.3 Baseline hypothesis for the mineral soil (MIN1)

This hypothesis assumes that there are only direct effects of environmental variables on C_BioR in the mineral soil. Climate and pH directly control the decomposition process. As the binding of organic matter with the mineral phase has been recognized as an important mechanism of C protection against decomposition (Doetterl et al., 2015; Kaiser et al., 2002; Porras et al., 2017), we assumed that there would be direct effects of soil texture and metal oxide contents on C_BioR. In the first years following a fire, the slow incorporation of charred residues from upper soil layers into the mineral soil could decrease the organic matter quality (Johnson and Curtis, 2001), resulting in a decrease of C_BioR with increasing TSF. Mn availability could directly modulate C_BioR (see O1), and exchangeable Al could impede microbial decomposition when in excess (Kunito et al., 2016).

2.3.4 Alternative hypothesis for the mineral soil (MIN2)

As an alternative to hypothesis MIN1, this hypothesis assumes that only TSF, pH, metal oxides, and Mn and Al have direct effects on C_BioR. Additionally, pH is assumed to decrease with TSF because of the imbalance in nutrient uptake caused by aggrading vegetation. Also, exchangeable cations are dependent on pH (Sanborn et al., 2011), and the decrease in pH favors the creation of organometallic complexes impeding microbial decomposition (Buurman and Jongmans, 2005; Porras et al., 2017). Contrary to hypothesis MIN1, which assumed direct effects of climate and soil texture on C_BioR, this hypothesis assumes that climate and soil texture have only indirect effects on C_BioR. The indirect effect of climate on C_BioR is mediated through its effect on mineral weathering (Doetterl et al., 2015) and the quantity of metal oxides leached from the upper soil layers (Schaetzl and Anderson, 2005). Compared to coarse-textured soils, fine-textured soils have more reactive surface sites that can bind additional Mn and Al ions (Petersen et al., 1996).

2.4 Calculations and data analyses

2.4.1 Index of moss dominance

In order to account for the effects of moss functional traits on C_BioR of the O layer, we differentiated between Sphagnum spp. and feather mosses, since they have different ecophysiological characteristics (Bisbee et al., 2001); e.g., feather mosses decompose faster than Sphagnum spp. (Fenton et al., 2010; Lang et al., 2009). Based on Nalder and Wein (1999), we calculated an index of moss dominance (IMD) using Eq. (3):

$$\begin{array}{}\text{(3)}& \text{IMD}={\displaystyle \frac{{\mathrm{O}}_{sph}}{{\mathrm{O}}_{sph}+{\mathrm{O}}_{pl}+{\mathrm{O}}_h+{\mathrm{O}}_{pt}}},\end{array}$$

where O is the sum of occurrence of a species in the 20 microplots (see Sect. 2.1), sph is Sphagnum spp., pl is Pleurozium schreberi (Brid.) Mitt., h is Hylocomium splendens (Hedw.) Schimp. and pt is Ptilium crista-castrensis (Hedw.).

Feather mosses dominate the moss stratum when the IMD tends toward 0, whereas Sphagnum spp. dominates the moss stratum when the IMD tends toward 1. Some sites (n=5) that recently had fires did not have any moss species regrowth at the time of the fieldwork. For these sites, we set the IMD to 0.

2.4.2 Soil C quality and bioreactivity

First, we wanted to estimate variation in the size of the bioreactive and recalcitrant soil C pools (C_fast and C_slow, respectively, expressed as stocks) with TSF. For each soil layer (O layer, the top 15 cm of mineral soil and 15 to 35 cm of mineral soil), we scaled up to the plot scale of the cumulative proportion of C mineralized at the end of the incubations and the proportion of acid-insoluble C using Eqs. (4) and (5):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\mathrm{C}}_{\text{fast}}={\displaystyle \frac{{\mathrm{C}}_{\text{BioR}}}{\mathrm{100}}}\times \mathrm{C}\times D_{\mathrm{B}}\times h,\text{(5)}& {\displaystyle }{\mathrm{C}}_{\text{slow}}={\displaystyle \frac{{\mathrm{C}}_{\mathrm{AI}}}{\mathrm{100}}}\times \mathrm{C}\times D_{\mathrm{B}}\times h,\end{array}$$

where C_fast and C_slow are the bioreactive and recalcitrant soil C pools (Mg ha⁻¹), C_BioR is the percentage of the initial mass of soil C lost through microbial respiration (%), C is soil C content (%), D_B is the bulk density (g cm⁻³), h is the soil depth (i.e., mean depth based on 20 measurements per plot for the O layer; cm) and C_AI is the acid-insoluble C fraction (%). Hereafter, the total C stock (C_tot), C_fast or C_slow pool size represents, within each plot, the sum of each C pool across all soil layers.

Secondly, in order to express the qualitative (relative) changes in soil C in relation to environmental variables, we used the proportion of the initial mass of soil C lost through microbial respiration as an index of C lability (see Sect. 2.2.4). In the statistical analyses, we considered the whole mineral soil (in the top 15 cm and from 15 to 35 cm) as a single soil layer by calculating the weighted mean by depth for all mineral soil variables.

2.4.3 Statistical analyses

First, we evaluated post-fire C stock changes in functional reservoirs (bioreactive versus recalcitrant) using the linear regression of C stocks against TSF. Preliminary analyses with generalized additive models and piecewise regressions did not show any significant nonlinear or segmented relationships. Secondly, we quantified direct and indirect causal relationships among variables and C_BioR using confirmatory path analysis with directional-separation tests (Shipley, 2000a), according to the set of alternative a priori hypotheses (Fig. 2). Path analysis was used together with Fisher's C test (Shipley, 2000b) as a simultaneous test of independence for a model basis set (i.e., all nonadjacent pairs of variables defined as claims of independence). This led us to quantify how our data supported each hypothetical DAG and to identify whether some hypotheses would be rejected based on a robust statistical test (Shipley, 2009). Fisher's C statistic was compared with a χ² distribution with 2 k degrees of freedom (where k is the number of claims of independence in a model basis set). We rejected a causal model at the significance level α=0.05 when p<α. Prior to analyses, we standardized (reduced and centered) all variables to quantify their relative contribution to the variability of soil C_BioR.

The fit of each DAG for every soil layer (O and mineral soil) was compared using a model selection approach together with the second-order Akaike information criterion (AICc) in order to account for small sample sizes (Shipley, 2013). For model selection, we used the relative AICc difference with the “best” model or relative weight (Symonds and Moussalli, 2010). To avoid having latent variables in the models and because we had no a priori knowledge about which specific climate, texture or exchangeable elements should be used for testing the validity of each hypothesis, we used the cross product of four climatic variables (mean annual temperature: MAT; growing degree-day above 5 ^∘C: GDD5; mean annual precipitation: MAP; and water balance: WB), three soil texture variables (sand percentage, silt percentage and clay percentage) and two exchangeable elements (Al and Mn, only for mineral soil). Therefore, we tested 12 and 24 model combinations for each hypothesis and DAG involving the O layer and the mineral soil, respectively. Given that each soil layer had two alternative causal hypotheses, we then compared 24- and 48-candidate DAG models using a model selection procedure for the O layer and mineral soil, respectively. Model-averaged estimates were calculated by multiplying each estimate within each model with the corresponding Akaike weight and by summing the resulting values across all models; this allowed all models to influence model-averaged estimates. By doing so, we guarded against making arbitrary decisions about which model should be considered. We used the “ggm” package to compute Fisher's C statistic (Marchetti et al., 2015). All calculations and statistics were made using the R software version 3.4.3 (R Core Team, 2017).

Table 2 Post-fire soil C pool size and accumulation rates.

^∗ TSF: time since fire (year). p values <0.05 in bold.

Figure 2 Path models for each of the multivariate causal hypotheses. Arrows indicate direct causal relationships. Climate: climate variable; texture: texture variable; TSF: time since fire; pH: pH of the O layer (O1 and O2) or of the top 35 cm of mineral soil (MIN1 and MIN2); IMD: index of moss dominance; Mn: exchangeable manganese of the O layer; Mpy: pyrophosphate-extractable metals; E_ex: exchangeable element (Al or Mn) of the mineral soil.

3 Results

3.1 Post-fire soil C pool size

Total soil C stock (C_tot, i.e., the sum of O-layer and mineral C stocks in the top 35 cm), the size of the recalcitrant C pool and the size of the bioreactive C pool (C_slow and C_fast, respectively) all increased linearly with TSF (Fig. 3a). A minimum C_tot value of 63 MgC ha⁻¹ was observed for a 100-year-old stand, which is close to the C_tot value of 66 MgC ha⁻¹ of the youngest (2-year-old) stand. A maximum C_tot value of 305 MgC ha⁻¹ was observed for a 283-year-old stand. On average, the C_slow size was 6-fold bigger than the C_fast size (Table 2). C_fast values ranged from 5 to 25 MgC ha⁻¹ for a 100-year-old stand and for a 91-year-old stand, respectively. C_slow values were 29 and 175 MgC ha⁻¹ for a 2-year-old stand and for a 91-year-old stand, respectively.

Figure 3 Carbon quality as a function of the time since fire (TSF). The upper panel shows the total soil C reservoir (C_tot), the recalcitrant C pool (C_slow) and the bioreactive C pool (C_fast) sizes as a function of TSF (a, on the left), and the kernel density for each pool (a, on the right). The lower panel shows the proportion of C_slow and C_fast relative to C_tot as a function of TSF (b, on the left) and their kernel density (b, on the right). ${}^{***}$ p<0.001; ${}^{**}$ p<0.01; n.s. p>0.05.

Using these simple linear trends, C_tot accumulated faster than C_slow and C_fast (F${}_{\mathrm{3},\mathrm{209}}=\mathrm{257.6}$, p<0.001; Table 2 and Fig. 3b). Our data indicate that the overall soil C quality did not vary quantitatively with TSF (R²<0.01, p≥0.81 for both C pools; Table 2) because the proportion of C_slow and C_fast remained constant over the time span of the fire chronosequence (Fig. 3b). These general trends are mostly influenced by the size of the C_slow and C_fast pools of the O layer, being 5 times and 2.4 times larger than the top 35 cm of the mineral soil one (i.e., sum of the two mineral soil layers), respectively (Table 2; Fig. S3). C_slow and C_fast decrease with soil layers from the surface soil horizon (O layer) to the deeper mineral soil, both in absolute size and proportion (Table 2; Fig. S3). Consistently with the whole dataset, the proportion of C_slow and C_fast do not vary quantitatively with TSF for all the soil layers analyzed separately (Table 2). The size of the C_slow pool increases linearly with TSF in the O layer only (R²=0.09, p=0.01), not in mineral soil layers (p>0.07 for both mineral soil layers). The size of the C_fast pool increases linearly with TSF in the O layer (R²=0.12, p=0.003) and in the top 15 cm of the mineral soil (R²=0.05, p<0.05), not in the deepest mineral soil layer from 15 to 35 cm (p>0.21).

Table 3 Model fitness to the data for a priori hypotheses for the O layer. Models are sorted by the increasing second-order Akaike information criterion (AICc).

Hypothesis: model name; climate: climate variable; texture: texture variable; C statistic: statistic for Fisher's C test; df: degree of freedom; p: p value for Fisher's C test (when p<0.05, the model is not supported by the data); K: number of free parameters; AICc: second-order Akaike information criterion; ΔAICc: relative AICc difference with the model that best fitted the data (in bold); W: Akaike weight. MAT: mean annual temperature; GDD5: growing degree-day above 5 ^∘C; MAP: mean annual precipitation; WB: water balance.

3.2 O-layer path analysis and model selection

The causal structure of the O1 baseline hypothesis, which assumes that there are only direct effects of covariates on C_BioR in the O layer, was rejected for all candidate models (Fisher's C statistic >45, p<0.05; Table 3). Instead, the data better supported the causal structure of the O2 alternative hypothesis for all models (Fisher's C statistic <31, p>0.25; Table 3), which indicates that indirect effects among covariates and C_BioR in the O layer need to be accounted for in order to properly assess the variation structure in the data. The model selection procedure revealed that the data were best explained by one leading model (hereafter, the best model; Fisher's C statistic = 22.3, p=0.67; Table 3); this model was associated with O2, with MAP as the climate variable and clay content as the texture variable. The Akaike weight for this model (68 %) was about 8 times greater than the weight of the second-most supported model (8 %). The model-averaging procedure revealed that exchangeable Mn and pH of the O layer were the two covariates that had the strongest direct and positive effects on C_BioR of the O layer (both with an averaged path coefficient of pc = 0.34, p<0.01; Fig. 4 and Table S1). TSF was the second-most important relative driver with a significant direct and positive effect on C_BioR (pc = 0.24, p<0.05; Fig. 4 and Table S1 in the Supplement). Moss dominance had no significant direct effects on C_BioR (pc = 0.08, p>0.41; Fig. 4 and Table S1). In addition, both TSF and moss dominance had indirect effects on C_BioR through their influences on pH (TSF→pH: pc = −0.32, p<0.01; moss dominance→pH: pc = 0.30, p<0.01; Fig. 4 and Table S1). In addition, the model contained an indirect effect of climate (MAP) on C_BioR through its direct and negative effect on moss dominance (pc = −0.25, p<0.05; Fig. 4 and Table S1). We detected no any effect of texture (clay content) of the mineral soil on moss dominance (−0.01 < pc < 0.05., p>0.43).

Figure 4 Model that best fitted the data to explain carbon bioreactivity (C_BioR) in the O layer. Arrows indicate direct causal relationships. The numbers are standardized averaged path coefficients obtained using model averaging (see Table S1 and text for further details). MAP: mean annual precipitation; clay: clay content in the top 35 cm of mineral soil; Mn: exchangeable manganese; pH: pH of the O layer; TSF: time since fire; IMD: index of moss dominance. ${}^{\ast }p<\mathrm{0.05}$; ${}^{\ast \ast }p<\mathrm{0.01}$.

By allowing all of the models (O1 and O2) to influence coefficient estimates, the model-averaging procedure indicated that the most important variables exerting a direct control over C_BioR of the O layer were as follows, by decreasing importance: pH and Mn, TSF, and moss dominance (Table S1). Moreover, we could not detect any direct effect of climatic and texture variables tested in this study on C_BioR in the O layer.

Table 4 Model fitness to the data for a priori hypotheses for the mineral soil. Models are sorted by the increasing second-order Akaike information criterion (AICc).

Hypothesis: model name; climate: climate variable; texture: texture variable; exchangeable element: exchangeable-element variable; C statistic: statistic for Fisher's C test; df: degree of freedom; p: p value for Fisher's C test (when p<0.05, the model is not supported by the data); K: number of free parameters; AICc: second-order Akaike information criterion; ΔAICc: relative AICc difference with the model that best fitted the data (in bold); W: Akaike weight. MAT: mean annual temperature; GDD5: growing degree-day above 5 ^∘C; MAP: mean annual precipitation; WB: water balance.

3.3 Mineral soil path analysis and model selection

The causal structure of the MIN1 baseline hypothesis, which assumed that there were only direct effects of covariates on the C_BioR of the mineral soil, was rejected (Fisher's C statistic >52, p<0.05; Table 4). Instead, the data best supported the causal structure implied by the MIN2 alternative hypothesis indicating that, similar to the O layer, indirect effects among covariates need to be accounted for in order to properly assess variation in the C_BioR of the mineral soil. The model selection procedure revealed that the data were best explained by one leading model (hereafter, the best model; Fisher's C statistic = 27.76, p=0.27; Table 4). This model was associated with MIN2, with WB as the climate variable, clay content as the texture variable and Al as the exchangeable-element variable. The Akaike weight for this model (47 %) was about 3 times greater than for the second-most supported model (14 %). The model-averaging procedure revealed that exchangeable Al had the strongest direct and negative effect on the C_BioR of the mineral soil (pc = −0.32, p<0.001; Fig. 5). Metal oxide content (pc = −0.27, p<0.05) and pH (pc = −0.25, p<0.05) were the second-most influential drivers with significant and negative direct effects on the C_BioR of the mineral soil. TSF had a small positive direct effect on the C_BioR of the mineral soil, but this relationship was not significant (pc = 0.05, p>0.36). In addition, pH induced two indirect effects on the C_BioR of the mineral soil, i.e., through its negative and direct effects on Al and Mpy (pH→Al: pc = −0.24, p<0.01; pH→Mpy: pc = −0.34, p<0.01). Clay content had an indirect effect on the C_BioR of the mineral soil, through its direct and positive effect on exchangeable Al (pc = 0.17, p<0.05). Also, water balance had a weak indirect effect on the C_BioR of the mineral soil through its direct effect on Mpy (pc = 0.11, p=0.07).

By allowing all the models (MIN1 and MIN2) to influence estimates, the model-averaging procedure indicated that the most important variables tested in this study and exerting direct control over the C_BioR of the mineral soil were as follows, by decreasing importance: exchangeable Al, metal oxide contents, pH and TSF (Table S2). Moreover, we failed to detect any direct effect of climate or mineral soil texture on C_BioR.

Figure 5 Model that best fitted the data to explain the carbon bioreactivity (C_BioR) in the top 35 cm of the mineral soil. Arrows indicate direct causal relationships. The numbers are standardized averaged path coefficients obtained using model averaging (see Table S2 and text for further details). Al: exchangeable aluminum in the top 35 cm of the mineral soil; pH: pH of the top 35 cm of the mineral soil; TSF: time since fire; Mpy: metal oxide content in the top 15 cm of the B horizon; WB: water balance; clay: clay content of the top 35 cm of the mineral soil. ${}^{\ast }p<\mathrm{0.05}$; ${}^{\ast \ast }p<\mathrm{0.01}$; ${}^{\ast \ast \ast }p<\mathrm{0.001}$.

4 Discussion

4.1 Post-fire soil C quality

Most of the studies on post-fire soil C have focused on immediate or short-term responses and found that fire affects soil C quality by creating profound changes in the structure of soil organic matter compounds through thermal oxidation (Certini, 2005; Gonzalez-Perez et al., 2004). By using a long-term chronosequence of TSF ranging from 2 to 314 years, our study provides new insights into the understanding of the trajectory of changes in soil C quality following a fire, over hundreds of years. Our estimates of the size of fast- and slow-cycling soil C pools and our results indicate that (i) both pools accumulate with TSF and (ii) the proportion of each C pool remains constant with TSF relative to total soil C stock (Fig. 2 and Table 1). These results do not necessarily imply that fire has no effect on soil C functional pools, because our chronosequence has a low resolution for the first few years following a fire, but they rather suggest that such changes, if present, are not long-lasting. Our results also highlight that the accumulation process of the bioreactive soil C reservoir does not reach an equilibrium, at least not in the first 3 centuries following a fire. Instead, environmental conditions limiting decomposition, such as cold temperatures under a thickening O layer which developed with TSF, could have slowed down labile C degradation and allowed its accumulation (Kane et al., 2005). Our results also emphasize that changes in the size of the soil functional reservoirs with TSF are stratified within the soil profile. This pattern is consistent with the fact that the impact of a fire on soil C stock is limited to surface soil horizons (Andrieux et al., 2018).

4.2 Control mechanisms of the soil carbon bioreactivity

This study shows that soil C_BioR is driven by several climatic and non-climatic variables, some being common both for the O layer and mineral soil and others not, suggesting that different mechanisms may be involved in the control of the decomposition process in the O layer and in the mineral soil (Shaw et al., 2015; Ziegler et al., 2017).

4.2.1 Soil carbon bioreactivity in the O layer

Our results suggested that pH and exchangeable Mn are important drivers of C_BioR in the O layer. Boreal evergreen coniferous species generate high-lignin litter and forest floor layers (Laganière et al., 2017). This is reflected in the high proportion of acid-insoluble C of O-layer samples (among all the O-layer samples, mean ± SD = 73±5 %; Fig. S3). Therefore, soil C cycling in boreal forests depends on the capacity of microbes to depolymerize lignin. Microorganisms in the acidic soils of this ecosystem are dominated by fungi that use metalloenzymes – such as Mn peroxidases – to metabolize lignin (Pollegioni et al., 2015) or are white-rot fungi (Basidiomycota) equipped with enzymes that oxidize lignin (Cragg et al., 2015). Soil C stocks in the boreal-forest humus layer have been found to be negatively correlated with exchangeable Mn availability (Stendahl et al., 2017). In our study, exchangeable Mn of the O layer was positively correlated with C_BioR, suggesting that increasing Mn availability stimulates organic matter breakdown and that an Mn bottleneck in soil C cycling may be present (Kranabetter, 2019). We also observed direct and positive causal relationships between pH and C_BioR of the O layer, indicating that acidic soil conditions limit soil C mineralization (Prescott et al., 2000). Bacterial respiration and microbial community composition were found to be strongly determined by soil pH in the forest soil (Bååth and Anderson, 2003). We found that pH of the O layer decreased with TSF. Alongside the direct and positive effect of TSF on the C_BioR of the O layer, our results indicate that dynamic processes constrained by chemical soil properties shifting with stand development after burning (e.g., pH) drive the nature of soil organic matter and potentially the rate of C losses by heterotrophic respiration from boreal-forest soils.

Altogether, these results emphasize the need to include both soil chemistry and biological mechanisms in models of soil C cycling to better anticipate the role played by the boreal forest in C cycle climate feedbacks. Soil C cycling in mechanistic models of forest C dynamics often assumes that climate drives decay and the transfer rate of and between soil C pools (see Deluca and Boisvenue, 2012). Based on our results, we argue that chemical drivers of soil organic matter decomposition, such as exchangeable Mn concentrations and pH, might also be used to modulate soil C dynamics in such models, and we especially advise accounting for temporal shifts in soil pH occurring with stand development.

We did not detect any direct effect of climate on soil C_BioR in the O layer. This finding is consistent with the results of unchanged soil C stocks with in situ experimental warming worldwide (van Gestel et al., 2018). Furthermore, when synthesizing data of in situ experimental warming, Carey et al. (2016) found no change in the soil respiration rate for warmed compared to control plots at the global scale, whereas changes were found to be significant for the boreal biome. The cumulative C mineralization of incubated soils in our study was not modulated by in situ temperature, which supports the results of Carey et al. (2016) for their entire dataset but not for the boreal biome-restricted dataset. However, Carey et al. (2016) did not study soils from the Canadian Boreal Shield.

4.2.2 Soil carbon bioreactivity in the mineral soil

As in the O layer, our results highlight the role of pH as a regulator of C_BioR in the mineral soil. In addition to having a direct effect on C_BioR, pH also had two indirect effects. The first indirect effect is through the stimulation of metal oxide production with increasing acidic conditions. We observed that low-pH conditions correlated positively with higher metal oxide contents, which in turn correlated negatively with C_BioR in the mineral soil. This result is consistent with previous findings showing the role played by pH in mineral weathering and the preservation of C from decomposition through organometal complexation (Andrieux et al., 2018). The second indirect effect of pH on C_BioR in the mineral soil is mediated through exchangeable Al only, not through Mn (Table S2). Microbes are vertically stratified within the soil column (Clemmensen et al., 2013; Ekschmitt et al., 2008; Hynes and Germida, 2013), with fungi populating the upper soil layers because of their higher metabolic-oxygen demand compared with bacteria, which can more easily dwell in the less-oxygenated deeper soil layers. Our results suggest that, contrary to the O layer, oxidative depolymerization of lignin compounds mediated by Mn peroxidases may not be a major process for C cycling in the mineral soil (see above). Instead, the negative effect of pH on exchangeable Al, together with the negative effect of exchangeable Al on C_BioR in the mineral soil, indicates that low-pH conditions favor a greater abundance of exchangeable Al, which in turn impedes organic matter decomposition. These findings are consistent with the observed pH-dependent Al toxicity that slowed microbial catabolic activities in acidic forest soils in Japan (Kunito et al., 2016) and in laboratory experiments (Wood, 1995). Our study goes one step further in that we show that the content of exchangeable Al is directly related to soil texture (especially clay content) in these podzolic soils. This supports the hypothesis that exchangeable Al bound to fine mineral particles, such as clay, might act as a source of stored Al that can be mobilized and complexed with C to impede decomposition.

Contrary to the O layer, we found that TSF was only weakly correlated with mineral soil C_BioR and pH, and these relationships were not significant. We also found that the indirect effect of climate (correlation between water balance and metal oxides) on C_BioR was marginal. These results indicate that effects of TSF (direct and indirect) and water availability (indirect) on C_BioR are restricted to surface organic horizons. Our results support the idea that properties of the organic layer are more likely to be affected by fire because they are directly exposed to surface heating (De Bano, 1990) and that the thick humus layer of boreal-forest soils (Table 1) protects deeper soil layers from shifts in environmental conditions. The fact that roots of black spruce mostly develop in the top soil (Yuan and Chen, 2010) could explain why we did not observe shifts in pH of the mineral soil with TSF.

4.3 Implication for carbon cycling and research needs

Theoretically, soil C dynamics can be predicted through a knowledge of soil C pool sizes, changes in inputs and sensitivity to environmental factors (Luo et al., 2016). Understanding the bioreactivity of the large boreal-forest soils C reservoir is key to predicting the future global C cycle in the face of global warming. As illustrated in our study, some factors may act as C stabilization agents of soil C (i.e., metal oxides binding organic matter), while others may contribute to accelerating or slowing down the rate of soil organic matter biological processing (i.e., exchangeable Mn as a co-metabolic compound of lignin degradation; when in excess, toxicity of exchangeable Al for microbes; soil acidity regulating the microbial activity), and finally, others are related to the quality of organic matter inputs to the soil (i.e., type of moss flora; low-quality organic materials left after a fire) or simply to the time required by the system to adjust and reach the steady state (i.e., causal relationships between TSF and pH or pH and metal oxides). Our best models showed that the climatic conditions experienced in situ, expressed here as temperature and water availability, had no direct effect on in vitro soil C_BioR. Moreover, the indirect effects of climate on soil C_BioR are limited to water supply factors, not to temperature. Best models also reveal direct and indirect effects on C_BioR of both site properties and factors that evolve with TSF. Understanding and predicting changes in the soil chemistry is therefore a key challenge that remains to be addressed in future works in order to improve our understanding of soil C balance with global change. Our results are in agreement with Davidson and Janssens (2006) and Davidson (2015), suggesting that improvements to Earth system models (ESMs) may arise from integrating the long-term effects of climate on soil properties with the environmental constraints on microbiological degradation of soil organic matter.

The results of this study identified new pathways for the control mechanisms of soil C_BioR that could help to predict the response of boreal-forest soils to global change. While Earth system models commonly focus on a temperature dependence of soil C decomposition (Bradford et al., 2016), our study showed, in agreement with Rasmussen et al. (2018), that key soil properties, because of their relationship to soil C bioreactivity, could improve ESMs for modeling soil C dynamics in relation to climate change. In particular, our study shows that predictive models need to include the direct effects of soil chemistry and the indirect effects of climate and soil texture on soil C_BioR. Moreover, while some factors (metal oxides and TSF) were found to affect both soil C_BioR (this study) and soil C stocks (Andrieux et al., 2018), at the same time, other factors did not have such effects (types of mosses and pH). For example, moss dominance had a direct effect on C stock (Andrieux et al., 2018) but not on C_BioR (this study) in the O layer.

The path analyses and model selection procedure used in our study have made it possible to distinguish direct from indirect effects of ecological drivers on soil C dynamics. We found that the local climate shaped soil C_BioR indirectly through effects on moss dominance and on metal oxides and that of the four climatic variables examined, only the variables related to water supply – and not temperature – significantly but indirectly affected soil C_BioR. This suggests that the forecasted increase of 11 % precipitation by the end of this century in eastern North America (IPCC, 2013) would indirectly modulate soil C stocks (Andrieux et al., 2018) and soil C_BioR (this study), together with the indirect effects of climate on the mechanisms of soil C stock and bioreactivity. How the boreal-ecosystem C balance will evolve in the context of global change might be assessed through further research focusing on the changes in soil physicochemical reactions pertaining to the mechanisms of soil organic matter decomposition and stabilization (Thornley and Canell, 2001).

5 Conclusions

Our study aimed to quantify the long-term post-fire changes occurring in the functional soil C pools and to disentangle the direct and indirect relative contribution of climate, moss dominance, soil particle size distribution and soil chemistry (pH, exchangeable Mn and Al, and metal oxides) on soil C_BioR. Using a chronosequence approach, we show that labile and recalcitrant soil C pools both increased continually from 2 to 314 years after a fire. Changes in the size of the bioreactive and of the acid-insoluble C pools with TSF were found to be stratified within the soil profile and reveal that the impact of a fire on soil functional C reservoirs is limited to surface soil horizons. The main drivers of C_BioR varied with the soil layer considered. The breakdown of organic matter in the O layer was constrained by pH and exchangeable Mn, while that of the mineral soil was dependent on the availability of exchangeable Al, metal oxide content and pH. Moreover, our results suggest that for both soil layers, the complex interplay among biogeochemical covariates needs to be accounted for to assess the variation structure of C_BioR, with climate (water supply parameters only and not temperature) having only indirect effects. We argue that the direct and indirect effects that covariates have on the bioreactivity of soil organic carbon need to be integrated in models simulating C dynamics. This is key to forecasting the response of the enormous boreal-forest soil carbon pool to global warming.