Implementation of mycorrhizal mechanisms into soil carbon model improves the prediction of long-term processes of plant litter decomposition

Ecosystems have distinct soil carbon dynamics, including litter decomposition, depending on whether they are dominated by plants featuring ectomycorrhizae (EM) or arbuscular mycorrhizae (AM). However, current soil carbon models treat mycorrhizal impacts on the processes of soil carbon transformation as a black box.

We re-formulated the soil carbon model Yasso15 and incorporated impacts of mycorrhizal vegetation on topsoil carbon pools of different recalcitrance. We examined alternative conceptualizations of mycorrhizal impacts on transformations of labile and stable carbon and quantitatively assessed the performance of the selected optimal model in terms of the long-term fate of plant litter 10 years following litter input.

We found that mycorrhizal impacts on labile carbon pools are distinct from those on recalcitrant pools. Plant litter of the same chemical composition decomposes slower when exposed to EM-dominated ecosystems compared to AM-dominated ones, and across time, EM-dominated ecosystems accumulate more recalcitrant residues of non-decomposed litter. Overall, adding our mycorrhizal module into the Yasso model improved the accuracy of the temporal dynamics of carbon sequestration predictions.

Our results suggest that mycorrhizal impacts on litter decomposition are underpinned by distinct decomposition pathways in AM- and EM-dominated ecosystems. A sensitivity analysis of litter decomposition to climate and mycorrhizal factors indicated that ignoring the mycorrhizal impact on decomposition leads to an overestimation of climate impacts on decomposition dynamics. Our new model provides a benchmark for quantitative modelling of microbial impacts on soil carbon dynamics. It helps to determine the relative importance of mycorrhizal associations and climate on litter decomposition rate and reduces the uncertainties in estimating soil carbon sequestration.

1 Introduction

Long-term soil carbon (C) sequestration is to a large extent determined by complex soil–plant rhizosphere and microbial interactions (Dijkstra and Cheng, 2007; Fernandez and Kennedy, 2016; Fontaine et al., 2007; Ostle et al., 2009). These interactions contribute to the atmospheric CO₂ balance (Ostle et al., 2009; Todd-Brown et al., 2012) and are increasingly recognized as processes that counteract climate change (Terrer et al., 2016). Plant associations with fungi, so-called mycorrhizas, are the most widespread symbiosis on Earth, featured by the majority of vascular plants including trees, shrubs, and herbs (Brundrett and Tedersoo, 2018). Mycorrhizae are hypothesized to play especially important roles in soil C sequestration, yet the actual mechanisms of mycorrhizal impacts on soil C dynamics are poorly understood.

Mycorrhizal fungi themselves are not capable of meaningfully obtaining carbon from decomposing plant litter (Bödeker et al., 2016; Lindahl and Tunlid, 2015). Instead, they receive carbon from their symbiotic host plants. However, the relation between mycorrhizal fungi and soil C dynamics is enabled through three potential pathways that likely complement each other (Frey, 2019): (i) provisioning of substrate for decomposition (Leake et al., 2004; Soudzilovskaia et al., 2015), (ii) mediating plant litter quality and amounts (Averill et al., 2019; Cornelissen et al., 2001; Phillips et al., 2013), and (iii) controlling the environment of plant litter decomposition, including mediation of the microbial community (Fernandez and Kennedy, 2016; Frey, 2019). Plant litter decomposition is an important component of soil C cycling and is affected by its chemical composition (Berg and McClaugherty, 2008; Cornelissen et al., 2007), which is generally grouped as labile and recalcitrant components. The fate of carbon originating from components of various chemical recalcitrance will ultimately determine the decomposition and sequestration dynamics (Aponte et al., 2012; Cusack et al., 2009; Kalbitz et al., 2003; McClaugherty et al., 1985). Among the three major pathways of mycorrhizal impacts on soil C dynamics, the pathway of mycorrhizal fungal control on this decomposition environment and the fate of carbon is arguably understood the least.

To understand mycorrhizal fungal impacts on soil C dynamics, we need to distinguish between arbuscular mycorrhiza (AM) and ectomycorrhiza (EM) types of symbiosis. Together, these types are possessed by over 80 % of plant species compromising the majority of terrestrial plant biomass (Brundrett and Tedersoo, 2018; Soudzilovskaia et al., 2019). While they are present in almost all ecosystems, it has been proposed that distinct mycorrhizal types are associated with specific ecosystems and soil attributes (Craig et al., 2018; Read and Perez-Moreno, 2003; Steidinger et al., 2019). Moreover, distinct mycorrhizal guilds differ in the pathways through which they affect the decomposition environment of plant litter. AM fungi (AMF) have limited or no ability to depolymerize organic macromolecules. They do not possess enzymes enabling nitrogen extraction and uptake from soil organic matter (Orwin et al., 2011; Treseder et al., 2016; Treseder and Allen, 2002) but primarily acquire inorganic nutrients mobilized by saprotrophic fungi and bacteria. Accordingly, plant litter subjected to AM fungi-dominated decomposition environment is likely to undergo a more balanced decomposition process, with both labile and recalcitrant components being degraded by saprotrophic decomposers. On the other hand, compared to AM fungi, most EM fungi (EMF) can produce enzymes involved in decomposing organic compounds of plant litter (Fernandez and Kennedy, 2015; Lindahl and Tunlid, 2015; Zak et al., 2019) and therefore have easier access to organic nutrients, especially so to nitrogen. It has been proposed that EMF increase the recalcitrance of decomposing litter, as their ability to uptake nitrogen while withholding carbon compounds from breaking down increases carbon-to-nitrogen ratios of decomposing plant litter (Nicolás et al., 2019; Read and Perez-Moreno, 2003; Rineau et al., 2013). This process of gradually increasing recalcitrance of plant litter subjected to an EM-dominated decomposition environment is further magnified by the suppression of saprotrophic decomposer activities, an effect known as the Gadgil effect (Fernandez and Kennedy, 2015; Gadgil and Gadgil, 1971; Smith and Wan, 2019). Yet the magnitude of the impacts induced by the differential roles of mycorrhizal types on the dynamics of decomposing plant litter is understood very poorly, especially so in quantitative terms.

