Model of Ectomycorrhiza Contribution in Forest Soil C and N Dynamics in a Frame of System of Models EFIMOD3

1. Introduction

The development of mathematical modelling in ecology, soil science, forestry and other biological disciplines began in the 20^th century with very simple functions, followed by a steady increase in the complexity of the model structure as new experimental data and ideas for quantification and parameterisation of processes in natural systems emerged.

The first equations for the dynamics of organic debris and soil organic matter (SOM) described the formation/accumulation of SOM, but then the focus shifted to the processes of SOM decomposition. At that time, the development of SOM models focused on SOM as an important factor for plant nutrition and improvement of edaphic conditions in terrestrial ecosystems [1,2,3]. Nowadays, the emphasis has shifted back to the sequestration of “soil organic carbon” (SOC instead of SOM) in the soil system due to the problem of climate change and the important role of soil in the global carbon cycle [4,5,6,7].

The number of recent models of SOM dynamics is quite large (see review in [8]). The effect of biological factors in the recent models is mostly limited by microbiological activity. SOM dynamics have not been directly linked to the whole soil biota. Modern models satisfactorily reproduce the processes of SOM transformation, decomposition and formation of SOM pools with relative stability or protection against its mineralisation within the framework of the classical concept of “plant – soil” interaction and ecosystem feedback through the “plant debris – soil nutrients” exchange. Some recent SOM models are integrated with plant/vegetation simulators and reproduce the organic carbon cycle at the ecosystem level [9,10,11]. This has led to increased interest in studying the finer biotic mechanisms of “plant – soil biota – soil” interactions.

Two widespread mechanisms of “plant – soil” biotic interactions are (i) the activation of rhizosphere soil microorganisms by plant roots, the so-called the “priming effect” [12,13], and (ii) the symbiosis of plant roots with mycorrhizal fungi [14,15,16]. Symbiotic microbial nitrogen fixation and nutrient transfer from the canopy to the soil through litterfall are also such “plant – soil” interactions [17,18,19]. The cumulative effect of these mechanisms of “plant – soil biota – soil” interactions leads to the optimisation of the edaphic environment in terrestrial ecosystems and increased plants growth and vitality.

The rhizosphere priming effect is well documented and has been studied in detail [12], which has allowed the mathematical simulation of this phenomenon [20,21]. However, the mycorrhizal symbiosis is a more powerful type of “plant – soil biota – SOM” interactions that has been underexplored for about 150 years. It is a widespread phenomenon in the world’s vegetation. There are five types of such symbiosis, but two are dominant. These are ectomycorrhiza in temperate and boreal forests and arbuscular mycorrhiza in forests and grasslands, mostly in warmer climates.

Ectomycorrhiza (EMR) forms a mycelial mantle (sheath) on the growing root tip that produces a large amount of extramatrical mycelium (EMM) that penetrates a whole volume of organic and organomineral soil horizons. Arbuscular mycorrhizal cells live inside the roots from where the EMM hyphae occupy a large soil volume. EMR explores the entire soil volume and its dense network of EMM “mines” the nutrients from SOM, thus covering their lack for fungal growth [22,23]. EMM accounts for one third of the total biomass of soil microorganisms [24].

Mycorrhizal fungi are obligatory biotrophs that exist on matter and energy obtained from the host plant in the form of root exudates (RE). The products of plant photosynthesis are the driving force behind the functioning of EMR. Trees spend 10–25% of their NPP for RE [23]. RE include carbohydrates, amino acids, organic acids and some other organic compounds [14,15,25]. However, RE do not cover all EMM growth requirements due to an imbalance of C and N in their composition [26]. First of all, RE in boreal and temperate forests are deficient in nitrogen, whereas the limiting nutrient for more southern vegetation zones is phosphorus. Therefore, EMR performs a rather specific SOM decomposition pattern, replenishing the missing N by “mining” from SOM [14,23,27,28]. EMM depolymerises large organic molecules of N-rich recalcitrant SOM for N extraction (“mining”), but does not mineralise all the “dissolved organic carbon” (DOC) obtained. This is a major difference between mycorrhiza and soil saprotrophic fungi: EMM produces a lot of DOC [24,29], but contributes about 2% to the total SOM mineralisation [22,30,31]. This DOC after N mining plus the rest of acid RE after mycorrhizal consumption of carbohydrates and amino acids are biochemical agents of soil minerals’ dissolution [32,33]. This resulted in an increase in the availability of phosphorus and other elements to plant roots.

EMM has a short lifespan, varying from a few days to several weeks [14,34,35]. This results in a high EMM turnover rate. Therefore, EMM produces a large pool of litter, which is larger than leaf litter and can reach 62% of the total annual above- and belowground litter from trees [36,37,38,39].

The decomposition of this very large and previously unknown pool of fungal litter is carried out by all soil biota: bacteria, fungi, micro- and mesofauna of food webs. However, the role of soil fauna in the fate of this specific litter is mostly missing in EMR studies [40], leading to an underestimation of the contribution of EMM litter to SOM formation. In the case of only microbial decomposition, the intermediate product of microbial exoenzymatic depolymerisation of SOM is DOC. It is consumed by microorganisms for their own biomass growth [41]. So, in effect, new microbial biomass, CO₂ and H₂O remain after complete SOM mineralisation by microorganisms. Nevertheless, there is a concept in soil science that the DOC produced microbial activity is an important source of new SOM formation by stabilising on a surface of soil minerals [42,43,44].

