Mixed Pine Forests in a Hotter and Drier World: The Great Resilience to Drought of Aleppo Pine Benefits It Over Other Coexisting Pine Species

Drought is an important driver of forest dynamics in the Mediterranean region. The forecasted increase in drought frequency and severity can notably influence tree growth, forest structure, composition and productivity. Understanding how coexisting tree species respond to drought is thus crucial to understand which are less vulnerable and will perform better in a warmer and drier world. To assess drought vulnerability, we used dendrochronology to study the radial growth trends and responses to a drought index of four pine species (Pinus halepensis, Pinus pinea, Pinus nigra, and Pinus sylvestris) coexisting in North-eastern Spain. We reconstructed the growth of each species and evaluated their short- and long-term growth response to drought for the common period 1980–2017. The growth of the four pine species depended on water availability and high early spring temperatures impacted the growth of P. nigra and P. sylvestris negatively. The occurrence of a severe drought between 2005 and 2007 lead to marked growth reductions in the four species, but it was greater in magnitude in P. pinea and P. halepensis in 2005, and in P. nigra in 2007. The results of basal area increment models at the individual tree level suggested that P. halepensis trees grow more than the rest of species. After accounting for age and drought effects, P. nigra and P. sylvestris displayed negative growth trends in the 2008–2017 period while P. pinea and P. halepensis displayed positive growth trends. P. sylvestris was the most resistant species and P. pinea the less resistant. Conversely, P. halepensis and P. pinea were slightly more resilient than P. sylvestris. Moreover, P. sylvestris was the species displaying the highest autocorrelation and the lowest coefficient of variation in ring-width indices. A marked drop in the autocorrelation of P. pinea ring-width index was observed in response to the 2005 drought. These results indicate that all study species are vulnerable to drought but in different degrees. The strong resilience capacity of P. halepensis suggests that it will better thrive in a drier future, but mixed pine forests, such as the one here studied, may contract or become rare due to the strong sensitivity of P. pinea to drought and the lower post-drought performance of P. nigra and P. sylvestris.