Traditional field experiments are typically too short to assess the full complexity of the mechanisms underpinning the potential difference of AM and EM impacts on plant litter decomposition processes over time. In addition, traditional field experiments have limitations in explicitly distinguishing the individual pathways of mycorrhizal impacts on the decomposition process, including the fate of litter fractions of different chemical recalcitrance. An alternative tool to progress in our understanding of mycorrhizal impacts on plant litter decomposition is testing different formulations of mycorrhizal impacts in process-based models of litter decomposition and examining how well the models fit the observations.

Current deterministic models of soil C decomposition (e.g. CENTURY, DAYCENT, DAISY, DNDC, NCSOIL, RothC, and Struc-C) do not explicitly account for mycorrhizas as a driver of plant litter decomposition processes. Instead, climate and litter quality, the well-acknowledged regulators of soil organic carbon (SOC) and litter decomposition (Cornwell et al., 2008; Coûteaux et al., 1998; Cusack et al., 2009; Parton et al., 2007; Zhang et al., 2008), are being modelled as primary drivers of all aspects of SOC dynamics. A body of recent studies have questioned the recognition of climate and litter quality as the only dominant regulators in SOC and litter decomposition (Bradford et al., 2016; García-Palacios et al., 2013; Wall et al., 2008) and plead for explicit inclusion of microbial and especially mycorrhizal impacts (Johnson et al., 2006; Shi et al., 2016) on SOC dynamics in biogeochemical models (Clemmensen et al., 2013; Craig et al., 2018; Todd-Brown et al., 2012; Wieder et al., 2013). However, so far, models assessing the role of mycorrhizas in SOC dynamics (e.g. Liang et al., 2017; Orwin et al., 2011; Shi et al., 2016) do not compare the relative impacts of mycorrhiza vs. climate on litter decomposition processes.

In this study, we aim to develop a framework allowing incorporation of mycorrhizal impacts on the decomposition of plant litter into a generic soil C model, specifically addressing one of the most poorly understood mechanisms of mycorrhizal impact on plant litter decomposition – the impact through controlling decomposition environment, separately from climate and other factors. Hereto we focus on answering the following four questions.

• What is the best conceptualization, and accordingly the best representation, in a soil C dynamics model to describe mycorrhizal impacts on the decomposition of plant litter labile and recalcitrant carbon compounds?

• To what extent does modelling mycorrhiza-associated impacts on the litter decomposition environment improve model performance, in terms of model errors, robustness, and temporal dynamics?

• What is the sensitivity of model predictions to the uncertainty of parameters and input describing the pathways of decomposition as affected by mycorrhiza vs. climate and other factors?

• How is the temporal dynamics of plant litter decomposition affected in AMF- vs. EMF-dominated decomposition environments in terms of both total C loss and loss of C from compounds of distinct recalcitrance?

2 Methods

Among available models of plant litter decomposition, the Yasso model (Tuomi et al., 2011a) provides an ideal framework for a mechanistic integration of mycorrhizal impacts into the modelling of plant litter decomposition processes. Yasso is among the models that underpin IPCC predictions of impacts of environmental change scenarios on global C cycles (IPCC, 2006; IPCC, 2019). In the Yasso model, plant litter is classified into five pools, characterized based on measurable chemical solubility of organic matter (Liski et al., 2005): compounds soluble in water (denoted with W); carbon compounds hydrolysable in acid (A); components soluble in a non-polar solvent, e.g. ethanol or dichloromethane (E); compounds neither soluble nor hydrolysable (N); and humus (H) (Berg and Agren, 1984; Palosuo et al., 2005). The W, A, and E pools together form the group of labile C fractions of soil organic matter, N is a recalcitrant but yet not mineral-bound C fraction, and the H pool represents the fraction of very stable soil C that remains in the soil for decades or centuries.

Yasso presents the litter decomposition process as a system of linear differential equations, and the total amount of carbon released from each pool is the result of fluxes between pools and C released to the atmosphere as carbon dioxide. Figure 1 presents the schematic representation of the Yasso model, with carbon flows quantified from results of the original Yasso model formulation (Tuomi et al., 2011a; Viskari et al., 2020). H pool-related flows are not specified in this figure because humus can only be produced in deeper soil accessible to mineral compounds; thus it is not considered in this study of 10 years of litter decomposition simulations. Detailed descriptions of the original Yasso model and the dataset used for its parametrization are provided in Appendices A and B, respectively.

The conceptualization of litter decomposition as a process of C conversion into pools representing measurable C fractions makes Yasso a particularly suitable model for incorporating new (in our case, mycorrhizal) pathways that are based on or affected by differences in litter decomposability.

Figure 1 Conceptual diagram of decomposition and mass flows between five carbon pools in Yasso. In the original Yasso model (Tuomi et al., 2011a, b), the fate of organic matter entering soil as plant litter material is represented as a series of carbon fluxes between carbon pools characterized by distinct decomposability (i.e chemical solubility) levels. Values in arrows show the percentage of C transformed between pools and leaving the pools per yearly time step (% yr⁻¹) according to the original Yasso formulation and parameterizations (Tuomi, et al., 2011a; Viskari et al., 2020). Small flows are dotted lines.

2.1 Implementation of mycorrhizal impacts on decomposition in Yasso: general principles and data

We modified the Yasso model by adding mycorrhiza as a factor controlling the plant litter decomposition environment. Our model focuses on explaining the fate of aboveground plant litter that is decomposed at the topsoil layer before entering into deeper mineral soil or subsoil. During this stage, decomposers pre-process plant litter, liberating carbon compounds which – in a later stage – contribute to the accumulation of mineral-associated organic matter (MAOM) and particulate organic matter (POM) through different pathways in deeper soils (Bradford et al., 2016; Cotrufo et al., 2013, 2015, 2019; Sokol et al., 2019).

We conceptualized the mycorrhizal environment as a driver of soil organic carbon dynamics in addition to the drivers already accounted for by Yasso (i.e. temperature, soil moisture, and litter chemical composition). We modelled impacts of the mycorrhizal environment on plant litter decomposition as the sum of impacts caused by the predominance of AM and EM fungal types. The AM and EM fungal impacts were assumed to depend on the fungal-type-specific ability to affect the litter decomposition process and its biomass. As there are no data currently available about the global distribution of mycorrhizal fungal biomass, we approximated the AM and EM fungal biomass to be proportional to AM and EM plant biomass (the latter estimated as the product of the proportion of AM and EM plant biomass and the total vegetation gross primary production (GPP, using MODIS product-MOD17 data) (Running et al., 2004; Zhao et al., 2005)).

