Tree Species Are Differently Impacted by Cumulative Drought Stress and Present Higher Growth Synchrony in Dry Places

The increase in frequency and intensity of droughts due to climate change might threaten forests under stress levels causing dieback and mortality episodes. Thus, deciphering how tree species from within a region respond to drought along environmental gradients should help us to understand forest vulnerability to climate change. To enlighten contrasting drought responses of dominant tree species, we reconstructed vegetation activity using Normalized Difference Vegetation Index (NDVI) and radial growth using tree-ring width series. We studied six tree species, three angiosperms (Fagus sylvatica, Quercus humilis, and Quercus ilex) and three gymnosperms (Pinus sylvestris, Pinus nigra, and Pinus halepensis), inhabiting a Mediterranean region in north-eastern Spain. We investigated if reduced growth resilience and increased growth synchrony after successive droughts (1986, 1989, 2005, and 2012): (i) were related to cumulative drought stress and (ii) preceded forest dieback in dry sites as compared to wet sites. In 2016, dieback affected Q. ilex and P. sylvestris stands in dry sites showing lower growth rates and NDVI. No dieback symptoms were observed in other species from dry (P. nigra, P. halepensis) or wet (F. sylvatica, Q. humilis, P. sylvestris) sites. Hot and dry summer conditions constrained growth and reduced NDVI. During 2005, a severe drought affected all species, but growth drops were more marked in dry places. All species were able to recover after extreme droughts, albeit angiosperms displayed lower than expected values of growth after the 2012 drought. Growth synchrony was higher in dry sites than in wet sites, and the differences were higher after the 2005 drought. This study reveals that the sensitivity of tree species to drought in species inhabiting the same region is species dependent, and it is contingent on local conditions with higher effects in dry sites than in wet sites. We describe how a cumulative impact of successive droughts increases growth synchrony and triggers the occurrence of dieback events in Mediterranean forests.


