The data were collected from an experimental plot situated in a natural maritime pine ( P. pinaster) forest located at 755–810 masl in a flat sandy region in the Castilian Plateau, central Spain (latitude 41º16’N and longitude 4º29’W). Rotation age has been approximately 100 years under a shelterwood system adapted to resin production. In this system, the regeneration period varies from 20 to 25 years in length. Climate is dry Mediterranean, with an annual rainfall of 401 mm and mean annual temperature of 11.2 ºC. Vegetation is dominated by Corynephorus canescens (L.) P. Beauv., Stipa spp., Retama sphaerocarpa (L.) Boiss., Lavandula stoechas L., Thymus mastichina L., and isolated trees of Pinus pinea L. A detailed description of the site was given by Nanos et al. (2004).

In 2000, a circular plot (radius=100 m) was installed in a mature stand within the study area where only dead trees had previously been removed during thinning operations. On all mature trees (380) the following data were collected: spatial position (polar coordinates, to nearest 0.1 m), ring count age, total height (to nearest 0.1 m), number of cones and diameter at breast height (DBH) to nearest 0.1 cm. Coordinates were measured with a total station. In 2006 and 2007, additional cone counts were conducted on the same trees. Number of cones was estimated on each tree by visual count (see details in Miguel et al., 2002). For each tree, BAL, i.e., the basal area in trees larger than the subject tree (Wykoff, 1990), was computed as an index of competition from above and basal area of neighboring trees (BA ₆) was calculated as an index of two-sided competition. Both BAL and BA ₆ for each subject tree were computed by considering only the six nearest trees. Crown dimensions of each tree were estimated with Lizarralde’s (2008) crown equations (see Fig. S1 [suppl] for details). Tree dimensions (stem and crown) were assumed to have experienced no change over the course of the study (from 2000 to 2007) because the stand studied is mature and well over its maximum growth period.

In 2011, 113 trees were selected by expert judgement to cover the range in cone yield and were cored to obtain radial increments. DBH for these 113 trees ranged from 25.3 cm to 62.7 cm while cone counts ranged from 0 to 78 (see more details in Table 1). Cores were polished and scanned at high resolution (2000 dpi) with an Epson Expression 1640 XL scanner with 0.01 mm accuracy. Tree rings were measured using WinDendro© V 6.5C version (Regent Instruments). COFECHA program 6.06P version (Grissino-Mayer 2001, www.ltrr.arizona.edu) was used to assess data accuracy by calculating correlation indices between ring-width series and identifying errors such as missing or false rings. At the end of this process, a 15-year growth series (from 1996 to 2010) was obtained for each tree. Tree-ring width from the cone count years and the three years previous to each cone count were extracted from each growth series. Only 105 trees per year (104 in year 2000) were used in the final analysis. Trees were classified into four social classes according to their relative BAL: dominant (if BAL was in the bottom quartile, i.e. lowest values), codominant (if BAL was in the second to smallest quartile), intermediate (if BAL was in the third to smallest quartile) and suppressed (if BAL was in the largest quartile, i.e., highest BALs). The cone count frequencies are presented in supplementary documentation (whole data set in Fig. S2 [suppl]). As a surrogate for growth efficiency different alternatives were computed from annual ring increment, including basal area growth per unit crown projection area (CP), growth per unit crown surface area (CSA), and growth per unit LAI _s.