We calibrated our new model using litter decomposition databases (Appendix B) used in Yasso modelling that included total mass loss and the dynamics of different chemical components over time (Tuomi et al., 2009, 2011a, b): CIDET with the measurements from Canada (Trofymow, 1998), LIDET with data from the USA and Central America (Gholz et al., 2000), and Eurodeco (ED) with data gathered from several European research projects (Berg et al., 1991). Chemical composition data consist of the initial composition of litter in terms of W, A, E, and N fractions which were measured for each site. These data, together with other environmental data, were used for initializing the model. In addition, for the ED dataset, W, A, E, and N components had been determined during the decomposition process and at the end of the decomposition. In addition, all datasets were supplemented with site-specific estimates on the fractions of AM and EM vegetation within total plant biomass, which was extracted from the global mycorrhizal distribution map of Soudzilovskaia et al. (2019). To avoid potential mismatches between the actual fractions of AM and EM plants within the total plant biomass and the (generalized) data of AM and EM fractions derived from the map of Soudzilovskaia et al. (2019), the ecosystem type of each site was carefully checked for consistency with the map.

2.2 Mycorrhizal impact on total decomposition

Figure 2 shows the general principle to include the impacts of the mycorrhizal environment on each decomposition pool. The total carbon outflux of each W, A, E, and N pool is controlled by two factors: climate (as in the original YASSO model) and mycorrhizal decomposition environment (the new factor added to the model). See Appendix A for details on the decomposition terms used in YASSO.

Figure 2 Carbon fluxes from and to each X pool of carbon, with X being W, A, E, or N, as represented by the modified Yasso model. The blue arrow and blue box show conceptualization of the added impact of the mycorrhizal environment on the litter decomposition process. While in the original version of the Yasso model plant litter decomposition process was represented as a function of climate and litter quality, in our model decomposition is a function of proportions of ectomycorrhizal and arbuscular mycorrhizal plants in vegetation, climate, and plant litter quality.

Accordingly, we modified the original form of the equations describing the decomposition rate of each W, A, E, and N element in the Yasso model. In the original model decomposition matrix (see Appendix A), only the climate was considered as a driver of decomposition K_i, where $i\in \mathit{\{}W,A,E,N\mathit{\}}$ through K_i (C) (see Eq. A3 in Appendix A). In the new model formulation we added a term M_i representing the mycorrhizal impact on the total C outflux of each W, A, E, and N pool by Eq. (1):

$$\begin{array}{}\text{(1)}& K_i=K_i\left(\mathrm{C}\right){}^{\prime }\cdot \left(\mathrm{1}+M_i\right),i\in \mathit{\left\{}WAEN\mathit{\right\}}.\end{array}$$

The M_i term is described by Eq. (2):

$$\begin{array}{}\text{(2)}& M_i=m_{i\mathrm{AM}}\cdot {\mathit{\lambda }}_{\mathrm{AM}}\cdot G_{\mathrm{pp}}+m_{i\mathrm{EM}}\cdot {\mathit{\lambda }}_{\mathrm{EM}}\cdot G_{\mathrm{pp}},i\in \mathit{\left\{}WAEN\mathit{\right\}},\end{array}$$

where m_iAM and m_iEM are the impacts of AM and EM mycorrhizas on C loss from pool i, λ_AM and λ_EM are the fractions of AM and EM vegetation within the total vegetation biomass, and G_pp is the gross primary production of mycorrhizal vegetation.

We compared four different conceptualizations of AM and EM impacts on distinct W, A, E, and N pools of decomposing litter, by evaluating the performance of four distinct model versions (Fig.3).

• Myco-Yasso.v1 is a model where the magnitude of mycorrhizal impact on carbon loss from each of the W, A, E, and N pools differs among the pools (Fig. 3a), based on the assumption that mycorrhiza impact each pool differently.

• Myco-Yasso.v2 is a model where mycorrhizal impacts on carbon loss from labile soil C pools (W, E, and A) are equal among the pools, while the mycorrhizal impact on carbon loss from the recalcitrant soil C pool (N) differs from the impact on C losses from labile pools (Fig. 3b), reflecting previous findings of the Yasso model that climate factors have similar impacts on W, A, and E pools but are different for the N pool (Tuomi et al., 2009; Viskari et al., 2020). We assume that mycorrhizal impacts are similarly differentiated.

• Myco-Yasso.v3 is a model where mycorrhizal impacts on carbon loss are equal for all pools (Fig. 3c), based on the assumption that the impact is the same for all pools.

• Myco-Yasso.v4 is a model where mycorrhiza affects only carbon loss from the recalcitrant soil C pool (N) (Fig. 3d), based on the assumption that mycorrhiza can only affect the most recalcitrant pool.

Figure 3 Four conceptualizations of the possible mechanisms of mycorrhizal impacts on litter decomposition, modelled with four versions of the Myco-Yasso model. (a) Myco-Yasso.V1: mycorrhizal impacts differ for each of the W, A, E, and N pools. (b) Myco-Yasso.V2: mycorrhizal impacts on W, A, and E pools have the same magnitude, but the mycorrhizal impact on N pool is different. (c) Myco-Yasso.V3: mycorrhizal impacts on W, A, and E pools and N pools are equal. (d) Myco-Yasso.V4: mycorrhizae impact only for the N pool.

