Predicting mushroom productivity from long-term field-data series in Mediterranean Pinus pinaster Ait. Forests in the context of climate change

Article Predicting Mushroom Productivity from Long-Term Field-Data Series in Mediterranean Pinus pinaster Ait. Forests in the Context of Climate Change Celia Herrero 1,2,*, Iosu Berraondo 3, Felipe Bravo 2,4, Valentín Pando 2,5, Cristóbal Ordóñez 2,4 , Jaime Olaizola 1, Pablo Martín-Pinto 2,4 and Juan Andrés Oria de Rueda 2,3 1ECM Ingeniería Ambiental. R + D + I department, C/Curtidores 17, 34003 Palencia, Spain; [email protected] 2Sustainable Forest Management Institute, University of Valladolid-INIA.ETS Ingenierías Agrarias, University of Valladolid, Avda, Madrid 44, 34071 Palencia, Spain; [email protected] (F.B.); [email protected] (V.P.); [email protected] (C.O.); [email protected] (P.M.-P.); [email protected] (J.A.O.d.R.) 3Departamento de Ciencias Agroforestales, E.T.S. Ingenierías Agrarias, Avda, Madrid 57, 34004 Palencia, Spain; [email protected] 4Departamento de Producción Vegetal y Recursos Forestales, E.T.S. Ingenierías Agrarias, Avda, Madrid 44, 34004 Palencia, Spain 5Departamento de Estadística e Investigación Operativa, E.T.S. Ingenierías Agrarias, Avda, Madrid 57, 34004 Palencia, Spain *Correspondence: [email protected]; Tel.: +34-979-741-006 Received: 31 January 2019; Accepted: 21 February 2019; Published: 26 February 2019   Abstract: Long-term field-data series were used to fit a mushroom productivity model. Simulations enabled us to predict the consequences of management and climate scenarios on potential mushroom productivity. Mushrooms play an important ecological and economic role in forest ecosystems. Human interest in collecting mushrooms for self-consumption is also increasing, giving forests added value for providing recreational services. Pinus pinaster Ait. is a western Mediterranean species of great economic and ecological value. Over 7.5% of the total European distribution of the species is found on the Castilian Plateau in central Spain, where a great variety of mushrooms can be harvested. The aim of this study was to model and simulate mushroom productivity in Maritime pine (Pinus pinaster Ait.) ecosystems in northern Spain under different silvicultural and climatic scenarios. A mixed model was fitted that related total mushroom productivity to stand and weather variables. The model was uploaded to the SiManFor platform to study the effect of different silvicultural and climatic scenarios on mushroom productivity. The selected independent variables in the model were the ratio between stand basal area and density as a stand management indicator, along with precipitation and average temperatures for September and November. The simulation results also showed that silviculture had a positive impact on mushroom productivity, which was higher in scenarios with moderate and high thinning intensities. The impact was highly positive in wetter scenarios, though only slightly positive and negative responses were observed in hotter and drier scenarios, respectively. Silviculture had a positive impact on mushroom productivity, especially in wetter scenarios. Precipitation had greater influence than temperature on total mushroom productivity in Maritime pine stands. The results of this paper will enable forest managers to develop optimal management approaches for P. pinaster forests that integrate Non-Wood Forest Products resources. Keywords: Maritime pine; modelling; mushrooms; non-wood forest products; SiManFor; simulation Forests 2019,10, 206; doi:10.3390/f10030206 www.mdpi.com/journal/forests Forests 2019,10, 206 2 of 18 1. Introduction Fungi play an essential ecological role in forest ecosystems. Mycorrhizal species influence the nutrient and water uptake, absorption [ 1 ], growth and survival [ 2 , 3 ] of symbiotic plants, the presence of most vascular plants [ 1 ], soil Carbon storage, structure, aeration and porosity [ 4 , 5 ] and plant resistance to pathogens, especially at the root level [ 6 ]. Saprotrophic fungi are essential decomposers of dead matter, and therefore crucial to nutrient cycling in forest ecosystems [ 7 ]. Among the most important non-wood forest products, forest mushrooms also provide recreational services and are marketable (potentially profitable) as edible, nutraceutical and medicinal products or components. In some areas, the value of mushroom resources greatly exceeds that of timber [ 8 , 9 ] and companies have been created to process edible and medicinal fungi. In the Mediterranean area, holistic management of forests is a great challenge because of the different ecosystem services provided. In light of expected climate change, it is necessary to develop sustainable forest management techniques and new silvicultural strategies to improve the resilience of forests and thus enable ongoing provision of goods and services, including mushroom resources. Silviculture can modify density, canopy cover, primary productivity, basal area, understory plant communities, soil conditions and soil microbial communities. All these micro-environmental conditions can strongly affect fungal taxa fructification and the subsequent overall yield and diversity of the forest system [10]. Because mushroom emergence and stochasticity vary significantly, models must be based on long-term mushroom inventory data to provide reliable estimates, especially when attempting to understand the effects of climatic correlations. Long historical dataset are difficult to obtain, which explains the limited number of mushroom yield models. However, the number of studies focused on researching the influence of different factors on mushroom productivity continues to increase. Studies of mushroom dynamics show the influence of site characteristics (altitude [ 11 ], slope, aspect), stand variables (tree species, stand density, stand age [ 12 ], basal area [ 13 ], dominant height, tree growth [ 14 , 15 ]), soil properties [ 16 , 17 ] and weather conditions (precipitation, temperature [ 18 , 19 ]). Different approaches have been used in mushroom modeling, including correlations [ 20 ], stepwise multiple regressions [ 21 ], classificatory simple model [ 22 ], linear and nonlinear models [ 23 , 24 ], mixed models [ 13 , 14 , 17 , 18 , 25 , 26 ] or the two-step model approach [ 27 , 28 ]. Most of the models incorporate the stochastic between-plot and between-year structure to account for nested structures in the data [ 29 ] and thereby avoid biased standard error estimators of the parameters. Incorporation of plot and year factors in the models has improved estimation efficiency by bringing more information to the model, independently of the other measurements [29]. Simulation tools are essential for sustainable forest management, since they allow foresters to explore different management and silvicultural alternatives. In this context, the SiManFor [ 30 ] online platform (www.simanfor.es) has been developed to facilitate simulation of sustainable forest management alternatives at the stand level. SiManFor provides utilities for assessing forest growth models and simulates silvicultural scenarios inherent to forest management. This tool could also be linked to climate scenarios, which are essential for understanding the potential effects of climate change in general, and mushroom productivity in particular [31,32]. Pinus pinaster Ait. is a common tree species, distributed throughout coastal regions in the Mediterranean basin of Europe and Africa, along the Atlantic coasts of Portugal, Spain and France, and in the mid-altitude mountain ranges and inland plains of the Iberian Peninsula. This conifer is widely used in artificial reforestation efforts to avoid soil loss and desertification of large areas and to recover the original forests. It is a species of great economic, ecological and aesthetic value, particularly on the Castilian Plateau in central Spain. There, P. pinaster forests cover more than 114,000 ha, which represents about 7.5% of the total European distribution of this species. Pinus pinaster forests are important Mediterranean ecosystems because of the fungi associated with them. Several edible and marketable species such as Hygrophorus latitabundus Britz, Lactarius deliciosus (L.) S.F. Gray, Macrolepiota excoriata (Schaeff.) M.M. Moser, Macrolepiota konradii (Huijsm. Ex Forests 2019,10, 206 3 of 18 Orton) Moser, Macrolepiota mastoidea (Fr.) Singer, Macrolepiota procera (Scop.) Sing, Suillus luteus (L.) Roussel, Tricholoma portentosum (Fr.) Quél or Tricholoma terreum (Sch.) Kumm can be found in these forests [33]. In this study, we looked at how stand, site and weather variables might be correlated with mushroom productivity. We hypothesized that stand, site and weather variables would influence mushroom productivity and that stand variables susceptible to forest management decisions could positively affect mushroom fruiting. In a more advanced approach linked to the predictive model, we also wanted to identify patterns associated with different silvicultural and climatic conditions. Since the Mediterranean climate pattern is characterized by high temperatures and low rainfall during summer, with less precipitation expected at higher temperatures, we hypothesized that climate change conditions would affect mushroom productivity and that predictive models and scenario simulations could be useful in the preservation and production of this interesting resource. Therefore, the specific aims of this study were (i) to develop a model to predict total mushroom productivity associated with Pinus pinaster forests in northern Spain and (ii) to simulate the effects of different silvicultural and climatic scenarios on these yields. 2. Materials and Methods 2.1. Data All sporocarps were collected on a weekly basis during successive autumn mushroom seasons (late September to late December) from 2003 to 2014. As the average duration of the fruiting bodies varies from 4 to 20 days depending on the species [ 20 ], it was difficult to choose a sampling frequency that suited all species and did not distort production. A weekly sampling interval has been used by several authors in previous works [ 34 , 35 ] and was adopted for this study. Mushroom sporocarps were harvested, transported to the laboratory and stored at 4 ◦ C. Fresh weight and identification characteristics were recorded within 24 h of collection. The sporocarps were identified at the species level whenever possible following several keys [ 36 – 44 ]. As in previous works [ 45 , 46 ], samples that could only be identified to genus level were classified into genus taxa. Mushroom taxa names and authors were obtained from the Index Fungorum database (www.indexfungorum.org). After identification, the fresh weight and number of sporocarps were measured. The sporocarps were then dried in air-vented ovens at 35 ◦ C up to constant weight and subsequently dry-weighed in order to obtain comparable biomass values. Dried sporocarps were stored and used to complete identification based on microscopic key characteristics when necessary. At the end of 2014, dendrometric and dasometric variables were recorded in plots established in the nine mushroom sampling sites, which were previously established in 2003, in order to estimate the main silvicultural characteristics linked to mushroom production. In each plot (80 × 30 m), the diameter at breast height (dbh (cm)) was measured with callipers (precision 1 mm) for all trees in two perpendicular directions and total height (ht (m)) was recorded for a sample of 15 trees per plot. The height of the remaining trees in each plot was estimated by the h-d model [ 47 ]. This model was selected as the best fit in terms of R 2 parameter, after fitting different models obtained from the literature. Cores at breast height were collected to measure annual radial growth for the previous 10 full years using the Windrendro program [ 48 ]. Plot tree basal area (BA) for each year was calculated from annual radial growth and expanded to the hectare. There were no thinning operations during the years of the study. The ratio between the BA of the sampled year iand tree density (N, trees ha −1 ) was calculated as a parameter for calculating stand management properties. This term was named SMI (Stand Management Indicator) and expresses the relationship between two of the most important stand variables for describing forest management in a given stand. Biomass equations from [ 49 ] were used to estimate total biomass and the biomass of the different tree fractions, particularly root biomass, due to its possible relationship with mycorrhizal species. Total volume and stem volume were calculated with the Cubifor application [50]. Forests 2019,10, 206 4 of 18 Climate data for the 2003–2014 period were sourced from the closest meteorological stations [ 33 ]. Annual precipitation, mean, minimum and maximum temperatures and total rainfall and average temperature for the first three months of the autumn mushroom season (September, October, November) and the preceding month (August) were included in the analysis. In each plot, site characteristics such as altitude (meters above sea level (m.a.s.l)), slope (%) and exposure were recorded (Table 1). Soil samples were also taken from each plot in mid-May 2014, before the fruiting season. In each plot, six soil samples were extracted using a cylindrical (2 cm radius, 20 cm deep, 250 cm 3 ) soil borer [ 51 ], after carefully removing the overlying litter and humus layer. The samples were taken at a minimum distance of 30 cm from any tree trunk [ 51 ]. Soil texture was determined according to the International Society of Soil Science [ 52 ] classification criteria. The pH (soil:solution ratio of 1:2.5) was measured using a pH meter. Exchangeable cations were extracted with 0.1 M BaCl 2 . Ca and Mg concentrations were determined by atomic absorption spectrophotometry and K and Na by emission spectrophotometry [ 53 ]. Base saturation (BS) (cmol c kg −1 ) was calculated as the sum of Ca, Mg, K and Na concentrations. Total organic matter was analyzed by the Walkey and Black method. Soil total Carbon content was determined from organic matter using a conversion factor of 1.74 [54] and total Nitrogen content was calculated by the Kjeldahl method [55]. Table 1. Stand, site and soil characteristics of the plots. Plot Stand Site Soil BA (m2ha−1) N (trees ha−1) QMD (cm) Hdom (m) V (m3ha−1) Broot (mg ha−1)Alt. S (%) Exp. pH SOM (%) C/N AP MT 1 52.9 408 40.6 19.1 1.082 0.1438 985 0.1 SE 5.2 5.28 25.09 573 10.1 2 40.3 333 39.2 18.5 0.987 0.1321 985 0.1 SE 5.4 4.86 28.88 573 10.1 3 45.3 379 39.0 18.2 0.971 0.1302 985 0.1 SE 5.1 4.39 22.02 573 10.1 4 12.1 133 33.9 19.8 0.788 0.1037 700 0.011 N 6.6 0.53 11.34 455 11.8 5 14.8 183 32.0 20.0 0.682 0.0908 700 0.011 N 6.5 0.68 19.57 455 11.8 6 10.5 146 30.3 18.1 0.564 0.0772 700 0.011 N 6.0 0.80 14.83 455 11.8 7 18.6 129 42.8 18.5 1.234 0.1624 890 0.01 N 6.6 0.82 18.60 429 12.1 8 18.1 146 39.8 18.5 1.029 0.1372 890 0.01 N 6.6 1.67 23.55 429 12.1 9 21.4 183 38.5 18.4 0.954 0.1278 890 0.01 N 6.6 1.24 19.30 429 12.1 Note: BA (m 2 ha −1 ) is the basal area of the stand; N (trees ha −1 ) is stand density; QMD (cm) is the quadratic mean diameter; Hdom (m) is the dominant height; V (m 3 ha −1 ) is tree volume; Broot (Mg ha −1 ) is the biomass of the roots; Alt. is the altitude in meters above sea level (m a.s.l.); S (%) is the slope; Exp. is the exposure (SE: southeast; N: north); SOM (%) is the percentage of soil organic matter; C/N is the relationship between soil Carbon and Nitrogen content; AP is the total annual precipitation (mm), MT is the annual mean temperature (◦C). 2.2. Statistical Analysis Fresh weight, dry weight and number of sporocarps were analyzed for all mushroom species combined, but also for mycorrhizal/saprotrophic and nonedible/edible/marketed species groups. A linear mixed model (LMM) was fitted to relate total mushroom productivity to stand, site, soil and climate variables. Statistical assumptions were assessed using descriptive and graphical analyses of the different variables. To characterize the stand for each year sampled, the following parameters were included in the model: stand basal area (BA (m 2 ha −1 )), stand density (N (trees ha −1 )), quadratic mean diameter (QMD (cm)), dominant height (Hdom (m)), total tree volume (V (m 3 ha −1 )), total tree biomass (BT (Mg ha −1 )), root biomass (Broot (Mg ha −1 )) and the SMI parameter (Stand Management Indicator, the ratio between BA and N). Site characteristics included in the model were altitude (meters above sea level (m.a.s.l)), slope (%) and exposure. The physical and chemical soil characteristics considered were: pH, texture, macro- and micronutrients, base saturation (BS, cmol c kg −1 ), total organic matter (%), total Carbon (%), total Nitrogen (%) and C/N relationship. Finally, the climate variables included were mean, minimum and maximum temperatures ( ◦ C), total annual precipitation (mm) and the average temperature (◦C) and precipitation (mm) for August, September, October and November. To obtain the final mixed model, the relationship between total mushroom productivity and all possible combinations of the independent variables were tested. A correlation study was also carried out to reduce the combinations among the variables. The proc corr procedure in the SAS 9.2 statistical program (Cary, NC, USA) was used for this [ 56 ]. To select the best model, we looked for one that combined the most interesting/most useful variables with a reasonable biological interpretation. Forests 2019,10, 206 5 of 18 Only the final selected model is presented here (Equation (1)). After obtaining the variance for each sampled year, LMM emerged as the best option for estimating total mushroom productivity, due to mushroom production variability in the sampled years: LnYij +1=β0+β1·LnSMIij +δi+β2·TsPsij +β3·TnPnij +εij (1) where Yij is the total fresh-weight mushroom productivity (kg fw ha −1 ) of plot iand year j. β0 is the intercept of the model. SMI ij is the stand management indicator (SMI) of the plot iand year j , where: SMIij =BAij(m2ha−1) Ni(trees ha−1) β1is the linear effect of the logarithm of the stand management indicator (SMIij). TsPsij is the September average Temperature * precipitation for plot iand year j. β2is the linear effect of the September average Temperature * precipitation variable. TnPnij is the November average Temperature * precipitation for plot iand year j. β3is the linear effect of the November average Temperature * precipitation variable. δiis the random effect of plot i, with: δi→N(0, σ2 δ) εij is the experimental error, verifying the following assumptions: εij →N(0, σ2 j) Covεij,εi0j0=(σjσj0ρ|j−j0|i f i =i0and j 6=j0 0i f i 6=i0 Therefore, the LMM model included 14 parameters of variance: σj is the variance between-plots; σ2 j with (j= 1–12) for within-plot error and ρ for the correlation between consecutive years. Restricted maximum likelihood (REML) was used to fit the mixed model, using PROC MIXED in the SAS 9.2 statistical program (Cary, NC, USA) [56]. The adequacy of the model was analyzed by testing the equation parameters simultaneously to compare the actual and predicted values (Equation (2)) of the total mushroom productivity variable. The bias correction factor was used to correct the predictions for the back-transformation bias to the original scale [ 57 ]. The simultaneous F-test of a = 0 and b = 1 was used to determine the adequacy of the model, as the regression between actual and predicted values should be a 45 ◦ line through the origin [58]. actual =a+b predicted (2) where actual is the observed value of the mushroom biomass, predicted is the value obtained from the model, and aand bare the parameters to be estimated. 2.3. Simulation Scenarios The fitted total productivity model was uploaded to the SiManFor platform to evaluate the effects of different silvicultural and climate scenarios on mushroom productivity. The silvicultural scenarios consisted of three thinning intensities: light (G10, 10% of basal area), moderate (G25, 25% of basal area) and heavy (G35, 35% of basal area). The type of thinnings was from below. The silvicultural scenarios were performed with the application of only one thinnings activity at the age of 63 years in the three cases (G10, G25 and G35) to compare among them. These were chosen to reflect criteria established in the silvicultural regulations for Pinus plantations in Spain [ 59 ] and because they coincided with the thinning intensities applied in a thinning experiment already installed in the area for the ‘Sustainable Innovative Mobilisation of Wood’ (SIMWOOD) research project [60]. The climatic scenarios tested were real (c1, c2, c3), wetter than real (c4, c5, c6), the wettest (c7, c8, c9 ), drier than real (c10, c11, c12), the driest (c13, c14, c15) and hotter than real (c16, c17, c18). So, c1 = precipitation and average temperature values in September (32 mm; 16 ◦ C) and in November Forests 2019,10, 206 6 of 18 (45 mm; 6 ◦ C); c2 = minimum precipitation and temperature values in September (1 mm; 14 ◦ C) and in November (12 mm; 4 ◦ C); c3 = maximum precipitation and temperature values in September (73 mm; 18 ◦ C) and in November (90 mm; 9 ◦ C). Precipitation for wetter scenarios (c4, c5 and c6) was 20 mm higher than average real scenario (c1) values for both months (Sep/Nov). Temperatures for the c4 scenario were equal to average real scenario (c1), while scenarios c5 and c6 represented the average temperature +2 ◦ C and − 2 ◦ C, respectively. The wettest scenarios were defined by adding 30 mm more precipitation in both months to c4 (c7), c5 (c8) and c6 (c9). Similarly, precipitation for drier scenarios (c10, c11 and c12) was 20 mm below the average real scenario (c1) values in both months (Sep/Nov). Temperatures for scenario c10 were equal to real scenario (c1), whereas scenarios c11 and c12 were the average temperature +2 ◦ C and − 2 ◦ C, respectively. The driest scenarios were defined by subtracting 20 mm precipitation from c10 (c13), c11 (c14) and c12 (c15). Hotter scenarios were characterized by precipitation of the c1 scenario and by adding from 1 ◦ C to 3 ◦ C to the average temperature in September and November: average +1 ◦C (c16), +2 ◦C (c17) and +3 ◦C (c18). 3. Results 3.1. Mushroom Productivity The average total mushroom fresh weight ranged from 2.25 kg fw ha −1 (0.17 kg dw ha −1 and 2833 ns ha −1 ) in 2005 to 350.51 kg fw ha −1 (36.13 kg dw ha −1 and 46,356 ns ha −1 ) in 2014 (Figure 1a). Most of the total productivity belonged to mycorrhizal ecological group, being more than 80% of the total mushroom fresh weight mycorrhizal fresh weight. In both groups, total productivity and mycorrhizal, maximum fw and dw values were observed in 2014, followed by 2013 and 2006, when almost 300 kg fw ha −1 of total productivity were harvested. In 2003, total and mycorrhizal fresh weight productivity exceeded 250 kg ha −1 (Figure 1a). Average total productivity and average mycorrhizal fresh weight were 167.50 kg fw ha −1 and 146.41 kg fw ha −1 , respectively, in the 12-year study period. In terms of economic value, average edible productivity and marketed were 46.10 kg fw ha −1 and 20.78 kg fw ha −1 , respectively. Edible productivity corresponded to 27.5% to the total productivity and marketed productivity corresponded to 45.1% of the edible mushroom productivity (Figure 1b). Results showed high values of different edible and marketed species such Lactarius deliciosus (L.) S.F. Gray, with more than 24 kg fw ha −1 in 2013. Mushroom productivity in terms of fresh and dry weight varied greatly among the sampled years and plots (Table 2). Forests 2019, 10, x FOR PEER REVIEW 7 of 18 (a) (b) Figure 1. (a) Average total and mycorrhizal fresh weight (kg ha−1) production by sampled year. (b) Average total, edible and marketed fresh weight (kg ha−1) production by sampled year. Table 2. Total fresh weight (fw kg ha−1), dry weight (dw kg ha−1) and number of sporocarps per hectare (ns ha−1) by plot and sampled year. Plot fw (kg ha−1) dw (kg ha−1) ns ha−1 fw (kg ha−1) dw (kg ha−1) ns ha−1 fw (kg ha−1) dw (kg ha−1) ns ha−1 2003 2004 2005 1 241.5 17.6 15000 73.4 10.4 12000 2.9 0.3 4000 2 392.7 27.3 15200 101.7 24.8 12600 0.8 0.0 1600 3 169.7 11.1 11200 181.5 21.1 10000 3.1 0.2 2900 2006 2007 2008 1 444.8 30.3 30400 - - - 13.4 1.5 80800 2 436.9 27.5 24800 - - - 22.6 2.4 63000 3 729.0 59.9 32900 - - - 26.9 2.8 75500 4 17.0 1.9 22800 62.1 9.5 30800 44.1 4.5 66500 5 6.8 0.9 10600 4.8 0.7 7800 5.0 0.6 17200 6 166.4 17.8 20700 118.7 18.3 27100 120.4 10.2 40700 7 391.9 38.8 50900 454.7 79.4 112900 201.1 17.7 88100 8 198.5 19.5 224500 216.2 39.7 110400 64.0 6.0 49500 9 290.0 29.2 26900 245.7 42.7 62700 92.5 9.7 33600 Figure 1. Cont. Forests 2019,10, 206 7 of 18 Forests 2019, 10, x FOR PEER REVIEW 7 of 18 (a) (b) Figure 1. (a) Average total and mycorrhizal fresh weight (kg ha−1) production by sampled year. (b) Average total, edible and marketed fresh weight (kg ha−1) production by sampled year. Table 2. Total fresh weight (fw kg ha−1), dry weight (dw kg ha−1) and number of sporocarps per hectare (ns ha−1) by plot and sampled year. Plot fw (kg ha−1) dw (kg ha−1) ns ha−1 fw (kg ha−1) dw (kg ha−1) ns ha−1 fw (kg ha−1) dw (kg ha−1) ns ha−1 2003 2004 2005 1 241.5 17.6 15000 73.4 10.4 12000 2.9 0.3 4000 2 392.7 27.3 15200 101.7 24.8 12600 0.8 0.0 1600 3 169.7 11.1 11200 181.5 21.1 10000 3.1 0.2 2900 2006 2007 2008 1 444.8 30.3 30400 - - - 13.4 1.5 80800 2 436.9 27.5 24800 - - - 22.6 2.4 63000 3 729.0 59.9 32900 - - - 26.9 2.8 75500 4 17.0 1.9 22800 62.1 9.5 30800 44.1 4.5 66500 5 6.8 0.9 10600 4.8 0.7 7800 5.0 0.6 17200 6 166.4 17.8 20700 118.7 18.3 27100 120.4 10.2 40700 7 391.9 38.8 50900 454.7 79.4 112900 201.1 17.7 88100 8 198.5 19.5 224500 216.2 39.7 110400 64.0 6.0 49500 9 290.0 29.2 26900 245.7 42.7 62700 92.5 9.7 33600 Figure 1. ( a ) Average total and mycorrhizal fresh weight (kg ha −1 ) production by sampled year. (b) Average total, edible and marketed fresh weight (kg ha−1) production by sampled year. Table 2. Total fresh weight (fw kg ha −1 ), dry weight (dw kg ha −1 ) and number of sporocarps per hectare (ns ha−1) by plot and sampled year. Plot fw (kg ha−1) dw (kg ha−1)ns ha−1fw (kg ha−1) dw (kg ha−1)ns ha−1fw (kg ha−1) dw (kg ha−1)ns ha−1 2003 2004 2005 1 241.5 17.6 15,000 73.4 10.4 12,000 2.9 0.3 4000 2 392.7 27.3 15,200 101.7 24.8 12,600 0.8 0.0 1600 3 169.7 11.1 11,200 181.5 21.1 10,000 3.1 0.2 2900 2006 2007 2008 1 444.8 30.3 30,400 - - - 13.4 1.5 80,800 2 436.9 27.5 24,800 - - - 22.6 2.4 63,000 3 729.0 59.9 32,900 - - - 26.9 2.8 75,500 4 17.0 1.9 22,800 62.1 9.5 30,800 44.1 4.5 66,500 5 6.8 0.9 10,600 4.8 0.7 7800 5.0 0.6 17,200 6 166.4 17.8 20,700 118.7 18.3 27,100 120.4 10.2 40,700 7 391.9 38.8 50,900 454.7 79.4 112,900 201.1 17.7 88,100 8 198.5 19.5 224,500 216.2 39.7 110,400 64.0 6.0 49,500 9 290.0 29.2 26,900 245.7 42.7 62,700 92.5 9.7 33,600 2009 2010 2011 1 115.3 10.9 128,900 133.6 15.2 75,100 117.8 11.0 81,300 2 95.7 7.4 128,100 97.0 8.0 57,300 104.7 8.1 101,100 3 26.2 2.3 76,500 65.4 10.5 57,800 102.8 8.1 102,700 4 159.2 15.3 28,600 - - - 45.5 4.3 54,900 5 2.8 0.3 3400 - - - 13.2 1.4 24,900 6 410.0 42.5 25,400 - - - 22.4 2.2 33,700 7 122.0 13.1 28,600 - - - 137.0 10.9 50,600 8 100.8 9.2 33,800 - - - 27.3 2.5 24,100 9 84.5 8.6 24,100 - - - 78.5 8.3 22,600 2012 2013 2014 1 70.3 5.1 18,800.0 497.7 51.9 113,500 684.5 61.5 94,300 2 22.0 2.0 16,200.0 193.9 20.9 138,400 396.8 43.5 66,400 3 37.7 2.3 20,400.0 67.4 7.0 65,800 341.2 48.5 54,300 4 57.3 5.6 45,500.0 150.8 22.3 14,800 124.8 17.8 36,400 5 170.4 19.1 50,600.0 224.3 28.6 18,300 102.2 14.9 12,500 6 65.8 6.9 53,000 354.6 53.7 17,500 162.9 24.2 46,800 7 467.0 47.7 101,100 478.3 65.3 34,300 648.5 56.8 52,400 8 179.0 17.6 52,500 298.3 38.0 35,900 460.6 38.8 34,700 9 98.5 10.1 48,500 421.4 61.7 26,800 233.2 19.1 19,400 Note: “-”indicates a lack of data for specific years. Forests 2019,10, 206 8 of 18 3.2. Mushroom Productivity Model Total fresh-weight productivity was related to the different stand, site, soil and weather variables in the sampled plots. Significant correlations were found among most stand variables, soil variables such as the Nitrogen/Potasium ratio and weather variables such as September and November precipitation and September average temperature (Table 3). The mixed model showed that the SMI parameter, the interaction between September precipitation and average temperature and the interaction between November precipitation and average temperature were significant at α < 0.05 (Table 4). Linear mixed-model results showed that higher SMI values, which corresponded to higher average basal area per tree, dealt with higher total mushroom productivity, after accounting for the weather variables. Thus, a doubling of SMI was associated with a 2 0.8797 = 1.84-fold increase in the median of total mushroom productivity. Higher rainfall during September and November, together with higher temperatures in early autumn, were also associated with an increase in total mushroom productivity. The fitted mixed model estimated a variance for each year sampled (Table 5) and included the variance between plots and the correlation between consecutive years. On the regression line between actual and predicted values, which describe the adequacy of the model, the intercept term was not significantly different from zero and the slope was not significantly different from one. The value of the bias correction factor [ 57 ] was 0.9736 and the pseudo R 2 statistics of the regression between actual and predicted values (pseudo R2= 0.5011) indicated good model performance. Table 3. Results indicating significant correlation levels between stand, site, soil and climatic variables and total mushroom productivity. Variable Correlation Correlation Level Stand N (trees ha−1)−0.21536 0.0713 QMD (cm) 0.43380 0.0002 Hdom (m) −0.27097 0.0223 V (m3ha−1)0.46478 <0.0001 Vs (m3ha−1)0.46474 <0.0001 Bbranches > 7 (Mg ha −1 ) 0.46548 <0.0001 Bbranches < 2 (Mg ha −1 ) 0.45570 <0.0001 Broot (Mg ha−1)0.45907 <0.0001 SMI parameter 0.47147 <0.0001 Soil pH 0.20622 0.0845 NMg −0.20905 0.0802 NK −0.25277 0.0334 Climatic SP (mm) 0.49268 <0.0001 NP (mm) 0.26051 0.0282 MT (◦C) 0.22344 0.0611 ST (◦C) 0.31817 0.0069 OT (◦C) 0.22684 0.0571 NT (◦C) 0.20052 0.0936 Note: N (trees ha −1 ) is stand density; QMD (cm) is the quadratic mean diameter; Hdom (m) is the dominant height; V (m 3 ha −1 ) is tree volume; vs. (m 3 ha −1 ) is stem volume; Bbranches > 7 (Mg ha −1 ) is the biomass of branches with diameter greater than 7 cm; Bbranches < 2 (Mg ha −1 ) is the biomass of branches with diameter smaller than 2 cm; Broot (Mg ha −1 ) is the biomass of the roots; SMI is the Stand Management Indicator (ratio between basal area of the stand (BA) and tree density (N)); NMg is the ratio between Nitrogen and Magnesium soil content; NK is the ratio between Nitrogen and Potassium soil content; SP and NP indicate September and November precipitation in mm, respectively; MT is the mean temperature ( ◦ C); ST, OT and NT indicate September, October and November average temperatures in ◦C, respectively. Forests 2019,10, 206 9 of 18 Table 4. Solution and Type 3 tests for fixed effects in the mixed model. Effect Estimate Standard Error DenDF F-Value Pr > F Incercept 10.5544 1.2807 8 8.24 <0.0001 Ln SMI 3.1786 0.5619 69 5.66 <0.0001 TsPs 0.000906 0.000173 69 5.24 <0.0001 TnPn 0.001926 0.000129 69 14.96 <0.0001 Note: DenDF is the denominator degrees of freedom, respectively; F-Value is the value of the Fstatistic; Pr > F is the P-value associated with the previous F-statistic. Ln SMI is the logarithm of the SMI parameter (Stand Management Indicator (ratio between basal area of the stand (BA) and tree density (N))); TsPs is the average Temperature * Precipitation variable for September; TnPn is the average Temperature * Precipitation variable for November. Table 5. Variance parameters of the linear mixed model. Estimated Variance Parameter Plot 0.1360 2003 0.3898 2004 0.7622 2005 15.0219 2006 0.7538 2007 0.5748 2008 0.9821 2009 1.7094 2010 0.02406 2011 0.4127 2012 1.3938 2013 3.6923 2014 0.01082 ρ0.6119 3.3. Simulations in SiManFor Platform Simulation results showed that total mushroom productivity before thinning was smaller when precipitation was lower (c2, c10–c15) and increased in wetter conditions (c3, c4–c9) and hotter scenarios (c16–c18) (Table 6). The highest simulated productivity occurred in the c8 scenario, which corresponded to the highest precipitation and temperature scenario studied. Forests 2019,10, 206 16 of 18 16. Castaño, C.; Lindahl, B.D.; Alday, J.G.; Hagenbo, A.; de Aragón, J.M.; Parladé, J.; Pera, J.; Bonet, J.A. Soil microclimate changes affect soil fungal communities in a Mediterranean pine forest. New Phytol. 2018 ,220, 1211–1221. [CrossRef] [PubMed] 17. Karavani, A.; De Cáceres, M.; de Aragón, J.M.; Bonet, J.A.; de Miguel, S. Effect of climatic and soil moisture conditions on mushroom productivity and related ecosystem services in Mediterranean pine stands facing climate change. Agric. For. Meteorol. 2018,248, 432–440. [CrossRef] 18. Hernández-Rodríguez, M.; de Miguel, S.; Pukkala, T.; de Rueda, J.A.O.; Martín-Pinto, P. Climate-sensitive models for mushroom yields and diversity in Cistus ladanifer scrublands. Agric. For. Meteorol. 2015 ,213, 173–182. [CrossRef] 19. Alday, J.G.; Bonet, J.A.; de Rueda, J.A.O.; de Aragón, J.M.; Martín-Pinto, P.; de Miguel, S.; Hernández-Rodríguez, M.; Martínez-Peña, F. Record breaking mushroom yields in Spain. Fungal Ecol. 2017,26, 144–146. [CrossRef] 20. Vogt, K.A.; Bloomfield, J.; Ammirati, J.F.; Ammirati, S.R. Sporocarp production by basidiomycetes, with emphasis on forest ecosystems. In The Fungal Community. Its Organization and Role in the Ecosystem, 2nd ed.; Dekker, M., Ed.; Marcel Dekker: New York, NY, USA, 1992; pp. 563–581. 21. De Aragón, J.M.; Bonet, J.A.; Fischer, C.R.; Colinas, C. Productivity of ectomycorrhizal and selected edible saprotrophic fungi in pine forests of the pre-Pyrenees mountains, Spain: Predictive equations for forest management of mycolocial resources. For. Ecol. Manag. 2007,252, 239–256. [CrossRef] 22. Vasquez Gassibe, P.; Fraile, R.; Hernández-Rodríguez, M.; de Rueda, J.A.O.; Bravo, F.; Martín-Pinto, P. Post-fire production of mushrooms in Pinus pinaster forests using classificatory models. J. For. Res. 2014 ,19, 348–356. [CrossRef] 23. Dahl, F.A.; Galteland, T.; Gjelsvik, R. Statistical modelling of wildwood mushroom abundance. Scand. J. For. Res. 2008,23, 224–249. [CrossRef] 24. Martínez-Peña, F.; de Miguel, S.; Pukkala, T.; Bonet, J.A.; Ortega-Martínez, P.; Aldea, J.; de Aragón, J.M. Yield models for ectomycorrhizal mushrooms in Pinus sylvestris forests with special focus on Boletus edulis and Lactarius group deliciosus.For. Ecol. Manag. 2012,282, 63–69. [CrossRef] 25. Bonet, J.A.; Palahí, M.; Colinas, C.; Pukkala, T.; Fischer, C.R.; Miina, J.; de Aragón, J.M. Modelling the production and species richness of wild mushrooms in pine forests of Central Pyrenees in north-eastern Spain. Can. J. For. Res. 2010,40, 347–356. [CrossRef] 26. Tahvanainen, V.; Miina, J.; Kurttila, M.; Salo, K. Modelling the yields of marketed mushrooms in Picea abies stands in eastern Finland. For. Ecol. Manag. 2016,362, 79–88. [CrossRef] 27. De Miguel, S.; Bonet, J.A.; Pukkala, T.; de Aragón, J.M. Impact of forest management intensity on landscape-level mushroom productivity: A regional model-based scenario analysis. For. Ecol. Manag. 2014,330, 218–227. [CrossRef] 28. Taye, Z.M.; Martínez-Peña, F.; Bonet, J.A.; de Aragón, J.M.; de Miguel, S. Meteorological conditions and site characteristics driving edible mushroom production in Pinus pinaster forests of Central Spain. Fungal Ecol. 2016,23, 30–41. [CrossRef] 29. Fox, J.C.; Ades, P.K.; Bi, H. Stochastic structure and individual-tree growth models. For. Ecol. Manag. 2001 , 154, 261–276. [CrossRef] 30. Bravo, F.; Rodríguez, F.; Ordóñez, A.C. A web-based application to simulate alternatives for sustainable forest management: SIMANFOR. For. Syst. 2012,21, 4–8. [CrossRef] 31. Agreda, T.; Águeda, B.; Olano, J.M.; Vicente-Serrano, S.M.; Fernández-Toirán, M. Increased evapotranspiration demand in a Mediterranean climate might cause a decline in fungal yields under global warming. Glob. Chang. Biol. 2015,21, 3499–3510. [CrossRef] [PubMed] 32. Büntgen, U.; Egli, S.; Galván, J.D.; Diez, J.M.; Aldea, J.; Latorre, J.; Martínez-Peña, F. Drought-induced changes in the phenology, productivity and diversity of Spanish fungi. Fungal Ecol. 2015 ,16, 6–18. [CrossRef] 33. Vásquez, P.; de Rueda, J.O.; Martín-Pinto, P.P. pinaster under extreme ecological conditions provides high fungal production and diversity. For. Ecol. Manag. 2015,337, 161–173. [CrossRef] 34. Ohenoja, E.; Koistinen, R. Fruit body production of larger fungi in Finland. 2. Edible fungi in northern Finland 1976–1978. Ann. Bot. Fenn. 1984,21, 357–366. 35. Baptista, P.; Martins, A.; Tavares, R.M.; Lino-Neto, T. Diversity and fruiting pattern of macrofungi associated with chestnut (Castanea sativa) in the Trás-os-Montes region (Northeast Portugal). Fungal Ecol. 2010 ,3, 9–19. [CrossRef] Forests 2019,10, 206 17 of 18 36. Moser, M.; Plant, S.; Kibby, G. Keys to Agarics and Boleti: Polyporales, Boletales, Agaricales, Russulales; Roger Phillips: London, UK, 1983; pp. 1–535. 37. Breitenbach, J. Champignons de Suisse Tome 1. Les Ascomycètes; Mykologia: Lucerne, Switzerland, 1984; pp. 1–310. 38. Breitenbach, J. Champignons de Suisse Tome 2. Champignons sans Lames; Mykologia: Lucerne, Switzerland, 1986; pp. 1–412. 39. Breitenbach, J. Champignons de Suisse Tome 3. Bolets et Champignons àLames: Première Partie; Mykologia: Lucerne, Switzerland, 1991; pp. 1–364. 40. Breitenbach, J. Champignons de Suisse Tome 4. Champignons àlames. Deuxième Partie; Mykologia: Lucerne, Switzerland, 1995; pp. 1–371. 41. Breitenbach, J. Champignons de Suisse Tome 5. Champignons àlames. Troisième Partie; Mykologia: Lucerne, Switzerland, 2000; pp. 1–340. 42. Breitenbach, J. Champignons de Suisse Tome 6. Russulaceae. Actaires. Russules; Mykologia: Lucerne, Switzerland, 2005; pp. 1–263. 43. Bon, M. Guía de Campo de los Hongos de Europa; Omega: Barcelona, Spain, 1987; pp. 1–352. 44. Knudsen, H.; Vesterholt, J. Funga Nordica: Agaricoid, Boletoid and Cyphelloid Genera; Nordsvamp: Copenhagen, Denmark, 2008; pp. 1–965. 45. Bonet, J.; Fischer, C.; Colinas, C. The relationship between forest age and aspect on the production of sporocarps of ectomycorrhizal fungi in Pinus sylvestris forest of the central Pyrenees. For. Ecol. Manag. 2004 , 203, 157–175. [CrossRef] 46. Martín-Pinto, P.; Pajares, J.; Diez, J. In vitro effects of four ectomycorrhizal fungi, Boletus edulis,Rhizopogon roseolus,Laccaria laccata and Lactarius deliciosus on Fusarium damping off in Pinus nigra seedlings. New For. 2006,32, 323–334. [CrossRef] 47. Meyer, H.A. A mathematical expression for height curves. J. For. 1940,38, 415–420. 48. Regent Instruments Inc. WinDENDROtm; Regent Instruments Inc.: Québec, QC, Canada, 2003. 49. Ruiz-Peinado, R.; Río, M.; Montero, G. New models for estimating the carbon sink capacity of Spanish softwood species. For. Syst. 2011,20, 176–188. [CrossRef] 50. Rodríguez, F.; Broto, M.; Lizarralde, I. CubiFor: Herramienta para cubicar, clasificar productos y calcular biomasa y CO2en masas forestales de Castilla y León. Rev. Montes 2008,95, 33–39. 51. De la Varga, H.; Agueda, B.; Martínez-Pena, F.; Parlade, J.; Pera, J. Quantification of extraradical soil mycelium and ectomycorrhizas of Boletus.For. Ecol. Manag. 2015,337, 161–173. [CrossRef] 52. ISSS-ISRIC-FAO. World Reference Base for Soil Resources. Draft; Wageningen: Rome, Italy, 1994; pp. 1–133. 53. MAPA. Métodos Oficiales de Análisis. Tomo III; Ministerio de Agricultura, Pesca y Alimentación: Madrid, Spain, 1994; pp. 1–662. 54. Nicholson, G. Methods of Soil, Plant and Water Analysis; New Zealand Forest Research Institute: Rotorua, New Zealand, 1984; pp. 1–24. 55. Benton, J.J.; Wolf, J.B.; Mills, H.A. Plant Analysis Handbook, a Practical Sampling, Preparation, Analysis and Interpretation Guide; Micro-Macro Publishing: Athens, GA, USA, 1991; pp. 1–213. 56. Sas Institute Inc. SAS/Stattm User’s Guide, Relase 9.4; Sas Institute Inc.: Cary, NC, USA, 2016; pp. 1–550. 57. Baskerville, G.L. Use of logarithmic regression in the estimation of plant biomass. Can. J. For. Res. 1972 ,2, 49–53. [CrossRef] 58. Huang, S.; Yang, Y.; Wang, Y. A critical look at procedures for validating growth and yield models. In Modelling Forest Systems; Amaro, A., Reed, D., Soares, P., Eds.; CAB International: Wallingford, UK, 2003; pp. 271–293. 59. Rio, M.; López, E.; Montero, G. Manual de Gestión para Masas Procedentes de Repoblación de Pinus pinaster Ait, Pinus sylvestris L. y Pinus nigra Arn. en Castilla y León; Consejería de Medio Ambiente, Junta de Castilla y León: Valladolid, Spain, 2006; pp. 1–102. 60. Bravo, F.; Herrero, C.; Ordóñez, A.C.; De La Parra, B.; Cuesta, J.; Santos, M.; Ruano, I.; Manso, R. Instalación de Ensayos de Gestión Forestal Adaptativa en Distintos Tipos de Pinares en el Bosque Modelo de Palencia; Consejería de Medio Ambiente, Junta de Castilla y León: Valladolid, Spain, 2015; pp. 1–15. 61. O’Dell, T.E.; Ammirati, J.F.; Schreiner, E.G. Species richness and abundance of ectomycorrhizal basidiomycete sporocarps on a moisture gradient in the Tsuga heterophylla zone. Can. J. Bot. 1999 ,77, 1699–1711. [CrossRef] Forests 2019,10, 206 18 of 18 62. Bonet, J.A.; de-Miguel, S.; de Aragón, J.M.; Pukkala, T.; Palahí, M. Immediate effect of thinning on the yield of Lactarius group deliciosus in Pinus pinaster forests in North-Eastern Spain. For. Ecol. Manag. 2012 ,265, 211–217. [CrossRef] 63. Straatsma, G.; Ayer, F.; Egli, S. Species richness, abundance and phenology of fungal fruit bodies over 21 years in a Swiss forest plot. Mycol. Res. 2001,105, 515–523. [CrossRef] 64. Straatsma, G.; Krisai-Greilhuber, I. Assemblage structure, species richness, abundance, and distribution of fungal fruit bodies in a seven year plot-based survey near Vienna. Mycol. Res. 2003 ,107, 632–640. [CrossRef] [PubMed] 65. Krebs, C.J.; Carrier, P.; Boutin, S.; Boonstra, R.; Hofer, E. Mushroom crops in relation to weather in southeastern Yukon. Botany 2008,86, 1497–1502. [CrossRef] 66. Berraondo, I.; Herrero, C.; de la Parra, B.; Olaizola, J.; Pando, V.; Oria de Rueda, J.A. Modelización de la producción de Tricholoma portentosum (Fr.) Quél. en masas de Pinus sylvestris L. In Proceedings of the V Congreso Forestal Español, Ávila, Spain, 21–25 September 2009. 67. Martínez-Peña, F. Producción y Aprovechamiento de Boletus edulis Bull.: Fr. en un Bosque de Pinus sylvestris L.; Consejería de Medio Ambiente. Junta de Castila y León: Valladolid, Spain, 2003; pp. 1–134. 68. Kauserud, C.L.; Stige, J.O.; Vik, J.O.; Okland, R.H.; Hoiland, K.; Stenseth, N.C. Mushroom fruiting and climate change. Proc. Natl. Acad. Sci. USA 2008,105, 3811–3814. [CrossRef] [PubMed] 69. Shaw, P.J.A.; Kibby, C.; Mayes, J. Effects of thinning treatment on an ectomycorrhizal succession under Scots pine. Mycol. Res. 2003,107, 317–328. [CrossRef] [PubMed] 70. Amaranthus, M.P.; Page-Dumroese, D.; Harvey, A.; Cazares, E.; Bednar, L.F. Soil Compaction and Organic Matter Affect Conifer Seedling Nonmycorrhizal and Ectomycorrhizal Root Tip Abundance and Diversity; Research Paper PNW-RP-494; USDA Forest Service: Portland, OR, USA, 1996; pp. 1–12. 71. Kropp, B.R.; Albee, S. The effects of silvicultural treatments on occurrence of mycorrhizal sporocarps in a Pinus contorta forest: A preliminary survey. Biol. Conserv. 1996,78, 313–318. [CrossRef] 72. Pilz, D.; Molina, R.; Mayo, J. Effects of thinning young forests on Chanterelle mushroom production. J. For. 2006,104, 9–14. [CrossRef] 73. Tomao, A.; Bonet, J.A.; de Aragón, J.M.; de Miguel, S. Is silviculture able to enhance wild forest mushroom resources? Current knowledge and future perspectives. For. Ecol. Manag. 2017,402, 102–114. [CrossRef] 74. Salerni, E.; Perini, C. Experimental study for increasing productivity of Boletus edulis s.l. in Italy. For. Ecol. Manag. 2004,201, 161–170. [CrossRef] 75. Nogués Bravo, D.; Araújo, M.B.; Lasanta, T.; López Moreno, J.L. Climate change in Mediterranean mountains during the 21st century. Ambio 2008,37, 280–285. [CrossRef] 76. IPCC. Summary for Policymakers. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change. In Climate Change 2013: The Physical Science Basis; Stocker, T.F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S.K., Boschung, J., Nauels, A., Xia, Y., Bex, V., Midgley, P.M., Eds.; Cambridge University Press: Cambridge, UK, 2013; pp. 1–1535. © 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).