A different situation occurs in the case of mycelial litter decomposition with the participation of fauna of soil food webs as active agent of litter transformation [45]. Soil micro- and mesofauna consume litter, transform organic matter and produce solid-phase metabolic products in the form of casts (coprolites). They have been known in soil ecology since the 19^th century as main products of SOM formation [46,47]. These casts, which are rich in N and not yet bound to the soil mineral matrix, sesquioxides and calcium, represent a previously unrecognised labile SOM fraction, which determines a positive role of EMR for plant growth and also a significant contribution of EMR to SOM recovery. Furthermore, the formation of available ammonia N in food webs is an additional factor in the positive effect of EMR on tree growth. It occurs because some excess N is produced during the consumption of microorganism (including EMM) by soil fauna, as a significant part of the consumed C is used for faunal respiration [48,49]. This pathway of EMM litter transformation and decomposition is consistent with the concept and experimental data on the prominent role of food webs in resource acquisition and consumption in soils of terrestrial ecosystems [21,48,49,50,51,52,53]. With regard to soil food webs, their active role has been confirmed by the oscillatory respiration on an hourly scale, which reflects the high activity of “predator – prey” type of relationships between the soil fauna of “microbial grazers” and microorganisms [54,55].

This role of EMR in SOM recovery and sequestration represents a feedback in the mycorrhizal symbiosis that determines a well-known positive effect of EMR on tree growth [56,57]. Thus, the increase in biomass increment, photosynthetic rate or respiration in tree seedlings with ectomycorrhiza was 1.2–2.6 times higher than the same in non-infected plants [14,58,59,60]. Currently, there is a growing interest in mathematical modelling of the functioning of EMR as a powerful phenomenon reflecting the interaction between tree and soil. Thus, there are studies on principles and approaches of EMR simulation [61,62,63], modelling of special patterns of EMM distribution in forest ecosystems [64], model assessment of nitrogen uptake by EMR fungi [4], quantification of factors of EMM decomposition [65], and well-parameterised simulation of EMM growth [66]. More promising is a successful attempt to integrate the EMR model into the whole ecosystem simulator, MYCOFON [67].

The majority of recent EMR models have two specific features: (i) they consider the dynamics of EMM litter decomposition as controlled only by microorganisms, without including the soil fauna of food webs; (ii) the models are not integrated into the whole ecosystem models. Therefore, the aim of this work was to develop an EMR model that can reproduce the role of all soil biota with special reference to soil food webs, SOM dynamics and plant nutrition as recently conceptualised [40]. This model has been developed for integration with the Romul_Hum model of SOM dynamics [8,53] and the most recent version of the forest ecosystem model, EFIMOD3 [68].

2. Model description

2.1. Basic postulates

The following postulates, based on previous studies on the functioning of EMM, were used to construct the model. In the mycorrhizal symbiosis, the fungi are very active and consume all N of the RE for their own growth (mainly for the EMM), thus making impossible the transfer of N of exudates to the soil through the dense Hartig net of the mycelial mantle.

EMM decomposes SOM to cover the N deficiency for EMM growth only. It is selective N mining from recalcitrant stable SOM with increased N concentration. In effect, the mycorrhizal fungi are de facto “working for themselves”.

EMM has a short lifespan. Therefore, EMM produces a large litter pool, which is supposed to be a dominant component in the structure of total litter input in forest ecosystems.

This large litter pool of EMM is consumed by soil micro- and mesofauna of different trophic levels in soil food webs. It ensures SOM replenishment by solid-phase products of soil faunal metabolism, namely labile SOM from faunal excreta.

Additional N in EMR occurs through N mining during selective decomposition of recalcitrant SOM and also after EMM consumption by soil fauna in food webs. This N becomes available to roots after microbial decomposition of rapidly decomposing N rich labile SOM originated from faunal excreta. This explains the positive role of EMR in tree nutrition and growth.

2.2. Initialisation

The model has been developed for steady-state conditions with developed mycorrhiza in a forest ecosystem, where mycorrhizal growth is in equilibrium with its die off and consumption by soil fauna (faunal consumption is a dominant process) and is linked to growth of fine roots. All calculations in the model are performed in terms of carbon stocks per unit area (e.g. in kg m−2), with a time step of 1 day. Hydrothermal conditions are assumed to be optimal with soil moisture at 60% of the specific water retention and a soil temperature of 20° C. Total SOC represents the sum of the carbon stocks of the organic and organomineral horizons (O and Ah/AhE) in the topsoil. The pools of SOC, SON, biomass of fine roots, a rate and composition of RE are obtained from the EFIMOD3 system of models as part of which the EMR model was developed. All parameters and variables are presented in Table 1 and Table 2.

The biomass of the ectomycorrhizal mantle reaches 30–60% of fine roots biomass [14,72,73]. We postulated it as a constant fraction of 0.4 of the fine root biomass of trees. To date, the fraction of fine roots determined experimentally in the field studies includes both roots and ectomycorrhizal mantle. It follows, that a fraction of the ectomycorrhizal mantle biomass can be isolated from the experimentally measured fine root biomass. Therefore, we used the fine roots biomass to determine the initial value and dynamics of the fungal mantle in the EMR model, which has been integrated with the EFIMOD3 model. It is included in the model formulation as part of the ectomycorrhizal biomass. The lifespan of the mantle is 2–3 weeks, while the lifespan of fine roots without mantle is about 3 years [14].

Ectomycorrhiza develops a dense network of EMМ, which also has the same lifespan as the mantle [14,34,35,74]. Its biomass represents 30% of the total biomass of all soil microorganisms [24], and this parameter was used to calculate the initial biomass of EMM. The biomass of soil microorganisms can be estimated in the absence of measured data using the equation derived from [75]:

