Effects of a warmer climate and forest composition on soil carbon cycling, soil organic matter stability and stocks in a humid boreal region

The maintenance of the large soil organic carbon (SOC) stocks of the boreal forest under climate change is a matter of concern. In this study, major soil carbon pools and fluxes were assessed in 22 closed-canopy forests located along an elevation and latitudinal climatic gradient expanding 4 ^∘C in mean annual temperature (MAT) for two important boreal conifer forest stand types: balsam fir (Abies balsamea), a fire avoider, and black spruce (Picea mariana), a fire-tolerant species. SOC stocks were not influenced by a warmer climate or by forest type. However, carbon fluxes, including aboveground litterfall rates, as well as total soil respiration (R_s) and heterotrophic (R_h) and autotrophic soil respiration (R_a), were linearly related to temperature (cumulative degree days >5 ^∘C). The sensitivity of soil organic matter (SOM) degradation to temperature, assessed by comparing Q₁₀ (rate of change for a T increase of 10 ^∘C) of soil respiration and R_s10 (soil respiration rates corrected to 10 ^∘C), did not vary along the temperature gradient, while the proportion of bioreactive carbon and nitrogen showed higher values for balsam fir and for warmer sites. Balsam fir forests showed a greater litterfall rate, a better litter quality (lower C : N ratio) and a higher R_s10 than black spruce ones, suggesting that their soils cycle a larger amount of C and N under a similar climate regime. Altogether, these results suggest that a warmer climate and a balsam fir forest composition induce a more rapid SOC turnover. Contrary to common soil organic matter stabilisation hypotheses, greater litter input rates did not lead to higher total SOC stocks, and a warmer climate did not lead to the depletion of bioreactive soil C and N. Positive effects of warming both on fluxes to and from the soil as well as a potential saturation of stabilised SOC could explain these results which apply to the context of this study: a cold and wet environment and a stable vegetation composition along the temperature gradient. While the entire study area is subject to a humid climate, a negative relationship was found between aridity and SOM stocks in the upper mineral soil layer for black spruce forests, suggesting that water balance is more critical than temperature to maintain SOM stocks.

1 Introduction

Of all the IPCC (Intergovernmental Panel on Climate Change) climate zones, the humid boreal region contains the largest terrestrial carbon (C) stocks (380 Pg). It is also the region where the proportion of carbon contained in the soil is the greatest (Sharlemann et al., 2014). Boreal regions are also experiencing the most intense warming (Hoegh-Guldberg et al., 2018). The consequences for the global climate of losing only a portion of the boreal soil organic matter (SOM) stocks are potentially disastrous as this could create a highly damaging positive feedback loop in which increasing warming would lead to greater soil C losses (DeLuca and Boisvenue, 2012). Recent assessments of global simulation models predicting the fate of soil C in northern regions with global change showed that the model predictions vary greatly between models and that they are prone to large uncertainties (Huntzinger et al., 2020; Wieder et al., 2019). These studies highlighted a great need to better understand the SOM dynamics and its relation to climate and the necessity of empirical observations to guide model development.

The use of climatic gradients to understand the impact of climate change on soil C dynamics complements the information generated by soil warming experiments (Crowther et al., 2016; Rustad et al., 2001). If vegetation types and soil conditions are held constant, climatic gradients allow for the observation of changes in soil fluxes with climate while avoiding short-term pulse effects that have been observed in manipulative warming experiments which may only be transient (Melillo et al., 2002). They can also be used to observe changes that takes long time to form, including the buildup of stable SOM (Lavallee et al., 2020).

Climate gradient studies generally show shorter organic matter turnover times under warmer climate (Tewksbury and Van Miegroet, 2007; Ziegler et al., 2017). However, because both C fluxes to and from the soil are accelerated by an increase in temperature, the net effect on soil C accumulation varies when the effect on one flux is greater than on the other. For example, a Finnish thermocline study showed greater accumulation of soil C in warmer climate (Liski and Westman, 1997) due to increased productivity, while the opposite was found for the SOM content of the soil organic layers of black spruce forests of Alaska (Kane et al., 2005; Kane and Vogel, 2009) and for both organic and mineral soil layers in the continental US (Garten, 2012). Other studies have revealed stable soil organic carbon (SOC) stocks along gradients of increased mean annual temperature in the Southern Appalachian region (Tewksbury and Van Miegroet, 2007) and in Newfoundland and Labrador (Ziegler et al., 2017). In the latter study, increased fluxes to the soil at warmer sites were associated with higher C stocks in the organic layer.

