Plasticity in reproduction and growth among 52 range‐wide populations of a Mediterranean conifer: adaptive responses to environmental stress

A plastic response towards enhanced reproduction is expected in stressful environments, but it is assumed to trade off against vegetative growth and efficiency in the use of available resources deployed in reproduction [reproductive efficiency (RE)]. Evidence supporting this expectation is scarce for plants, particularly for long‐lived species. Forest trees such as Mediterranean pines provide ideal models to study the adaptive value of allocation to reproduction vs. vegetative growth given their among‐population differentiation for adaptive traits and their remarkable capacity to cope with dry and low‐fertility environments. We studied 52 range‐wide Pinus halepensis populations planted into two environmentally contrasting sites during their initial reproductive stage. We investigated the effect of site, population and their interaction on vegetative growth, threshold size for female reproduction, reproductive–vegetative size relationships and RE. We quantified correlations among traits and environmental variables to identify allocation trade‐offs and ecotypic trends. Genetic variation for plasticity was high for vegetative growth, whereas it was nonsignificant for reproduction. Size‐corrected reproduction was enhanced in the more stressful site supporting the expectation for adverse conditions to elicit plastic responses in reproductive allometry. However, RE was unrelated with early reproductive investment. Our results followed theoretical predictions and support that phenotypic plasticity for reproduction is adaptive under stressful environments. Considering expectations of increased drought in the Mediterranean, we hypothesize that phenotypic plasticity together with natural selection on reproductive traits will play a relevant role in the future adaptation of forest tree species.


