The present study was conducted with samples taken from a monospecific spacing experiment of Pinus sylvestris L., established in 1935 and situated in the north-east of Scotland (57° 36′ 23″ N, 4° 16′ 50″ W) at 170m a.s.l and two spacing (density) treatments were selected for use: high density (ρ_H) and low density (ρ_L).

Sample plots were randomly established so that ten 0.02 ha plots were present in both ρ_H and ρ_L. Measurements of tree diameter at breast height (DBH – 1.3 m) were collected for all trees within the sample plots. Sample trees were then selected from within these sample plots using the diameter distribution of measured trees to represent three dominance classes (Dominant, Co-dominant, and Sub-dominant). Dominant trees were selected as those closest to, but above the upper quartile diameter value, co-dominant trees were selected as those closest to, but above the median diameter value and sub-dominant trees were those closest to, but above the lower quartile diameter value across all trees in the 10 sample plots in both ρ_H and ρ_L where present. This approach meant 14 dominant, 14 co-dominant and 6 sub-dominant trees in each treatment were felled in the winters of 2002 and 2003, resulting in 34 sample trees from both ρ_H and ρ_L. Cross sectional discs were taken from the base of each tree at 0.3m ± 30cm and scanned, resulting in 29 (ρ_H) and 27 (ρ_L) usable discs images for analysis (Table 1). Mortality assessments had been conducted in these plots approximately every 5 years from 1955 to 1990 (Supplementary Table 1).

Individual tree annual basal area increments (BAI) were calculated from raw ring width data following Eq.1 where R is the radius of the tree in years t and t-1 using the dplR package (Bunn et al., 2019) in R version 3.6.1 (R Core Team, 2019). See Ovenden et al. (2021) for further details of the site and dendroecological methods used to process the collected samples.

To identify notable drought years, we used interpolated climate data at 1 km resolution, obtained from the climate hydrology and ecology research support system (CHESS) meteorology dataset for Great Britain (Robinson et al., 2017) for the study period (1961–2002) to calculate monthly values for both the standardized precipitation evapotranspiration index (SPEI) (Vicente-Serrano et al., 2010) and the climatic water deficit (CWD). SPEI was calculated for August using the five previous months data (i.e., using a 6 month integration period of March–August) (SPEI_Aug6) across the study period using the SPEI package in R (Beguería et al., 2014). Increasingly negative SPEI values indicate increasingly severe drought conditions. Huang et al. (2015) identified an SPEI threshold of – 1.64 in other species of pine (Pinus edulis and Pinus ponderosa), below which drought was linked to significant growth decline. As such, –1.64 was the thresholds adopted here for defining a drought event using the SPEI (Figure 1).

Climatic water deficit was calculated for the site at monthly time steps between 1961 and 2002 following Lutz et al. (2010) using R code developed by Redmond (2019). The approach used by Lutz et al. (2010) to calculating CWD applies a Thornthwaite-type water-balance model (Thornthwaite, 1948) which is considered most appropriate when climatic data are limited (e.g., when data on historic wind speed, humidity, and net radiation are absent). As such, CWD was calculated here using slope (degrees), latitude (decimal degrees), folded aspect (degrees) (McCune and Keon, 2002), monthly total precipitation (mm), mean monthly temperature (°C), and the soil available water capacity (mm) in the top 200 cm of the soil. Soil available water capacity for the study site of 143.63 mm was obtained from the James Hutton Institute’s Available Water Capacity (AWC) dataset (Gagkas et al., 2019). Monthly CWD values between March and August were then summed annually to give a single annual CWD value for each year during the study period (Figure 1).

Only years that showed a clear drought signal in the calculations of both SPEI and CWD were selected for analysis, resulting in the identification of three drought events in 1976 (SPEI_Aug6 of – 1.8, CWD of 60.2 mm), 1984 (SPEI_Aug6 of –2.0, CWD of 66.0 mm) and 1995 (SPEI_Aug6 of – 1.9, CWD of 58.5 mm), meaning that both SPEI and CWD indicate that 1984 was the most severe drought in the study period (Figure 1).