A conceptual model of SOM dynamics, the Microbial Efficiency-Matrix Stabilization (MEMS) framework developed in Cotrufo et al. (2013), suggests that a larger input of labile plant material leads to larger and more stable SOM stocks. Thus, if warming enhances plant productivity, it could favour the accumulation of larger SOM reservoirs. Stable SOM makes up the bulk of SOM. For example, Andrieux et al. (2020) found that only 10 % of total SOM could be considered fast or bioreactive C. To the contrary, several studies that considered the decomposability of SOM have shown a higher content of bioreactive or poorly stabilised soil C in colder soils (Fissore et al., 2009; Ishikuza et al., 2006; Laganière et al., 2015; Norris et al., 2011), suggesting that a cold temperature maintains reservoirs of easily decomposable SOM and consequently that cold soils may be more susceptible to C losses upon warming.

The interpretation of climate gradient studies is often complicated by changes in vegetation composition and in soil properties. The objective of the present study was to assess changes in SOC stocks, quality and C fluxes to and from the soil along a climatic gradient, while maintaining vegetation composition and soil properties as constant as possible for two important boreal forest types: (1) forests dominated by balsam fir (Abies balsamea (L.) Mill.), a fire avoider that becomes abundant in landscapes with a long fire return interval, and (2) forest dominated by black spruce forests (Picea mariana (Mill.) B.S.P.), a fire-adapted species that nevertheless can persist in the landscape in the absence of fire (Couillard et al., 2018). Wildland fires are projected to increase in number, size and intensity, due to climate warming and fire frequency could increase by 1.5 to 4 times before 2100 in the Canadian boreal forest (Boulanger et al., 2014), a situation that would likely lead to a reduction in the area covered by balsam fir forests.

We hypothesised the following:

• i.

  Stocks. Balsam fir sites, presumably because they produce greater amounts of litter of higher quality, and sites at the warmest end of the gradient should accumulate more stable SOM in line with the Microbial Efficiency-Matrix Stabilization (MEMS) framework developed in Cotrufo et al. (2013), which suggests that a larger input of labile plant material leads to larger and more stable SOM stocks. A larger stock of stable SOM should be apparent in total C stocks because it represents the largest share of total SOM (Andrieux et al., 2020).

• ii.

  Cycling. Warmer sites should exhibit enhanced C cycling to and from the soil.

• iii.

  SOM bioreactivity. Warm sites should contain less bioreactive SOM because higher microbial activity at high soil temperature leads to its degradation (Laganière et al., 2015). The proportion of bioreactive SOM can be evaluated with the use of long-term incubation. A lower proportion of bioreactive SOM should also be reflected in a higher Q₁₀ for soil respiration rate according to the “C quality-temperature” (CQT) hypothesis (Laganière et al., 2015). Hence, balsam fir sites and warm sites should have greater soil C stocks and higher Q₁₀.

In comparison to other transect studies dealing with SOC dynamics, this study considered two forest types. Forest composition can have major impacts on SOC storage and cycling (Mayer et al., 2020). In addition, studying the impact of climate on SOC within vegetation type is important in the light of studies pointing out that vegetation shifts with climate change were found to be much more modest than expected. This has been related to non-climatic factors including seed predation, pathogens and soil-related limitations (Brown and Vellend, 2014). In fact, land-use change and disturbances were found to have much greater effects on forest composition than climate change (Danneyrolles et al., 2019). These studies suggest that forest composition shows a strong inertia to climate change and hence that it is timely to study the impact of climate change within a vegetation type.

2 Material and methods

