Phosphorus regulates ectomycorrhizal fungi biomass production in a Norway spruce forest

Ectomycorrhizal fungi (EMF) are important components of soil microbial communities, and EMF biomass can potentially increase carbon (C) stocks by accumulating in the soils as necromass and producing recalcitrant structures. EMF growth depends on the C allocated belowground by the host trees, and the nutrient limitation on tree growth is expected to influence this allocation. Therefore, studying EMF production and understanding the factors that regulates it in natural soils are important to understand C cycling in forests.

Fungal mycelium collected from ingrowth mesh bags is commonly used to estimate EMF biomass, but these measurements might not reflect the total EMF production since turnover rates of the hyphae are not considered. Here we estimated EMF production and turnover in response to P fertilization (applied as superphosphate) in a Norway spruce forest where nitrogen (N) deposition has resulted in phosphorus (P) limitation of plant production by using a combination of mesh bags with different incubation periods and with Bayesian inferences. To test how localized patches of N and P influence EMF production and turnover we amended some bags with a nitrogen source (methylene urea) or P source (apatite). Additionally, the Bayesian model tested the effect of seasonality (time of mesh-bag harvesting) on EMF production and turnover.

We found that turnover of EMF was not affected by P fertilization or mesh-bag amendment. P fertilization had a negative effect on EMF production in all the mesh-bag amendments, suggesting a reduced belowground C allocation to the EMF when P limitation is alleviated. Apatite amendment significantly increased EMF biomass production in comparison with the pure quartz bags in the control plots but not in the P-fertilized plots. This indicates that P-rich patches enhance EMF production in P-limited forests, but not when P is not limiting. Urea amendment had a generally positive effect on EMF production, but this was significantly reduced by P fertilization, suggesting that a decrease in EMF production due to the alleviated P limitation will affect N foraging. Seasonality had a significant effect on EMF production, and the differences registered between the treatments were higher during the warmer months and disappeared at the end of the growing season.

Many studies highlight the importance of N for regulating belowground C allocation to EMF in northern coniferous forests, but here we show that the P status of the forest can be equally important for belowground carbon allocation to EMF production in areas with high N deposition.

1 Introduction

In terrestrial ecosystems forest soils are important reservoirs for carbon (Falkowski et al., 2000). Boreal forests contribute approximately 50 % of the total forest carbon stock from which around 85 % is stored in the soil (Malhi et al., 1999). At least half of the carbon stock in boreal soils originates from belowground carbon allocation through roots (Clemmensen et al., 2013), and a large portion of boreal forest primary production is allocated belowground by the trees (Gill and Finzi, 2016). The carbon dynamics in forest soils are highly dependent on soil microbial communities that either enhance C losses by degrading organic matter or increase C stocks by immobilizing C (Clemmensen et al., 2013). Filamentous fungi forming mycorrhizal associations, for example, play an important role for C fluxes since some species have the capability to degrade a great variety of organic compounds, while others can contribute to soil organic matter formation by releasing exudates that promote soil aggregation (Rillig, 2010) or produce slowly decomposing and highly melanized hydrophobic tissues (Almeida et al., 2022). The effect of EMF on soil microbial communities might not be trivial since up to 20 % of the net primary production is allocated belowground to support the symbiosis (Hobbie, 2006). Therefore, ectomycorrhizal mycelium is expected to be a significant part of the soil fungal biomass, and its production and turnover play an important role in forest carbon cycling and organic matter formation (Ekblad et al., 2013). For that reason, the development of methods that allow us to quantify EMF growth in forests natural soils is of paramount importance (Fernandez, 2021).

Therefore, understanding the factors that regulate the growth rates of filamentous fungi like EMF is important to understand carbon dynamics in soils. Growth rates of free-living fungi from natural soils have been studied in the laboratory by measuring labeled acetate incorporated in the fungal membrane component ergosterol (Zheng et al., 2022; Rousk and Bååth, 2007) or labeled water incorporated into DNA (Schwartz et al., 2016). Quantifying growth (production) of EMF natural communities, on the other hand, is more complicated since EMF is dependent on plant roots (Smith and Read, 2008) and such measurements must be performed when the fungi are living in symbiosis. Many studies have attempted to quantify EMF production in situ in forest soils by using ingrowth mesh bags and fungal biomarkers like ergosterol phospholipid fatty acids (PLFAs) (Wallander et al., 2013). In those studies, EMF production has been estimated based on the standing fungal biomass measured in mesh bags after a specific time of incubation in the soil (Ekblad et al., 2013; Wallander et al., 2013, 2001). However, the standing biomass does not necessary reflect growth since the standing biomass is the result of the interaction between fungal growth and the residence time of the fungal mycelium in the mesh bag (Ekblad et al., 2016). In order to overcome these shortcomings, some studies have estimated fungal production and mycelium turnover by repeated harvests of mycelial mesh bags, applying ergosterol as a marker of mycelial biomass, and mathematical models to estimate the production and turnover of EMF mycelium biomass (Hagenbo et al., 2021, 2017) or, combined with analyses of chitin, to enable estimates of production and turnover of both biomass and necromass (Ekblad et al., 2016). In these studies, the standing biomass and necromass were analyzed in bags incubated over periods varying in length, combining several shorter periods, one after the other, with overlapping longer periods. Common assumptions in these studies were that EMF growth occurs at a constant rate and that biomass and necromass were lost at constant exponential rates (Ekblad et al., 2016).