We applied a set of exclusion criteria to SPEI_Aug6 values for the study site across all years preceding each of the three drought events to remove years that indicated unusually dry conditions or potential drought recovery years. This was necessary to ensure that only BAI values from non-drought years or drought recovery years were used to calculate pre-drought growth averages, against which drought responses would be subsequently compared. We then used the corresponding BAI values for the retained SPEI- constrained years to calculate a cumulative moving average (CMA) BAI value at an annual time step (BAI_con). In this way, BAI_con changes over time as additional SPEI-constrained growth years are included in the average. Initially, we identified all years during the study period (1961–2002) where the SPEI_Aug6 values were between –1 and 1 (Figure 1). This threshold is commonly used in drought studies to indicate non-drought conditions (Hoffmann et al., 2018; Vanhellemont et al., 2018; Slette et al., 2019; Bose et al., 2020). In order to avoid including BAI values that might represent growth during a recovery year [i.e., the year(s) after SPEI values < –1] in the calculation of an SPEI-constrained growth average, we applied thresholds for data exclusion by omitting 1 year of BAI after an observed SPEI value of ≤–1 but >–1.5, 2 years after a SPEI value ≤ –1.5 but > –2.0 and 3 years if the SPEI value was ≤–2, reflecting commonly used thresholds of drought severity (Hoffmann et al., 2018; Vanhellemont et al., 2018). This approach was designed to be conservative (by excluding growth years that may themselves be a recovery year) and reflect the fact that higher stress levels may be more likely to result in physical damage such as hydraulic failure (Adams et al., 2017) and potentially slower recovery (Gessler et al., 2020). Thus, BAI_con was calculated annually from the CMA of all prior years that had not been excluded using one of the above criteria. In this way, SPEI-constrained growth years were defined as having an SPEI_Aug6 value of > –1, with CMA BAI calculations also omitting any BAI values where growth was potentially recovering from a previous dry year. This approach meant that for each tree separately, the BAI1976_con for the 1976 drought was calculated using nine prior annual growth records, the BAI1984_con for the 1984 drought was calculated using 13 prior growth records and the BAI1995_con for the 1995 drought was calculated using 21 prior growth records from the study period (Figure 1).

Resistance (Rt) is a measure of the immediate impact of drought on tree radial growth relative to a pre-drought growth average calculated over a pre-defined period, often between 2 and 5 years before drought (Lloret et al., 2011). As a result, Rt is calculated following Eq. 2:

where, PreDr is average pre-drought growth and Dr is growth during the drought year, providing a measure of drought impact as the ratio of drought year growth performance relative to a pre-drought growth average. We made a simple modification to Eq. 2 which allowed us to express Rt as PGC following Eq. 3:

where, the pre-drought growth average in the Rt calculation is replaced with individual tree BAI_con values (BAI1976_con, BAI1984_con, and BAI1995_con) and Dr is replaced with BAI_t which represents the observed BAI in any year (during or after drought). The use of Eq. 3 over Eq. 2 enabled us to assess BAI_t annually against an SPEI-constrained pre-drought growth average for individual drought episodes. The use of Eq. 3 also allowed an assessment of recovery, with an annual PGC value of 0% indicating growth rates have returned to BAI1976_con, BAI1984_con, or BAI1995_con levels. We calculated PGC for the three drought years and for the 4 years following each drought but could not calculate PGC for the years immediately preceding each drought event due to the inclusion of some of these years in the calculations of the pre-drought growth averages (BAI1976_con, BAI1984_con, and BAI1995_con) (Figure 1). We also calculated PGC annually over the same period for all three drought events using ring-width data that had been detrended using a cubic smoothing spline with a 30-year cut-off to ensure our results derived from BAI data were robust.

We fit a mixed-effects model using the nlme package in R (Pinheiro et al., 2020) following Eq. 4 to investigate at an annual resolution whether P. sylvestris trees show differences in growth (indexed by BAI) after droughts of differing severity, at different stand densities and to ascertain how long any differences might persist post-drought.

In Eq. 4, BAI_ij is the jth measure (representing measurement at each j year) of BAI for the ith tree, X is an n x p matrix of fixed effect variables, where n = the number of observations and p = the fixed effect variables, including drought event, stand density and timepoint (year), with an interaction between drought event and timepoint (year), β is a p x 1 column vector of regression estimates, b0_i represents the random effect of tree, where b0i ∼ N(0, σ²₀) and ε represents the residual error term, where ε_ij ∼ N(0, σ²) (where the error terms are assumed to be normally distributed with a mean of 0 and standard deviation σ²).