We used a Bayesian framework and a differential evolution Markov Chain with a snooker updater (DEzs; Braak and Vrugt, 2008) algorithm and Markov chain Monte Carlo (MCMC) (Haario et al., 2001) for calibrating all relevant parameters following the original Yasso framework (Tuomi et al., 2011a; Viskari et al., 2020). Essential parameters from the original Yasso and newly derived mycorrhizal dependencies with corresponding symbols and units are explained in Table 1. We allowed m_iAM and m_iEM to vary from negative to positive values. The only control on priors of m_iAM and m_iEM is limiting M_i>−1 in Eq. (1) to make the algorithm meaningful. The other parameter priors were adopted from previous Yasso research (Tuomi et al., 2009). We performed cross-validation for each model, using 80 % of the decomposition time series randomly drawn from the dataset for calibration and the remaining 20 % of the decomposition time series for validation. After parameterization, all model versions were examined for Pearson's r and RMSE values of the correlation between the predicted and observed data for both the validation dataset and the full dataset. To account for the fact that the data in the different datasets varied in measurement uncertainty and the number of observations, we opted to compare the Pearson r and RMSE values of models separately for each dataset. We use the root mean square error (RMSE) from the 20 % validation dataset and Akaike information criterion (AIC) and Bayesian information criterion (BIC) based on the 80 % data used for calibration as the criteria for comparing the relative quality of the models. The conceptualization with the lowest RMSE, AIC, and BIC was selected as the optimal model with the best performance.

Table 1 Parameters calibrated for each model version.

2.3 Performance of the selected best mycorrhizal model of soil C sequestration

2.3.1 Model residuals and uncertainty analysis

We examined the residuals (differences between measurements and model-predicted litter decomposition) as a function of AM and EM fractions in the biomass of mycorrhizal vegetation. In addition, the uncertainty of the selected Myco-Yasso model was assessed in two aspects.

• a.

  Variability in estimating total C mass loss through litter decomposition. The variability in the percentage of C mass remaining after 10 years of litter decomposition, as revealed by the original Yasso model and the selected Myco-Yasso model was examined by conducting Monte Carlo simulations for a hypothetic site. In line with previous sensitivity tests of Yasso (Liski et al., 2005), we chose the following input data to represent the conditions of decomposition: mean annual temperature of 5.2 ^∘C and annual precipitation of 840 mm. For the Myco-Yasso model, the mycorrhizal impact in Eq. (2) was quantified by assuming an AM mycorrhizal plant biomass proportion of 38 %, EM mycorrhizal plant biomass proportion of 36 %, and a GPP of 1516 g m⁻² yr⁻¹. We used the following global mean values for the chemical composition of litter: W fraction −20.6 %, A fraction −43.0 %, E fraction −8.7 %, and N fraction −27.7 %. We ran 1000 simulations using parameter values randomly selected from an even distribution of the input parameters within their uncertainty ranges.

• b.

  Sensitivity to parameters and input. With environmental conditions and chemical composition of the litter being the same as used in part (a), we evaluated the sensitivity of litter decomposition by separately increasing model parameters by 1 % and input values by 1 % of variations across the dataset in 10 years of model runs. This test was conducted for both the original Yasso model and the selected Myco-Yasso model.

2.3.2 Temporal dynamics of the model

• 1.

  Model performance over time.

  We examined the ability of the selected model in predicting litter decomposition at different times following litter input, by comparing model predictions of C mass remaining to real measurements of the remaining C at the same time moment. The time slots were classified according to the datasets' characteristics, as litter decomposition measurements for different datasets were taken at different months in a year.

• 2.

  Mycorrhizal impact on labile and recalcitrant litter pools.

  To analyse total litter decomposition and the different litter pools, simulations of 10 years of litter decomposition were conducted in different mycorrhizal environments with varying AM : EM vegetation biomass proportions. Input values in terms of environmental factors and chemical fractions were set consistent with the standard conditions as used in the sensitivity analysis. To evaluate model performance consistency, the analysis was done using chemical fractions of typical root and leaf litter as contrasting litter types (Appendix D).

3 Results

3.1 Model comparison and selection

For all four model versions examined, the calibration based on all three decomposition datasets showed a high correlation between measurements and model predictions, with Pearson's r being 0.84–0.86 for CIDET, 0.67–0.68 for LIDET, and highest with 0.90–0.91 for ED, which also contained information on carbon pools through time. Small differences occurred between individual versions of Myco-Yasso models. However, the model RMSE comparisons revealed that the Myco-Yasso V2 provided the strongest RMSE decrease among all Myco-Yasso models compared to the original Yasso15 model. This pattern was consistent through all datasets (Table 2). The AIC and BIC confirmed that the Myco-Yasso V2 has the best performance. Based on the RMSE, AIC, and BIC, we selected Myco-Yasso V2 as a model representing the optimal conceptualization of mycorrhizal impact on plant litter decomposition. In this model V2, the mycorrhizal impact is similar among labile C compounds (WAE) but different for the recalcitrant C compound (N). Hereafter, this optimal model is referred to as Myco-Yasso. Scatterplots showing model improvement in terms of observed vs. predicted values for the Myco-Yasso model compared to the original Yasso15 model are provided in Appendix C. Details of parametrization outcomes of the Myco-Yasso model are provided in Table C1.

Table 2 Model performance based on RMSE for the 20 % validation dataset and AIC and BIC for the 80 % data used for calibration. Model predictions are based on the total mass remaining in plant litter of different mycorrhizal model versions for different litter decomposition datasets. The selected model and its performance indexes are in bold.

3.2 Model performance across the range of mycorrhizal plant biomass fractions in vegetation

The standardized residuals for the litter decomposition measurements (percent of C decomposed from initial plant litter) as a function of AM and EM fractions in the biomass of mycorrhizal vegetation are shown in Fig. 4. Within the 95 % probability density covered by 2σ intervals, model predictions agreed well with measurements across the entire range of fractions of biomass AM and EM plants in vegetation. The model had relatively large negative residuals at low values of the AM fractions (AM <10 %) and high values of EM fractions (EM >85 %) but relatively large positive residuals at low values of EM fractions (EM <10 %), which suggest a lower predictive power for these conditions.

Figure 4 Standardized residuals (predictions – measurements, P–M) of decomposition, expressed as percent mass loss, modelled by Myco-Yasso as a function of the abundance of AM mycorrhizal plants (a) and EM mycorrhizal plants (b) in vegetation. The circle in the middle of each line is the mean value of the residuals. The intervals contain residuals within 95 % probability intervals.

3.3 Variability in litter decomposition estimations

The 1000 simulations of the Yasso15 model ran for the conditions of a hypothetical site with prescribed environmental conditions revealed a normally distributed dataset with μ = 22.56 % and σ = 1.81 % mass remaining. The same simulations conducted by the Myco-Yasso model yielded a dataset with a lower μ (16.90 %) and lower σ (1.19 %), indicating a lower total sensitivity of the Myco-Yasso model to variation in input parameters. The best-fit normal distributions of these two model predictions are shown in Fig. 5.