In total, 12 balsam fir and 10 black spruce sites were selected along a climatic gradient caused both by latitude and by elevation (Table 1; Fig. 1). All sites were mature closed-canopy stands with no signs of recent human or natural disturbances. The climatic gradient extended over a 4 ^∘C mean annual temperature difference, which represents mean annual temperature (MAT) increase expected for this region (2011–2080) according to IPCC scenario RCP4.5 (Ouranos, 2015). The study region has a very humid climate with an aridity index comparable to that of the rainforests such as the Pacific coast temperate rainforests of North America or the Amazon (Trabucco and Zomer, 2018). Aridity was assessed using the Penman–Monteith equation calculated daily from May to October over a 30-year period (Table 1). Post hoc, we found out that temperature (degree-days: DD) and aridity were strongly correlated for balsam fir sites only (R²: 0.86; p<0.0001) and that this relationship is not significant for black spruce sites; P=0.26 (Supplement Fig. S1). For balsam fir, warmer sites are drier. This relationship does not hold for our black spruce sites because cold sites included both wet (high elevation) and drier sites (high latitude) (Fig. 1). Sites were at least a few kilometres apart from each other and are the same sites as those described in Larocque et al. (2014), with the addition of three sites north of Lac St-Jean that were added to increase the number of cold black spruce sites. These sites were Lac Tirasse (TIR), RESEF (RES) and eastern old black spruce (EOBS), a flux tower site that was intensively studied during the Fluxnet Canada and the Canadian Carbon program from 2002–2011. Only the closed-canopy section of this latter site (plots FCM 17–24) was considered and not the open lichen woodland portions (FCM 36–39), which would have make it too different from other black spruce sites along the climatic sequence. Four black spruce and three balsam fir sites were intensively studied. These sites are marked in bold in Tables 1 and 2 (FM, FR, LS, J2, TIR, PR, OEBS). On these sites, root trenching was conducted to distinguish the contribution of autotrophic from that of heterotrophic soil respiration, and soil samples were collected for long-term laboratory incubations to determine the proportion of bioreactive soil C. Methods are described below. At each site, five 2 m × 2 m sample plots were located on the perimeter of a 400 m² circular plot where forest tree mensuration measurements were conducted (see Larocque et al., 2014). In addition, on the seven intensive sites, five 2 m × 1 m root trenching plots (Tr) were installed, paired with each control plots. Field fluxes were measured during the snow-free period for 5 years 2002–2006 at a monthly interval during the snow-free period, from early May to late October.

Figure 1 Location of the study sites showing how the temperature gradient is related to elevation and latitude. Data from Abrams et al. (2020).

Table 1 Site characteristics. Intensively studied sites are marked in bold characters.

Note that temperature, precipitation and aridity were estimated daily over a 30-year period using BioSIM (Régnière and Saint-Amant, 2008). Aridity was estimated from May to October. B.A. – basal area; Sa, Si, Cl – sand, silt and clay respectively; Humo-Ferric Podzol (HF-P); Ferro-Humic Podzol (FH-P) or Humic Podzol (H-P); textural class – sandy loam (SL), loamy fine sand (LFS), silty loam (SiL), loam (L) Podzol (H-P). MAT – mean annual temperature; MAP – mean annual precipitation.

2.1 Soil C pools

At each sample plot (five per site), three volumetric soil samples of the organic layer (FH layers) were sampled using a 400 cm² template, and three samples of the top 0–20 cm mineral soil layer were sampled with a 4.70 cm diameter corer. We discarded the loose portion of the top litter because it is short-lasting and varies during the season. It represented a very small portion of the O layer. This generated 15 soil samples per layer and site. In addition, one 20–40 cm volumetric mineral soil sample was taken per sample plot. The exact depth of the sampling was noted in the cases when it was not possible to sample to these depth because of rocks. Samples were kept in a cold room at 3 ^∘C until they were processed. Then they were oven-dried at 70 ^∘C and sieved at 2 mm to remove roots and coarse fragments. They were then weighted, and C content was determined by dry combustion using a LECO CNS-2000 analyser (Leco Corp., St. Joseph, MI, USA).

2.2 Litterfall

Eight square litter traps (2862 cm²) were installed on each site. Litterfall was collected at least twice a year during 3 years. Litterfall per trap was divided by the number of days for the full sampling period and then multiplied by 365 to get annual estimates. Samples were kept at 3 ^∘C, before being sorted to distinguish foliar material from other fine material (small branches, mosses, insects, insect frass, cones, flowers and unknown materials) dried at 70 ^∘C. They were ground on a Wiley mill and analysed for C and nitrogen (N) content as above.

2.3 Soil respiration in situ

In each sample plot, both in control and trenched plots, two 5 cm deep polyvinyl chloride (PVC) collars measuring 10 cm in diameter were inserted into the soil at a 1 m distance from each other (10 collars per site) to measure total soil respiration (R_s). On the seven intensive sites, trenches were dug to a depth of approximately 30 cm on a 2 × 1 m area near each sample plot to sever roots. Trenches were lined with landscaping fabric and backfilled. Herbs and shrubs were removed from trenched plots periodically. Two PVC collars were installed in each plot to measure heterotrophic soil respiration (R_h). Installation was conducted in June, and measurements started the following year to let the severed roots decompose and avoid a method-induced CO₂ peak. CO₂ efflux measurements were performed once a month from early May to mid-late-October. They were measured by placing the soil chamber (LI-6400-09) of a portable infrared gas analyser (LI-6400, LI-COR Inc., Lincoln, NE) above the collar. During soil CO₂ efflux measurements, soil temperature at a depth of 10 to 15 cm was recorded as well as soil moisture (i.e. volumetric water content, % $v/v$) using a TDR-300 time domain reflectometry moisture probe (Spectrum Technologies Inc., Plainfield, IL). This relatively shallow depth was chosen because soil temperature at this depth fluctuates daily along the growing season. Soil temperature was measured continuously on the intensive sites with thermocouple wires attached to a Campbell Scientific data station, while WatchDog 100 Series Water Resistant Button Loggers were installed on the other sites. As detailed in Table 5, 5093 CO₂ efflux soil collar measurements were performed on control plots and 2661 on trenched plot over the duration of the study. Measurements were performed over 4 years for most sites; three sites were sampled for only 3 years and one site for 2 years at a rate of five to six measurement days per year (details in Table 5). Several parameters were derived from the respiration vs temperature relationship. The relationship of soil CO₂ efflux vs temperature was parameterised with a simple exponential equation (Buchmann, 2000):