Basal area increments values were square root transformed prior to analysis to improve model fit. In order to compare the three drought events, the year of the drought was assigned timepoint “0” (1976, 1984, or 1995). The next four recovery years were then described as years 1–4. Timepoint (year) (zero to four, where zero is the drought year and one to four are the post-drought years) was fitted using third order orthogonal polynomials to accommodate for non-linearities in BAI over time, with the optimal degree of polynomials selected by minimizing AIC values. Tree ID was fitted as a random effect and interaction terms between drought event and year were included. To correct for temporal autocorrelation, the correlation structure was modeled using a corARMA correlation structure of p = 2, q = 0 and year as a time covariate, which also served to detrend the BAI data.

We subsequently compared the association between BAI and the three drought events over time using estimate marginal means calculated using the “emtrends” function in the “emmeans” R package and used the Tukey HSD method to correct for multiple comparisons (Lenth, 2020). Estimated marginal means are the mean response for each factor, averaged across the other variables in the model. When the adjusted marginal means were extracted from the model, the estimates were back-transformed to the original scale and these are the values presented in Figure 2. Adjusted marginal means and unadjusted 95% confidence intervals were also obtained using “emmeans” for each year and drought event. All calculations were performed using R version 3.6.1 (R Core Team 2019).

We then used a generalized linear mixed model using the “lme4” package to test whether the proportion of trees growing at >2SD below their pre-drought growth average was different during and following the 1976, 1984, and 1995 droughts. Drought event and year were fit as fixed effects along with an interaction between these two variables, while tree ID was fit as a random effect. A 2 degree polynomial was applied to Year. Adjusted marginal means were extracted from the model using the “emmeans” package and the proportions of trees with reduced growth compared at each year. A Bonferroni correction was applied to adjust for multiple comparisons.

To assess the extent to which drought induces shifts in stand structure, we calculated three complementary indices annually for each treatment using all BAI records; G_d, S_i, and S_a. G_d was calculated manually in R following Method 3 outlined in West (2018) by plotting cumulative proportional tree basal area (BA) (χ) against the corresponding cumulative proportional tree BAI (γ)and then fitting multiple polynomial functions using ordinary least squared regression and calculating the area under the fitted curve. We selected the best fitting polynomial function using the lowest AIC value with a maximum polynomial order of eight for each year.

Growth dominance measures the degree to which larger or smaller trees are contributing to overall stand growth relative to their size and ranges from –1 to 1, with a hypothetical value of 0 indicating no size class is dominating growth and all trees are growing at rates directly proportional to their size (West, 2018). The closer G_d is to 1, the more larger sized trees are dominating total stand growth, while the closer G_d is to –1, the more smaller trees are dominating total stand growth (West, 2014, 2018).

Size inequality reflects the variability in tree sizes and ranges from 0 to 1, with inequality being zero (i.e., total equality) if all trees are identical in size. S_a is intimately linked to S_i, with values ranging from 0 to 2. S_a measures how much larger trees (S_a values ranging from 1 to 2) or smaller trees (S_a values ranging from 0 to 1) in the population contribute to the observed levels of S_i. Both S_i and S_a were calculated using the Gini coefficient and the Lorenz asymmetry coefficient, respectively, in the ineq package in R (Zeileis and Kleiber, 2014) by plotting cumulative proportional tree size (BA) against cumulative proportional tree frequency.

