Growth plasticity of conifers did not avoid declining resilience to soil and atmospheric droughts during the 20th century

1.   Introduction

Understanding how forests respond to environmental drivers is crucial for projecting impacts of climatic changes worldwide (IPCC, 2013). Drought stress is a major driver of plant distribution (Engelbrecht et al., 2007), demography (Anderegg et al., 2016) and physiology (McDowell, 2011). Drought severity is predicted to increase in many regions as a result of climate change (Ault, 2020; Cayan et al., 2010; Vicente-Serrano et al., 2020; Vicente-Serrano et al., 2020a). Even in the absence of changes in rainfall regime, increased temperatures will lead to increased atmospheric water demand and increased vapour pressure deficit (VPD). Other drivers, such as nitrogen deposition and increased atmospheric CO₂ may lead however to significant water savings for plant leaves and canopies (Levesque et al., 2017; Craine et al., 2018; Guerrieri et al., 2020). Because of the interacting nature of many environmental drivers and plant physiological processes, it is difficult to predict whether long-term chronic effects of repeated severe droughts become amplified during subsequent events (Reichstein et al., 2013). By contrast, it is possible that past drought events coupled to elevated CO₂ concentrations and nitrogen deposition, elicit positive adjustments to future severe droughts (Drake et al., 2018; Franklin et al., 2020).

Trees can be plastic in their responses to severe droughts, in the sense that the sensitivity of the growth response to an environmental driver (i.e., the slope of the relationship between growth and that factor) can vary over time (Peltier and Ogle, 2020), i.e., lack of stationarity (Wilmking et al., 2020). This plasticity may result in adjustments to the gradually changing environmental conditions, such that impacts of subsequent severe droughts on growth may either increase or decrease. Overall, the relative balance of these positive and negative effects caused by antecedent events and their legacies is central in understanding the long-term impacts of extreme climatic events, and especially drought, on tree responses (Schwalm et al., 2017; Yin and Bauerle, 2017; Serra-Maluquer et al., 2021).

Adjustments of trees to drought are thought to depend on the existence of a stress memory (Fleta-Soriano and Munné-Bosch, 2016; Julio Camarero et al., 2018; Mu et al., 2022). Adjustment mechanisms that are adaptive during periods of severe stress are costly to maintain during subsequent periods favourable to growth (DeWitt et al., 1998; Valladares et al., 2007), but they may reduce the negative impacts of the stressful events in the face of their increased frequency and/or intensity. The alternative hypothesis is that accumulated physiological damage caused by recurrent and severe droughts will lead to growth decline and, ultimately, to catastrophic loss of function and tree death, as stress severity increases.

Whether soil drought severity (understood in terms of reductions in the balance between rainfall and evaporation) has changed systematically across most biomes is debated (Trenberth et al., 2014; Ault, 2020; Padrón et al., 2020). However VPD, which mostly depends on temperature, has increased globally since about the 1970s–1980s (Yuan et al., 2019). Increased VPD amplified not only the frequency of years of high-VPD, but also of events when high-VPD and soil drought occurred jointly (De Boeck et al., 2010; Zscheischler and Seneviratne, 2017; Alizadeh et al., 2020). The frequency of combined hot drought events is predicted to increase further in the future (IPCC, 2013; Allen et al., 2015). Because xylem growth is mechanistically linked to plant water potential and turgor (De Schepper and Steppe, 2010; Cabon et al., 2020a), this may lead to reduced resistance (i.e., growth occurring during the severe drought) and resilience (post-drought growth relative to growth during non-drought years). Here we take advantage of long-term records of growth in tree rings during the 20th century to examine whether limitations to growth caused by two classes of potentially very different droughts, i.e., anomalies in site water balance between precipitation and potential evapotranspiration (i.e., soil drought) and anomalies in VPD (i.e., atmospheric drought), resulted in plastic, i.e., non-stationary, growth responses over time. We then examine the consequences of these plastic responses. We acknowledge, however, that our tree-level analysis does not explicitly account for compensating processes at the stand level, such as stand density reductions resulting from drought-induced mortality (Lloret et al., 2012; Giuggiola et al., 2016).

Because of land-atmosphere feedbacks, strong spatial and temporal interactions are expected to occur between soil and atmospheric droughts (Seneviratne et al., 2006; Zscheischler and Seneviratne, 2017; Pagán et al., 2019; Zhou et al., 2019). This correlation is especially strong at multi-daily to seasonal time scales, as a result of weather patterns. Here, we exploit the low correlation existing between cumulative 12-month SPEI (Standardized Precipitation-Evapotranspiration Index), a widely used drought index accounting for changes in site water balance, and spring-to-summer VPD to explore individual and combined effects of these two types of severe droughts during the 20th century. Our approach to exploring the effects of severe droughts based on SPEI and VPD is explained in Fig. 1. Because both SPEI and VPD are standardized on a Gaussian scale (i.e., site-specific anomalies), we employ identical thresholds to select events of equivalent severity, following prior work, e.g., Anderegg et al. (2020), Li et al. (2020), and Senf et al. (2020). Combining 1-year SPEI and spring-to-summer VPD to determine growth responses to severe climatic events is important for several reasons: 1) temporal variability in SPEI and other drought indices is driven more by precipitation than by atmospheric evaporative demand for many biomes (Vicente-Serrano et al., 2020b), 2) atmospheric evaporative demand and soil drought conditions have very different effects on plant physiology (Novick et al., 2016; Kannenberg et al., 2019; Brodribb et al., 2020; Grossiord et al., 2020) and both can contribute to growth reductions (Williams et al., 2010; Anderegg et al., 2013) and 3) tree growth shows strong memory effects (i.e., antecedent conditions can affect drought responses and drought effects can last for several years post-drought) at different time scales (Peltier et al., 2018). Although, VPD and SPEI provide only a partial representation of atmospheric and soil drought conditions, their joint examination provides complementary perspectives on the responses of growth to severe drought events (Liu et al., 2020).

Figure  1.  Bivariate distribution of anomalies in 1-year cumulative SPEI and mean spring-to-summer VPD for 1,660 sites. Both SPEI and VPD are standardized (centred around zero and normalized to standard deviation of 1) for each site independently. As severe drought years, we select those site-years when either a) spring-to-summer VPD ≥1.5 and the sum (VPD - SPEI) < 3.0 (pale blue points), b) 1-year SPEI ≤ −1.5 and the sum (VPD - SPEI) < 3.0 (green points), or c) the sum (VPD - SPEI) ≥ 3.0 (red points on or above the red bisector line of slope +1). Cases a) and b) represent situations when the severe event is defined based on the occurrence of either severe spring-to-summer atmospheric or 12-months soil water balance limitations but not both (according to an OR scenario of a Survival Kendall multivariate distribution, i.e., either one of the severe events occurs (Zscheischler and Fischer, 2020)), and case c) when the two types of events co-occurred (AND scenario (Zscheischler and Fischer, 2020)), with a variable contribution of the two components. Values exceeding a threshold of 1.5 on a standardized scale equate to a probability of occurrence of ~7% for univariate distributions. On average, ~11% of all site-year events are classed as severe (i.e., less than the expected 14% if the two univariate distributions did not overlap). Of these events, 20% of all SPEI and 27% of all VPD are compound events. Grey points in otherwise coloured areas indicate events in time series with <5 events overall, which are excluded from further analysis, resulting in 1,660 series retained (cf., Tables S1–2 for further details on species and sites sampled). (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

DownLoad: Full-Size Img PowerPoint

Following the reasoning above, here we test the following hypotheses: 1) The sensitivity of conifer tree growth to VPD and soil moisture balance during droughts changed during the 20th century for the three types of droughts (high-VPD, low-SPEI or compound droughts). 2) As a consequence of these shifts in sensitivity throughout the century, drought-year and post-drought growth differed from the growth obtained if sensitivities had not changed over time. Because high-VPD events became globally more frequent towards the end of the 20th century, we hypothesise that the increased severity of stressful events elicited a memory response and that the consequent endogenous adjustments reduced the negative effects on growth for VPD and compound events. 3) We finally assess whether this plasticity in growth responses was sufficient to prevent long-term trends of growth reductions during or after the severe droughts of the 20th century (complete adjustment) or whether there is evidence of incomplete adjustment, depending on the nature of the drought (high-VPD, low soil moisture or compound events).

