Early Differentiation of the Phenotypic Space and Performance of Juniperus thurifera Across Woodland-Expanding Areas

Forest expansion after rural abandonment is changing landscapes, bringing new challenges and opportunities. Juvenile performance has an effect in the establishment of new forests and the subse-quent acclimation of adult trees. However, most studies have focused on the performance of adult trees. Here, we analysed differences in the phenotypic trait space of juveniles across gradients of Juniperus thurifera woodland expansion and explored which traits shape juveniles’ performance, seeking speciﬁc drivers of such differentiation and the role of the ontogeny. Thus, we surveyed functional leaf and plant-size traits in 180 juveniles growing in 18 plots at three different tree density stages (mature woodland, recently established area and a transition zone between both) and recorded local environmental conditions around each individual. We also estimated their radial growth and resilience to drought (in terms of resistance, recovery and resilience) as performance indicators. We found no relevant effects of the study environmental variables, but we did ﬁnd a shift in the phenotypic space among juveniles across the gradient, with taller crowns and higher amount of leaf N and d 13 C at the expanding edge. This trait differentiation across the gradient implied higher performance of young trees growing in old agricultural ﬁelds, both in growth rate and drought resilience terms. Our results showed a relevant land-use legacy on early trait differentiation that shapes more efﬁcient phenotypes in the expanding edge compared to those in mature woodlands, suggesting that the future of these new forests is not compromised under current conditions

local environmental conditions around each individual. We also estimated their radial growth and resilience to drought (in terms of resistance, recovery and resilience) as performance indicators. We found no relevant effects of the study environmental variables, but we did find a shift in the phenotypic space among juveniles across the gradient, with taller crowns and higher amount of leaf N and d 13 C at the expanding edge. This trait differentiation across the gradient implied higher performance of young trees growing in old agricultural fields, both in growth rate and drought resilience terms. Our results showed a relevant land-use legacy on early trait differentiation that shapes more efficient phenotypes in the expanding edge compared to those in mature woodlands, suggesting that the future of these new forests is not compromised under current conditions. Key words: forest expansion; functional traits; growth; junipers; land-use legacy; resilience.

INTRODUCTION
In the last decades, the expansion of woodlands over former agricultural fields associated to rural exodus is occurring pervasively in temperate Northern Hemisphere (Song and others 2018) and particularly in the Mediterranean areas (Chauchard and others 2007; Valladares and others 2014). Colonisation of abandoned fields by tree species requires the presence of pioneers, whose seedlings are able to thrive on open areas under intense radiation. In Mediterranean areas, increasing temperatures during the growing season frequently match with rainfall drop, which increase evaporation and reduce soil water availability, triggering summer droughts. Therefore, pioneer species in these areas do not only tolerate high irradiance, but also have to cope with the excess of intercepted solar radiation when carbon assimilation is limited by stomatal closure due to water stress (Joffre and others 2008).
Functional traits determine how individuals use resources from their environment (McGill and others 2006). Thus, regardless of the potential light/ shade tolerance of each species (Niinemets and Valladares 2006), plants are able to adjust their traits within a species-specific range in order to optimise the photosynthetic capacity in response to the available irradiance. This plastic response can therefore boost intra-specific and intra-individual variability within the canopy across horizontal and vertical gradients of light (Niinemets 2007;Valladares and Niinemets 2007). Among others, foliage acclimation to high light typically involves less nitrogen allocation to light harvesting, smaller leaf areas, increased photoprotective compounds, thicker mesophyll, better-developed sclerenchymatic tissues or denser stomata (Givnish 1988;Hanba and others 2002;Robakowski and others 2004). At the whole plant level, plant shape and size can also be modified as a function of the light environment, for example, plants are capable of orienting their leaves or modify the branching pattern or leaf arrangement (that is, phyllotaxis, angle, orientation), to reduce the fraction of leaf area directly exposed to the sun and protect the photosynthetic apparatus (Ellison and Niklas 1988;Stenberg and others 1998;Valladares and Niinemets 2007).
Plastic trait responses accommodate phenotypes to environmental heterogeneity, pursuing to keep up performance under different conditions (Pigliucci 2001). In fact, plant performance can be considered the association between phenotypes and the ecological success (Koehl 1996). In Mediterranean areas, performance assessed as growth rate and resilience to drought may greatly determine plants' chances to thrive (fitness) in these ecosystems, especially in juveniles that are more susceptible to extreme climatic events (Castro and others 2005;Herrero and Zamora 2014). Thus, the higher secondary growth saplings experience, the more likely those individuals overcome the first summer droughts, and consequently survive and recruit (Zamora and others 2001;Benavides and others 2015). However, relationships between traits and performance vary across species (Granda and others 2018), plant age (Gibert and others 2016) and resource availability gradients (Gimeno and others 2012a;Kimball and others 2012;Laughlin and others 2018). Despite emerging works delving into the ecological processes underlying the expansion of woodlands in Mediterranean ecosystems (Gimeno and others 2012a, b; Alfaro-Sá nchez and others 2019; Acuñ a-Míguez and others 2020; Martín-Foré s and others 2020; Villellas and others 2020), there are still important gaps about patterns and drivers of the phenotypic space and performance of juveniles across the colonising gradient (but see experiments developed by Gimeno and others (2012a, b) that need to be explored. Increasing our understanding about traits that mediate performance, particularly in relation to the response to droughts, will help us to predict their future chances in a global change scenario.
Here, we analysed the phenotypes, in terms of traits and performance, of Juniperus thurifera juveniles across gradients of woodland expansion, and we tested whether colonising phenotypes, those juveniles in the expanding edge (that is, recently established on former agricultural fields), differ from those growing in mature woodlands. Based on previous knowledge, same genetic pool prevails across such close areas (Villellas and others 2020), and therefore, phenotypic differences found among adult trees mainly derived from acclimation in response to different local environments. However, the understanding of this differentiation process, namely when it is expressed, the factors involved, and the consequences on individuals' performance, are still unsolved. We hypothesised that trait differences respond to contrasting conditions (for example, higher radiation, soil nutrients or less compacted soil in the expanding front determined by the land-use legacy) that young colonising individuals experience compared to those growing under canopy. This assumption implies an early trait differentiation within the individual development, opposed to a posterior differentiation expressed as individuals grow older and compete with conspecific trees for resources. We also expect that differences in traits would shape differential performance, being those with more competitive and acquisitive traits the trees that grow faster (Reich 2014) but accounted higher vulnerability to drought events (Martínez-Vilalta and others 2012; Serra-Maluquer and others 2018).