Figure 5 The probability density distribution of litter mass remaining as predicted by Yasso15 compared to Myco_Yasso. Compared to predictions of Yasso15, Myco-Yasso reduces the variation in the predictions of C mass remaining from decomposing litter after 10 years of decomposition. The probability density is based on 1000 model runs for conditions of a hypothetical site with the prescribed environmental conditions (see descriptions in Sect. 2.3.1).

3.4 Sensitivity of litter decomposition to parameters and input values

The sensitivity of the Myco-Yasso model to the individual litter decomposition parameters is shown in Fig. 6. The Myco-Yasso model showed the highest sensitivity to the impact of arbuscular mycorrhizal vegetation of the N pool (m_{N_AM}) out of the four mycorrhizal impact parameters (Fig. 6). This implies that an AM environment has a much stronger stimulating effect on the decomposition of the recalcitrant pool compared to an EM environment. In contrast, the decomposition from the labile pools was only a bit more stimulated by an EM environment than by an AM environment. Concerning the decomposition rate parameters, the overall carbon loss in the Myco-Yasso model has a considerably lower sensitivity to the total decomposition rate of the N pool (α_N) and a slightly increased sensitivity to the decomposition rate of the A pool (α_A) compared to Yasso15. However, the total impact of all α terms together on carbon loss predictions is generally similar in Myco-Yasso compared to Yasso15 (Fig. C3).

Figure 6 Model sensitivity to 1 % increase in individual litter decomposition parameters.

Analysis of the environmental dependencies of the Myco-Yasso (Fig. 6) revealed that the new model is less sensitive to the overall variability in temperature parameters (b₁, b₂, b_N1, and b_N2) than the original Yasso15 model, although the overall effect of temperature sensitivity decrease (Fig. C3) is mostly driven by the decreased sensitivity of N pools to temperature (b_N1 and b_N2). The sensitivity to precipitation parameters (g and g_N) of Myco-Yasso is generally similar to Yasso15, with a slight increase in sensitivity of the W, A, and E pools to precipitation (g parameter) and a slight decrease in sensitivity of the N pool to precipitation (g_N). The resulting impacts on carbon transformation in the respective pools are provided in Fig. C4.

Figure 7 shows the model sensitivity to an increase in each input parameter by 1 % including effects of initial plant litter chemistry, climate parameters, and the mycorrhizal environment of decomposition. With an increase in biomass of mycorrhizal plants, less carbon will remain of the decomposing plant litter after 10 years of decomposition. This impact is similar in magnitude to the impact of temperature increase. An increase in EM dominance leads to a slight increase in carbon accumulation, while AM dominance speeds up decomposition to a similar extent as temperature increases. The Myco-Yasso shows a slight decrease in sensitivity to climate variables compared to Yasso15, confirming our supposition that potential mycorrhizal impacts were partly accounted for by climate variables in the original Yasso15. The magnitude of sensitivity of plant litter decomposition to the mycorrhizal environment is comparable to the sensitivity to climate (Figs. 6, 7, and C3).

Figure 7 Model sensitivity to 1 % increase in model input values. Impacts of input parameters are shown in terms of the relative change in total C remaining after 10 years of decomposition. The “AM” bar shows the impact of an increase in AM plant biomass by 1 %, while EM plant biomass remains unchanged. The “EM” bar shows the impact of an increase in EM plant biomass by 1 %, while AM plant biomass remains unchanged. The “Biomass of mycorrhizal plants” bar shows the impact of an increase in the combined biomass of AM and EM plants by 1 % while the AM and EM distribution within the vegetation remains unchanged.

3.5 Model predictions of temporal dynamics of plant litter decomposition

3.5.1 Model performance over time

Figure 8 illustrates how the model predictions of the Myco-Yasso improve the modelled decomposition over time compared to the original Yasso15 model using the full dataset. The differences in models' prediction accuracy (RMSE of the Yasso15 predictions minus RMSE of Myco-Yasso predictions) has a trend of increment over time, indicating an increasing impact of mycorrhizas on litter decomposition dynamics across 10 years. The examination with only validation comparing the original model Yasso15 and Myco-Yasso is provided in Appendix C, Fig. C5.

Figure 8 Improvement of performance comparing the Myco-Yasso model to the original Yasso model over the decomposition period. Bars represent the relative RMSE differences between Yasso15 and Myco-Yasso per period. The line with dots shows the absolute value of the RMSE differences (Yasso15- Myco). Predictions examined from the full dataset simulation.

3.5.2 Mycorrhizal impact on labile and recalcitrant pools

Assessments of the dynamics of total litter mass decomposition under the dominance of AM and EM vegetation with Myco-Yasso (Fig. 9a) revealed that, at the 10th year of decomposition, plant litter (with equal initial chemical composition) will have ca. 15 % less carbon remaining if decomposed in an AM-dominated environment compared to an EM-dominated environment. During the first decomposition year, litter subjected to AM or EM decomposition environments decomposes with a similar rate, while at the later stages (after 1 year), litter subjected to an AM environment decomposes faster. The difference in the total mass remaining in an AM vs. EM dominant environment increases during the decomposition period from 2–10 years.

Figure 9 Dynamics of plant litter decomposition in AM-dominant vs. EM-dominant environments. (a) Decomposition of total carbon mass from plant litter. (b, c, d) The dynamics of C remaining from labile carbon components (W – water-soluble C fraction, E – ethanol-soluble C fraction, A – acid-hydrolysable C fraction). (e) Dynamics of carbon remaining from recalcitrant C component (N – non-hydrolysable fraction).

Examining the dynamics of carbon loss from distinct individual decomposition pools (Fig. 9b–e) shows that labile carbon components of plant litter (W, A, E) decompose with a similar rate in AM and EM environments. Recalcitrant carbon litter compounds (N) tend to accumulate during the first 2 years. After that, C loss starts to take place in an EM-dominant environment, promoting the accumulation in the N pool compared to an AM-dominant environment. Comparison among distinct litter types reveals that this pattern is not affected by initial litter quality (Figs. D1 and D2).

