Yield models for predicting aboveground ectomycorrhizal fungal productivity in Pinus sylvestris and Pinus pinaster stands of northern Spain

Abstract

Background

Predictive models shed light on aboveground fungal yield dynamics and can assist decision-making in forestry by integrating this valuable non-wood forest product into forest management planning. However, the currently existing models are based on rather local data and, thus, there is a lack of predictive tools to monitor mushroom yields on larger scales.

Results

This work presents the first empirical models for predicting the annual yields of ectomycorrhizal mushrooms and related ecosystem services in Pinus sylvestris and Pinus pinaster stands in northern Spain, using a long-term dataset suitable to account for the combined effect of meteorological conditions and stand structure. Models were fitted for the following groups of fungi separately: all ectomycorrhizal mushrooms, edible mushrooms and marketed mushrooms. Our results show the influence of the weather variables (mainly precipitation) on mushroom yields as well as the relevance of the basal area of the forest stand that follows a right-skewed unimodal curve with maximum predicted yields at stand basal areas of 30–40 m²∙ha^− 1.

Conclusion

These models are the first empirical models for predicting the annual yields of ectomycorrhizal mushrooms in Pinus sylvestris and Pinus pinaster stands in northern Spain, being of the highest resolution developed to date and enable predictions of mushrooms productivity by taking into account weather conditions and forests’ location, composition and structure.

Background

The importance of wild forest fungi as a key component of both ecosystem processes and services at scales ranging from local to global have been widely reported (Boa 2004). Among them, ectomycorrhizal fungi are especially relevant because of their ecological and socioeconomic importance (Smith and Read 2008; Bonet et al. 2014). Mycorrhizal species form symbiotic associations with their host plants that directly influence nutrient and water availability for trees (Brunner 2001; Hartnett and Wilson 2002; Guidot et al. 2003; Smith and Read 2008), facilitating seedling establishment and supplying and recycling soil nutrients (Egli 2011; van der Heijden et al. 2015), protecting plant hosts from soil pathogens and environmental extremes (Smith and Read 2008), and playing an important role in the sequestration of C in soil and trees (Treseder and Allen 2000; Egli 2011). From a socioeconomic point of view, the commercialization of the fruitbodies of ectomycorrhizal fungi provides important economic benefits to collectors and rural communities (Samils et al. 2008; Martínez de Aragón et al. 2011; Bonet et al. 2014). In fact, the value of edible epigeous fungi may surpass the value of timber in certain regions such as the Mediterranean (Palahí et al. 2009) or the northwest part of Turkey under some restrictive conditions based on forest management (Mumcu Küçüker and Başkent 2017b), and emerging activities such as mycotourism can further contribute to increasing this value (Büntgen et al. 2017).

Within this context, it is important to have accurate estimations of ectomycorrhizal fungal production at national level, not only for forest and land managers to integrate them into forest management planning or for the industrial sector to establish business strategies, but also to comply with international reporting requirements for Forest Europe (Forest Europe 2015) and FAO (Forest Resource Assessment-FRA) (MacDicken 2015). Empirical models can contribute directly to these tasks by providing both qualitative understanding and quantitative predictions of the impact of various management practices and climatic scenarios on forest ecosystem behavior over different spatio-temporal scales, allowing to integrate the production of different type of products in existing management planning systems such as stand simulators and Decision Support Systems (DSS) tools (Sánchez-González et al. 2015; Mumcu Küçüker and Başkent 2017a, 2017b).

In Spain, different edible mushroom yield models at local and regional levels have been developed so far, mainly for pine forest ecosystems. Bonet et al. (2008) developed a model for predicting ectomycorrhizal mushroom yield and species richness as a function of site and forest stand variables in Scots pine forests of north-eastern Spain. Later on, Bonet et al. (2010) developed similar models for pine forests (Pinus sylvestris L., Pinus halepensis Mill., Pinus nigra J.F.Arnold) in south-central Pyrenees. Martínez-Peña et al. (2012) developed models for predicting the productions of ectomycorrhizal mushroom in P. sylvestris forests located in north-central Spain with a special focus on the most valuable species. De-Miguel et al. (2014) developed a model-based scenario analysis for predicting the effect of forest management intensity on mushrooms productivity in pine forests (P. sylvestris, P. halepensis, P. pinaster and P. nigra) in Catalonia region (north-eastern Spain). More recently, Taye et al. (2016) predicted edible mushroom yield in Pinus pinaster forests of central Spain under different meteorological and site conditions. However, there is a lack of models enabling accurate enough prediction of mushroom yield on a larger scale. The reason behind this lack of national-level models is the difficulty to obtain large quantities of data over several years, especially if the models are intended to account for the effect of changing meteorological conditions (Taye et al. 2016; Alday et al. 2017a). In Spain, these kind of data are only available in the northern part of the country, which represents one of the largest existing spatio-temporal data series on aboveground fungal yield worldwide. All the aforementioned mushrooms yield models already used part of the data used in the present study, but none have used the whole dataset.