$$C_{mic}=3.91-0.30\times C_{SOM},$$

(1)

where $C_{mic}$ is the C of all microorganisms (all fungi and all bacteria, but excluding microfauna), % of total SOC ($C_{SOM}$).

The biomass of EMM ($C_{EMM}$) can be calculated as

$$C_{EMM}=0.3\times C_{mic}.$$

(2)

The C:N ratio of the EMM ($C{N}_{EMM}$) is also required to calculate the EMM dynamics in the model. It can be calculated using the following equation [53]:

$$C{N}_{EMM}=0.85\times C{N}_{SOM},$$

(3)

where $C{N}_{SOM}$ is the C:N ratio of soil organic matter.

2.3. Root exudate input

The input of root exudates ($C_{RE}$, kg [C] m⁻² day⁻¹) is determined in the photosynthesis submodel of EFIMOD3 as part of the gross primary production (GPP) of the trees, or this input can be obtained from the available measured data. The actual RE consumption (${∆C}_{RE}^{cons}$) can be partially reduced due to competition with priming ($k_{pr}^{EMM}$) and RE consumption efficiency ($k_{RE}^{cons}$):

$${∆C}_{RE}^{cons}=C_{RE}\times k_{pr}^{EMM}\times k_{RE}^{cons}.$$

(4)

The second and third terms of this equation are assumed to be 1.0 in this pilot version of the model due to the lack of experimental data.

2.4. EMM growth

Consumed RE are used for EMM growth (${∆C}_{EMM}^{RE}$) taking into account (i) the N pool in RE ($N_{RE}$) and (ii) the need for RE carbon for EMM growth respiration. It can be expressed by two equations:

$${∆C}_{EMM}^{RE}=N_{RE}\times {CN}_{EMM},$$

(5)

where ${CN}_{EMM}$ is the C:N ratio of the EMM biomass;

$$R_{EMM}={∆C}_{EMM}^{RE}\times {(k}_{EMM}^R/k_{EMM}^{Gr}),$$

(6)

where $R_{EMM}$ is RE respiration as a cost of EMM growth; $k_{EMM}^{Gr}$ is a portion of consumed C spent on EMM growth; $k_{EMM}^R$ is a portion of consumed C spent on respiration (with the condition that $k_{EMM}^{Gr}+k_{EMM}^R=1.0$).

The use of RE for EMM growth and respiration does not utilise the entire RE pool due to insufficient nitrogen content (the C:N ratio of RE is significantly higher than the C:N ratio of EMM biomass). The rest of the RE (${rest}_{RE}$) is calculated as follows:

$${rest}_{RE}={∆C}_{RE}^{cons}-R_{EMM}-{∆C}_{EMM}^{RE}.$$

(7)

This pool contains a significant amount of organic acids as RE proteins and carbohydrates are consumed by EMM for biomass production and respiration.

2.5. N mining

Complete assimilation of RE carbon for EMM growth requires additional nitrogen, which fungi “mine” from SOM by decomposing it to soluble compounds (dissolved organic carbon, DOC) during the process of enzymatic depolymerisation of large organic molecules [15,28]. It is therefore necessary to determine N deficiency for the assimilation of the rest of RE for EMM growth (${∆N}_{mng}$):

$${∆N}_{mng}={\displaystyle \frac{{rest}_{RE}}{{CN}_{EMM}}}.$$

(8)

It is possible to calculate the amount of SOC, which is “mined” from the total pool of SOC ($∆DOC$), using the following equation:

$$∆DOC={∆N}_{mng}\times {CN}_{SOM},$$

(9)

where ${CN}_{SOM}$ is the C:N ratio of total SOM. This SOC represents a depolymerised soluble organic matter with a strongly reduced N content after its “mining” for the assimilation of the remaining C of RE for EMM growth:

$${∆C}_{EMM}^{mng}={∆N}_{mng}\times {CN}_{EMM},$$

(10)

where ${∆C}_{EMM}^{mng}$ is the additional growth of EMM using “mined” N and the rest of RE.

The processes of RE assimilation for EMM growth and N mining from total SOM lead to the production of two fractions of nitrogen-poor and acidified soluble organic matter. This mixture represents “dissolved organic carbon” (DOC) at a given time step i, which can be described by the following equation:

$$DO{C}^i={DOC}^{i-1}+{rest}_{RE}-{∆C}_{EMM}^{mng}+∆DOC.$$

(11)

2.6. EMM consumption and mortality

The biomass of EMM decreases under the influence of natural mortality and consumption by the soil fauna of food webs:

$${∆C}_{EMM}^{mort}=C_{EMM}\times k_{EMM}^{mort},$$

(12)

$${∆C}_{EMM}^{cons}=C_{EMM}\times k_{EMM}^{cons},$$

(13)

where ${∆C}_{EMM}^{mort}$ is the EMM mortality; ${∆C}_{EMM}^{cons}$ is the consumption of EMM by soil micro- and mesofauna in food webs; $k_{EMM}^{mort}$ and $k_{EMM}^{cons}$ are the rates of natural mortality and consumption by soil fauna, respectively.

The total carbon budget of EMM at a given time step i ($C_{EMM}^i$) can be expressed as the difference between the corresponding value at a previous time step ($C_{EMM}^{i-1}$) and all incoming and outcoming fluxes, taken with the corresponding sign:

$$C_{EMM}^i=C_{EMM}^{i-1}+{∆C}_{RE}^{cons}-R_{EMM}-{∆C}_{EMM}^{mort}-{∆C}_{EMM}^{cons}.$$