Introduction
The timing of the onset of reproduction and the number of offspring produced by an individual are two fundamental life-history traits closely linked to fitness in an environment (Stearns, 1992;Braendle et al., 2011). According to life-history theory, individuals that start reproducing earlier in life tend to be favoured under harsh environments, due to reduced life expectancy (Roff, 1992). The initiation of reproduction in plants is often related to size rather than age (De Jong & Klinkhamer, 2005). For example, individuals should build a large vegetative body and invest all available resources in reproduction just before death, that is, a bang-bang strategy (King & Roughgarden, 1982). But uncertainty about the moment of death, for example, due to disturbances will tend to favour reproduction at smaller sizes (and younger ages), a graded reproductive investment and bet-hedging strategies (Childs et al., 2010). Thus, it is expected that plants, particularly long-lived perennials, will delay reproduction in favourable environments until they reach an optimal size for reproduction both by means of genetic change and phenotypic plasticity provided that selective forces act at local and broad scales (Kozłowski, 1992;Roff, 1992).
Experiments on herbaceous plants demonstrate that varying environmental factorsnamely resource availability and competitioninduce plasticity in reproductive strategies (Sultan, 2000;Weiner et al., 2009b;Anderson et al., 2011;Nicholls, 2011). In addition, plant populations are often genetically differentiated along environmental clines for size at reproduction and reproductive allometry, that is, the relationship between reproductive output and vegetative size (Lacey, 1988;Alexander et al., 2009;Guo et al., 2012). However, phenotypic plasticity of reproduction is driven to an important extent by size effects, as a strong positive relationship between vegetative and reproductive size is typically found and vegetative traits commonly respond plastically to environmental conditions. In comparison, phenotypic plasticity of the relationship between vegetative and reproductive size has been claimed to have a minor contribution to reproductive output, but this is still debated (Weiner et al., 2009a,b).
Long generation time in long-lived perennials implies that the same genotypes cope with year-to-year changing environmental conditions. On the other hand, populations of annuals or short-lived perennials can undergo genetic changes in shorter periods (Franks & Weis, 2008). Therefore, plasticity might be of greater importance as an adaptive strategy in trees and woody plants compared with short-lived plant species (Willson, 1983) such that long-lived species might exhibit plasticity in both vegetative (Chambel et al., 2005) and reproductive traits like size at reproduction and reproductive investment. The few studies published on long-lived species highlight strong selection on the threshold size at first reproduction and the allometry of reproduction, leading to genetic differentiation at large spatial scales (Thomas, 1996;Matziris, 1997;Niklas & Enquist, 2003;Climent et al., 2008;Santos-del-Blanco et al., 2010) and promoting phenotypic plasticity in life histories at local scales (Fang et al., 2006). Despite consistent predictions of plasticity in the threshold size of reproduction, little is known about the costs of plasticity in terms of final reproductive output relative to vegetative size (Roff, 2000). Reproductive efficiency (RE) can be defined as the slope of the reproductive-vegetative size developmental trajectory that connects threshold size for reproduction with reproduction at a given developmental stage or at the onset of senescence (Bonser & Aarssen, 2009). It is expected that early reproduction will imply lower RE, modifying reproductive allometries and, in turn, reproductive variability within and among populations. However, this has rarely been tested even in shortlived semelparous species (but see Bonser et al., 2010).
The Mediterranean pine Pinus halepensis Mill. (Aleppo pine) is a suitable model species for testing hypotheses on the evolution of reproductive strategies in long-lived perennials. It is precocious, bearing female cones from as early as 3 to 6 years of age, and commits heavily and regularly to reproduction, most notably female reproduction (Ne'eman et al., 2004). Pinus halepensis is widespread over a large circum-Mediterranean distribution area, and low population differentiation in neutral markers has been reported in the Iberian Peninsula due to recent range expansion (Soto et al., 2010). Pinus halepensis shows a wide ecological breadth among populations and is putatively adapted to a large range of abiotic stressors and perturbations, particularly fire and drought (Ne'eman et al., 2004), although intense drought episodes might be detrimental to reproduction (Girard et al., 2011). However, information regarding among-population variation in phenotypic traits in this species remains scarce.
Previous works described significant ecotypic differentiation for size at maturity in P. halepensis (Climent et al., 2008;Santos-del-Blanco et al., 2010). In this study, we focus on phenotypic plasticity and among-populations genetic variation in plasticity for reproductive allometry in range-wide P. halepensis populations assessed in a common garden experiment replicated in two contrasted sites (low and high environmental stress). Our objectives are to (i) assess the existence of phenotypic plasticity for size at maturity and the reproductive-vegetative size (R-V) relationship in range-wide populations subject to contrasting field conditions; and (ii) to compare genotype 9 environment patterns for vegetative and reproductive traits and correlations between both sets of traits representing trade-offs that might describe adaptive strategies. First, we expect that similar environmental cues defining favourable or unfavourable growth conditions will act in the same direction considering genetic differentiation and plasticity (Anderson et al., 2012;Chevin et al., 2012). Based on life-history theory, this would imply that the more stressful the environment (both at the origin of populations and at the trial site), the greater amount of resources would be devoted to reproduction. Specifically, we expect that environmental stress will induce reproduction at smaller sizes, associated with higher slopes of the R-V relationship. Finally, we also expect reproductive strategies to be governed by trade-offs between precocity and lifetime fitness, so that individuals that reproduce late benefit from a higher lifetime reproductive investment relative to their size. eastern and central Spain. The trial includes 52 native populations from continental Spain, Balearic Islands (Spain), France, Italy, Greece and Tunisia, as well as four non-native populations (see Climent et al., 2008 for details) (Fig. 1, Table S1), thus covering most of the species' range. Only data relative to native populations were used in the present study. The minimum requirements for assessing plasticity in our experiment were, first, data measured at same age and identical protocols between sites and, second, contrasted enough environments. Only two of the six sites fulfilled both requirements.
Summary data of environmental conditions at both trial sites obtained from a functional model (Gonzalo-Jim enez, 2010) are shown in Table 1. Valdeolmos trial site (hereafter 'low-stress site') has sandy loam deep soil, whereas Rinc on de Ademuz trial site (hereafter 'high-stress site') has shallow and rocky soil. In addition, mean annual rainfall is ca. 25% higher in the low-stress site, and winters are slightly warmer compared with the high-stress site. As a result of combined effects of poorer soil, lower rainfall and slightly colder winters, the high-stress site is much more limiting for P. halepensis vegetative growth compared with the lowstress site. This constant environmental difference between sites should not be confounded with withinsite year-to-year meteorological variation that has been previously described in this species (Girard et al., 2011).
Population seedlots were obtained by bulking openpollinated seeds from a subsample of 20 to 30 trees spaced at least 100 m apart in each population. At both sites, 832 one-year-old seedlings from native populations were planted in 1997 in a row-column design on the intersections of a 2.5 9 2.5 m grid, with four replicates and four contiguous plants per population and replicate (16 trees per population). One replicate in the low-stress site was lost due to rabbit herbivory and was not included in this study (624 trees remaining). Due to other causes of mortality, final sample size for this study was 589 in the low-stress site and 633 in the high-stress site.