4 Discussion

Mycorrhizal vegetation types are widely recognized to have a strong impact on plant litter decomposition processes and soil carbon pool dynamics. Yet, the mechanisms of mycorrhizal impacts on the soil C cycle are not well-understood, and available data of the relationship between soil C pools and dominance of distinct mycorrhizal types of vegetation are often contrasting each other at both the local (Craig et al., 2018; Phillips et al., 2013) and global scales (Soudzilovskaia et al., 2019; Steidinger et al., 2019). The matter is additionally complicated by the fact that mycorrhizas affect C cycles via three distinct pathways of (1) provisioning dead mycelium as substrate for decomposition, (2) mediating plant litter quality and amounts, and (3) controlling the environment of plant litter decomposition. Earlier works did not explicitly differentiate between these pathways (Johnson et al., 2006) or focused mainly on the second pathway (Brzostek et al., 2014). Our study is the first attempt to incorporate the impacts of different mycorrhizal environments on litter decomposability, i.e. the third pathway, into a plant litter decomposition model. Herewith, we explicitly focus on impacts of the mycorrhizal environment on the plant litter decomposition process in topsoils, where plant litter is transformed into soil organic matter and carbon compounds are pre-processed for further potential incorporation into particulate organic matter or mineral-associated (i.e. stable) organic matter. We assessed a full range of concepts representing mycorrhizal impacts on labile and stable components of decomposing litter across a wide range of eco-environmental conditions varying in plant species, litter types, and climate variables (Table B1). Overall, the model Myco-Yasso fits the litter decomposition datasets well, especially considering the high level of noise in some of the data and the environmental variation among the sites, in terms of geology, soil quality, and other alike parameters not described by the model. Based on this assessment, we provide insights into the role of distinct mycorrhizal types in long-term decomposition processes of labile and recalcitrant components of plant litter.

4.1 Improved representation of temporal dynamics of litter C

There are still many uncertainties and unknowns in the temporal dynamics of litter decomposition, even though it is an essential process with the soil C cycle. Decomposition encompasses changes in both the composition of soil and litter C and in their breakdown (García-Palacios et al., 2016). This duality, in combination with the long-term nature of the processes involved, makes experimental assessments of temporal dynamics of SOC formation extremely difficult. This is especially true for measuring flows between pools. This supports using modelling approaches, although models may lead to misinterpretations when lacking a theoretical basis. Incorporation of mycorrhizal impacts into Yasso improved the overall model predictions of topsoil C across 10 years, indicating that mycorrhizal impact is a vital factor to be accounted for in analyses of long-term litter decomposition processes, at least in the topsoil layer. The mycorrhizal impacts are likely less visible in the short-term (<3 years), and detectable effects of the mycorrhizal environment on litter decomposition should be assessed over a longer period. This is in agreement with earlier studies (e.g. Paterson et al., 2008) that have shown in short-term ¹³C-labelling experiments that labile and recalcitrant plant litter fractions are utilized by distinct microbial communities, while in the short-term, these communities are not shaped by the presence or activity of mycorrhizal fungi.

4.2 Explicit separation of climate vs. mycorrhizal impacts

Our model allows explicit quantification of mycorrhizal impacts on the decomposition environment and separates these impacts from those of climatic factors. In the original Yasso model, soil C pools are controlled by litter quality and climate, with the “climate” factor implicitly accounting for all global variation in environmental conditions. That original model had high predictive power, especially so for short-term decomposition processes, to which our re-formulation could provide only an incremental improvement. However, the oversimplification of the role of climate without considering microbial factors hinders the ability of models to examine future impacts of alterations in the climate on soil C dynamics (Pongratz et al., 2018). Such a lack of mechanistic and quantitative representation of belowground processes is recognized to be a principal source of uncertainty in our quantification of global terrestrial biogeochemical cycles (Nyawira et al., 2017; Pongratz et al., 2018; Todd-Brown et al., 2013; Trumbore, 2006). There have been recent efforts to incorporate microbial impacts to better represent soil processes in models, such as CORPSE (Sulman et al., 2014), MIMICS (Wieder et al., 2015), and Millennial (Abramoff et al., 2018). Our study is among the first attempts to enable quantification of the impacts of the mycorrhizal environment and to explicitly model mycorrhizal impacts on litter decomposition processes and topsoil C dynamics.

Compared to the original Yasso15 model, the Myco-Yasso model has a lower sensitivity to variation in temperature. The decrease in decomposition sensitivity to temperature suggests that the impact of temperature on decomposition could have been overestimated in previous global modelling attempts that did not consider mycorrhizae as a driving factor. Undoubtedly, the temperature regime controls soil and litter respiration (Hobbie, 1996), making the sensitivity to temperature in a soil C cycle model an essential issue for better estimating future soil C stocks change and its feedback to climate. While modelling approaches allow us to distinguish these mechanisms, separation of these two factors from global field observations is extremely difficult because of a tight correlation of mycorrhizal distributions to gradients in temperature (Barceló et al., 2019; Soudzilovskaia et al., 2015).

4.3 Mycorrhizal impact on labile litter pools is distinct from that on recalcitrant litter pools

We tested four principally distinct concepts on the impact of the mycorrhizal environment on plant litter decomposition. The selected model imposes distinct impacts in terms of both magnitude and direction on labile vs. recalcitrant carbon pools. This finding supports the theory that the turnover of litter depends largely on its composition and recalcitrance of biopolymers (Baldrian, 2017; Berg and McClaugherty, 2008; Cornelissen et al., 2007; Gui et al., 2017), while distinct mycorrhizal types differ in strategies with respect to processing simple organic compounds and recalcitrant compounds (Rajala et al., 2011). This translates into an accumulation of recalcitrant C components of not-yet-decomposed plant litter material in EM-dominated environments, while AM-dominated environments generally promote decomposition. While the presence of AM does not directly affect decomposition, the theory that AMF can exert an indirect influence on this process through regulating free-living groups of decomposers in the soil is well supported. AM fungi alter the physicochemical environment for the microbial community and modify soil bacterial communities (Gui et al., 2017; Nuccio et al., 2013; Offre et al., 2007). AMF stimulate the activity of particular bacteria (Franco-Correa et al., 2010), which are known to be capable of catalysing an efficient degradation of labile and recalcitrant plant litter (Bayer et al., 1998; Kersters et al., 2006). Furthermore, AMF has been shown to prime the decomposition of organic matter by supplying plant-derived labile C to saprotrophic fungi and bacteria (de Vries and Caruso, 2016), which results in higher microbial turnover and respiration and decreasing the soil C pool.