To identify whether changes in G_d following drought were indicative of a threshold-type response, we initially tested for structural change in G_d over time (separately for high (ρ_H) and low (ρ_L) density stands) using the efp function in the strucchange package in R (Zeileis et al., 2002) and specifying “OLS-CUSUM” which runs an empirical fluctuation process of OLS residuals to test whether the null hypothesis of no structural change over time is supported. For G_d in both ρ_L and ρ_H the null hypothesis was rejected, indicating the presence of significant structural change at some point(s) in both G_d timeseries (p < 0.05 for G_d in both ρ_H and ρ_L). We subsequently used change point analysis [where the term “change point” is synonymous with the term “threshold” (Andersen et al., 2009)] to identify the number and temporal location of any thresholds for G_d in both ρ_H and ρ_L. This method has been previously used to assess abrupt decreases in tree growth linked to drought (Vanoni et al., 2016a,b) and does not require the number of potential thresholds to be pre-determined. Similarly, this method can detect both positive (abrupt increases in G_d) or negative (abrupt decreases in G_d) change points. The optimum number of change points was simultaneously estimated using the Bayesian Information Criterion (BIC) and the location of each change points along with their 95% CI were estimated as calendar years using the breakpoints function from the strucchange package in R (Zeileis et al., 2002). G_d values for ρ_H were also compared to ρ_L in all three drought years (1976, 1984, and 1995) by permuting the difference between values 10,000 times and using the distribution of these G_d differences to determine the probability of observing a particular value by chance.

During the 1976 drought year, mean PGC was positive in both the high density (ρ_H) (+21%) and low density (ρ_L) (+5%) stands but slightly negative in the 1995 drought year for both ρ_H (−3%) and ρ_L (−1%) (Figure 3). I contrast, mean PGC was −23% in ρ_H and −25% in ρ_L in 1984 (Figure 3) and continued to decline in both stand densities in the 2 years following the 1984 drought, so that by 1986, mean PGC was −43% and −46% in ρ_H and ρ_L respectively, meaning that average tree BAI was 43% and 46% lower than BAI1984_con in 1986. While the impact of all three drought events was more pronounced for spline detrended ring width data, the general patterns in PGC were the same as those derived from the BAI analysis (Supplementary Figure 1). The annual pairwise comparison between droughts showed that the 1984 drought year was associated with significantly more trees growing at >2SD below average than in the 1995 drought year and in all three post-drought years after both the 1976 and 1995 droughts (p < 0.05 in all cases Supplementary Table 4). By 1986, 16 trees in ρ_L (55%) and 13 trees in ρ_H (48%) were growing at rates > 2SD below BAI1984_con levels (Supplementary Table 5). Following 1986, the pattern of continued growth decline reversed and mean PGC became positive in 1988 in both ρ_L and ρ_H (Figure 3) and the number of trees with BAI > 2SD below BAI1976_con, BAI1984_con or BAI1995_con levels was not statistically different between all three drought events (p > 0.9 in all cases) (Supplementary Table 4).

Mixed-effects model analysis showed a significant difference in BAI among drought events, but this difference changed over time (p < 0.001) (Figure 2 and Table 2). Post-hoc analyses of estimated marginal means for BAI were used to compare differences between droughts at an annual resolution (Figure 2). BAI differed significantly after the 1984 drought compared to the 1976 and 1995 droughts for 3 years following the event (p < 0.0001 in all cases, Figure 2 and Supplementary Table 3). However, differences in BAI between the 1976 and 1995 drought were only significant (p < 0.05) in the first post-drought year (p = 0.014), but not in the drought year itself or any other post-drought years (p > 0.06 in all cases) (Figure 2 and Supplementary Table 3). As a result, the slightly more extreme drought of 1984 was associated with a greater absolute change in BAI in the drought year and in the following 3 years relative to the less severe 1976 or 1995 droughts.

Interannual variability in G_d was higher in ρ_H than ρ_L across the study period (Figure 4A). Between 1961 and 1973, G_d in ρ_H indicates that smaller trees were dominating stand growth (small trees were contributing more to total BAI than the same trees were contributing to total BA), but by 1975 larger trees began to dominate growth (larger trees started growing at rates disproportionately fast for their size) (Figure 4A). During the 1976 drought year in ρ_H, all trees were growing at rates roughly proportional to their size (G_d), however, in the years following 1976 larger trees contributed proportionately more to total stand BAI than their BA contributed to total stand BA (Figure 4A). In contrast, for ρ_L, G_d remained relatively stable, oscillating around 0 from 1961 to 1985, indicating that all trees were growing at rates roughly proportional to their size (Figure 4A).