We focus on Gymnosperm, and specifically coniferous, trees for two reasons. First, there is evidence that Gymnosperms, and more extensively the Pinaceae, have been severely affected by extreme events in recent decades (Williams et al., 2010; Gazol et al., 2018; Cailleret et al., 2019; Anderegg et al., 2020; DeSoto et al., 2020). For example, Gazol et al. (2018) showed that evergreen gymnosperms of semi-arid Mediterranean forests showed the lower resistance to drought episodes compared to deciduous angiosperms, while Cailleret et al. (2019) showed a gradual increase in inter-annual growth variability and a decrease in growth synchrony in conifers in the final part of the 20th century. Finally, DeSoto et al. (2020) showed that conifers show stronger and longer legacies of drought on tree radial growth relative to angiosperms. Second, significant historical trends were reported in both growth resistance and resilience of Gymnosperms for the second half of the 20th century (Li et al., 2020; Zheng et al., 2021). Although Angiosperms are also known to have been affected by severe droughts at several locations around the world, the availability of a large tree core collection for Gymnosperms provides the opportunity for a robust test of these hypotheses.

2.   Materials and methods

2.1   Tree-ring width data

We employed the site-averaged quality controlled tree ring width dataset produced by Zhao et al. (2019) and downloaded from https://doi.org/10.5061/dryad.kh0qh06. This dataset is a subset of the International Tree-Ring Data Bank (ITRDB), a global repository for tree-ring data, mainly from extra-tropical biomes (Grissino-Mayer and Fritts, 1997) and focusing on living trees. To obtain site-specific chronologies, individual ring width data were detrended and standardized using the dplR library (Bunn, 2008). Starting from the raw files, we used a cubic smoothing spline with a 50% frequency cut-off at 30 years to eliminate long-term trends in all time-series (2,195 series) (Klesse, 2021). Detrending and standardization involved dividing ring-width values by the corresponding fitted values to generate a series of normalized growth anomalies with mean close to 1 (Cook, 1985; Cook and Kairiukstis, 2013). Finally, series were averaged onto site chronologies using Tukey's Biweight Robust Means. We cut all datasets to the 20th century. A first level of filtering was carried out to retain only those sites in which a) gymnosperm species were present (where sites represent sampling locations within species, 1,830 sites). We then retained only those sites in which b) at least five severe SPEI or VPD events occurred during the 20th century, to reduce the uncertainty in the detection of temporal trends (1,707 sites). Additionally, c) we only retained sites where, for all decades in the 20th century, series age >90 years (1,686 sites) and finally, d) in which (see later) VPD estimates were extracted from data available from the local cell as opposed to interpolation from surrounding cells (1,660 sites, genera and spatial distributions are shown in Table S1 and Fig. S1). Severe events occurring at short time distances from one another were considered separate events. Separately, as a cross-check, we also considered the subset of 429 series that began before 1900 and ended after 2000. Tree-ring data for non-surviving trees were not used, as they do not allow assessing resilience. We acknowledge that this may bias our results towards higher resilience, as stand level responses (e.g., mortality) may result in lowered competition after the disturbance event.

2.2   Site-level climate, drought indices and additional variables

To characterise climatic differences across the 1,660 sites, elevation-corrected monthly temperature and precipitation data were downloaded from the CRU database v.4.05 on a 0.5° latitude by 0.5° longitude grid (https://crudata.uea.ac.uk/cru/data/hrg/). To estimate soil water balance at each site, we used the Standardized Precipitation Evapotranspiration Index (SPEI) (Vicente-Serrano et al., 2010). SPEI is based on the difference D between precipitation PRE and potential evaporation PET. D values are cumulated over pre-defined periods of time up to a specified month. The standardization of site-level SPEI allows comparing drought anomalies across sites. SPEI thresholds of drought severity have a probabilistic meaning, e.g., SPEI ​= ​−1.0, −1.5 and −2.0 imply a 16%, 7% and 2% probability or a one in six-, 14- and 50-year return intervals, respectively. SPEI has been validated to detect moisture anomalies for agricultural (Zipper et al., 2016; Páscoa et al., 2017) and environmental applications (Pasho et al., 2011; Greenwood et al., 2017; Zhang et al., 2017). We downloaded 6-month and 12-month-long SPEI data ending in the months from May to August for the Northern hemisphere (November to February for the Southern hemisphere) from the CSIC SPEI global database available at http://spei.csic.es/database.html. We compared different SPEI thresholds (i.e., < −1.0, < −1.5 and < −2.0, respectively) and used different end-months for both datasets. Similar results were obtained with the various combinations (data not shown). August was preferred since it cumulates D through the entire growing season and the previous 12 months. From the same website, we also downloaded time series of 3-months SPEI values ending in the months of February, May, August and November for each site and year. The season with the lowest SPEI value (variable Season; spring: March–May, summer: June–August, autumn: September–November, winter: December–February; reversed classification for the southern hemisphere) was determined for each site.

We also considered severe atmospheric droughts based on high values of vapour pressure deficit VPD spanning the period 1900-present. We downloaded monthly gridded temperature and vapour pressure data from the CRU time series (TS4.04 Dataset, Harris et al., 2020). The VPD time series has much lower spatial coverage for the second than for the first half of the 20th century, because of lower density of vapour pressure measurements. In the CRU dataset, VPD is estimated either from primary vapour pressure measurements or from empirical regressions against daily temperature ranges, a variable more widely measured than vapour pressure. In version TS4 of CRU, these synthetic estimates of VPD are calculated from primary station observations as opposed to gridded data and they therefore reflect the intrinsic spatial dependency of daily temperature range as opposed to the mean regional climatology (with a mean correlation decay distance of about 750 ​km). We used the available quality flag to exclude all site-year combinations where VPD interpolation was carried out from the regional climatology as opposed to the local gridded cell. We calculated absolute vapour pressure for the months of May to August, then obtained saturated vapour pressure using formulas given in (Yuan et al., 2019), corrected for atmospheric pressure and site elevation, then calculated VPD as saturated minus actual vapour pressure. We calculated temporal anomalies of VPD for each site and forced a normalization of the data onto a Gaussian curve with mean of zero and standard deviation of one, similar to SPEI. We chose a threshold of +1.5 to detect years of severe high-VPD values.

Because the same probability threshold was chosen for both types of droughts, we could examine differences in tree responses to the two classes of events. To detect years when strongly negative values of SPEI and strongly positive values of VPD co-occurred, we summed years when (VPD ‒ SPEI) ​≥ ​3.0. Under those conditions, either a) the two types of severe events co-occurred (VPD ≥1.5 and SPEI ≤ −1.5), or b) one of the two drivers did not reach the threshold for severe values but the second driver more than exceeded its threshold thereby compensating the other. This approach is an approximation to modelling compound severe events using a survival Kendall distribution (Zscheischler and Seneviratne, 2017), a frequently applied bivariate copula for the study of, e.g., the return periods of extreme temperature and precipitation and their effects on crop yield variability (Zscheischler et al., 2017). We similarly determined the sum of the 1-year SPEI and of the spring-to-summer VPD values that occurred during the four years prior to the severe event (SPEI_prev and VPD_prev) as well as during the four years following the event (SPEI_post and VPD_post). Overall, each chronology had between 5 and 23 severe events (median of 10 events). On average, 10.8% of all site_year events were classed as severe. Of these events, 20% of all SPEI events and 26.6% of all VPD events were compound events. Overall, we examined >16,000 events. Before use in the regressions, spline-frequency detrending was applied also to SPEI and VPD data, as above for tree ring data.

