The study was carried out in the arboretum of the Iturrarán botanical garden (hereafter Iturrarán, 43.249°N, 2.155°W, 70 m a.s.l., Gipuzkoa province, Basque country, Spain).¹ The botanical garden is about 25 ha large, and its arboretum holds a living collection of around 300 oak species (Quercus spp.) with several mature individuals per species (typically 3–5 individuals). Previous studies have been done in this site to evaluate the interspecific variation in oak leaf morphology and hydraulic traits (Sancho-Knapik et al., 2021; Martín-Sánchez et al., 2024a,b). The site elevation varies between 120 and 210 m a.s.l. The soils are mainly composed of clay and their pH is neutral. The climate conditions are softened due to the influence of the Atlantic Ocean. The elevated precipitation of the site, between 1,500 and 1,800 mm per year, and mild temperature conditions (range of temperatures 5.7–22.5 °C) lead to oceanic climate conditions (data from Hernani station, 43.265°N, 1.976°W, 44 m a.s.l.). The potential tree species of the study area is Quercus robur L.

Sampling was carried out in May 2025 with the invaluable help of the former arboretum manager (Francisco Garin). We selected a total of 20 oak species represented for at least three adult individuals. Most sampled species were winter-deciduous (n = 14), and the rest were evergreen species (n = 6) (Supplementary Table S1). We measured the diameter at breast height (dbh) of each selected tree at 1.3 m and extracted two cores perpendicular to the maximum slope of the terrain using 5-mm-wide Pressler increment borers (Haglof, Sweden).

The extracted cores were processed and analyzed using dendrochronological methods including cross-dating and detrending of ring-width series (Fritts, 1976). Briefly, cores were air dried after collection, glued onto wooden supports, and sanded with progressively finer-grain sandpaper until ring boundaries were clearly visible. After that, we used a high-resolution scanner (Epson 13000 XL) to scan the processed wood samples at 2,400 dpi resolution, and the resulting images were used to measure ring widths (0.01 mm resolution) and cross-date the samples using the CooRecorder-CDendro software (Larsson and Larsson, 2022). The last complete ring we measured corresponded to the year 2024. The quality of cross-dating was checked using the COFECHA software. The cross-dating was done by calculating moving correlations between individual series of ring-width values and the mean series of the site (Holmes, 1983).

To evaluate trends in radial increment and compare growth trajectories between species, we converted ring-width series into mean basal area increment (BAI) series assuming a concentric stem growth (Biondi and Qeadan, 2008). To calculate BAI we used the formula:

where RWLt and RWLt-1 are the ring-width in the years t and t − 1, respectively. To evaluate year-to-year variations in radial increments we transformed ring-width measurements into dimensionless ring-width indices (RWI). Individual ring width measures were detrended using a smoothing spline with a frequency response of 0.5 and a wavelength of 15 years. Thus, BAI series contain low- and high-frequency variations, whilst RWI series mainly contain high-frequency variations. Thus, we used BAI as a measure of radial increment magnitude and trajectory, while RWI was used to study the relationship between climate variables, including a drought index, and growth variability.

For each species, we quantified the mean ring width and its standard deviation, the Gini coefficient and the first-order autocorrelation (AR1), which measure inequality and year-to-year similarity, respectively. These statistics were obtained by averaging the individual measures for each species (see Supplementary Table S1).

As mentioned above, the site has been previously used as a quasi-common garden to study interspecific differences in leaf and hydraulic traits (Sancho-Knapik et al., 2021; Martín-Sánchez et al., 2024a). Thus, we take advantage of the data published in these studies to study the relationship between xylem and leaf traits and growth variables at the interspecific level. Particularly, Martín-Sánchez et al. (2024a,b) quantified the interspecific variation in leaf area (LA, in cm²), leaf dry weight (LW) and dry leaf mass per area (LMA, in g m⁻²) for 141 oak species in Iturrarán, and 85 of these samples were previously studied by Sancho-Knapik et al. (2021). Small leaves with high LMA, representative of slow resource use strategy (Reich, 2014) are more prevalent in oaks from arid regions while larger leaves with low LMA appear in cold-wet temperate regions (Martín-Sánchez et al., 2024a). We gathered the values of LA, LW and LMA that accompanied the publication of Martín-Sánchez et al. (2024a) and were publicly available in Zenodo (Martín-Sánchez et al., 2024b). Note that data was available for 19 out of the 20 species studied in our study since Quercus phellos L. was not considered by Martín-Sánchez et al. (2024a). We also obtained the Index of Lobation (ILB) from Martín-Sánchez et al. (2024a,b) because this morphological trait, used to measure the degree of indentation of leaf margins, also characterized oak tolerance to heat stress since narrow and deeply lobated leaves are better adapted to warmer conditions climates due to thinner boundary layers and enhanced hydraulic efficiency (Sisó et al., 2001).