By using this approach, Ekblad et al. (2016) tested the effect of nitrogen (N) fertilization on EMF turnover and growth in a Pinus taeda forest. They reported that fertilization significantly decreased both EMF standing biomass and growth, but turnover rates of biomass and necromass were not affected. It was suggested that the decrease in EMF growth was regulated by changes in carbon allocation as a result of an increase in soil fertility. These results are in line with evidence indicating that the relative amount of carbon allocated to EMF is sensitive to plant nutrient status and soil fertility (Gill and Finzi, 2016). Thus, in boreal forests where N is the nutrient that limits tree growth (Högberg et al., 2017), high amounts of carbon are invested belowground to support ectomycorrhizal symbiosis to facilitate N uptake (Gill and Finzi, 2016).

The role of N as a limiting nutrient in high-latitude forested ecosystems and its effect on EMF are well known and have been described in several studies (Binkley and Högberg, 2016; Hedwall et al., 2013; Gill and Finzi, 2016). However, there is some evidence suggesting that anthropogenic N deposition can potentially change the forest's nutrient requirements and push the system toward phosphorus (P) limitation (Almeida et al., 2019; Jonard et al., 2015; Talkner et al., 2015; Prietzel et al., 2020; Du et al., 2021). In fact, in a region with high N deposition in southwestern Sweden, Almeida et al. (2019) reported that P fertilization had a stronger effect on tree growth than N fertilization, subverting the expectation that N is the main nutrient regulating plant growth in northern forests. The effect of the transition from N to P limitation on the belowground C allocation and EMF growth has not been studied in natural soils, but P deficiency is expected to increase EMF biomass to improve P foraging and uptake (Rosenstock et al., 2016; Ekblad et al., 1995; Wallander and Nylund, 1992). In a field study, Rosenstock et al. (2016) reported an increase in root and standing biomass in a Norway spruce (Picea alba) forest limited by P compared to forests with sufficient P. In the field study performed by Almeida et al. (2019), however, no effect on EMF standing biomass was found in mesh bags incubated for 133 d. Yet, since only the standing biomass was measured and the turnover rates and production were not estimated, we cannot exclude the possibility that P fertilization had an effect on EMF production, an effect that cannot be detected by studying the standing biomass alone.

In this study, we aimed to improve our understanding of EMF production and turnover in natural soils by testing how fungal biomass collected from ingrowth mesh bags is affected when P is limiting tree growth. In the forest described by Almeida et al. (2019) we estimated fungal production (which is assumed to be dominated by EMF production) and turnover using the mathematical model of Ekblad et al. (2016) with Bayesian inferences. Our first hypothesis was that P fertilization will decrease EMF biomass production in this P-limited forest as a result of the limitation being alleviated.

In addition, because EMF growth is subsidized by the host, in exchange for N and P, EMF production in the mesh bags should be affected by the nutrients found at the hyphal front. Indeed, EMF biomass in P-poor forests is stimulated around localized patches of the P-rich mineral apatite (Rosenstock et al., 2016; Berner et al., 2012; Hagerberg et al., 2003). Therefore, besides purely sand-filled mesh bags, we incubated mesh bags amended with apatite or methylene urea (referred to as urea throughout the paper) in order to simulate soil N and P nutrient patches, respectively. We expected that the nutrient patches would increase EMF biomass production depending on fertilization. In particular, apatite amendment will increase EMF biomass production in the control plots but not in P-fertilized plots (second hypothesis), and urea amendment will increase EMF biomass production in the P-fertilized but not in the control plots (third hypothesis).

Finally, since belowground C allocation follows the three phenological cycles (Endrulat et al., 2016), EMF production is likely to vary with season, peaking in autumn (Hagerberg and Wallander, 2002; Wallander et al., 2001; Hagenbo et al., 2021), we performed a more extensive incubation scheme and more frequent harvests of bags than in Ekblad et al. (2016). This allowed us to not only test effects of treatments (P fertilization) and of mesh-bag amendments, but also to estimate possible seasonal effects. Therefore, our fourth hypothesis was that EMF biomass production will be higher in autumn than in summer.

2 Materials and methods

2.1 Field site and fertilization treatments

This study was performed at Tönnersjöheden forestry research station (56^∘41^′ N, 13^∘6^′ E, 80 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$) with a mean annual temperature of 6.4 ^∘C and a mean annual precipitation of 1064 mm (Högberg et al., 2013). Soils are podzols developed in a glaciofluvial parent material with a pH (in H₂O) of 4.05 and a $\mathrm{C}/\mathrm{N}$ of 25.1 in the mor layer (Hansson, 2011; Högberg et al., 2013). The forests consist of managed Norway spruce (Picea abies) planted on former pastureland in 1979. The site is in southwestern Sweden with an N deposition of 14.5 $\mathrm{kg}\mathrm{N}{\mathrm{ha}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Rosenqvist et al., 2007), which is high in comparison with most other forests in the country (Akselsson et al., 2010; Högberg et al., 2013) and exceed the N critical loads in which negative changes in the function and composition of an ecosystem are expected (Kuylenstierna et al., 1998; Pardo et al., 2011; Pihl Karlsson et al., 2017).

The total experimental area comprised 2.1 ha. The experiment consisted of six plots (30–40 m × 25 m): three control and three P-fertilized plots. Since availability of P can be very low in soils, and to make sure there is enough available P for the trees to get a growth effect if P is limiting, an excess of P was added. Thus, the P-fertilized plots received 200 kg P ha⁻¹ of superphosphate (100 kg ha⁻¹ applied twice in September 2011 and July 2012).