2.3   Growth time series modelling

The auto-correlated nature of tree growth must be taken into account in the analysis of growth rates during (drought-year growth) or following severe events (post-drought growth). A frequently applied procedure is to divide growth during or after a severe event by growth during a reference condition (e.g., growth prior to the severe event, Lloret et al., 2011), thereby obtaining estimates of resistance and of resilience ratios. This is problematic because in many widely employed approaches, resistance and resilience strongly and spuriously co-vary with one another because of growth year-to-year persistence, an issue that was only recently identified (Li et al., 2020; Zheng et al., 2021; Serra-Maluquer et al., 2022). In Supplementary Materials, we show that spurious covariation between resistance and resilience occurs in almost all of the pairs of resilience (R_s)/resistance (R_t) indices found in the literature. Consistent with this, we employ the approach by MacGillivray et al. (1995) for all our subsequent analyses, because this approach ensures that resilience and resistance are mathematically independent from one another, albeit they may still show temporal covariation as a result of autocorrelation in the time series, reflecting a real biological property as opposed to a mathematical artefact. Using this definition, resistance R_t is obtained dividing observed drought-year growth by the expected growth that would have taken place had the drought not occurred in that year. Conversely, resilience R_s is obtained dividing observed post-drought growth by the growth expected in the hypothesis that the drought had not occurred (Martínez-Vilalta et al., 2012). Expected drought-year growth and post-drought growth are derived by employing the autocorrelation properties of the time series and forecasting the time series for those years.

Growth expected to occur in the absence of the severe event and the subsequent four years was estimated for each time series independently using ARIMA (autoregressive moving average) modelling. The period of four years over which to assess growth is a compromise between the recovery usually taking up to 4 years after a severe drought and the need to minimise the overlap with subsequent drought periods (Anderegg et al., 2015; Wu et al., 2018). For each time series, we determined the four parameters of the best non-seasonal ARIMA model for the entire time series, i.e., intercept c, autoregressive term p, moving average term q and degree of first differencing d, using the ARIMA function of the fable package in R (Hyndman and Khandakar, 2008) applied to log-transformed tree ring indices (to avoid occasional negative predictions). The ARIMA function uses a variation of the Hyndman-Khandakar algorithm, which combines unit root tests, minimisation of Akaike Information Criterion AIC and maximum likelihood estimation to automate a stepwise search procedure for the best fitting ARIMA model (Hyndman and Athanasopoulos, 2018). We did not include environmental variables as predictors because environmental drivers of growth may have changed over time during the 20th century (Brienen et al., 2012; Duchesne et al., 2019), a finding that our own results support (cf., Table S9 for a summary of all the ARIMA models). In addition, the analysis in the Supplementary Materials demonstrates that modelling based only on ARIMA properties replicates the main features of the time series. After excluding the values corresponding to severe SPEI or VPD years (plus subsequent four years), the fitted ARIMA parameters were used to derive expected growth values for each event using a Kalman smoother, obtained using the average of the forward and backward forecasts for the missing section corresponding to the severe event. The interpolation routine used the interpolate function in fable. The median root mean square error of the forecasts across all models was 0.19 (min ​= ​0.04, max ​= ​0.76). A summary of the workflow employed for the analysis is given in Fig. 2. In the case of partially overlapping severe events, the interpolation was carried out for the entire sequence of excluded years. We also considered two more ARIMA models. One in which the autocorrelation properties were determined after excluding all drought episodes in each non-gapped filled time series and one in which ARIMA parameters were derived sequentially after excluding each drought at a time. All ARIMA methods gave very similar results. Finally, all log-transformed values were exponentiated.

Figure  2.  Summary of the workflow employed for the construction of the resilience/resistance indices. Further details are given in the text.

DownLoad: Full-Size Img PowerPoint

2.4   Temporal changes in resilience indices