In 1984 G_d shifted away from larger trees toward smaller trees in ρ_H, meaning all trees contributed to total stand BAI at rates roughly proportional to their BA. This directional shift in G_d continued into 1985 as smaller trees showed increased G_d, contributing proportionally more to stand level BAI than they did to total stand BA (Figure 4A). G_d briefly returned to roughly zero (all trees were again growing at rates proportional to their size) in 1986, however, in contrast to values of G_d prior to 1984, smaller trees contributed proportionally more to stand level growth from 1987 to 1991. Larger trees again began to dominate growth from 1992 in ρ_H, a pattern that appeared to temporarily reverse following 1995, where smaller trees dominated growth until 1998. In ρ_L, G_d became consistently negative from 1986 until 1993 when G_d returned to roughly zero, though G_d never became positive again after 1986 in ρ_L (Figure 4A). Despite these shifts in G_d, the ten largest trees for all 3 years prior to the 1984 drought were still the same ten largest trees in 1988.

Across the study period (1961–2002), Si was consistently higher for all years in the higher density (ρ_H) than the lower density (ρ_L) treatment (Figure 4B) meaning that there was a greater range of tree diameters in the higher density stand. There was also a slight decrease in S_i in both ρ_H and ρ_L across the study period but no discernible change in S_i during or after any of the three drought events (Figure 4B).

In general, S_a remained stable throughout the study period in ρ_H, while ρ_L showed a brief depression in S_a between 1978 and 1981, indicating larger trees were temporarily contributing less to S_i than before 1978 (Figure 4C). However, S_a returned to pre-1978 levels in ρ_L by 1982. No immediate change in S_a followed the 1984 drought in either ρ_H or ρ_L, however, between 1990 and 1993, S_a shifted to a higher and stable level in ρ_L, indicating that from 1990 onward in this lower density stand, larger trees were consistently contributing more to size inequality than previously (Figure 4C).

Change point analysis of G_d revealed three change points (i.e., thresholds) for ρ_H (1973, 1983 and 1990) and two change points for ρ_L (1970 and 1985) (Supplementary Figure 2). Since the date assigned to each breakpoint is allocated to the last year of the previous period (i.e., 1 year before a change to a lower or higher G_d level) a change point in 1983 indicates a change occurred in 1984. As a result, the extreme drought event of 1984 aligns with a threshold change in G_d in the same year in ρ_H, while a threshold change in G_d for ρ_L occurred 2 years later in 1986 (Supplementary Figure 2). Neither the severe drought of 1976 or 1995 fell within the 95% CI of any of the other G_d thresholds identified in either ρ_H or ρ_L (Figure 4A and Supplementary Figure 5). Change point analysis was not conducted for S_i or S_a due to apparent insensitivity of these indices around all three drought years (Figure 4) and no significant differences in G_d were noted between ρ_H and ρ_L in either 1976, 1984 or 1995 (p > 0.613 in all cases).

Many studies calculate resistance, recovery and resilience relative to a pre-drought growth average derived from the years immediately preceding a drought event (Zang et al., 2014; Gazol et al., 2017, 2018). However, recent work has shown that the subjective choice of the number of years to include when calculating these pre-drought growth averages can bias estimates of drought resilience (Schwarz et al., 2020; Ovenden et al., 2021). Here we developed an approach which uses the SPEI to climatically standardize the selection of these pre-drought years to avoid the need for this subjective selection and minimize the risk of including pre-drought year BAI values that may themselves be the product of an abnormally dry year in the calculation of pre-drought average growth.

