Intraspecific responses to climate reveal nonintuitive warming impacts on a widespread thermophilic conifer

Many ecologically important forest trees from dry areas have been insufficiently investigated for their ability to adapt to the challenges posed by climate change, which hampers the implementation of mitigation policies. We analyzed 14 common-garden experiments across the Mediterranean which studied the widespread thermophilic conifer Pinus halepensis and involved 157 populations categorized into five ecotypes. Ecotype-specific tree height responses to climate were applied to projected climate change (2071–2100 AD), to project potential growth patterns both locally and across the species’ range. We found contrasting ecotypic sensitivities to annual precipitation but comparatively uniform responses to mean temperature, while evidence of local adaptation for tree height was limited to mesic ecotypes. We projected intriguing patterns of response range-wide, implying either height inhibition or stimulation of up to 75%, and deduced that the ecotype currently experiencing more favorable (wetter) conditions will show the largest inhibition. Extensive height reductions can be expected for coastal areas of France, Greece, Spain and northern Africa. Our findings underline the fact that intraspecific variations in sensitivity to precipitation must be considered when projecting tree height responses of dry forests to future climate. The ecotype-specific projected performances call for management activities to ensure forest resilience in the Mediterranean through, for example, tailored deployment strategies.


Introduction
Trees are crucial organisms in many terrestrial ecosystems, and their responses to climate variations could influence the global carbon and water cycles considerably and affect the ecosystem services that forests provide (Viglizzo et al., 2016;Sheil, 2018). Conifers, in particular, show strongly differentiated populations (Neale & Wheeler, 2019) that result from long-term selective processes, the history of niche recolonization after glaciations, and the trade-offs derived from patterns of covariation among fitness-related traits (Kawecki & Ebert, 2004). A high degree of intraspecific genetic variability is essential for the resilience of forests (Alberto et al., 2013). However, the importance of climatic adaptation of forest trees, which involves both standing genetic variation and differential phenotypic plasticity of populations to buffer environmental changes (M aty as, 1996), has been largely overlooked when forecasting responses to climate change. This omission may result in unrealistic projections for the status and distribution of forest species in the forthcoming decades (O'Neill et al., 2008;Reich & Oleksyn, 2008;Peterson et al., 2019).
The analysis of adaptive strategies has focused primarily on forest trees of high economic relevance in the context of breeding programs (e.g. Pinus contorta; Montw e et al., 2016). However, the ability of many ecologically important species to adapt to the challenges posed by climate change remains insufficiently studied (M aty as et al., 2009). This lack of knowledge precludes the implementation of policies and measures for climate change mitigation and adaptation. Besides, these species tend to occupy peripheral, arid areas that are subject to severe human impact. One such example is the circum-Mediterranean conifer Aleppo pine (Pinus halepensis). This thermophilic pine is the most widely distributed tree species in the Mediterranean basin. It shows an extensive ecological breadth and is seemingly adapted to a broad range of abiotic stressors and disturbances (Ne'eman et al., 2004;Schiller & Atzmon, 2009;Klein et al., 2013;Santini et al., 2019). Previous studies have shown that populations in xeric sites exhibit different characteristics than those in mesic areas, indicating intricate anatomical, morphological, and physiological adaptations of the species (Santos del Blanco et al., 2013;Voltas et al., 2015;David-Schwartz et al., 2016). However, a comprehensive analysis of the climate factors underlying intraspecific variation is currently lacking for P. halepensis, which exemplifies the significant gaps in our knowledge on the adaptive patterns of tree species and populations inhabiting drought-prone regions (M aty as et al., 2009).
Local adaptation, which reflects a genetic adjustment of individuals to a particular environmental niche (i.e. higher fitness of resident than of foreign genotypes; Williams (1966)), is widespread in forest trees (Alberto et al., 2013;Valladares et al., 2014;Peterson et al., 2019;S aenz-Romero et al., 2019). The first step towards detecting local adaptation is to find evidence that a trait is relevant for fitness, which can be suggested by linking trait variability of populations with the historical climate at their location of origin. This has been extensively shown for tree height and diameter in the case of P. halepensis (Climent et al., 2008;Voltas et al., 2008;Esteban et al., 2010), which suggests ecotypic specialization for either resource-rich habitats favoring vegetative growth, or resource-poor habitats favoring early reproduction. However, local adaptation is challenging to characterize, because the differential allocations of resources to reproduction, storage, defence, and above-and below-ground growth throughout ontogeny trade-off for maximizing fitness, as evidenced for P. halepensis (Sbay & Zas, 2018;Santini et al., 2019). Despite these difficulties, height growth is the most widely used fitness-related trait to date for understanding adaptation in forest trees (Peterson et al., 2019;Fr ejaville et al., 2020), mainly because of its ample availability in extensive networks of common gardens. In particular, tree height can be considered a potential proxy for fitness in P. halepensis because this species is shade-intolerant, thus showing light-dependent growth patterns, and presents phenotypic plasticity subject to resource availability (Zavala et al., 2011).
The analysis and interpretation of intraspecific adaptive patterns can be accomplished using common-garden experiments and thorough statistical modeling of genotype by environment (G 9 E) interactions (M aty as, 1996;Rehfeldt et al., 2002). Within this context, the analysis of climate-response-functions, which describe the relationships between the performance of a population across test sites and the climate of those sites (Rehfeldt et al., 2002), can be useful in forecasting the effects of climate change on tree species (Peterson et al., 2019). Thus, the response of populations to the new environments at the test sites can be interpreted as a simulation of their response to climate change. In this study, climate-response-functions were obtained using a mixed factorial regression model approach. This approach is often used in plant breeding given the relevance of G 9 E for the optimal deployment of new releases (van Eeuwijk et al., 2005;Malosetti et al., 2013;Heslot et al., 2014), but is less frequently applied in forest genetics (Sixto et al., 2016). Factorial regression describes G 9 E interactions as differential genotypic sensitivities to specific environmental factors (i.e. differences in the slope of response). This methodology has two main advantages in gaining insight into the environmental causes underlying adaptation. First, it allows for straightforward statistical testing of the drivers and relevance of differential plasticity linked to sizeable G 9 E interactions. Second, it is flexible in accounting for population structure (i.e. ecotypic variation) in a multi-environment setting, which is valuable for highly unbalanced datasets denoting incomplete replication of populations across testing sites.
We hypothesized that the combined responses to temperature and precipitation have shaped the intraspecific adaptive strategies of the widespread thermophilic pine, P. halepensis. We used tree height to unravel the potential of populations to perform across a range-wide environmental gradient. Out of the components of above-ground growth, tree height is under the strongest genetic control and, as such, is commonly used as an indicator of population and site differences in productivity (Skovsgaard & Vanclay, 2007;Moles et al., 2009). Our goal was to delineate ecotypic responses to climate and obtain projections of future growth under climate change, which are related to the carbon storage capacity of dry forests and their role as net carbon sinks (Zhu et al., 2018). Under climate fluctuations, the productivity of species should decrease from the 'leading-edge' (recently favorable habitats where populations are expanding) to the 'trailingedge' (where populations are contracting from recently unsuitable habitats), following a climate gradient.
Specifically, this study addresses the following questions: What is the relative importance of precipitation and temperature as factors that could explain adaptive variation in a thermophilic species subject to varying degrees of water shortage? Does accounting for G 9 E interaction improve the estimation accuracy of tree height for the current climate and modify the projections of tree height under climate change? And are ecotypes from dry environments more likely to decline in growth by climate warming than those occupying wetter areas, hence following the commonly assumed 'trailing-edge, leading-edge' range-shift pattern? To answer these questions, we analyzed a broad yet highly unbalanced trial network in which five major P. halepensis ecotypes were evaluated across the Mediterranean, including independent data validation (new population-trial combinations). Our study is likely to be the first range-wide examination of adaptive variation for an ecologically relevant forest species native to warm, dry environments which uses independent data to validate predictions. Our study further analyzes the observed ecotypic patterns to forecast the species' response to future conditions (period 2071-2100 AD).

Field trials and plant material
Nine Aleppo pine (Pinus halepensis Mill.) provenance trials (hereafter, training trials) were used to obtain intraspecific response functions to climate. The training trials were distributed across the Mediterranean basin (Israel, Italy, and Spain; Supporting Information Fig. S1), in areas representative of the climatic conditions where Aleppo pine is found naturally. In particular, MAT (mean annual temperature) ranged from 11.1 to 19.5°C, and MAP (mean annual precipitation) from 291 to 757 mm (highresolution WorldClim database v.1.4; Hijmans et al., 2005). The climate of the training trials (Table S1) closely matched the species' core climate envelope, with an average difference between the midpoint of training trials and the midpoint of the species distribution range of +1.3°C (MAT) and À55 mm (MAP). The training trials (Table S2) had been previously used to characterize the extent of ecotypic differentiation in tree height plasticity and survival for the species, with no significant ecotype by trial interactions for survival (Voltas et al., 2018). Five additional trials with climate characteristics within the range of those of the training dataset were used as test trials to evaluate the predictive ability of response functions (see Tables S1, S2) and were from Australia (Spencer, 1985), Italy, Morocco (Sbay & Zas, 2018) and Spain (Santos del Blanco et al., 2013).
The training trials evaluated seed sources from 82 populations, representing most of the species' natural distribution across the Mediterranean basin. The test trials evaluated seed sources from 142 populations, with 67 populations being also present in the training trials ( Fig. S1; Table S3). Both datasets were extremely unbalanced (Methods S1). A total of 157 different populations were considered. For each population, climatic variables at their origin were obtained for the period 1961-1990 (Climatic Research Unit, CRU TS 3.24; Harris et al., 2014), at a coarse resolution (0.5°9 0.5°grid-box basis) following Voltas et al. (2018), and were retrieved from each population's nearest grid. The climate variables were selected based on previous identification of climate factors underlying ecotypic variation (Tapias et al., 2004;Voltas et al., 2008). These included mean annual temperature (MAT o ), the maximum temperature of the warmest month (i.e. July; TMX o ), minimum temperature of the coldest month (i.e. January; TMN o ), the temperature annual range (i.e. TMX o -TMN o ; TAR o ), mean annual precipitation (MAP o ), and mean summer precipitation (MSP o ) at the origin.
Using MAT o , MAP o , TAR o and MSP o , we obtained linear discriminant functions of the five ecotypes identified through clustering. The discriminant functions were employed to allocate each of the 75 new populations of the test dataset to a particular ecotype and to obtain subsequent probabilities of membership (Table S6).

Multi-environment trial analysis and factorial regression modeling
The nine training trials were combined and analyzed as multi-environment mixed model analysis of variance for population-trial means with the fixed trial (T), ecotype (E) and trial by ecotype interaction effects, and random effects for the population (P) nested to ecotype and population nested to ecotype by trial interaction (random terms are underlined): where Y ijk is the observation (tree height) of the j th population of the k th ecotype in the i th trial, l is the general mean, T i is the fixed effect of the i th trial, E k is the fixed effect of the k th ecotype, P(E) jk is the random effect of the j th population nested to the k th ecotype, (TE) ik is the fixed effect of interaction between the k th ecotype and the i th trial, and TP(E) ijk is the random (residual) effect of interaction between the j th population nested to the k th ecotype and the i th trial. For the explanation of ecotype by trial effects, we used a factorial regression modeling approach by which regression terms for the interpretation of genotype by environment (G 9 E) interaction are incorporated in the form of environmental covariates to the levels of the environmental factor. Factorial regression helps to detect genotypes that are differentially sensitive to changes in identified environmental components (van Eeuwijk, 1995;Malosetti et al., 2013).
We fitted mixed factorial regression models (Denis et al., 1997) to ecotype by trial effects in the model described in Eqn 1 (model 1, hereafter), searching for climate variables underlying differential plastic responses at the ecotype level (Voltas et al., 2018). In factorial regression, an explicit environmental covariable Z can be included in the levels of the trial factor to describe the trial main effect plus the interaction term of interest, (TE) ik : where a stands for the general sensitivity of all populations to the environmental variable Z i , b k stands for the specific sensitivity of the k th ecotype to the same variable, and T i and (TE) ik are the residual effects of the i th trial and the interaction between the k th ecotype and the i th trial, respectively. The model described in Eqn 2 (model 2) allows identification of a climate variable that can explain absolute changes in tree height among trials, differential ecotypic reactions across trials, or both simultaneously. This model can be extended to the multiple covariables case: where Z il denotes the values for l climate variables underlying variation in T i , (TE) ik , or both. The terms a l Z il and b kl Z il are fixed terms because Z il values contain known, explicit environmental information.

Characterization of ecotypic responses to climate: response surfaces
We tested, for the selection of climatic factors in the model described in Eqn 3 (model 3), all possible combinations of one temperature-and one precipitation-related variable because temperature (or precipitation) variables were significantly correlated and, hence, partly redundant for explaining differential ecotypic performances: MAT and MAP, MAT and MSP, TAR and MAP, and TAR and MSP. For each combination, a quadratic response surface was incorporated into the levels of the trial main effect and the ecotype by trial interaction simultaneously in model 3, as follows: Eqn 5 where Z il denotes the values for l = 2 climate variables of temperature and precipitation at the trial i. We refer to this model as a factorial regression response surface model (hereafter FRSM). The heterogeneity in the b k 's for Z 1 and Z 2 accounts for the interaction between ecotypes and trials, while the sum of multiplicative terms P 2 l ¼1 b kl Z il approximates this interaction (van Eeuwijk, 1995). To facilitate interpretation of the b k 's, the climate variables were standardized to a mean value of 0 and a standard deviation of 1. For the selection of the best combination of climate variables in the FRSM, we used the marginal R 2 (Nakagawa & Schielzeth, 2013), which informed on the variance explained by the fixed (climate) factors alone, that is, the predictive ability of each model based on the (nonstandardized) climate variables, and omitting any trial-related effect from model 3. The data met the underlying assumptions of normality and homogeneity of error variances ( Fig. S3; Methods S2).
The best FRSM model included MAT and MAP as climate variables to the detriment of other combinations (Table S7). They provided noncorrelated information about the climatic characteristics of the training trials (r = À0.11; P = 0.78). Based on this model, response surfaces to climate were obtained at the ecotype level, and for the whole dataset (i.e. general performance model; see Methods, S2).

Future forecasts of climate responses
For the current distribution area of P. halepensis (EUFORGEN, 2009), height growth models were projected at 30 arcsecond resolution (c. 1 km 2 ) using climatic maps for 'reference' conditions (i.e. the period 1961-1990 AD) from WorldClim (Hijmans et al., 2005) and for 'future' conditions (i.e. the period 2071-2100 AD). For future conditions, we used the climate projections from S aenz-Romero et al. (2017) under two representative concentration pathways: a 'moderate' scenario (RCP 4.5 W m À2 ; hereafter RCP 4.5) and a 'pessimistic' scenario (RCP 8.5 W m À2 ; hereafter RCP 8.5). Based on the general performance model (i.e. disregarding ecotypic variation), and also for each of the ecotypic models independently, we projected tree height for every pixel of the current distribution of P. halepensis. We then produced maps of height differences between reference conditions and height projections, as a percentage of reference height of each pixel (c. 306 000 occurrence points; EUFORGEN, 2009). Maps were generated under the RCP 4.5 and 8.5 scenarios, and for the general performance and each ecotype model (see Methods S3). To delineate the spatial patterns of ecotypic performance, we compared the relevance of differential plasticity on tree height. We first assumed that each ecotype could be potentially distributed over the entire distribution range of the species (hereafter, regional distribution). In this way, we made the comparison easier among ecotypes, such that the implications of these results could be better assessed for management purposes (e.g. deployment or assisted migration strategies).
We then produced a map of the 'tallest ecotype' by projecting the ecotype with the highest height for every pixel of the distribution under each scenario. This map was filtered at the pixel level for tree height of the tallest ecotype, which exceeded the secondbest ecotype by at least 0.26 m (least significant difference of ecotype means at P = 0.05 in model 1). Moreover, we created histograms to summarize the height frequency distribution under reference and future climates for general performance and each ecotype separately and independently for both regional distribution (i.e. maximum occupancy) and local conditions (i.e. local occupancy). Nonparametric statistics were used to assess changes in frequency distribution (local vs regional, reference vs future climate; Methods S3).

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Results
Evidence of differential ecotypic plasticity linked to climate The thermophilic P. halepensis showed phenotypic plasticity for the total height that was genetically structured, as revealed by significant ecotype-by-environment interactions (likelihood-ratio test = 15, P < 0.001, for models with and without ecotype by trial interaction term). According to the selected FRSM, the nature and magnitude of differential tree height responses were adequately described by simple climate factors, MAT and MAP in particular (marginal R 2 = 43.1%; Table S7). For reference climate conditions (period 1961(period -1990, and ignoring ecotypic variation, tree height at age 15 yr exhibited a 2.5-fold increase with increased precipitation and temperature (from 4 to 10 m approximately), and the largest tree height was mainly found for locations with MAP > 75% of potential evapotranspiration (Fig. 1a). FRSM, however, pointed to variable ecotypic responses to precipitation (MAP, MAP 9 MAP) and its combined effect with temperature (MAT 9 MAP) (Table S8), and allowed quantification of specific sensitivities to climate for each ecotype (Table S9). These sensitivities produced height variability underlying G 9 E interactions (i.e. differential ecotypic plasticity) and also ecotypic differentiation in mean height across trials (i.e. standing ecotypic variation, with up to 13% relative difference between extreme ecotypes; Table 1). The FRSM-derived ecotypic equations (Table 1), which are solved in Fig. 1 for the climate envelope defined by the training trials, had, in most cases, coefficients of determination equal to or higher than 0.75, indicating a reasonable fit. The relationship between tree height and climate, therefore, varied if ecotypes were considered ( Fig. 1b-f). The standard deviations of predicted tree height, obtained through bootstrap of regression estimates, are shown in Fig. 2. Most ecotype predictions had deviations below 0.3 m across MAT and MAP combinations, but the weakest prediction power corresponded to ecotype DHT, followed by WHC. Also, the weakest spots for prediction purposes across ecotypes corresponded to the four combinations of extreme MAT and MAP values, especially those of high MAP (Fig. 2).
The MAP variation across trials (from 300 to 800 mm) underlay varying ecotypic height performances. The Wet-summer/sub-Humid/Cool (WHC) ecotype (mainly from northeast Spain and southeast France) displayed the sharpest dependence on precipitation, as indicated by significant above-(linear) and below-average (quadratic) sensitivities to MAP (Table S9). Conversely, the Dry-summer/Semiarid/Temperate (DST) ecotype (from coastal southeast Spain, inland Algeria and Tunisia, and Israel) and the Dry-summer/sub-Humid/Temperate (DHT) ecotype (from Greece and coastal Tunisia) showed a subtle sensitivity to precipitation, for either linear (DST, DHT) or quadratic (DHT) sensitivities (Table S9). Contrary to precipitation, responses to temperature were approximately uniform across ecotypes (Table S8) and showed a positive reaction to warming for the whole species (Table 1). In some instances, however, the abovementioned responses to MAP were contingent on the simultaneous increase in MAT and MAP (Table S9), which strongly benefitted DHT while penalizing WHC under the same conditions. Based on the ecotypic equations, the combined sensitivities to climate factors translated into height differences between extreme ecotypes (DHT and DSC) of 0.8 m at the shortest trial (Ademuz, Spain) and 2 m at the tallest trial (Oria, Italy) (Table S1).
Cross-validation on the training dataset (Table 1) indicated that the performance of ecotypic models was reasonably good, with mean absolute errors of c. 1 m (height range of populationtrial combinations = 10.2 m) and with correlations equal to or higher than 0.75 between observed and predicted values. Using independent validation, the predictive ability of ecotypic models was superior (WST) or similar (DSC, DST, WHC) to that of the general performance model (ASE statistic; Table 1). The only exception was DHT, whose heights were most often underestimated by both the ecotypic and general models (Table S10). Overall, these results demonstrated the good predictive ability of climate functions.

Description of ecotypic responses to climate
The optimum height trajectories for the current average climate of each ecotype (curved black arrows in Fig. 1) indicated the direction and strength of maximum response and provided information on differential plasticity. In particular, DST responded positively to MAT up to a mean temperature of c. 17°C, beyond which it showed a mild sensitivity to MAP (Fig. 1b). By contrast, the Dry-summer/Semiarid/Cold (DSC) ecotype (from inland Spain and Morocco) showed positive responses to increasing temperature (up to 1-2°C beyond their current mean climate) when accompanied by increasing precipitation (Fig. 1c). DHT showed positive, strong responses to temperature and precipitation (Fig. 1d), attaining the maximum height of all ecotypes within its climate envelope (%12 m). Finally, the Wet-summer/ Semiarid/Temperate (WST) ecotype (from the Balearic Islands, Italy, and eastern Spain) and the WHC ecotypewith the coldest and wettest originmildly benefitted from increased temperature and precipitation (Fig. 1e,f). Most ecotypes displayed saddle-shaped response surfaces (Fig. 1), with no obvious combination of temperature and precipitation being optimum for tree height (Table S9). The exception was WHC, which was projected to reach a climate optimum at 1500 mm MAPwell above its current local climate in southeast France and northeast Spain (Table S9).

Spatial patterns of current performance
The climate scope of populations at their origin matched the species' current climate envelope well, as shown in Fig. 3(a) (left panel, average climate differences between midpoint of populations and midpoint of species distribution range: MAT = +1.3°C, MAP = À62 mm). Comparing the differential plasticity of ecotypes to climate, DHT and WHC should attain, at present, the highest and lowest median heights, respectively, for the entire region (M R = 6.82 m and 6.58 m for DHT and WHC; Fig. 4). These results resembled the ecotypic ranking across trials (Table 1), which supported the effectiveness of the training trials for modeling purposes. The results also showed that that DHT and WHC should perform best, at present, in extensive areas of the Mediterranean basin (Fig. 3a, central panel). We observed much larger median differences among ecotypes locally, that is, considering only occurrence of each ecotype at its present location. These differences amounted to c. 3 m (M L values, Fig. 4).   The number of population-trial combinations per ecotype (N), estimated ecotype means across trials (and their standard errors), full quadratic models (response surfaces) at the whole species (general performance) and ecotype levels as a function of mean annual temperature (MAT) and mean annual precipitation (MAP), and predictive performance of ecotype models derived from a factorial regression response surface model for tree height at age 15 yr. The predictive performance of ecotype models is evaluated by k-fold cross-validation based on k = 10 replicate samples (training dataset) and by the predictive ability of the response surfaces using independent data (test dataset). For cross-validation, the Mean Absolute Error (MAE) and the correlation (r) between observed and predicted values are shown. For predictive ability using independent data, the Average Squared Error (ASE, the sum of squared errors divided by the number of cases, N) of the corresponding ecotype model and the correlation (r) between observed and predicted values using the entire test dataset and a simplified dataset containing the subset of populations also evaluated in the training dataset (r*) are shown. For comparison with ecotypic models, ASE values for the training and test datasets using the general performance model are also presented. In the table, E k is the fixed effect of the k th ecotype, a (a') is the linear (quadratic) sensitivity of the whole species to a given environmental variable (or to the cross-product of variables MAT and MAP), and b (b 0 ) is the ecotype-specific linear (quadratic) sensitivity for the same variable (or cross-product of two variables MAT and MAP).

New Phytologist
Indication of local adaptation was limited to the mesic ecotypes DHT and WHC, which often showed a significantly greater height at origin (> 40% of their local area) when compared with the other ecotypes (Fig. 3a, right panel).

Spatial forecasts under future conditions
Projections for the period 2071-2100 AD indicated a widespread increase in tree height across the species' current distribution ( Fig. 5; more details in Figs S4-S7; Table S11). This increase may occur despite climate change is projected to push most ecotypes beyond their current temperature range for both RCP 4.5 (1.8-3.5°C warmer and 5-111 mm drier on average) and 8.5 (3.6-5.7°C warmer and 73-188 mm drier on average) (Figs. 3b,c, left panels, and Table S12). However, wet summer ecotypes were the exception, for which growth inhibition was often projected both regionally (WHC) and locally (WST and WHC). New Phytologist (2020) Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com

New Phytologist
For RCP 4.5, the general performance model, disregarding ecotypic variation, projected height enhancement in 51% of the current species' distribution (northern, northwestern, and southeastern Mediterranean, and continental areas of northern Africa; Fig. 5a). However, notable differences were observed among the ecotypes (Fig. 5). Regionally, the largest area of height enhancement (65%) was projected for DHT, followed by DST (59%). Inhibition at over 40% of the distribution was projected for WHC (see also Fig. S6a). Overall, height inhibition of up to 1.4 m and enhancement of up to 2 m was expected (5% and 95% quantiles of height differences, respectively, for general performance; Figs S6a, S7a). Locally, dry summer ecotypes (DST, DSC, and DHT) were projected to respond positively to warming, whereas wet summer ecotypes (WST and WHC) were mainly projected to experience stability or inhibition (Fig. 5). These responses translated into significant height differences among ecotypes both regionally (difference in median values between extreme ecotypes = 1.4 m) and locally (Fig. S4). Locally, however, differences were reduced from c. 3 m (for reference conditions) to 1.7 m, and the ecotypes showing superiority at origin changed compared with reference conditions, with the growth of DHT being favored and WHC inhibited under moderate warming (Fig. 3b, right panel).
For RCP 8.5, the general performance model projected height enhancement for a slightly smaller fraction of the current distribution (48%), compared to RCP 4.5 (Fig. 5b). Overall, height , spatial projections of tallest ecotypes for tree height (m) at age 15 yr across the entire species' distribution range in the Mediterranean basin, considering maximum occupancy of each ecotype (center), and percentage of the area where each ecotype shows height superiority (P < 0.05) at the origin, that is, under local occupancy (right). (a) Reference conditions (period 1961(period -1990 based on WorldClim data (Hijmans et al., 2005). (b) Future conditions (period 2071-2100 AD) using the climate layers for RCP 4.5. (c) Future conditions using the climate layers for RCP 8.5. In left panels, trials are depicted as black triangles (upward-pointing triangles, training trials; downward-pointing triangles, test trials) and populations as colored squares (populations and trials represented in (b) and (c) are based on reference climate). In the central panels, the pie charts represent the frequency of each ecotype significantly outperforming (P < 0.05) other ecotypes for each climate condition considering maximum occupancy. The species range is derived from the EUFORGEN distribution map (EUFORGEN, 2009) for a resolution of 10 0 . DHT, Dry-summer/sub-Humid/Temperate ecotype; DSC, Dry-summer/Semiarid/Cold ecotype; DST, Dry-summer/Semiarid/Temperate ecotype; WHC, Wet-summer/sub-Humid/Cool ecotype; WST, Wet-summer/Semiarid/Temperate ecotype.
inhibition of up to 2.8 m and enhancement of up to 2.5 m was expected (5% and 95% quantiles of height differences, respectively, for general performance; Figs S6b, S7b). The dry summer DST and DSC ecotypes still showed superior growth compared to growth under reference conditions (Fig. 5), both regionally and locally. Conversely, DHT and WST ecotypes mainly experienced enhancement regionally, and WHC largely exhibited both regional and local inhibition (Fig. 5b). Only the dry summer DSC and DHT ecotypes exhibited superiority at origin (Fig. 3c,  right panel). These contrasting responses exacerbated the height differences among ecotypes projected for RCP 4.5 regionally (difference in median values between extreme ecotypes = 2.3 m), but not locally (difference = 0.7 m; Fig. S5). While DST, DSC, and WST showed a significantly higher median height regionally than locally under both climate scenarios, the reverse situation was observed for ecotypes DHT and WHC (Figs S4, S5).
From these results, it could be concluded that climate change would lead to smaller height differences among ecotypes, assuming there is local occupancy (i.e. no dispersal), especially for the most stringent RCP 8.5 scenario. This was mainly because the dry summer ecotypes benefitted from warming at their origin. Conversely, height differences were amplified under maximum occupancy. Accordingly, the climate models forecasted that WHC would be outperformed by DHT and also by DSC, across the whole distribution range (Fig. 3b,c).

Discussion
We developed, for the widespread thermophilic conifer P. halepensis, ecotype-specific vertical growth responses to climate which were then projected for reference and future climate conditions (2071-2100 AD) across the Mediterranean. Tree height inferences based on a factorial regression response surface model (FRSM) revealed good predictive ability of ecotypic climate functions. Importantly, the array of provenance trials used closely approximated the species' climate envelope, which is critical for (b-f) ecotype level, considering maximum occupancy in the Mediterranean basin. The histograms depict the class frequencies (%) of tree height at the ecotype and general levels for the entire species' distributions. The class frequencies of tree height at the geographic origin of each ecotype, estimated for a 25 km radius around each ecotype's populations, are also shown in the histograms as hatched bars. A Median test (Z) for comparison between local (L) and regional (R) performance, together with their median values (M, in m), is included. Positive Z values indicate that local level is superior to regional level performance.

Research
New Phytologist establishing realistic projections (Wang et al., 2006). Additionally, the climate-based ecotypic models were validated using an independent dataset, which reduced the uncertainty in our projections of climate responses for the species. Our approach, based on ecotypic classification, can provide parsimonious and, thus, hands-on information on the adaptive structure of a species to forest managers and end-users of genetic material. Other methods of investigating intraspecific variability may cluster populations based directly on their trait responses (DeLacy et al., 1996), or fit 'universal response functions' as an interaction between the trial's climate and the climate of a population's origin (Wang et al., 2010). However, these methods are sensitive to unbalanced data and may require intricate procedures for estimating population performance at sites with extreme climate (Wang et al., 2010).
Despite the importance of gene flow in P. halepensis (Steinitz et al., 2012), the FRSM model indicated that differentiation by climate is strong throughout the entire species' range (Aitken & Bemmels, 2015). Simple climate factors such as mean annual temperature and precipitation adequately described the nature and magnitude of intraspecific vertical growth responses. Therefore, when ecotypes are considered, the relationship between tree height and climate varies considerably. However, evidence of local adaptation was limited to mesic ecotypes. For the remaining ecotypes, the absence of local adaptation for tree height could be related to complex trade-offs among traits maximizing fitness for this species (Santini et al., 2019), or to the existence of adaptation lags related to historical climate (e.g. Rehfeldt et al., 2002). Importantly, the climate-based grouping of ecotypes matched faithfully a recent definition of genetic groups of the species based on molecular (putatively neutral) SNP information (Serra-Varela et al., 2017), which identified populations with similar evolutionary histories. The observed match between classifications is supportive (but not conclusive) of the impact of selection processes under different environmental conditions, as some SNP markers used by Serra-Varela et al. (2017) were likely influenced by adaptive selection. As such, our results suggest that selection by climate may have favored different genotypes in populations subjected to distinctive bottlenecks during the post-glacial westward Mediterranean expansion of P. halepensis (Grivet et al., 2009).
Indeed, the relevance of climate in shaping adaptive patterns for P. halepensis has been well recognized for many morphophysiological traits, including growth and reproduction (e.g. Calamassi et al., 2001;Climent et al., 2008;Schiller & Atzmon, 2009;Klein et al., 2013;Voltas et al., 2015;David-Schwartz et al., 2016). Altogether, our results support rangewide G 9 E regarding climate in P. halepensis (Savolainen et al., 2007;Peterson et al., 2019). This thermophilic pine is known to exhibit a general adaptive syndrome by which less plastic populations are characterized by higher survival and larger investment in reproduction and reserves than in growth (Climent et al., 2008;Voltas et al., 2018;Santini et al., 2019). In particular, intraspecific variation in tree height was mainly associated with precipitation, suggesting that variation in natural selection is more related to water availability than to temperature. However, this is not unexpected for a conifer subject to steep gradients of water availability (Richardson et al., 2014). Therefore, the observed height responses, deduced from ecotype-specific climate sensitivities, are fundamental for interpreting intraspecific adaptation (Benito-Garz on & Fern andez-Manjarr es, 2015). and reference conditions (period 1961(period -1990 across the distribution range of Pinus halepensis. Right panels. Pie charts summarizing the relative importance of height differences (enhancement (> 10%), inhibition (< À10%), stability (À10%, +10%)) at the local (charts on the left) and regional (i.e. maximum occupancy of each ecotype, charts on the right) levels for general performance and for each ecotype. DHT, Dry-summer/sub-Humid/Temperate ecotype; DSC, Dry-summer/Semiarid/Cold ecotype; DST, Dry-summer/Semiarid/Temperate ecotype; WHC, Wet-summer/sub-Humid/Cool ecotype; WST, Wet-summer/Semiarid/Temperate ecotype. See Supporting Information Fig. S7 for absolute tree height differences.
The projected climate responses presented an evident northeast-southwest cline, with mesic ecotypes outperforming ecotypes of drier climates across a large proportion of the species' distribution area. In particular, it indicated that ecotypes WHC and DHT should grow best, at present, in extensive areas of the Mediterranean basin. These ecotypes from northern and eastern Mediterranean areas allocate fewer resources to reproduction (Climent et al., 2008) and transpire more in summer (Santini et al., 2019), but they do not show reduced survival than their drier counterparts from southern and western regions (Voltas et al., 2018). In particular, DHT (native to the eastern Mediterranean basin) exhibits high genetic diversity and, importantly, encompasses most of the functional variation of P. halepensis (Grivet et al., 2009). Along with this suite of features, the high degree of climate-related height plasticity of DHT suggests that this ecotype may be decisive for the persistence of species through climate change and, thus, important for the implementation of deployment or assisted migration strategies. Although taller trees may be more vulnerable to conduction-blocking embolisms (Olson et al., 2018) and, therefore, to warming-induced drought stress, David-Schwartz et al. (2016) have shown that eastern populations of the species are not more embolism-susceptible than their western counterparts (David-Schwartz et al., 2016). The latter suggests that fast vertical growth and drought tolerance are compatible within P. halepensis.
The larger height differences observed among ecotypes locally (relative to range-wide differences) could mainly be a consequence of more favorable conditions occurring at the wettest spots of the species' range for a conifer limited by drought stress (S aenz-Romero et al., 2019). However, those height differences at the local level also reflected differential ecotypic sensitivities to climate, as revealed in this study. Notably, the lack of superiority at the origin for most ecotypes suggests that differential plasticity is not necessarily locally adaptive, at least for vertical growth, despite the existence of intraspecific plastic variation for tree species with long lifespans (Benito-Garz on et al., 2019;G arate-Escamilla et al., 2019). In this regard, while specialization in harsh environments is likely linked to limited plasticity and a conservative resource-use strategy, the opposite is to be expected for specialization in favorable environments (Valladares et al., 2007).
Under future climate (i.e. 2071-2100 AD), and despite the projected warmer and drier conditions, projections indicated widespread height enhancement of P. halepensis, but smaller height differences among ecotypes assuming there is local occupancy (i.e. no dispersal), especially for the most stringent RCP 8.5 scenario. Height enhancement becomes particularly evident for the most xeric populations of the western Mediterranean, which are forecasted to be favored even under extreme warming. On the other hand, the projected height inhibition of wet summer ecotypes, either regionally (WHC) or locally (WST and WHC), suggests that these ecotypes could be at risk. The future height decline of wet summer ecotypes may lead to forest fragmentation as a result of maladaptation to rising drought (Montw e et al., 2016). Importantly, forecasted trends for RCP 4.5 and RCP 8.5 were similar, but with more inhibition under the more severe RCP 8.5 for wet summer ecotypes. Nevertheless, our dataset lacked test trials at sites warmer than 20°C MAT, which represents 17% and 37% of the distribution of the species under RCP 4.5 and 8.5, respectively. Thus, the reliability of projections beyond this temperature threshold might be compromised.
Prevalent positive responses to warming are supported by studies reporting radial growth stimulation of old Aleppo pine forests even in semi-arid/arid and warm environments such as the pre-Saharan Algerian steppe (Choury et al., 2017) or Israel (Dorman et al., 2015). These responses could involve a shift of the actual growing season towards earlier months in winter-spring (Bigelow et al., 2014). This shift would make growth more reliant on water availability during the very early growing season (Choury et al., 2017), a period progressively more suitable for tree performance as the climate becomes warmer and drier. Also, elevated CO 2 (eCO 2 ) stimulation effects on growth cannot be discarded, although uniformity of ecotypic height responses to eCO 2 can be expected (Resco de Dios et al., 2016). Anyhow, tree height responses to eCO 2 are indirectly accounted for in our projections (see Notes S1 on eCO 2 effects). Other studies, however, have reported drought-induced declines in radial growth for P. halepensis under severe dry spells (de Luis et al., 2013;Gazol et al., 2017). In this regard, weather extremes can have a substantial impact on Aleppo pine performance (Santos del Blanco et al., 2013), which could only be approximated by the climate averages used in our study. The sequence of consecutive extreme weather events and, also, the dynamics of biotic damage will, almost certainly, tip the balance towards survival or mortality (Sang€ uesa- Barreda et al., 2015), as future conditions will likely produce a mismatch between the species' fitness limits and its current distribution. On the other hand, our projections agree with alternative method-testing scenarios for assisted migration in P. halepensis, in which increased volume is projected for 2050 AD, except in southeast France and coastal Spain (Benito-Garz on & Fern andez-Manjarr es, 2015).
In our approach, we substituted space for time in forecasting height growth responses to climate (but see Meril€ a & Hendry, 2014, for potential problems). Nevertheless, the evolutionary history of P. halepensis also involves fitness traits other than height, as well as trade-offs mainly related to drought and fire events (Aitken & Whitlock, 2013;Pausas, 2015;Santini et al., 2019), that could be further used to describe adaptive patterns. Alternative methodologies that remain to be tested could involve treering records and FRSMs by untangling the effects of the site (space) and year (time) in the explanation of differential ringwidth responses in trial networks. Indeed, dendrochronological approaches applied to common-garden experiments are emerging that aim to forecast responses to climate fluctuations (Montw e et al., 2015;Suvanto et al., 2016;George et al., 2019). Also, the combination of height records with radial growth modeling could improve estimates of productivity and carbon sequestration potential (Klesse et al., 2018). For a holistic perspective of our short-term projections, future studies should compare these forecasts to long-term dispersal rates and gene flows, which are considered to be high in P. halepensis (Steinitz et al., 2012;Serra-Varela et al., 2017).
This study demonstrates that phenotypic changes in tree height, for a widespread and ecologically diverse conifer such as P. halepensis, are, to a large extent, genetically based and not exclusively the result of uniform plastic responses. We show that this variation is related to climate and precipitation in particular, and indeed, accounting for intraspecific sensitivity to climate is a crucial step for forecasting tree responses (Peterson et al., 2019). The latter is exemplified here for P. halepensis, a species native to warm, dry environments of substantial ecological, cultural and economic importance which may have to confront very different conditions in the near future. Changes in climate may, therefore, drive nonintuitive shifts in range, particularly affecting the edge of the species' range in wetter regions. We do not yet know to what extent plasticity of the integrated phenotype (or the exhibition of particular arrays of traits in a changing environment) may be crucial for the resilience and persistence of forests in the near future. Nevertheless, the patterns we report here are crucial in moving towards a comprehensive understanding of prevalent climatedriven trade-offs. Further insights are urgently needed to conclusively predict the evolutionary response of this exemplary Mediterranean conifer to climate change and the reactions of forest trees typical of dry environments in general.

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article.     Spatial projections (future RCP 8.5 scenario) for tree height (m) at age 15 yr at the species (general performance) and ecotype level for Pinus halepensis in the Mediterranean basin assuming maximum occupancy.

Fig. S6
Cumulative distribution (frequency classes) and density plots of the predicted height differences at age 15 yr between future and reference climate conditions. Methods S1 Study sites and populations assignment to ecotypes.
Methods S2 Ridge analysis and validation of response surface functions.
Methods S3 Quantification of ecotypic responses to future climate.
Table S1 Geographic and climatic characteristics of the nine trials used as the training dataset and of the five trials used as the test dataset.

Table S2
Mean survival, number of populations per ecotype at each trial, characteristics of the experimental design, and layout of the nine trials used as the training dataset and of the five trials used as the test dataset.

Table S3
Geographical and climatic (CRU TS 3.24, 1961-1990 period; 4) characteristics of the 82 populations of Pinus halepensis used as the training dataset for factorial regression modeling.

Table S4
Main climate characteristics of ecotypes (AE standard deviation of population means) based on the WorldClim dataset.

Table S5
Information criteria used to select the best classification of populations into ecotypes for tree height at age 15 yr.

Table S6
Probability of assignment to each ecotype of 44 new populations evaluated in test trials using linear discriminant functions.

Table S7
Factorial regression response surface models of tree height at age 15 yr partitioning the ecotype 9 trial interaction term.

Table S8
Full factorial regression response surface model (FRSM) of tree height at age 15 yr using mean annual temperature (MAT) and mean annual precipitation (MAP) at the trial level.

Table S9
Ecotypic sensitivities (b MAP , b MAT9MAP , b 0 MAP ) to selected response surface terms of Eqn 5 from the main text, related to mean annual temperature (MAT) and mean annual precipitation (MAP) for tree height (m) at age 15 yr.
Table S10 Predictive performance of quadratic models at the ecotype level for tree height at age 15 yr of Pinus halepensis based on mean annual temperature and mean annual precipitation.
Table S11 Comparison of frequency distributions of tree height between reference and future climate conditions (RCP 4.5 and RCP 8.5) at the species (general performance) and ecotype levels.
Table S12 Summary of climatic variables (MAT = mean annual temperature; MAP = mean annual precipitation) within the distribution of Pinus halepensis for reference and predicted future conditions under future RCP 4.5 and RCP 8.5 scenarios.
Please note: Wiley Blackwell are not responsible for the content or functionality of any Supporting Information supplied by the authors. Any queries (other than missing material) should be directed to the New Phytologist Central Office.
New Phytologist is an electronic (online-only) journal owned by the New Phytologist Trust, a not-for-profit organization dedicated to the promotion of plant science, facilitating projects from symposia to free access for our Tansley reviews and Tansley insights.
Regular papers, Letters, Research reviews, Rapid reports and both Modelling/Theory and Methods papers are encouraged. We are committed to rapid processing, from online submission through to publication 'as ready' via Early View -our average time to decision is <26 days. There are no page or colour charges and a PDF version will be provided for each article.
The journal is available online at Wiley Online Library. Visit www.newphytologist.com to search the articles and register for table of contents email alerts.
If you have any questions, do get in touch with Central Office (np-centraloffice@lancaster.ac.uk) or, if it is more convenient, our USA Office (np-usaoffice@lancaster.ac.uk) For submission instructions, subscription and all the latest information visit www.newphytologist.com New Phytologist (2020) Ó 2020 The Authors New Phytologist Ó 2020 New Phytologist Trust www.newphytologist.com

New
Phytologist