2.2 Experimental design

To estimate EMF mycelial production, ingrowth mesh bags (Wallander et al., 2001) were incubated in the plots. The mesh bags were cylindrical, 2 cm wide and 10 cm long. They were made of 50 µm nylon mesh and filled with approximately 40 g of acid-washed quartz sand. The mesh size is used to exclude fine roots but allow mycelium hyphae colonization (see Wallander et al., 2001).

Three different amendments in the mesh bags were used: quartz-only (pure sand), apatite-amended (quartz and 1.5 %($w/w$) crushed apatite mineral with a grain size of 50 to < 250 um) and urea-amended (quartz and 0.5 %($w/w$) granulated methylene urea). The mesh bags were vertically installed into holes made with a soil corer (2 cm diameter) with the upper end of the bag level with the soil surface. The amount of apatite used is similar to the amount used by Rosenstock et al. (2016) and was chosen to provide enough mineral to the fungi to sustain growth for the whole duration of the experiment. The apatite grain size was chosen based on other studies showing that the respiration of EMF colonizing the mineral became higher when the grain size was small compared to when it was larger (Leake et al., 2008). The 50 µm mesh of the used bags is small enough to avoid losses of the apatite material.

More apatite than urea was provided since apatite is a more recalcitrant source and has a relatively low content of P (18 %) compared to the N content (42 %) of methylene urea.

To calculate turnover rates and biomass production as done by Ekblad et al. (2016), sequential mesh-bag incubations were performed for a 5-month period starting in July 2015 and ending in November 2015 (150 d). This period was chosen since it covered the expected productive period for EMF growth from summer to late autumn. The mesh bags were incubated for variable periods of time (30, 60, 90, 120 or 150 d; Fig. 1). The shorter period (30 d) was chosen as it was reported that the mean residence time of the fungal biomass in mesh bags in a previous study was around 28 d (Ekblad et al., 2016).

Figure 1 Overview of the incubation design. Different color bars represent the incubation time periods: yellow corresponds to 30 d, light green to 60 d, dark green to 90 d, purple to 120 d and blue to 150 d of incubation. The arrows represent the points in time when the same holes from the previous incubation were used to incubate the next set of mesh bags.

There were five different 30 d incubation periods: four 60 d incubation periods each overlapping with two 30 d incubation periods, three 90 d incubation periods each overlapping with three 30 d incubation periods, two 120 d incubation periods each overlapping with four 30 d incubation periods and one 150 d incubation period overlapping with all 30 d incubation periods.

The bags incubated over 30 d were incubated sequentially, and when one set of bags was collected, a new set of bags was directly installed using the same holes as the ones just emptied (Fig. 1).

In each plot, a quartz-only mesh bag for each of the incubation periods described above was placed along a 15 m long transect. The distance between each mesh bag was approximately 1.5 m. The apatite-amended and urea-amended bags were placed 10 cm (perpendicular to the long transect) at each side of the quartz-only mesh bags. Three 15 m long transects were done to have three sub-replicates (for each set of bags) that were pooled before further analysis to give one sample from each incubation period and amendment (quartz-only, apatite and urea) per plot.

Each incubation period consisted of 54 mesh bags (two treatments $\mathrm{C}/\mathrm{P}$, three replicated plots, three sub-replicates, three amendments: 2 × 3 × 3 × 3 = 54). In total, 810 mesh bags were installed and collected according to their incubation period.

Upon harvest, the mesh bags were kept in an icebox until arrival at the laboratory, up to 10 h in summer and up to 7 h in autumn, where they were stored at −20 ^∘C.

The fungal cell membrane compound ergosterol was used as a proxy for fungal biomass to infer EMF growth as this compound has been used as a quantitative estimation of living EMF biomass (Ekblad et al., 2016; Wallander et al., 2013; Hagenbo et al., 2018, 2021). Other markers like fungal DNA sequences can be also used to estimate EMF abundance, although this method relies on amplicon relative abundances and is semiquantitative. Ergosterol was extracted and measured from 5 g of the pooled samples as per Bahr et al. (2013) using high-pressure liquid chromatography (auto-sampler L2130 with UV detector L2400 by Hitachi, Japan). It was assumed that after incubation in the soil the mesh-bag contents were dominated by EMF as has been shown by metabarcoding (Almeida et al., 2019; Rosenstock et al., 2016; Berner et al., 2012; Wallander et al., 2010; Hedh et al., 2008) and isotopic studies (Wallander et al., 2001). Therefore, the fungal biomass collected was expected to be of EMF origin.

2.3 Mathematical models

The turnover rates and EMF biomass production were estimated by applying the mathematical model used in Ekblad et al. (2016). In this paper, however, the mathematical model was tested under two assumptions: EMF production was dependent on the treatments alone (Model 1), or EMF production was dependent on treatments and sampling season (Model 2), allowing us to test for the interactions between treatment and seasonal effects.

Model 1

This model works under the assumptions that EMF production occurs at a constant rate (in units of biomass over time) and biomass is lost at a rate which depends on the number of microbial cells times a linear mortality coefficient (see Hagenbo et al., 2017; Ekblad et al., 2016). EMF growth is assumed to be constant because growth is limited by transport of nutrients from outside, which is its kinetic limitation.

