IntroductionTop
Sustainable management of the forests require accurate tools for forecasting the different ecosystems services provided under different management schedules, and different environmental and economical scenarios. While growth and timber production models have experienced a huge advance in the last years for the main part of the Mediterranean forest species (see e.g. Bravo et al., 2011), the development of models focusing on non-wood forest products has not been so intense (Calama et al., 2010). Nevertheless, important research has been recently carried out on some non-wood products as mushrooms (Bonet et al., 2010), cork (Sánchez-González et al., 2007; Vázquez-Piqué & Pereira, 2008) and especially on pine nuts.
Pine nuts from stone pine (Pinus pinea L.) constitute the most important edible seed and one of the most important non-wood products obtained from the Mediterranean forests. Pine nuts are a high valuable product, with current retail prices exceeding 100 €/kg, and with an expanding market out of the region. In many forest areas in the Mediterranean basin, pine nuts represent the most important income for forest owners, even surpassing the profits obtained from timber and fuelwood (Finat & Gordo, 2014). In the last decades, a huge effort has been carried out in Spain, Turkey and Portugal aiming to the domestication of the species, by means of intensive plantation, reproductive propagation by grafting of selected clones for cone production (Mutke, 2013) and application of irrigation and fertilization techniques (Calama et al., 2007). However, the main part of the cones collected and marketed are still cropped from trees growing on natural forests or mature afforestations (Mutke et al., 2012).
Several tools are nowadays available aiming to carry out predictions of cone and pine nut production at different spatial and temporal scales. While the first efforts in modelling cone production were oriented to develop stand level models (Castellani, 1989; García-Güemes, 1999), tree level models have performed better (Cañadas, 2000; Calama et al., 2008), since they allow taking into account both between and within-tree variability in production. However, all these models were pure spatial, allowing to predict average annual cone production at different tree, stand and forest scales, but neglecting the strong masting habit of the species. Mutke et al. (2005a) identified the most important climate factors ruling interannual variability in cone production, establishing the main basis for the annual tree level model for cone production in the Spanish Northern Plateau, by Calama et al. (2011). This spatiotemporal model, constructed using data from 750 trees with cone series covering ten years (1996-2005), constitutes up to the present the only available tool for predicting annual cone production at tree and forest scales in Spain, and has been implemented within the process based model PICUS (Pardos et al., 2015).
Despite the wide use of this model, severe limitations are identified. The first one is to the narrow spatial range of applicability of the model, limited to the Northern Plateau region in inner Spain (60,000ha). Direct application of the model to other relevant stone pine regions (such as Andalusia or the Central Range) has resulted in biased predictions, thus specific models for these regions are highly required.
A second main limitation is that the model was constructed using a limited series of data covering ten years (1996-2005), and then validated over an extended series of three additional years (2006-2008). Besides this, due to the absence of an independent series for validation, available data set was split out into a fitting and a validating data set, which resulted in a significant reduction in the size of the fitting data set. The accuracy of the model was then successfully evaluated over different years out of the initial temporal series. However, the constructed model has failed to predict the observed yield decline in the three most recent crops (2012-2013, 2013-2014 and 2014-2015) in the inner part of Spain, a process common to other regions (Mutke et al., 2014). Two main hypotheses have been postulated to explain such decline: a constraining climate factor not considered in the current formulation of the model; and the effect of a new biotic agent, the exotic seed bug Leptoglossus occidentalis (Heid.) (Bracalini et al., 2013; Lesieur et al., 2014). Seed predation by this insect has been described as a cause of the abortion of unripe conelets, as well as for provoking a significant reduction of the seed yields in mature cones (Mutke et al., 2014; Calama et al., 2015).
Whatever the main cause is, our current model tend to overestimate the production in recent years, thus the construction of a new model using an extended series of cone production covering up to the most recent crops will give insight on the potential cause for this decrease on the production.
Another factor that might be taken into account for modeling the masting pattern of stone pine, additionally to the predominance of drought-prone climatic covariates (Calama et al., 2011), is the negative autocorrelation with lag 3 observed in the 40-year cone-yield time series at regional scale (Mutke et al., 2005a). This effect has been interpreted as reduction in the number of female flowers induced by the high physiological cost during the ripening of heavy crops (Mutke et al., 2005b). This inhibitory effect of depletion has been largely described in other species (e.g. Sork, 1993; Sala et al., 2012). The original model at tree and stand scale did not consider it, because it would have reduced the ten years series to only seven.
In relation to the original formulation and structure of the model, it included as explanatory covariates stand age and site index, which prevents the sound application of the model in multi-aged stands, where these attributes are not interesting. Apart from this, while common management inventories on even-aged stands may include an estimate of stand age (which can be then applied to compute site index), regional or national scale inventories, as National Forest Inventory, do not include such information, thus limiting the potential use of the model in predicting cone production at these extended spatial scales. Additionally, the annual model does not produce period interval predictions of cone production. Due to this, it is difficult to couple the model with either the existing distance independent tree-level model for the species PINEA2 or the associated stand-level simulator, which makes out predictions at five years steps.
In consequence, the main aim of the paper is to present enhanced and new tools for forecasting cone production in Pinus pinea forests in inner Spain aiming to overcome the previously exposed limitations. To tackle with this we propose to attain three specific objectives:
- To build up an enhanced version of the model developed by Calama et al. (2011) for annual tree cone production for the Northern Plateau, using new available data series. The constructed model will: (i) be compatible with common input data from different sources (NFI, forest management inventories), (ii) be flexible enough to consider different changing climate scenarios, (iii) allow compatibility among annual and periodic estimates of cone production; and (iv) take into account resource depletion effects.
- To build up a new model for annual tree cone production with validity in the forests of the Central Range, sharing similar characteristics and structure as those previously described.
- To present an independent validation on the constructed models based on an approach for upscaling predictions from tree to forest level.
Our main hypothesis are that: (i) the observed recent decay in the cone production can be explained by including climate factors not adequately considered in the previous versions; (ii) the consideration in the model of the delayed inhibitory effect of large cone crops will improve the annual estimates of cone productions; and (iii) the processes controlling cone production in both regions are similar, thus main regional differences are related with the timing (phenology) intensity of the factors rather than on genetic differences.
MaterialTop
Studied regions
Pinus pinea forests in the inner part of Spain occupy two main regions, defining two different genetic provenances (Prada et al., 1997): Northern Plateau (from now on NP) and Central Range (from now on CR), with different environmental conditions and management tradition and practices (Figure 1).
|
|
The NP region is a plain defined by the basin of Duero river, owing two main differentiated units: sandy areas, at an average altitude of 700 – 750 m and limestone plains, at an altitude over 800-850 m. Within the NP Pinus pinea covers more than 60,000 ha, mainly in the province of Valladolid. Lithological differences have resulted in different soil types. Sandy soils show a very high content of sand (> 90%) and very low water holding capacity (WHC<100 mm), while soils in limestone area , with percentage of clays and limes over 40-50% , reach values of WHC > 250 mm. With respect to climate conditions, NP is a quite homogeneous territory, characterized by a Mediterranean continental climate, with very low precipitations (average annual rainfall: 450 mm), summer drought (<50 mm between July-September) and cold temperatures (average annual temperature 10.5–13.4 ºC, minimum absolute temperature below –10 ºC). These forests have been managed since the end of 19th century, aiming to guarantee soil protection and optimize both cone and timber production, mainly resulting in pure, even-aged stands. However, more complex forests, showing uneven aged structures and mixtures with other conifers and broadleaves, are commonly found in limestone areas.
Within the CR region stone pine occupies the rocky granitic slopes, which define the basins of the Alberche and Tiétar rivers, between the provinces of Ávila and Madrid. These valleys conform a mountainous area, with an orography of hills and small valleys, with steep slopes defining an abrupt landscape. In this territory, stone pine grows at an altitudinal strip between 600 and 1,000 m occupying more than 30,000 ha. Soils are in general poor, mainly composed by sands (70%-80%), though the main differentiating characteristic is the high presence of granite rocks in surface layers, accounting for 50% of the total soil volume. Concerning climate, CR is a moist Mediterranean region, with average annual rainfall exceeding 700 mm but summer precipitation below 50 mm, and warmer temperatures (average temperature: 14.0 ºC with absence of extreme freezing events). The complex orography of these stands as well as the recurrent frequency of forest fires have oriented their main production towards pine nuts and pastures, resulting in mixed, open forests with multi-aged structure, where overmature great cone producers are kept standing and no thinnings are commonly applied (Montero et al., 2003).
Fitting data set: INIA permanent plots
In 1996, INIA-CIFOR, in cooperation with the forest services of Valladolid, Ávila and Madrid, installed a net of permanent plots in pure even-aged stands of Pinus pinea within the two studied regions. The net included 141 plots in NP and 52 plots in CR. The plots are circular in shape, with variable radius and include a fixed number of trees: 20 in NP, 10 or 20 in CR. Plots were selected aiming to cover the whole range of site conditions, stand stocking and age identified within each region, searching for a uniform spatial distribution, and were always located in public forests.
At plot installation, diameter at breast height, total height, crown diameter, height to crown base and tree coordinates were measured in all the trees within the plot. In a subsample of trees, total age was determined by extracting cores at stump height with a Pressler increment borer. Plots were inventoried in 2001, 2008 and 2015. Main differences among the studied regions in stand and tree level attributes are related with the overmaturity and lower stocking of the stands in CR (Table 1).
Table 1. Mean dendrometric characteristics of the fitting data set plot and trees in each region (NP: Northern Plateau, CR: Central Range) and between-regions differences
In the five trees closer to the plot center, mature cones were collected every year by specialized climbers. Cones from each tree were classified as healthy or damaged by pests (larvaes from moth Dyorictria mendacella (Stgr.) or the beetle Pissodes validirostris (Gyll.), and cones from each category were counted and weighted separately. Between 1996 and 2005 (10 crops years) cones were collected in all the plots within the regions, comprising a total of 750 trees in the NP and 260 trees in CR. From 2006 to 2008 cones were collected in a subsample of 30 plots (150 trees) in NP and 15 plots (75 trees) in CR. Finally, in 2009, 2010, 2013 and 2014 cone collection was maintained in the 30 plots of the NP. Due to budget restrictions, no cones were collected in 2011 and 2012. According to this design, the theoretical number of tree x year observations of cone production should result in 8550 for the Northern Plateau and 2825 for the Central Range. However, due to illegal cone collection, and/or tree natural or catastrophic mortality, total number of tree x year observations resulted in 8074 in NP and 2671 in CR.
Validation data set: provincial series
As an independent data set for validating the model we used the data from the annual series of cone production at forest scale managed by the forest services. Annual cone harvesting rights in public forests are sold every year in public auctions, and at the end of the campaign, enterprises should give a report including total weight of cones collected within the forest. These post-crop data series at forest scale extend from 1962 up to the present in NP (Valladolid province) while in the CR there are available estimates of total cone production covering from 2004 to 2012 campaigns.
MethodsTop
Response variable: annual weight of healthy cones at tree level
As in the previous spatiotemporal model, the proposed response variable was the weight (expressed in kg) of sound cones annually collected in each tree wc, which is a better indicator than cone number of allocated resources to the reproductive effort, and takes into account the processes of female floral induction and pollination as well as subsequent cone growth, maturation and effect of common endemic cone pests. As stated in Calama et al. (2008, 2011), this variable shows severe departures from normality assumptions, due to asymmetry in the distribution, skewness, and excess of zeros. With respect to the excess of zeroes, for the analyzed series average percentage of null observations (Table 2) reach 54.2% in NP (ranging between 18.8% in 2000 and 85.4% in 2014) and 27.89 % in CR (5.5% in 2000 and 44.3% in 2008).
Table 2. Mean anual value of cone production at tree level for the studied series (1996-2014 in Northern Plateau; 1996 – 2008 in Central Range)
Zero –inflated models
The observed excess of zeros in the data set largely prevents from using classical statistical approaches. Zero-inflated models (Lambert, 1992) conform an interesting approach for dealing with this kind of data. These models account for the zero observations and the positive outcomes separately. Distributional particularities are dealt with by combining a distribution that explains the binary nature of absence or occurrence of zeros, with a second distribution that models the response variable, conditioned by its occurrence. In this case, the phenomenon under examination is considered the result of two different processes: the first, with a binary outcome, is the occurrence of the event (in our case, fruiting); and the second, determining the intensity of the event (in our case, total weight of the tree cones), conditioned by the occurrence of the event.
The binary part is commonly modelled using a logistic function assuming a binomial distribution. Given the continuous nature of the abundance part, we propose to use a log-normal distribution, since it avoids zero-truncation (meaning that zero responses can only occur in the binomial part of the model) and offers greater potential for modelling highly skewed distributions (Tu, 2002). The so-called Zero-Inflated Log-Normal (ZILN) model can then be expressed in a bietapic form:
Where θ represents the probability for having a non-zero value of wc (i.e., obtaining a non null crop), and μ represents the expected value of the logarithmic transformation of wc, conditional to a non-null crop; xand z represent vectors of known possible explanatory covariates for modelling θ and μ, respectively; while α and β are vectors of unknown but estimable parameters. For more information on ZILN distributions and modelling see Calama et al. (2012).
Explanatory covariates and model construction
Given the huge differences observed in individual tree cone production between NP and CR, the ecological differences observed among the two analyzed regions, and the different length of the analyzed series, we proposed to fit independent ZILN models for each region. However, we attempted to maintain a similar structure of the model in both regions. To do that, we evaluated the same covariates acting at different spatial and temporal scales as potential predictors of cone production.
Spatial attributes (tree and plot scale)
- Tree size: diameter at breast height, section at breast height, total height.
- Tree level competition: basal area of trees larger than subject tree, ratio dbh / mean squared diameter, ratio tree section / basal area.
- Stand attributes: stocking density (stems/ha), basal area, mean squared diameter, Reineke’s stand density index , dominant height, crown cover factor.
Logarithmic, root and inverse transformation of these variables were also evaluated as potential predictors. Stand age and site index, which were predictors in the previous version of the model, were not more considered as potential predictors.
Gordo (2004) proposed a stratification for the territory of the NP into Natural Units (1,400 – 7,200 ha of stone pine forest each), based on climate, soil, lithology, orography and management characteristic which was included as explanatory in Calama et al. (2008, 2011). In the present study, we proposed a similar stratification for the CR region into four regions (2,000 – 4,700 ha), evaluating the entrance of these stratifications as categorical dummy variables (Table 3, Figure 1).
Table 3. Natural Stratification of the territory in the two analyzed regions
Temporal attributes
- Rainfall: monthly precipitation, starting from the month of January four years before the year of cone maturation and ending up in October of the year of maturation. Apart from single month precipitation, cumulative precipitations in periods covering up to twelve months were also evaluated.
- Temperature: number of freezing days per month.
- Lagged response: the possible inhibitory delayed effect three years after a bumper crop observed, was taken into account by using data series at regional scale. Average regional value for cone production per ha one, two, three and four years before the current crop were evaluated as potential predictors. This data were obtained from the provincial series of annual cone production.
Rainfall and temperature monthly series were obtained from the most complete meteorological station located within each region: Valladolid meteorological station (41° 39’ 8’’ N - 4° 43’ 24’’ W, 690 m a.s.l.) in NP, and El Tiemblo station (40°24’56” N - 4°30’02” W, 685 m a.s.l.) for CR. The use of a single reference station per region was preferred, because masting pattern had been found to be regionally synchronized rather than locally (Mutke et al., 2005a; Calama et al., 2011).
To ensure the compatibility among annual and periodic predictions, climate variables and the lagged response were standardized by substracting to the observed annual value the mean value for the series and dividing by the standard deviation. By doing this, periodic predictions over an average year can be carried out by fixing a value of zero for all the temporal attributes. Additionally, the standardization allows to give insight on the relative importance of each temporal predictor over the response variable.
Covariate selection and model construction were carried out following a sequential procedure. In an initial phase, two separate models for the occurrence of the event (fruiting /no fruiting) and the abundance of the event (total weight of cones, in logarithmic scale, only considering nonzero observations) were separately fitted. In this preliminary step, models were fitted as generalized linear mixed models, including random effects at plot, tree, year and plot x year effects. Starting from a basic model that only includes the intercept, the empirical Bayesian estimates for random effects at different levels were correlated with the evaluated covariates acting at those scales, i.e., year random effects vs temporal covariates, plot and tree random effects vs stand and tree level attributes. The entrance of the covariates largely explaining observed variability at each level was then evaluated for both the occurrence/abundance submodels. The process was sequentially repeated until no more significant covariates entered the submodels. This preliminary set of explanatory covariates was then used as predictors to fit a simultaneous zero-inflated lognormal model. The final model only includes random effects at plot level in both the occurrence and the abundance submodels.
Goodness of fit
Model accuracy was first evaluated over the fitting data set. To do that, marginal predictions (assuming unknown random effects
equaling zero) over the response variable, expressed in original weight of cones scale, were computed. Antilogarithmic transformation
of the response variable was carried out by applying exponential transformation and then multiplying by the factor,
represent the variance associated with random plot effect and residual term, respectively.
Thereafter we computed mean error (E), level of significance of the error (p-value), root mean square error (RMSE) and modelling efficiency (MEF) at different scales:
- Tree x year: comparing observed and predicted values of the response variable wc.
- Tree: comparing the observed and predicted total production at tree level for the whole period.
- Plot x year: comparing the observed and predicted annual production at plot level.
- Plot: comparing the observed and predicted total production at plot level for the whole period.
- Year: comparing observed and predicted annual average production at tree level.
Independent validation
The constructed models for each region were validated against the available data from the series of cone production at forest scale previously described. The model for the Northern Plateau was validated over the annual series 2004 – 2014 from five forests (numbers 35, 40, 44, 39 and 79, covering more than 5500 ha, 10% of the total stone pine area within the region). We validated the model for the Central Range over the annual series 2004-2012 from forest 75, the most important in the region, covering 1000 ha. Upscaling the predictions from tree level to forest level was carried out by using data from the forest management inventory. Sampling plots are circular, with a fixed radius of 15 m, and are systematically located in a grid of 200 m x 200 m within each forest. Collected attributes include tree dbh for all the trees within the plot and tree height and age in a subsample of trees. The constructed model was then applied to each tree within the plot, and plot annual production was computed as the sum of the predicted productions for each tree. Forest annual production was then upscaled by using common sampling estimation procedures (Mandallaz, 2007). Annual predictions at forest scale were then contrasted with annual post-crop estimates by means of E, RMSE and MEF. In the case of forest 75 (CR) we took into account the observation by the forest service indicating that uniquely in 25% of the total area of the forest cones are collected with commercial purposes due to steep slopes preventing mechanization, isolated areas and distance from forest roads.
ResultsTop
Annual cone production in the studied regions
There is a significant difference in cone production at tree level between both regions (Table 2). CR almost doubles the average production of NP for the studied series in both the weight of cones (4.575 kg/tree in CR vs 2.472 kg/trees in NP, p-value <0.0001) and the number of healthy cones per tree (15.064 cones/tree in CR vs 8.952 cones/tree in NP, p-value <0.0001). This difference is also patent in the percentage of “zero” observations (trees not bearing a single cone in a given year), reaching 54.21% in NP (ranging between 18.88% in 2000 and 85.42% in 2014) and 27.89 % in CR (5.49% in 2000 - 44.28% in 2008).
For the common thirteen-years series (1996-2008), no significant differences between regions are found in 2005, 2006 and 2008 (p-value according to non-parametric Kruskall-Wallis test 0.1442, 0.0563, 0.1660 respectively). Uniquely in 1996 and 1997 average crop in NP was higher than in CR. Constant differences in these values point out to huge and permanent spatial differences in the process of fruiting among regions.
Large interannual variability in cone production is found in the two studied regions, with no significant pattern of synchrony between them (Tau´b kendall: 0.0512, p-value: 0.8072; Figure 2). This could indicate (i) that different climate attributes or phenological timing are acting at each region, or (ii) that even though the same climate factors ruled annual variability in cone production in the two analyzed regions, the intensity of these factors show interregional differences. However, it is patent that some bumper years are common in both regions (e.g. 1996, 2000), which evidences common factor triggering huge cone crops.
|
|
Model fitting
After an independent sequential procedure for model fitting, a similar set of covariates was selected for each region, revealing large interregional similarities in the factors ruling cone production (Table 4). Concerning spatial variability, main predictors explaining cone production at stand and tree level are common in NP and CR: the inverse of the mean squared diameter (1/dg), Reineke’s stand density index (SDI) and the ratio between tree diameter (dbh) and dg (dbh/dg); while tree dbh only enters the model in the abundance part for NP. Parameter sign and magnitudes are similar in both regions. The proposed natural stratification of the territory is also highly significant in both regions, pointing out to the importance of this type of stratification to synthetize the effect soil, lithology, orography and other spatially explicit covariates. Noteworthy, once the proposed set of covariates enter the models, the entrance of either stand age or SI did not improve the models.
Table 4 Parameter estimates for the ZILN models for cone production for Northern Plateau and Central Range
In relation with temporally explicit covariates, results show that in both regions the most important factor ruling cone production is the amount of rainfall from spring to fall during the year previous to flowering, together with the rainfall during the months just before and after flowering, as well as the rainfall during the last year of maturation. While general similitudes are observed, the exact season do not match exactly between the two regions, probably due to phenological shifts along termal gradients (Mutke et al., 2003). However, two main differences between the regions are found: (i) the observed negative impact of the number of days with severe freezing temperatures (< -5ºC) during the first winter after flowering in NP, and (ii) the positive effect of the rainfall from September to December during the 2nd year of maturation in CR. Finally, a negative effect of the cumulative sum of the cone crops two and three years before cone maturation is found in both regions. Taking into account that all the temporally explicit covariates are standardized over the mean and standard deviation values for the climate series 1980 – 2010, further use of the model will require these values (Table 5).
Table 5. Mean annual and standard deviation values of the temporally explicit covariates entering the models (required for standardization)
Model goodness of fit
Accuracy of the fitted models over fitting data set was similar for both regions. Predictions were significantly unbiased at all the analyzed spatial and temporal scales, with relative error at scale of prediction always below 6% (Table 6). The fitted models explained over 66% – 69% of the total variability in cumulative cone production during the studied period between plots, and 50% of the variability between trees. Concerning spatiotemporal dynamics, 50% of the total annual variability between plots and 30%-40% of the residual variability at tree x year scale is explained by the models. Finally, modelling efficiency reach values close to 90% in explaining interannual variability in cone production at region scale, clearly mimicking the observed masting pattern (Figure 3).
Table 6. Goodness of fit statistics of the marginal predictions (in real scale) of the fitted models over fitting data set
|
|
We compared the constructed model with a new refit of the previous model by Calama et al., (2011) to the currently available data set for the Northern Plateau, identifying that the new proposed set of covariates explain better the observed variability in the response variable (AIC: 19279 vs 19593). Apart from this, the previous set of covariates failed to explain the observed decay in the crops following 2012-2013 campaign, a process which is now mimicked by the new model (Figure 3).
Independent Validation
Upscaling the predictions by the model to the forest scale led to unbiased results for the annual cone production of the six analyzed forests over the studied series 2004-2014 (Table 7). Mean error ranges from –14 to 45 kg.ha.year–1. Modelling efficiency at forest scale range from 0.61 to 0.90, except for the forest 40, where due to the exceptional bias in 2005-2006 (predicted crop at forest scale over 205000 kg, while collected crop uniquely reached 66000 kg), MEF is reduced down to 0.23. Root mean square error for the analyzed forests varies from 20 to 95 kg.ha.year–1. The model attains to mimic the observed pattern of interannual variability at forest scale, mainly in differentiating among bumper and non-bumper crops (Figure 4), and especially predicting adequately the reduced crops observed after the 2010-2011 exceptional bumper year.
Table 7. Prediction accuracy for the six forests included in the independent set of evaluation
|
|
DiscussionTop
The constructed model for the Northern Plateau means an improvement over the previous model (Calama et al., 2011), making use of the extended data series of cone production, now covering up to 17 years in some plots. From the original model, some variables difficult to measure or time and money consuming, as the stand age, or necessarily estimated from other models, as the site index, were successfully removed by using alternative approaches, based on stand density, mean squared diameter, and the stratification in Natural Units. In the case of site index this result indicates that a sound ecologically based stratification of the territory can be as useful as site index in explaining differences in quality, as previously stated in other works by the species (Calama et al., 2008). Stand age has been substituted by the inverse of the mean squared diameter, a predictor commonly measured in forest inventories that synthetizes the effect of both stand maturity and stand stocking. Stand density, defined as the number of stems per ha, was now integrated in the Reineke’s stand density index, which among other advantages considers both stocking and tree size, being much more robust to immediate changes just after a low thinning. Finally, as in the previous version of the model, both tree size (dbh) and competitive status (dbh/dg) entered the model.
With respect to the predictors explaining interannual variability the main difference with the previous model is the entrance of the model of the delayed effect occurring two and three years after a bumper crop. This inhibition of flower induction by resource depletion due to ripening crops has been previously described in stone pine (Mutke et al., 2005a) and is related with the proposed hypothesis of resource depletion (Sork, 1993; sagi et al., 1997I). Our result indicate that after a bumper year reduced crops can be expected with a delay of two and three years, pointing that the allocation of resources to enlarge cones in a bumper crop results in simultaneous inhibition of both bud and female floral induction (Mutke et al., 2005b, Crone et al., 2009, Sala et al., 2012).
The evaluation of a more complete set of climate covariates and an enlarged series of years resulted in the entrance within the model of the cumulative sum of the rain between May and November the year before flowering, instead of the separated May-June and October-November precipitations that entered in the previous version. This indicates a continuous influence of water availability during the whole vegetative period the year before flowering, which can directly enhance flowering by affecting bud development and differentiation (which was the main previous theory), but also indirectly by means of an increment on needle length, leaf area index, annual net assimilation and accumulated carbohidrates (Mutke et al., 2003). This finding agrees with the theories postulating complex trade-offs between flowering-fruiting process and other key life-history variables (Knops et al., 2007). Another finding of the present work is the influence of rainfall events just before flowering, pointing that the process of bud differentiation into vegetative or reproductive shoots might be delayed with respect to previous hypothesis. Finally, rainfall occurring the summer immediately after pollination, during the last months before maturation as well as the occurrence of extreme freezing events influence conelet survival and cone enlargement during the last period of maturation (Calama et al., 2011).
In the case of the Central Range, the constructed model conforms the first approach to modelling cone production in the region based on climate, stand, tree and site attributes. Previous modelling efforts, as the one presented by Calama & Montero (2007), mainly focused on identifying and quantifying the different spatiotemporal sources of total variability in cone production at different spatiotemporal scales, not aiming to obtain predictions at tree of forest scale. Predictors entering the model for CR are mainly common with those for NP. In the case of spatial variability the unique difference is that dbh is missing as a predictor in CR, probably due to the wider range of diameters and the particular irregular structure of some of the stands observed in the region (Montero et al., 2003). A stratification of the territory based on soil, climate and altitude attributes, similar to that of NP, successfully substituted site index. In the temporal explicit variables differences are found in the timing of the periods of seasonal rainfall influence (Figure 5), indicating an anticipated phenology and an extended period of vegetative growth in CR with respect to NP due to warmer spring and raindalss (Montero et al., 2008). Other main differences are related with the null influence of extreme freezing events in CR, where are quite uncommon.
|
|
As expected, our results confirm that the factors affecting cone production in both regions are similar, giving support to the hypothesis of correlated environmental variables rulling spatial synchrony at large scales, as stated by other authors (Koenig & Knops, 2013). Thus the main differences detected in cone production between the two regions are more related with the observed differences in the timing and intensity of these influential factors rather than with regional differences in the response to these factors. In this sense it is noteworthy to indicate that larger cone productions in Central Range are mainly related with the larger amount of annual rainfall and the absence of extreme freezing events, together with intrinsic stand factors, as is the larger ageing and lower stocking of the stands (Table 1). This overmaturation of the stands in CR, linked with the lack of application of a strict silviculture can result, in the next decades, in a significant reduction in cone crops within the region associated with the already observed processes of natural mortality and decay of overaged trees (Spathelf et al., 2014).
The predictions by the models in both region mimicked the observed patterns of spatial and interannual variability within the fitting data set, accurately differentiating good and bad crop in both regions. Model application lead to unbiased estimates of cone production in all the different spatiotemporal scales analyzed, resulting in goodness of fit statistics similar in both regions. In the case of the NP the new fitted model is able to explain the recent decay in cone production during the campaigns extending from 2012 to 2015. This result supports the hypothesis of climate causes of the reduction, aggravated by the inhibitory effect of the exceptional 2010 cone crop, against the hypothesis focusing on biotic agents like Leptoglossus occidentalis. Expected decrease in cone production associated with new scenarios of rainfall reduction has been also predicted by means of complex process based models (Pardos et al., 2015). However, the effect of this exotic bug over kernel yield per cone weight has already been described in the regions (Mutke et al., 2014, Calama et al., 2015), and the potential effect on immature conelets should be evaluated by contrasting observed and predicted cone crops in years where larger cone production is expected from climatic conditions. Moreover, specific ongoing experiments studying the evolution of female conelets through the ripening cycle under controlled levels of the insect will give insight on this topic (Strong, 2015).
Upscaling the predictions of the model from the tree to the forest scale permits to obtain unbiased estimates of annual cone production for the forest management units. Application of the models using input data from forest management inventories and updated climate data from regional meteorological stations will allow obtaining reliable estimates of cone production by the end of July of the year of cone maturation, four months before the beginning of cone collection, and by the time where public auctions for cone production are carried out. Given the standardized character of the climate variables, it will be possible to make even one-year advanced predictions of cone production, defining a confidence interval for the prediction based on different scenarios of rainfall during the last months of maturation. All this advanced information is essential for a sound planning of all the activities involving the different stakeholders, from the ownerships to the cone collection enterprises, forest services and pine nut industries.
Although the fitted models requires climate predictors extending up to three-years before cone collection, given the standardized character of these predictors it is possible to make average predictions for a given period by assuming a zero value, and even consider the effect of future climate scenarios by assuming variations on the climate standardized predictors within the proposed range. In this sense, the constructed models will be easily implemented within a stand level simulator, as PINEA2 (https://sites.google.com/site/regeneracionnatural/pinea2), allowing to compare and define optimal silviculture under different scenarios. In addition, the required input dataset for the model is compatible with the data from the Spanish National Forest Inventory, permitting to obtain annual estimates of cone production at regional scale.