Measurement of traits and environmental variables
We measured height for each tree at ages 7, 11 and 13 years for both sites (2003, 2007 and 2009, respectively). Diameter at breast height was measured at both sites at ages 11 and 13 years and used to infer biomass from allometric models (Montero et al., 2005) ( Table 2).
The onset of female and male reproductive functions in P. halepensis is decoupled, with trees generally starting reproduction as females (protogyny) and male reproduction being delayed for up to several years (Shmida et al., 2000). Thus, we focused on the study of P. halepensis early investment in female function to estimate threshold sizes for first reproduction and reproductive investment.
Female cones in P. halepensis remain attached to the branches even after dehiscence (normally also delayed several years, Ne'eman et al., 2004). Differences in size and colour allow discrimination of several cohorts within tree crowns (Ne'eman et al., 2011) and therefore enable retrospective record of female reproduction. Up to three successive cohorts of female cones were counted at ages 7 and 13 years (2003 and 2009) therefore dating back to the very first reproductive events up to the generalization of reproduction at both sites.
We defined the cumulative reproductive investment (CRI) as the sum of all counted female cones produced by an individual until last measurement at age 13 (Table 2).
We collected data for six meteorological and three spatial variables describing the environmental conditions found in the natural populations (Table 1, Table  S1). Meteorological data for Iberian populations were obtained from a functional model (Gonzalo-Jim enez, 2010), and data for other populations (i.e. Balearic Islands, France, Tunisia, Italy and Greece) were obtained from WorldClim-Global Climate Data at 5′ resolution (Hijmans et al., 2005). To test hypotheses of local adaptation, we also calculated the Gower's ecological distance for each population at both trial sites (Rutter & Fenster, 2007). This adimensional index informs about the environmental distance between the native environment of each population and the environment where they were grown in the common garden. The analysis was limited to Iberian and Balearic populations due to the unbalanced number of eastern Mediterranean populations in the experiment.

Data analysis
All reported models and tests were implemented in R (R Development Core Team., 2012) using packages lme4 (Bates et al., 2011) and MCMCglmm (Hadfield, 2010).

Survival and vegetative size
Mixed linear models for size (height and biomass) at age 13 and generalized linear mixed models (logit link, binomial error) for survival and proportion of reproductive individuals at age 13 were fitted. In all models, site effect was treated as fixed. Population, site-by-population interaction and replicate within site were treated as random. A common interpretation of model parameters is as follows: significant differences among populations indicate intraspecific genetic variability; significant differences between sites reflect phenotypic plasticity, and significant site 9 population interaction indicates genetic variation for plasticity among populations (Schlichting, 1986). However, deviations from that framework need also to be considered. For example, environmental factors can significantly affect seeds  *7 corresponds to ages 5 and below, 6 and 7 years; 13 corresponds to 11, 12 and 13 years.
during development, causing epigenetic changes in gene expression (Johnsen et al., 2005). Also, a significant site 9 population interaction can indicate local adaptation if populations have a better performance in the site most similar to the conditions of their site of origin (Vergeer & Kunin, 2013). To test the significance of site, population and site 9 population terms, we performed likelihood ratio tests (LRTs) comparing full models containing all terms with those lacking the relevant term to be tested. Variance components and adjusted means for size at final measurement of each population were derived from analogous models fitted for each trial site.