Using these climatically constrained pre-drought growth averages to address our first objective, we show that the extreme drought of 1984 was associated with a protracted period of low growth and changes in growth dominance at the stand level, consistent with a threshold-type drought response noted in other pine species (Huang et al., 2015; Asbjornsen et al., 2021). However, we did not find any evidence of a similar response following the lower severity droughts on 1976 and 1995, nor did we document any significant differences between the high or low density stands (Table 2), the latter result being in keeping with other recent work on the growth response of P. sylvestris to drought which also documented no effect of tree density in this species (Bello et al., 2019). As such, we found that a small increase in drought severity [measured in both absolute (CWD) and relative (SPEI) terms] in 1984 compared to the 1976 and 1995 droughts was associated with a significant decrease in average tree BAI and a significant increase in the number of trees growing a > 2SD below their pre-drought average for up to three post-drought years (Figure 2 and Supplementary Tables 3, 4). These results reflect other studies that have documented multi-year legacies on tree growth following extreme drought (Anderegg et al., 2015b; Wu et al., 2018). Huang et al. (2015) reported an SPEI threshold of −1.64 for Pinus edulis and Pinus ponderosa, after which progressively more negative SPEI values cause significant declines in forest growth. This value reported by Huang et al. (2015) is slightly less negative than the SPEI recorded for all three droughts in this present study (1976 SPEI_Aug6 of – 1.8, 1984 SPEI_Aug6 of –2.0 and a 1995 SPEI_Aug6 of – 1.9) (Figure 1). As such, our results support other recent work that suggests that thresholds of drought severity likely vary among species (Cavin et al., 2013; Kolb, 2015; Gazol et al., 2020) but also highlights how setting generic, climatologically defined drought thresholds too low can risk including events that may not have been strong enough to elicit an organismal or systemic response (Smith, 2011), potentially clouding our understanding of how different species respond to extreme events. Similarly, the characterization of drought purely in terms of severity using indices (e.g., SPEI or CWD) may also obscure important differences between droughts stemming from the inherent multi-dimensionality of these events (e.g., timing, duration, and intensity). Due to the limited number of drought events in the climate data at this study site, it was not possible to ascertain whether the timing, duration or intensity of drought at this site were linked to patterns in tree growth response. While none of the post-drought years in this study appeared to be abnormally dry (Figure 1), we acknowledge that post-drought climate likely plays an important role in regulating forest recovery dynamics but was not explicitly considered in this study, however, where possible future studies on forest resilience to drought should aim to included post-drought climate to ensure a more complete understanding of drought legacy.

At the stand level, mean PGC in both the high and low density stands remained relatively constant during and following both the 1976 and 1995 droughts with some trees showing reductions in radial growth while others showed radial growth increases. The net effect of this inter-tree variability may have partly buffered against stand-level changes in BAI during these two less severe droughts. While individual tree growth variability is negatively linked to individual drought resilience (Bose et al., 2020), inter-tree variability in pre-drought growth may act as a form of response diversity (Mori et al., 2013) and reflects other recent studies that highlight the potential for intraspecific differences to partly buffer against the impacts of climate change (Oney et al., 2013; Taeger et al., 2015; Moran et al., 2016). In contrast, after the 1984 drought stand growth reduced to levels substantially below BAI1984_con (Figure 3), temporarily reduced inter-tree growth variability (Figure 3) and shifted growth dominance (Gd) toward smaller trees (Figure 4), indicating that a threshold for intraspecific variability to buffer against stand level growth loss may have been exceeded.

The drought severity required to cross such thresholds is likely to be linked to the ecophysiological limits of a species hydraulic system to drought stress (Choat et al., 2012; Adams et al., 2017) and pre-drought conditions, which can lead to phenomenon such as structural overshoot (Jump et al., 2017). While it is possible that observed shifts in Gd in this study are partly the result of the death of neighboring trees, assessments of both the high density (ρ_H) and low density (ρ_L) treatments in 1985 and 1990 (1 and 6 years after the 1984 extreme drought) show mortality was not abnormally high during these post-drought periods (Supplementary Table 1). While we acknowledge that in some cases mortality might be delayed for many years following drought (Bigler et al., 2007), we show that the crossing of drought thresholds does not necessarily need to result in widespread or elevated tree mortality to be associated with detectable changes in growth dominance and stand dynamics (Stuart-Haëntjens et al., 2015; Batllori et al., 2020), addressing our second research question. Persistent, interspecific shifts in competitive dominance have been documented following extreme drought in other species (Cavin et al., 2013), but such shifts at the population level were not observed in a much larger scale study of P. sylvestris (Martínez-Vilalta et al., 2012).

Collectively, our results demonstrate how small increases in drought severity can be associated with changes in tree growth at different tree densities and may influence stand dynamics to a degree that is not observed under less severe drought conditions. While we did not detect any large changes in Si or Sa in response to drought, this study exclusively looked at even aged, monospecific stands of P. sylvestris. Future work should investigate whether drought induces changes in Si, Sa, and Gd in more structurally diverse, species rich stands of variable age classes to further our understanding of the role of drought in driving novel developmental trajectories and structuring community composition. Similarly, understanding which combinations of drought intensity, timing, duration, and frequency are associated with differences in forest growth and stand development will be key.