NAO data were obtained from the Climate Analysis Section, National Center for Atmospheric Research, NCAR (https://climatedataguide.ucar.edu/climate-data/hurrell-north-atlantic-oscillation-nao-index-station-based). Monthly and Seasonal NAO indexes are based on the difference of normalized sea level pressures (SLP) between Ponta Delgada, Azores and Stykkisholmur/Reykjavik. As of April 2003, the Ponta Delgada station stopped reporting, so to continue this time series NCAR staff substituted it by the nearest National Centers for Environmental Prediction/National Center for Atmospheric Research, NCEP/NCAR reanalysis grid point (see details at http://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.html). The correlation between a monthly time-series constructed with the Ponta Delgada observational data and one constructed with NCEP/NCAR reanalysis data was, according to NCAR information, 0.99 for the period from 1948 to 2003. NAO data from the cone count years (t=0) and the three previous years (t=1, t=2 and t=3) were extracted from the full dataset. NAO _winter was computed as the average of the consecutive December, January and February NAOs).

Determining correct variance-covariance structure

Because the data included three non-consecutive measurements from each tree, potential autocorrelation within trees could lead to biased estimates of variances and incorrect p-values for statistical tests. Four different variance-covariance structures, typically used in ecological studies (Barnett et al., 2010), were tested in a model for cone occurrence (binomial regression) and in a model for cone intensity (Poisson regression): independent, exchangeable, unstructured and autoregressive of order one. The models with alternative variance-covariance structures were fitted by SAS Proc GENMOD and compared by the quasi-likelihood information criterion (QIC) proposed by Pan (2001). No evidence could be found for significant autocorrelation in the pooled data, so further data analysis was pursued assuming independence of observations.

A zero-inflated Poisson (ZIP) model (eq. 1) proposed by Lambert (1992) was fitted to test the five hypotheses proposed in the Introduction. The ZIP model consisted of two parts: one to estimate the occurrence of cone production by using a logistic model and other to estimate the cone set or intensity of the reproduction event (number of cones) by using a Poisson model. The ZIP model fitted by SAS Proc COUNTREG had the following form:

where y_i is the observed cone count θ,represents the probability of observing no cones (y_i=0) modeled via a logit link function (eq. 2 below) and f(y) indicates the probability mass function for a Poisson function with expected mean μ (eq. 3) and variance equal to μ (1 + αμ) ².

where x _i are the covariates representing climate, tree size, stand density, and inter-tree competition and β _i and γⁱ are parameters estimated from the data. Expected cone occurrence (θ) and expected cone intensity (μ _i) are not assumed to be affected by the same covariates.

To test the first set of alternative hypotheses (see end of Introduction) about response of cone occurrence and cone intensity to tree size, stand density, inter-tree competition and climate, five corresponding models were developed (Table 2). The best among the five candidate models was selected on the basis of AIC.

To test the second general hypothesis about marginal effects of crown attributes, radial growth and tree growth efficiency, another set of alternative models was developed (Table S1 [suppl]). New variables were included only in the intensity part of the model. Different variables were initially tested in the occurrence part of the model, but no improvements were achieved, suggesting that the occurrence model should include neither crown variables nor any radial growth or growth efficiency effects.

The final step in the analysis was to identify the best full model ( i.e. the best model for explaining cone occurrence and cone intensity as a function of tree size, stand density, inter-tree competition, and climate), with the intensity portion of the model augmented by the best variables representing crown, growth and growth efficiency. The final model was selected on the basis of percentage reduction in AIC over the basic model (best model from Table 1). Percentage reduction in AIC (PRAIC) was computed as:

Additionally, Akaike weight ( W _i) was computed following the procedure presented by Johnson & Omland (2004). W _i provides a relative weight of evidence for each model and can be interpreted as the probability of model i is the best model for the observed data set. Equal values of W _i means that the performance of the models is equal.

Winter NAO three years prior to cone maturation was the only variable tested that consistently had a significant influence on occurrence of cones. In contrast, tree size, inter-tree competition, and local stand density variables all contributed to explaining variation in cone intensity, along with winter NAO either during the year of cone primordia initiation (year t-3 in models 2, 3, and 4; Table 3) or during year of flowering (year t-2 in model 5; Table 3).

In model 5, all variables were significant and AIC achieved the lowest value (Table 3). Model 5 was therefore further considered as the base model for testing hypothesis 2 with regard to marginal effects of crown attributes, growth efficiency, and radial growth on cone count (ZIP full models, see below), starting with the following estimated forms for occurrence and intensity, respectively:

Models with variables representing crown dimensions, radial growth and growth efficiency showed an improvement between 0.62 and 4.45 % in AIC reduction over the basic model (Table S1 [suppl]). The variables identified as best from each group (crown dimensions, radial growth and growth efficiency) were tested in different combinations to identify a best model that potentially included one variable from each group (Table 4).

Three of the combined models improved the performance of the basic model by including at least crown length (CL) as an additional explanatory variable in the basic model (Table 4). The best combined model (Basic+ CL + EFI_LAI _{s_t} in Table 4) included the two of the three variables that had the largest marginal effect on cone count when considered as the only variable added to the base model, i.e., CL, Growth ₀ and EFI_LAI _{s_0} (Table S1 [suppl]). AIC was reduced over 4.7% by this full model (Table 4). We conclude that cone occurrence and intensity were driven by climate (NAO _winter in the bud formation and flowering year), tree size (DBH), tree competition (BAL), tree density (BA ₆), and crown length (CL) and growth efficiency (EFI_LAI _{s_1}), as indicated in below in eqs. [7] and [8]:

Cone occurrence

When NAO _winter is relatively high, the winter is drier in the Iberian Peninsula. The logistic model for cone occurrence showed that a year with relatively high NAO _winter will lead three years later to a year with at least some cone production. In general, dry winters stimulate the development of reproductive buds in P. pinaster. This result is consistent with known relationships between environmental conditions and the life cycle of pine species. Mutke et al. (2005) found that rainfall during the fourth and the third years before cone maturation are key variables in cone production in P. pinea, a species that takes three years to produce mature cones (one year more than in P. pinaster).

Intensity of cone production

Four variables affected intensity of fruiting in our model: winter weather of the flowering year, tree size, stand density and competition.

Wet winters (indicated by lower NA _Owinter) during the flowering year (2) lead to higher cone yields. These results suggest that the necessary conditions for high cone counts include a relatively dry winter followed by a relatively wet winter. Dry winters ensure the occurrence of cones by favoring development of cone primordia, and a subsequent wet winter ensures a high cone count by enhancing primordia survival, flowering and pollination success (Fig. 1a).

Trees with larger diameter are predicted to produce a larger number of cones than smaller trees (Fig. 1b). This result has been previously documented in other species such as Pinus ponderosa (Krannitz & Duralia, 2004) and Pseudotsuga menziesii (El-Kassaby & Barclay, 1992), and is probably explained in part by the generally greater vigor of larger trees in even-aged stands. Higher tree vigor has also been shown to produce a higher number of recruits in next generation (González-Martínez et al., 2006). As thinning promotes diameter growth, an increase in cone production is expected as a longer-term benefit of thinning (Krannitz & Duralia, 2004).

As stand density declines the number of cones produced by a given tree increases (Fig. 1c). This result is consistent with previous observations of increased cone production after thinning (Verkaik & Espelta, 2006; Ruano et al., 2013), and can probably be explained in part by the greater resources available for each tree in less dense stands, allowing more carbon to be allocated to reproductive efforts.

Trees that are experiencing higher levels of competition were predicted to produce higher average cone counts than more dominant trees (Fig. 1d). This result seems to contradict the previously described effect of stand density, by suggesting that trees under more competitive stress from above may be induced to invest more resources in reproductive effort. Some previous studies suggest the opposite (Fowells & Schubert, 1956; Larson & Schubert, 1970), so in P. pinaster we similarly expected that dominant trees would produce the largest number of cones. The positive effect of BAL must be considered with caution, however, because P. pinaster stands in Central Spain are typically of very low density and tree vertical differentiation is weak, perhaps resulting in similar cone production for trees in different social classes (Fig. 1d). Furthermore, large trees (note dominant effect of tree diameter) may benefit from some surrounding tall trees if, for example, a larger proportion of cone primordia survive frost events with some protective cover.

No tradeoffs were detected between vegetative growth and cone production. However, a reduction of vegetative growth when cone production is high has been repeatedly found in different pine species (Linhart & Mitton, 1985; Krannitz & Duralia, 2004; Climent et al., 2008). In some cases, these tradeoffs have been proposed to have a genetic basis (Schmidtling, 1981; Linhart & Mitton, 1985). In natural stands of P. pinaster (this study), trees with higher growth efficiency also produced higher cone counts on average (Fig. 2), but the frequency distribution of cone counts was also wider than for trees with lower growth efficiency. In a study of Pinus contorta, Dick et al. (1990) found that presence of male cones reduced tree photosynthetic potential while presence of female cones did not. Tradeoffs between secondary growth and defense (Sampedro et al., 2011) and secondary growth and cone production (Climent et al. 2008) have been previously reported for Mediterranean pines evaluated under controlled conditions. Sampedro et al. (2011) found negative tradeoffs between growth and defensive investment in a greenhouse experiment with P. pinaster seedlings from different provenances. However, these latter tradeoffs emerged only under phosphorus limiting conditions. Climent et al. (2008) found that cone production in a Pinus halepensis seed orchard was incompatible with high vegetative growth while Almqvist et al. (2001) found no statistically significant correlation between cone production and height growth in Picea abies. Stands of P. pinaster in Central Spain are located on sandy soils of various depths and are characterized by high heterogeneity in available water and nutrients. Therefore, sites with good conditions for growth may also foster reproduction, while in more homogenous conditions (as those of greenhouse or common garden experiments), tradeoffs between growth and reproduction could emerge or simply be more detectable under the more uniform conditions.

In a detailed fruiting study in natural populations of Quercus lobata in southern California, Sánchez-Humanes et al. (2011) found that the cost of seed production differed at different levels of organization within the tree (tree canopy, individual branch and shoot within branch). Likewise the cost of seed production also varied among trees with different average seed production level. Only trees with consistently low seed production potential supported the trade-off theory, and only at the branch level. Their results suggest that any tradeoffs in reproductive allocation in natural populations would be detected at lower levels within trees, such as leaves or branches rather that at the level of the entire tree (Lovett-Doust & Lovett-Doust, 1988; Obeso, 1997). Our results for P. pinaster in a natural stand in Central Spain agreed with a lack of tradeoffs at the tree level, as also reported by Knops et al. (2007).

As an alternative possible microsite effects in extreme Mediterranean landscapes (see above), lack of apparent growth/reproduction tradeoffs in P. pinaster may be controlled by more general climatic factors. Cone production and growth interact with climate and this fact could have hidden tradeoffs in our study, as growth and cone development were driven by the same variables describing climate (Despland & Hoyles, 1997). Winter NAO is negatively related with radial tree growth over the Iberian Peninsula. However, Bogino & Bravo (2008) found that the growth variance explained by atmospheric indices at a finer scale among P. pinaster forests within the Iberian Peninsula was low (8.95 to 37.46%).

Climate, stand density and tree conditions (size and vigor, competition and growth efficiency) influence significantly both cone occurrence and intensity of fruiting as shown by a ZIP model. As the climate variables included in the model (based on NAO) are general and easily obtained, the proposed model has practical applicability to predicting Pinus pinaster cone production in the Iberian Peninsula. Our results are limited because we used just one monitoring site during three years, so specific studies are needed at local scales and in other geographical ranges of the species before the method can be applied at operational level. Judging from this model, no strong tradeoffs occurred between growth and reproduction in P. pinaster from Central Spain, so that phenotypic selection of those trees that have larger cone production as seed parents should not forfeit future stand growth in this region. However, trade-offs with other uses of photosynthates such as for defensive chemicals or structures (studied by Moreira et al., 2015 or Di Matteo & Voltas, 2016) should be also considered. Our model also points to silvicultural actions ( e.g. changes in stand density) that could increase cone production in natural forests of the species, which may prove useful in the face of impending environmental threats such as climate change