Size at first reproduction
A generalized linear mixed model (logit link, binomial error) was fitted for cumulative female reproduction (either present or absent) data at ages 7 and 13. We included height as covariate and height 9 site interaction as fixed term. Then, independent models were fitted for each population and site by Markov Chain Monte Carlo (MCMC) methods (see Santos-del-Blanco et al., 2012 for further details), and median threshold size for first reproduction (TSFR) was defined as the size at which the probability for a tree to have reached sexual maturity was 50% (Wesselingh et al., 1997) and computed by dividing slope by intercept estimates. We also calculated the size of the smallest reproductive individual (SRI) at each population and used this information to classify nonreproducing trees into juvenile (smaller than SRI) or vegetative (larger than SRI) (Mendez & Karlsson, 2004).

Fecundity and reproductive-vegetative size relationships
Generalized linear mixed models (log link, Poisson error) were fitted to the CRI. Juveniles were removed from the data set prior to analysis. The models also included an individual-level random effect to model additive overdispersion (Elston et al., 2001). Reproduction in plants is typically size dependent (Niklas & Enquist, 2003;Weiner, 2004). We accounted for size-dependent differences in reproductive allocation by calculating reproductive allocation per population first as the mean value across individuals (CRI/biomass) and second as the expected reproductive value based on fitted reproductive-vegetative size (R-V) models and then divided into average size. Thus, each approach represents the mean reproductive allocation per population and reproductive allocation of an average-sized individual in a population, respectively. Similar values for both indexes would indicate that the estimation of reproductive allocation is robust, although issues remain about spurious correlations with size.
Reproductive-vegetative size models describe the relationship between reproductive and vegetative allocation using two parametersan intercept and a slope (Weiner et al., 2009a;Guo et al., 2012). However, depending on the range of sizes used to fit the models, those two parameters might not be independent in a set of populations due to collinearity (Pinheiro & Bates, 2000). Thus, to summarize reproductive output and compare R-V relationships while accounting for tree vegetative size, we fitted generalized linear mixed models (log link, Poisson error) with independent intercept and slopes to CRI data according to: where g i is the linear predictor linked to the expected value of natural logarithm of CRI [ln (l i )]. x′ represents the design matrix containing the values for the fixed size effects. b is a vector containing the fixed intercept and slope associated with size, to be estimated. z′ is the design matrix for the random populations effects. b is the vector of random coefficients that follow a normal distribution. e is the vector containing the errors that follow a Poisson distribution. Two models were fitted per site, the first one with b containing random effects for intercepts and the second containing random effects for the slopes. AIC values from both models at each site were very close, indicating that either random intercept or random slope models had similar explanatory power. Population-adjusted intercepts and slopes were derived from MCMC models fitted at each site and used as fecundity indicators; this allowed us to compare general estimates from both sites. Random intercepts associated with population reflect constant deviations across sizes from the general model, that is, a constant higher or lower commitment to reproduction across sizes. Random slopes associated with population represent deviations on reproductive output proportional to vegetative size, that is, enhanced or decreased commitment to reproduction along vegetative size. Our analysis of size at first reproduction and reproductive output divided in two steps (binomial and Poisson submodels) was thus similar to a hurdle model (Brophy et al., 2007;Haymes & Fox, 2012).

Reproductive efficiency
We estimated RE as the slope of the size-reproduction developmental trajectory, linking vegetative size at first reproduction and vegetative and reproductive output at final development (age 13) (Bonser & Aarssen, 2009). RE was estimated at the population level for both sites. We tested whether there were significant correlations between the threshold size for reproduction and RE at each site and whether RE was affected by the environment, comparing RE values between both sites with a paired t-test.

Local adaptation patterns
Pearson's correlation tests at each trial site were used to test the relationship between Gower's distance and fitness. We used CRI and female TSFR, as the variables most closely related to fitness but also explored the correlation between Gower's distance and vegetative growth traits (Leimu & Fischer, 2008). We also tested whether increased environmental distances were correlated with changes in trait means.