(14)

According to the model formulation, the EMM necromass becomes a resource for microbiological decomposition, as postulated in microbiological models of SOC dynamics [5,20].

2.7. Food webs dynamics

The main consumers of EMM biomass are the soil micro- and mesofauna of the food webs. A specific feature of soil fauna physiology is the formation of solid excreta as a major metabolic by-product. This process is important for understanding the role of micro- and mesofauna in the soil carbon balance. Cast excretion (${∆C}_{MEF}^{excr}$), respiration ($R_{MEF}$) and growth (${∆C}_{MEF}$) of micro- and mesofauna can be described by operating on the value of EMM biomass consumed, cast production rate ($k_{MEF}^{excr}$) and respiration rate ($k_{MEF}^R$):

$${∆C}_{MEF}^{excr}={∆C}_{EMM}^{cons}\times k_{MEF}^{excr},$$

(15)

$$R_{MEF}={∆C}_{EMM}^{cons}\times k_{MEF}^R,$$

(16)

$${∆C}_{MEF}={∆C}_{EMM}^{cons}-{∆C}_{MEF}^{excr}-R_{MEF}.$$

(17)

These three equations describe the functioning of the soil fauna in relation to the use of EMM as a basic level of food webs, although in other cases the basic level is represented by saprotrophic microorganisms.

In this case, part of the EMM consumed has been used for respiration, while some N remains excess and is excreted as ammonium (${∆N}_{MEF}^{min}$) [48,49]. This can be calculated using the following equation ([48,21], adopted):

$${∆N}_{MEF}^{min}={\displaystyle \frac{{∆C}_{EMM}^{cons}}{C{N}_{MEF}}}-{\displaystyle \frac{{∆C}_{EMM}^{cons}}{C{N}_{EMM}}},$$

(18)

where $C{N}_{MEF}$ is the C:N ratio of micro- and mesofauna. The proportion of the faunal population ($C_{MEF}$) that die-off is described as follows:

$${∆C}_{MEF}^{mort}=C_{MEF}\times k_{MEF}^{mort},$$

(19)

where ${∆C}_{MEF}^{mort}$ is the faunal necromass production (in terms of carbon) and $k_{MEF}^{mort}$ is the faunal mortality rate.

Micro- and mesofauna are a food resource for higher trophic levels with larger predators in soil food webs. In this case, the consumption of micro- and mesofauna can be represented as follows:

$${∆C}_{MEF}^{cons}=C_{MEF}\times k_{MEF}^{cons},$$

(20)

where ${∆C}_{MEF}^{cons}$ is the consumed biomass (in terms of carbon) of micro- and mesofauna and $k_{MEF}^{cons}$ is the consumption rate.

An entire carbon budget of the trophic level of micro- and mesofauna in food webs at a given time step i can be represented as follows:

$$C_{MEF}^i=C_{MEF}^{i-1}+{∆C}_{EMM}^{cons}-R_{MEF}-{∆C}_{MEF}^{excr}-{∆C}_{MEF}^{mort}-{∆C}_{MEF}^{cons}.$$

(21)

The representation of the higher trophic levels of food webs (soil macrofauna) in the current version of the model does not include aboveground fauna and earthworms. The macrofauna consume the biomass of the micro- and mesofauna. The fluxes associated with the soil macrofauna are described in the same terms as for the micro- and mesofauna:

$${∆C}_{MAF}^{excr}={∆C}_{MEF}^{cons}\times k_{MAF}^{excr},$$

(22)

$$R_{MAF}={∆C}_{MEF}^{cons}\times k_{MAF}^R,$$

(23)

$${∆C}_{MAF}={∆C}_{MEF}^{cons}-{∆C}_{MAF}^{excr}-R_{MAF},$$

(24)

$${∆N}_{MAF}^{min}={\displaystyle \frac{{∆C}_{MEF}^{cons}}{C{N}_{MAF}}}-{\displaystyle \frac{{∆C}_{MEF}^{cons}}{C{N}_{MEF}}},$$

(25)

$${∆C}_{MAF}^{mort}=C_{MAF}\times k_{MAF}^{mort},$$

(26)

$${∆C}_{MAF}^{cons}=C_{MAF}\times k_{MAF}^{cons},$$

(27)

$$C_{MAF}^i=C_{MAF}^{i-1}+{∆C}_{MEF}^{cons}-R_{MAF}-{∆C}_{MAF}^{excr}-{∆C}_{MAF}^{mort}-{∆C}_{MAF}^{cons}.$$

(28)

The outputs of the model are shown in Table 1. This list can be modified according to the specifics of the simulation tasks. The general scheme of the model is shown in Figure 1.

Figure 1. Flowchart of the model of EMR functioning in forest soils. The solid lines represent C fluxes and the dotted lines represent N fluxes. The numbers in the circles represent the processes of EMR functioning reproduced in the model.

Figure 1. Flowchart of the model of EMR functioning in forest soils. The solid lines represent C fluxes and the dotted lines represent N fluxes. The numbers in the circles represent the processes of EMR functioning reproduced in the model.

3. Sensitivity analysis