To examine the first hypothesis, i.e., that 20th century environmental changes modified the sensitivity of growth to its drivers during severe events, we applied hierarchical models separately for each of the decades from 1900 to 2010 and for each event type. We incorporated descriptors of concurrent soil drought and VPD severity (${\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{d}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{d}}$, and Season of strongest drought) and the pre-drought and post-drought 4-year cumulative SPEI and average VPD (${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{v}}$, ${\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{v}}$ and ${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$, ${\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$).

where $R_{\mathrm{s}}$ and $R_{\mathrm{t}}$ are the resistance indices. Resistance is defined as the ratio of observed growth during the drought year to the growth in the same year predicted from the autocorrelation property of the site-specific time series. Equivalently, resilience is defined as the ratio of observed growth in the four years following the drought to the predicted growth in the same four years. ${\gamma }_{0,i}$ is the random intercept effects of the i site nested within species, while ${\beta }_j$ are also random coefficients varying by decade j, crossed with species/site. Three separate models were fitted, one for each type of severe event (low-SPEI, high-VPD, compound). Random effects were tested with the ranova command library lmerTest, and were found to be always very highly significant (i.e., at least P ​< ​0.0001). Trends of changes in the random ${\beta }_j$ coefficients of Eqs. 1 and 2 across decades were examined with weighted least square linear models (with or without a quadratic term for decade to examine curvilinear responses) as well as with weighted generalized additive models (GAM). The weights were set to the inverse of the variance of the coefficients for each decade, because variance accounts both for unequal sample size as well as unequal precision over time. The GAM models were employed primarily to confirm that the quadratic models were able to capture the trends in the data. Finally, we considered the potential confounding effect caused by varying contributions from different sites during the 20th century, because different chronologies started and ended in different years. To account for this effect, we repeated these analyses by subsetting the database only to those 429 series that began before 1900 and ended after 2000.

To examine the second hypothesis, i.e., that the effects of sensitivity shift on resistance and resilience differ depending on the type of the event considered, we employed Eqs. 1 and 2 to predict resistance $R_{\mathrm{t}}$ and resilience $R_{\mathrm{s}}$ separately for each decade in the century. As above, models were run separately for each event type. We used two complementary approaches to isolate the effects of changes in sensitivity to drivers on resilience indices from the changes in drivers themselves. In the first approach, because drought severity varied from decade to decade, we calculated the value of the differences of resilience indices estimated using decade-specific sensitivities (as per in the first hypothesis) minus resilience indices estimated using reference-decade sensitivities, while letting VPD_prev/d/post and SPEI_prev/d/post vary from year to year. As for the reference decades, we take the sensitivities of 1901–1910 for the period 1911–1960, and those of 1961–1970 for the period 1971–2010. For each decadal period, we test whether resilience indices predicted from variable sensitivities are larger or smaller than predicted from constant sensitivities. We split the century at 1960 and chose the first decade of each period as the reference decade based on the results of the previous analyses (Fig. 3), showing that many sensitivities to the environmental factors changed following a parabolic function with an inflection point at or after the middle of the century. This also matches trends of increased VPD which began to become apparent from the subsequent decade (1970s).

Figure  3.  Changes in the sensitivities of growth to VPD and SPEI by decade of the 20th century for resistance (top half, a–d) and resilience (bottom half of the graph, e–j). Sensitivities are the slope of standardized resilience indices to the two environmental drivers. P values refer to the significance of the regression of each sensitivity coefficient against decade for the three severe event types, using a linear or a quadratic model (cf., Fig. S4 for GAM models results), with the linear model preferred if the AIC of the quadratic model remains within two units of the linear model. Bars are 95% confidence intervals for each point. For each response variable, suffixes prev, d and pp refer to variables measured in the 4 years prior to the event, during or in the 4 years following the event year.

DownLoad: Full-Size Img PowerPoint

In the second approach, mean decadal values of the environmental drivers (SPEI_prev/d/post and VPD_prev/d/post) were calculated separately by site within species, year of the decade and type of severe event, and eventually averaged across sites and species to correspond with the decadal estimates of the sensitivity coefficients. We then selected one reference decade and recycled its values of the environmental drivers for the entire century. We chose to recycle the decade 1991–2000 because it is representative of a period where increases in VPD and compound events are already apparent (see Fig. S6, calculations were also made using the mean climatic conditions of the 1901–1910 and the 1951–1960 decades). We then used all the decadal-values of the sensitivity coefficients from 1901 to 2010, together with the recycled 1991–2000 climatic factors to calculate century-long resilience indices at constant climate.

Therefore, in the first analysis, climatic data vary through the century and one compares growth under constant versus variable sensitivities. Conversely, in the second analysis climate is constant and only the decadal variations in sensitivities affect resilience indices. Overall, the first approach asks the question of how varying sensitivities affect resilience indices relative to constant sensitivities under the varying climatic conditions of the 20th century. The second approach asks the question of how resilience indices varied at constant climate in response to varying sensitivity.

For the first approach, we compared the paired differences between the two growth rates across all sites with paired Kruskal-Wallis tests using library coin and the command independence_test, using a conditional distribution for the null against an asymptotic approximation of the exact distribution of the statistics. Analysis of the same data after aggregation of all events to species level led to the same conclusions. We confirmed this analysis by running a linear mixed model of predicted growth based on variable sensitivities against predicted growth based on fixed sensitivities for the reference decade, decade, event type and the interaction between these last two variables. For the second approach, we examined whether growth at constant climate changed through the century using linear, quadratic and GAM regression models as described above for the first hypothesis.

To examine the third hypothesis, i.e., whether plasticity in growth responses was sufficient to prevent long-term trends of growth reductions during or after severe droughts, we employed a mixed model framework similar to the one of Eqs. 1 and 2. In this analysis we pooled all the data for the century and described the linear responses of resistance and resilience to time (Year), with and without the effects of environmental drivers. The following two models were therefore employed:

where $R_{\mathrm{x}}$ stands for either $R_{\mathrm{s}}$ or $R_{\mathrm{t}}$ (both centred to mean of 0 and normalized to a standard deviation of 1), Year is the centred and normalized year in the century and Severe is a three-level variable (VPD, SPEI or compound events).

where ${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{x}}$ and ${\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{y}}$ stand for ${\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{d}},{\mathrm{V}\mathrm{P}\mathrm{D}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{v}},{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{d}},{\mathrm{S}\mathrm{P}\mathrm{E}\mathrm{I}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{v}}$, in the case of regressions of $R_{\mathrm{t}}$, whereas they stand for $\mathrm{V}\mathrm{P}{\mathrm{D}}_{\mathrm{d}},\mathrm{V}\mathrm{P}{\mathrm{D}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{v}},\mathrm{V}\mathrm{P}{\mathrm{D}}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}},\mathrm{S}\mathrm{P}\mathrm{E}{\mathrm{I}}_{\mathrm{d}},\mathrm{S}\mathrm{P}\mathrm{E}{\mathrm{I}}_{\mathrm{p}\mathrm{r}\mathrm{e}\mathrm{v}}\mathrm{o}\mathrm{r}\mathrm{S}\mathrm{P}\mathrm{E}{\mathrm{I}}_{\mathrm{p}\mathrm{o}\mathrm{s}\mathrm{t}}$ in the case of $R_{\mathrm{s}}$.

In Eqs. 3 and 4, ${\beta }_1,{\beta }_{\mathrm{x}}\mathrm{a}\mathrm{n}\mathrm{d}{\beta }_{\mathrm{y}}$ are now global coefficients for the century, not decade-specific. Eq. 3 provide an estimate of overall linear time trends, Eq. 4 of time trends net of the effects of year-specific environmental conditions. In these models, ${\gamma }_{0,i}$ is the random intercept of sites nested within species. Trends of changes in the resilience indices of Eqs. 3 and 4 were also examined with hierarchical GAM to confirm the linearity assumption inherent in Eqs. 3 and 4. All analyses were conducted in R (v. 3.5.1) (R Core Team, 2018). Mixed models corresponding to Eqs. 1 and 2 were fitted with lmer, those corresponding to Eqs. 3 and 4 with nlme. Library arm was used to obtain estimates of standard deviation of the random effects. Linearity in responses, model convergence, singularity, multicollinearity, normality of fixed and random effects, violations of homogeneity of variance, homoscedasticity and overall goodness-of-fit were examined using library performance (Lüdecke et al., 2020). We employed library BestNormalize (Peterson and Cavanaugh, 2019) and choose the Ordered Quantile Normalization (OQN) to normalize all variables (including the response variables of growth) to a mean of 0 and standard deviation of 1 (cf., further discussion on why this was required in Supplementary Materials). The OQN is a one-to-one transformation technique that is guaranteed to convert any arbitrary distribution to a normal distribution. It can be shown that if a dataset fits a pre-defined theoretical distribution of any kind, the OQN will produce a transformation that will approximate that distribution (Peterson and Cavanaugh, 2019). Significance of each variable was tested using the anova command using type III sums of squares. Marginal (fixed factors) and conditional (fixed and random factors) proportions of explained variance were obtained using the r.squaredGLMM function in the MuMIn package (Barton, 2009). The packages effects (Fox, 2003) and emmeans (Lenth, 2016) were employed to plot interactions and test differences across levels of categorical factors.

3.   Results