In this study, we have used data from P. sylvestris and P. pinaster stands because these pine species are monitored for mushroom yield from permanent plots established in different locations of northern Spain, which is a consequence of the importance of both pine species in terms of their wide distribution range and mushroom production in the Iberian Peninsula. Previous works have shown that climatic factors are key for controlling the yields of mushrooms (Kauserud et al. 2008; Ágreda et al. 2016; Alday et al. 2017b; Karavani et al. 2018). Wet and warm autumns seem to promote ectomycorrhizal fungi yields in Spanish pine forests (Martínez-Peña et al. 2012; Taye et al. 2016; Alday et al. 2017a;). However, climatic factors alone do not fully explain the emergence of fungal sporocarps: topographical factors such as slope, aspect or altitude (Egli 2011), soil characteristics (i.e., pH, texture) (Martínez-Peña et al. 2012) and stand structure variables (Tahvanainen et al. 2016) seems to be also influential in mushroom fructifications (Tomao et al. 2017). Among the latter group of variables, stand basal area has been reported to be correlated with mushroom production (Bonet et al. 2008; Bonet et al. 2010; Martínez-Peña et al. 2012; de-Miguel et al. 2014; Tahvanainen 2014). Those studies shown that forest stands with too low or too high basal areas can be less productive in terms of mushroom yield than stands with intermediate values of stand basal area.

The main aim of this study is to develop yield models for predicting the annual productivity of ectomycorrhizal mushroom in P. sylvestris and P. pinaster stands in northern Spain, taking into account temperature, precipitation, stand basal area and local site characteristics as explanatory variables. Models were fitted for the following groups of fungi separately: all ectomycorrhizal mushrooms, edible mushrooms and marketed mushrooms, based on data from 90 sample plots from different locations in northern Spain. We hypothesized that mushroom production will follow certain biogeographic patterns and, therefore, the variability among plots in relation to mushroom production, climate and site conditions was also studied.

Methods

Sampling design and data collection

The study area is located in three different regions of Northern Spain. In total, seven monitoring sites amounting 90 sampling plots (100 m² each plot) have been considered in this study. Of those plots, 39 plots correspond to pure P. sylvestris forest stands whilst 51 plots represent pure P. pinaster forest stands.

The P. sylvestris plots are located in Catalonia region (North-eastern Spain, 19 plots measured between 1997 and 2015), in Soria Province, Eastern Castilla y León region (North-Central Spain, 17 plots measured between 1995 and 2015) and in Palencia Province, Western Castilla y León region (North Western Spain, 3 plots measured between 2008 and 2015). The P. pinaster plots are located in Catalonia region (North-eastern Spain, 28 plots measured between 2008 and 2015), in Soria Province, Eastern Castilla y León region (North-Central Spain, 14 plots measured between 1997 and 2015) in Palencia Province, Western Castilla y León region (North Western Spain, 3 plots measured between 2003 and 2015) and in Valladolid Province, Western Castilla y León region (North Western Spain, 6 plots measured between 2006 and 2015).

All the considered plots have been weekly monitored during the autumn season, with data recorded for at least 9 consecutive years or even more. The permanent mushroom plots are representative of the heterogeneity of the forest areas of both pine species in northern Spain, being located in a wide range of ecological conditions and subjected to different forest management treatments (mainly forest thinning). The data were characterized by high interannual variability in mushroom occurrence and yield associated with differences in weather conditions between years and between plots, and representing different past forest management practices Additional file 1: Table S1). Weather data was obtained from the nearest weather station to each study area. More detailed descriptions of the study areas as well as sampling methodology can be found in Vásquez-Gassibe et al. (2016) and Hernández-Rodríguez et al. (2015) for Western-Castilla y León (Palencia and Valladolid provinces), Martínez-Peña et al. (2012) and Taye et al. (2016) for Eastern-Castilla y León, and de-Miguel et al. (2014) and Alday et al. (2017b) for Catalonia.