The analysis of model sensitivity to uncertainty in parameter estimation was carried out in the same way as for the priming model [21]. The approach of testing the sensitivity to each parameter in isolation is often criticised for its inability to account for the complexity arising from the possible non-linear interactions. Therefore, we used the approach proposed in [76]. According to this approach, we generated the uniform distributions for each input parameter in a range of ±30% from the default values. A dataset containing all possible combinations of parameter values was then compiled. The model was run for one year (365 daily time steps) for each record in the dataset, and the sum of EMM, soil micro-, meso- and macrofauna respiration (${R_{EMM}+R}_{MEF}+R_{MAF}$) was calculated as the target variable, whose calculation involves (either directly or indirectly) all analysed parameters. Then, the multiple linear regression analysis of the standardised values ($x_i^S$)

$$x_i^S=\left({\displaystyle \frac{x_i-\overline{x}}{S\left(x\right)}}\right)$$

(29)

of the output variable with respect to the set of standardised values of the input parameters (I₁ … I_n) was performed, where $x_i$ is the input parameter or target variable, $\overline{x}$ is the mean value, and $S\left(x\right)$ is the standard deviation. The standardised regression coefficients c₀ … c_n (the subscript 0 refers to an intercept) can be considered as a measure of the sensitivity of the model to the variation of a given parameter, while the value of the coefficient of determination, R², indicates the non-linearity of the model (the lower the R², the higher the non-linearity).

The results of the sensitivity analysis (Table 3) showed that the target variable is most sensitive to C:N ratios, especially those of EMR fungi. Other influential parameters include the respiration rates of the soil meso- and macrofauna. The analysis also showed that the model is almost linear.

Table 3. Analysis of model sensitivity to uncertainty in input parameters.

Parameter	Parameter Name	Standardised Coefficient of Linear Regression
C:N ratio of root exudates	${CN}_{RE}$	−0.245***
C:N ratio of ectomycorrhiza	${CN}_{EMM}$	0.723***
C:N ratio of soil mesofauna	${CN}_{MEF}$	−0.344***
C:N ratio of soil macrofauna	${CN}_{MAF}$	−0.164***
respiration rate of ectomycorrhiza	$k_{EMM}^R$	0.238***
mortality rate of ectomycorrhiza	$k_{EMM}^{mort}$	−0.0553***
consumption rate of ectomycorrhiza	$k_{EMM}^{cons}$	0.0730***
mortality rate of soil micro- and mesofauna	$k_{MEF}^{mort}$	−0.0126***
consumption rate of soil micro- and mesofauna	$k_{MEF}^{cons}$	0.0167***
respiration rate of soil micro- and mesofauna	$k_{MEF}^R$	0.326***
respiration rate of soil macrofauna	$k_{MAF}^R$	0.142***
Intercept		~0
R2		0.921

Note: ^*** indicates the significance level (p < 0.001).

Table 3. Analysis of model sensitivity to uncertainty in input parameters.

Parameter	Parameter Name	Standardised Coefficient of Linear Regression
C:N ratio of root exudates	${CN}_{RE}$	−0.245***
C:N ratio of ectomycorrhiza	${CN}_{EMM}$	0.723***
C:N ratio of soil mesofauna	${CN}_{MEF}$	−0.344***
C:N ratio of soil macrofauna	${CN}_{MAF}$	−0.164***
respiration rate of ectomycorrhiza	$k_{EMM}^R$	0.238***
mortality rate of ectomycorrhiza	$k_{EMM}^{mort}$	−0.0553***
consumption rate of ectomycorrhiza	$k_{EMM}^{cons}$	0.0730***
mortality rate of soil micro- and mesofauna	$k_{MEF}^{mort}$	−0.0126***
consumption rate of soil micro- and mesofauna	$k_{MEF}^{cons}$	0.0167***
respiration rate of soil micro- and mesofauna	$k_{MEF}^R$	0.326***
respiration rate of soil macrofauna	$k_{MAF}^R$	0.142***
Intercept		~0
R2		0.921

Note: ^*** indicates the significance level (p < 0.001).

4. Model validation

Despite a very large number of publications on mycorrhizal symbiosis (cited in [14,15]), a number of papers on EMR dynamics with a complete dataset of environmental and physiological parameters of EMR functioning in forest ecosystems are not so abundant. We found 4 papers with sufficiently detailed data on mycorrhizal respiration, soil climate, EMM biomass production and soil properties. All papers were carried out on the same experimental plot at the same time. This allowed us to reconstruct all the input parameters and data necessary for the validation of the EMR model.

The comprehensive EMM respiration data from [77] were used to validate the model. In this paper, the average EMM respiration values for the spring, summer and autumn months are given for Larix gmelinii var. principis-rupprechtii plantations of 3 different ages (11, 25 and 45 years) located in northeastern China. Data from [78] were used to reconstruct the seasonal dynamics of mycorrhizal respiration. In this paper, seasonal dynamics of heterotrophic and autotrophic respiration (including mycorrhizal respiration), soil temperature and soil moisture are provided with a time step of 2 weeks. Data from [78] were interpolated to daily resolution using a cubic spline method [79]. The seasonal dynamics of mycorrhizal respiration was calculated in proportion to the seasonal dynamics of autotrophic respiration:

$$R_{EMM}^{DOY}={\displaystyle \frac{R_A^{DOY}}{\sum R_A^{DOY}}}\times \sum DO{Y}^0\times \overline{R_{EMM}^{DOY}},$$

(30)

where $DOY$ is the Julian date ($\sum DO{Y}^0$ denotes the total number of days for which the measured data are available), $R_{EMM}^{DOY}$ is the measured ectomycorrhizal respiration at DOY, $R_A^{DOY}$ is the measured autotrophic respiration at $DOY$, $\overline{R_{EMM}^{DOY}}$ is the mean daily ectomycorrhizal respiration.