INTRODUCTION
The occurrence of extreme climate events such as droughts and heat waves among others has been identified as a key factor reducing tree growth and leading to forest dieback (Allen et al., 2015;Senf et al., 2018;Resende et al., 2020). In Mediterranean forests, tree growth is strongly limited by water availability (Pasho et al., 2012) and drought-adapted tree species can present a bimodal growth, maximizing resource use during spring, early summer, and continuing in autumn, when water availability is not limited (e.g., Camarero et al., 2010;Gutiérrez et al., 2011). However, if recurrent drought and heat waves occur, even dieback in drought-adapted Mediterranean tree species happen Sánchez-Salguero et al., 2017;Sánchez-Salguero and Camarero, 2020). Given the prediction of higher frequency and intensity of drought and heat waves (IPCC, 2018), detecting which tree species and populations are more vulnerable to these climatic constraints emerges as a fundamental question (Brodribb et al., 2020).
To understand the temporal and spatial occurrence of hotter forest drought-induced dieback and mortality, and which species will suffer more, multidisciplinary approaches are needed (Anderegg et al., 2013). Drought induced dieback causes bud, leaf and shoot death as well as growth reductions (Dobbertin, 2005). However, leaf shedding and canopy dieback often happen during and after drought occurrence , and methods allowing the retrospective quantification of droughts impact on long-term growth and vigor are required. In this sense, tree rings are natural archives on how trees responded to past droughts (Fritts, 1976;Ogle et al., 2000). Moreover, combining tree ring growth with estimates of primary productivity and vegetation activity (e.g., Correa-Díaz et al., 2019;Vicente-Serrano et al., 2019) enhances our ability to upscale growth patterns from individuals or stands to regional scales (Babst et al., 2018;Gazol et al., 2018;Kannenberg et al., 2019). This information, combined with high-resolution climate data, improves our understanding on how forests respond to drought at different spatial scales and allow us to understand drought-driven forest mortality (Hartmann et al., 2018).
The response of tree species to drought can vary at the intraspecific level depending on both, regional and local factors (Orwig and Abrams, 1997;Gazol et al., 2017Gazol et al., , 2018Aragão et al., 2019;Baker et al., 2019;Crouchet et al., 2019). It is expected that tree species will be more vulnerable to drought near the rearedge (equatorward) and xeric limits of their distribution ranges (Sánchez-Salguero et al., 2017) where stronger legacy effects of drought on growth can be expected (Anderegg et al., 2015;Gazol et al., 2020). However, some studies have also reported greater sensitivity of tree growth to drought at the core and even the cold distribution margins (Muffler et al., 2020). In any case, the cumulative impact of successive droughts on forests can result in more pervasive drought legacies (i.e., lagged effects that can last for years despite drought being meteorologically alleviated), more marked growth reductions (Peltier and Ogle, 2019), and loss of resilience capacity (i.e., a reduced ability to recover growth values as those before drought occurrence). Interestingly, local factors as soil water holding capacity (Baker et al., 2019) and neighbor's competition for water (Crouchet et al., 2019) can also determine tree responses to drought (Galiano et al., 2012;Anderegg et al., 2013). However, it remains unclear the way in which all these local factors interact and modulate drought-induced dieback (Hartmann et al., 2018).
At the interspecific level, local environmental constrains interact with functional differences in plants response to drought among co-occurring species, such as differences in belowground and aboveground hydraulic conductance (McDowell et al., 2008;Johnson et al., 2018). Isohydric trees will be more vulnerable to carbon starvation because they regulate stomata closure to avoid xylem cavitation, whereas anisohydric trees are more vulnerable to hydraulic failure because they keep open stomata during drought (McDowell et al., 2008). Long-lasting growth reductions preceding tree death after drought have been found in several species, usually gymnosperms, and agree with the deterioration of their hydraulic system (Pellizzari et al., 2016). In contrast, carbon starvation has been considered as the main driver in declining angiosperms (Galiano et al., 2012;Klockow et al., 2018). These results question the validity of the isohydry-anisohydry dichotomy as a framework to determine dieback susceptibility because other local (e.g., Sánchez-Salguero and Camarero, 2020) and regional factors can limit growth (Körner, 2015). Recent studies suggest that factors such as oxidative stress, fire-drought interactions and root mortality among others can play important roles on tree mortality caused by drought (Parra and Moreno, 2017;Ameye et al., 2018). In this respect, drought timing arises as an important factor in Mediterranean regions since species inhabiting dry areas are more constrained by water availability and show responses to longer droughts, whereas those species from wet regions depend more on temperature and respond to short droughts (Pasho et al., 2011(Pasho et al., , 2012. Severe winter droughts or summer dry spells impact growth of coexisting tree species in different ways . In 2016 a severe forest dieback episode characterized by abundant leaf shedding (defoliation) and rise in mortality rates was observed in lowlands of the Montsec massif situated in Catalonia, north eastern Spain (Supplementary Figure S1). The dieback mainly affected Holm oak (Quercus ilex L.) and Scots Pine (Pinus sylvestris L.) (DeBosCat, 2017;Michel and Seidling(eds), 2017). Furthermore, trees in these regions were still recovering from the damages caused by the severe 2012 drought (DeBosCat, 2017). Motivated by this event, we started a dendrochronological study to evaluate the growth response to drought of six tree species (three gymnosperms of the Pinaceae family, and three angiosperms of the Fagaceae family) present along an altitudinal gradient (from 400 to 1300 m a.s.l.) in the Montsec massif. Furthermore, we combined this information with estimations of past vegetation activity at the ecosystem level derived from remotely sensed images (Normalized Difference Vegetation Index, NDVI; Myneni et al., 1995). Thus, we can obtain a better understanding on the relationship between radial growth and vegetation activity, and upscale the consequences of forest dieback from population to ecosystem levels. Specifically, we expect that: (i) all species are negatively impacted by drought resulting in severe growth reductions; (ii) growth and NDVI will be positively linked since wood formation is an important component of vegetation activity; and (iii) the cumulative impacts of successive droughts on growth and their legacies will be higher in dry low-elevation sites resulting in a stronger growth synchrony between species and stands, a lower vegetation activity and a higher risk of drought-induced dieback. We argue that the retrospective quantification of growth responses to past droughts can help us to understand why some tree species and stands are more sensitive to drought than others guiding future research on the mechanisms that lead to a greater drought vulnerability.

MATERIALS AND METHODS
The study area is the Montsec massif, situated in Catalonia, north eastern Spain (42.03 • N; 1.02 • E). The massif is 40 km long and 5 km wide with an elevation range from 400 to 1676 m a.s.l. With a Mediterranean climate it presents cold winters and hot summers characterized by a marked dry season from July to September. During the last years, the air temperature has increased in the region leading to more frequent and intense drought events since the 1980s (Supplementary Figure S2).
Due to the strong elevation and bioclimatic gradient there is a marked change in forest species composition with Aleppo pine (Pinus halepensis Mill.) and Holm oak dominating the lowland drier parts of the study area, and Scots pine and European beech (Fagus sylvatica L.), dominating the high-elevation sites. Other sub-Mediterranean species such as downy oak (Quercus humilis Mill.) and black pine (Pinus nigra J. F. Arn.) are present at intermediate elevations. Soils are rocky and basic.
We selected two different regions: a mid-to high-elevation site without evident signs of dieback (hereafter wet site) and a lowland site with marked signs of dieback (hereafter dry site; Supplementary Figure S1). In each site several tree species were selected (see Table 1). Beech and oak were sampled in the wet site, Aleppo pine (a natural forests and a plantation subjected to natural climate conditions) and Holm oak were sampled in the dry site, and Scots pine and black pine were sampled in both sites. For those sites in which the two species showing dieback were dominant (Q. ilex and P. sylvestris), 50-m long transects were established to quantify stand characteristics and dieback incidence ( Table 2).
For each species in each site, 10-25 dominant trees, usually separated by at least 5 m, were sampled for dendrochronological analyses. The diameter at breast height of each tree (Dbh) was measured at 1.3 m using a metric tape. In addition, we annotated whether the tree was dead, i.e., if it showed a crown without living leaves and branches, or alive. In the case of sprouting species as oaks we tagged them and checked that dead trees did not sprout 2 years after sampling. In the case of the two most negatively impacted stands (Scots pine and Holm oak populations from the dry site) we sampled eight (pine) and 10 (oak) individuals recently dead, respectively, by taking cross sections.

Tree-Ring Processing
Wood cores were extracted at 1.3 m using 5-mm increment borers (Haglöf, Sweden), and perpendicular to the slope. After that, wood samples were air-dried, glued and polished using a series of sand-paper grids until tree-ring boundaries were clearly visible. Cores were visually cross-dated and measured to the nearest 0.01 mm using a LINTAB measuring device (Rinntech, Heidelberg, Germany). Cross-dating accuracy was checked by using the software COFECHA (Holmes, 1983).
Individual tree-ring width series were detrended using a cubic regression spline with a frequency response of 0.5 at a wavelength of 32 years to obtain yearly values of detrended ring-width indices (RWI). The resulting standardized RWI tree-level series (RWIstd) were pre-whitened (fitting an autoregressive model to the time series) to obtain series of pre-whitened ring-width indices (RWIres). Finally, the RWIstd and RWIres individual series were averaged yearby-year using a biweight robust mean to obtain a standard and residual site chronology, respectively. All the processes of tree-ring series detrending and chronology computation were performed using the ARSTAN software (Cook and Krusic, 2007). Finally, we calculated three statistics to characterize the chronologies: the first-order autocorrelation of RWIstd (AR1), the mean sensitivity (MS) of RWIres to quantify the year-to-year relative changes in growth, and the mean inter-series correlation with the master chronology (Briffa and Jones, 1990).
We calculated the drought effect on tree growth following Kannenberg et al. (2019): (1) Since the average RWIstd (RWI std ) is equal to one this is a measure of the percentage of growth reduction due to the occurrence of drought.
Legacies of drought were calculated as the difference between standardized growth (RWIstd) and detrended growth variation (RWIres) predicted by drought following Anderegg et al. (2015). The RWIres predicted by the drought index (SPEI) represents the expected impact of drought on growth in the absence of autocorrelation. Negative legacies indicate that the growth is lower than expected by the effect of drought.

Climate and NDVI Data
We used 1 km 2 -gridded climatic data assembled at a weekly timescale based on values of daily mean air temperature, sunshine duration, wind speed, relative humidity, and precipitation coming from the Spanish Meteorological Agency (AEMET). The dataset was homogenized and quality-checked prior to the estimation of gridded values for the period 1962-2016 (see Vicente-Serrano et al., 2017). Based on this gridded dataset, we calculated the monthly water balance as the difference between precipitation (P) and potential evapotranspiration (PET). Further, the water balance (P-PET) was used to calculate the Standardized Precipitation-Evapotranspiration index (SPEI), which is a drought index that allows for spatial and temporal comparability of drought severity (Vicente-Serrano et al., 2010). The SPEI is a multiscalar drought index that may take negative and positive values, indicating The geographical location of the population sampled (site number, latitude, longitude, and elevation) is shown together with the average precipitation for the site according to the data presented in Supplementary Figure S2, number of sampled trees, its diameter at breast height (Dbh), estimated stand age (in years), and mean tree-ring width. In addition, the mean inter-series correlation with the site chronology (C), the mean sensitivity (MS) of ring-width indices, and the first-order autocorrelation (AC1) of ring-width data are shown. *Plantation. All woody plants (height > 25 cm) were identified, measured (Dbh) and their defoliation (%) visually estimated. For each site, the dominant species, the total density (D), the density of the dominant species (D dom.), the diameter (Dbh) and defoliation of the dominant species are shown. Dbh and defoliation values are means ± SD. These sites are the same presented in Table 1. dry and wet periods, respectively (Vicente-Serrano et al., 2010). Previous studies have found that the radial growth of Mediterranean tree species responds to drought at long timescales (Pasho et al., 2011), thus we calculated the drought index at 12-month long time-scales for the entire study period . We generated five NDVI time series, one per study site (see Table 1), by combining the MODIS Vegetation Index Products MOD13Q1 and MYD13Q1 (data available at https: //modis.gsfc.nasa.gov/data/dataprod/mod13.php). Since NDVI time series may contain spurious data such as missing values or cloudy pixels, data were filtered using the quality metadata provided by MODIS products. In this way, time series for the period 2003-2019 at a resolution of 16-day were generated by combining the 9 × 9 NDVI images obtained at high-resolution (250 m × 250 m) for each study site. Missing values were filled as the average between the previous and the posterior value. Yearly values of NDVI were calculated by summing the values from February to November of the year in which the tree-ring is formed, as the relationship between RWI and NDVI reaches maximum values during this period in Mediterranean forests (Vicente-Serrano et al., 2019).

Statistical Analyses
Bootstrapped correlations were used to test for the relationship between tree-ring width chronologies (RWIres) and mean maximum and minimum air temperature, precipitation and the 12-month long SPEI. Analyses were performed from September in the year before the tree-ring was formed to October of the year of tree-ring formation based on previous studies in these species (Andreu et al., 2007;Pasho et al., 2011). Significance was assessed by resampling 1000 times using the exact bootstrapping method in the Seascorr package (Meko et al., 2011). Mean maximum and minimum air temperature and precipitation data were detrended using a cubic regression spline as done for RWI. In this study all statistical analyses were performed in the R environment (R Development Core Team, 2019).
To assess the impact of severe drought on growth, we selected four drought events (1986, 1989, 2005, and 2012) which present SPEI values below −1.5 (Supplementary Figure S2c) and that have impacted the growth of trees in the Iberian Peninsula (Vicente-Serrano et al., 2017;Gazol et al., 2018). Further, we quantified the cumulative impact of successive drought episodes (1986, 1989, 2005, and 2012) on the growth of each species by summing the drought effects which were individually computed for each drought episode using Eq. 1.
Legacies of drought were calculated as the difference between observed RWIstd and the RWIres predicted by the 12-month long August SPEI for all the years in the period 1983-2016 to reflect those years in which growth was lower than the value predicted by the drought effect.
Linear mixed-effect models (LME; Pinheiro and Bates, 2000) were used to evaluate whether the relationship between tree-ring width chronologies (RWIres) and drought varies between dry and wet sites. This was done by testing the significance of interactions between the 12-month long August SPEI and site. Chronology identity was used as random factor (random intercept only) since each chronology represents repeated measures over the same tree population. We run all potential models combining site and SPEI and selected the most parsimonious one based on information theory (Burnham and Anderson, 2002). First, we listed the models according to their corrected Akaike Information Criterion (AIC) values. Second, we selected the model with the lowest degrees of freedom among those models with a difference lower than 2 AIC units from the lowest AIC model. The final selected model was evaluated graphically, and its fit was quantified with the marginal and conditional pseudo-R 2 values which account for the variation explained by the fixed only and fixed and random factors, respectively (Nakagawa et al., 2017).
We quantified synchrony across ring-width chronologies for the common period 1983-2016 (period in which all chronologies were robust and well replicated). Particularly, we tested for the existence of within-and between-site synchrony comparing dry vs. wet sites, as well as angiosperms vs. gymnosperms. Separate analyses were performed. Briefly, we fitted seven variance covariance (VCOV) models relating RWIres at site level and RWIres at functional group level. Similar analyses were performed with RWIstd and we decided to use RWIres which produced stronger results. The existence of within-and betweengroup synchrony in homocedastic or heterocedastic models was tested based on information criterion (Burnham and Anderson, 2002). Particularly, we compared the AIC for each model and we selected the most parsimonious one (i.e., the model with the lowest degrees of freedom among those models with a difference lower than 2 AIC units from the lowest-AIC model). After that, we also tested for the existence of temporal trends in synchrony by analyzing models of 20-year moving windows with 1-year lag.
Generalized Additive Mixed Models (GAMM; Wood, 2017) were used to study the temporal trend in NDVI and whether it varies between wet and dry sites. NDVI was modeled using a cyclic cubic regression spline to represent intra-annual variations and a thin plate regression spline to explain long-term NDVI trends. Further, we included interactions with site to decipher if the intra-annual variation and long-term trend varied between wet and dry sites. Site identity was included as random factor since NDVI series represent repeated measures for every location.
Finally, the degrees of freedom for each smooth term were set to a maximum of 6.
Linear mixed-effect models were also used to evaluate the relationship between tree-ring-width data and NDVI. We tested for the relationship between RWIstd and the NDVI (10-month long cumulative NDVI for November). Chronology identity was used as random factor as RWIstd series represent repeated measures over the same population. In addition, we tested for the influence of site (i.e., dry vs. wet site) and functional group (gymnosperms vs. angiosperms) in the relationship between RWIstd and NDVI for the period 2003-2016. Different models were created and ranked according to their AIC. Finally, we selected the most parsimonious model based on the AIC and the degrees of freedom of the model.
Bootstrapped correlations were calculated using the treeclim package (Zang and Biondi, 2015), Synchrony between ringwidth chronologies was calculated using the DendroSync package (Alday et al., 2017(Alday et al., , 2018, and LMEs were fitted and compared, in that order, with the nlme (Pinheiro et al., 2018) and MuMIn (Barton, 2018) packages. The mgcv package (Wood, 2017) was used to fit the GAMMs and the visreg package was used to visualize the regression models (Breheny and Burchett, 2017).

RESULTS
Patterns in growth variability (RWIstd) varied between sites but growth reductions in dry years were common and marked for all species (Figure 1 and Supplementary Figure S3). All species, excepting Aleppo pine and black pine from dry sites, showed significant positive correlations with June or July precipitation and negative correlations with the mean maximum temperature in June (Supplementary Figure S4). Significant positive correlations were also found between growth of most species and the precipitation in the previous winter (Supplementary Figure S4). Tree growth was strongly limited by drought as demonstrated by the strong relationship between growth and the 12-month SPEI (Supplementary Figure S4).
Tree growth reductions because of drought were more manifest in dry than in wet sites (Figure 1 and Supplementary Figure S5). The strongest growth reductions were observed in response to the 2005 drought at the dry sites (Figure 1 and Supplementary Figure S5). The cumulative drought-induced growth reductions were higher for the species in dry sites, and particularly for Scots pine and black pine (Supplementary Figure S5b). The results of the LMEs confirmed that the relationship of tree growth with the 12-month August SPEI was stronger in dry than in wet sites (Supplementary Tables S1, S2 and Figure 2).
The difference between the observed growth and that predicted by drought (i.e., expected drought impact on growth) showed that the 12-month August SPEI was not able to fully predict the impact of the 2005 drought on growth ( Figure 3A). Lower than expected growth values were found after the 2005 ( Figure 3B) and 2012 droughts ( Figure 3C). These drought legacies were more intense and lasted longer after the 2012 drought, particularly for beech and oaks. Most of the dead trees  Figure S6).
Within-site tree growth synchrony varied between wet and dry sites (Figure 4 and Supplementary Table S3). In the case of the comparison of dry vs. wet sites, we selected the homocedastic unstructured model, i.e., allowing for heterogeneous variances and covariances. Similar results were observed for RWIstd and RWIres (Supplementary Table S3). In the case of RWIres, the selected model (homoscedastic unstructured model) provided support for the existence of lower synchrony in wet sites (0.44 ± 0.06) than in dry sites (0.61 ± 0.06). There was also a notable within-group variation (0.48 ± 0.06). In the case of the comparison between gymnosperms and angiosperms, the best VCOV model proved to be a broad evaluation model suggesting the lack of synchrony. The temporal trend in synchrony showed an increase in dry sites as compared to wet sites (Figure 4).
We found that NDVI peaked during May and winter (December-January) and decreased during summer (July and August) but the pattern was different in dry and wet sites (Figure 5 and Supplementary Figures S7-S9). The GAMM accounted for 50% of the variation in NDVI and showed that the intra-annual variation varied more than the temporal trend (Supplementary Figure S8) and that this variation was stronger in dry sites (Supplementary Table S4). The longterm trend of NDVI showed reductions in the dry years  Figures S7, S8), which were more marked in dry sites (Supplementary Figure S9). Lower NDVI values were found in dry places suggesting low productivity (Supplementary Figure S9). In this sense, we found that RWIstd and NDVI were positively related (0.37 ± 0.03; p < 0.05). The LME also showed that RWIstd depended on site conditions (0.86 ± 0.15; p < 0.05; see Supplementary Table S5) but not its interaction with NDVI, and so the slope of the relationship between NDVI and RWIstd was the same in dry and wet sites but the intercept varied ( Figure 5). The model accounted for 79% of the variation in RWIstd and 55% of the variation was explained by the effects of NDVI and site.

DISCUSSION
We found that the radial growth of the six studied tree species was greatly constrained by early summer drought according to our first hypothesis. This phenomenon was fully evidenced during the 2005 extreme drought that lead to growth reductions higher than 40% in all species. We also found that growth reductions were more marked in dry places due to a stronger dependency on precipitation and leading to an increasing trend in withinsite spatial synchrony between species, a phenomenon that may be considered an early warning of the impacts of climate change (Shestakova et al., 2016;Tejedor et al., 2020). Drought legacies were particularly intense and durable for Fagaceae species after the 2012 drought and suggest that the long-term deterioration of the hydraulic system and/or the depletion of carbon reserves can explain drought-induced mortality in Holm oak in the studied sites (Galiano et al., 2012). However, we cannot discard than the causes that lead to tree mortality can depend on other factors  , 1986, 1989, 2005, and 2012).
such as oxidative stress (Ameye et al., 2018). The strongest legacy, i.e., the deviation between growth observed in 2005 and droughtpredicted growth, pointed again to the largest impact of that drought compared to previous dry events. Our results point out to cumulative effects of more frequent hotter droughts on tree growth and vigor, dieback, and increased mortality of vulnerable species in dry sites from the studied region. Further, droughtinduced growth reductions also resulted in vegetation activity reductions because greater growth sensitivity to drought was strongly linked with a remarkably lower vegetation activity as indicated by the NDVI drops. Hot and dry conditions leading to increased evapotranspiration rates and drought stress in spring and summer (Supplementary Figure S10) can have pervasive and negative consequences on the growth recovery of tree species in sites with low soil water-holding capacity (Sánchez-Salguero and Camarero, 2020), as the one studied here. Whether this result can be extrapolated to other regions remains unclear due to the mainly local extent of our study. However, we demonstrate that FIGURE 4 | Temporal trends in within-group tree-ring width (RWIres) synchrony calculated for dry (red symbols and areas) and wet (blue symbols and areas) sites. The synchrony was calculated for the best model considering 20-year moving intervals lagged by 1 year over the period . The x-axis shows the central year of the moving intervals. Lines and points represent observed values and shadows are standard errors. Synchrony increased sharply between the second and the third periods (1985-2005 and 1986-2006), and was higher in dry sites than in wet sites during the following periods. The results correspond to a homocedastic variance covariance (VCOV) unstructured model. accounting for growth sensitivity to drought helps to identify those species and sites most impacted by drought.
We found that drought is an important factor constraining the growth of the studied species as previously found in comparable Mediterranean ecosystems (Pasho et al., 2012;Gazol et al., 2018). All species, excepting Aleppo pine and black pine, responded positively to precipitation and negatively to maximum temperature during early summer (Supplementary Figure S4). Aleppo pine is a species adapted to drought that may present a bimodal growth under xeric conditions (e.g., Camarero et al., 2010) and its growth strongly relies on precipitation during winter (De Luis et al., 2013). It is also plausible to think that under dry circumstances, the growth of black pine also depends on water availability before the growing season which recharges soil water pools in winter (Sangüesa- Barreda et al., 2019). Indeed, these two pine populations responded significantly to the 12month long SPEI indicating that their growth depends on winter conditions. In this sense, the impact of SPEI was only nonsignificant in the case of Scots pine and black pine in wet sites. Pasho et al. (2011) found that gymnosperms from wet sites responded to drought at short-time scales in contraposition to other species from dry sites. Furthermore, the black pine population in the wet site was particularly young and thus it can be less vulnerable to drought occurrence (Martínez del Castillo et al., 2018). These results also indicate that local conditions can influence how tree growth responds to drought, and the FIGURE 5 | Relationship between the standard tree-ring width chronologies (RWIstd) and the 10-month long cumulative NDVI for November in dry (brown circles) and wet (blue squares) sites. The lines show the effect of NDVI on RWIstd according to the linear-mixed effect models (see Supplementary  Table S5) and the dashed lines show the 95% confidence intervals for the predictions. Red and orange symbols show the values observed in the 2005 and 2016 droughts, respectively.
increased sensitivity to drought of the two pine species in the driest sites suggests a higher drought-vulnerability toward drier sites (Sánchez-Salguero et al., 2017). This conclusion is supported by the stronger impact of SPEI on growth in dry sites (Figure 2), and it is further corroborated by the enhanced growth synchrony in those sites (Figure 4).
The extraordinary growth reductions observed during the 2005 drought (Figure 1 and Supplementary Figure S5) were not captured by the effect of SPEI on RWI (Figure 4), in contrast to what occurred in response to previous drought events (Figure 4). That is, the observed growth in 2005 was much lower than the growth predicted by drought (Figure 4). Furthermore, the marked increase in growth synchrony between species in dry sites started when the 2005 drought was considered in the time window (i.e., from the 1985 to 2005 period onward). These results point out to the extraordinary dimension of that drought regarding its impact on forests (González-Hidalgo et al., 2018). In this sense, widespread foliage loss, dieback and mortality of Mediterranean tree species have occurred since the early 1990s and aggravated in response to the strong 2005 drought in the region (Carnicer et al., 2011). This drought was characterized by elevated temperatures and evapotranspiration rates during the growing season (Supplementary Figure S10). Probably, the negative effects of the 2005 drought on growth lasted over the last decade, and its negative impacts have reactivated by recent warming trends and ended in a strong dieback event in 2016 (DeBosCat, 2017). This is further reinforced by the aggravation in drought legacies, i.e., the difference between observed and drought-predicted growth, from the 2005 to 2012 droughts (Figure 3). These results confirm that the impacts of drought on growth can last for several years (e.g., Ogle et al., 2000;Anderegg et al., 2015;Peltier et al., 2016), and that these drought legacies can aggravate with the accumulation of severe droughts resulting in more pervasive growth reductions (Peltier and Ogle, 2019). This is particularly notable in the case of Holm oak, the most affected species, which showed more pervasive drought legacies after the 2012 drought, which was characterized by warm-dry conditions during the 2011-2012 winter and the 2012 spring, than after the 2005 drought. However, deviation between observed and drought-predicted growth were in general lower in response to the 2005 drought than after drought occurrence. Deciphering whether those growth deviations are a consequence of prolonged drought impacts may require further research over longer time periods. In any case, the results points that the cumulative impact of hotter droughts can have negative consequences on some Mediterranean tree species and forests.
There might be a trade-off between growth resistance (i.e., how growth is impacted by drought) and growth recovery (i.e., capacity to recover after drought) that might depend on species adaptations. Gymnosperms are known to have marked growth reductions (i.e., low resistance) in response to drought but in some cases, they can recover fast. Such is the case of Mediterranean pines (Gazol et al., 2020). Thus, growth resistance and recovery alone provide useful but incomplete information as many other factors can ultimately determine how tree species tolerate drought. For instance, the occurrence of favorable conditions at some point may facilitate the performance of a species that can suffer strong mortality rates when climate conditions worsen (Jump et al., 2017). In this respect, combining different sources of information might be useful and the use of permanent plots that allow quantifying stand characteristics and recruitment capacity are of great value. This evidences the utility of the intensive, long-term monitoring of those sites showing dieback to understand its causes and consequences and how it varies between species.
Angiosperms were more severely defoliated in 2016 than gymnosperms in the region (DeBosCat, 2017). Interestingly, GAMM showed a stronger impact of the 2005 drought on NDVI in wet than in dry sites (Supplementary Figure S8) which might indicate that primary production was more constrained in those sites dominated by Q. humilis, F. sylvatica, and P. sylvestris. However, growth reductions in response to drought were much more intense in gymnosperms, particularly those in dry sites during the 2005 drought (Figure 1 and Supplementary Figure S5), than in angiosperms. This is in line with the higher prevalence of abrupt growth reductions in gymnosperms (Schweingruber, 1986), and with previous studies indicating that pines from the Mediterranean region present a plastic growth response to drought (Camarero et al., 2010). According to our results the capacity to recover from drought of gymnosperms from dry sites is faster than that of angiosperms in general (Gazol et al., 2018). However, it is also true that most studies have found longer and more lasting legacies of drought in gymnosperms than in angiosperms (e.g., Anderegg et al., 2015;Gazol et al., 2020). Recently, DeSoto et al. (2020) found that those trees that die because of drought presented lower growth resilience than surviving ones. Our results indicate that the growth of Holm oak was lower than expected in the 4 years after the 2012 drought occurrence. However, these long-lasting growth reductions were also common in the rest of tree species, particularly in Fagaceae species from wet sites, and did not result in dieback and elevated mortality.
The potential occurrence of pervasive drought legacies cannot explain why some Holm oak and Scots pine trees showed more dieback and mortality rate despite the mechanisms can be partially inferred (e.g., Cailleret et al., 2017). It is unlikely to think that differences in morphological, structural, or hydraulic traits drive differences in drought-induced mortality in the Montsec sites. Most likely, a feasible explanation for the decline of several holm oak and Scots pine trees can be on the influence of local site conditions in general (Kannenberg et al., 2019;Stralberg et al., 2020), and belowground conditions in particular. In these drought-prone and low productive sites, tree survival during drought may depend on the capacity to access to rock bedrock groundwater (McDowell et al., 2019;Mackay et al., 2020). Several studies have found that small trees from different oak species can be more vulnerable to drought than larger individuals (Colangelo et al., 2018), and this could be due to a smaller root system and a reduced ability to use deep soil water sources (Ripullone et al., 2020). Along this, the diameter of the trees was the smallest in the two sites impacted by drought (see Table 1), which can suggest a greater sensitivity of smaller individuals. However, we lack enough data to provide a statistically sound comparison between dead and living trees. Future studies will benefit from a random sampling of trees (Nehrbass-Ahles et al., 2014) together with a characterization of tree height and stand density of all dead and living individuals within permanent plots, which can improve the interpretations of the effects of tree size on growth response to drought. If the causes of dieback and tree death depend on the ability to access belowground water resources, the size and architecture of trees can also determine their vulnerability by regulation of non-structural carbohydrates allocation (Galiano et al., 2012). On the other hand, the dieback and mortality of Scots pine could be more related to hydraulic failure (Pellizzari et al., 2016). However, more research is needed to fully understand whether long-term hydraulic failure or carbon depletion contributed to explain Holm oak dieback and mortality, which seems to be mainly triggered by local soil factors.
How drought-induced growth reductions translate to changes in primary productivity and vegetation activity at the ecosystem level remains unclear (Gazol et al., 2020). Interestingly, we found a positive relationship between RWIstd and NDVI, which is in accordance with previous results (Vicente-Serrano et al., 2019). We also found that NDVI was lower in dry sites indicating that those sites that are more susceptible of drought-induced growth reductions are also less productive. This confirms that carbon uptake and growth are limited by water availability (Körner, 2015). Moreover, seasonal patterns in NDVI were very different between sites indicating that phenology varies between sites according to the species dominating the forests (Aragonés et al., 2019). However, recent studies have identified that NDVI variations depend on the interaction of multiple factors cautioning about the interpretation of the impact of single factors (Piedallu et al., 2019). These results indicate that combining tree-and ecosystem level measures help to understand how dieback events translate into changes in vegetation activity patterns at the ecosystem level.

CONCLUSION
We combined tree-ring measures and satellite-derived vegetation indexes to show that droughts severely constrain growth and that this has an impact on forest productivity. As hotter droughts become more frequent and intense their negative consequences on terrestrial ecosystems are exacerbated. However, the impacts of drought on growth are clearly site-and species-specific thus making difficult to drive universal rules on how trees respond to drought. Accordingly, additional local studies combining different approaches are required to improve our understanding on how forests respond to drought and what can be expected in the near future. Such is the case of this study, which despite being local in geographic extent offers information on how six tree species varying in their functional characteristics respond to drought in a region.
The most important conclusion of our study is that dry sites are more vulnerable to drought but that the species will respond in different ways. We found that growth reductions in response to drought were strong in Mediterranean pines, but mortality was stronger in Holm oak and in Scots pine inhabiting dry sites. Thus, understanding how tree growth responds to drought is not enough to understand the resilience capacity of forests to drought. Further research is required to corroborate if the drought-induced dieback and growth decline of some species is determined by their ability to access groundwater resources. The importance of such soil water and nutrients in driving dieback and tree mortality episodes in dry ecosystems, and how they interact with the physiological behavior of tree species, should be considered in further studies.

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

AUTHOR CONTRIBUTIONS
AG led the analyses and writing of the manuscript with substantial contributions by JC. GS-B and JC were responsible of tree-ring analyses. All authors carried out field sampling, edited and discussed the final draft. FUNDING JC, AG, and RS-S acknowledge the support of the FORMAL (RTI2018-096884-B-C31) and LESENS (RTI2018-096884-B-C33) projects from the Spanish Ministry of Science, Innovation and Universities; and VULBOS project (UPO-1263216) from the ERDF funds Andalusia regional government 2014-2020. GS-B was supported by Spanish Ministry of Economy, Industry and Competitiveness Postdoctoral grant (FJCI 2016-30121; FEDER funds) and from the European Union through the CANOPEE Interreg POCTEFA project.