Preliminary analysis

An exploratory principal component analysis (PCA, Legendre and Legendre 1998) was performed to explore the variability among plots with respect to mushroom production and climate variables. The PCA were run with PRINCOMP procedure available in SAS version 9.4. (SAS Institute Inc. 2016).

Next, we clustered the plots based on mushroom production, climatic and site conditions. Hierarchical clustering is a method of forming clusters iteratively, starting with each object in its own cluster and then proceeding by combining the most similar pairs of clusters step by step, thus forming a hierarchy of clusters (e.g. Everitt et al. 2011). We performed hierarchical clustering on the mean values of the yields over a study period, using Euclidean distance as a measure of similarity and Ward’s minimum variance method as the clustering method. The number of clusters was selected based on the dendrogram and the cubic clustering criterion (CCC) (Milligan and Cooper 1983; Yeo and Truxillo 2005). The cluster analysis was run with CLUSTER procedure available in SAS version 9.4. (SAS Institute Inc. 2016).

Modeling mushroom production

The annual mushroom yields were modelled for the following groups of fungi separately: all ectomycorrhizal mushrooms, edible mushrooms (those ectomycorrhizal fungi considered as edible in the available fungal literature) and marketed mushrooms (ectomycorrhizal edible fungi usually sold in markets) (de-Miguel et al. 2014; Alday et al. 2017a), as a function of location and variables representing different meteorological conditions (i.e., monthly total rainfall and mean temperature), stand characteristics (i.e., stand basal area) and thinning treatment. Since, in Spain, wild edible mushrooms are usually commercialized on a fresh weight basis, fresh mushroom biomass in kg∙ha^− 1∙yr^− 1 was selected as the response variable by pooling the yield data of all fungal species according to the three levels of grouping. Different combinations and transformations of predictors were tested. Weather variables were aggregated in different ways, e.g., the accumulated precipitation during August and/or September (late summer) or during September and/or October (early autumn), to further test their combined effect on the response variables in addition to testing the influence of the disaggregated monthly rainfall and temperature. The different combinations of meteorological variables also aimed at testing hypothetical delayed responses of mushroom yield to the combined effect of different predictors (e.g., previous research has reported a delay of several weeks in the combined effect between rainfall events and favorable temperatures) (Martínez de Aragón et al. 2007; Martínez-Peña et al. 2012).

Since the available data are based on repeated measurements of the same sampling plots during several years, measurements taken on a given plot are likely to be more correlated than measurements taken from different plots. Similarly, measurements taken closer in time on the same plot (i.e., in a given year) are likely to be more correlated than measurements taken further apart in time. Such autocorrelation patterns implies that assumptions about error variance being independent are no longer valid (Wolfinger 1996; Littell et al. 2000). The analysis of repeated measurements requires that correlations between the observations made on the same sampling unit must be taken into account as well as possible heterogeneous variances among observations on the same plot over time. Second, data are unbalanced because the number of sample plots varied among groups and each plot there were no available data for the same years. Third, mushroom yield has a stochastic nature that coupled with the rather small size of sample plots results in the occurrence of “zero” production in many sample plots.

To deal with these characteristics of the data, hurdle models within a generalized linear mixed-effects modeling framework (GLMM) were used. Hurdle models model the zeros and non-zeros as two separate processes (Hamilton and Brickell 1983) which, as compared with single-model functions, can also provide further insight into mushroom dynamics by analyzing those factors driving mushroom occurrence and abundance separately (de-Miguel et al. 2014). Therefore, we applied this approach for the three considered groups of fungi although the percentage of zeros differs in each one, being 4.25% in the all ectomycorrhizal mushroom group, while in the edible and marketed groups the proportion of observations with zero mushroom production was 8.86% and 31%, respectively. The first part of the hurdle models aimed at predicting the probability of occurrence of mushroom production based on binomially distributed data (i.e., absence or presence) using logistic regression (Eq. 1) along with a logit link function (Eq. 2). The second part of the hurdle models aimed at predicting mushroom yield conditional on the probability of mushroom occurrence by means of Gamma regression (Eq. 3) along with a log link function (Eq. 4). Finally, the expected mushroom yield was obtained by multiplying the estimates provided by Eq. 1 and Eqs. 3 and 5.