We also obtained xylem traits, including wood density, of the selected species considering anatomical traits related to mechanical strength, sap transport and xylem vulnerability to embolism (Fontes et al., 2022). Specifically, we gathered data on hydraulic diameter (Dh), percentage conductive area (Cond.A) and the Huber Value (HV, i.e., conductive area per leaf area ratio). Stem wood density was measured using the water-displacement method for volume and oven drying at 101 °C for 72 h for mass (Fajardo, 2025). This was done in the laboratory using an additional core extracted at 1.3 m from 2 to 4 individuals per species. In the case of xylem variables, we obtained three to five 1-year-old shoots from different individuals of each species and cut 20-to-30-μm-thick cross sections using a sledge microtome (AS 200; Anglia Scientific Instruments, Cambridge, UK). Sections were mounted on glass slides, stained with safranin (0.5% in distilled water), Astra Blue (2%) and fixed with Eukitt®. Images of sections were taken at 20–40× magnification with a digital camera mounted on a light microscope (Olympus BH2). They were stitched with the ICE software (Microsoft©) and analyzed using ImageJ-Fiji4 (Schindelin et al., 2012). The same image analyses software as used to measure the total leaf area in order to calculate HV following García Saganta (2006). Lastly, the Dh was calculated as.

where D is the mean diameter (ellipses) of the n vessels measured in the whole ring (Scholz et al., 2013).

We gathered the species occurrence data from Martín-Sánchez et al. (2024a) dataset, since they evaluated how macroclimate conditions in species distribution range affects leaf trait variation, and obtained occurrence data for each species from GBIF.² We complemented the dataset by including the geographical distribution coordinates of Quercus phellos L. from GBIF (see text footnote 1; 1 October 2025) (GBIF.org, 2025). The coordinates associated with the occurrence of each species were used to download bioclimatic variables (Supplementary Table S2) as done by Martín-Sánchez et al. (2024a) (Figure 1).

Finally, we used the previous study done by Bose et al. (2024) to obtain the growth characteristics and growth response to drought of oak species in their regions. We retained the values in that study to check the coherence between them and the findings presented here. Only eight out of the 20 species studied were shared with previous dendrochronological records and allowed comparing resilience measured in Iturrarán and drought legacies in the field at large scale. Bose et al. (2024) quantified drought legacies as the departure of observed tree growth from expected tree growth based on the relationship between mean site series of ring-width indices and monthly climate variables. This is methodologically different but conceptually similar to the resilience components quantified in this study. Resilience measured in the field at large scale.

We quantified the radial increment responses of each species to drought using resilience indices (Lloret et al., 2011; Thurm et al., 2016; Schwarz et al., 2020). We focused on five indices that account for the magnitude of radial increment reductions following drought and the capacity and time required to recover after drought. We focused on the 2022 drought as it was the most intense (see Figure 2), and calculations were based on BAI and considering 2 years before and after drought occurrence as pre- and post-drought periods, respectively. We calculated: (i) drought resistance (Rt), which is the ratio between drought and pre-drought BAI; (ii) recovery (Rc), i.e., the ratio between post-drought and drought BAI; and (iii) resilience (Rs), which is the ratio between post- and pre-drought BAI (Lloret et al., 2011). Additionally, we also considered the recovery period (Rp) and the relative growth reduction (Gr) which accounts for the time required to recover pre-drought BAI (in years) and the percentage of BAI loss during the drought and recovery period (Thurm et al., 2016).

The average climate conditions and additional bioclimatic variables (19 bioclimatic variables; see Supplementary Table S2) were obtained for each site from WorldClim database.³ The bioclimatic variables are derived from monthly temperature and precipitation records and provide information on climate annual trends, seasonality, and extremes (Fick and Hijmans, 2017; Hijmans et al., 2005). We used the moderate resolution global data (2.5 min) to represent the bioclimatic conditions for each species occurrence. Additionally, mean values were calculated for each bioclimatic variable across species to represent the bioclimatic distribution for the species (Figure 1).