Biomass and necromass losses are both just organic matter decomposition, so we expect them to follow the same kinetics as any other organic matter (linear first derivative and exponential integrated over time).

Briefly, the sum of the biomass during two sequential short incubation periods is expected to exceed the biomass in an overlapping longer incubation period due to an on average older mycelium and hence larger turnover in bags with a longer incubation period.

The model in its differential form is defined as

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\mathrm{d}B}{\mathrm{d}t}}=P-\mathit{\mu }\cdot B,\end{array}$$

where P is the production of new mycelium (in mass units), B is the mycelium biomass (also in mass units) and μ represent the mortality, the fraction dying over a specified time period (adimensional). This equation is solved over time as

$$\begin{array}{}\text{(2)}& B\left(t\right)={\displaystyle \frac{P_k}{{\mathit{\mu }}_k}}\cdot (\mathrm{1}-e^{{\mathit{\mu }}_{k}t}).\end{array}$$

In our case we assumed that both P_k and μ_k are influenced by the fertilization treatments, denoted here by k, and we therefore assigned a specific (unknown) P and μ to each treatment in the Bayesian model.

Model 2

Equation (2) has been utilized in other publications (Hagenbo et al., 2021, 2017; Ekblad et al., 2016), and one of the main assumptions of this model is that EMF production occurs at a constant rate. However, EMF production can vary depending on the time of the year (Coutts and Nicoll, 1990; Hagenbo et al., 2021), so we tested a modification of the model by introducing an additional degree of freedom into the model represented by the term β_k,j, dependent on sampling seasons (j) and their interactions with treatments (k) so that the calibration can apply a correction for seasonality to each treatment (independent from the other treatments). When ${\mathit{\beta }}_{k,j}=\mathrm{1}$, then the model is equivalent to what described in Eqs. (1) and (2). We utilized this model to decompose P into two components, defining a new term P^′.

$$\begin{array}{}\text{(3)}& P_{k,j}^{\prime }=P{\mathrm{0}}_k\cdot {\mathit{\beta }}_{k,j}\end{array}$$

$P_{k,j}^{\prime }$ corresponds to P_k (if the distributions were perfectly symmetric the average for P and P^′ should converge to the same value), but the predicted biomass production is now the result of the interactions between sampling season and treatments.

Equation (3) is then substituted into Eq. (2) by substituting P with P^′. The resulting model is equivalent to the one described by Eq. (2) for certain parameter combinations and describes the same curve. The only difference is that now two components are used to decompose the variance explained by the calibrated model in two separate terms: P0_k, which expresses the production variable with treatments only (k), and β_k,j, which expresses the effects of seasonality and their interactions with treatments. P0_k is now equivalent to the production normalized by the seasonality effect $\frac{P_{k,j}^{\prime }}{{\mathit{\beta }}_{k,j}}$. By letting P0_k and β_k,j vary independently (therefore describing each point as a combination of k and j) we avoid making any strong assumption on the effect of seasonality (since we are not imposing a parametric function of time to describe it, but we leave it free to vary for each time point) or on its interactions with treatments (which are still free to vary depending on the treatment), while on the other end we maximize the information we can extract from the data by representing the interactions between the terms in one single model calibration. If we had instead relied on fully independent calibrations within each subset of seasons × treatments we would have had to divide the data into j×k subsets for which we would calibrate each model parameter independently, limiting each calibration to a smaller number of samples.

2.4 The calibration

The model was calibrated within a formal Bayesian framework, developed with the Stan toolbox (Stan Development Team, 2021). This approach is based on a numerical implementation of Bayesian statistics, which allows for a continuous update of the knowledge while new data are developed, based on stochastic principles (through a modification of the Metropolis–Hastings sampler). The main assets of the method are that (a) we can integrate and utilize previous information in the calibration, defining it as prior probability distributions of model parameters (from now on, “priors”), and (b) such information is combined with the statistical information contained in the data to determine the posterior distributions of model parameters and consequently predictions; such a distribution is nonparametric (so not assuming any specific shape but determined only by the available information). The methodology is therefore extremely useful to combine multiple sources of information and very valuable when information is scarce, and at the same time it is quite robust given that it estimates detailed posterior probability distributions (which can be examined closely).

In our case the methodology allows us to draw information from previous studies.

In particular, we used information from an EMF production study in a conifer forest by Hagenbo et al. (2017). This information is considered probabilistically. It does add information to our final results (our posterior distributions), but such information is combined with the information contained in our data. The chosen statistical approach updates the old information with new data, and old and new information can be therefore compared.

We calibrated both a model with only Eq. (2) (so considering only treatment effects; Model 1) and one considering Eqs. (2) and (3) (considering treatments × seasonality effects; Model 2).

Priors for P_k and μ_k were derived from the mean EMF biomass production and turnover for a forest of similar age as the forest in the current study and estimated by Hagenbo et al. (2017) after unit conversion. Both priors were expressed as normal distributions with deviation prudentially estimated as 25 % of the mean (please note that this does not mean that the prior was limited within this range due to the tails of the normal distributions).

P_k was expressed as

$$P_k\sim N(\mathrm{0.099},\mathrm{0.099}\cdot \mathrm{0.25}),$$

while μ_k is

$${\mathit{\mu }}_k\sim N(\mathrm{0.009},\mathrm{0.009}\cdot \mathrm{0.25}).$$

The Bayesian system was run considering one independent P_k and μ_k for each treatment.

When we also considered Eq. (3), priors for P0_k were defined as the priors for P_k, while priors for β_j were set as uniform between 0 and 5.