$$ p\left({y}_{ij}=1|x\right)=\uppi (x)=\frac{1}{1+{e}^{-\left[\left({\alpha}_0+{v}_{1j}\right)+\alpha {\boldsymbol{X}}_{ij1}\right]}} $$

(1)

$$ g(x)=\log \left[\frac{\pi (x)}{1-\pi (x)}\right]=\left({\alpha}_0+{v}_{1j}\right)+\alpha {\boldsymbol{X}}_{ij1}\kern0.5em $$

(2)

$$ {\mathrm{yield}}_{cij}={\mathrm{e}}^{\beta_0+{v}_{2j}}{\boldsymbol{X}}_{ij2} $$

(3)

$$ g(x)=\log \left({\mathrm{e}}^{\beta_0+{v}_{2j}}{\boldsymbol{X}}_{ij2}\right)=\left({\beta}_0+{v}_{2j}\right)+\beta \log \left({\boldsymbol{X}}_{ij2}\right) $$

(4)

$$ {\mathrm{yield}}_{ij}=p\left({y}_{ij}=1|x\right)\bullet {\mathrm{yield}}_{cij} $$

(5)

where p(y_ij = 1| x) is the probability of occurrence of all ectomycorrhizal, edible and marketed mushrooms in plot i and year j, yield_cij is all ectomycorrhizal, edible and marketed mushrooms yield conditional on mushroom occurrence in plot i and year j (kg∙ha^− 1∙yr^− 1), yield_ij is the predicted all ectomycorrhizal, edible and marketed mushrooms yield in plot i and year j (kg∙ha^− 1∙yr^− 1), α₀ and β₀ denote fixed-effects, v_1j and v_2j denote year random effects which were specified as crossed effects, and X_ij1 and X_ij2 denote vectors of predictor variables in plot i and year j.

The resulting groups of plots from the cluster analysis were included as fixed dummy variables within the models. In addition, to test whether differences among groups and years were statistically significant a repeated measures ANOVA was performed. The differences were examined using pairwise comparisons according to the Tukey method using MIXED procedure available in SAS 9.4 (SAS Institute Inc. 2016).

Several site and climatic variables, as well as their transformations were included as potential predictors in the model. Models were fitted adding 2 year random effects in the intercept of each part of the hurdle model, v₁ and v₂. These effects are distributed under a normal distribution with mean zero ${\upsigma}_{\mathrm{b}}^2,{\upsigma}_{\mathrm{s}}^2$. The unstructured covariance structure was used to describe the variance-covariance structure of the random effects (Littell et al. 2000). Plot random effect were not included in the models because the variance of the plot random effects were practically zero. All the models were fitted using the NLMIXED procedure available in SAS version 9.4 (SAS Institute Inc. 2016).

Model selection and evaluation

The models were selected and evaluated according to the following criteria: biological sense, goodness-of-fit and predictive ability. The biological sense was evaluated considering whether alternative models behaved logically, i.e., whether they represented biologically or ecologically consistent relationships between predictors and the response variables according to current scientific and expert knowledge. Only those models whose coefficients were statistically significant (p < 0.05) were further considered in the analysis.

The goodness-of-fit of the models was analyzed according to the mean bias, which reflects the deviation of model predictions against observed values, and the root mean square error (RMSE), which account for the precision of the estimates. The relative values of these statistics (BIAS%, RMSE%) were calculated by dividing the mean bias and RMSE by the mean of the predicted mushroom yield. The Akaike’s information criterion (AIC) and likelihood ratio tests were also used to further guide the selection of predictor variables to prevent overfitting by accounting for the trade-offs between model parsimony and goodness of fit. Furthermore, uncertainty was assessed also using resampling techniques, namely bootstrapping based on 2000 bootstrap samples with replacement, to ensure the stabilization of the estimates, and by computing prediction and confidence intervals accounting for the residual variance, the uncertainty in the fixed coefficients, and the uncertainty in the variance parameters of the year random effects. In addition, receiver operating characteristic (ROC) curves and the corresponding area under the ROC curve (AUC) along with their bootstrapped confidence intervals were computed for the logistic models for the probability of mushroom occurrence.