To examine our three hypotheses on the significance of plasticity in growth responses of Gymnosperm trees to drought, we compare growth rates measured during or in the four years following droughts (defined based on high-VPD, low-SPEI or compound events) against expected growth from the ARIMA model. In other words, the prediction is obtained from the autocorrelation properties of the time series (cf., Methods). We purposely avoid traditional resistance/resilience indices based on measured growth rates prior to droughts because our earlier results showed that these metrics are often not independent, preventing the interpretation of the two indices in isolation (Zheng et al., 2021; Serra-Maluquer et al., 2022). We comprehensively examine the literature on resilience and find that, almost universally, its definition implies non-independence from resistance (Supplementary Materials, Sections 1 to 3, Figs. S2 and S3; Table S3‒S5). Because resistance and resilience define conceptually different dimensions of vegetation responses to stresses (Nimmo et al., 2015), avoiding non-independent definitions is essential. The use of ARIMA models (Martínez-Vilalta et al., 2012) achieves this goal perfectly. Using autocorrelation features to derive estimates of ‘expected’ growth is consistent with the definition of resilience (MacGillivray et al., 1995) and with our own earlier work on growth responses to drought events (Martínez-Vilalta et al., 2012; Ovenden et al., 2021).

3.1   Plasticity of the sensitivities to environmental drivers

We first examine hypothesis 1, that during the 20th century, growth sensitivities to environmental factors are not stationary for the three types of high-VPD, low-SPEI and compound events. To test this hypothesis, we calculate the response functions of resistance and resilience as a function of VPD and SPEI as predictor variables for each decade of the 20th century in linear mixed models (Eqs. 1 and 2). For resistance (Fig. 3a‒d), three out of the four sensitivities to VPD and SPEI vary significantly (the one for VPD_d only marginally) throughout the 20th century, suggesting plasticity in growth sensitivity to some environmental drivers. Specifically, the significant quadratic trend for SPEI_d (the annual SPEI during the drought year, Fig. 3b) suggests a positive effect on resistance during severe events for the majority of the 20th century with a minimum around 1960s, i.e., other things being equal, positive SPEI_d coefficients correlate with increased resistance. The signs of VPD_d (the summertime VPD during the drought year) are mostly negative for high-VPD and compound events (Fig. 3d, averages of −0.11, −0.09 and 0.03 for compound, high-VPD and low-SPEI events, respectively), indicating that higher VPD_d reduced resistance and drought-year growth, while the signs for the period preceding the severe events (VPD_prev) are expected to vary in the opposite direction. Overall, the magnitudes of the changes in sensitivity are substantially larger for compound than high-VPD or low-SPEI events (mean ranges, i.e., max-min, of 0.56, 0.45, 0.36). To relax assumptions of linearity, we also employed weighted GAM (generalized additive models) regressions of sensitivities against decade by extreme event type instead of linear regression. Similar results are obtained (Fig. S4a‒d).

For resilience (i.e., observed divided by expected post-drought growth), four out of six of the sensitivities to VPD and SPEI significantly vary throughout the 20th century (Fig. 3e‒j), suggesting stronger plasticity in growth sensitivity to individual environmental drivers following severe events. In two instances, sensitivities vary nonlinearly over time with an inflection point around the middle of the 20th century (AIC of quadratic versus linear regression « −2). The magnitudes of these changes are substantially larger for compound events relative to high-VPD and low-SPEI events (mean ranges of 0.73, 0.37, 0.35, respectively). For all three kinds of droughts, negative correlations are found between SPEI_prev and R_s at early and late 20th century, however, the effect all turn positive during the middle of the century (Fig. 3e). Decreasing trends are found for the effect of SPEI_d on R_s in low-SPEI and compound events, while its effect in high-VPD events is almost unchanged (Fig. 3f). The values of slope for R_s to SPEI_post in all drought types are always above 0, indicating that higher post-drought SPEI always leads to higher post-drought growth, yet the sensitivity in high-VPD and compound events decline with the time and that for low-SPEI events is stable (Fig. 3g). The majority of the slopes of R_s to VPD_d are below 0, expressing that high VPD_d usually result in low post-drought growth, the sensitivity of R_s to VPD_d in three kinds of droughts all shift considerably (Fig. 3i). Most significantly, for the four variables where significant trends are seen (SPEI_prev, SPEI_d, SPEI_post and VPD_d), the directions of the changes during the second half of the century are always towards a reduction in the sensitivity of resilience to changes in the drivers, especially so for high-VPD and compound extreme events. Overall, considering both resistance and resilience, the sensitivities to SPEI metrics are substantially more plastic than those of VPD metrics (Fig. 3a‒j).

In all these analyses, the random effects of species (and sites within species) are almost always highly significant (often P ​< ​0.0001) for both growth metrics, suggesting that the magnitude of the temporal trends vary systematically from species to species and site to site. When linearity assumptions are relaxed using GAM models, similar changes in the sensitivity to environmental drivers are also seen (Fig. S4 bottom). When the analysis is repeated using only the 429 series that started prior to 1990 and ended after 2000, similar results were also found (Fig. S5), albeit with lower significance level. This confirms that responses to VPD and SPEI during and following severe events changed significantly and often curvilinearly throughout the century.

3.2   Effects of shifts in sensitivities to environmental drivers

Our second hypothesis predicts that, as a consequence of shifts in sensitivities, resilience indices vary, with the variation dependent on the nature of the severe event, i.e., whether high-VPD, low-SPEI or compound event. For this analysis, we couple two approaches that allow isolating the plasticity effect (i.e., eliminating the climatic effects) in two complementary ways (as described in 2.4).

Using the first approach, for the period 1901–1960, both resistance (Fig. 4a top line) and resilience (Fig. 4b top line), showed distributions that are significantly skewed towards negative values (except R_t in low-SPEI events). These negative distributions are more accentuated during high-VPD and compound years than during low-SPEI years, especially for resilience (P ​< ​0.0001 using paired t-tests within event types, Fig. 4b). The predominant negative departures during drought events suggest that with changes in growth sensitivity, both resistance and resilience with varying sensitivities was lower than it would have been with the constant 1901–1910 sensitivities. Keeping in mind the results of R_t in Fig. 3, considering the relatively larger changes of sensitivities to environmental drivers in high-VPD/compound events taken place during 1901–1960, the most likely interpretation is that the additive effects of the changes in sensitivities 1–4 years prior to and during high-VPD/compound events increased the negative effects of drought on growth, which result in the reduction of resistance. For resilience, in addition to the previous reasons, a significant decline in the sensitivity of SPEI_post may also account for the decrease (Fig. 3g). Finally, for low-SPEI events, changes are small for resistance, but a larger negative shift in the distribution is seen for resilience.

Figure  4.  Analysis of the effects of plasticity on the sensitivity parameters (pars) of growth rates. The time-varying sensitivity parameters are calculated for each decade as in Fig. 3. The x axis plots the value of the differences of resilience indices predicted based on decade-specific sensitivity parameters ($\beta$ calculated from Eqs. 1 and 2 as we did in subtitle 3.1) minus the same indices predicted based on a reference-decade sensitivity parameters ($\beta$ fixed at the reference decade), for each severe event. The same values of environmental variables are employed in both cases to scale the resilience indices. Results are given for resistance R_t (panel A) and resilience R_s (panel B). Within each panel, the first line refers to the 1901–1960 period (reference decade 1901–1910) and the second line to the 1961–2010 period (reference decade 1961–1970). From left to right, results are given for high-VPD, low-SPEI and compound events, respectively. Values to the left of the vertical dashed line at 0 indicate that decade-specific sensitivities led to smaller resilience indices compared to reference-decade sensitivities. The Delta and P values in each figure give the mean difference and probability of the Kruskal-Wallis test of the differences between the two growth estimates by pooling all events of each period. Results are almost identical when site-level data are averaged by species.