To calculated climate-growth relationships we used the 4-km gridded Terraclimate dataset (Abatzoglou et al., 2018). We downloaded the following monthly climate variables: mean maximum and minimum temperatures, precipitation, vapor pressure deficit and downward radiation. These variables were obtained for 2000–2024, a period which included most tree-ring series, but were available since 1958 (Supplementary Figure S1). We also used the Standardized Precipitation Evapotranspiration Index (SPEI) to represent drought severity in the same period. The SPEI is a multi-scalar drought index that measures drought severity through the quantification of its intensity and duration (Vicente-Serrano et al., 2010). The SPEI is computed from the difference between potential evapotranspiration (PET) and precipitation and thus it accounts for the effect of both temperature and precipitation on drought (Vicente-Serrano et al., 2015). Positive and negative values of the SPEI correspond to cumulative water surplus and water deficit, respectively (Vicente-Serrano et al., 2010). We obtained weekly SPEI data at 1-, 3-, 6-, 9-, 12-, and 24-month long temporal scales for 2000–2024 from the 1.1-km² Spanish SPEI dataset (Vicente-Serrano et al., 2017) to analyze RWI-SPEI relationships.

The relationship between radial increment variability (RWI), climate variables and drought (SPEI) were studied by means of Pearson correlation (r) analyses considering the best-replicated period (2000–2024). A mean site chronology was calculated using the individual chronologies for all the species (20 species, 62 radii). We studied the relationship between RWI, climate variables and SPEI at different scales and for different months. In the case of climate variables, correlations were calculated from the previous October to the current September.

We used Pearson correlation, Principal Component Analysis (PCA) and path analysis to study the relationship between leaf (LA, LMA, LW, ILB), xylem traits (Dh, Cond.A, HV), wood density (WD), and size or radial growth traits (diameter, mean ring width, AR1, Gini coefficient) and resilience components (Rt, Rc, Rs, Gr, Rp). Correlation analyses were used to study bivariate relationships between all the variables. We used Principal Component Analysis (PCA) to decipher major axes of size and radial growth variables, and resilience components variation and the coordination between them (Legendre and Legendre, 2012). Finally, linear regression analyses and path analysis were used to study how climate variables and functional traits (leaf, xylem and wood traits) affect the interspecific variation in resilience components (Rt, Rc, Rs, Rp and Gr). In this sense, we tested how resilience was related to climate-of-origin, and whether this relationship was mediated by functional traits. Due to the reduced sampled size, we restricted the path analyses to the climate variables with the strongest correlation with resilience (Rs) and the functional trait significantly correlated to Rs. We proposed two path models: (1) climate-of-origin affects Rs through its effect on the functional trait; (2) as the previous one but including a direct effect of climate-of-origin. The two models were compared based on their Akaike Information Criterion (AIC).

All statistical analyses were done in the R statistical environment (R Core Team, 2025). The BIOMEplot package was used to create the Whittaker biomes figure (Kunstler, 2014). The package “dplR” (Bunn, 2008, 2010; Bunn et al., 2025) was used to calculate dendrochronological statistics, quantify BAI and detrend the ring-width series. The geodata (Hijmans et al., 2024) and raster (Hijmans, 2024) packages were used to download bioclimatic data and extract information for species occurrences. The pointRes (van der Maaten-Theunissen et al., 2015, 2021), corrplot (Wei and Simko, 2024), vegan (Oksanen et al., 2025), lavaan (Rosseel, 2012) and visreg (Breheny and Burchett, 2017) R packages were used to quantify resilience indices, perform statistical analyses and create graphs.

Q. shumardii and Q. rubra displayed the largest diameter values while Q. phillyraeoides was the species with the lowest average diameter (Supplementary Table S1). All ring-width series covered the period 2007–2014 (Supplementary Figure S2). The BAI trajectory for all oak species together showed an increase until 2011 and a subsequent stabilization (Supplementary Figure S2). A marked BAI drop was observed in the dry years 2013 and 2022 (Figure 2). The mean ring-width index for the year 2022 was the lowest in the period 2007–2024 (Figure 2a), and summer maximum temperatures were extremely warm during that year (Supplementary Figure S1).