Plant trait correlations and ecotypic trends
We calculated Pearson's correlations among plant traits at the population level and among those traits and environmental conditions found in the natural populations. Correlations among plant traits can be interpreted as genetic correlations modified by common environmental effects. Correlations were conducted at each trial site separately to check whether trait-trait correlations and ecotypic trends of variation were site dependent. We also obtained the site-to-site correlations for phenotypic traits, as a double-check of site-by-population interaction (Pigliucci, 2001).

Vegetative traits
We found that plants in the high-stress site had lower biomass and height compared with those in the lowstress site, thus confirming that overall environmental differences between both sites had an effect on vegetative growth (Table 3, Fig. 2). In addition, tree survival was significantly lower in the high-stress site compared with the low-stress site (v 2 1 ¼ 76:2, P < 0.001) ( Table 3). All populations attained larger sizes in the low-stress than in the high-stress site. However, there was no evidence for population effect alone, but differences between populations were site specific, and a significant site-by-population interaction was found for all vegetative traits (Table 4, Figs 2 and 3), consistent with genetic variation in plasticity for vegetative traits among populations and, possibly, local adaptation. Among-population variance was larger for biomass in the low-stress site [45.0 (31.6-69.2)] compared with the high-stress site [4.7 (3.3-7.2)], but no significant differences were found for height [829 (582-1275) lowstress site; 797 (560-1225) high-stress site]. Population means for vegetative traits at the low-and high-stress sites can be accessed in Tables S2 and S3, respectively.

Reproduction and threshold sizes
Mean size of reproductive individuals was greater than that of nonreproductive ones at both sites and both years (low-stress site: biomass v 2 1 ¼ 21:0, height v 2 1 ¼ 34:7; high-stress site: biomass v 2 1 ¼ 85:8, height v 2 1 ¼ 242:1, all tests P < 0.001). At the early measurement date (age 7 years), the proportion of reproductive individuals was slightly greater in the low-stress site than in the highstress site (v 2 1 ¼ 58:6, P < 0.001). However, at 13 years of age, 96% of trees were reproductive in the stressed environment, whereas only 84% were in the more favourable environment (v 2 1 ¼ 34:7, P < 0.001). As a result, at final measurement, the number of vegetative individuals was higher in the low-stress than in high-stress site (34 vs. 5).
We found a significant effect of both site and population on the threshold size for reproduction, as shown by the significant site and population terms (Table 4). Thus, threshold size for reproduction is both highly plastic and variable among populations (Figs 2 and 3). By contrast, the site-by-population interaction term was not significant, indicating that there was no significant genetic variation for plasticity in the threshold size for reproduction among populations. The probability of reproducing at a given size was significantly smaller in the low-stress site than in high-stress site, evidenced by a reduced slope of the model (data not shown).
We were able to fit independent threshold models for all but three populations in the low-stress site and all populations but one in the high-stress site. For all but two populations, the point estimate of the threshold size for reproduction was higher in the low-stress site than in the high-stress one (Figs 2 and 3).
For CRI, site and population effects were significant, but not, although marginally, population-by-site interaction (Table 4). When tree biomass was included as a covariate in the R-V model, site (indicating a different R-V relationship in both sites), and population terms were significant, but not site 9 population interaction (Table 4). When height was used as a covariate, similar results were obtained although site was not significant ( Table 4). The subsequent GLMM models fitted by MCMC aimed at estimating fecundity at the population level while controlling for size effects revealed an enhanced reproductive allocation in the high-stress site with respect to the low-stress one, defined by a larger intercept and slope (Table 3). Here, a positive intercept must not be regarded as biologically implausible, as it represents a population, not an individual developmental trajectory. Mean reproductive allocation per population and expected reproductive allocation of an average-sized individual per population yielded similar results (data not shown), so we used the first index to describe reproductive allocation (RA) as it was derived primarily from the data. Consistently with fecundity descriptors, RA was larger in the high-stress than in the low-stress site. Population means for reproductive traits at the low-and high-stress sites can be accessed in Tables  S2 and S3, respectively.