Net primary production (NPP) data of different tree organs for these plantations from [80] were used to create an input dataset on the amount of root exudates entering the soil. The global root exudate flux is estimated to be about 9% (7–14%) of the global annual gross primary productivity (GPP) for mature trees [81]. At the same time, a number of studies estimate NPP to be 0.47±0.04 of GPP [82]. As stated in [83], the production of root exudates is proportional to the proportion of total NPP used for root growth. Therefore, the following equation was used to calculate root exudate (RE) production for 3 scenarios of exudate production rate and 3 stand ages:

$$RE={\displaystyle \frac{NP{P}_{tot}^{age}}{0.47}}\times R{E}_{prt}\times {\displaystyle \frac{NP{P}_{roots}^{age}\times NP{P}_{tot}^{mature}}{NP{P}_{roots}^{mature}\times NP{P}_{tot}^{age}}},$$

(31)

where $NP{P}_{tot}^{age}$ is the total NPP for stands aged 11, 25 or 45 years; $R{E}_{prt}$ is the proportion of root exudates in total GPP for 3 different scenarios of exudate production (0.07, 0.09 and 0.14, respectively); $NP{P}_{roots}^{age}$ is the NPP of roots for stands aged 11, 25 and 45 years; $NP{P}_{tot}^{mature}$ is the total NPP for stand aged 45 years; $NP{P}_{roots}^{mature}$ is the NPP of roots for stand aged 45 years.

To reconstruct the seasonal dynamics of root exudate supply, it was assumed that the daily NPP is proportional to the amount of mineralised soil N. The dependence of the mineralisation rate on soil temperature and moisture from the Romul_Hum model [8] was used to calculate the seasonal dynamics of N mineralisation and, consequently, root exudate supply. Finally, ectomycorrhizal respiration calculated by the model was compared with measured data from [77] using Theil’s U index of inequality (0…1; [84]):

$$U={\displaystyle \frac{\sqrt{\frac{1}{n}\times {\sum }_1^n{\left(o_t-m_t\right)}^2}}{\sqrt{\frac{1}{n}\times {\sum }_1^n{o}_t^2}+\sqrt{\frac{1}{n}\times {\sum }_1^n{m}_t^2}}},$$

(32)

where n is the number of observations, o_t is the observed EMM respiration at day t, m_t is the simulated EMM respiration at day t. The case where U = 0 corresponds to full identity of measured and simulated values.

Figure 2. Comparison of measured and simulated seasonal dynamics of ectomycorrhizal respiration for stands of different ages. U is the Theil’s U index of inequality.

Figure 2. Comparison of measured and simulated seasonal dynamics of ectomycorrhizal respiration for stands of different ages. U is the Theil’s U index of inequality.

The validation results (Fig. 2) showed that the model underestimates EMM respiration in very young plantation with low stand density. However, the model overestimates respiration in dense 25-year-old forest. The best agreement between experimental and simulated data was obtained for 45-year-old larch forest. The low Theil’s U index in all cases indicates the satisfactory results of this validation.

5. Testing the EMR model

The EMR model was tested using data from the 1 ha (100×100 m) permanent sample plot in a mixed uneven-aged forest on Albic Luvisol dominated by Pinus sylvestris L., Betula spp. and Picea abies L. The forest community is located in the Prioksko-Terrasny Biosphere Reserve, south of the Moscow region, Russia (54.88876° N, 37.56273° E). Information on the vegetation and soils of this site are presented in publications [68,85].

The vegetation composition and soil conditions vary considerably within the plot, which was taken into account in the simulations. According to the formulation of the EFIMOD3 model, the computational experiments took into account the spatial heterogeneity of the distribution of tree roots, depending on species-specific features, location, size and age of the trees. The spatial distribution of root biomass in turn determined spatial differences in the autotrophic respiration rates calculated in EFIMOD3 and the amount of root exudates input into the EMR model. In addition, the spatial structure, species composition and size class distribution of the simulated stand determined spatial differences in the input of species-specific fractions of above- (needles and leaves, branches, stems) and belowground (fine and coarse roots) litter, as well as soil hydrothermal conditions, which influenced the pools and mineralisation rates of SOM in EFIMOD3.

We compared the values of autotrophic (tree roots) and heterotrophic (bacteria and non-mycorrhizal fungi) soil respiration calculated by EFIMOD3 and ectomycorrhizal respiration calculated by the EMR model. In addition, we compared the amount of available nitrogen from SOM decomposition calculated by EFIMOD3 and the amount of ammonium nitrogen (N-NH₄) calculated by the EMR model as a pool of micro- and mesofauna excreta.

All calculations were performed on a simulation grid of 100×100 m, divided into 10,000 cells of 1×1 m, for a period of 5 years with a daily time step. This allowed us to assess the possible variation of the indicators (components of soil respiration and nitrogen formed in FWs) both in time and in space. The first (temporal) is reflected in the variation of values between different dates of specific years or terms of the growing season. The second (spatial) is related to the distribution of the fine roots of the trees over the simulation plot and is reflected in the ranges between minimum and maximum values for each specific day of the simulation period.

The calculation of the amount of root exudates for the EMR model was carried out similarly to the validation and the same C:N ratio of RE (equal to 40) was used. Temperature and moisture dynamics of the organic layer and mineral soil were simulated using the observed data on air temperature and humidity and precipitation for 2016–2020 measured at the EMEP station “Danki” located in the reserve [86,87].