$$\mathit{\beta }\sim U(\mathrm{0},\mathrm{5})$$

Please note that β_j=1 means no seasonality effect, and β_j=5 means a 5-fold increase in production due to seasonality, while β_j=0 means a complete halt of production due to a seasonal effect.

2.5 Statistical analysis and probability distribution comparisons

The standing biomass data were tested for homogeneity of variances and normal distribution using Levene and Shapiro–Wilk tests, respectively. Analysis of variance (ANOVA), Tukey's post hoc test and Dunn analyses were performed on the data to check for statistical differences between the fertilization treatments and mesh-bag amendments. The Levene and Shapiro–Wilk tests, as well as ANOVA and Dunn analyses, were done by using R (R Core Team, 2014).

The stochastic approach of the Bayesian method produces Markov chain Monte Carlo (MCMC) that represents a probability distribution with as many discrete parameter values as iterations in the chains (in our case 10 independent chains of 10 000 iterations, so a total of 100 000 iterations), with a histogram that approximates a continuous distribution (probability distribution). Thus, the predicted EMF production and turnover for each treatment (fertilization regime and mesh-bag amendment) are represented by a probability distribution.

The means of the probability distributions were calculated and the highest density intervals of the estimated parameters were interpreted as confidence intervals at 95 % and 90 % (Kruschke and Liddell, 2018). To test the significance of the treatments (fertilization regime, mesh-bag amendment and season), the confidence intervals of the probability distributions were compared.

3 Results

3.1 Mycelial standing biomass

The standing biomass of mycelia in the mesh bags was significantly affected by incubation period (time of the year) (Kruskal–Wallis, p < 0.0001, X² = 116.4).

Figure 2 Standing EMF biomass in the mesh bags that have been in the soil for a period of 2 months (a) and a period of 1 month (b). The error bars represent the standard error of the mean (n=3). Lowercase letters represent statistically significant (P < 0.05) differences between the incubation periods according to Dunn's test.

Figure 3 Standing EMF biomass in the three mesh-bag amendments (quartz-only, apatite and urea) and in the control plots (red symbols) as well as the P-fertilized plots (blue symbols) and control plots during different incubation times (30, 60, 90, 120 and 150 d). The error bars represent the standard error of the mean.

Biomass in 1-month incubation mesh bags from July, August and September was significantly higher than the biomass collected in October and November for both control plots and P-fertilized plots (Dunn's test, p < 0.001, X² = 26.1) (Fig. 2). Biomass in 2-month incubation mesh bags from July–August and August–September was significantly higher than the biomass collected in September–October and October–November for both control plots and P-fertilized plots (Dunn's test, p < 0.001, X² = 27.7; Fig. 2). Fertilization significantly affected the standing biomass in the quartz-only, apatite and urea-amended mesh bags (Kruskal–Wallis, p < 0.05 and X² = 6.5, p < 0.0001 and X² = 18, p < 0.0001 and X² = 15.5, respectively). Phosphorus fertilization reduced the standing biomass at all the incubation times (numbers of incubation days) for the apatite- and the urea-amended mesh bags (Fig. 3). Apatite amendment significantly increased the standing biomass in comparison with the quartz-only bags in the control plots after 60 and 150 d of incubation (Dunn's test, p < 0.05 and X² = 18, p < 0.05 and X² = 11.2, respectively), and the effect of apatite was stronger after 150 d of incubation when the biomass in the apatite bags was on average 3-fold higher than the biomass in the quartz-only bags. Apatite amendment did not increase biomass in the P-fertilized plots at any incubation time, while urea amendment increased biomass at most of the incubation times and for both C- and P-fertilized plots (Dunn's test, p < 0.05) (Fig. 3).

Table 1 Mean of the EMF production in different treatments (P_k) estimated by Model 1. The highest density intervals (HDIs; Kruschke and Liddell, 2018) represent the boundaries of each estimate at different degrees of confidence.

3.2 EMF production and turnover rates (Model 1)

The predicted EMF biomass production varied between the P-fertilized plots and the control plots and between the mesh-bag amendments (Fig. 4a). P fertilization significantly decreased EMF production in all the mesh-bag amendments (urea and apatite and quartz-only) (Table 1). In the P-fertilized plots the EMF production was reduced to a third in the apatite and quartz-only bags in comparison with the prior used to set the model (0.099 g m² d⁻¹). P fertilization caused a reduction on average of 43 % in the quartz bags, 60 % in the apatite bags and 39 % in the urea bags in comparison with the control plots.

Figure 4 (a) Probability distribution of the predicted EMF biomass production (P_k) (g m² d⁻¹) for the different fertilizer treatments (control and P fertilization) and mesh-bag amendments (quartz-only, apatite and urea). (b) Probability distribution of the turnover rates (d⁻¹) for the different fertilizer treatments (control and P fertilization) and mesh-bag amendments (quartz-only, apatite or urea).

The mesh bags amended with urea had the highest predicted biomass production in both control and P-fertilized plots (Fig. 4). Relative to the quartz bags, the urea amendment doubled the production in both fertilizer treatments. The apatite amendment, in contrast, gave no significant change in production relative to the quartz bags in the P-fertilized plots, while a 35 % increase was found relative to the quartz bags in the control plots (Table 1).

According to the mathematical modeling, the biomass turnover rates were not affected by P fertilization or mesh-bag amendment (Fig. 4b).

