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

. Following a wildﬁre, 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-ﬁre C stock changes in functional reservoirs (bioreactive and recalcitrant) us-ing the proportion of C mineralized in CO 2 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 ﬁre. The bioreactive and acid-insoluble soil C pools increased at a rate of 0.02 and 0.12 MgC ha − 1 yr − 1 , respectively, and their proportions relative to total soil C stock remained constant with the time since ﬁre (8 % and 46 %, respectively). We quantiﬁed 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 bioreac-tivity, followed by the time since ﬁre (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 ﬁre (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 .yr -1 and 0.12 MgC.ha -1 .yr -1 , respectively, and their proportions relative to total soil C stock remained constant with time since 15 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 characterise 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 20 both: 0.34) factors directly related to soil organic C bioreactivity, followed by 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 (modelaveraged estimator: -0.32), followed by metal oxide (-0.27), pH (-0.25), 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 indirect effect (O 25 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.

30
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 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 35 (Davidson and Janssens, 2006). Hence, the C-quality temperature hypothesis states that more "recalcitrant" soil organic matter should have 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 time scales relevant to global warming. Thus, quantifying the size of the bioreactive soil C reservoir and understanding the 40 controlling factors of soil C bioreactivity (CBioR) is key to inform models of C cycle in 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 fire (e.g., soil acidity) and some not (e.g., particle size C stock should decline as total soil C stocks increases with TSF; and 3) soil CBioR is primarily controlled by TSF and moss dominance in the O layer, and by soil physico-chemistry in the mineral soil (see section 2.3 for detailed hypotheses).

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) 90 and mesic drainage conditions. The soils that develop under these cool and humid climate with acidic litter inputs typically belong to the Podzolic order (Table 1). In boreal forests, podzolic soils often have 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; 95 Le Goff et al., 2007;Le Goff et al., 2008;Portier et al., 2016) to establish the chronosequence. We assumed that the black spruce canopy composition did not change significantly with time and that the forest cyclically returned to a black spruce dominance after 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 100 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 (NFI, 2016).
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 centre 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 105 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). 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 (ø = 4.7 cm, height = 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 110 podzolic B horizon or to the bedrock when possible, for 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 (ø = 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 days for each region. They were kept in the dark at 2°C until processing during the fall.

115
The fieldwork was conducted in 2015, from 15 June to 8 September. The sampling effort covered 72 sites in black spruce stands where 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/) 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).

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 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 top 15 cm of the 125 mineral soil, mineral soil from 15 to 35 cm or 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. 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 sub-samples (< 0.5 mm) were used for C concentration, pyrophosphate extractable Fe and Al (top 15 cm of the B horizon only) and acid hydrolysis 130 analyses.

Soil physico-chemistry
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.,

135
Waltham, MA, USA). Pyrophosphate extractable Fe and Al (i.e., organically complexed metals; Mpy, hereafter defined as metal oxides) were extracted with a 0.1N Na4P2O7 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).

Incubation settings
Soil incubation followed the method described in Paré et al. (2011). We prepared a total of 215 microcosms (72 sites x 3 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 145 moisture adjustments) and pre-washed with HCl (0.1 M) followed by deionized water. The microcosms were saturated with deionized water, drained for 24 hours 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 soil moisture to 85% of the water-holding capacity. Except during CO2 production measurements, each microcosm was stored in a 500-mL Mason jar kept open to 150 maintain aerobic conditions and to prevent CO2 accumulation to toxic levels. A rubber septum was installed on the metal lid for gas sampling when measuring CO2 production.