In contrast, efficient nutrient uptake by EM fungi promotes immobilization of soil nitrogen in complex organic molecules of high recalcitrance and therewith promotes the formation of microbial communities, mostly consisting of saprotrophic fungi, able to decompose such recalcitrant organic substrates (Fernandez and Kennedy, 2016; Langley and Hungate, 2003). While multiple studies examining the genetic potential of ectomycorrhizal fungi have shown that EM fungi are capable of producing enzymes degrading complex C and humus (Nicolás et al., 2019), the abundance of such genes is generally low compared to saprotrophic fungal guilds.

Yet, the question of which direction EM impacts soil C the most in the long term has remained unanswered. Similarly, the long-term impacts of AM fungi on saprotrophic communities have, to our knowledge, never been evaluated quantitatively. Our modelling exercises provide a quantitative examination of the long-term consequences of the differential AM and EM impacts on topsoil C across 10 years and suggests that more C is conserved in an EM-dominant environment than an AM environment, particularly due to the accumulation of recalcitrant carbon compounds (independent of the associated litter quality). Moreover, we show that the long-term impacts of both types of mycorrhizas on labile carbon components are similar.

4.4 Future improvements of mycorrhizal impacts of SOC modelling

Our model improves the accuracy of predictions of SOC dynamics even though we assessed the litter decomposition processes in topsoil profiles across 10 years only. Formation of the most recalcitrant soil pool, defined by the Yasso model as “humus” (Tuomi et al., 2011a), was not examined in our study, because we assumed that a 10-year period of litter decomposition for which we had detailed data for model calibration was not long enough for forming humus. Future work should aim at including mycorrhizal impacts on humus formation, linking short- and medium-term decomposition processes to the ultra-long SOC dynamics.

Furthermore, our current work examines the dynamics of topsoil C in terms of labile and stable compounds, yet not addressing the fate of stable, mineral-associated soil C, the ultimate pool of soil-sequestered C. During the last decade, the question of whether mineral-associated soil C originates from labile C components, possibly undergoing microbial transformation (Cotrufo et al., 2015, 2019; Mambelli et al., 2011), or develops through direct sorption of poorly decomposed plant compounds was intensively debated (Bradford et al., 2016; Sokol et al., 2019). Recent research (Sokol et al., 2019) has proposed that both pathways are possible, depending on the capability of the environment to support the release of labile C compounds. While our work does not address the pathway of formation of mineral-stabilized carbon, it provides insights into the important processes preceding C mineral stabilization, as we examine the long-term processes in labile C pools that are potentially available for microbial uptake and the development of recalcitrant plant litter pools that potentially form MAOM by binding to mineral particles. Our study suggests that an EM-dominated decomposition environment tends to promote the accumulation of poorly decomposed plant compounds, supporting the pathway of mineral-associated soil C from undecomposed plant material, which suggests that EM- and AM-dominated ecosystems differ in POM and MAOM fractions contributing to the process of further SOC decomposition. The question of to what extent this pathway dominates the entire flux of soil C into the pool of mineral-associated C needs to be further evaluated. Such evaluation should additionally consider the processes omitted in this study such as fluxes of labile C from the root and fungal exudates and C fluxes originating from the decomposition of dead mycelium of mycorrhizal fungi (Baskaran et al., 2017; See et al., 2021).

5 Conclusions

Our study is the first attempt to model the impacts of the mycorrhizal environment on litter decomposition in topsoil profiles based on differences in carbon release from specific soil chemical pools within pathways of respiration and mass transformation. While mycorrhizae are widely recognized as important factors controlling SOC dynamics, the quantification of these impacts has not been possible thus far. Our work creates a benchmark in such quantifications and enables explicit separation of mycorrhizal impacts from impacts by climate factors in determining topsoil carbon formation processes, which can be applied to a broad range of ecosystems.

The dynamics of decomposition and accumulation of labile and recalcitrant litter compounds is shaped by the abundance of arbuscular and ectomycorrhizal plants in vegetation: if plant litter is decomposed in EM-dominated vegetation, the accumulation of recalcitrant components of that litter in soil is twice as high as in soils of ecosystems dominated by AM vegetation. This difference is likely to affect pathways of accumulation of soil C. We conclude that mycorrhizal traits are an important driver of soil carbon dynamics the impacts of which should be examined quantitatively when estimating future terrestrial carbon storage and predicting impacts of climate change.

Appendix A: Methodological details of the Yasso model structure

The Yasso model represents the decomposing plant litter as five pools of soil carbon, W, A, E, N, and H, varying in recalcitrance. Each pool has its specific decomposition rate (independent from litter type and the initial amount of composition) (Liski et al., 2005; Tuomi et al., 2011a). The model presents the litter decomposition process as a system of linear differential Eq. (A1):

$$\begin{array}{}\text{(A1)}& \mathit{x}{}^{\prime }\left(t\right)=\mathbf{A}\left(D\right)\mathit{x}\left(t\right)+b\left(t\right),\mathit{x}\left(\mathrm{0}\right)={\mathit{x}}_{\mathrm{0}},\end{array}$$

where x(t) is a vector describing the mass of individual carbon pools as a function of time (t); x(0) = x₀ represents the initial amount of each carbon fraction; b(t) is the litter input; and A(D) is a matrix describing the total decomposition as a function of climatic conditions (D), where the diagonal values represent the fraction of C being removed from the pool and the non-diagonal terms specify the amount of C transferred to other pools (Viskari et al., 2020).

The total amount of carbon released from individual W, A, E, N, and H pools is the result of two fluxes: (1) carbon transformation from and to other pools and (2) carbon that is not transferred to other pools but instead released to the atmosphere as carbon dioxide. The mass fluxes between the different pools and to outside the system are accordingly determined by two parameter sets: p_ij represents the mass transportation between pools, and α_i represents the total decomposition rate of each pool, i.e. the C mass leaving the pool (the sum of C transfer to other pools and C released into the atmosphere). The total mass flux between two pools is thus a product of these two parameters; e.g. the mass flux from pool A to pool W is α_A⋅p_AW.