The growth index showed a significant positive relationship with the drought index (SPEI) at short temporal scales (Table 1). In particular, the 1-month SPEI from mid-July to early-August was strongly related with RWI (r ≥ 0.5) as well as the 3-month SPEI from late-September to early-October (Table 1). The strongest relationship with RWI was found for the 1-month SPEI in the first week of August (r = 0.700; p = 0.01) and for the 3-month SPEI in the third week of September (r = 0.775; p = 0.01). Oak growth was enhanced by warm and sunny April conditions, whereas elevated temperatures in the previous December and the current August reduced growth (Table 2). Sunny July-to-August conditions were also related to lower growth indices.

The first two axes of the PCA accounted for 56% of the variation in radial growth traits (diameter, mean ring width, AR1, Gini coefficient) and resilience components (Figure 3). The first axis separated species with large diameter at breast height and wide rings from (e.g., Q. pagoda) those species showing lower growth rates but presenting high values of the Gini coefficient (growth variability) and first-order autocorrelation (e.g., Q. acutissima). Accordingly, the tree diameter was positively related with the mean and the median ring width and also with the ring-width standard deviation (Figure 4). The second axis was mainly related with the separation of species with a high ring-width kurtosis coefficient and Rc (e.g., Q. phillyraeoides) and those with high values of standard deviation of ring-width (e.g., Q. palustris).

The Gini coefficient was negatively related with mean and median growth and positively with the standard deviation in mean growth rate and the Gini coefficient. The 2022 recovery (Rc) and resilience (Rs) indices were positively correlated between them. The 2022 recovery period (Rp) and the growth reduction (Gr) were also positively correlated, and negatively correlated with resistance (Rt) and resilience (Rs), indicating that larger growth reductions need more time to recover pre-drought radial increment.

We found that radial growth traits were related with leaf and wood traits (Supplementary Table S3). Particularly, tree diameter and the standard deviation in ring-width were negatively related to LMA. Similarly, the Dh was negatively related with mean ring-width, whereas the HV and the relative growth reduction (Gr22) presented a positive correlation. A negative correlation was found between Rs and LW.

Regarding functional traits, leaf area and LMA were negatively correlated, but LMA was positively correlated with the HV (Supplementary Figure S3). In contrast, the ILB was negatively correlated with LMA and the HV. We found that radial growth traits were related with leaf and wood traits (Supplementary Table S3). Particularly, tree diameter, mean ring-width and its standard deviation were negatively related to LMA. Similarly, the Dh and conductive area and the WD were negatively related with mean ring-width, and its standard deviation, whereas the HV and the relative growth reduction (Gr22) presented a negative correlation.

The resilience components varied between species, but they were rarely related to functional traits. Tree diameter was positively correlated with growth rates and their standard deviation (Figure 4). Interestingly, we found that the 2022 recovery (Rc) and resilience (Rs) indices were positively correlated between them. The 2022 recovery period (Rp) and the growth reduction (Gr) were also positively correlated indicating that larger growth reductions need more time to recover pre-drought radial increment.

Tree diameter and growth rates showed positive correlations with BIO2, the mean diurnal range, but tree diameter was negative correlated with BIO12, BIO16 and BIO18, i.e., tree size increased with wetter conditions (Supplementary Table S4). However, resilience (Rs) decreased as climate conditions were wetter, but Rp and Gr increased. Resilience displayed the strongest negative correlation with the precipitation of the wettest quarter (BIO16) while Rp and Gr were positively related to this variable. Rc and Rs showed negative correlations with the minimum temperatures of the coldest month (BIO6). Regarding functional traits, LW, LMA and HV increased in response to temperate or mesic (e.g., higher BIO3) and wetter conditions (e.g., higher BIO13 and BIO16) (Supplementary Table S5), whereas LA and WD decreased. Additionally, LA increased in climates with wider thermal ranges and colder winters. Conductive area was negatively correlated with the mean diurnal range (BIO2) and Dh with the mean temperature of driest quarter (BIO9).

The minimum temperature of the coldest month (BIO6) was negatively correlated with the recovery (Rc) and resilience (Rs) indices (Table 3). The temperature annual range showed a negative association with the recovery period (Rp), and the mean temperatures of the wettest (BIO8) and coldest (BIO11) quarters showed negative correlations with Rs and Rc and Rs, respectively, In general, bioclimatic variables indicating wet conditions (e.g., BIO12, BIO13, BIO15, BIO16, BIO18) were negatively correlated with Rs, but positively with Rp and the growth reduction (Gr) The precipitation of the driest month (BIO14) showed negative associations with Rp and Gr.