$$\begin{array}{}\text{(1)}& R_{\mathrm{s}}=a{e}^{(b\cdot (T_{\mathrm{s}}-\mathrm{10})},\end{array}$$

where R_s is soil respiration (in μ moles C m⁻² s⁻¹), a and b are coefficients, and T_s is soil temperature (in ^∘C) at a depth of 10–15 cm. b1 is the estimate of respiration at 10 ^∘C or R_s10. Q₁₀ expresses the rate of change of R_s with a temperature change of 10 ^∘C. It can be calculated from b:

$$\begin{array}{}\text{(2)}& Q_{\mathrm{10}}=e^{\left(\mathrm{10}b\mathrm{2}\right)}.\end{array}$$

The parameters of the equations were estimated for each site using all measurements from all years. We combined all the measurements per site because we noted that the fit was better when the range in soil temperatures used was wider and because we noted the trends were comparable between years.

2.4 Assessing annual soil respiration fluxes

We estimated soil respiration from 1 May to 31 October (183 d) period. In the first measurement period, soil temperature was close to 0 ^∘C, and in the last measurements, it varied from 2 to 4.5 ^∘C depending on local climate. This period is close to the snow-free period. Ecosystem respiration during the snow-covered period has been estimated at about 10 % of total ecosystem respiration for one of our sites (Bergeron et al., 2007). The methodologies outlined in Lavigne et al. (2003) were used to scale CO₂ efflux over the year (183 d) and to derive CO₂ effluxes between measurement dates and for a few days at the end or beginning of the season. A full description of this method and its rationale is presented in Lavigne et al. (2003). This method accounts for seasonal trends in R_s that are not related to temperature change, for example, changes that may be due to a seasonal flush of labile C with fine root development or with fresh litter fall inputs. Briefly, a R_s10 value was calculated for each site and measurement date assuming a Q₁₀ value of 2. This value was chosen because most models use this value (Meyer et al., 2018) and because we have no knowledge if the Q₁₀ value that was estimated for the whole growing season is appropriate for specific periods within the season. A R_s10 value was then linearly interpolated for each calendar day between each measurement day. This daily R_s10 value was used, together with the recorded soil temperature at 10–15 cm depth, to estimate daily soil respiration per site with Eq. (3). This equation converts the units from µmol CO₂ s⁻¹ m⁻² into kg C ha⁻¹ d⁻¹ assuming a Q₁₀ of 2. Daily measurements were summed from 1 May to 31 October to obtain yearly estimates of R_s.

$$\begin{array}{}\text{(3)}& R_{\mathrm{s}}=R_{\mathrm{s}\mathrm{10}}\cdot {\mathrm{2.7183}}^{(\mathrm{0.0693}\cdot (T_{\mathrm{s}}-\mathrm{10}\left)\right)}\cdot \mathrm{3600}\cdot \mathrm{24}\cdot \mathrm{12}\cdot {\mathrm{10}}^{-\mathrm{5}}\end{array}$$

We compared the results from this calculation method with one using Q₁₀ values that vary with site but that were the same for the whole season for the intensively studied sites. The overall difference between the two methods varied between sites and was only 1 % on average in yearly soil respiration (more details provided in Table S1)

2.5 Partitioning into heterotrophic and autotrophic respiration

On the intensive sites, the same calculation as above was performed for the trenched plots to obtain an annual soil efflux per site for the trenched plot. The same parameters were fitted using Eqs. (1) to (3). The parameters derived from these equations were analysed and discussed. However, because the values were higher than those generally reported in the literature (R_h averaged 77 % of R_s) although close to the values reported by Lavigne et al. (2003) for cold sites, we suspected that the trenching was incomplete. We therefore used the often used (Naidu and Bagchi, 2021) relationship of Bond-Lamberty et al. (2004):

$$\begin{array}{}\text{(4)}& \ln \left(R_{\mathrm{h}}\right)=\mathrm{1.22}+\mathrm{0.73}\ln \left(R_{\mathrm{s}}\right),\end{array}$$

where R (h or s) is respiration (in g C m⁻² yr⁻¹). This value was subtracted from R_s to estimate R_a, autotrophic respiration.

2.6 Lab incubation

For the microcosm experiment, soil samples from four of the five plots per site were pooled by depth to yield composite samples for four intensively studied sites (mild fir – FR; war fir – LS; cold spruce – J2; warm spruce – PR; see Table 1 for site codes). Organic and mineral layer samples, 6 and 40 g dry weight, respectively, were placed on a layer of glass wool in 120 mL plastic containers (28 cm² surface area), wetted at field capacity and placed in 500 mL glass jars (Mason type). These microcosms were left to stabilise for 1 week at 2 ^∘C after handling and were then incubated in the laboratory at 3, 10, 15 and 22 ^∘C at field capacity for 354 d. These temperatures were chosen to cover the range of soil temperatures that decomposers would experience in the field. Soil respiration rates were measured monthly during periods of 4 to 24 h during which the lid was closed. The length of these periods depended on rates and CO₂ evolution in order to keep concentrations within the calibration range of the IRGA, a LI-6400 portable photosynthesis system (LI-COR, Lincoln, NE, USA). The methodology used is described in Andrieux et al. (2020). At each measurement period, samples were flushed with distilled water, and NO₃–N and NH₄–N were analysed by FIA (Quickchem 8500, Lachat Instruments, Loveland, Colorado). Flushing maintains the soil humid and prevents the accumulation of metabolic products that may interfere with the decomposition process. Using linear interpolation between measurement dates, the total amount of C and N mineralised during the full incubation was estimated and considered as bioreactive C or N.

2.7 Statistical analysis

The effects of DD, species and their interaction on the dependent variables (soil carbon pools and fluxes as well as soil respiration parameters (R_s10, Q₁₀) and litterfall properties) were tested using PROC GLM (SAS, 2013), DD being a continuous variable and species a categorical one. Sites are considered as the experimental units. For laboratory incubation, the effects of incubation temperature, climate and species were tested using linear regression models fitted separately on each soil layer. Due to the small sample size, interactions between independent variables could not be evaluated. Effects were considered significant at the 0.05 threshold level.

3 Results

3.1 Soil carbon stocks

Soil C stocks varied widely among sites, ranging from 108 to 202 C t ha⁻¹ for fir sites and from 76 to 310 t ha⁻¹ for spruce sites (Table 2). This variability could not be attributed to species, DD or their interaction (Table 2). On balsam fir sites, the organic layer contained an average of 25 % of estimated soil C stocks, while this proportion was 37 % on spruce sites. Excluding the two sites where the 20–40 cm layer could not be sampled, the latter value falls to 33 % and is thus close to that found for the balsam fir sites. No significant effects of tree species, DD and their interaction on soil C stocks were found in either O.L. or the top 0–20 cm mineral soil layer (Table 3; Fig. 2a).

Table 2 Soil carbon stock (t C ha⁻¹) mean and coefficient of variation (CV) for the study sites. O.L. stands for organic layer (FH) and 0–20 and 20–40 cm for the soil depth increments of mineral soil samples. ^* CWD: coarse woody debris. Total C: total of all listed pools. Intensively studied sites are marked in bold characters. NA: not available.

Table 3 Results of statistical analyses (p values) from a GLM testing for forest type (tree species), degree-days (DD) and their interaction for C fluxes, soil C pools and estimated parameters from the in situ soil respiration–soil temperature relationships.

3.2 Carbon fluxes to and from the soil

Aboveground litterfall was linearly related to DD as well as to tree species, but the interaction of these two factors was not significant (Table 3; Fig. 2b, c), indicating that the slope of litterfall with degree-days did not differ between the two forest types. On average, fir forests produced 71 % more foliar litter and 37 % more total litter than black spruce ones. The N concentration of foliar litter was higher for balsam fir than for black spruce, and this is reflected in lower C : N values for the former (Table 3; Fig. 2d). This variable was remarkably stable for balsam fir across the temperature gradient, while it showed greater variability for black spruce but no trend with DD (Fig. 2d).

Cumulative soil respiration from 1 May to 31 October was linearly related to DD and was not significantly affected by tree species (Table 3, Figs. 3, S1). Autotrophic, heterotrophic and total soil respiration rates were strongly and linearly related to degree-days (Fig. 3). The slope of autotrophic respiration with DD was greater than that of heterotrophic respiration (Table S2), although this difference was not significant at the 5 % threshold (p=0.0832), suggesting that a higher share of autotrophic respiration to total soil respiration for warm sites as observed by Lavigne et al. (2003). For two of the coldest sites, the values obtained from trenched plots (indicated by stars in Fig. 3) were similar to the value estimated with the Bond-Lamberty et al. (2004) equation. However, the field values obtained for the warm sites were much higher than the modelled values. This could indicate that the trenching of the warm sites was incomplete.

Figure 2 Relationships between degree-days (DD) and (a) soil carbon stocks (O.L. + top 0–20 cm mineral soil) (b), soil annual foliar litterfall, (c) total annual aboveground litterfall, (d) foliar litter C : N ratio, (e) Q₁₀ and (f) R_s10 derived from field soil respiration measurements. Only the relationships with litterfall (b, c) were significant (see Table 3).

Figure 3 Seasonal soil respiration (183 d: 1 May–31 October) with cumulative degree days for R_s (total soil respiration), R_h (heterotrophic soil respiration) and R_a (autotrophic soil respiration); stars indicates heterotrophic respiration assessed for trenched plots. R_h and R_a are otherwise calculated using the equation of Bond-Lamberty et al. (2004). All slopes are significant at the 0.001 level. More information in Table S2.

3.3 Reactivity of SOM to temperature and tree species

Soil respiration measured in situ was dependent on soil temperature, but the fit of the Q₁₀ relationships varied greatly between sites. R² values ranged from 0.17 to 0.62 and averaged 0.43. The Q₁₀ values ranged from 1.55 to 3.12 and averaged 2.51. R_s10 varied from 2.41 to 3.87 and averaged 3.09 µmol m² s⁻¹ (Table S3). The effect of DD on both Q₁₀ and R_s10 was not significant (Table 3; Fig. 2e, f). This indicates that colder sites did not show a greater reactivity of soil respiration to temperature and that at a similar temperature, cold sites would not show a greater soil respiration than warm soils. R_s10 was significantly affected by tree species, while the interaction with DD was not significant (Table 3). These results suggest that balsam fir soils would show a greater soil respiration rate than black spruce ones under a similar climate. On trenched plots, the values of R_s10 ranged from 1.94 to 3.36 and averaged 2.53 and were, as expected, systematically lower as compared to their respective control plots (Table 3, Fig. 2f: empty symbols). The Q₁₀ values of trenched plots ranged from 1.55 to 3.12 and averaged 2.51 across sites and were not systematically lower nor higher than the values of their respective control plots (Table 3, Fig. 2e).

Laboratory incubation indicated a higher proportion of bioreactive C and N in the organic layer than in the mineral soil (Table S2). Incubation temperature had an overall significant effect for both soil layers and elements (Fig. 4, Table S4). The effects of DD and species on bioreactive C and N were only significant in the OL (Table S4). Warm and balsam fir sites showed consistently higher N and C mineralisation rates (Fig. 4). This effect was significant at the 0.05 threshold for N but barely missed that mark (p=0.06) for the species effect on C mineralisation (Table S4).

4 Discussion

4.1 Soil carbon stocks

No significant changes in SOC stocks along the climate gradient were found, and we refute hypothesis (i). Our result do not give support to the MEMS framework (Cotrufo et al., 2013). A warmer climate increases litterfall inputs as well as soil respiration rates without causing a net effect on C stocks. However, we must be cautious with this interpretation because soil C stocks estimates have a high level of uncertainty, and it is not possible, given our sampling effort, to detect small changes (Yanai et al., 2003). Climate gradient studies in the boreal region have shown that warmer conditions can increase (Liski and Westman, 1997), have no effect (Ziegler et al., 2017) or lower soil carbon stocks (Norris et al., 2011; Kane and Vogel, 2009) despite a higher productivity. These latter two studies were conducted in regions with a much drier climate than the one of this study and of that of the other studies cited above. This suggests that a negative effect of climate warming on SOM stocks appears above a certain threshold of aridity. When heat coincides with drought, which is more likely in drier climates, these events strongly reduce gross primary production (GPP) but yield smaller reduction in ecosystem respiration, leading to strong reduction in net ecosystem productivity (NEP) (von Buttlar et al., 2018).

We tested the relationship between aridity and SOM stocks (Fig. S2). No relationship was found for the organic layer. However, we found a negative relationship between aridity and SOM stocks in the upper mineral soil for all species combined (R²=0.29; p: 0.015) and a stronger relationship for black spruce sites (R²=0.44; p: 0.03) but no significant relationship for balsam fir (Fig. S2). Because aridity is not related to temperature along our black spruce site gradient, these results indicate that aridity may have a more important role to play in maintaining boreal SOC stocks than temperature. The influence of these two factors may often be confounded in climate gradient studies.

Under the cold and wet conditions at our sites, microbial decomposition as well as plant productivity is strongly temperature-limited, and greater litter inputs due to warmer conditions are likely offset by a greater microbial decomposition. In addition, our study area does not appear to suffer from water deficit, and tree growth is not sensitive to summer soil moisture (Girardin et al., 2021). Furthermore, the size of the SOM stock is not only controlled by climate or net primary production (NPP) but is also influenced by soil types and drainage (Andrieux et al., 2018; Dalsgaard et al., 2016). While we tried to maintain these conditions as constant as possible in our site selection, they may still have an important control on soil C stocks overriding that of C fluxes (Dalsgaard et al., 2016). We were not able to explain the large variability in soil C stocks across sites. This property is highly variable at small scale and has notoriously been difficult to map (Paré et al., 2021). A much larger dataset would be required to test for the effect of environmental factors on soil C storage and should include gradients of texture, drainage, vegetation composition and climate including aridity with as much independence as possible between driving variables.

Figure 4 Fitted relationships for organic layer cumulative soil CO₂–C and N mineralisation with incubation temperature, site climate (warm vs. cold) and tree species (balsam fir vs. black spruce). Red colour identifies warm sites and blue colour cold sites, with points showing observed values. Shaded areas indicate 95 % confidence intervals.

4.2 Carbon fluxes in relation to forest type and a warmer climate

Carbon fluxes to and from the soil were all linearly related to the cumulated degree-days in agreement with hypothesis (ii). The litterfall rate was the only parameter that also showed a clear species effect. This may be due to the fact that fir trees do not keep their needles for as long as black spruce, about 5 years (Fleming and Piene, 1992) compared to 5 to 30 years (Grossnickle, 2000). Both species also showed a stable litter C : N ratio along the climate gradient, suggesting that the stoichiometry of C to N is not affected by climate. The difference in litter C : N ratio between species is not surprising because black spruce typically shows lower foliar N concentration than balsam fir (Paré et al., 2013). Despite this difference for N, and perhaps greater N limitation in black spruce forests, the C litterfall rate of both species reacted positively and with the same slope to a warmer climate. Litterfall has been found to be a good predictor of net primary production (Malhi et al., 2011). Our results thus suggest an enhanced tree productivity under warmer conditions. This result contrasts with the results of Larocque et al. (2014) conducted on approximatively the same sites, which did not show any significant differences in radial tree growth along the climate gradient. This lack of effect could be explained by the fact that the forest stands were mature, at a stage where growth is slow. D'Orangeville et al. (2016), using a much larger dataset, did observe a positive effect of warming on tree growth in Quebec for both balsam fir and black spruce. They predicted that enhanced tree growth would be observed for an increase of 4 ^∘C even without a change in precipitation for black spruce but would require a 15 % increase in precipitation for balsam fir. Interestingly, we found a positive relationship between aridity and balsam fir litterfall and total soil respiration for all species (Fig. S2). Although a strong correlation between aridity and DD for balsam fir precludes disentangling the driving factors, the enhanced fluxes with warming suggest that water limitation is not an important driver of tree productivity and soil respiration in this wet climate.

The linear trends that we observed between carbon fluxes and temperature using a relatively small number of sites indicate that the assessment of biogeochemical fluxes, including soil respiration measurements or foliar litterfall may be more suited to study changes in NPP than the often used tree ring analysis. These fluxes are less affected by stand maturity and are more directly connected to carbon capture by photosynthesis than stemwood biomass, which only represents from 10 % to 30 % of photosynthate allocation and depends on many physiological processes, including allocation to below- and aboveground parts, reproduction, defence and energy storage (Litton et al., 2007). Our results showed clear and linear trends with warming for litterfall and soil respiration, while tree ring analysis has shown divergent results in the Canadian boreal forest (e.g. Girardin et al., 2016 vs. Hember et al., 2019).

Increased litterfall rate with a warmer climate was accompanied by a linear and positive response of soil respiration, indicating a greater turnover rate of SOC as the climate warms. These results are coherent with the observations of Ziegler et al. (2017) for a climatic gradient in balsam fir forests in Newfoundland and Labrador and suggest that soil C stocks can be maintained despite a greater turnover rate. These results are in line with an acceleration of C turnover rate observed at the global scale over the recent decades (Carvalhais et al., 2014). Interestingly, the maintenance of SOC pools, despite a faster turnover, has also been observed for a tropical mountain climatic gradient (Giardina et al., 2014), suggesting that within a biome where water limitation to forest production is of minor importance, the same pattern of accelerated C cycling and maintenance of pools may occur.

4.3 Bioreactive SOM in relation to a warmer climate

Several studies have concluded that a warmer climate would lead to a depletion of bioreactive or poorly stabilised SOC (Fissore et al., 2009; Ishikuza et al., 2006; Laganière et al., 2015; Norris et al., 2011). Our results indicated, to the contrary, that the bioavailability of SOC is either similar or higher with a warmer climate. Field results revealed that Q₁₀ and R_s10 were stable along the climate gradient, while results from the laboratory incubation indicated a greater proportion bioreactive C and N for balsam fir and for warmer climate but only for the O layer. Both results support the rejection of hypothesis (iii). The stability of R_s10 along the climate gradient suggests that bioreactive SOM does not accumulate more in cold than in warm soils. The higher R_s10 values found for balsam fir soils, irrespectively of climate, are congruent with the greater litterfall rates observed for this species in this study.

According to the C quality-temperature (CQT) hypothesis, decomposition of organic matter (OM) exhibiting a low decomposition rate should have a greater temperature sensitivity than that of faster-turnover substrates (see Laganière et al., 2015), while we observed no change in Q₁₀ with climate. Our results diverged from those of Laganière (2015) in that we did not find a lower reactivity to temperature, nor a greater amount of bioreactive OM in the soils of the cold sites. The increase in bioreactive organic matter on cold sites, observed in this latter study, could be due to a greater abundance of bryophytes in the colder sites of this transect study, favouring the accumulation of bioreactive organic matter (Kohl et al., 2018). We did not observe a trend for a greater abundance of bryophytes as a ground cover along our climatic gradient, and this could explain the absence of a change in soil bioreactive OM in our study. Bryophyte abundance is linked to canopy openness (Pacé et al., 2016), and our sites all had a closed canopy. The impact of the development of a bryophyte cover on soil C dynamics can be important (Simard et al., 2007).

A reduction of bioreactive soil C with warming had been observed in a manipulative experiment on a nearby site. D'Orangeville et al. (2013) as well as Marty et al. (2019) found a reduction in labile C with artificial soil warming, while they observed no change in SOM composition congruently with our observations. The difference in the results from these studies and ours regarding changes in soil respiration reactivity to temperature (Q₁₀) and bioreactive C suggests that artificial soil warming may not have the full effect that long-term climate warming has, both on fluxes to and from the soil. A balance between C inputs and outputs may contribute to the maintenance of soil organic matter stocks and properties despite warming. Climate gradient studies may also be more suited than warming experiments to observe changes that take a long time to appear.

Stabilised MAOM (mineral-associated organic matter) reservoirs may reach a saturation point, especially in non-recently disturbed soils (Lavallee et al., 2020). Boreal regions show the highest concentrations of dissolved organic carbon (DOC) in the surface soil globally (Langeveld et al., 2020), indicating that the capacity of these soils to immobilised DOC as water percolates through the soil column is limited. Cotrufo et al. (2021) suggested that under cold and wet conditions, it is not MAOM but poorly stabilised particulate organic matter (POM) that dominates the dynamics of SOC cycling. If indeed MAOM reservoirs have reached saturation, and POM dominates the SOC cycling, our results suggest that warming, while accelerating SOC cycling, does not lead to changes in the stocks of either POM or MAOM. However, more research is needed to determine how the different fractions of SOM are impacted by changes both in aridity and in temperature and to identify thresholds in aridity from which SOC stocks become vulnerable to warming.

5 Conclusions

Our results showed no evidence of net SOM losses or a reduction of the most active SOM fraction with a warmer climate and contrast with the conclusions that were reached in other studies derived from climate gradient studies or from direct soil warming experiments. Our results are applicable to the context of the study: a cold and wet climate where available moisture does not limit productivity or decomposition, closed-canopy coniferous forests that were not recently disturbed and the absence of change in forest composition along the climate gradient. Other studies conducted under a wet climate corroborate these results. Altogether, results from this study and from others suggest that a certain aridity threshold needs to be reached before a warmer climate has implications for SOM stocks. Recognising that the effects of climate change on SOM are not linear, identifying these thresholds would greatly improve our capacity to predict the fate of the important SOM pool in boreal ecosystems under climate change. Our results also suggest that aridity may be more important than temperature in controlling SOM stocks.

Forest composition had effects on both SOM quality and the rate of C cycling. However, warming did not show interacting effects with forest composition. This suggests that differences in C cycling with forest composition remained stable with projected changes to temperature. These results indicate that the impact of climate change on SOM storage and dynamics needs to be studied both within and among forest ecosystem types to separate the direct effects of climate change from that of vegetation change.