Results

Exploratory analysis of data

The principal component analysis (PCA) applied for exploring the variability among plots with respect to mushroom yield, climate and site variables showed that the differences among groups of plots are related to higher precipitation levels in summer (Fig. 1). The first two principal components (PC1 and PC2) accounting for 70.78% of the data variance (PC1 46.83%, PC2 23.95%) were assumed as the most principal components. The sum of all remaining (5) PCs accounted for the 29.22% of the data variance. As an expected result, PCA results indicated that high water availability in late summer favored good mushrooms yields while temperature is a less determining factor for mushrooms production in the study area. Besides, the score plots (Fig. 1) suggests some groups of plots which were verified by a cluster analysis.

Fig. 1

Scatterplots of the loadings (above) and the scores (below) obtained by PCA applied on plots. In parentheses are shown the percentage of variability accounted by the two first components. “total” denotes all ectomycorrhizal mushrooms yield, “edible” denotes edible mushrooms yield, “marketed” denotes marketed mushrooms yield, “altitude” denotes the plot altitude, “latitude” denotes the plot latitude, “longitud” is denotes the plot longitude, “slope” denotes the main slope of the plot, “aspect” denotes the plot aspect, “slopeasp” is a synthetic variable created by multiplying slope by aspect, “P_annual” is the annual accumulated precipitation, “P_ag” is the accumulated precipitation of August, “P_set” is the accumulated precipitation of September, “P_oct” is the accumulated precipitation of October, “P_nov” is the accumulated precipitation of November, “P_as” is the accumulated precipitation of August and September, “P_aso” is the accumulated precipitation of August, September and October, “T_annual” is the annual mean temperature, “T_ag” is the mean temperature of August, “Tm_set” is the mean temperature of September, “T_oct” is the mean temperature of October, “T_nov” is the mean temperature of November, 21 is Pinus sylvestris and 26 is Pinus pinaster, Cat is Catalonia region, Pal is Palencia province, Sor is Soria province and Val is Valladolid province

The Ward’s hierarchical clustering method was performed for the mushroom yield, climate and site variables. The number of clusters was determined based on the cubic clustering criterion (CCC) method. The value obtained for five clusters was 4.53, which is higher than 2 indicating good clusters (Milligan and Cooper 1983). The plots formed the following five different clusters according to mushroom yield, climate and site characteristics (Fig. 2):

• 1) CAT1: made up of the 9 plots of Pinus sylvestris in Catalonia with the highest precipitations in summer and located at higher altitude.

• 2) CAT2: made up of 5 plots of P. sylvestris in Catalonia with higher precipitation levels in summer but located at lower altitude that the previous group.

• 3) CAT3: made up of all the plots of P. pinaster (28) in Catalonia and 5 plots of P. sylvestris in Catalonia with higher temperatures than the other plots of the same species in that region.

• 4) CyL1: made up of all the plots of P. sylvestris (20) in Castilla y Leon region (Palencia and Soria provinces) and three plots of P. pinaster (3) located in Palencia Province.

• 5) CyL2: made up of all the plots of P. pinaster (20) in Castilla y Leon region (Soria and Valladolid provinces).

When these five groups of plots were considered in the modelling as a fixed dummy variables, likelihood ratio tests indicated a significant improvement in the fit of the mushroom yield models.

Fig. 2

Clustering of 90 plots based on mushrooms yield, weather and site variables. Different colors indicate the five different clusters identified as follows: dark orangeCAT1, blue CAT2, green CAT3, light orange CyL1, and purple CyL2

A summary of the weather variables in each group of plots is shown in Additional file 1: Table S1. Additional file 2: Figure S1 and Additional file 3: Figure S2 of show mushroom yields by each group of plots considered and years. The average annual yield of all ectomycorrhizal mushrooms was 117.83 ± 4.53 kg·ha^− 1·yr^− 1. The highest yields amounting 500 kg·ha^− 1 were obtained in the plots of the group called CAT2 in 2014. The mean yield of all ectomycorrhizal mushrooms in 2006 was 330.76 ± 36.71 kg·ha^− 1 and in 2014 was 271.93 ± 21.62 kg·ha^− 1.