DownLoad: Full-Size Img PowerPoint

Conversely, for the 1961–2010 period, resistance (Fig. 4, second line of panel a) shows only very small shifts in the distribution for VPD events and moderate negative changes for SPEI and compound events. This suggests that, despite changes in sensitivities after 1960, plasticity did not improve growth responses, and drought-year growth remained even lower than at the reference decade for all types of drought events. Resilience shows instead positive shifts in the distributions for VPD and compound, but not SPEI, years (Fig. 4b, second line). These results imply that plasticity in growth sensitivities increased growth especially following high-VPD and compound events towards the end of the century. Because we pooled data within each period, we also ran mixed model regressions against Decade and Event type to confirm these conclusions (Table S6).

We then use our second approach, to further isolate the effects of temporal changes in the sensitivity coefficients from climate variability. We predict resilience indices at constant climatic conditions (taken from the 1991–2000 decade) but with varying sensitivities. To facilitate interpretation, we draw horizontal lines in our figures corresponding to the unnormalized values of R_t ​= ​1 (corresponding to observed growth equals expected growth during the drought) or R_s ​= ​1 (corresponding to observed growth equals expected growth following the drought), but it needs to be remembered that resilience indices for evaluating the effects of plasticity have been normalized. Temporal changes in the sensitivity coefficients at constant climatic conditions do not cause systematic trends in predicted resistance (P ​> ​0.05 for the effect of decade and the interaction between decade and event type, Fig. 5a). The lack of a temporal trend in predicted resistance at varying sensitivities and constant climate is consistent with the limited plasticity in the sensitivity coefficients for resistance through the century (Fig. 4a).

Figure  5.  Changes in the resistance (a, R_t) and resilience (b, R_s) by decade of the 20th century. In this analysis, climate was kept constant at the severity values for all sites using the data from the decade 1991–2000, whereas the sensitivity coefficients were allowed to vary by decade (cf., Fig. 3). P values refer to the significance of F-tests for decade from the regression ANOVAs, and separate lines are given for the three types of severe event years. The combined regressions against decade × Severe event type had R² = 0.13 (AIC of parabolic relative to linear model = 1.6) and R² = 0.50 (AIC of parabolic relative to linear model = −7.7), for resistance (a) and resilience (b), respectively. Bars are 95% confidence intervals for each point. To facilitate legibility, data points are dodged 3 units to the right side from the beginning year of the decade. Black points are averages across all event types. R_t and R_s are standardized at mean = 0 and standard deviation = 1. The thin black lines give the values at which unnormalized R_t = 1 and R_s = 1. Similar analyses were carried out using climatic data from the 1901–1910 and the 1951–1960 decades (Fig. S6).

DownLoad: Full-Size Img PowerPoint

In contrast, for resilience, plastic responses in the sensitivity coefficients at constant climatic conditions are sufficient to cause substantial increases for the final third of the century (Fig. 5b). Interestingly, using a quadratic regression model, the temporal recovery is identical in magnitude for all three event types, without a significant interaction term (P > 0.05) and despite a significant difference in absolute values across event types (P ​< ​0.05). Fig. 5 demonstrates that what might seem as relatively minor trends in growth sensitivity coefficients (Fig. 3) result in substantial differences when the effects of individual coefficients are added up. As expected from Fig. S5, the substantial increases in resilience at constant climate is also confirmed when data from only the 429 century-long series are considered (Fig. S6).

Considering together the results of Figs. 4 and 5 suggests a substantial difference between low-SPEI events on the one hand, and high-VPD and compound events, on the other, in how resilience varied in responses to both climate and plasticity in sensitivities. For low-SPEI events, the increase in resilience at constant climate post-1960 (Fig. 5b) entirely disappears when decadal-related climatic differences are considered (Fig. 4b, the middle of second line). In other words, the plasticity of many of the sensitivity coefficients (Fig. 3, bottom panel) translates into positive adjustments in resilience only under constant climatic conditions for low-SPEI events. By contrast, for compound and even more high-VPD events, the plasticity in growth remains whether resilience is compared at the constant climate of the 1991–2000 decade (i.e., the second approach, Fig. 5b) or by employing actual post-1960 climate. Furthermore, we should notice that the adjustments of sensitivities can occur not only at individual, but importantly also at stand scale. In this case, the historical trends in resilience during the 20th century would be biased towards zero, when considering the fact that dead trees were not involved, i.e., the true effects of historical shifts in sensitivities would be larger than we report.

3.3   Long-term trends in resilience indices

The evidence above addresses hypotheses 2 and 3, i.e., that physiological adjustments reduced, but did not eliminate, the negative effects of severe droughts on resilience towards the end of the 20th century. We further examine hypothesis 3 by looking at the long-term trends in resistance and resilience throughout the 20th century, while taking both trends in environmental conditions as well as potential plasticity in growth responses into account. In this analysis, we pool data from all decades and fit linear mixed models as a function of year, or year plus environmental conditions. Resistance declines throughout the century only for high-VPD years (P ​< ​0.05 for the significance of the linear decline, and P ​< ​0.0001 for the difference in slope between low-SPEI and high-VPD years) (Fig. 6a; Table S7), while trends for low-SPEI and compound events are flat. The magnitude of the response varies across event types (P ​< ​0.0001), with lower values for compound and low-SPEI relative to high-VPD years (P ​< ​0.0001, Tukey post-hoc comparisons significant level ​= ​0.05, for both cases). Interestingly, while resistance predicted during high-VPD years is substantially higher than for the other event types during the earlier part of the 20th century, the difference disappears towards the end of the century. The inclusion of environmental covariates in the model does not alter these time trends (Fig. 6b).

Figure  6.  Results of linear mixed models of the changes of resistance (top line) or resilience (bottom line) as a function of year in the 20th century. In each line, the panels give the trends as a function of year only (model I, left-hand-side), or as a function of year plus SPEI and VPD drivers (model II, right), separately for the three types of severe events. Each panel gives the ANOVA results of the F-tests for year, event type and the interaction between the two, ns is short for not significant. Bands around each line are one standard error.

DownLoad: Full-Size Img PowerPoint

Resilience declines throughout the century for all drought types (P ​< ​0.0001 for the significance of the linear decline, Fig. 6c and d; Table S8). Even though compound events tend to have a significantly steeper slope, the difference across event types is not significant, i.e., P ​= ​0.35). The magnitude of the response varies across event types, with higher values for low-SPEI relative to high-VPD and compound years (P ​< ​0.0001, Tukey post-hoc comparisons significant level ​= ​0.05, for both cases). The difference between time trends is consistent with results in Figs. 4 and 5, in showing a flat or marginal declining trend for resistance and a consistently declining trend for resilience. Similar results are obtained using GAM models (Fig. S7). We note again, however, that these results focus on surviving trees, and therefore result from a combination of individual-and stand-level adjustments.

4.   Discussion

4.1   Effects of shifts in sensitivities