Plant trait correlations and ecotypic trends
Differences in CRI were related to the Gower's environmental distance at both trial sites. That is, there was a significant positive relationship between the environmental similarity between each population with respect to the common garden and the number of cones it produced as revealed by Pearson's correlation tests (Table 5). Female threshold size for reproduction was also significantly negatively related to Gower's distance in the low-stress site but only marginally in the high-stress site. By contrast, for vegetative traits, even when correlations were only marginally significant, they showed opposite patterns at each trial site. Correlations at the high-stress site were negative and at the low-stress site were positive ( Table 5). The threshold size for reproduction was the only variable significantly correlated with a change in environmental distance. Closer distances were related to larger thresholds for reproduction. Correlation with CRI was negative but nonsignificant (Table 5).

RV Slope
High stress site Low stress site Fig. 2 Site-site graphs for the interpretation of phenotypic plasticity in several vegetative and reproductive traits measured in two Pinus halepensis provenance common gardens subject to contrasting environmental conditions (low and high environmental stress). Points represent mean values per source population. Units are as defined in Table 2. Site-site correlations of population-adjusted means were not significant for vegetative traits, but strong positive correlations were found for reproductive traits (Table S4). This corroborates the high site-by-population interaction for vegetative traits vs. nonsignificant site-by-population interaction for reproductive traits seen by previous analyses (Table 4).
Within sites, correlations among traits representing potential trade-offs between growth and reproduction (e.g. height and reproductive allocation) showed, in general, stronger correlations in the high-stress site than in the low-stress one (Table S4). We found significant positive correlations between height and TSFR in the high-stress site, but not in the low-stress one. All other correlations between vegetative and reproductive traits were nonsignificant. We also found extensive correlations within vegetative traits and reproductive traits (Table S4).
We found no significant difference in RE between sites (t 47 = 1.442, P = 0.156), nor for growth above the site-specific median threshold size for reproduction (t 47 = À1.037, P = 0.305). RE was negatively correlated with TSFR in the low-stress site but not in the highstress site (Table S4).
Correlations between environmental factors and plant traits, reflecting ecotypic trends of variation, were higher in reproductive traits compared with vegetative traits. In turn, they were higher in the low-stress site, but in both sites, the sign of the correlation was the same. Traits indicative of more precocious or abundant reproduction were related to higher altitude, and warmer summers and colder winters (therefore higher continentality index). However, no significant correlations were found between rainfall and any of the phenotypic traits measured at either site (Table S4).