Figure 5 Seasonality effect on biomass production expressed by the β parameter for the different months of the growing season.

3.3 Seasonal effect (Model 2)

The effect of seasonality as described by β had a positive effect on the predicted EMF production, and this effect was highest in July and decreased over time. Moreover, the effect of β on EMF production differed depending on the fertilization and on the mesh-bag amendment (Fig. 5).

For example, in July the model suggests a seasonal effect increasing the predicted EMF production by up to 5 times in the quartz mesh bags from the P-fertilized plots and up to 2.5 times in the urea mesh bags in the control plots in comparison with the apatite bags from the P-fertilized plots where season had no effect on EMF production. The positive effect of sampling season on the EMF production, as identified by the model, decreased in general with time, and at the end of the growing season (October and November) β had the same effect on all the samples independently from the treatment (fertilization and mesh-bag amendment).

Figure 6 Probability distribution of $P_k^{\prime }$ (g m² d⁻¹) for the different months of the growing season.

Even though the β probability distributions of the different treatments were not significantly different, the effect of the season on biomass production was important, and when we decompose EMF production by seasonality ($P_k^{\prime }$), the differences in EMF production between P fertilized and control plots and between the mesh-bag amendments are present only early in the season (July, August) and disappear in September, October and November (Fig. 6).

4 Discussion

4.1 Effect of P fertilization on EMF biomass production and turnover

In support of our first hypothesis, EMF biomass production declined in response to P fertilization in all mesh-bag amendments (Fig. 4a). This reduction in EMF production was not trivial, and P fertilization decreased the predicted EMF production to a third in comparison with the EMF production of a forest of similar age estimated by Hagenbo et al. (2017) (0.099 g m² d⁻¹). These results contrast with those of Almeida et al. (2019), who tested the effect of P fertilization on the EMF standing biomass in the same plots as in the present study. This contrast does not depend on variation in turnover rates between control and P-fertilized plots since mortality was not significantly affected by fertilization as shown indirectly in the current results. In the present study, P fertilization had a negative effect on the EMF standing biomass in most of the incubation periods (Fig. 3). Thus, the standing biomass of one given incubation time might not truly reflect the effect of fertilization on EMF growth. The use of the sequential incubation method and the mathematical model allowed us to have a more robust estimate of the effect of P fertilization on the extramatrical mycelium in this forest.

Fertilization experiments have largely been used to evaluate the effect of soil fertility and nutrient status of the trees on carbon allocation and EMF production (Bahr et al., 2015; Ekblad et al., 2013). However, studies on the effect of nutrient additions on EMF in boreal forests have predominantly focused on N fertilization (Leppälammi-Kujansuu et al., 2013), probably because N is the most common limiting nutrient in boreal forests (Högberg et al., 2017). Therefore, the effects of P fertilization alone on boreal forests have not been widely tested despite evidence that the steep increase in anthropogenic C and N inputs can lead to unbalanced nutrition and push forested ecosystems to P limitation (Jonard et al., 2015; Peñuelas et al., 2013; Talkner et al., 2015; Prietzel et al., 2020; Du et al., 2021). Indeed, in the study performed by Almeida et al. (2019) in the same experimental plots as the current experiment, it was reported that P fertilization enhanced tree growth. Moreover, the authors reported that the foliar N:P ratios measured in the unfertilized control plots corresponded to suggested tipping points where the ecosystem shifts towards P limitation (see Suz et al., 2021, and van der Linde et al., 2018). The results of the current paper suggest that this shift is linked to changes in EMF growth as shown by the reduction of EMF biomass production when P fertilization alleviates the nutrient limitation. We propose that the decreased EMF production in the P-fertilized plots in our study is a result of a decrease in belowground C allocation due to reduced tree dependency on EMF for P foraging and acquisition. Fine root production and root tip colonization by EMF could be advisable as an independent second method to confirm that the decrease in EMF growth in the P-fertilized plots was an effect of reduced C allocation by the trees.

A potential decrease in belowground C allocation is also expected to alter EMF community composition, selecting for C efficient species when the ecosystem has crossed the nutritional tipping point thresholds (Suz et al., 2021). Indeed, in the soil EMF survey performed in the same experimental plots as the present study, Almeida et al. (2019) reported that the relative abundance of Tylospora asterophora was significantly increased after P fertilization. This species has been reported to extensively occupy ingrowth mesh bags while colonizing a relatively low number of tree root tips, which might suggest either a high C efficiency or lower turnover rates (Jörgenssen, 2021). The lack of difference in turnover rates between fertilized and unfertilized plots in the present study might suggest the former.

4.2 Effect of nutrient amendment on biomass production and turnover

Both nutrient amendments (urea and apatite) increased EMF production in comparison with the quartz-only mesh bags in the control plots. This is consistent with mesocosm experiments that have shown that when organic (Wallander and Pallon, 2005; Leake et al., 2001; Bending and Read, 1995) and mineral nutrient patches (Smits et al., 2012; Leake et al., 2008) are colonized by EMF, mycelial branching and proliferation increase to explore the nutrient patch. In support of our second hypothesis, apatite amendment increased EMF production in comparison with the quartz-only bags but only in the control plots. Our results are consistent with the view that trees in the control plots are P-limited and that they allocate more resources to the EMF when exploring a P source like apatite. When P limitation is alleviated by fertilization however, there is probably a decrease in C allocation to the root symbionts, which could cause the reduced EMF colonization in the apatite bags. This is supported by other studies reporting that apatite amendment increases EMF standing biomass in mesh bags under P-poor conditions (Rosenstock et al., 2016; Berner et al., 2012; Hedh et al., 2008; Hagerberg et al., 2003). In a fertilization study in nearby plots in the same forest, Bahr et al. (2015) showed that apatite addition stimulated EMF standing biomass in mesh bags in control and in N-fertilized plots, but when N was added in combination with P, on the other hand, no significant differences were found between apatite-amended and quartz-only bags. All together these results provide evidence that EMF growth is responsive to P nutrient patches, but this response is dependent on the P demand of the host.