We show that the sensitivities of tree growth to several environmental variables during severe drought years has changed throughout the 20th century in a direction that reduced the sensitivities to external drivers during high-VPD and compound events. Besides, these changes have also altered the resilience indices independent of climatic changes. Finally, the adjustments in these responses in the final third of the century substantially reduced the negative trends in resilience of Gymnosperm trees. The finding of plastic responses to environmental factors was the same, independent of which specific ARIMA imputation approach was employed. It also remained the same, independent of whether we included all the time series or we subset only those series that spanned the entire century. Both of these checks suggest that it is a robust finding. The consequences on growth of these adjustments in sensitivities were visible during all types of severe events, if climate conditions were kept constant at 1991–2000 levels. However, under realistic climatic conditions of the 20th century, these adjustments in sensitivities did not alter resilience following years of low-SPEI events. Altogether, our results are consistent with the hypothesis that more frequent high-VPD years affected growth rates negatively after the 1970s, but also elicited compensatory responses by trees that reduced these negative effects. These compensatory responses are largely driven by changes in how drought legacies affect responses to drought.

Our comparative analysis of resilience indices highlights a conceptual problem whose significance had been overlooked. Despite their widespread use and repeated invitations to examine the different components of resilience (Orwin and Wardle, 2004; Lloret et al., 2011; Nimmo et al., 2015; Ingrisch and Bahn, 2018), spurious correlations between resilience and resistance indices were almost universal. Spurious correlations complicate the analysis of the different physiological or ecological processes leading to resilience. We show that it is possible to avoid spurious correlations in the analysis of the components of resilience via a reference null model. Conceptually, our approach is equivalent to (MacGillivray et al., 1995), originally intended for experimental comparisons against a control. Our null model maintains the temporal memory of the process via its autocorrelation structure (Peltier et al., 2016; Ovenden et al., 2021) and differs therefore from (Anderegg et al., 2015; Jiang et al., 2019), where null expectations of growth were obtained from static growth-climate relationships.

Up to the decade of the 1960s, plasticity in responses following severe VPD and SPEI years (Fig. 3) was associated with declines in resilience and partly resistance (Fig. 4, Table S6). Large declines in resilience were also seen when we replicated the 1991–2000 climate conditions across the entire 20th century (Fig. 5), confirming the casual link between plasticity in sensitivity and growth responses. Although the rate of decline in resilience predicted under the constant climatic conditions of the 1991–2000 decade did not vary significantly across all three event types, values were generally lower for high-VPD and compound events, suggesting substantial greater impacts of VPD compared to SPEI on post-drought growth across all sites and species. The causes of this substantial decline in growth caused by altered sensitivities need to be investigated further. One can speculate that this may have resulted from canopy overshoots, one possible mechanism being CO₂-driven increases in canopy size, leading to wrong-way responses and amplification of competition for water during/after droughts (Jump et al., 2017; Zhang et al., 2021).

These effects on resilience dramatically reversed during the second half of the 20th century, especially but not only for high-VPD and compound years, when changes in sensitivity to environmental variables prevented further losses in post-drought growth. Changes during the second half of the century support the hypothesis that a higher frequency and/or severity of high-VPD events after the 1970s–1980s elicited a positive memory response, and therefore that physiological adjustments took place. This interpretation is consistent with the results of both analyses where the effects of decadal changes in sensitivity on growth were examined, i.e., either by varying the climate through the 20th century, while comparing sensitivities against a reference decade (Fig. 4), or by varying sensitivities at constant climate (Fig. 5). Our mixed models assumed linearity in the responses of growth to environmental drivers. To test whether responses may have varied through the range of possible conditions, we tested the constant climate scenario of Fig. 5 by replicating instead the climate data from the 1901–1910 and the 1951–1960 decades. We found almost identical patterns, both with regard to absolute differences across event types as well as to opposing trends during first and second parts of the century (Fig. S8). In additional preliminary analyses, we also employed GAM models to directly explore the linearity in the responses to environmental drivers and found very little evidence of nonlinearities. We are therefore confident that our analysis captures a genuine shift in the responses of trees to environmental and endogenous drivers.

Table 1 summarises all estimates of resilience losses avoided thanks to plasticity obtained in this analysis. These results suggest substantial resilience adjustments in the last decades of the 20th century, implying that plasticity in sensitivity made a marked difference to growth during high-VPD and compound years. Altogether, these findings suggest that acclimation to more frequent and severe high-VPD years may have occurred during the final part of the 20th century.

Table  1.  Summary of the estimates of the role of plasticity in avoiding resilience losses during the second half of the 20th century. Estimates are obtained based on different assumptions using the results given in the corresponding figures and tables. The first line in the table provides an estimate of plasticity having accounted for climatic conditions in different ways. The second line gives the estimate of avoided losses in resilience by relating the plasticity changes of the first line to the observed net temporal losses of Fig. 5c. Data from Fig. 3 are obtained from the model in Table S6 as the linear slope of Decade for compound events for the post-1960 period. Estimates of plasticity effects from Figs. S6 and 4 are obtained as the linear slope of Decade for compound events for the post-1960 period for the 1901–1910, 1951–1960 and 1991–2000 constant-climate scenarios, respectively. Because resilience indices are standardized with mean of 0 and standard deviation of 1, changes due to plasticity are given as % of 1 standard deviation.

Index	Source
	Fig.3/Table S6	Fig. S6/1901–1910 climate	Fig. S6/1951–1960 climate	Fig.4/1991–2000 climate	Mean effect/year ±standard error
Post-drought	1.96%/decade	0.2%/decade	0.7%/decade	2.3%/decade	1.29% ± 0.50%
Plasticity effect on growth	87% avoided losses	41% avoided losses	70% avoided losses	89% avoided losses	72% ± 11%

 | Show Table

DownLoad: CSV

4.2   Drivers of the plasticity in sensitivity

Although our data allow to analyse the variations over time in the sensitivity of resilience indices to VPD and SPEI during the 20th century, identifying the factors that might have driven this plasticity is more complex. Several processes may explain adjustments to severe drought years during the latter part of the 20th century, i.e., lower stomatal sensitivity to soil drought or VPD (Poyatos et al., 2007), higher sapwood-leaf area ratios (Rosas et al., 2019; Trugman et al., 2019), direct temperature responses in growth (Cabon et al., 2020b) and altered whole-tree carbon balance (Dietze et al., 2014). Considering the directions of the changes in the sensitivity coefficients during the second half of the century and the fact that the most important shifts in the sensitivity coefficients occurred for SPEI, rather than VPD, metrics (Fig. 3), one might conclude that trees have become less sensitive to changes in soil drought drivers under high VPD conditions. Responses to changes in VPD over time may be confounded by covariation between temperature and VPD. We believe that this is unlikely, because of the strong mechanistic link between growth and soil and atmospheric drought (De Schepper and Steppe, 2010; Cabon et al., 2020a). Even when temperature has been shown to control cambial sink activity, this appears to be mediated by cell turgor effects (Peters et al., 2021). It cannot be excluded that long-term declines in nitrogen availability (globally evident for most of the 20th century in some cases (Mason et al., 2022)) contributed to reduced nitrogen concentration in leaves and therefore helped with responses to extreme events. In any case, our central argument remains, that plasticity in growth sensitivity during severe VPD and compound-event years reduced but did not entirely prevent reductions in resilience during the second half of the 20th century.