Sites, Sampling Design and Environmental Characterisation
We conducted the study in the Alto Tajo Natural Park and surroundings (Guadalajara, central Spain) (Figure 1). The climate is continental-Mediterranean with cold snowy winters and hot dry summers. The mean annual rainfall is 477.1 mm and the mean annual temperature is 10.4°C (Molina de Aragó n 40°50¢40¢¢ N, 1°53¢07¢¢ W, Spanish National Meteorology Agency). Soils are poorly developed on limestone and dolomite substrates and frequent outcrops. This area is largely covered by pines (Pinus nigra Arn. subsp. salzmannii (Dunal) Franco, which is replaced by P. sylvestris in colder areas), oaks (mainly Quercus ilex L. and Quercus faginea Lam.) and juniper woodlands (J. thurifera L.) which appear in the most continental parts (Ferrero and others 2006). Juniperus thurifera L., target species of this study, is a dioecious wind-pollinated tree, which forms low-density woodlands on poor, shallow, rocky soils. It tolerates a broad temperature range and can thrive on moderate to low precipitation typical of Mediterranean-climate regions (Montesinos Torres 2007).
We used 14 plots from a previous study (Villellas and others 2020) located in three sites close to towns Ribarredonda, Maranchó n and Huertahernando (Figure 1), where land abandonment has favoured the woodland expansion of J. thurifera over old agricultural fields. In particular, using satellite images they identified areas of well-conserved mature woodlands (MA), adjacent areas recently established (expanding edge or area, EA) and transition zones (TR) between two former stages. We completed the design in order to have two plots in each stage within the gradient (EA, TR, MA) and location, 18 plots in total (3 sites 9 3 stages 9 2 plots), following the same methods described in Villellas and others (2020). In particular, we used satellite images to identify the required areas (EA in Maranchó n, MA in Huertahernando and TR in Ribarredonda) and, later in the field, we framed the new plots considering the density of adult trees, that is, including at least 50 adult trees ( Figure 1; Supplementary Information Table S1).
In April 2018, we randomly selected and georeferenced 10 juveniles of J. thurifera in each plot among the censused juveniles in a preliminary characterisation, totalling 180 juveniles, 60 per stage (EA, TR, MA). We surveyed recruited juveniles over 50 cm height (to avoid small, vulnerable saplings) and diameter at breast height under 3 cm, considered the minimum size for reproductive individuals in the area (around 30 years old) (Montesinos Torres 2007; Acuñ a-Míguez and others 2020).