The results of the estimation of the 5-year cycle of soil C-CO₂ emission fluxes (autotrophic, mycorrhizal and heterotrophic) are presented in Fig. 3. All calculated fluxes show a remarkable temporal and spatial variation. The maximum of all fluxes corresponds to the middle of the growing season (July). The main contribution to the total soil respiration is made by heterotrophic respiration, whose share varies from 45% in summer to 95% in winter, against the background of a corresponding decrease in the average total flux of C-CO₂ from 3.5 to 0.4 g m⁻² day⁻¹ from July to December (Fig. 4).

In the first half of the growing season (April–June), ectomycorrhizal respiration is estimated by the EMR model at an average of 0.25–0.3 g m⁻² day⁻¹ (in terms of C), and at 0.35–0.4 g m⁻² day⁻¹ in the second half (July–September). The dynamics of ectomycorrhizal respiration rates correlate with autotrophic (stand) respiration rates, which follows from the general statements of the EMR model. The 5-year average estimates for the growing season (April–October) show that EM respiration can reach 35–60% of the autotrophic C-CO₂ flux and up to 15% of the respiration of soil bacteria and non-mycorrhizal fungi.

The results of the comparison of the N pools available to plants calculated by the EFIMOD3 and EMR models are shown in Fig. 5. The obtained patterns of dynamics of these indicators are similar to the dynamics of C-CO₂ fluxes mentioned above. The amount of ammonium nitrogen (N-NH₄) formed in soil food webs during the summer months (June–August) is estimated by the EMR model to be 5–10 mg m⁻² day⁻¹ on average (in terms of N) with a range of variation from close to zero values to 40 mg m⁻² day⁻¹. The 5-year average estimates for the growing season (April–October) show that the pool of ammonium N formed as by-products of micro- and mesofauna activity in soil food webs can reach 10–12% of the total pool of available N formed during SOM mineralisation (Fig. 6). It should be noted that another noticeable fraction of the EMR-derived N comes to trees via the decomposition of the labile SOM pool of faunal excreta. In addition, we found that the total C of the EMM-derived casts of micro- and mesofauna of the food webs varied between 0.1 and 0.3 t ha⁻¹ for 3 months of summer.

Figure 3. Simulated soil fluxes of C-CO₂ for the simulation plot (10,000 cells, light grey): Rh is the heterotrophic respiration, Ra is the autotrophic respiration (tree roots), Rm is the ectomycorrhizal respiration. Bold lines indicate mean values over the simulation plot.

Figure 3. Simulated soil fluxes of C-CO₂ for the simulation plot (10,000 cells, light grey): Rh is the heterotrophic respiration, Ra is the autotrophic respiration (tree roots), Rm is the ectomycorrhizal respiration. Bold lines indicate mean values over the simulation plot.

Figure 4. Simulated intra-annual dynamics of total soil respiration (monthly averages for 5 years) and the ratio of mycorrhizal (Rm), autotrophic (Ra) and heterotrophic (Rh) fluxes of C-CO₂.

Figure 4. Simulated intra-annual dynamics of total soil respiration (monthly averages for 5 years) and the ratio of mycorrhizal (Rm), autotrophic (Ra) and heterotrophic (Rh) fluxes of C-CO₂.

Figure 5. Simulated dynamics and variation of the total pool of available nitrogen calculated by EFIMOD3 (Total) and N-NH₄ calculated by the EMR model (Ectomycorrhizal FW) for the simulation plot (10,000 cells).

Figure 5. Simulated dynamics and variation of the total pool of available nitrogen calculated by EFIMOD3 (Total) and N-NH₄ calculated by the EMR model (Ectomycorrhizal FW) for the simulation plot (10,000 cells).

Figure 6. Intra-annual dynamics of the total pool of available N in the soil (monthly averages for 5 years) and the ratio of available (Total) and ammonium N fluxes (Ectomycorrhizal FW) calculated by EFIMOD3 and EMR models, respectively.

Figure 6. Intra-annual dynamics of the total pool of available N in the soil (monthly averages for 5 years) and the ratio of available (Total) and ammonium N fluxes (Ectomycorrhizal FW) calculated by EFIMOD3 and EMR models, respectively.

6. Discussion

Currently, there are two main concepts of the symbiotic interaction between plant and EMR [28,88]: (i) mutualistic with direct exchange of resources between plant and fungi (“carbon by nitrogen”), and (ii) indirect effect of EMR on the host plant with a specific pathway of N return from EMM to plant roots via EMM decomposition. The proposed model is an attempt to quantify the second approach, with EMR as an obligate biotroph that consumes all C, N and, moreover, SOM-derived additional “mined” N from RE for its own growth only, thus actually acting as a commensal. N return to the plant occurs through the fast cycle of EMM transformation and decomposition by the biota of the soil food webs, including the activity of microorganisms, micro- and mesofauna [40]. The majority of recent EMR models describe EMM growth without this N feedback to the plant, with the exception of [4], which calculates different pathways of N acquisition.

The structure of the proposed model has been developed for integration into the forest ecosystem model EFIMOD3. Therefore, the model is focussed on N production for tree growth and C turnover and sequestration in the soil. Both tasks were solved by incorporating food webs into the model structure by establishing a previously missing N-rich “fresh and free” SOM pool of faunal excreta with a solid phase organic matter not yet bound to R₂O₃, Ca and clay minerals.

The sensitivity analysis revealed a leading role of C:N ratio in root exudates, EMM and soil fauna biomass. In principle, the same effect has already been described in the simulation model of the rhizosphere priming effect [21].