The resilience index was correlated with the climate of origin, specifically precipitation (Figure 5). The results of the path analysis showed that the model including direct and indirect, through leaf weight (LW), effects of BIO16 on Rs was more parsimonious (lower AIC values; Table 4). According to the mode, the direct effect of BIO16 on Rs was significant, while the indirect effect through LW was not significant.

Considering the eight oak species sampled in Iturrarán and also in the field, we found that drought legacies previously reported for field oak populations were negatively related with Rs, but this relationship was not significant (r = −0.652; p = 0.08; Figure 6). That is, the species that were more resilient to drought in the studied arboretum presented lower drought legacies in the field.

The resilience of temperate oak species, growing under similar environmental conditions, in response to the unprecedented 2022 drought was related to the climate of origin, specifically the precipitation of the wettest quarter (BIO16; Figure 5). This inference is of correlational nature given that both traits and growth resilience depend on climate conditions and our approach does not allow disentangling these influences. Oak species from regions with drier conditions showed higher resilience and vice versa. Neither wood nor leaf traits showed significant correlations with growth resilience. A major novelty of our study is the use of retrospective growth data which allows framing the temporal context and impact of the 2022 drought in a common garden experiment.

First, the cross-species correlation between resilience and BIO16 could be explained because traits associated with a conservative strategy and persistence under resource limitation increased with aridity, whereas traits associated with resource use increased in wet regions (Schwinning et al., 2022). An increase of tree size and leaf area with wetter origins and a decrease of WD and LMA have been already reported (Perez-Martinez et al., 2025). Second, a site-level mechanism of reduced recovery could be due to soil-water recharge. In Mediterranean climates, soils are mainly recharged by autumn–winter rainfall. After the 2022 drought, insufficient soil recharge may have limited growth recovery in the following spring differently impacting the study oak species. It is noteworthy that the climatic variables related to resilience did not correspond to climate variables enhancing growth in the study arboretum (Table 2), where wood production increased in response to cold winter and summer conditions, linked to reduced carbohydrate consumption and evapotranspiration rates, respectively, warm springs leading to an earlier and more rapid xylem resumption, and wet conditions in July–August (Table 1) enhancing latewood production (see García-González and Eckstein, 2003; Alla and Camarero, 2012).

The post-drought growth recovery could be mechanistically interpreted as an enhanced radial growth, which is linked to a higher HV, i.e., a greater production of conductive sapwood area relative to canopy leaf area, thus improving drought resilience (Rowland et al., 2023). This would explain why functional traits (e.g., LMA, hydraulic diameter) are related to radial growth rate (Supplementary Table S3). A dense wood is formed by xylem cells with thick walls and abundant fibers, which enhance mechanical support and resistance to embolism, thus providing drought tolerance (Chave et al., 2009; Hacke et al., 2001; Li et al., 2018) and reducing the drought-induced mortality risk (Anderegg et al., 2016). Moreover, the drought legacies were lower for more resilient oak species indicating that an improved recovery capacity lead to smaller carryover effects, which can be explained by a higher concentration of non-structural carbohydrates (NSC) used as source of osmolytes or as carbon sink for meristems (Anderegg et al., 2015). Future studies could relate resilience with the concentrations of NSCs in leaves, stem sapwood and roots (Wang and Wang, 2015).

Growth resilience measured in the arboretum was negatively related to drought legacies measured in the field, as expected, but the relationship was not significant (Figure 6). This lack of significance may be explained by several reasons. First, we sampled a few individuals of each species growing under similar environmental conditions, whereas field studies are based on many individuals and populations occupying a wide range of environmental conditions (e.g., Anderegg et al., 2015; Bose et al., 2024; Gazol et al., 2018). Second, our analyses only considered a recent hot drought, whereas most studies dealt with several drought episodes of diverse intensity or timing differently impacting tree growth (Kannenberg et al., 2019). Third, the inconsistency between arboretum and field resilience data may be due to phylogeny. However, this was partly considered because most study species belong to the Lobatae and Quercus sections of the subgenus Quercus whose ancestors are inferred to have evolved at higher latitudes under short growing seasons (Hipp et al., 2018; Jiang et al., 2019).