Regarding edible mushrooms, the highest yields, above 396.75 kg·ha^− 1, were obtained also in plots belonging to group CAT2 in 2014. Accordingly, this group of plots also produced the highest yields of marketed mushrooms in that year. In edible and marketed mushrooms, the average annual yield was 66.83 ± 2.87 and 36.18 ± 2.15 kg·ha^− 1·yr^− 1 respectively, while in 2014 the production was 171.84 ± 16.21 and 78.98 ± 12.78 kg·ha^− 1. In 1997 the production of edible mushrooms was 121.17 ± 21.44 kg·ha^− 1 and of marketed mushroom was 78.28 ± 16.90 kg·ha^− 1.

Models for the probability of mushrooms occurrence

The information about the selected models for predicting the probability of occurrence of all ectomycorrhizal, edible and marketed mushrooms is presented in Table 1. The logistic regression analysis showed that the mean temperature in November had a significant effect on the probability of occurrence of all ectomycorrhizal mushroom. The inclusion of this variable in the probability of occurrence of all ectomycorrhizal model decreased the random between-year variation. The probability of occurrence of mushrooms was different among the five groups of plots considered.

Table 1 Parameters estimations (Est) for the hurdle models for the three groups of fungi. α_i are the parameters of the part of the hurdle model that predicts the probability of occurrence of mushroom production, β_i are the parameters of the part of the hurdle model that predicts yield conditional on the probability of mushroom occurrence, v_i are variance of year random effect in each part of the hurdle model, Cat1 is a dummy variable set to one for CAT1 group of plots and zero for the rest of the groups, Cat2 is a dummy variable set to one for CAT2 group of plots and zero for the rest of the groups, Cat3 is a dummy variable set to one for CAT3 group of plots and zero for the rest of the groups, CyL1 is a dummy variable set to one for CyL1 group of plots and zero for the rest of the groups, CyL2 is a dummy variable set to one for CyL2 group of plots and zero for the rest of the groups, G is stand basal area (m²·ha^− 1), Pag is the accumulated precipitation of August, Pset is the accumulated precipitation of September, Poct is the accumulated precipitation of October, Pnov is the accumulated precipitation of November, Tnov is the mean temperature of November, the term “ln” refers to the natural logarithm. LowB and UpB are the lower and upper bounds of the bootstrapped 95% confidence intervals

The bootstrapped values of the area under the ROC curves (Fig. 3) indicated an excellent performance of the logistic models for the probability of all ectomycorrhizal mushroom and edible mushroom, with AUC values ranging between 89%–96%, and 83%–89%, respectively. For the marketed mushroom, the performance of the logistic model can be considered acceptable with AUC values ranging between 76%–82%. Although we recognized that these values may be overestimate due to having used the same data for model fitting and for calculating the AUC values (Copas and Corbett 2002), we think that the three models can be considered good enough for detecting the occurrence of the three types of fungi.

Fig. 3

The performance of the probability of occurrence models are depicted by the receiver operating characteristic (ROC) curves and the area under the ROC curve (AUC). In parentheses are shown the AUC estimates bootstrapped 95% confidence intervals

Models for mushrooms yield conditional on the probability of mushrooms occurrence

The information about the selected models for predicting all ectomycorrhizal, edible and marketed mushrooms yield conditional on the probability of mushrooms occurrence is presented in Tables 1 and 2. The yield of ectomycorrhizal, edible and marketed mushrooms was significantly different among all the considered groups of plots. All ectomycorrhizal and edible mushrooms yield conditional on the probability of mushroom occurrence were also influenced by stand basal area, the rainfall from August to November, and the mean temperature of November. The effect of basal area on the amount of mushrooms production was the same in the three groups of fungi considered, following a right-skewed unimodal curve (Fig. 4). Increasing values of the selected meteorological predictors increases ectomycorrhizal mushrooms yields, except for the rainfall in November indicating that high values of precipitation in this month might cause a decrease in the mushroom yields when also the previous months had high values of precipitation. Similar effects of these variables were also observed for the production of edible mushrooms. Regarding the marketed mushrooms yield, the yield of this type of fungi was only influenced by the total rainfall of September, indicating that higher precipitations in that period increases the yield of marketed mushrooms.