Of the two nutrient amendments, urea had the highest effect on EMF growth in both the control and P-fertilized plots, partially confirming our third hypothesis. From a phytocentric point of view it could be expected that EMF growing on a P-rich source like apatite is rewarded with more C from the P-limited trees than EMF colonizing N bags. The stronger response of EMF growth to the N nutrient patches than to P nutrient patches in the P-limited control plots suggests that even though the forest is limited by P, N still has an important effect on the growth of the extramatrical mycelium. It is possible that P limitation results in a general increase in C allocation to the root symbionts and the C invested by the tree is delivered indiscriminately among its fungal symbionts, independently of the nutrient patch they are colonizing.

This is probably not surprising since N is needed by the fungus and plant alike, and in order to produce biomass to forage for P and enzymes to mineralize it, EMF requires N. Thus, N uptake can improve the P nutrition of the mycorrhizal system, and positive feedback between the plant and fungus might happen.

Despite the strong effect of N patches on EMF growth, P fertilization decreased growth in all mesh bags independent of the amendment. EMF communities in forests are diverse and composed of species with different abilities to mineralize the different nutrients present in the soils (Lilleskov et al., 2011). By amending the mesh bags with different nutrient types, EMF communities are selected depending on the nutrient added (Almeida et al., 2019; Rosenstock et al., 2016). The consistent effect of P fertilization on both nutrient patches and even in the barren quartz-only bags suggests that P fertilization affects growth of different EMF communities alike and reduces nutrient foraging for both N and P. This is consistent with the idea that alleviated P limitation results in a general decrease in C delivered to the roots and the mycorrhizal symbionts.

Previous studies on EMF growth have focused on EMF biomass collected from mesh bags filled with acid-washed sand (see Hagenbo et al., 2021; Hagenbo et al., 2017; Ekblad et al., 2016). However, since the quartz-only mesh bags are devoid of nutrients (except probably for dissolved organic material entering the bags during incubation), they might underestimate EMF production in soils. Moreover, in soils most N and P is heterogeneously distributed in nutrient patches (Hodge, 2006). For this reason, amending the mesh bags made it possible to imitate the soil nutrient conditions that influence EMF growth in forests and to understand how the nutrient regimes (both as inorganic nutrient fertilization and as nutrient patches) affect EMF production. In fact, the EMF growth in this study was influenced by both the nutrient at the hyphal front (N and P amendment) and by the C provided by the roots (as shown by the effect of P fertilization).

There were no differences in mycelium turnover between the different mesh-bag amendments. This contrasts with previous studies showing that the nature of a nutrient patch could also affect hyphal turnover (Ekblad et al., 2013; Jansa et al., 2011). Mineral substrates like feldspar have been shown to maintain EMF growth for up to 15 weeks (Rosling et al., 2004), while organic nutrient patches have been shown to sustain EMF growth for around 5 weeks (Bending and Read, 1995). Therefore, organic substrates like urea are expected to be quickly depleted in soils. As a result, the EMF hyphae are expected to autolyze and transfer the nutrients to other locations of the exploring mycelium faster than during the slow weathering of mineral substrates like apatite (Ekblad et al., 2013; Jansa et al., 2011). Therefore, it should be expected that the apatite bags show lower turnover rates than the urea bags. In the present study, however, we could not detect differences between the two nutrient patches. The material used to amend the urea mesh bags in this study is methyleneurea, which is a slow-N-release molecule. Thus, methylene urea is hydrolyzed to ammonium at a slower rate than the urea molecules (Högberg et al., 2021). Therefore, even if there is evidence that some EMF species can directly consume urea (Morel et al., 2008; Yamanaka, 1999), these slow-releasing nutrient sources might require a more persistent mycelium than other organic sources.

Additionally, previous mesocosm experiments have shown that when EMF mycelium grows on sand, longevity is enhanced in comparison with EMF growing on nutrient patches (Wallander and Pallon, 2005). Nutrient patches enhance growth and metabolic activity of EMF, which may enhance turnover rates. For example, Bidartondo et al. (2001) tested ectomycorrhizal growth response to apatite and ammonium in growth chambers with EMF-colonized Pinus muricata seedlings. It was found that apatite and ammonium addition increased the respiration rates of EMF, which could be taken as an indication of higher metabolic activity and probably higher mortality. Thus, it can be expected that EMF growing on the quartz bags has lower turnover than the mycelium colonizing the nutrient amendments, but this was not the case in this study. These discrepancies relating EMF turnover rates between the current and previous studies might be caused by shortcomings of the sequential incubation method used for the model in this paper. This method relies on the premise that the sum of the biomass from mesh bags incubated for short continuous periods should exceed the biomass from mesh bags incubated for a long incubation time. However, in a number of cases the mycelial biomass from a long incubation period was greater than the sum of the consecutive shorter intervals. This could be caused by a delay or a lag phase in EMF colonization inside the bags. It is possible that when a mesh bag was collected and the same hole was used to replace a new bag (Fig. 1) there was a lag phase before the hyphae could colonize the newly placed mesh bag due to disturbance of the mycelial connections (Wallander et al., 2013). Thus, those data points could have created noise in the data, making the turnover estimates less robust. In any case, if turnover in the EMF communities colonizing the nutrient amended bags is higher (as suggested by previous studies) and was underestimated in the current study, then the high standing biomass measured in the urea and apatite bags can only be explained by even higher EMF production than predicted in these results.