CO2 production measurements
CO2 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), 155 using a LI-6400 portable photosynthesis system (LI-COR ® , Lincoln, NE, USA) connected to an N2 carrier gas (LI-COR ® Application Note # 134). The flow rate of the carrier gas was set to 100 mL.min -1 using the gas flow meter FMA1812A (Omega Engineering, INC., Norwalk, CO, USA). Initial CO2 measurements were taken directly after hermetically sealing a jar with a metal lid and final CO2 measurements were taken after 4 h to 24 h, depending on the soil layer and the progress of the experiment; these measurements were carried out using a 2.5-mL or 10-mL (for O layer or mineral soil, respectively) air 160 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. CO2 concentration (µmol.mol -1 ) was predicted using the linear regression of a sample's measured CO2 peak against calibration curves obtained from benchmark gas (CO2 at 800 ppm and 3,000 ppm). The first measurement (initial CO2 concentration) accounted for the CO2 concentration of the ambient air when closing the jars.
This value was subtracted from the final CO2 concentration to account solely for the CO2 produced by the microcosm.

165
All data were subsequently standardized to a 24-hour 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 CO2 concentration (µmol.mol -1 ) 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):

170
where Mt (mg CO2-C) is the cumulative mass of C mineralized at time t, Rt (mg CO2-C.d -1 ) is the daily respiration rate at time t, and t is the Julian day (d). Mt was divided by the initial C mass (g) of each sample to compute the specific respiration rate (Rs; mg CO2-C.g -1 Corg). Then, dividing Rs by 10 gave the percentage of initial mass of soil C lost through microbial respiration, or CBioR.

175
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 slow-cycling soil C pool . Hydrolysis was carried out by refluxing 2 g of soil with 50 mL HCl (6M) brought to the boiling point on a hot plate. We used a two-hour reaction time because the majority of soil organic matter is hydrolyzed during the first two hours 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 180 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 section 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 ) of a sample was calculated based on the following Eq. (2): where CAI is the percentage of the acid-insoluble C (%), [CAI] and [Ci] are the C concentration of the acid-insoluble residues and of the initial soil (%), and MAI and Mi are the mass of the acid-insoluble residues and of the initial soil sample (g), respectively.

Ecological a priori hypotheses
According to our third general hypothesis and to address the complex interplay among climatic and non-climatic factors, we 190 first selected the following environmental variables documented in the literature as being important drivers of soil CBioR and pedogenesis of Podzolic soils: climate (temperature and water supply), soil texture, TSF, dominance of the moss functional type, soil pH, and 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 195 Germida, 2013), suggesting different drivers of the decomposition process in O-and mineral soil layers. Hence, for each of the O-and mineral soil 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 CBioR (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 200 a falsifiable form as regards the underlying mechanisms of soil CBioR in the two soil layers. For each soil layer, the first hypothesis assumed only direct relationships between environmental variables and soil CBioR (hypotheses 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 205 both direct and indirect effects among variables and soil CBioR were explicit (hypotheses O2 and MIN2 in Fig. 2). Justifications for each a priori ecological hypothesis are listed below.

Baseline hypothesis for the O layer, O1
This hypothesis assumes that soils that have developed under cooler conditions limiting microbial activity should have a greater CBioR (Laganière et al., 2015) once temperature constraints have been removed. Rainfall under good drainage conditions (such 210 as in this study) should promote greater decomposition rates, and hence lower CBioR. 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 CBioR 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 215 al., 2008). We expected decreasing CBioR 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 CBioR with more Sphagnum. Also, Mn availability has been shown to be a good predictor of boreal soil C stocks (Stendahl et al., 2017); hence, Mn being a co-metabolic compound of lignin degradation, we assumed that Mn has a direct positive effect on CBioR. 220

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 CBioR 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 225 (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 electro-neutrality of the vegetation (Driscoll and Likens, 1982). Contrary to hypothesis O1, which postulated that climate and soil texture had direct effects on CBioR, hypothesis O2 assumed that these drivers had only indirect effects on CBioR through their influence on moss dominance. Indeed, we expected that Sphagnum spp. would dominate over feathermosses under wetter conditions induced by greater precipitation, 230 fined-texture soils holding more water, or both, because of their greater dependence to high soil water content.

Baseline hypothesis for the mineral soil, MIN1
This hypothesis assumes that there are only direct effects of environmental variables on CBioR 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., 235 2017), we assumed that there would be direct effects of soil texture and metal oxide contents on CBioR. In the first years following 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 CBioR with increasing TSF. Mn availability could directly modulate CBioR (see O1), and exchangeable Al could impede microbial decomposition when in excess (Kunito et al., 2016).

Alternative hypothesis for the mineral soil, MIN2
240 As an alternative to hypothesis MIN1, this hypothesis assumes that only TSF, pH, metal oxides and Mn/Al have direct effects on CBioR. 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 favours 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 CBioR, this hypothesis assumes that climate 245 and soil texture have only indirect effects on CBioR. The indirect effect of climate on CBioR 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).
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 260 the fieldwork. For these sites, we set the IMD to 0.