Table 2 Estimates of the fitting statistics of the hurdle models for the three groups of fungi. The estimates of the fitting statistics of the yield models not conditional on the probability of mushroom occurrence are shown between brackets. BIAS is the average error, RMSE is the root mean square error, AIC denotes the Akaike’s information criterion, LowB and UpB are the lower and upper bounds of the bootstrapped 95% confidence intervals

Fig. 4

Effect of basal area on mushroom yield. The values assigned to the predictors correspond the mean values of each groups of plots

The performance of the selected models for mushroom yield conditional on the probability of mushroom occurrence for each considered group of fungi is shown in Fig. 5. The best performance was shown by the model for marketed mushrooms, in which the expected relationship between observed and predicted values almost overlapped the perfect equality line. A good performance was shown as well by the model for edible mushrooms whose predictions also practically overlapped the equality line through the range of observed values. The performance of the model for all ectomycorrhizal mushrooms was also acceptable along the range where most observations fall, although tended to overestimate very low yield and underestimate maximum yields.

Fig. 5

Performance of the selected models for mushroom yield conditional on the probability of mushroom occurrence (Eqs. (3) and (4)). The solid line represents perfect fit (1:1 equality line) whereas the dashed line indicates the regression line between the measurements and the back-transformed conditional predictions (including random year effects) of the selected model along with their confidence intervals

Hurdle models predictions

The estimated variance of year random effect in both parts of the hurdle model was lower in the marketed mushrooms models than in the other two groups of mushrooms (Table 1), where the between-year variation was very similar. In the three selected models that between-year variation were higher in the first part of the hurdle model that predict the probability of occurrence of mushrooms production. This is explained by the fact that none meteorological variable was finally included in that part of the model. Although the inclusion of different precipitation- and temperature- related variables were tested, no one resulted statistical significant or the parameter sign was not biologically meaningful.

Table 2 shows fitting statistics along with their bootstrapped confidence intervals of the hurdle models for predicting the yield of the three considered groups of fungi. In terms of the deviation of the models with respect to observed values (bias), the selected models for each group of fungi were almost unbiased. While in terms of the precision of the estimates (RMSE), the best performance was shown by the model for all ectomycorrhizal mushrooms yield, followed by the edible mushrooms yield model.

For illustrating how the hurdle models behave related to mushrooms yield predictions, in Table 2 the values of the fitting statistics of the selected yield models not conditional on the probability of mushroom occurrence are shown between brackets. As expected, the greatest improvement achieved by the hurdle models appeared when the percentage of observations with zero mushroom production was higher, which was the case of the marketed mushrooms group (31%). But even when the percentage of zeros was rather low, such as for the edible mushrooms group (9%) and all ectomycorrhizal mushrooms group (4%), the goodness-of-fit was better or very similar for the hurdle models. Indeed, the hurdle models tended to enhance model performance especially in terms of mean bias reduction.

Discussion

Mushrooms production variability

Mushroom yields are characterized by a high temporal (i.e., interannual) and spatial (i.e., between-location) variability (Alday et al. 2017a). The interannual variability is mainly due to changes in the meteorological conditions between years, while spatial variability can be mainly due to differences in soil properties, forest management, stand composition and structure, and site characteristics (Martínez-Peña et al. 2012; Primicia et al. 2016; Taye et al. 2016; Collado et al. 2018). The mushroom yield models developed in this study account for both sources of variability which will provide forest managers with further knowledge about the patterns of mushroom production in the north of the Iberian Peninsula.

The production of ectomycorrhizal mushrooms were exceptional high in the plots of the group called CAT2 in year 2014. Alday et al. (2017a) also identified that same year as exceptionally productive when studying mushrooms production from 1997 to 2014 in different forest stands in Spain.

Effect of weather on mushrooms occurrence and productivity

Our results show a stronger influence of the weather conditions on mushroom yields than on the appearance of fungal fruit bodies (Table 1). Although previous studies have reported that precipitation and temperature are important drivers for both mushroom occurrence and yield (Martínez-Peña et al. 2012; Hernández-Rodríguez et al. 2015; Mumcu Küçüker and Başkent 2015; Tahvanainen et al. 2016; Karavani et al. 2018), in our models for predicting the probability of occurrence of all ectomycorrhizal mushrooms only the mean temperature in November appeared as a significant predictor and no precipitation variable appeared as a significant predictor in the probability of occurrence of any group of fungi considered. This result can be explained by the introduction of dummy variables related to the group of plots and the year random effects considered in model fitting. Both of them contributed to explaining part of the between-year variation arising from annual changes on the meteorological conditions of each considered regions, and this effect is more patent in occurrence models than when predicting the yield.