In the model validation, the underestimation of simulated EMM respiration (compared to observed values) in sparse tree plantations (aged 11 years), where soil patches without tree roots may occur, can be a result of an overestimation of respiration for the entire experimental plot. The other reasons could be both high metabolic activity of young trees and differences between actual and simulated soil properties on the study plot, as the soil was assumed to be the same in all simulation scenarios. The overestimation of simulated respiration in the 25-years-old stand was apparently related to an unaccounted effect of a colder microclimate in dense stands of this age [89,90]. The overestimation of simulated respiration in early spring and late autumn may be related to the unaccounted influence of short daylight hours.

The results of the test simulations in a frame of the EFIMOD3 system also showed the same order of magnitude as the experimental data [91,92,93].

The excess carbon of root exudates forces EMM to mine N for biomass growth from protected stable SOM with a low C:N ratio. The model allows the calculation of the amount of N mined for the growth of short-living EMM. It provides an opportunity for a more accurate assessment of the EMR contribution to the N budget of the forest ecosystem.

EMM produces a lot of DOC during N mining [24]. This occurs because EMM first depolymerises large polymer molecules into soluble OM with low molecular weight. Then, EMM extracts and uses only N-rich OM, and does not consume the rest as other saprotrophic microorganisms do. The proposed model is, to our best knowledge the first to include this specific EMR-derived DOC pool in the calculation of SOM dynamics.

A significant detail is that the concept prevailing now assumes that microbial-derived DOC is a significant source of SOC sequestration on minerals surfaces [94]. However, soil microorganisms produce DOC for their own use at a first stage of nutrient consumption, i.e. during SOM decomposition [41]. Therefore, microorganisms have a limited capacity for free DOC production and consequently for SOM formation (in a process of SOM decomposition). Logically, DOC in forest soils can be mainly represented by the mycorrhiza-derived substances obtained from the specific SOM depolymerisation during N mining.

Furthermore, if the activity of food webs fauna is excluded from the processes of EMM transformation and decomposition (as in the vast majority of biogeochemical and microbiological SOM decomposition models), the capacity of EMR for SOM sequestration will be significantly reduced. In this case, the intermediate product of decomposition will be the small amount of DOC that is bound only by clay minerals. Soil minerals have a limited sorption capacity with the possibility of “carbon saturation” [94] compared to the unlimited potential of the coprolite formation (excretes, casts) of micro- and mesofauna. These are micromorphologically visible structures [95,96] that can be further protected in soil by R₂O₃ or Ca²⁺.

The calculated total production of faunal casts (in terms of С) can reach 0.1–0.3 t ha⁻¹ for 3 summer months in a temperate broadleaved-coniferous forest, as was found in the model experiments. This value is only three times lower than the cast production by earthworms, which varies from 0.1 to 1.0 t ha⁻¹ (in terms of С) during the growing season (adopted from [97,98]). The casts from EMM-derived litter are processed by the gut microbiomes of soil micro- and mesofauna. Apparently, they play the same role as the “carbon trap” mediated by earthworms [99].

EMM was rather rapidly consumed by soil biota [45,100], with C being returned to SOM as the above-mentioned solid-phase pool of faunal metabolites. They are used for organic C sequestration and also to increase the amount of N available to roots in the case of this young N-rich organic matter mineralisation. The balance between these two processes depends on biological and environmental factors. This cycle has not yet been considered in mycorrhizal models and is therefore a task for future research.

The formation of available ammonia N in food webs is an additional factor in the positive effect of EMR with the significant increase in tree growth. This occurs because some excess N is produced during the consumption of microorganisms (including EMM) by soil fauna, as a significant proportion of food C is used for faunal respiration [48,49,51].

From an ecological point of view, the results of this model validation and test prove that the representation of mycorrhiza as a “commensalist” can be a main pattern to describe the functioning of EMR in forest ecosystems. EMR “works for itself”, but its positive effect on tree growth can be explained by the rapid transformation and decomposition of EMM by all soil biota, including microorganisms and fauna of soil food webs. In other words, the mycorrhizal symbiosis can be represented as a triad of biological actors: “roots – mycorrhiza – all soil biota” [40]. They are functionally and perhaps evolutionally linked in the forest ecosystem. New data are needed to specify the processes of N mining, DOC formation, food web activity and natural mortality of EMM.

In the described model, there are no pathways for direct resource exchange between host plant and ectomycorrhiza (carbon from plant for nitrogen from ectomycorrhiza). Nevertheless, the satisfactory results of the model validation show that the positive effect of EMR on the plant with the production of additional nitrogen and SOM sequestration by faunal casts can only be achieved during the fast cycle of transformation and decomposition of EMM litter by all soil biota.

7. Conclusion

The proposed model of ectomycorrhiza in forest ecosystems is designed to be incorporated into the EFIMOD3 forest ecosystem model and possibly other forest models. Therefore, the role of EMR in the nitrogen budget and SOM dynamics in forest ecosystems is emphasized. The model is based on the concept of an indirect positive effect of EMR on plant growth due to the rapid decomposition of short-living EMM, which represents a dominant belowground litter pool in forest soils. The model allows the assessment of EMM growth, N mining and release for plant nutrition, and SOM formation. The growth of EMM is determined by the input and the C:N ratio of root exudates. Mycorrhizal N mining of SOM is calculated in the model only for EMM own use, with a significant amount of dissolved organic matter as a by-product. The decomposition of EMM by all soil biota of food webs allowed the estimation of (i) the available N-NH₄ that can be directly used by plants and (ii) the solid-phase by-product of faunal metabolism, which is a pool of fresh N-rich organic matter that has never been represented in EMR models. The validation and testing of this model demonstrated the adequacy of the theoretical background used in this model.