Soil C quality and bioreactivity
First, we wanted to estimate variation in the size of the bioreactive and recalcitrant soil C pools (Cfast and Cslow, respectively, where Cfast and Cslow are the bioreactive and recalcitrant soil C pools (Mg.ha -1 ), CBioR is the percentage of initial mass of soil C lost through microbial respiration (%), C is soil C content (%), DB is the bulk density (g.cm -3 ), h is the soil depth (i.e., mean 270 depth based on 20 measurements per plot for the O layer; cm) and CAI is the acid-insoluble C fraction (%). Hereafter, the total C stock (Ctot), Cfast or Cslow 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 initial mass of soil C lost through microbial respiration as an index of C lability (see section 2.2.4). In the statistical analyses, we considered the whole mineral soil (in the top 15 cm and in the 15-to 35-cm layer) as a single soil layer 275 by calculating the weighted mean by depth for all mineral soil variables.

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 non-linear or segmented relationships. Secondly, we quantified direct and indirect causal relationships 280 among variables and CBioR 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 non-adjacent 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 285 with 2k 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 CBioR.
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

305
Total soil C stock (Ctot, 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 (Cslow and Cfast, respectively) all increased linearly with TSF (Fig. 3a) Using these simple linear trends, Ctot accumulated faster than Cslow and Cfast (F3,209 = 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 Cslow and Cfast remained constant over the timespan of the fire chronosequence (Fig. 3b).

315
These general trends are mostly influenced by the size of the Cslow and Cfast 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 (

O layer path analysis and model selection
The  Fig. 4 and Table S1). TSF was the second most important relative driver with a significant direct and positive effect on CBioR (pc = 0.24, p < 0.05; Fig. 4 and Table S1). Moss dominance had no significant 335 direct effects on CBioR (pc = 0.08, p > 0.41; Fig. 4 and Table S1). In addition, both TSF and moss dominance had indirect effects on CBioR 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 CBioR 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).

340
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 CBioR 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 CBioR in the O layer.

345
The causal structure of the baseline hypothesis MIN1, which assumed that there were only direct effects of covariates on the CBioR 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 alternative hypothesis MIN2 indicating that, similar to the O layer, indirect effects among covariates need to be accounted for assessing in order to properly assess variation in the CBioR of the mineral soil. The model selection procedure revealed that the data were best explained by one leading model (hereafter, "best" model; Fisher's C statistic = 350 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 three 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 CBioR 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 355 effects on the CBioR of the mineral soil. TSF had a small positive direct effect on the CBioR 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 CBioR 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 CBioR 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 CBioR of the mineral soil through its direct effect on 360 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 a direct control over the CBioR 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 CBioR.

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 two to 314 years, our study 370 provides new insights into the understanding of the trajectory of changes in soil C quality following 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 fire, but rather suggest that such changes, if present, are not long-lasting. Our results 375 also highlight that the accumulation process of the bioreactive soil C reservoir does not reach an equilibrium, at least not in the first three centuries following fire. Instead, environmental conditions limiting decomposition, such as cold temperatures under a thickening O layer developed with TSF, could have slowed down labile C degradation and allow 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 fire impact on soil C stock is limited to surface soil horizons 380 (Andrieux et al., 2018).

Control mechanisms of the soil carbon bioreactivity
This study shows that soil CBioR is driven by several climatic and non-climatic variables, some being common both for 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).

Soil carbon bioreactivity in the O layer
Our results suggested that pH and exchangeable Mn are important drivers of CBioR in the O layer. Boreal evergreen coniferous species generate high-lignin litter and forest floor layers . 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 390 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 CBioR, suggesting that increasing Mn availability stimulates organic matter breakdown and that an Mn bottleneck in soil C cycling may be present (Kranabetter, 395 2019). We also observed direct and positive causal relationships between pH and CBioR 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 CBioR 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 400 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 into models of soil C cycling to better anticipate the role played by 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 405 exchangeable Mn concentrations and pH, might also be used to modulate soil C dynamics in such models, and we especially advise to accounting for temporal shifts in soil pH occurring with stand development.
We did not detect any direct effect of climate on soil CBioR 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 soil respiration rate for warmed compared to control 410 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.

415
As in the O layer, our results highlight the role of pH as a regulator of CBioR in the mineral soil. In addition to having a direct effect on CBioR, 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 CBioR 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 organo-metal complexation 420 (Andrieux et al., 2018). The second indirect effect of pH on CBioR 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 greater need for metabolic oxygen 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 425 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 CBioR in the mineral soil, indicates that low pH conditions favor a greater exchangeable Al abundance, which in turn impedes organic matter decomposition. These findings are consistent with the observed pHdependent 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 exchangeable Al content is directly 430 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 to mineral soil CBioR 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

435
CBioR was marginal. These results indicate that effects of TSF (direct and indirect) and water availability (indirect) on CBioR 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 black spruce roots 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.

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 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 445 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 fire) or simply to the time require by the system to adjust and reach 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 450 here as temperature and water availability, had no direct effect on in vitro soil CBioR. Moreover, the indirect effects of climate on soil CBioR are limited to water supply factors, not to temperature. "Best" models also reveal direct and indirect effects on CBioR of both site properties and factors that evolve with TSF. Understanding and predicting changes in 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 455 improvements to 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 CBioR that could help to predict the response of boreal forest soils to global change. While earth system models (ESMs) 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 460 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 CBioR. Moreover, while some factors (metal oxides, TSF) were found to affect both soil CBioR (this study) and soil C stocks (Andrieux et al., 2018), at the same time, other factors did not have such effects (types of mosses, pH). For example, moss dominance had a direct effect on C stock (Andrieux et al., 2018), but not on

465
CBioR (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 CBioR 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 CBioR. This suggests that the forecasted increase of 11% 470 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 CBioR (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 physico-chemical reactions pertaining to the mechanisms of soil organic matter decomposition and stabilization (Thornley and Canell, 2001).

Conclusion
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 CBioR. Using a chronosequence approach, we show that labile and recalcitrant soil C pools both increased continually from 2 to 314 years after fire. Changes in the size of the bioreactive and of 480 the acid-insoluble C pools with TSF were found to be stratified within the soil profile and reveal that fire impact on soil functional C reservoirs is limited to surface soil horizons. The main drivers of CBioR varied with the soil layer considered. The breakdown of organic matter in the O layer was constraint by pH and exchangeable Mn, while that of the mineral soil was dependent on exchangeable Al availability, 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 CBioR,

485
with climate (water supply parameters only, 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 forecast the response of the enormous boreal forest soil carbon pool to global warming.

Data availability
All data presented in this paper can be accessed online free of charge on Canada's Open Government website 490 (https://open.canada.ca, DOI to be determined).

Author contribution
BA, DP, YB and PG participated in writing the funding application and in designing the study. BA and PG did the fieldwork.
BA and DP contributed to the lab work. BA and JB analyzed the data. BA prepared the manuscript with contributions from all co-authors.

Competing of interest
The authors declare that they have no conflict of interest.        Table S1 and text for further details).   Table S2 and text for further details).