As already reported by Karavani et al. (2018), total and edible mushrooms were related to the same set of climate predictors, which is an expected result because most of the annual ectomycorrhizal mushroom yields can be considered as edible. The cumulated precipitation in late summer and early autumn, from August to October, was an important predictor of total and edible mushroom productivity, while for marketed mushrooms the precipitation of September was the main meteorological predictor. These results agree with previous studies (Tahvanainen et al. 2016; Alday et al. 2017a) indicating that high water availability before or at the beginning of the fruiting season would favor good mushrooms yields. However, to maintain high total and edible mushrooms yields during the whole fruiting season it seems to be necessary lower rainfall and higher temperature in November. The negative effect of precipitation on mushroom yields in excessive wet conditions was already reported by Boddy et al. (2014) and Karavani et al. (2018). These last authors also found that cold temperatures at the end of the season can also have a negative influence on mushroom yields. In contrast, marketed mushroom yields was not related to precipitation from October to November or to temperature in November. This result suggest an earlier phenology of the marketed species (Karavani et al. 2018).

Effect of stand structure on mushrooms occurrence and productivity

The composition of the mycorrhizal fungal community is strongly influenced by forest stand structure, being stand basal area one of the most used variables to describe the stocking of forests stands (Tomao et al. 2017). In our models, stand basal area had a significant effect on the yield of all the groups of fungi considered. The optimal stand basal area at which mushroom production was maximized was about 40 m²·ha^− 1 for total and edible mushrooms and 30 m²·ha^− 1 for marketed fungi, which is in accordance with previous research. In Spain, Martínez-Peña et al. (2012) reported that values around 40 m²·ha^− 1 maximized the Boletus edulis yields in the P. sylvestris forests of Central Spain; and de-Miguel et al. (2014) reported values between 30 and 40 m²·ha^− 1 as the optimum basal area for edible and marketed mushroom production in pure P. pinaster stands growing in good conditions. In Finland, Tahvanainen et al. (2016) found that the most suitable basal area for Lactarius sp. and all marketed mushrooms was 30 m²·ha^− 1.

Basal area’s effect follows a right-skewed unimodal curve, although the decreasing trend is not very pronounced (Fig. 4), indicating that high mushrooms yields can be found from 20 to 30 m²·ha^− 1 over a wide range of basal areas. This pattern seems to indicate that mushrooms yields may be maximized when the forest stand is not too dense or too sparse. The effect of thinning was found not significant in the modeling process. This result can be explained by the fact that thinning was applied only once in 31 out of the 90 plots included in the dataset, and also because the groups of fungi considered included many species. A review of the literature shows different and sometimes contradictory results about the effect of thinning on mycorrhizal ecosystems, because such forestry operation may modify fungal succession patterns aboveground and provide favorable conditions for certain species in detriment to others (Bonet et al. 2012; Mumcu Küçüker and Başkent 2017b; Tomao et al. 2017), although recent research did not find any thinning effect on the fungal community belowground (Castaño et al. 2018).

Conclusions

The mushroom yield models developed in this study are the first empirical models for predicting the annual yields of ectomycorrhizal mushrooms in Pinus sylvestris and Pinus pinaster stands in northern Spain. These models are of the highest resolution developed to date and enable predictions of mushrooms productivity by taking into account weather conditions and forests’ location, composition and structure. The modeling approach used, i.e., hurdle mixed models, allowed a better insight into the processes driving fungi emergence and abundance, and therefore into the ecological system. In addition, these models can be very useful for forest and land managers, since the integration of them into forest management planning can facilitate the decision making process in multifunctional forests because provide further knowledge about the patterns of mushroom production and associated ecosystem services in the north of the Iberian Peninsula. Since the basal area has been demonstrated as a significant factor influencing the mushroom yields, this open the range of options for forest managers, who are usually managing the forests only considering timber as an economic option. The consideration of both timber and mushrooms in the definition of the optimal management schedule through optimization techniques may increase the profitability of the northern Iberian Peninsula forests as demonstrated by Palahí et al. (2009).