Predicting Mushroom Productivity from Long-Term Field-Data Series in Mediterranean Pinus pinaster Ait . Forests in the Context of Climate Change

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.


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.
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.

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][37][38][39][40][41][42][43][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 i and 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].
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].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).

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.
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: where Y ij is the total fresh-weight mushroom productivity (kg fw ha −1 ) of plot i and year j.β 0 is the intercept of the model.SMI ij is the stand management indicator (SMI) of the plot i and year j, where: β 1 is the linear effect of the logarithm of the stand management indicator (SMI ij ).
TsPs ij is the September average Temperature * precipitation for plot i and year j.
β 2 is the linear effect of the September average Temperature * precipitation variable.TnPn ij is the November average Temperature * precipitation for plot i and year j.
β 3 is the linear effect of the November average Temperature * precipitation variable.δ i is the random effect of plot i, with: δ i → N(0, σ 2 δ ) ε ij is the experimental error, verifying the following assumptions: 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].
where actual is the observed value of the mushroom biomass, predicted is the value obtained from the model, and a and b are the parameters to be estimated.

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].

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).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.

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 R 2 = 0.5011) indicated good model performance.

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.The effect of thinning was assessed for each silvicultural scenario (G10, G25 and G35) in the c1 climatic scenario (Figure 2).Mushroom production decreased after thinning in G10 (Figure 2a), but not in G25 and G35 (Figure 2b,c), indicating a positive effect from thinning at moderate or high thinning intensity.
The effect of thinning was assessed for each silvicultural scenario (G10, G25 and G35) in the c1 climatic scenario (Figure 2).Mushroom production decreased after thinning in G10 (Figure 2a), but not in G25 and G35 (Figure 2b,c The effect of thinning intensity was evaluated in different climatic scenarios (Figure 3).Comparison of the three silvicultural scenarios (G10, G25, G35) in the c1, c2 and c3 climatic scenarios revealed that total mushroom productivity increased after thinning in G25 and G35 and also increased with thinning intensity (Figure 3a,b,c).The highest increment appeared in the climatic conditions of the c3 scenario (the highest precipitation with the hottest conditions) (Figure 3c).A similar pattern was found in wetter (c4 (Figure 3d) and c7 (Figure 3e)), drier (c10 (Figure 3f), c13 (Figure 3g)) and hotter climatic scenarios (c16 (Figure 2h), c18 (Figure 3i)).More specifically, the The effect of thinning intensity was evaluated in different climatic scenarios (Figure 3).Comparison of the three silvicultural scenarios (G10, G25, G35) in the c1, c2 and c3 climatic scenarios revealed that total mushroom productivity increased after thinning in G25 and G35 and also increased with thinning intensity (Figure 3a-c).The highest increment appeared in the climatic conditions of the c3 scenario (the highest precipitation with the hottest conditions) (Figure 3c).The effect of thinning was assessed for each silvicultural scenario (G10, G25 and G35) in the c1 climatic scenario (Figure 2).Mushroom production decreased after thinning in G10 (Figure 2a), but not in G25 and G35 (Figure 2b,c), indicating a positive effect from thinning at moderate or high thinning intensity.The effect of thinning intensity was evaluated in different climatic scenarios (Figure 3).Comparison of the three silvicultural scenarios (G10, G25, G35) in the c1, c2 and c3 climatic scenarios revealed that total mushroom productivity increased after thinning in G25 and G35 and also increased with thinning intensity (Figure 3a,b,c).The highest increment appeared in the climatic conditions of the c3 scenario (the highest precipitation with the hottest conditions) (Figure 3c).A similar pattern was found in wetter (c4 (Figure 3d) and c7 (Figure 3e)), drier (c10 (Figure 3f), c13 (Figure 3g)) and hotter climatic scenarios (c16 (Figure 2h), c18 (Figure 3i)).More specifically, the A similar pattern was found in wetter (c4 (Figure 3d) and c7 (Figure 3e)), drier (c10 (Figure 3f), c13 (Figure 3g)) and hotter climatic scenarios (c16 (Figure 3h), c18 (Figure 3i)).More specifically, the effect of higher thinning intensity was very positive in wetter scenarios, slightly positive in hotter scenarios (Figure 3h,i) and very negative in drier scenarios (Figure 3b,f,g).
In contrast, comparison of the climatic scenarios within one silvicultural treatment, for instance G35, showed that the highest total mushroom productivity was obtained in wetter-and hotter-than-average conditions (Figure 4a,b).In drier scenarios, c12 showed the smallest decrease in production due to the effect of temperature and production was lowest in the driest conditions (Figure 4c).In hotter scenarios (Figure 4d), total mushroom productivity increased slightly with respect to the c1 scenario.
Forests 2019, 10, x FOR PEER REVIEW 12 of 18 effect of higher thinning intensity was very positive in wetter scenarios, slightly positive in hotter scenarios (Figure 3h,i) and very negative in drier scenarios (Figure 3b,f,g).
In contrast, comparison of the climatic scenarios within one silvicultural treatment, for instance G35, showed that the highest total mushroom productivity was obtained in wetter-and hotter-thanaverage conditions (Figure 4 a,b).In drier scenarios, c12 showed the smallest decrease in production due to the effect of temperature and production was lowest in the driest conditions (Figure 4c).In hotter scenarios (Figure 4d), total mushroom productivity increased slightly with respect to the c1 scenario.

Discussion
In this study, mushroom production in Pinus pinaster stands in northern Spain was analyzed.A predictive model was fitted to estimate total productivity, using empirical data for mushroom production.The model was also uploaded to the SiManFor platform to simulate the effects of different thinning intensities and climatic scenarios on mushroom production.
The mushroom productivity model indicated that mushroom production increased as the average tree basal area (SMI) of the stand increased.Many mushroom species establish mycorrhizal symbioses with trees; thus, modifications in stand structure will likely affect mushroom production [8,9,13,25].Stand variables have been included as regressors in previous fitted models.Previous researchers showed that the production of Lactarius deliciosus was influenced by dominant height and stand basal area [24].In a separate study, they also found that the latter variable heavily influenced the yield of Boletus edulis Bull.Fr.As well as this, [13] demonstrated the importance of basal area in mushroom production.They found maximum mushroom yields associated with basal areas of approximately 20 m 2 ha −1 in planted Pinus sylvestris L. stands in north-eastern Spain.Other authors [27] reported optimal basal area ranging between 10-15 m 2 ha −1 in Pinus halepensis Mill.(growing under rather arid conditions) and between 35-40 m 2 ha −1 in Pinus pinaster stands growing under good

Discussion
In this study, mushroom production in Pinus pinaster stands in northern Spain was analyzed.A predictive model was fitted to estimate total productivity, using empirical data for mushroom production.The model was also uploaded to the SiManFor platform to simulate the effects of different thinning intensities and climatic scenarios on mushroom production.
The mushroom productivity model indicated that mushroom production increased as the average tree basal area (SMI) of the stand increased.Many mushroom species establish mycorrhizal symbioses with trees; thus, modifications in stand structure will likely affect mushroom production [8,9,13,25].Stand variables have been included as regressors in previous fitted models.Previous researchers showed that the production of Lactarius deliciosus was influenced by dominant height and stand basal area [24].In a separate study, they also found that the latter variable heavily influenced the yield of Boletus edulis Bull.Fr.As well as this, [13] demonstrated the importance of basal area in mushroom production.They found maximum mushroom yields associated with basal areas of approximately 20 m 2 ha −1 in planted Pinus sylvestris L. stands in north-eastern Spain.Other authors [27] reported optimal basal area ranging between 10-15 m 2 ha −1 in Pinus halepensis Mill.(growing under rather arid conditions) and between 35-40 m 2 ha −1 in Pinus pinaster stands growing under good conditions.Understanding the influence of basal area in mushroom production is very important and can be taken into consideration in management schedules for planning thinning operations.Our results would be of great interest to foresters seeking to direct forest management towards mushroom production, as it introduces the SMI parameter, which can inform strategies involving thinning from above or from below.Further data series and fitting process could advance in the response of the different mushroom groups (mycorrhizal/saprotrophic and nonedible/edible/marketed) in order to know their response, and moreover, the pattern of mushroom species with a clear social and/or economic value.
Our results also linked total productivity to autumn weather characteristics, as we hypothesized.The 12-year field-data series showed a strong correlation of autumn rainfall and temperature with mushroom production.In the literature, weather is mentioned as the most critical factor affecting the emergence of sporocarps [21,61,62].Studies based on long-term mushroom inventories, such as those conducted in Switzerland [63], Austria [64], Canada [65] and Spain [19], also indicate a strong positive correlation between increased rainfall and greater yields for all mushroom species in general.Moreover, rainfall and temperature were significant predictors in previous fitted models [14,21,22,24].Total production positively correlated to mean annual precipitation, and correlation improved when precipitation during the fruiting season and the month prior was considered [21].September precipitation showed significant positive relationship with the yield of ectomycorrhizal species in Pinus pinaster in Northeastern Spain [14].On the other hand, while [63] reported a correlation between productivity and precipitation from June to October in a mixed forest in Switzerland, [25] found a correlation between mushroom production and precipitation from August to November in pine forests in north-eastern Spain.In a research focused on Tricholoma portentusum emergence, the authors [66] showed that this species required a minimum initial amount of rainfall, followed by appropriate temperatures.Moreover, accumulated precipitation from late summer was found to be a driving factor in the occurrence of edible, marketed and L. deliciosus mushrooms in Pinus pinaster forests in Central Spain [28].Regarding temperature, [67] found a negative influence between Boletus production and the maximum temperature two weeks before the sampling week.In contrast, [17] found that high temperatures limited mushroom yield at the beginning of the fruiting season, but tended to enhance it towards the end in Pinus pinaster in Northeastern Spain due to an elongation of the fruiting season.Similar to our results, [19] found a slight correlation with November temperature, suggesting that high temperatures in mid-to-late-autumn can contribute to expansion of the fruiting period.In the current context of climate change, our model can support management decision-making, since changes in climatic conditions (increasingly irregular precipitation accompanied by extreme rainfall events [68]) are expected in the Mediterranean area.
In a more advanced approach linked to our predictive fitted model, the simulations described in this paper provided additional information about mushroom emergence.Our findings showed that moderate and high thinning activities increased total mushroom productivity.The greatest increment in mushroom production occurred at the high thinning intensity.Although in this paper the direct thinning effects on mushroom productivity were not considered, our results have shown that in G10 silvicultural scenarios, the reduction in mushroom productivity did not balance by opening the canopy, a pattern that could explain the positive increment in moderate and high thinning intensity.Several studies have evaluated how thinning intensity affects total mushroom productivity.While [69] did not find clear positive effects, other authors [14,62,[70][71][72][73] observed that mushroom productivity increased with thinning.A literature review carried out by previous researchers [73] indicated that the influence of forest management practices depended on the specific fungal ecological needs, reproductive strategies, forest type and management regime.Moderate thinning intensity could immediately boost mushroom productivity as it was shown in Maritime pine forests in northeast Spain [62] in Italy [74].Our simulations corroborated these findings, indicating that thinning operations adjusted to encourage maximum sustained yield within the framework of continuous cover forestry may also enhance the supply of wild mushrooms.Further studies could research about the impact of thinning in a bigger range of stand basal area.In this paper, the results have not shown a saturation effect or a parabolic curve, which would reflect a more realistic pattern, because of the values of the studied basal area stands.Therefore, future research could be carried out in order to know the response in more mature stands.
In addition, our simulations showed that total mushroom productivity increased in warmer weather conditions.Similar results were found in a simulation for 2016-2100 due to the elongation of the fruiting season [17].Current predictions indicate that temperatures will increase between 1.4 • C and 5.1 • C by 2055 [75].However, our findings showed that precipitation was found to have a greater effect on total mushroom productivity than temperature.Total annual precipitation projections show a tendency towards less precipitation (10%-15% reduction, principally during the summer months) combined with more extreme rainfall events [76].Simulations of drier conditions pointed to a critical situation for mushroom productivity in this type of forest ecosystem.In the three silvicultural scenarios simulated in this study, total mushroom productivity dropped to near zero in the driest conditions and the effect of silviculture was minimal due to the dominance of climatic effects.
Among silvicultural activities, forest thinning has been suggested as a forest management option for mitigating the impacts of climate change on Mediterranean forests.Thinning has the potential to increase water availability and water use efficiency in the remaining trees.This would also alter soil microclimatic conditions.
In this study, a combined analysis of silvicultural and climatic scenarios showed that under favorable climatic conditions, thinning increased total mushroom productivity.This supports our call for further research to evaluate the effect of different forest management strategies in this ecosystem services.

Conclusions
Forestry is progressing toward more complex and multifunctional silvicultural schemes that seek to optimize tree and mushroom production in the current context of climate change.Long-term studies provide valuable data for fitting relevant models that can inform forest management decisions.The estimated model presented here, based on our long-term field-data series, showed the influence of SMI, precipitation and temperature on mushroom production.The fitted model and simulations enable us to predict and quantify the consequences of management and climate scenarios on potential mushroom productivity in Maritime pine stands.
Results of the simulation study showed that wetter-and hotter-than-average scenarios increased total mushroom productivity.The best mushroom productivity occurred in scenarios with higher temperatures and precipitation.The simulation results also showed that scenarios with moderate and high thinning intensities had a positive impact on mushroom productivity.The impact was highly positive in wetter scenarios, whereas only slightly positive and negative responses were found in hotter and drier scenarios, respectively.Thus, the simulations indicated that precipitation had greater influence than temperature on total mushroom productivity.These results are very important in light of climate forecasts that predict drier situations in the Mediterranean basin in the coming decades.Understanding how silvicultural and climate scenarios affect mushroom productivity can enable us to design adaptive forest management strategies in the current context of climate change.

Figure 1 .
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.

Author Contributions:
Conceptualization, J.A.O.d.R.; methodology, J.A.O.d.R., V.P. and F.B.; software, C.O.; formal analysis, C.H. and I.B.; data curation, P.M.-P.; writing, C.H.; review and editing, C.H. and P.M.-P.; supervision, P.M.-P.and J.O. Funding: This study was made possible through the Simwood Project [Sustainable Innovative Mobilisation of Wood, co-funded by the European Union Seventh Framework Programme FP7 under Grant Agreement n 613762]; the Formixing Project [Mixed Forest complexity and sustainability: dynamic, silviculture and adaptive management tools, funded by the Spanish Ministry of Economy and Competitiveness under Grant Agreement n 345151982-51982-45-514] and the Torres Quevedo programme [grant PTQ-12-05409, co-funded by the Spanish Ministry of Economy and Competitiveness and the European Social Fund].

Table 1 .
Stand, site and soil characteristics of the plots.

Table 3 .
Results indicating significant correlation levels between stand, site, soil and climatic variables and total mushroom productivity.

Table 4 .
Solution and Type 3 tests for fixed effects in the mixed model.

Table 5 .
Variance parameters of the linear mixed model.

Table 6 .
Total sporocarp productivity before thinning by climatic scenarios.