Discussion
Our experiment showed that at the more stressful site, P. halepensis trees started reproducing at smaller sizes and completed female reproductive maturity earlierboth in size and in timethan at the least stressful site, therefore following theoretical expectations (Roff, 1992). By definition, threshold size for reproduction accounts for differences in size, so a plastic response in this trait should be considered as a true plastic response and not driven solely by plasticity in size (Sugiyama & Bazzaz, 1998;Weiner, 2004). Hypotheses regarding plasticity of threshold size for reproduction have been addressed in plants only in few cases (Bonser & Aarssen, 2009;Kagaya et al., 2009;Bonser et al., 2010), in contrast with predictions of the high relevance of this type of plasticity (Burd et al., 2006). However, evidences pointing at this phenomenon are common through the literature (Nagy & Proctor, 1997;Fang et al., 2006). We relied on two natural environments to test our hypothesis, which also allowed us to study local adaptation patterns. Nonetheless, a more precise control of environmental stress could be achieved by artificially inducing drought or watering or by setting the experiment at contrasting soil depths and/or nutrient levels, for example, leading to more general conclusions.
The adaptive value of reproduction at larger sizes in favourable conditions relies on a positive relationship between fecundity and size at reproduction, so that Table 4 Results of general and generalized linear mixed models for Pinus halepensis vegetative and reproductive traits measured in two experimental sites with contrasting environmental conditions. Full models were fitted including all terms. Site, population and population-by-site models were fitted excluding the relevant terms to test plasticity, genetic variation and genetic variation for plasticity. log-likelihood (logLik) is given for each model. Chisquare statistic (Chisq), degrees of freedom (d.f.) and P-value are given for likelihood ratio tests between the full model and reduced models. attaining a larger size implies an increased lifetime reproductive output (Metcalf et al., 2003). Although this relationship is clear in semelparous species, as they only have one reproductive event in their life (bangbang strategy), in iteroparous species like trees the relationship is not straightforward because individuals allocate significant amounts of resources to maintenance each season throughout their lives (De Jong & Klinkhamer, 2005). For example, two trees with a similar adult size might differ in reproductive output due to differential investment in maintenance along their lives, whereas in annual plants, no differences in reproductive output are expected for similar-sized individuals (Weiner et al., 2009a). In plants, favourable environments for growth generally favour lower reproductive investment relative to size (Matyas & Varga, 2000;Ortiz et al., 2011;Haymes & Fox, 2012). As these conditions are typically associated with increased competition (Grime, 1977), delayed reproduction in these environments could be driven both by a positive relationship between size at reproduction and lifetime reproductive investment, but also by a persistent allocation to growth and maintenance in crowded stands (Zhang, 2006).
At the population level, a low threshold size for reproduction was correlated with steeper slopes for the R-V relationship at both trial sites. Thus, genetic and  Table 2.  environmental factors promoting allocation to reproduction acted consistently along the development of the trees in our experiment. Trade-offs between reproduction and growth are predicted to be more relevant in limiting conditions (Karlsson & Mendez, 2005). Accordingly, in the more stressful site, we found a negative correlation between vegetative traits and reproductive precocity (hence positive with TSFR), and correlation coefficients between reproductive allocation and size among populations were higher than in the less stressful site. However, contrary to our expectations, we found no differences for RE between sites; that is, delayed onset of reproduction in the low-stress site was not rewarded with a proportionally higher reproductive output. We did not extensively test for adaptive plasticity by inducing plants to express an inappropriate phenotype in a given environment (see Sultan, 2000). Nonetheless, our findings regarding enhanced reproduction at the more stressful site highlight the adaptive value of plasticity for reproduction (Anderson et al., 2012) as they support theoretical expectations (Pigliucci, 2001). Indeed, if reproduction at the high-stress site had followed the same allometric trend as in the low-stress site, the risk of becoming locally extinct after a severe disturbance would be dramatically higher. However, at the same time, our results raise uncertainty about what the benefits of delayed reproduction are in environments favourable to vegetative growth. Also, as we did not consider male reproduction in our analysis as the trees were still young, we cannot rule out a possible trade-off between female and male reproduction, so that female reproduction is reduced in favourable environments, but male reproduction could be enhanced. To gain better insight into these uncertainties, data covering a longer time period for both sexual functions would be needed. In our experiment, tree size was a weak predictor of reproductive output in both trial sites. Although our trees were young and the relationship may strengthen with time (Weiner et al., 2009a), a loose, although significant, relationship between size and reproduction is common in perennial species, notably in trees (Climent et al., 2008;Haymes & Fox, 2012). Population 9 site interaction was important for vegetative traits, but we found no evidence for larger plants corresponding to shorter environmental distances (Vergeer & Kunin, 2013). Actually, for the low-stress site, the trend was opposite, being larger plants from ecologically distant populations. Instead, reproductive output did show a negative relationship with the environmental distance of the original populations to each trial site. As reproductive output is closely linked to fitness, this suggests that populations have adapted to local climate conditions, and climate is important in controlling the expression of reproduction (Leimu & Fischer, 2008). Thus, we advise against using tree size as a proxy for fitness and encourage the use of reproductive output in tree evolutionary ecology studies.
Several additional factors might interact with raw reproductive output to define individual fitness (Braendle et al., 2011). Within populations, some individuals remained nonreproductive well above their population TSFR, a phenomenon also described in biennials (Wesselingh & Klinkhamer, 1996), and the highest reproductive output was typically achieved by mediumsized individuals in consistency with other experiments in this species (Climent et al., 2008). This pattern was more evident in the low-stress site, where a higher number of trees remained vegetative. A likely explanation for this observation is a diversifying bet-hedging strategy (Simons, 2007), with individuals reproducing according to a genetically determined allocation curve and others situated below that curve (Weiner, 2009a). If a disturbance occurred at either trial site, population and individual would be the most important factors determining the number of available seeds for the next generation. This would imply that if the primary reason for delaying reproduction was a larger future reward through increased size and greater potential future reproduction, many individuals would be making a nonprofitable investment. However, an enhanced allocation to growth would also increase fitness through an increased likelihood of survival (Zhang, 2006). The relative importance of these nonexclusive explanations deserves more attention to better understand adaptive responses in trees.
Contrary to expectations, reproductive output for the whole set of populations was very similar between the two contrasting environments. Plasticity for cumulative cone production was much lower (up to twofold) than that for biomass (up to 10-fold) (Fig. 2). Reproductive output emerged from a combination of plastic responses in growth (larger in the less stressful site) and allometry (higher reproduction for a given size in the more stressful site). An ecotypic trend of enhanced reproduction towards higher altitudes and more extreme temperatures, already described in Climent et al. (2008), was not related to population differentiation in plasticity. Interestingly, we found plasticity for both reproductive allometry and vegetative traits, but only genotypeby-environment interaction for vegetative traits.
Phenotypic plasticity is expected to arise in environments that change in a predictable fashion (Van Kleunen & Fischer, 2005). Within species, higher plasticity is generally expected in populations subject to greater interannual variance in precipitation and extreme temperatures and also those living in more patchy environments (Sultan & Spencer, 2002;Baythavong, 2011). In addition, traits might differ in their sensitivity to the environment, or may be constrained resulting in some being more plastic than others (Matesanz et al., 2010), as is the case in our experiment, where not only phenotypic plasticity but also its variation among populations was higher for growth than for reproduction, consistent with findings in different plant genera (reviewed in Weiner et al., 2009a). In our experiment, the lack of population differentiation for plasticity for reproductive allometry could be due to nonexclusive causes such as (i) a strong stabilizing selection for plasticity of reproductive allometry among populations, (ii) a canalization or total dependence of reproductive traits on vegetative traits, like internal cues, or (iii) the perception of environmental heterogeneity, selecting for plasticity, differing between reproductive and vegetative traits. For example, vegetative traits may be more dependent on fine-grain variability of soil depth or nutrient availability, but reproductive traits depend more on factors acting at a larger scale like climate and severe perturbations. In the high-stress site, variation in responses for vegetative traits was constrained, whereas in the low-stress site, among-population differences were neatly expressed, revealing cryptic genetic variation (Schlichting, 2008). Reproductive traits, however, displayed similar levels of variation at both trial sites, so in the environment with most limiting conditions, variation for reproductive traits was more relevant than that for vegetative traits.
Considering expectations of increased drought in the Mediterranean due to climate change (Lindner et al., 2010) and assuming a high heritability of reproductive allometry (Santos-del-Blanco et al., 2010 and in prep;Wesselingh & De Jong, 1995), we hypothesize that phenotypic plasticity coupled with subsequent natural selection on this trait (Anderson et al., 2012;Chevin et al., 2012) will play a relevant role in future adaptation of forest species.

Supporting information
Additional Supporting Information may be found in the online version of this article: Table S1 List and location of the Pinus halepensis populations comprised in the present study and planted at two trial sites with contrasting environmental conditions (high and low environmental stress). Table S2 Vegetative and reproductive traits in 52 natural Pinus halepensis populations grown in a common garden placed at Valdeolmos (Madrid, Spain), referred as low stress site. Table S3 Vegetative and reproductive traits in 52 natural Pinus halepensis populations grown in a common garden placed at Rinc on de Ademuz (Valencia, Spain), referred as high stress site. Table S4 Pearson correlation coefficients for the correlation at the population level between climatic variables and Pinus halepensis traits in each site, the low stress site above diagonal and the high stress site below the diagonal.