The approach we took here is similar to the comparison of observed growth versus growth predicted under reference null scenarios (Anderegg et al., 2015; Jiang et al., 2019), although previous studies did not examine temporal plasticity. Several hypotheses have been put forward to explain variable growth-climate sensitivity (Peltier and Ogle, 2020). The hypothesis of changes in climate change drivers (Levesque et al., 2017; Peñuelas et al., 2017) is consistent with the significant directional shifts in sensitivities (Fig. 3). However, as mentioned above, plastic changes in the sensitivity to VPD and SPEI do not necessarily reflect a direct response to historical changes in these two abiotic drivers. Even so, the intensity of both low-SPEI and high-VPD events varied through the 20th century, but in different directions (towards less strong SPEI, but stronger VPD events, Fig. S9), while their frequency increased on average (by 4% and 18% for SPEI and VPD, respectively), between 1901‒1950 and 1951–2000. Attributing observed plastic changes to individual environmental or biotic drivers will require a separate analysis, but we note that some regional studies have attributed recent mortality episodes of gymnosperm-dominated forests to temperature and VPD (Stovall et al., 2019). A recent global synthesis identified the imprint of hotter-drought events on reported mortality episodes (Hammond et al., 2022). The hypothesis of directional changes in plant physiological processes (i.e., hydraulic damage, carbon source-sink dynamics, structural changes (Peltier and Ogle, 2020)) is also consistent with previous diachronic and synchronic analyses of the autocorrelation term in paired chronologies from surviving and dead gymnosperm trees (Cailleret et al., 2019; DeSoto et al., 2020).

In our analyses, the random effect of species was almost always highly significant, suggesting that the magnitude of these temporal trends varied from species to species. In a separate analysis using a slightly larger dataset (1931 site series), we showed that species having a stronger relationship between ring width index and SPEI drought episodes, and having low resistance and high recovery from low-SPEI years, tended to have higher wood density and more negative leaf water potentials, consistent with the interpretation provided above and suggesting that species hydraulic traits underlie some of these temporal responses (Serra-Maluquer et al., 2022). One should note that the vast majority of the tree growth records from the ITRDB come from natural stands, not plantations. We recognize that our conclusions would be stronger if site-specific information was available on stand conditions (e.g., type of management, species composition, density) of the stands from which the tree cores were taken. Unfortunately, this information is not available for the vast majority of the sites in ITRDB. This likely adds considerable noise (and potential biases, cf. below) to the results and prevents a more in-depth interpretation. We also recognize that additional insights into species-specific drought responses could be obtained by a more detailed examination of resilience patterns of genera known to differ for their average rooting depth (e.g., deep-rooted pines versus shallow-rooted spruce). Although this analysis was not possible for this study, we mention it as an obvious future expansion of the present project.

4.3   Long-term trends in resilience indices

Despite evidence of compensatory mechanisms, our data show that resilience did decline during the 20th century even when changes in VPD_prev/d/post and SPEI_prev/d/post were accounted for. This historical decline in resilience is consistent with an earlier analysis using a different resilience metric (Zheng et al., 2021), but contrasts with the conclusion using an impact-normalized index that quantifies recovery (Li et al., 2020), by taking the obtained level of drought stress as a reference, as opposed to resilience per se (Ingrisch and Bahn, 2018). Our results help to clarify the distinction between trends in recovery and resilience. We show that growth sensitivities changed plastically over time (Fig. 3) and that these plastic adjustments reduced the magnitudes of time-related declines in resilience post-1960 (Figs. 4 and 5). We therefore posit that these findings explain the increases in impact-normalized post-1950 recovery reported by (Li et al., 2020), i.e., increased recovery capacity via plasticity did not elicit complete adjustment and did not prevent long-term declines in actual resilience. This conclusion is robust whether calculated relative to measured growth prior to drought (Lloret indices (Lloret et al., 2011) in (Zheng et al., 2021)) or to expected growth that would have occurred in absence of droughts (this paper). We also show that this conclusion applies not only to low-SPEI but even more to high-VPD and compound events.

Publicly available tree ring data remain spatially limited (Zhao et al., 2019) and biased towards old or dominant trees of a limited number of species often growing under marginal conditions, for which strong climatic sensitivities can be expected (Klesse et al., 2018). The distribution of tree ages remained roughly constant through the 20th century, with a prevalence of relatively older trees (Fig. S10). Also, tree age was not found to be a significant variable in our models. To avoid potential confounding age effects, we excluded series with age lower than 90 years for each decade. However, when we did not exclude these data, similar trends were found. Similarly, the number of study sites and sample species (Table S2) only showed declines in frequencies in the decade after the year 2000 (Zhao et al., 2019). Because temporal adjustments in sensitivity can only be documented in trees that survive extreme events, it is possible that our analysis over-estimates the significance of plasticity and consequently under-estimates the decline occurring in trees that eventually lose resilience and succumb (Brienen et al., 2012; DeSoto et al., 2020). Improved information on stand-level responses to drought, which are not usually documented in the ITRDB, would facilitate a better understanding of the drivers of the adjustments observed at the tree-level. Overall, we show that the combination of negative effects caused by climate change forcing and subsequent adjustment of Gymnosperm trees during the 20th century can be conceptualised in terms of the net balance between the physiological damage caused by increased severity of extreme droughts, and trait plasticity and potentially acclimation preventing higher levels of damage, especially in response to high VPD. Understanding the nature and magnitude of these compensatory mechanisms in relation to current warming trends is crucial to predict whether trees will continue to head towards more challenging conditions during the 21st century.

5.   Conclusions

We studied the sensitivity of resilience indices of conifer trees to climatic factors (VPD and SPEI) during three types of droughts (high-VPD, low-SPEI, and compound events) in the 20th century. Our findings demonstrate that the sensitivities of resistance and resilience to VPD and SPEI changed remarkably throughout the 20th century, with the directions of these changes often reversing in the second half of the century. Changing sensitivities had positive effects on resilience, especially following years of high-VPD and compound events, avoiding growth losses that would have occurred if sensitivities had remained constant. Despite sensitivity changes, resilience was still lower at the end of the 20th century compared to the beginning of the century. Future adjustments to low-SPEI and high-VPD events are likely to continue to compensate for the trends in climate only partially, leading to further generalized reductions in tree growth of conifers during the 21st century. An improved understanding of these plastic adjustments and their limits, as well as potential compensatory processes at the stand level, is needed to project forest responses to climate change.

Author's contributions

TZ: Funding acquisition, Conceptualization, Data curation, Formal analysis, Methodology, Writing - original draft. JMV: Conceptualization, Methodology, Writing - review & editing. RGV: Conceptualization, Data curation, Formal analysis, Methodology, Writing - review & editing. AG: Conceptualization, Methodology, Supervision, Writing - review & editing. JJC: Conceptualization, Methodology, Supervision, Writing - review & editing. CM: Funding acquisition, Supervision, Writing - review & editing. MM: Conceptualization, Methodology, Formal analysis, Writing - review & editing.

Availability of data and materials

Data and code for analyses are available at: https://github.com/mmencuccini/ITRDB-data-analyses.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

We thank all dendrochronologists for providing data to the ITRDB. We thank China Scholarship Council and professor Changcheng Mu for providing the funding and opportunity.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.fecs.2023.100107.