Forest Systems
25(3), e079, 17 pages (2016)
eISSN: 2171-9845
http://dx.doi.org/10.5424/fs/2016253-09671
Instituto Nacional de Investigación y Tecnología Agraria y Alimentaria (INIA)
RESEARCH ARTICLE
OPEN ACCESS
Enhanced tools for predicting annual stone pine (Pinus pinea L.)
cone production at tree and forest scale in Inner Spain
Rafael Calama1,2*, Javier Gordo3, Guillermo Madrigal1, Sven Mutke1,2, Mar Conde1, Gregorio Montero1,2
and Marta Pardos1,2
1
INIA-CIFOR. Ctra. A Coruña km 7.5 28040 Madrid. Spain. 2 iuFOR. Sustainable Forest Management Research Institute UVa-INIA. Spain.
3
Servicio Territorial de Medio Ambiente de Valladolid. C/ Duque de la Victoria, 8. 47001 Valladolid. Spain
Abstract
Aim of the study: To present a new spatiotemporal model for Pinus pinea L. annual cone production with validity for Spanish
Northen Plateau and Central Range regions. The new model aims to deal with detected deficiencies in previous models: temporal
shortage, overestimation of cone production on recent years, incompatibility with data from National Forest Inventory, difficulty
for upscaling and ignorance of the inhibitory process due to resource depletion.
Area of study: Spanish Northern Plateau and Central Range regions, covering an area where stone pine occupies more than
90,000 ha.
Material and methods: Fitting data set include 190 plots and more than 1000 trees were cone production has been annually collected from 1996 to 2014. Models were fitted independently for each region, by means of zero-inflated log normal techniques.
Validation of the models was carried out over the annual series of cone production at forest scale.
Results: The spatial and temporal factors influencing cone production are similar in both regions, thus the main regional differences in cone yield are related with differences in the phenological timing, the intensity of the influent factors and forest intrinsic
conditions. A significant inhibition of floral induction by resource depletion was detected and included into the model. Upscaling
the model results in accurate prediction at forest scale.
Research highlights: [1] The new model for annual cone production surpass the detected deficiencies of previous models, accurately predicting recent decay in cone production; [2] Regional differences in cone production are due to phenological and seasonal climatic differences rather than to between provenances genetic differences.
Keywords: zero-inflated models; pine nut; conelet losses; Leptoglossus occidentalis; forest upscaling.
Citation: Calama, R., Gordo, J., Madrigal, G., Mutke, S., Conde, M., Montero, G., Pardos, M. (2016). Enhanced tools for predicting annual stone pine (Pinus pinea L.) cone production at tree and forest scale in Inner Spain. Forest Systems, Volume 25, Issue
3, e079. http://dx.doi.org/10.5424/fs/2016253-09671.
Received: 18 Mar 2016. Accepted: 15 Sep 2016.
Copyright © 2016 INIA. This is an open access article distributed under the terms of the Creative Commons Attribution-Non
Commercial (by-nc) Spain 3.0 Licence, which permits unrestricted use, distribution, and reproduction in any medium, provided the
original work is properly cited.
Funding: This work was carried out within the budgetary framework of the national research projects RTA-2013-00011-C2.1
and AGL2010–15521; the European project FP7-KBBE-2012-311.919-STARTREE and the regional grant PROPINEA, from the
Deputation of Valladolid.
Competing interests: The authors have declared that no competing interests exist.
Correspondence should be addressed to Rafael Calama: rcalama@inia.es
Introduction
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 nonwood 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
2
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
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 tempo-
Forest Systems
ral 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.
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
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.
Material
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
Forest Systems
3
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
December 2016 • Volume 25 • Issue 3 • e079
4
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
Natural Units
1:700,000
Legend
N
Natural Units
1
21
2
22
3
23
4
24
5
6
7
8
9
10
1:700,000
Figure 1. Studied regions and Natural Units definition (uniquely for public forests) overlapped with the total area with presence of
Pinus pinea L.
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).
Forest Systems
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
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
5
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
Plot level
Tree level
Variable
NP
CR
p-value
Variable
NP
CR
p-value
BA
N
Age
Ho
SI
SDI
FCC
Dg
18.3
320
62
10.0
14.6
337
0.51
29.9
19.2
219
73
10.9
14.3
327
0.48
35.3
ns
0.0173
0.0422
ns
ns
ns
ns
0.0031
dbh
g
h
cw
CR
BAL
30.4
0.084
9.5
5.2
0.565
0.203
35.0
0.096
10.1
5.6
0.566
0.578
<0.0001
<0.0001
0.0137
0.0045
ns
<0.0001
Where BA: basal area (m2/ha), N: number of stems/ha, Age: years, Ho: dominant height (m), SI: site index (m), SDI: Reineke´s stand
density index, FCC: crown coverture factor, Dg: mean squared diameter, dbh: breast height diameter (cm), g: breast height section (m2),
h: tree height (m), cw: crown diameter (m), CR: crown ratio, BAL: basal area of trees larger than subject tree (m2). P-value refers to the
level of significance for the null hypothesis of interregional homogeneity.
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.
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).
Validation data set: provincial series
Zero –inflated models
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.
Methods
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,
Forest Systems
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-
December 2016 • Volume 25 • Issue 3 • e079
6
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
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)
Northern Plateau
Central Range
YEAR
wc
nc
% zeroes
n
wc
nc
% zeroes
n
5.951
1.165
2.225
3.697
13.778
4.066
3.092
5.572
4.501
1.986
4.757
5.509
1.624
18.580
5.735
5.119
12.596
40.459
14.514
11.826
18.905
15.922
8.140
19.320
15.627
4.886
14.3
56.2
42.3
24.7
5.5
26.1
26.7
18.2
22.9
40.7
14.7
28.0
44.3
245
260
260
259
255
249
247
231
231
214
75
75
70
4.575
15.063
27.8
2671
1996
1997
1998
1999
2000
2001
2002
2003
2004
2005
2006
2007
2008
2009
2010
2011
2012
2013
2014
6.433
3.451
1.462
0.701
5.489
3.214
0.984
0.913
0.150
1.873
4.910
1.200
1.526
2.661
4.988
23.325
13.770
4.233
2.813
17.871
9.256
4.304
2.725
0.595
9.454
20.255
5.193
4.942
12.658
17.909
32.1
40.1
49.7
67.7
30.5
64.4
72.0
66.9
82.7
36.8
31.0
79.3
58.7
30.6
18.9
751
750
724
722
698
711
701
680
681
672
145
150
155
111
143
0.185
0.196
0.794
0.944
77.9
85.4
136
144
AVERAGE
2.472
8.952
54.2
8074
Where wc: annual average weight of healthy cones per tree (kg), nc: annual average number of healthy cones per tree, % zeroes: percentage of trees bearing no cones in a given year, n: number of sampled trees.
Normal (ZILN) model can then be expressed in a bietapic form:
⎛ θ ⎞
= xα
logit θ = ln ⎜
⎝ 1− θ ⎟⎠
μ = E [ln(wc) | wc > 0] = zβ
E(wc | x, z) = θ μ =
exp ⎡⎣ xα ⎤⎦
1+ exp ⎡⎣ xα ⎤⎦
zβ
Where θ represents the probability for having a nonzero 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;
x and 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
Forest Systems
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
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
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).
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.
7
— 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).
Table 3. Natural Stratification of the territory in the two analyzed regions
Annual
rainfall
(mm)
Water
retention
(mm)
Nº plots
500
350
2
480
313
17
480
295
18
385
107
17
385
175
6
400
160
12
400
186
19
380
111
17
420
166
10
16.8
Alluvial
meadow soils
Tertiary
sediments
Tertiary
sediments
Tertiary
sediments
Wind
deposits
Wind
deposits
Wind
deposits
Alluvial –
Wind
Wind
deposits
Wind
deposits
Alluvial
deposits
400
158
22
925
13.6
850
150
17
2010
875
15.3
780
200
14
4713
773
15.4
Rocky
Wind
granitic
deposits
slopes
Quartz sands Alluvial –
& limes
wind
Quartz sands Alluvial –
& limes
wind
700
200
14
3980
740
12.2
Granitic boles
650
170
7
Region
Unit
Name
Area
(ha)
Altitude
(m)
SI (m)
NP
1
Torozos
1834
847
14.0
NP
2
4757
845
14.9
NP
3
2284
854
14.4
NP
4
Limestone
plain W
Limestone
plain E
Valladolid
1908
687
13.1
NP
5
Nava del Rey
1480
710
11.8
NP
6
Viana de Cega
7211
707
15.1
NP
7
Iscar
2945
745
16.7
NP
8
Medina
4044
746
12.5
Quartz sands
- Clays
Quartz sands
NP
9
1401
712
12.6
Quartz sands
NP
10
Tudela de
Duero
River terraces
1664
688
CR
21
Hoyo Pinares
2720
CR
22
Cofio Valley
CR
23
CR
24
Tiétar&North
Alberche
Valleys
South
Alberche
Valley
Geology and
texture
Soil
Origin
Limestone –
Marl
Limestone Marl
Limestone –
Marl
Quartz sands
Quartzitic
gravels
Quartz sands
Wind
deposits
Where SI is site index, as stated by Calama et al. (2003); total area is uniquely referred to public owned forests.
Forest Systems
December 2016 • Volume 25 • Issue 3 • e079
8
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
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
1⎛
factor exp ⎜
2⎝
∑σ
2
v
⎞
+ σ 2e ⎟ , where σ 2v and σ 2e 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:
Forest Systems
— 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.
Results
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
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
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 nonparametric 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. How-
9
ever, 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.
16
Central Range
Northern Plateau
14
Weight of healthy cones (kg.tree–1)
12
10
8
6
ns
4
2
ns
ns
0
1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 2011 2012 2013 2014
Figure 2. Average tree annual cone production in Central Range and Northern Plateau. Years with non-significant differences between
regions are indicated.
Forest Systems
December 2016 • Volume 25 • Issue 3 • e079
10
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
Table 4. Parameter estimates for the ZILN models for cone production for Northern Plateau and Central Range
Northern Plateau
Central Range
Covariate
Logistic
Part (α)
Log normal
Part (β)
Covariate
Logistic
Part (α)
Log normal
Part (β)
Intercept
0.3997
-0.8654
Intercept
–1.0177
0.4824
nv_7_3_std
mr_3_2_std
sp_4_2_std
jl_3_0_std
jl_7_0_std
N_hel_std
0.7697
0.6209
0.3389
0.3956
0.4772
0.0992
0.2780
0.9125
0.5316
0.0874
0.1529
0.2496
dc_10_3_std
my_2_2_std
sp_2_2_std
dc_4_1_std
jn_4_0_std
-0.3237
Lagged effects
lag_3_2_std
-0.5422
-0.3067
lag_3_std
-0.2029
Natural Unit
1
2
3
4
5
6
7
8
9
10
2.0233
0
-0.0375
–1.6412
–1.5405
–1.1506
-0.8612
–1.7904
-0.8290
-0.3490
0.4904
0
-0.2106
-0.8235
–1.0489
-0.3565
-0.1578
-0.6489
-0.7304
0.3302
21
22
23
24
0.2914
1.5677
1.3401
0
-0.3915
0.4458
0.1174
0
Stand covariates
1/dg
SDI
-58.0943
-0.0020
-25.4095
-0.0008
1/dg
SDI
-43.5477
-0.0020
-76.6994
-0.0010
Tree covariates
d/dg
dbh
2.3315
1.0931
0.0322
d/dg
3.0211
2.4468
Random components
σ2u (plot)
σ2v (plot)
σ2 (residual)
ρ (correlation
term)
0.8512
σ2u (plot)
σ2v (plot)
σ2 (residual)
ρ (correlation
term)
0.8544
Climate covariates
0.2541
1.0880
0.3307
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
Forest Systems
0.4972
0.2337
0.3244
0.1951
1.0118
0.2702
series 1980 – 2010, further use of the model will require these values (Table 5).
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).
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
11
Table 5. Mean annual and standard deviation values of the temporally explicit covariates entering
the models (required for standardization)
NP
Mean
STD
CR
Mean
STD
nv_7_3 (mm)
mr_3_2 (mm)
sp_4_2 (mm)
jl_3_0 (mm)
jl_7_0 (mm)
N_hel (days)
Cones (kg/ha)
239.5
102.8
89.2
92.6
241.0
3.41
194.526
78.8
59.7
45.0
43.5
76.0
3.91
164.689
dc_10_3 (mm)
my_2_2 (mm)
sp_2_2 (mm)
dc_4_1 (mm)
jn_4_0 (mm)
520.9
105.8
44.8
327.0
169.7
181.5
42.6
32.0
160.9
57.6
Cones (kg/tree)
4.605
3.512
Note: average annual cone production and standard deviation for NP are given in kg/ha, while in CR
are referred to kg/tree.
Table 6. Goodness of fit statistics of the marginal predictions (in real scale) of the fitted models over fitting data set
Northern Plateau
Tree
year–1
E (kg)
0.016
p-value (t-test) 0.8245
RE (%)
0.64%
wc_pred (kg) 2.472
wc_obs (kg)
2.491
EF (%)
31.37%
n
8041
RMSE (kg)
6.706
Central Range
Plot
year–1
Tree
Plot
Year
Tree
year–1
Plot
year–1
Tree
Plot
Year
-0.043
0.9540
-0.33%
12.879
12.836
49.14%
1555
29.548
-0.171
0.8944
-0.65%
26.284
26.112
49.86%
763
35.564
-0.477
0.9717
-0.34%
142.042
141.565
69.14%
141
158.709
0.026
0.8736
1.10%
2.346
2.372
89.47%
17
0.639
-0.264
0.0610
-5.77%
4.840
4.575
40.59%
2671
7.282
–1.312
0.2375
-5.76%
24.072
22.759
52.17%
537
25.736
-2.701
0.2813
-5.77%
49.528
46.826
56.71%
261
40.431
–13.558
0.4918
-5.70%
248.592
235.034
66.32%
52
140.485
-0.202
0.3884
-4.53%
4.658
4.455
93.14%
13
0.809
Where E: mean error, p-value: level of significance for null hypothesis E=0, RE: relative error; wc_pred and wc_obs: predicted and
observed weight of cones, EF: modelling efficiency, n: number of observations, RMSE: root mean square error.
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
Forest Systems
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.
Discussion
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
December 2016 • Volume 25 • Issue 3 • e079
12
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
8
PRED
OBS
NORTHERN PLATEAU
7
kg cones . tree–1 . year–1
6
5
4
3
2
1
0
1995
1997
1999
2001
2003
2005
2007
2009
2011
2013
2015
16
PRED
OBS
CENTRAL RANGE
14
kg cones . tree–1 . year–1
12
10
8
6
4
2
0
1995
1997
1999
2001
2003
2005
2007
2009
Figure 3. Observed (solid line) vs predicted (dotted line) value of average tree annual cone production, for the studied series in the two regions (1996-2014 in NP, 1996-2008 in CR). Predictions are the
marginal ones (assuming zero values for the random plot fixed effects). In 2011 and 2012 no cones
were collected in NP.
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
Forest Systems
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;
Isagi et al., 1997). 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 fe-
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
13
Table 7. Prediction accuracy for the six forests included in the independent set of evaluation
Region
Forest
Area (ha)
Obs (kg.
year-1)
E (kg.year–1)
p-value
RE (%)
RMSE
(kg.year–1)
EF
NP
NP
NP
NP
NP
CR
35
39
40
44
79
75
1275
1000
1477
1057
1078
811*
208545
50081
43863
157872
112350
69355
22314
6629
-20688
44995
28380
-22846
0.4475
0.2884
0.1390
0.1428
0.0941
0.0941
10.7%
13.2%
-47.2%
28.5%
25.3%
-32.9%
92028
19840
45652
100144
56226
45978
85.1%
90.1%
23.6%
71.1%
86.7%
60.5%
Where Obs: average annual cone production in the forest, E: mean error, p-value: level of siginifanca for the null hypothesis E=0, RE:
relative error, RMSE: root mean square error EF: modelling efficiency. *indicates that both observed and predicted values in forest 75
refers uniquely to the area commonly collected, which represents 25% of the total forest area.
male 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
Forest Systems
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,
December 2016 • Volume 25 • Issue 3 • e079
14
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
Forest 35
Forest 40
800000
250000
200000
600000
500000
400000
obs = 0,95 pred + 31604
R² = 0,86
300000
200000
Observed cone production (kg. year_1)
Observed cone production (kg. year_1)
700000
150000
100000
50000
obs = 0,61 pred + 476'
R² = 0,68
100000
0
0
0
100000
200000
300000
400000
500000
600000
700000
0
50000
100000
Predicted cone production (kg.year )
Forest 39
250000
Forest 44
500000
150000
100000
obs = 1,05 pred + 4481
R² = 0,91
50000
Observed cone production (kg. year_1)
200000
Observed cone production (kg. year_1)
200000
600000
250000
400000
300000
obs = 1,11 pred + 32625
R² = 0,78
200000
100000
0
0
0
20000
40000
60000
80000
100000
120000
140000
160000
0
180000
50000
100000
–1
150000
200000
250000
300000
350000
400000
450000
Predicted cone production (kg.year–1)
Predicted cone production (kg.year )
Forest 79
Forest 75 (Central Range)
600000
300000
500000
250000
400000
300000
200000
obs = 1,29 pred + 4097
R² = 0,95
100000
0
Observed cone production (kg. year_1)
Observed cone production (kg. year_1)
150000
Predicted cone production (kg.year–1)
–1
200000
150000
100000
obs = 0,9539 pred + 26042
R² = 0,76
50000
0
0
50000
100000
150000
200000
250000
300000
350000
–1
Predicted cone production (kg.year )
0
50000
100000
150000
200000
250000
300000
Predicted cone production (kg.year–1)
Figure 4. Agreement between observed and predicted values of annual cone production at forest scale, for the six studied forests of
the independent validation data set. Solid line represents the observed = predicted 1:1 functions. Dotted line represents the fitted model
observed = a + b . predicted, whose parameters and R2 are shown in the plot.
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
Forest Systems
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
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
15
Year T - 3
Year T - 2
Year T - 1
Year T (cones are collected in
November)
M A M J J A S O N D
J F M A M J J A S O N D
J F M A M J J A S O N D
J F M A M J J A S O N
NP_Occ
NP_Ab
CR_Occ
CR_Ab
Process
Bud
formation
Bud differentitation
Female
Female
flowering
flowering
survival
Damage by
freezing
Cone
enlargement
Cone
crop
Figure 5. Rainfall (dark gray) and freezing (pale gray) events affecting cone production through the whole flowering, fruiting and
maturation cycle in Northern Plateau and Central range. Physiological related processes are matched with the months of occurrence.
NP refers to Northern Plateau, CR to Central Range, Occ refers to submodel for occurrence of fruiting, Ab to submodel for abundance
of fruiting.
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.
Acknowledgements
The authors thank to the Regional Forest Services
of Valladolid, Ávila and Comunidad de Madrid their
continuous support in maintaining the experimental
Forest Systems
trials. The authors also wish to thank Rodrigo Gandía,
Luis Finat, David Cubero, Rebeca Martín, Enrique
Garriga and Ángel Bachiller for their appreciated help
in carrying out this research.
References
Bonet JA, Palahí M, Colinas C, Pukkala T, Fischer CR, Miina
J, Martínez De Aragón J, 2010. Modelling the production
and species richness of wild mushrooms in pine forests
of Central Pyrenees in northeastern Spain. Can J For Res
40: 347-356. http://dx.doi.org/10.1139/X09-198.
Bracalini M, Benedettelli S, Croci F, Terreni P, Tiberi R,
Panzavolta T, 2013. Cone and Seed Pests of Pinus pinea
L.: Assessment and Characterization of Damage. For
Entomol 106: 229-234.
Bravo F, Alvarez-Gonzalez JG, Del Rio M, Barrio M,
Bonet JA, Bravo-Oviedo A, Calama R, Castedo-Dorado F, Crecente-Campo F, Condes S, et al., 2011. Growth
and yield models in Spain: historical overview, contemporary examples and perspectives. Forest Systems
20(2): 315-328. http://dx.doi.org/10.5424/fs/201120211512.
Calama R, Montero, G, 2007. Cone and seed production from
stone pine (Pinus pinea L.) stands in Central Range
(Spain). Eur J For Res 126(1): 23-35. http://dx.doi.
org/10.1007/s10342-005-0100-8.
Calama R, Madrigal G, Candela JA, Montero G, 2007. Effect
of fertilization on the production of an edible forest fruit:
stone pine (Pinus pinea L.) nuts in SW Andalusia. Forest
Systems 16(3): 241-252. http://dx.doi.org/10.5424/
srf/2007163-01013.
Calama R, Mutke S, Gordo J, Montero G, 2008. An empirical ecological-type model for predicting stone pine (Pinus
pinea L.) cone production in the Northern Plateau (Spain).
For Ecol Manage 255 (3/4): 660-673.
Calama R, Tome M, Sánchez-González M, Miina J, Spanos
K, Palahi M, 2010. Modelling Non-Wood Forest Products
in Europe: a review. Forest Systems 19: 69-85. http://
dx.doi.org/10.5424/fs/201019S-9324.
December 2016 • Volume 25 • Issue 3 • e079
16
Rafael Calama, Javier Gordo, Guillermo Madrigal, Sven Mutke, Mar Conde, Gregorio Montero and Marta Pardos
Calama R, Mutke S, Tomé JA, Gordo FJ, Montero G, Tomé
M, 2011. Modelling spatial and temporal variability in a
zero-inflated variable: the case of stone pine (Pinus pinea
L.) cone production. Ecol Model 222: 606-618. http://
dx.doi.org/10.1016/j.ecolmodel.2010.09.020.
Calama R, Manso R, Tomé JA, 2012. ¿Cómo modelizar datos
con exceso de ceros? Métodos y aplicaciones a la investigación forestal. Cuadernos SECF 34: 55-65.
Calama R, Pardos M, Conde M, Madrigal G, Mutke S,
Montero G, Gordo FJ, 2015. Pérdidas de rendimiento en
piñón en piñas de Pinus pinea L.: análisis interregional.
III Reunión Científica de Sanidad Forestal SECF. Madrid
7-8 octubre 2015.
Cañadas MN, 2000. Pinus pinea L. en el Sistema Central
(valles del Tiétar y del Alberche): desarrollo de un modelo de crecimiento y producción de piña. Ph D Thesis.
Universidad Politécnica, Madrid.
Castellani C, 1989. La produzione legnosa e del fruto e la
durata economico delle pinete coetanee di pino domestico
(Pinus pinea L.) in un complesso assestato a prevalente
funzione produttiva in Italia. Annali ISAFA 12: 161-221.
Crone E, Miller E, Sala A, 2009. How do plants know when
other plants are flowering? Resource depletion, pollen
limitation and mast-seeding in a perennial wildflower.
Ecology Letters 12: 1119-1126. http://dx.doi.org/10.1111/
j.1461-0248.2009.01365.x.
Finat L, Gordo J, 2014. Sensibilización acerca del papel de
la propiedad pública y gestión forestal sostenible en la
provincia de Valladolid. Jornada proyecto PROPINEA,
Pedrajas de san Esteban, noviembre 2014. http://propinea.
es/wp-content/uploads/Sensibilizacion-acerca-del-papelde-la-propiedad-publica.pdf.
García Güemes C, 1999. Modelo de simulación selvícola
para Pinus pinea L. en la provincia de Valladolid. Ph D
Thesis. Universidad Politécnica, Madrid.
Gordo FJ, 2004. Selección de grandes productores de fruto
de Pinus pinea L. en la Meseta Norte. Ph D Thesis. Universidad Politécnica, Madrid.
Isagi Y, Sugimura K, Sumida A, Ito H, 1997. How does masting happen and synchronize? J Theor Biol 187: 231-239.
http://dx.doi.org/10.1006/jtbi.1997.0442.
Knops JMH, Koenig WD, Carmen WJ, 2007. A negative
correlation does not imply a trade-off between growth and
reproduction in California oaks. PNAS (USA) 104: 1698216985. http://dx.doi.org/10.1073/pnas.0704251104.
Koenig WD, Knops JMH, 2013. Large-scale spatial synchrony and cross-synchrony in acorn production by two
California oaks. Ecology 94(1): 83-93. http://dx.doi.
org/10.1890/12-0940.1.
Lambert D, 1992. Zero-inflated Poisson regression, with an
application to defects in manufacturing. Technometrics
34 : 1-14. http://dx.doi.org/10.2307/1269547.
Lesieur V, Yart A, Guilbon S, Lorme P, Auger-Rozenberg
MA, Roques A, 2014. The invasive Leptoglossus seed
bug, a threat for commercial seed crops, but for conifer
diversity?. Biol Invassions 16(9): 1833-1849. http://dx.
doi.org/10.1007/s10530-013-0630-9.
Mandallaz D, 2007. Sampling techniques for Forest inventories. Chapman & Hall/CRC. Applied Environmental
Forest Systems
Statistics. Boca Ratón, 256 pp. http://dx.doi.
org/10.1201/9781584889779.
Montero G, Cañadas N, Yagüe S, Bachiller A, Calama R,
Garriga E, Cañellas I, 2003. Aportaciones al conocimiento de las masas de Pinus pinea L. en los Montes de Hoyo
de Pinares (Ávila - España). Revista Montes 73.
Montero G, Calama R, Ruiz Peinado R, 2008. Selvicultura
de Pinus Pinea L. En Montero G, Serrada R, Reque J
(Eds.) Compendio de Selvicultura de Especies, pp 431470. INIA - Fundación Conde del Valle de Salazar. Madrid, España.
Mutke S, Gordo J, Climent J, Gil L, 2003. Shoot Growth and
Phenology Modelling of Grafted Stone Pine (Pinus pinea L.)
in Inner Spain. Ann For Sci 60(6): 527-537. http://dx.doi.
org/10.1051/forest:2003046.
Mutke S, Gordo J, Gil L, 2005a. Variability of Mediterranean stone pine cone production: yield loss as response
to climatic change. Agric For Met 132: 263-272. http://
dx.doi.org/10.1016/j.agrformet.2005.08.002.
Mutke S, Sievänen R, Nikinmaa E, Perttunen J, Gil L, 2005b.
Crown architecture of grafted stone pine (Pinus pinea L.):
shoot growth and bud differentiation. Trees-Structure and
Function 19(1): 15-25. http://dx.doi.org/10.1007/s00468004-0346-7.
Mutke S, Calama R, González-Martinez S, Montero G, Gordo
J, Bono D, Gil L, 2012. Mediterranean Stone Pine: Botany and Horticulture. Horticultural Reviews 39: 153-202.
Muttke S, 2013. Stone pine in Mediterranean forests. In:
Besacier C, Briens M, Duclercq M, Garavaglia V (coord.)
State of Mediterranean Forests 2013 (SoMF 2013), FAO
Silva Mediterranea, Rome, pp. 83-87.
Mutke S, Martínez J, Gordo J, Nicolas JL, Herrero N, Pastor
A, Calama R, 2014. Severe seed yield loss in Mediterranean stone pine cones. medPINE5, 5th International
Conference on Mediterranean Pines. Solsona, Spain.
Pardos M, Calama R, Maroschek M, Rammer W, Lexer MJ,
2015. A model-based analysis of climate change vulnerability of Pinus pinea stands under multi-objective management in the Northern Plateau of Spain. Annals of
Forest Science 72(8): 1009-1021. http://dx.doi.
org/10.1007/s13595-015-0520-7.
Prada MA, Gordo FJ, de Miguel J, Mutke S, Catalán-Bachiller G, Iglesias S, Gil L, 1997. Regiones de procedencia
de Pinus pinea. Ministerio de Medio Ambiente, Madrid,
Spain.
Sala A, Hopping K, McIntire EJ, Delzon S, Crone E, 2012.
Masting in whitebark pine (Pinus albicaulis) depletes
stored nutrients. New Phytol 196: 189-199. http://dx.doi.
org/10.1111/j.1469-8137.2012.04257.x.
Sánchez González M, Calama R, Cañellas I, Montero G,
2007. Variables influencing cork thickness in Spanish cork
oak forests: a modelling approach. Ann For Sci 67(2):
301-312. http://dx.doi.org/10.1051/forest:2007007.
Sork VL, 1993. Evolutionary ecology of mast seeding in
temperate and tropical oaks (Quercus spp). Vegetatio
107/108: 133-147. http://dx.doi.org/10.1007/978-94-0111749-4_9.
Strong W, 2015. Lodgepole pine seedset increase by mesh
bagging is due to exclusion of Leptoglossus occidentalis
December 2016 • Volume 25 • Issue 3 • e079
New models for cone production
(Hemiptera: Coreidae). J Entomol Soc Brit Columbia 112:
3-17.
Spathelf P, van der Maaten E, van der Maaten-Theunissen
M, Campioli M, Dobrowolska D, 2014. Climate change
impacts in European forests: the expert views of local
observers. Ann For Sci 71: 131-137. http://dx.doi.
org/10.1007/s13595-013-0280-1.
Forest Systems
17
Tu W, 2002. Zero-inflated data. In: El-Shaarawi AH, Piegorsch
WW (Eds.) Encyclopedia of Environmetrics. John Wiley
and Sons, Chichester, United Kingdom, pp. 2387-2391.
Vázquez-Piqué J, Pereira H, 2008. What to take into account
to develop cork weight models?: review and statistical
considerations. Forest Systems 17(3): 199-215. http://
dx.doi.org/10.5424/srf/2008173-01035.
December 2016 • Volume 25 • Issue 3 • e079