4.3 Seasonal effects on EMF growth

The general assumption of Model 1 is that fungal growth occurs at a constant rate. However, this approximation has some limitations, since seasonality usually affects the amount of C allocated to the roots (Endrulat et al., 2016) and consequently EMF root colonization (Taniguchi et al., 2018). Indeed, the standing EMF biomass in the mesh bags peaked in July and decreased over autumn, contradicting our fourth hypothesis (Fig. 2). In this paper Model 2 allowed the predicted fungal growth to vary both with seasonality and with the treatments (P fertilization and mesh-bag amendment). The introduction of these different dependencies in the model allowed us to test for the interactions between treatment and seasonal effects. It must be noted that the predicted EMF growth resulting from Model 1 is not incorrect and truly reflects the EMF growth differences between the treatments. However, by including seasonality in Model 2, we could determine that those differences predicted earlier were highly dependent on the season. Indeed, EMF growth not only increased early in the season, but the magnitude of this increase also depended on the treatments (Fig. 5). Therefore, the differences in biomass production between the fertilization regime and mesh-bag amendments were significant only early in the season (Fig. 6).

In contrast with our fourth hypothesis, the EMF biomass production peaked in summer and decreased in autumn. This contrasts with previous studies that have reported that the standing biomass in mesh bags collected from Pinus sylvestris (Hagenbo et al., 2021; Wallander et al., 2001), Pinus pinaster (Hagenbo et al., 2021) and Picea albies (Wallander et al., 2001) forests was higher during the autumn season. The unexpected growth seasonal patterns could have been caused by year-to-year variation in climatic conditions. However, according to climate data, the temperature and precipitation differences between summer and autumn in the year of sampling were not particularly different from other years. In a study performed in the same experimental area as the present study, Wallander et al. (2013) found that the standing biomass in September–October incubations was lower than the standing biomass in July–August incubations. It has been reported that different EMF species have different seasonal peaks (Castaño et al., 2017; Iotti et al., 2014; De la Varga et al., 2013), which could explain the differences in EMF growth between previous studies and the current experiment. Our results are also consistent with those from Coutts and Nicoll (1990), who found that the mycelium extension of Laccaria proxima and Telephora terrestris inoculated in Picea sitchensis peaked during July and decreased in autumn. The mycelial extension was associated with soil temperature, which peaked early in the growing season.

4.4 Potential non-mycorrhizal growth in the mesh bags

It could also be possible that non-mycorrhizal fungi contributed to the fungal growth detected in the current study. The main assumption that the ergosterol in this experiment comes mostly from EMF relies on previous evidence that the mesh-bag system favors the growth of EMF over non-mycorrhizal fungi (Almeida et al., 2019; Rosenstock et al., 2016; Berner et al., 2012; Wallander et al., 2010, 2001; Hedh et al., 2008). However, it has been shown that the shorter the time period a mesh bag remains underground the higher the proportion of non-mycorrhizal fungi inside the bags (as measured by the proportion of non-mycorrhizal DNA in Hagenbo et al., 2018).

Thus, non-mycorrhizal fungi growth could partially explain the seasonal effect detected as this fungal guild has been reported to respond positively to temperature (Pietikäinen et al., 2005). Unfortunately, the current study lacks non-mycorrhizal biomass controls (ie: fungal biomass from ingrowth bags collected in a trenched root-free area) that can be used to estimate the contribution of non-mycorrhizal fungi. Therefore, we cannot rule out the possibility that some of the ergosterol measured in the bags came from non-mycorrhizal fungi (i.e., methylotrophic yeasts in the urea-amended bags that could use methylene urea as both C and N sources). Even so, the significant negative effect of P fertilization on all the mesh-bag types suggests that the decrease in fungal growth might be related to a potential reduction in C allocation by the trees as discussed earlier. Moreover, the effects of the P fertilization and mesh-bag amendment on fungal growth were higher early in the season, which might imply that the seasonal effect seen in the current study is explained mostly by EMF.

It must be noted nevertheless that a potential reduction in belowground C allocation could decrease root activity and possibly root exudates, which might reduce labile sugars in the soils, affecting saprotrophic fungi as well. Further studies are necessary to evaluate the effect of P limitation on root dynamics and other members of soil microbial communities.

In conclusion, EMF production was strongly reduced when the P fertilizer was added to the forest, suggesting a decline in belowground C allocated by the trees to EMF when the P limitation was alleviated. This decline affected the colonization of the apatite and urea mesh bags, which might indicate that a potential decrease in belowground C allocation affected foraging for P but also foraging for N patches. The strong negative effect of P fertilization on EMF production suggests a central role of P in regulating EMF biomass production in N-rich forests. Moreover, the effect of the reduced belowground C allocation and the nutrient patches on EMF growth was significant only in the warmest months of the growing season, suggesting an important effect of seasonality on EMF growth dynamics and nutrient uptake.