The total decomposition represented by matrix A(D) within the whole system can be represented as a mathematical equation with a mass flow matrix, where parameters $p_{ij}\in [\mathrm{0},\mathrm{1}]$ denote the flows between each pair of W, A, E, N, and H pools i and j, and K(D) represents the impact of climate on decomposition rate Eq. (A2).

$$\begin{array}{}\text{(A2)}& \mathbf{A}\left(D\right)=\left[\begin{array}{rrrrr}-\mathrm{1}& p_{WA}& p_{EA}& p_{NA}& \mathrm{0}\\ p_{AW}& -\mathrm{1}& p_{EW}& p_{NW}& \mathrm{0}\\ p_{AE}& p_{WE}& -\mathrm{1}& p_{NE}& \mathrm{0}\\ p_{AN}& p_{WN}& p_{EN}& -\mathrm{1}& \mathrm{0}\\ p_H& p_H& p_H& p_H& -\mathrm{1}\end{array}\right]\cdot \mathbf{K}\left(D\right)\end{array}$$

In the matrix K(D), each element k_i(D) describes the decomposition of W, A, E, N, and H as a function of temperature (T) and annual precipitation (P) modelled according to Eq. (A3):

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{k}_i\left(D\right)& ={\displaystyle \frac{{\mathit{\alpha }}_i}{J}}{\sum }_{j=\mathrm{1}}^J\exp \left(b_{i\mathrm{1}}{T}_j+b_{i\mathrm{2}}{T}_j^{\mathrm{2}}\right)\\ & \left(\mathrm{1}-\exp \left(g_{i}P\right)\right),i\in \mathit{\left\{}WAENH\mathit{\right\}},\end{array}\end{array}$$

where α_i represents decomposition rate parameters; b_i1 and b_i2 are parameters describing the dependency of heterotrophic respiration on temperature, assessed through a Gaussian model (Tuomi et al., 2009); and g_i is a parameter describing the dependency of heterotrophic respiration of precipitation, assessed through an exponential function (Tuomi et al., 2009). Systematic error in the litter decomposition resulting from litter leaching out of the litter bags was corrected by leaching parameters.

Appendix B: Methodological details of calibration and databases of litter decomposition data used

There are three main litter decomposition databases used in both the original Yasso modelling (Tuomi et al., 2011a) and our new model parameterization: CIDET dataset with the measurements from Canada (Trofymow, 1998), the LIDET dataset with data from the USA and Central America (Gholz et al., 2000), and the Eurodeco (ED) dataset with data gathered from several European research projects (Berg et al., 1991). The distributions of these experimental sites are shown in Fig. B1. Details of these datasets used to parametrize our new model are shown in Table B1.

The original Yasso model also uses a dataset with information of SOC accumulation along thousands of years at sites in Finland (Liski et al., 2005) and a large global soil C stock measurements dataset (Zinke et al., 1986) to infer the dynamics of the most stable carbon – humus pool in soil (see Fig. 1). However, given our focus on impacts of mycorrhizas on the dynamics of chemical compounds during plant litter decomposition, and given that the LIDET, CIDET, and ED databases of litter decomposition, used for calibration of our modified YASSO model, store data of 0–10.5 years of decomposition, we assumed no measurable amounts of humus being formed during this time frame. Therefore, the mycorrhizal impacts on H-related decomposition terms (α_H and p_NH; see Fig. 1) were set to zero.

Figure B1 The distribution map of litter bag experiment sites.

Table B1 Dataset description with general environmental conditions.

Appendix C: Details of Myco-Yasso parameters and performance

Table C1 Posterior and 95 % Bayesian credibility intervals (confidence limits) for the Yasso_Myco parameters.

Figure C1 Scatter plot of predictions for C loss from plant litter made by the original Yasso15 model (grey circles) and predictions made by the Myco-Yasso model (blue dots) compared to experimental measurements.

Figure C2 Correlations between parameters of the Myco-Yasso C model. The gradient from the most intensive blue colours to the most intensive red colours indicate correlations from completely negative (−1) to completely positive (1).

Figure C3 Model sensitivity of the Yasso15 and Myco-Yasso models to individual groups of parameters. The impact of increases by 1 % of each parameter is shown.

Figure C4 Sensitivity of Yasso15 and Myco-Yasso models to 1 % increase in mass flow parameters. The two leftmost bars show the sensitivity to the joint impact of all p-term parameters being increased by 1 %. The other bars show the impact of 1 % increase in individual p terms: p_WA – C flux from the W to A pool, p_EN – C flux from the E to N pool, and p_WN – C flux from the W to N pool.

Figure C5 Improvement of performance comparing the Myco-Yasso model to the original Yasso model over the decomposition period. Bars represent the relative RMSE differences between Yasso15 and Myco-Yasso per period. The line with dots shows the absolute value of the RMSE differences (Yasso15 – Myco).

Appendix D: Experiment on litter decomposition dynamics of different initial litter quality in AM-dominated vs. EM-dominated environments

Figure D1 Dynamics of plant root litter decomposition in AM-dominated vs. EM-dominated environments. (a) Loss of total carbon mass from root litter. (b, c, d) Dynamics of loss of labile carbon components (W – water-soluble C pool, A – acid-hydrolysable C pool, E – ethanol-soluble C pool). (e) Dynamics of loss of recalcitrant (non-hydrolysable) carbon (N pool). The initial W, A, E, and N composition of root material is 17 % – W, 55 % – A, 9 % – E, and 20 % – N (typical for plant roots).

Figure D2 Dynamics of plant foliage (leaf) litter decomposition in AM-dominated vs. EM-dominated environments. (a) Loss of total carbon mass from foliage litter. (b, c, d) Dynamics of loss of labile carbon components (W – water-soluble C pool, A – acid-hydrolysable C pool, E – ethanol-soluble C pool). (e) Dynamics of loss of recalcitrant (non-hydrolysable) carbon (N pool). The initial W, A, E, and N composition of leaf material is 25 % – W, 45 % – A, 12 % – E, and 18 % – N (typical for plant foliage).