Sample Collection and Trait Assessment
We chose plant traits from two dimensions important to define the global plant spectrum (Díaz and others 2016). The trait dimensions included plant size, which reflects the ability to compete for light and occupy space (in both the vertical and hori-Early Differentiation of the Phenotypic Space and Performance zontal planes); and leaf traits, morphological and chemical characteristics linked to trade-offs between acquisition and conservation of resources and water use (Wright and others 2004;Chave and others 2009;Díaz and others 2016). Specifically, in the field, we measured individual height (H, m); diameter at ground level (D, cm); crown size, calculated with two perpendicular diameters of the crown projection (CP, m 2 ) and height at which the crown is larger (Hmax, m), whose comparison to total H (Hmax:H) indicates if an individual has a more conical (Hmax in the lower third of H), oval (mid third) or vase (upper third) shape. For multistem individuals, we measured all diameters and assessed the quadratic diameter (the square root of the sum of squared diameter of each stem) (Stewart and Salazar 1992).
In addition, we collected one branch from a fully exposed part of the crown to assess leaf traits, which were stored in cool and fresh conditions until they were taken to the laboratory (within the next 48 h). There, branches were rehydrated for at least 12 h in the darkness (Garnier and others 2001). Juniperus thurifera has scaly leaves adpressed to the branches, hindering the standard estimation of leaf traits (Perez-Harguindeguy and others 2013). Therefore, we selected five of the youngest green twigs as units analogue of leaves (Paula and Pausas 2006) and weighed them for the saturated weight. They were then oven-dried at 60°C for 72 h and weighed for dry mass, and we estimated leaf dry matter content (LDMC, the oven-dry mass of an individual leaf-like unit divided by its watersaturated fresh mass, mg 1 g -1 ). This trait is sensitive to variations in nutrient and water supply, similarly to others more frequently used based on leaf area (for example, specific leaf area or leaf mass per area) but more reliable for non-flat leaves (Wilson and others 1999).
Then, we ground oven-dried samples, and took between 3.5 and 4.0 mg, weighed them and sealed them in tin foil capsules for chemical traits, namely C and N content (LCC and LNC, %) and the isotopic signatures of stable isotopes 13 C (d 13 C, &) and N (d 15 N; &). These samples were shipped to UC Davis Stable Isotope Facility, where the chemical analyses were carried out with an elemental analyser interfaced to a continuous flow isotope ratio mass spectrometer (IRMS) (for details see https:// stableisotopefacility.ucdavis.edu/). Leaf N correlates with the photosynthetic capacity (Field and Mooney 1986) in a way that nutrient-rich environments favour N allocation to the photosynthetic machinery, instead of to supporting tissue which is more typical in poor environments (Berendsen and others 2007). Leaf C content is the structural basis and its relationship with other nutrients (leaf stoichiometry) that can be indicative of responses to stress (Zhang and others 2017). For instance, an increase in the C:N ratio in plant tissues can reflect a shift from photosynthetic to structural tissue allocation under harsher conditions (Å gren 2008). The abundance of 13 C in leaves is linked to plant's water use efficiency (WUE). Plant photosynthesis discriminates against the stable 13 C isotope (Farquhar and others 1989) until stomata close and intercellular CO 2 concentration drops. Then, higher levels of d 13 C are achieved under prompt stomatal closure. The isotopic signature of stable isotope 15 N in leaves reveals spatial and temporal patterns of N cycling (Amundson and others 2003) and varies among species, ecosystems and climatic gradients (Craine and others 2009 and references therein).
Finally, in each individual, we also collected one wood microcore or core at ground level (depending on the individual stem size, totalling 155 cores and 25 microcores) to assess performance using a Trephor tool (2 mm diameter (Rossi and others 2006)) or increment borers coupled to a drill (5 mm diameter; Pressler Haglö f), respectively. The sampling was entirely conducted in April 2019.

Performance Measurements
We used the core extracted per tree to count and measure annual ring widths. We glued the cores into wooden mounts, air-dried, and polished them using sandpapers of progressively finer grain until tree rings were visible. We scanned cores (2400 d.p.i.), dated them with Cdendro v9.1. (Cybis Elektronik & Data AB, 2018) and estimated individual age by counting rings until the pith or using Cdendro software (Table 1). We measured ringwidths to an accuracy of 0.01 mm using the software CooRecorder v9.1 (Cybis Elektronik & Data AB, 2018) and checked cross-dating of individual series using CooRecorder and COFECHA programs (Holmes 1983) (Table 1). We then calculated the annual basal area increment for each tree (BAI, cm 2 AEy -1 ), as it provides a better proxy for biomass growth even for juvenile growth (Rozendaal and others 2010; Brü llhardt and others 2020). We assessed the average of the annual BAIs as a measure of secondary growth that in our case is uncorrelated with age (see Supplementary Information, Figure S2).
To study juvenile drought resilience, we first selected an outstanding drought event within the life span of most of our study individuals (on average 18 years old). We used the 6-month SPEI (standardized precipitation-evaporation index) from 2006 to 2016, that is, a 10-year period prior the sampling year. It represents a multiscale drought index that accounts for both temperature and evapotranspiration effects on the water balance. Negative SPEI values indicate dry periods, whereas positive SPEI values indicate wet periods (Vicente-Serrano and others 2010). These values were downloaded from the SPEIbase (http://sac.c sic.es/spei/database.html, SPEIbase v2.6, Beguería and others 2018), and we chose 2012 as drought year because the two annual values, both in spring and autumn (Camarero and others 2010) showed the lowest SPEI values of all the study period (SPEI < -2; Figure S3, Supplementary Information).
Then, we calculated the components of resilience, that is, resistance, recovery and resilience (Rt, Rc and Rs, respectively), described by Lloret and others (2011) using the basal area increments in 2012, BAI 2012 , our drought year reference, and averaged BAI in preceding and posterior years (BAI pre , BAI post ). In short, resistance (Rt) is the ability to endure a stress situation induced by a drought event without reducing growth (equation 1); recovery (Rc) reflects the growth increase after the drought event (equation 2), and resilience (Rs) indicates the ability of a tree to reach a predrought growth level (equation 3).
We chose a five-year period to represent preceding (2007-2011) and posterior (2013-2017) periods to the drought event. We also obtained the indices with a 3-year period (considering 2009-2011 and 2013-2015), in order to increase the sample size (the youngest trees are not discarded), and had similar results (Gonzá lez de André s and others 2021).

Statistical Analyses: the Phenotypic Space
We analysed the variation of juveniles' phenotypes across the woodland expansion gradient (or stages), namely EA, TR, MA. First, using linear mixed models, we explored the distribution of the variance of each trait across spatial scales, that is, among sites, stages within the expansion gradient, plots within stage and among individuals within plots represented by the residual variance (equation 4). Before, trait values were log transformed, as appropriate, to optimise normality of the residuals.
Then, we explored the variation of phenotypes with an integrated (multivariate) perspective and estimated the functional trait hypervolume (functional trait space) of juveniles in each stage across the woodland expansion gradient following the methodology described by others (2014, 2018). Before running the hypervolume estimations, we reduced the number of traits (axis) to five (see Supplementary Information for details about trait selection), so the ratio between observations and number of dimensions was over 10 (Blonder and others 2014). We normalised the data (z-scores) of selected traits (LDMC, LNC, d 13 C, CP, Hmax:H) to have axes with comparable units (standard deviation). Then, we obtained the size of hypervolumes, which represents the variation of the phenotypic space (expressed by fifth power of standard deviations), and their centroids at each stage across the forest expansion gradient, which represent the (standardised) mean trait values (Blonder and others 2014). We tested whether the size of functional space (hypervolumes) enlarges or shrinks (that is, whether traits diverge or converge in each stage), and whether the functional space shifts across the woodland expansion gradient using null models (Gotelli and Graves 1996). We built up hypervolumes based on 999 null communities (null expectations) by randomly selecting 999 times 60 individuals from the total pool, no matter their actual stage within the woodland expansion gradient, and we obtained from each realisation the hypervolume size and the centroid coordinates that allowed us to calculate pairwise distances among centroids of realisations. In biological terms, these null communities represent communities where the functional space occupied by T. thurifera juveniles does not vary across a woodland expansion gradient (no land-use legacy). Then, we compared distances among centroids of observed hypervolumes along the woodland expansion gradient to the distribution generated with those assessed from null communities. We obtained a P value (equation 5) that represents the probability of getting the observed distance or lower, and a statistically significant shift was therefore considered when (1 -P) < 0.05.
Similarly, we tested whether the size of the observed hypervolumes were significantly different to those obtained for the null communities; thus, we calculated again the probability P (equation 5) that the observed value is equal or lower than expected in a null distribution. We standardised the metric to range from -1 to 1 obtaining the effect size ES (equation 6) (Bernard-Verdier and others 2012). We considered that the observed hypervolume is statistically significant greater than by chance (trait divergence) when ES > 0.95 (at the 0.05 significance level), and smaller than by chance (trait convergence) when ES < -0.95. We performed a SEM for each performance metric (BAI, Rt, Rc, Rs), where performance was modelled by extrinsic and intrinsic factors with a random intercept for plots within sites (plot nested in site) to consider the autocorrelation within locations. We also run a joint SEM for all perfor-mance metrics to consider any interaction among them. We allowed extrinsic factors to directly affect the performance, as well as indirectly through shaping functional traits. Parameters reflecting the extrinsic factors were the biotic and abiotic environmental variables recorded at the individual level together with the stage (woodland expansion gradient: EA, TR, MA) that may capture environmental differences due to the underlying land-use legacy. Parameters representing the intrinsic factors were the age and functional traits. Note that plantsize traits (H, D, CP) were highly correlated (q > 0.7) (Supplementary Information, Figure S4); therefore, we kept CP because in this low-density forest horizontally space occupation is more meaningful of individuals' ability to compete for resources than vertically. Similarly, C:N and d 15 N were strongly correlated with LNC; therefore, they were removed from further analyses, using the five traits used to define the phenotypic space. We expected that shrub and tree covers would affect architectural traits (CP and Hmax:H) and water use-related traits (d 13 C) given the nursery effect, frequent in Mediterranean areas (Gó mez-Aparicio and others 2004). Whilst the competition with herbaceous and shrub vegetation (Benavides and others 2016) and micro-topography and slope would affect traits related to water and resources Coverherb, covershrub and covertree represent the cover occupied by herbaceous vegetation, shrubs and the projection of adult trees' crowns around each juvenile (area delimited by the projection of the crown + 2 m). LDMC: leaf dry matter content; LNC: leaf nitrogen content; d 13 C: isotopic signature in leaves of 13 C; CP: crown projection; Hmax:H: ratio between height at which the crown is horizontally larger and total height; BAI: average of the annual basal area increments; Rt, Rc and Rs represent the resistance, recovery and resilience of juveniles against a drought event.
use (d 13 C, LDMC and LNC). We had previously checked whether the environmental factors varied across stages (Supporting Information, Figure S5) and verified that only tree cover increased significantly from the expanding area towards mature woodlands (correlation between cover tree and stage). This result agrees with a previous work in the same plots that showed a significant difference in tree density between EA and MA (c.a. 28 trees/ ha and 116 trees/ha on average, respectively (Alfaro and other 2021).
Given that stage (woodland expansion gradient) was categorical, and no standardised estimates can be assessed, we recoded the values of the factor, namely EA = 1, TR = 2, MA = 3, and then model as numeric. Thus, the coefficient represented the expected change in the response variable as going from the expanding area to mature woodlands (Lefcheck 2021). When our full model had an inadequate goodness-of-fit (P value < 0.05 of Fisher's C), we reconsidered our model and checked that all relevant pathways were included. Then, we looked for variables with no relationship to other response variables and removed them, one by one, and compared the model to their respective full model (nested model) through their Fisher's C statistic according to Shipley (2013) and Lefcheck (2021), similarly to the use of AIC when comparing nested models (Zuur and others 2009).
All statistical analyses were performed in R (version 4.0.2, R Core Team 2020), using the stats package for correlation analysis and linear models, lme4 package (Bates and others 2015) and MuMIn (Barton 2020) for linear mixed models. We also used hypervolume (Blonder 2019) and piecewiseSEM (Lefcheck 2016), emmeans (Lenth 2021) and nlme packages (Pinheiro and others 2020) to, respectively, conduct the functional space estimation and the piecewise SEM analysis.

RESULTS
Sampled juveniles were on average 21 years old (SD = 9.35), with younger individuals in the expansion area compared to those in mature woodlands, 15.7 y (SD = 5.94) and 26.3 y (SD = 9.61), respectively (P value = 0.03) (Supporting Information, Figure S5). Their phenotypes across the forest expansion gradient only differed in some leaf traits (LNC, C:N, d 13 C) and crown shape (Hmax:H), capturing around 10% of the total variation in each one (Figure 3). However, important differences appeared among sites in all traits (going from 4 of LDMC to 54% of d 15 N of the total variance), within stages (plot accounted for 20-30% of the variance of LDMC, LCC and d 15 N,) and within plots with over 50% of the variance found among coexisting individuals (except for LCC and d 15 N which accounted for 43 and 22%, respectively) (Figure 3).
Results from the integrated phenotypic space (hypervolumes) showed that juveniles in mature woodlands were different to those in opener places, with greater distance between hypervolumes' centroids than expected by chance (Table 2; Supplementary Information Figure S6). Juveniles in woodlands had leaves with lower dry matter content, d 13 C and LNC, and had smaller and squatter crowns (conical shape). Moreover, juveniles growing in both extremes of the gradient, in woodlands and expanding areas, had trait values that tended to converge (ES = -0.83, ES = -0.89, respectively) although they did not reach the significance threshold -0.95, which would indicate significant convergence of traits.
The overall model showed no relationship between radial growth (BAI annual ) and any of the resilience metrics, therefore, results are shown in independent networks. The piecewise SEM for radial growth (BAI annual ) suggested that stage (landuse legacy) has both direct and indirect effects on juveniles' growth (Fisher's C = 11.84 with P value = 0.921, Figure 4). Juveniles with greater increments in secondary growth appeared in the woodland-expanding area (EA), where individuals had larger amount of LNC and d 13 C in leaves. Although the isotopic signature of 13 C in leaves did not directly affect growth, it was positively correlated with crown size (CP) and the amount of N in leaves, which in turn were positively associated to basal area increments (Figure 4). Age did not affect directly BAI, but through CP, that is, older individuals had bigger crowns and grew more. Besides, conical architectures (lower values of Hmax:H) also performed slightly better, with bigger mean annual BAI.
The final piecewise SEM for resilience components showed that juveniles in woodland-expanding areas (EA) resist more (that is, growth drop during a drought event was less intense in the expanding area) and tend to recover quicker, resulting in a higher resilience in the EA compared to those juveniles in mature woodlands ( Figure 5). Indirectly, stage was connected to resistance through d 13 C in leaves (proxy of water use efficiency), appearing higher resistance in the expanding area where juveniles had higher d 13 C in leaves. Resistance to drought was also positively associated to smaller crowns and younger individ-uals (Fisher's C = 5.440, P value = 0.86; Figure 5a).
Stage also affected indirectly the recovery and resilience (Fisher's C = 13.192 with P value = 0.659 Figure 5b, and Fisher's C = 11.148 with P value = 0.516, Figure 5c, respectively) through effects on age and LCN, being more resilient young individuals. Similarly to BAI, recovery and resilience were also affected by the amount of N in leaves (positively) and the size and shape of the crown, this time individuals with bigger and squatter crowns (that is, larger part of the crown at the base: lower Hmax:H) yielding lower performance after a drought event.
Finally, the comprehensive piecewise SEM including the three resilience metrics (shown in Supplementary Information Figure S8) showed that resistance against drought and recovery after a drought event were positively correlated to resilience, in agreement to their definition (partial correlations to resilience of 0.457 and 0.519, P value < 0.001, respectively, Fisher's C = 21.508 with P value = 0.920). At the same time, they were themselves negatively correlated showing a tradeoff between them (partial correlation -0.455, Pvalue < 0.001).  Hypervolumes represents the phenotypic space expressed by fifth power of standard deviations (SD 5 ), based on five functional traits: leaf dry matter content, leaf N content, isotopic signature of 13 C in leaves, crown projection area and ratio between height at which crown is larger and total height. P value represents the probability of an observed value in a null distribution where juveniles' phenotypes are not led by the environmental gradient across the woodland expansion. SES (standardised effect size) shows whether the observed hypervolume is significantly larger (SES > 0.95, divergence of traits), lower (SES < -0.95, convergence of traits) or not different to what is expected by chance. For distances between centroids, 1 -P represent the significance threshold. Bold values show significant different to those expected by chance.

DISCUSSION
We found relevant land-use legacies on the phenotypic space of J. thurifera juveniles that drove early differentiation of key functional traits. This early phenotypic differentiation rendered into juveniles colonising old-agricultural fields with higher performance in terms of growth and resilience against droughts, compared to those growing in mature woodlands. Our results suggested that the functional differentiation previously observed in J. thurifera mature woodlands (Acuñ a-Míguez and others 2020; Villellas and others 2020) is therefore the result of early trait acclimation of traits, especially on traits related to water and resource use like leaf isotopic signature of 13 C and N content that remains over time.

Environmental Conditions and Trait Differentiation
Establishment in open areas under Mediterranean climates implies great abiotic stress in summer, triggered by high irradiance and rainfall scarcity.  case, we did not find any trace of facilitation as shown in a previous experiment with transplanted saplings in the area (Gimeno and others 2012a), as a higher amount of shrubs around juveniles in the expanding area would have suggested. Nevertheless, our survey does not imply a lack of nursing effect on saplings either, but rather a lack of clear relationship between shrubs and juveniles that recruited relatively long ago.
In addition, we did not find differences in the study local factors across the expansion gradient, except for a higher tree cover in mature woodlands. Neither did our data capture significant direct effects of extrinsic factors (local factors) on juveniles' phenotypes, likely because they were either coarsely recorded or non-representative of the main environmental variability. Note that, despite canopy cover varied across the woodland expansion Figure 5. Final piecewise SEM representing the factors that affect juveniles' resistance (Rt), recovery (Rc) and resilience (Rs) against a drought event (see Appendix Figure S7, for the overall model). Only significant relationships are shown, with their path coefficient and significance (***p < 0.001, **p < 0.05, AEp < 0.1). Solid lines represent positive associations and dashed lines negative associations. Traits in white background and thinner lines are traits without direct effect on the response variable. R2 of the response variables of the different pieces of the models are included: R2 marginal (R2m) considering fixed effects and R2 conditional (R2c), fixed and random effects. Grey lines and coefficients represent non-significant trends (AEp < 0.1). LNC: leaf nitrogen contents, d 13 C: isotopic signature in leaves of 13 C; CP: crown projection: Hmax:H: ratio between height at which the crown is horizontally larger and total height; covertree: cover occupied by the projection of adult trees' crowns; coverherbs: cover occupied by the herbaceous vegetation around each juvenile (area delimited by the projection of the crown + 2 m). gradient, it did not exert direct effect on any trait nor performance metric, likely due to the low tree density within the plots (ranging from c.a. 28-116 trees/ha, see Alfaro and other 2021), not imposing a relevant differential environmental pressure and suggesting that the expansion gradient encapsulates other extrinsic factors.
However, we did find a shift in the phenotypic space in juveniles, with squatter (maximum crown size closer to the ground) and their leaves with less N content and d 13 C in mature woodlands than in the expanding edge. Previous studies in the area showed a similar pattern in leaves of adult trees, with higher d 13 C and lower C:N (inversely correlated to our LNC) in the expanding edge (Gimeno and others 2012a; Acuñ a-Míguez and others 2020). These authors suggested that higher stoniness (implying less soil water-retention capacity) and radiation (higher evapotranspiration) in the recently established area would have triggered this leaf d 13 C differentiation (proxy of the WUE), namely to adjust individuals' water use under harsher conditions (greater d 13 C-WUE-with less water availability). In our sites, high irradiance in the expanding edge can lead this increment in the leaf proportion of 13 C (Aranda and others 2007). However, we did not find higher stoniness there, at least closely around our surveyed juveniles (recorded as % cover of stones/rocks within each juvenile influence area, data not shown), neither did another study in the area others 2012a, 2012b).
Leaf N content is related to the photosynthetic capacity (Field and Mooney 1986) and, at the interspecific level, normally decreases under increasing irradiance-light-demanding species versus shade-tolerant species-and in nutrientpoor soils (Givnish 1988;Meziane and Shipley 2001;others 2001, 2005;Hanba and others 2002), while it usually increases in species from drier sites others 2001, 2004) likely due to a larger volume fraction of mesophyll cells within the leaf (Roderick and others 1999). However, at the intra-specific level, that is, comparing individuals within species across gradients, contrasting results turned out (Siefert and others 2015;Wright and Sutton-Grier 2012;Niinemets 2015), reflecting the complex interactions among environmental factors and individuals, determined by local conditions (limiting factors) and physiological trade-offs. Across the study gradient, we found higher LNC in the expanding edge in agreement with a typical response to water scarcity caused by a higher evapotranspiration or/and a better nutrient status of soils there, due to the le-gacy of past agricultural practices (Flinn and Marks 2007; Alfaro-Sá nchez and others 2019; Acuñ a-Míguez and others 2020).

Extrinsic and Intrinsic Factors Shaping Performance
Aforementioned trait differentiation across woodland expansion gradients triggered differences in juveniles' performance. Hence, we found that juveniles in the expanding edge grew more, as we had expected (see Figure S9, SI), whereas they were more resilient after a drought than those in mature woodlands, refuting our hypothesis. These results agreed with previous studies in the area, which showed an unexpected higher survival of saplings in recently established areas than in 'safer sites' within mature woodlands (Gimeno and others 2012a) and higher growth rates in adult trees attributed to a lower competition in opener areas (Alfaro-Sá nchez and others 2021).
Together with the extrinsic factors encapsulated within woodland expansion gradient, our results showed a surprising positive effect of the herbaceous layer on juveniles' recovery. Although a nursing effect exerted by shrubs has been often reported (Gó mez-Aparicio and others 2004; Castro and others 2006; Valladares and others 2014; Benavides and others 2016), herbs have been usually seen as direct competitors against seedlings and saplings for water (Davis and others 1999;Benayas and others 2005;Caldeira and others 2014). Therefore, the positive correlation with drought recovery might indicate a competition release with juveniles that recruited long ago and whose root systems are now deeper than herbs' (Redmond and others 2018), and a potential role slowing down soil evapotranspiration while increasing nutrients' concentration (Gilliam 2007).
Regarding the intrinsic factors, that is, functional traits shaping juveniles' performance, LNC and d 13 C were positively related to juvenile performance. In particular, LNC favoured radial growth and recovery after a drought, similarly as (Gonzá lez de André s and others 2021) found for other tree species in North Spain and as could be expected due to its relationship to the photosynthetic capacity (Field and Mooney 1986). Meanwhile, d 13 C was positively associated to a higher resistance to drought, and indirectly to radial growth, as it is an integrative estimation of the ratio between photosynthesis rate and stomatal aperture (WUE) over a determined growth period (Farquhar and Richards 1984). Interestingly, other intrinsic factors drove performance, namely age and crown size (both positively correlated among themselves and with plant height and diameter at ground level) with contrasting directions depending on the performance metric. They correlated positively with radial growth and inversely with resilience expressed in any of their components. It has been long recognised that during juvenile phase, tree growth increases exponentially with age, and for a given amount of available resources, juveniles with larger crowns can harvest more light, increasing their carbon gains and growing more than smaller ones until a size threshold is reached. For our study species, this peak occurs c.a. 30-50 cm diameter breast height (Gó mez-Aparicio and others 2011), which is far from our study size range. However, when water is scarce, bigger crowns consume more water (higher respiration costs) and taller trees have longer hydraulic paths and experience more gravitational potential opposing water ascent (Mencuccini and others 2005;Ryan and others 2006), which can reduce resilience to drought. In the same line, crown shape of juveniles also affected individual performance, with squatter or conical crowns (lower Hmax:H) favouring the individual recovery and resilience compared to those more oval (taller crowns), supporting the hypothesis of hydraulic limitation.
A limitation of our study is that we did not consider any belowground trait. In water-limited regions, like the Mediterranean area, survival of juveniles after drought events is related to the size and length of roots able to search for water in deep layers, that is, to a higher investment in root growth (Padilla and Pugnaire 2007;Ramírez-Valiente and others 2018). Thus, belowground traits indicating a higher investment in the root system that would allow the access to deeper soil moist layers could also underlie our pattern, favouring a greater annual growth in expanding areas, despite a lower water availability due to higher irradiance and evaporation.

Relationship Between Performance Approaches Within a Climate Change Scenario
The Mediterranean climate demands certain growth resilience to drought, which is directly associated to individual vulnerability and longterm survival (DeSoto and others 2020). Although previous studies have shown that fast-growing trees are commonly less resistant to drought than slowing-growing trees, the former usually have a faster recovery (Martínez-Vilalta and others 2012; Serra-Maluquer and others 2018; Alfaro-Sá nchez and others 2019; Fang and Zhang 2019). However, in this study, we did not find any relationship between total growth and any component of drought resilience, showing no trade-off between growth rate and resilience at this early life stage (based on a comprehensive SEM not shown, which showed no connection). A possible explanation is that at an early-life stage, when reproduction is not an issue and individuals are rather vulnerable to physical and drought physiological damages, fast secondary growth may involve a quick escape from this vulnerable stage, bestowing an advantage for surviving (Zamora and others 2001;Benavides and others 2015) that might cancel out the negative relationship found on adult trees.
In agreement with Lloret and others (2011) definition of resilience used in this study, the capacity to reach pre-drought growths (resilience) was positively influenced by both individual's resistance (reversal of the reduction in ecological performance during disturbance) and recovery (ability to recover relative to the damage), two components that are complementary. However, different strategies of forests to cope with drought imply different weights of both components, which depend on the species, biomes and climatic conditions (Gazol and others 2017). Here, we found a trade-off between resistance and recovery (Lloret and others 2011). This trade-off previously described in literature others 2017, 2018; Hoffmann and others 2018) entails that those individuals that resist better the water scarcity, maintaining growth levels during the drought, have usually slower recovery in relation to the predrought mean of radial growth. We also found higher correlation and more similar SEM structures between resilience and recovery ( Figure 5), compared to those between resilience and resistance, which suggest a higher impact of recovery on final resilience of juniper juveniles. This outcome is in the line of a comprehensive study conducted by Gazol and others (2018), who analysed resilience of 11 tree species (including our study species) across biomes within the Iberian Peninsula. They showed that resilience in evergreen gymnosperms dominating semi-arid forests (our study case) mainly lay in a high recovery rather than resistance, which is more typical in species dominating humid forests.
To conclude, alarming data are evidencing a progressive and substantial drying of the Mediterranean land surface since last century, with Early Differentiation of the Phenotypic Space and Performance drought and heat wave durations and frequencies increasing more than six-fold since the 1960s (EEA 2017). Consequently, studies tackling forest resilience through time, despite scarce, have shown a drop in the last decades of any of the resilience components, as a response to the increasing extreme events occurrence (Serra-Maluquer and others 2018; Fang and Zhang 2019). On the other hand, recent studies have posed that forests in drought-prone areas may develop an ecological drought memory (Walter and others 2011) that might partially offset the negative effects of recurrent droughts, through a process named priming (Hilker and Schmü ulling 2019), which entails acclimation of phenotypes that allows them to improve their response to future stress. In our study area, outcomes about juveniles (ours and in Gimeno and others 2012a), together with other studies based on adult J. thurifera trees (Acuñ a-Míguez and others 2020; Villellas and others 2020; Alfaro-Sá nchez and others 2021), have showed phenotypic adjustments that are favouring higher radial growth and drought resilience. This suggests that these new J. thurifera forests are not particularly vulnerable under current conditions and may be outcompeting coexisting species which are more affected by climate changes (Gimeno and others 2012a). Further research is needed to unveil whether these adjustments are preparing these new woodlands for future stress, in which including other coexisting species and soil variables would help us to understand better these new forests and forecast their future.

AC KNOWLEDGEMENT S
We are in debt to all the people collaborating with the field and laboratory work collecting and assessing data, especially to David Ló pez-Quiroga, Antonio Mas, Borja Ruiz-Valdepeñ as Martin, Ana Sanz Barciela, Cristina Nevado Luna and André s Bravo.

FUNDING
Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Funding was provided by Spanish projects REMEDINAL TE-CM (Comunidad Autó noma de Madrid, S2018/ EMT-4338) and COMEDIAS (CGL2017-83170-R, Ministerio de Ciencia, Innovació n y Universidades). RB was funded during the conductance of the study partially by the Contrato Programa CEEC Individual, Fundaçã o para a Ciê ncia e Tecnologia, Portugal.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/li censes/by/4.0/.