INTRODUCTION
As climate warms, a rearrangement of plant communities looms (García-Valdés et al., 2015;Buras and Menzel, 2019) impacting the structure, composition and functioning of forests (García-Valdés et al., 2020). Currently, climate warming has boosted the occurrence of severe droughts (IPCC, 2022). It has been recognized that drought is a major constraint of tree radial growth worldwide (Babst et al., 2019) and a recurrent cause of forest dieback and mortality both directly and indirectly (Allen et al., 2015;Hartmann et al., 2022). The Mediterranean Basin represents a climate change hotspot in which forests are suffering the direct and indirect consequences of increasingly warmer and drier conditions (revised in Peñuelas and Sardans, 2021). Directly, a warmer climate and more intense droughts may cause tree death by inducing xylem hydraulic failure among other responses (Hereş et al., 2014;Camarero et al., 2015a;Pellizzari et al., 2016;Hevia et al., 2019). Indirectly, severe droughts can also enhance tree mortality by weakening trees, reducing their carbon pools and making them more vulnerable to pathogens (McDowell et al., 2008;Gaylord et al., 2015;Stephenson et al., 2019;Trugman et al., 2021). A key question emerging here is: when coexisting, which tree species will be more vulnerable to drought?
In drought prone regions such as most Spain, differences in radial growth rate, wood traits and xylogenesis are important factors determining how trees respond to water shortage separating angiosperms from gymnosperm species (Carnicer et al., 2013;Camarero, 2018). Under drier scenarios, forests will shift toward less productive and more drought tolerant communities (García-Valdés et al., 2021). In the long term, drought tolerant angiosperms such as anisohydric oaks will be benefited against drought vulnerable gymnosperms such as isohydric pines in the case of mixed pine-oak stands (Quero et al., 2011;Gea-Izquierdo et al., 2014. However, differences exist within pines which may become more important when coexisting (Camarero et al., 2015a;Marqués et al., 2016;Salazar-Tortosa et al., 2018;Andivia et al., 2019). Pines represent an important part of Spanish forests accounting for ca. 40% of the total forested area which is ca. 15 million ha (Astigarraga et al., 2020). Therefore, assessing how they respond to drought is fundamental to advance in our understanding of how future Spanish forests will look like in the near future (Andreu et al., 2007;Camarero et al., 2021a;Herrero et al., 2021). Poleward and upward migrations of Mediterranean tree species are forecasted in the long-term (Buras and Menzel, 2019), but these projections should also consider the role played by drought as constraint of performance in some pine species (Andreu et al., 2007).
It is expected that Mediterranean pines with lower lethal water potential, i.e., more resistant to xylem embolism and the loss of hydraulic conductivity, and a facultative bimodal growth pattern, i.e., potentially growing in spring and autumn, could better tolerate drought than drought-vulnerable temperate pines with unimodal growth patterns showing a major growth peak in spring (Camarero et al., 2010). That is, Mediterranean species with either a facultative or obligated bimodal growth pattern such as Pinus halepensis Mill., Pinus pinaster Ait. or Pinus pinea L. (e.g., Castagneri et al., 2018;Campelo et al., 2021) may have advantages over other species from mesic sites with unimodal growth patterns such as Pinus nigra J.F. Arn. or Pinus sylvestris L. (Rossi et al., 2006). Experimental evidence indicates that the weaker stomatal control of P. halepensis as compared to P. sylvestris and P. nigra may allow this species to maximize carbon uptake during drought (Borghetti et al., 1998;Salazar-Tortosa et al., 2018). These species also differ in the amount of biomass invested in belowground organs, with species such as P. halepensis and P. pinea creating long-lasting large root systems which may potentially reach deep water sources (Sarris and Mazza, 2021). Besides, common garden experiments suggest that Mediterranean pine species have also larger root systems than pine species from more mesic sites (Andivia et al., 2019).
It is unclear how pine species will respond to drought when they coexist under similar climate conditions. However, finding sites where several pine species coexist is challenging (e.g., Granda et al., 2018;Hernández-Alonso et al., 2021), and thus those regions offer valuable and informative settings to test which species are move vulnerable to drought. Understanding how trees responded to past droughts through the study of tree rings may provide clues in this respect (Andreu et al., 2007;Eilmann and Rigling, 2012;Camarero et al., 2015a). Past droughts leave imprint information in tree-ring records (Fritts, 1976), which allows understanding how surviving individuals responded to them (Ogle et al., 2015). Droughts cause abrupt growth reductions in tree growth (Schweingruber, 1986) that can last for several years creating carryover effects or drought legacies (see Anderegg et al., 2015) and impacting drought resilience capacity (Lloret et al., 2011). Importantly, tree growth responses to past droughts may determine the likelihood of dieback and mortality occurrence in the future and be used as early warning signals of tree vulnerability (Camarero et al., 2015a(Camarero et al., , 2021aCailleret et al., 2017;De Soto et al., 2020;Keen et al., 2022).
Here we studied growth patterns and responses to drought of four co-occurring pine species (P. sylvestris, P. nigra, P. pinea, and P. halepensis) coexisting in mixed pine forests located in Northeaster Spain. We take advantage of a natural experiment that evaluated the impact of a Diplodia shoot blight (Diplodia sapinea (Fr.) Fuckel (Syn: Sphaeropsis sapinea (Fr.) Dyko and B. Sutton).) on the four co-occurring pine species (Caballol et al., 2022). We hypothesize that recent warming trends and the occurrence of severe droughts may have left an imprint in the radial growth series of those species offering valuable information to infer their vulnerability to water shortage. Thus, we investigate differences in growth and its response to drought between species over the common period 1980-2017 (prior to the Diplodia outbreak). We hypothesize that: (i) the growth of the four pine species will respond to drought and all species will present drought-induced growth declines, but (ii) P. nigra and P. sylvestris, species found outside the Mediterranean Basin in other Eurasian regions, will be more affected by climate warming and recent droughts than species restricted to the Mediterranean Basin (P. pinea and P. halepensis). This is supported on the recently observed declines of P. sylvestris in its equatorward distribution limit or rear edge (e. g. Gea-Izquierdo et al., 2014;Camarero et al., 2015a) and its poor adaptation to withstand drought (Camarero et al., 2010). Along this, recent droughtinduced dieback and mortality episodes of P. nigra and P. pinea in Mediterranean dry regions (Linares and Tíscar, 2010;Petrucco et al., 2017) suggest that they will be more vulnerable to drought than P. halepensis which will be the species more able to couple with drought among those coexisting due to its above-and belowground growth plasticity (Borghetti et al., 1998;Pacheco et al., 2016;Salazar-Tortosa et al., 2018;Andivia et al., 2019;Campelo et al., 2021;Sarris and Mazza, 2021).

Study Site and Climate Data
The study was conducted in mixed pine forests located near the village of Oristà, central Catalonia, North-eastern Spain (41 • 55 58.4 N, 2 • 3 35.7 E, 468 m a.s.l.). The studied stands were composed by mixed pine forests of P. sylvestris, P. nigra, P. pinea, and P. halepensis. In the year 2020, eight stands in the area were sampled for dendroecological purposes (Caballol et al., 2022). Two stem cores were taken at 1.3 m using Pressler increment borers (Haglöf, Sweden) for trees from the four pine species to reconstruct their growth patterns and estimate their age. Here we analyze a total of 123 trees sampled in six of these stands ( Table 1).
The climate in the region is humid subtropical or warm temperate (Cfa) according to Köppen-Geiger classification (AEMET-IMPA, 2011). The mean annual temperature in the region is 12.6 • C, the coldest and warmest months are January (4.0 • C) and July (22.0 • C), respectively (Supplementary Figure 1). The total annual precipitation is 655 mm (Supplementary Figure 1). The soils are thin and acid, developed on shales, sandstones and conglomerates. The vegetation is composed by mixed conifer forests with accompanying angiosperms such as Quercus ilex L. The presence of a Mediterranean climate with spring and autumn precipitation peaks facilitates the co-occurrence of the four pine species (Figure 1).
We obtained monthly mean temperature ( • C) and precipitation (mm/day) data for the common 1980-2017 period from the 0.25 • -gridded E-OBS climate dataset (Cornes et al., 2018). We used the Standardized Precipitation  Evapotranspiration Index (SPEI; Vicente-Serrano et al., 2010) to characterize drought severity at temporal resolutions from 1, 3, 6, 9, 12, and 24 months (Supplementary Figures 2, 3). SPEI values were obtained at ∼1.1 km 2 spatial resolution for the study region (Oristà) from the Spanish SPEI database (Vicente-Serrano et al., 2017). 1 We selected the 12-month long SPEI in further analyses since it captured the main drought episodes in the study area.
To describe the climate niche of each species we first downloaded distribution maps across Europe from the European Forest Genetic Resource Programme website. 2 These maps were created by combining different sources of information across Europe (Caudullo et al., 2017). Mean annual temperature and total annual precipitation were downloaded from the WorldClim database (Fick and Hijmans, 2017). Polygons from the distribution for each species were converted into points and the climate conditions were obtained from each point. The climate space for each species was described by the climate conditions in each point (see Figure 1).

Tree-Ring Analyses
A total of 123 trees were sampled for dendrochronological purposes ( Table 1) and their diameter at breast height (DBH) was measured in the field. The cores were prepared following standard dendrochronological procedures. Cores were air dried, glued onto wooden supports and sanded until the tree rings were clearly visible. The processed samples were visually cross dated (Fritts, 1976) and the ring widths were measured to a 0.001 mm resolution using scanned images (Epson Expression 10,000 XL, resolution 2,400 dpi) and the CDendro software (Larsson and Larsson, 2018). The age of each sampled tree at 1.3 m was estimated by counting the maximum number of rings in the two cores extracted per tree. The COFECHA software (Holmes, 1983) was used to validate the cross dating by calculating moving correlations between the individual series and the mean series of each species.
First, to study growth responses to climate at the species level after removing biological growth changes, we detrended tree-ring series to eliminate the influence of changes in tree size and age or disturbances (Figure 2). Each series was detrended with a cubic smoothing spline that has a 50% cut-off response in 20 years. We used the same detrending for all samples to facilitate the comparison of the chronologies and their response to climate (Klesse, 2021). Individual detrended series were also subjected to pre-whitening to remove temporal autocorrelation and then averaged using bi-weight robust averages to obtain residual chronologies. To assess within-species growth coherence, we calculated the mean correlation between individual series (Rbar) and the Expressed Population Signal (EPS), which is a measure of the robustness of the chronology (Wigley et al., 1984).
To study growth trends and responses to climate at the individual tree level, we converted tree-ring widths measures into basal area increments (BAI; Figure 2), assuming a circular shape of stems and using the equation: where r t and r t−1 are the tree radius in the year of tree-ring formation (t) and the year before tree-ring formation (t-1), respectively. BAI adequately reflects growth changes in the early stages of tree life and retains high-frequency variations. BAI was calculated for each individual tree after averaging the tree-ring width measures of the two series from the outside inwards.
Pointer year analysis (Schweingruber, 1986) was used to identify and quantify the presence of abrupt growth enhancements and reductions for each species. Analyses were done at the species level considering all trees of each species (tree level not series level) and using BAI. For each tree, growth enhancements and decreases were quantified as growth increase of 40% and decreases of 60% as compared to the three preceding years, respectively. According to Schweingruber (1986) pointer years are those growth enhancements or decreases that occur in 75% or more of the individual series.
Resilience growth components were quantified according to Lloret et al. (2011). Further, we considered two additional metrics proposed by Thurm et al. (2016): recovery period and total growth reduction. Resistance and Resilience were quantified for the 2005-2007 severe drought ( Supplementary  Figures 2, 3). Resistance (Rt) was considered as the ratio in growth (BAI) between the first dry year (2005) and the three preceding years. Resilience was quantified as the ratio in growth between the pre- (2002-2004) and post-drought (2008-2010) periods. Both the recovery period and total growth reduction were calculated for a maximum period of 11 years after 2007.
Finally, early warning signals of critical transitions were quantified by quantifying the coefficient of variation (CV) and the first-order autocorrelation (AR1) of the ring-width indices (RWI) of each species considering 25-year moving windows shifted by 1year intervals (Dakos et al., 2012). For these analyses, we used the standard chronology (i.e., before pre-whitening) of each species as it maintains short-term temporal autocorrelation but has no effect of long-term biological trends (in contrast to BAI series).

Statistical Analyses
Correlation analyses were used to test for the relationship between ring-width indices and, temperature, precipitation, and drought (SPEI) for the period 1980-2017. Bootstrapped correlation analyses were performed considering the period 1980-2017. Significance was assessed by resampling 1,000 times using the classical bootstrapping method (Zang and Biondi, 2015). Analyses were performed at the species level using the residual chronologies for each species (i.e., after pre-whitening) and the mean temperature, precipitation and 12-month long SPEI from the month of September in the year before growth to the month of September in the year of growth.
We used linear mixed-effect models (LMM; Pinheiro and Bates, 2000) to study growth trajectories and response to climate for each species at the tree level. The model was of the form: Where Y is the log-transformed BAI (log(x+1)) of each individual s in the year t, f(X) is the set of fixed effects, us represents the tree nested within stands random effects, vt is a normally distributed random effect for calendar year t (year effect), and est is the normally distributed residual for tree s at year t. We used this model following Mehtätalo et al. (2014) to study the temporal variation in BAI while accounting for the fixed effects together with unspecified tree-, and year-level factors.
As fixed factors, we included the age trend, the drought index and species identity. The age trend was represented by the log-transformed age of each tree s at year t. To represent drought conditions, we selected the 12-month long SPEI for the month of September (Supplementary Figure 2). We included potential interactions between both variables and species identity. The analyses were carried out for the period common to all species . Fixed factors were standardized to have zero mean and unit variance prior to the analyses. Model selection was based on the Akaike Information Criterion (AIC) (Burnham and Anderson, 2002). All potential models combining the variables listed above were created and the one with the lowest AIC was selected. The model was evaluated by graphically inspecting the residuals (Zuur et al., 2009). For the selected model, we calculated the conditional (R 2 c; variance explained by fixed plus random effects) and marginal R 2 (R 2 m; variance explained by fixed effects) coefficients of determinations according to Nakagawa et al. (2017).
We studied if growth trends after removing age-and droughteffects differed between species before and after the dry period. To this end, we obtained the residuals of the models proposed above and studied how they varied between species and periods using LME's. The residuals were used as response variable and species identity, calendar year and their interaction were used as fixed factors using tree nested within stand as random factor. Separate models were created for the period 1980-2004 and 2008-2017. To test for the significance of year, species and / or their interaction, we run model selection based on AIC.
We used linear models to compare the resilience components between species. The corresponding indices (Rt, and Rs) were used as response variable and species identity, tree age and DBH at sampling time were used as covariates. The resilience components were log-transformed prior to the analyses (log(x+1)). In the case of the recovery period and the absolute growth reductions, those trees that did not recover pre-drought growth levels were not analyzed. In all cases, model selection was applied and the models with the lowest AIC were selected. When significant effect of species or interactions between species and year were found, we performed least-square means based on the Tukey's honest significant difference (HSD) to test for the differences in growth trends and resilience components between species (Bretz et al., 2010).
All analyses were performed in the R environment for statistical computing (R Core Team, 2021) using Rstudio (RStudio Team, 2021). The package dplR (Bunn, 2008;Bunn et al., 2021) was used to manage tree-ring width series, detrend them and calculate BAI; the package pointRes (van der Maaten-Theunissen et al., 2015) was used to calculate pointer years and resilience components; climate-growth analyses were done with the treeClim package (Zang and Biondi, 2015); the package early warnings (Dakos et al., 2012) was used to quantify early warning signals for critical transition; the package lme4 (Bates et al., 2015) was employed to create linear models and LMMs, the package MuMIn to perform AIC model selection (Barton, 2018) and the package emmeans (Length, 2022) was used to perform trends and factor comparisons using estimated marginal means; and finally the package effects (Fox and Weisberg, 2019) was used to visualize regression graphs.

Growth Trends and Response to Climate
Mean chronologies (pre-whitened RWI) were similar between species (Figure 2) suggesting that the four species are driven by similar climate conditions (Figure 3). In this line, correlations between residual chronologies and climate showed that P. halepensis, P. pinea, and P. sylvestris growth was enhanced by March precipitation (Figure 3). All species except P. pinea responded positively to the precipitation in June and July and P. sylvestris responded also positively to the precipitation in May. High temperatures during March reduced the growth P. nigra and P. sylvestris (Figure 3). Finally, all species showed positive significant correlations with the 12-month long SPEI from the months of May to September (Figure 3).
The models of BAI at the tree level also supported the strong impact of the SPEI on growth ( Table 2 and Supplementary  Table 1). The selected model accounted for 69% of the variation in the data (conditional pseudo-R 2 ) which the fixed factors accounting for 33% of it (marginal pseudo-R 2 ). BAI increased as trees aged, and it varied between species (Figure 4). The relationship between BAI and SPEI was the strongest for P. nigra and the lowest for P. sylvestris ( Table 2). In general, P halepensis  For each variable and their interactions, the estimated regression coefficient and standard error are shown together with the t-value and associated probability. Species' codes: PN, P. nigra; PP, P. pinea; PS, P. sylvestris.
had a higher age-corrected BAI than the rest of species (i.e., fitted values for individuals with age fixed at its mean across species, Figure 4). The residual variation in BAI after accounting for age-and SPEI-effects showed different results in the two periods studied. In the first period , we found no differences in the residual trend between species (Supplementary Table 2). Conversely, we found differences in the trends between species in the second period (2008-2017; Figure 4 and Supplementary Table 2). In this second period, the model accounted for 56% of the variation in the data although only 2% was because of the interaction between calendar year and species. The growth trends differed between species markedly because while P. halepensis and P. pinea had positive trends, P. sylvestris and P. nigra had clearly negative trends (Figure 4 and Supplementary Table 4).

Resilience Components
The pointer year analyses showed similar patterns in the occurrence of positive and negative pointer years between species although with subtle differences (Figure 5). Growth enhancements higher than 40% and occurring in 75% or more trees were only recorded by P. halepensis in 2004. Growth decreases higher than 60% and occurring in at least 75% of trees were observed in P. pinea and P. halepensis in the years 2005, and P. nigra in the year 2007 (Figures 2, 5). P. nigra also displayed important growth decreases in 1998-1999 and 2016. The strongest growth reductions in P. sylvestris were observed in 2007 and particularly in 2016. The 2004-2007 period was characterized also by substantial growth enhancements in 2004 affecting all species, except P. sylvestris, and followed by subtle and severe growth reductions in 2005 and 2007.
The model of Rt accounted for 31 % of its variation and showed differences between species as well as effects of DBH and tree age (Supplementary Figure 4 and Table 3 and  Supplementary Table 3 Table 5). The models of resilience (Rs) showed small differences between species being P. halepensis and P. pinea the most resilient species and P. nigra and P. sylvestris the less resilient (Figure 5 and Table 3). These differences were not significant according to the post-hoc analyses (Supplementary Table 5). The models accounted for 6% of the variation in Rs. No differences between species were found in the average time to recover and the absolute growth reduction (Supplementary Table 3). However, there were 24 trees that were not able to recover pre-drought BAI levels and corresponded to 9 P. nigra (26%), 7 P. sylvestris (27%), 5 P. pinea (12%), and 3 P. halepensis (12%) individuals.

Early Warning Signals for Critical Transitions
Early warning analyses showed differences in the first-order autocorrelation and the coefficient of variation (CV) of standard chronologies (i.e., RWI) between species (Figure 6). In general, P. sylvestris displayed the highest autocorrelation and the lowest coefficient of variation in RWI while P. nigra displayed the highest coefficient of variation in RWI followed by P. halepensis. P. pinea displayed a marked drop in autocorrelation coinciding with the 2005 drought (less intense in P. halepensis) and had a positive trend in autocorrelation after it (Figure 6). P. halepensis displayed an autocorrelation close to zero, excluding the effect of the 2005 drought, during most of the studied period. The trend in autocorrelation and the coefficient of variation of standard chronologies was strong and positive in P. sylvestris and P. pinea according to the non-parametric Kendall tau correlation (Figure 6). In the case of P. nigra the trend was strong for autocorrelation only while no clear trends were observed for P. halepensis.

DISCUSSION
The results presented in this study corroborates that drought is a major factor determining the growth of pines in Mediterranean forests (Pasho et al., 2011;Peñuelas and Sardans, 2021). What is novel here is that we have been able to compare how four coexisting pine species with different growth behavior and biogeographical origins respond to drought when coexisting. Previous studies have compared how two (Martín-Benito et al., 2013;Marqués et al., 2016) or even three of these species responded to drought (Granda et al., 2018), or used some of them to understand their physiological responses to drought in controlled conditions (Andivia et al., 2019). Here we show how they vary in their growth response to drought and post-drought resilience capacity. Our results indicate that an extreme drought event, such as that occurring in 2005−2007 (Supplementary Figures 2, 3), differently impacted each tree species and had disparate effects on their growth. The growth of the Mediterranean P. pinea and P. halepensis was more impacted by the first drought event (2005) while P. sylvestris and P. nigra growth was more impacted by the second part of the dry period (2007) (see Figures 2, 4, 5). In this sense, P. sylvestris and P. nigra were more resistant to the 2005 drought than P. pinea and P. halepensis. However, P. sylvestris and P. nigra trees were slightly FIGURE 4 | Expected and observed growth (BAI) for each pine species. Panels (A,C,E,G) show the predicted growth pattern (solid and dashed lines represent main effect and 95% confidence intervals, respectively) including the SPEI*species interaction while keeping age constant (mean age standardized = 0; ∼30 years). Panels (B,D,F,H) show the observed growth trajectories considering the common period 1980-2017. Gray lines represent the observations for each individual tree. Straight lines in the period 2008-2017 (secondary y-axis) represent the expected growth trends after accounting for age-and SPEI-effects (solid and dashed lines represent main effect and 95% confidence intervals, respectively). The positive trend in the residuals means that the model underpredicts observed growth trends during the studied period.
less resilient, presented negative growth trends after accounting for age-and SPEI-effects and were less able to retrieve predrought growth levels (Figures 4, 5). The analyses of climategrowth couplings revealed the sensitivity of all species to drought (i.e., SPEI; Figure 3C) but also suggested a greater sensitivity of P. sylvestris and P. nigra to summer precipitation (particularly, P. sylvestris, Figure 3A). These results at the species level contrasted with the results of BAI at the tree level that showed a lower sensitivity of P. sylvestris to drought. Early warning signals indicated differences between species and an important impact of the 2005 drought in P. pinea growth autocorrelation (Figure 6). All together, these results suggest that all species are vulnerable to drought occurrence but that P. halepensis might be in an advantageous position to thrive in a drier world and outperform the other species, at least in the study site or in other sites with similar climate conditions. It has been suggested that negative growth trends are common before drought-induced dieback in conifer species (e.g., Camarero et al., 2015aCamarero et al., , 2018Camarero et al., , 2021aCailleret et al., 2017;De Soto et al., 2020;Keen et al., 2022). We observed that growth trends for P. sylvestris and P. nigra were negative in the period 2008−2017 after accounting for age-and SPEI-effects (Figure 4). Moreover, the analyses of the resilience components suggested that 9 out of 34 (26%) and 7 out of 26 (27%) P. nigra and P. sylvestris trees were not able to recover pre-drought BAI levels. These age-corrected growth trends contrasted with the positive trends found in P. halepensis and P. pinea. According to this, we would interpret that the more Mediterranean species (P. halepensis and P. pinea) have a better capacity to thrive with the current climate conditions than the other two species (P. nigra and P. sylvestris). Thus, it can be interpreted that P. halepensis and P. pinea have a better capacity to recover after the 2005-2007 dry period than P. nigra and P. sylvestris in the studied forests. However, the analyses of the resilience components were not so clear (Figure 5 and Supplementary Figure 4). First, we found a higher resistance to the 2005 drought of P. nigra and P. sylvestris and differences in resilience between species were weak. Tree species differ in their resistance to hydraulic damage, which can make some species more resilient against drought than others (e.g., Choat et al., 2018). Across biomes, these species vary in their growth resilience capacity to drought being drought severity one of the most important factors explaining such variations but displaying similar trade-offs in resistance and recovery between species (e.g., Gazol et al., 2018). However, when different tree species co-habit in a region, and are thus exposed to comparable drought conditions, only species differences, microsite conditions and trees' features can determine growth responses to drought (Baker et al., 2019;Gazol et al., 2020;Férriz et al., 2021). Overall, the four species were resilient to the drought event (2005)(2006)(2007), but the negative growth trend of P. sylvestris and P. nigra may suggest their greater vulnerability to drought and their lower capacity to recover.
The capacity of P. halepensis and P. pinea species to present facultative bimodal growth patterns (e.g., Camarero et al., 2010;Pacheco et al., 2016;Castagneri et al., 2018;Campelo et al., 2021) can explain their resilience to drought, which contrast to the low capacity of P. sylvestris to resume growth (Oberhuber et al., 2021). The occurrence of drought often results in low photosynthesis rates in gymnosperms due to stomatal closure (Forner et al., 2018), which may give a competitive advantage to those species that maintain growth rates during dry periods (e.g. Salazar-Tortosa et al., 2018). In this sense, the weaker stomatal control of P. halepensis and P. pinea (Manzanera et al., 2016) can help these species to maintain growth and tolerate dry periods better than P. nigra or P. sylvestris. Another potential explanation for the greater resilience capacity of P. halepensis and P. pinea can be found in their greater capacity to invest more in belowground organs in response to drought (Sarris and Mazza, 2021). It has been suggested that the capacity to access bedrock water may serve as a "hydraulic refugia" allowing trees to avoid drought damage (McDowell et al., 2019). Andivia et al. (2019) found that seedlings of these species grew belowground organs faster than those of P. sylvestris and P. nigra allowing them to access deeper water sources. In fact, Férriz et al. (2021) found that larger P. pinea trees were less vulnerable to drought than small individuals and pointed to relationship between aboveground size and the root system. However, it is important to consider that P. sylvestris can be more vulnerable to pathogens than P. halepensis and P. pinea in the studied forest. Caballol et al. (2022) found that this species was more affected by Diplodia shoot blight than the rest of Pine species. We do not have information on previous hailstorms and subsequent Diplodia outbreaks that may explain P. sylvestris declining trends, but the greater sensitivity of the species to this pathogen and its low resilience capacity points to its vulnerability to climate change in sight of the importance of combined impacts of drought and pathogen attacks (McDowell et al., 2008(McDowell et al., , 2019Trugman et al., 2021).
Another aspect that deserves further attention is the differences in the response of growth to the dry period across species (Figures 2, 3): while Mediterranean pines responded more to the first part of the dry period (the year 2005; Figure 5A), pines from mesic sites responded more to the second part (the year 2007). This is particularly true in the case of P. pinea which shows a very marked growth reduction in response to the 2005 drought that has also marked effects on growth autocorrelation (Figure 6). Interestingly, Férriz et al. (2021) found that the growth of declining and non-declining P. pinea trees in central Spain started diverging after the 2005 drought. The growth and xylogenesis of P. pinea is strongly determined by water availability during the growing season and before (Castagneri et al., 2018;Calama et al., 2019;Mechergui et al., 2021). In fact, Piraino (2020) found that planted P. pinea forests in central Italy were sensitive to spring rather than to summer drought. This is also true for P. halepensis which tends to respond to drought at large temporal scales (Pasho et al., 2011). In our study, P. pinea was the species with a strongest response to the precipitation in March and showed significant correlations with the 12-month SPEI from January to September. In P. halepensis we also observed a significant response to the precipitation in March, but its growth was also dependent on summer precipitation (Figure 3). It has been argued that P. nigra growth is enhanced by May precipitation and reduced by June temperatures (Martín-Benito et al., 2008;Sánchez-Salguero et al., 2012) as well P. sylvestris (Andreu et al., 2007) and they respond to drought at shorter time scales. Our results partially confirm these findings as we found negative impacts of March and May temperature in their growth and a strong linkage with summer precipitation. This was particularly true in the case of P. sylvestris. However, P. sylvestris was the species displaying the least intense relationship with SPEI at the tree level. It should be noted that P. sylvestris displayed the lowest correlation between series (Table 1) which can explain their lower correlation with SPEI at the tree level because of the different response of different trees. Thus, it is difficult to explain why these species responded differently to the 2005 and 2007 years ( Figure 5). Alternatively, the capacity of P. halepensis and P. pinea to maximize carbon uptake during drought (Borghetti et al., 1998;Salazar-Tortosa et al., 2018) can explain their ability to cope with the 2007 dry year. In any case, these results evidence that drought timing and its interaction with species phenology are important factors controlling tree growth responses to drought (Camarero et al., 2015b).
Both, P. halepensis and P. pinea are Mediterranean pines species adapted to cope with summer drought (Castagneri et al., 2018;Campelo et al., 2021). In our study, differences between these two species in drought resistance and resilience were minimal as well as in growth trends after accounting for ageand SPEI-effects. We only found small differences suggesting a greater sensitivity of P. pinea to the 2005 drought in the form of a sudden drop in its growth autocorrelation (Figure 6). In addition, we found clear differences in projected growth trajectories between species with P. halepensis showing greater growth than P. pinea (Figure 4). However, important age differences between species makes difficult to assure that these results indicate that P. halepensis tolerates better current drought conditions in the region than P. pinea. Interestingly, Caballol et al. (2022) found linkages between P. pinea defoliation and growth, indicating that those trees showing higher defoliation values grew less in the past. This may indicate that this species is suffering from drought-induced dieback and points to its drought vulnerability in the studied forests. A recent review indicates that P. pinea is a species very sensitive to climate change whose habitat will suffer strong reductions in a drier future (Mechergui et al., 2021). It is important to note that we are not considering other aspects that can influence the persistence of the species in the forests such as regeneration (e.g., Férriz et al., 2021). Further studies should investigate whether P. pinea perform worse than P. halepensis in the studied area testing potential mechanistic explanations such as a higher sensitivity to warmer conditions and increased vapor pressure deficit (Mechergui et al., 2021) or different root systems (Sarris and Mazza, 2021).
Drought-induced growth declines and canopy dieback episodes of P. sylvestris been widely reported since the late 20th century across southern Europe (e.g., Hereş et al., 2014;Caudullo and Barredo, 2019) and process-based model projections forecast a worsening of the growth conditions for the species in the region (Sánchez-Salguero et al., 2017). Despite being considered a drought-tolerant species, drought-induced dieback episodes have been also observed in P. nigra (Sánchez-Salguero et al., 2012;Camarero et al., 2018), and Circum-Mediterranean studies suggest an enhanced vulnerability of its growth to drought . Even in the case of the droughtadapted Mediterranean P. pinea and P. halepensis droughtinduced growth declines, dieback and mortality episodes have been observed in some regions after very warm and dry conditions (e.g., Camarero et al., 2015a;Mechergui et al., 2021). These observations indicate that all the study pine species and other Mediterranean species such as P. pinaster (Gea-Izquierdo et al., 2019), are potentially vulnerable to dry spells, but P. halepensis will be the species performing better in terms of growth when coexisting with P. pinea, P. nigra, and P. sylvestris. However, these results should be interpreted with caution due to age differences between sampled individuals of each species.

CONCLUSION
The climatic characteristics of the Oristà forest facilitated the coexistence of four pine species with contrasting drought resistance and resilience. Results indicate that P. halepensis has a greater capacity to recover pre-drought growth levels after drought than P. sylvestris and P. nigra when coexisting. We found that P. halepensis grew better and is slightly less sensitive to drought than P. pinea in the studied forest. Thus, growth patterns and response to climate suggest that P. halepensis has a greater capacity to thrive in a drier and warmer world than the other three species studied. Further studies could address if this capacity of P. halepensis to withstand drought as compared to other pines is based on above-(e.g., stomatal conductance, regulation of hydraulic conductivity through needles and stems) or belowground (e.g., root architecture and conductivity, fine root longevity, sources of soil water) adjustments.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JO designed the study. MC, JC, and CV collected and processed the wood samples. AG analyzed the data. AG and JC drafted the manuscript and all authors commented on it. All authors contributed to the article and approved the submitted version.

ACKNOWLEDGMENTS
We thank Maria Caballol and Maia Ridley for their help in the field.