Delayed and altered post-fire recovery pathways of Mediterranean shrubland under 20-year drought manipulation

drought ( (cid:0) 30% rainfall) on the pathways of vegetation metrics related to species richness, community composition and abundance dynamics for an early-successional Mediterranean shrubland. The results indicate that the pathways of vegetation metrics were differently affected by experimental drought. The abundance of Globularia alypum follows pathway 1 (altered mature state). Simpson diversity and abundance of Erica multiflora follow pathway 2 (delayed succession) while species richness, community abundance and shrub abundance follow pathway 3 (alternative stable state). There were no significances for the resilience to extremely dry years (the ratio between the performance after and before severe events) between control and drought treatment for all vegetation metric. But, their resilience for the metrics (except Simpson diversity) to extremely dry years in 2016 – 17 were significantly lower than that of 2001 and of 2006 – 07, possibly caused by the severe water deficits in 2016 – 17 at mature successional stage. Principal component analysis (PCA) shows that the first two principal components explained 72.3 % of the variance in vegetation metrics. The first axis was mainly related to the changes in community abundance, shrub abundance and species richness while the second axis was related to Simpson diversity and abundance of G. alypum and E. multiflora . Principal component scores along PC1 between control and drought treatment were significantly decreased by long-term experimental drought, but the scores along PC2 were not affected. Further research should focus on successional pathways in more water-deficit conditions in Mediterranean ecosystems and the consequences of changes in vegetation recovery pathways on ecosystem functions such as biomass accumulation and soil properties.


Introduction
Anthropogenic disturbances (e.g. land-use change, climate and frequent fires) present an unprecedented risk for natural ecosystems, having significantly altered 75% of terrestrial ecosystems, leading to abrupt biodiversity loss and ecosystem function decline (Moreno-Mateos et al., 2017;Díaz et al., 2019). Considerable efforts have been dedicated to understanding how species diversity and composition, and ecosystem functions recover over time (Anderson-Teixeira et al., 2013). However, increasing water deficits characterized by lower soil moisture and higher atmospheric aridity are posing greater challenges to the disturbancerecovery due to the resultant decreases in new establishment and increases in plant mortality (Angeler and Allen, 2016;Schwalm et al., 2017;Stevens-Rumann et al., 2018;Astigarraga et al., 2020). Increased frequency and intensity of extreme episodes (e.g. heatwaves and severe droughts) may trigger greater limitation on plant growth and survival, leading to the abrupt losses in local species diversity, composition shifts and further declines in the carbon sink for terrestrial forest ecosystems (Ciais et al., 2005;Nolan et al., 2018;Trisos et al., 2020). Therefore, accurately understanding the dynamics of ecological recovery in waterstressed ecosystems is of great importance to conserve the local diversity and maintain ecosystem functions in a rapidly changing climate for coming decades.
There are substantial knowledge gaps on the dynamics of ecological recovery in response to water deficit. First, which vegetation metrics could be affected that is not well-addressed. Vegetation recovery is mainly focused on the ecosystem functions such as vegetation biomass and productivity (Anderson-Teixeira et al., 2013). Multiple vegetation metrics such as species richness, community composition and abundance (or density) are important indicators of vegetation recovery dynamics that are less studied (Anderson-Teixeira et al., 2013). Second, the trajectories of vegetation recovery over time in a changing climate are not clear. Increasing summer water deficit can limit the species occurrences of reseeding (regenerating from seed) rather than resprouting species, potentially shifting community composition over time (Parra and Moreno, 2018; van Blerk et al., 2021a). However, the rates of vegetation recovery may be strengthened if some drought-adapted species become the dominant and reach a different mature state (Anderson-Teixeira et al., 2013). In addition, severe droughts may delay or shift vegetation post-recovery via the alteration in growth resilience (the capacity to reach pre-disturbance performance) (Parra and Moreno, 2018;Lloret et al., 2011Lloret et al., , 2012. Third, effective approaches to investigate ecological recovery at stand level is greatly limited. Ecosystem modelling has been applied to estimate the recovery time of global forests after disturbances such as biomass harvest, large fires and insect outbreaks (Pugh et al., 2019), but the effect of increasing drought stress on recovery was not investigated. Interestingly, Baudena et al. (2020) have predicted that increasing aridity may reduce the post-fire recovery (e.g. plant cover) of Mediterranean forests based on a modelling approach, with large areas shifting towards open shrublands. However, it is difficult to generalise modelling approaches on the dynamics of vegetation metrics for local communities, since the recovery processes greatly depend on the vegetation types, soil characteristics and biotic competition (Capitanio and Carcaillet, 2008;Anderson-Teixeira et al., 2013;van Blerk et al., 2021b). Therefore, finding effective ways to capture the recovery dynamics of multiple attributes of community in response to climate change is essential.
Manipulative field experiments have the potential to verify the recovery dynamics of species diversity, community composition and ecological functions and biogeochemical processes by conducting a continuous climate change in natural ecosystems (Wu et al., 2011;Beier et al., 2012;Peñuelas et al., 2018). In particular, long-term (decadal) experiments can offer unique possibilities to detect the lagged and accumulated impacts, which include plant physiological, demographical and evolutionary adjustments (Harte et al., 2015;Zhang et al., 2017;Liu et al., 2020) and biotic interactions (Wright et al., 2017;Liu et al., 2018b). However, most of these climatic experiments have been performed in mature or equilibrium ecosystems where the structure and function of plant communities are relatively complete (Beier et al., 2012;Kröel-Dulay et al., 2015;Liu et al., 2020). Thus, vegetation responses, e. g. species richness, community composition and productivity were relatively resistant to climate treatments to the point that only minor and even no net changes were detected for arctic tundra (Hudson and Henry, 2010), temperate grassland (Grime et al., 2008) and Mediterranean shrubland (Tielbörger et al., 2014). In contrast, considerable vegetation changes in response to climate manipulation were detected for the recovering communities (or early-successional status) since the species assembles were unstable Liu et al., 2017;Parra and Moreno 2018;van Blerk et al., 2021a). It fits well with hierarchical theory that abrupt and nonlinear changes in community structure would likely occur over the long term by species reordering and competition outcomes (Smith et al., 2009;Luo et al., 2011). In addition, the impact of extreme natural droughts on community and population dynamics could be further intensified, due to changes in the resilience in the climate-manipulated ecosystems Liu et al., 2015Liu et al., , 2018aParra and Moreno 2018). But it has also been Fig. 1. Conceptual diagram hypothesising four successional pathways in which climate change may impact the dynamics of vegetation recovery. The blue line represents the value of a vegetation metric following fire disturbance in the absence of experimental drought (control). Metrics include species richness, community composition and abundance dynamics (for community, shrub and species levels). The yellow lines indicate the four pathways simulated by experimental drought (− 30% rainfall), representing a possible impact of climate change (Fig. 1a). The vertical dashed line separated the vegetation recovery into two stages (young and mature) in control, the peak year was selected based on the highest value of vegetation metrics. The right figures (Fig. 1b,c,d,e) show the percent changes of vegetation metrics (calculated as: 100*(drought -control)/control), with the horizontal lines indicate the values at the beginning of experiment. Therefore, long-term drought experiment and extreme droughts could lead to four different succession pathways for the early successional ecosystem. 1) Altered mature state (pathway 1): the species can obtain resource benefits under drought conditions and perform better than the control (Fig. 1b). 2) Delayed succession (pathway 2): the new establishment and growth are transiently limited, and the recovery of vegetation metrics may take a longer time to reach its expected mature state (Fig. 1c). 3) Alternative stable state (pathway 3): some species are permanently filtered out due to species reordering, plant community will achieve an alternative stable status (Fig. 1d). 4) Advanced degradation (pathway 4): species establishment and growth are strongly limited and tend to degrade over time due to the accumulated and legacy effects of drought (Fig. 1e). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) reported that vegetation may have the resilience after extreme droughts via the stabilisation process of increasing new recruitment and regrowth (Seifan et al., 2010;Lloret et al., 2012;Ploughe et al., 2019). To date, few long-term drought experiments have been performed in recovering and unstable habitats; detailed knowledge on the succession pathways of species diversity, composition and abundance dynamics, is highly limited.
The successional recovery is strongly affected by drought in the water-limited Mediterranean-type ecosystems. The Mediterranean Basin is one of the most highly disturbed regions, due to a combination of anthropogenic and natural activities, and it has been demonstrated that these disturbed (or successional) communities are sensitive to ongoing and future climate change (especially drought) (Myers et al., 2000;Doblas-Miranda et al., 2015;Cramer et al., 2018). Negative impacts on ecological recovery (e.g. forest degradation) are expected to be further exacerbated by the projected 2-4 • C summer temperature increase and continuous 30% summer precipitation decrease by the end of this century (Cramer et al., 2018;Lionello and Scarascia, 2018). Experimental results have indicated that drought has already delayed the post-fire successional process by decreasing the species richness and density of recruited seedlings (Lloret et al., 2004(Lloret et al., , 2009del Cacho et al., 2013). Once the communities reached a more stable status, increases in species richness under drought can happen, due to the establishment of droughtadapted species (e.g. G. alypum) (Liu et al., 2018b). But abundance changes may vary greatly with functional groups (e.g. herbs vs shrubs). For instance, in an early-successional Mediterranean shrubland ecosystems, both long-term warming and drought had no significant influence on the abundance of shrubs . Indeed, the shrubs and the dominant species (G. alypum and E. multiflora) may acquire resource benefits due to their long life-history, physiological adaptions (e.g. increased photosynthesis and water-use efficiency) (Liu et al., 2016;Zhang et al., 2017) and delayed leaf-unfolding and flowering (Peñuelas et al., 2004;Prieto et al., 2008). However, the successional pathways of species diversity, community composition and abundance dynamics affected by long-term drought, are still unclear.
We conducted a drought (− 30% rainfall) manipulation in an earlysuccessional Mediterranean shrubland to study the recovery dynamics over the period of 1999-2018. The study site had suffered a wildfire 4 years before the drought manipulation. The amount of rainfall reduction in our experiment fell into the range of rainfall changes in the projection of climatic models (up to − 30% with the scenario of 2-4 • C) for the Mediterranean Basin during the coming decades (Cramer et al., 2018;Lionello and Scarascia, 2018). We hypothesized that long-term (20year) climatic manipulation, combined with severe droughts (a 10th percentile threshold of annual precipitation), would shift the pathways of vegetation recovery and reduce the vegetation resilience to severe events in this Mediterranean early-successional shrubland. The objectives were to 1) study the impacts of long-term drought on species diversity, community composition, abundance dynamics (for community, shrub and two dominant species; respectively) and to identify the pathways (altered mature state, delayed succession, alternative stable state and advanced degradation) explained in Fig. 1; 2) compare the likely differences in the resilience to severe droughts between drought treatment and control; 3) analyze the effect of experimental drought on the changes in spatial variation of the vegetation metrics, which can be investigated by the change in the main principal components.

Study site
The study site is located in the Garraf Natural Park south of Barcelona, northeastern Spain (41 • 18Ń, 1 • 49É; 210 m a.m.s.l.). The climate is typically Mediterranean-type, characterised with dry and hot summers (June to August) and wet springs (March to May) and autumns (September to November). Mean annual precipitation (MAP) during the period of 1951-2018 was 571.2 mm, with nearly 62 % occurring in spring and autumn (354.1 mm). Mean annual temperature (MAT) during the same period was 15.3 • C, and the hottest season occurred in summer (mean temperature is 22.8 • C) (Fig. S1). The soil type is calcareous and composed of marlstone and limestone. The soil depth above the bedrock is about 10-40 cm. The study site has suffered two severe wildfires, one in summer 1982 and the other in spring 1994. The conifer forest (dominated by Pinus halepensis) was damaged by the two fires. Thus, most of the current vegetation, mainly Mediterranean shrubs, has sprouted from the underground organs after the fires. The vegetation in this site has the potential to recover to become a secondary forest after the two fires. But tree seedlings (e.g. Pinus halepensis) were cleared irregularly to maintain the vegetation as shrubland. The height of shrubland vegetation is<1.5 m. The canopy vegetation cover is about 40-60%. There are nearly 30 Mediterranean shrub and herbaceous species in this study site. The dominant shrubs are Globularia alypum (GA) and Erica multiflora (EM). The abundance of GA and EM relative to total abundance is larger than 50% in our study site. Moreover, there are other shrub species, such as Rhthtosmarinus officinalis (RO), Ulex parviflorus (UP) and Dorycnium pentaphyllum (DP). The undergrowth also includes some herbaceous species (mostly Poaceae species (Poaceae spp.)).

Field manipulations
The treatments of control and drought experiments were initiated in 1999 (1998 was a pre-treatment year). Three blocks with two plots (4 m × 5 m) each were selected along a south-facing slope (ca. 13%). Treatments were randomly distributed in each block (Beier et al., 2004;Liu et al., 2017Liu et al., , 2020. Scaffolding (1.2 m height) was constructed in all of the plots. But there were no curtains or covers in control plots. Future drought condition was simulated by intercepting and removing rainfall with transparent waterproof covers during the wet seasons (spring and autumn). The drought curtains were activated by a rain sensor and covered the drought plots during rainfall events larger than 5 mm and were retracted once the rain ended. The automatic system for drought was also damaged by lighting in 2015 and then was redesigned to passive continuous rain interception using transparent V-type plastic over the scaffolds designed to intercept nearly 30% of the annual rainfall ( Fig. S2), the similar percentage that was intercepted with the automatic system. The cover system had minor influence on the soil and surface soil temperature. Soil moisture at 0-15 cm depth was measured biweekly using Time Domain Reflectometry (TDR, Tektronix 1502C; Tektronix, Beaverton, OR, USA) in each plot. Throughout the study period of 1999-2018, the 30% rainfall reduction in drought plots reduced the soil moisture at 0-15 cm depth by 16.7%, compared to control plots (15.50 ± 0.7% and 18.59 ± 0.6% for drought and control respectively, n = 20).

Vegetation measurements
Vegetation was measured by a nondestructive pin-pointing method conducted annually during the non-growing season in late July and the beginning of August from 1998 (pre-treatment) to 2018. With the pinpoint method, we successfully recorded every hit within the vegetation profile that allowed us to estimate total number and density frequencies for every species and that were used to follow the changes over time in the plant community structure and composition. We conducted vegetation measurements since the pre-treatment year in five parallel 3m transects (0.8 m apart) marked permanently in the central area of each plot (3 m × 4 m). Each transect had a total of 61 sampling points that were separated every 5 cm, totaling 305 points per plot. For the sampling, we lowered a 2 m long steel pin (3 mm diameter) through each sampling point and recorded, for each plant hit, the organ (e.g. leaf, branch and reproductive of flower and fruit) and its status (alive or dead). Our measurements were performed during the non-growing season, the dead hits included both those resulting from naturally drought and those caused by the climate manipulation. Since the dead hits occurred in both control and drought plots.

Vegetation metrics and climate variables
Species richness and abundance dynamics were measured yearly in each plot during the study period of 1998-2018. The abundance (hits) of all species in the plot was recorded and used to calculate the index of abundance dynamics and Simpson diversity. Here abundance metric was calculated as the number of alive hits or frequency. The number of alive hits for each species were naturally log-transformed and then summed to get the equal weight of community abundance. The shrub abundance each plot was obtained by summing all the log-transformed hits of shrub species according to the growth-form identified in our study site. The abundance of the two dominant species G. alypum and E. multiflora was analyzed separately. Community composition was estimated by the Simpson diversity index for each plot throughout the study period. It was calculated as: where n is the number of alive hits in a specific species and N is total number of alive hits for all species. Simpson diversity takes into account the number of plant species and the abundance of each species in a community. The values of its range from 0 to 1, with higher values indicating higher diversity. Simpson diversity were calculated with the vegdist function in the vegan 2.5-6 package (Oksanen et al., 2017) in R version 3.3.0.
Climate variables of monthly mean temperature and total precipitation for 1951-2018 were compiled from a nearby meteorological station during the period of 1951-2010 and our study site during the study period of 1998-2018 (linear regression for the two sites for 1998-2010 of R 2 = 0.99 for monthly temperature and R 2 = 0.83 for monthly precipitation). Based on this profile, we calculated the historical climate variables of annual mean temperature and total precipitation for the study site. The natural extremely dry years were obtained according to the lowest 10th percentile precipitation year during the period of 1951-2018 (Knapp et al., 2015). During the experimental period of 1998-2018, extremely dry years occurred in 2001, 2006-07 and 2016-17 (Fig. S1).

Statistical analyses
We selected generalized additive mixed models (GAMMs) to analyze the dynamics of vegetation metrics with a smoothing function using the mgcv package (version 1.8.3) (Wood, 2017), because vegetation metrics of post-fire recovery process are nonlinear and flexible. The bam function was implemented in the GAMM analysis. The smoothing functions could be used for testing potentially different trends over time for different treatments (control, drought). The model formula was performed with ordered factors that included intercept difference and reference smoothing: vegetation metric ~ treatment + s(year) + s(year, by = treatments). Here, treatment that had two levels, drought and control, was considered as a reordered factor, and year as a numeric variable from 1998 (pre-treatment year) to 2018. The plot nested within block was treated as the random factor in the analyses. In order to explain the temporal autocorrelation in the time series, we specified the  Table S1).
We conducted the assessments on the resilience to extremely dry years of 2001, 2006-07 and 2016-17. Resilience is the capacity to reach pre-disturbance performance, which was calculated as the values after an extremely dry year versus those in the pre-drought year in each plot (PostDry / PreDry) (Lloret et al., 2011) (2010). Two-way ANOVA (treatment and extreme dry years) was used to test the significant influences. Once the significance was obtained, we further analysed the ANOVA with Tukey's HSD (honest significant difference) post-hoc test.
Principal component analysis (PCA) was conducted to explore and visualize the relationships between the multiple quantitative vegetation metrics and to explore which metrics are the most crucial in distinguishing the functional differences for treatments. The vegetation metrics of species richness, Simpson diversity, community abundance, shrub abundance, abundance of G. alypum (GA) and E. multiflora (EM) were selected. We applied Horn's parallel analysis (paran package of version 1.5.2) to determine how many axes of vegetation metrics should be retained in principal component analysis. The result of the parallel analysis indicated that all of the 6 metrics should be retained in the PCA (Fig. S3). The treatment (control and drought) was considered as a categorical variable. The relationship between the vegetation metric and study years were analyzed using the Performance Analytics package (version 2.0.4) in R (Peterson et al., 2019). Here, we used the packages FactoMineR (version 1.34) for the analysis (Husson et al., 2013) and factoextra (version 1.0.7) for ggplot2-based visualization (Kassambara, 2016). The scores of the vegetation variables in each plot during the period 1999-2018 were obtained according to the score function. Here we only selected the scores along the first two dimensions. Due to the non-independence of the measurements in control and drought over time, we included an autocorrelation structure in the model using nlme package. The significant differences for the interaction of treatments (control and drought) and year were analysed. Plot was selected as the random factor in the best model. Repeated measure for a two-way analysis of variance (ANOVA) was applied to test the significance. The mean score of each treatment × year combination was calculated by groupwiseMean function from the rcompanion package, and the confidence interval for the mean scores which were presented in the Fig. 6b, c.
All the analyses were performed in R (version 4.0.5).

1 Species richness and Simpson diversity in response to long-term drought
Species richness and Simpson diversity were significantly decreased by drought (both p < 0.001) during the period 1999-2018 based on the GAMM test ( Fig. 2a and Fig. 2b; respectively) ( Table S1). The smoothing in control for species richness and Simpson diversity were nonlinear, and they were significantly different from 0 in the model analysis. We do not observe a significant interaction of treatment and year for these two vegetation metrics. Based on the result of model analysis and percentage difference between control and drought, we can conclude that species richness and Simpson diversity follow pathway 3 (alternative stable state) and pathway 2 (delayed succession), respectively.

Abundance dynamics affected by drought treatment
Both community abundance and shrub abundance were significantly decreased by drought treatment over the study period (Fig. 3a,b; Table S1). The GAMM tests suggest that the smoothing was significantly increasing over the study period (both p < 0.001). The changes in abundance at community and shrub levels by drought follow the pathway 3 (alternative stable state) (Fig. 3c, d).
Drought treatment had different effects on the abundance for the two dominant species (Fig. 4a, b; Table S1). Drought treatment significantly increased abundance of G. alypum (p < 0.001) and there was an interactive effect of drought and year (p < 0.05). Yet it did not affect abundance of E. multiflora during the study period. The changes in abundance of G. alypum and of E. multiflora support pathway 1 (altered mature state) and pathway 2 (delayed succession); respectively.

The resilience to extreme droughts
We observed that there were no significant differences for the resilience of any of the vegetation metrics between control and drought treatments based on the ANOVA test ( Fig. 5; Table S2). However, there were significant differences among the three extreme events for the resilience of vegetation metrics (except Simpson diversity). Moreover, the values of resilience were similar between the first two extreme events (2001 and 2006-07), which occurred during the early successional stage, and the values at the mature successional stage  were lower than the values at the early stage.

The variation in vegetation metric space over time
The first two principal components (PC; axes) of the PCA explained 72.3 % of the variance in vegetation metrics, PC1 (47.4 %) and PC2 (24.9 %) (Fig. 6a). The first dimension (PC1) was related to the changes in community abundance, shrub abundance and species richness, that were highly positively associated in the first dimension (contributed by 56.3 %, 52.9 % and 48.1 %, respectively). However, the second dimension (PC2) was linked to the Simpson diversity (contributed by 47.3 %) and abundance of G. alypum and E. multiflora (contributed by 74.2 % and a 37.2 %, respectively). Simpson diversity was negatively correlated with the changes in the abundance of G. alypum and E. multiflora (Fig. S4). In addition, the scores for the two PCs varied differently with the treatments over time (Fig. 6b). PC1 scores were significantly decreased by experimental drought (anova, p < 0.001). PC2 scores were not affected by treatment, but the changes varied significantly over time (p < 0.001).

Changes in the species richness and community diversity
Long-term drought had significant impacts on species richness and composition. Species richness and Simpson diversity were significantly decreased by drought treatment but the smoothing differences between control and drought did not change over time. The changes in species richness under drought fit well with the pathway 3 (altered maturity state) that reach a relatively stable state over time (Fig. 2a), possibly due to the loss of drought-sensitive species (e.g. Galium lucidum, Ulex parviflorus and Lithodora fruticosa). Simpson diversity in drought appears to follow pathway 2 (delayed succession) that took a longer time to reach its expected mature state. This can be attributed to increases in abundance of drought-resistant species (e.g. G. alypum, Polygala rupestris and Fumana thymifolia) in dry condition (Liu et al., 2018b). However, we observe the resilience for species richness to the extreme dry years of 2016-17 (mature stage) was lower than to those of 2001(young stage). This result suggests more negative influences by extreme events during the mature stage of succession, although this response could also be related to the cumulative effect of multiple droughts on the recovering vegetation. This contrasts with the response reported for mature Mediterranean shrubland in Israel where no net changes in species richness and biomass accumulation were reported under decadal drought (− 30% rainfall), partly due to high resistance of mature communities (Tielbörger et al., 2014). Moreover, it was reported that few sites had significant changes in shrubland species richness by experimental drought over time for study sites across the climate gradients of Europe (Kröel-Dulay et al., 2015).

Experimental drought affected abundance dynamics over time
Community abundance and shrub abundance were significantly decreased by experimental drought. These two metrics under drought fit with the pathway 3 (alternative stable state) that being restricted to recover as the mature state in control plots and reach a relatively stable state. Less abundance (density) for both community and shrub levels under drought might be associated with the limitation on occurrence of new organs (e.g. leaves, branch and flowers) Liu et al., 2017). However, the significant decrease in both community abundance and shrub abundance over time, especially in the extremely dry years of 2016-17. One reason may be related to self-thinning once the plant communities were approaching mature status (Yu et al., 2019). The resilience of the community abundance and shrub abundance to extreme events were decreased over time; the values were lower in mature stage (in 2016-17) than in young stage (in 2001 and 2006-07) (Fig. 5c, d). Another reason that the rainfall of these two years was relatively low (average 402 mm), which was much less than the average rainfall (586 mm across 1951-2018).
The drought impacts on abundance for the dominant species differed greatly. Drought significantly increased the abundance of G. alypum compared to control, and this drought difference was increased over time. But experimental drought decreased abundance of E. multiflora in the first few years and then had almost no influences after 2007. Therefore, the abundance of G. alypum follows pathway 1 (alternative stable state), but E. multiflora follows the pathway 2 (delayed succession). Previous studies have reported that the dominant species (especially G. alypum) have advantages to endure warmer and drier conditions by adjusting their photosynthesis and water-use efficiency (Liu et al., 2016;Zhang et al., 2017) and phenological alteration (e.g. delaying or advancing the leaf appearance) (Peñuelas et al., 2004;Prieto et al., 2008). However, the abundance of G. alypum abruptly declined in both control and drought plots in extremely dry years of 2016-17, which may be caused by accumulatively water deficits. We also observed a lower resilience in the period of 2016-16 than that of extremely dry years of 2001 and 2006-07 (Fig. 5g, h). Thus, more severe extreme drought in future may decrease the growth resilience of these two species. Thus, the two dominant shrub species may lose competitive advantages in hot and/or dry conditions if other more drought-adapted species (e.g. Polygala rupestris and Fumana thymifolia).

Change in vegetation metric space over time
The first two principal components clearly reflect the variation of the vegetation metric space over time. In the analysis, the first axis was largely driven by the changes in community abundance, shrub abundance and species richness (Fig. 6a). High positive correlations were observed for these three metrics in the positive direction of first dimension (Fig. S4). The second axis was mainly driven by the Simpson diversity and abundance of G. alypum and E. multiflora, and Simpson diversity was negatively correlated with the changes in abundance of two dominant species (Fig. 6a). However, the scores along PC1 were significantly higher in control than drought, indicating the vegetation space was shifted by experimental drought over time. The assessment of variation in vegetation metrics could be useful to understand the metric space associated with plant ecological strategies (Díaz et al., 2016;Bruelheide et al., 2018).

Greater structure and composition shift for future climate scenarios
The changes of vegetation metrics in drought condition varied greatly over time, it is still challenging to definitively attribute to pathway types. Such attribution should become easier as the timeseries is extended, highlighting the importance of conducting multi-decadal manipulations to understand the real types of succession pathway over time. For pathway 4 (advanced degradation) to occur, the drought resistance and plasticity are likely overwhelmed (Tielbörger et al., 2014), which may take many years or decades to become apparent. Thus, it is plausible that given a longer drought perturbation or moderately more severe extremes, these systems may yet shift towards pathway 4, with more severe consequences on ecosystem functions in future. Thus, continued long-term drought manipulation studies are essential to provide an early warning of how climate change affects species diversity, community composition and abundance dynamics for recovering communities.
These impacts on recovery also need to be put in the context of changing disturbance regimes. Global warming and intensified water stress may increase the occurrence and severity of fires in the Mediterranean ecosystems (Díaz-delgado et al., 2002;Doblas-Miranda et al., 2015;Cramer et al., 2018). This would shift larger areas of the Mediterranean Basin into an earlier successional stage, exposing them to the altered successional processes that we have described herein, and thus accelerating a shift in species diversity and community composition across the region. Changes in disturbance types and intensity may also be important. Whilst Mediterranean-type ecosystems typically have a fast recovery after fire disturbance (Pausas et al., 2008;Kelly et al., 2020), many studies have indicated that high burn severity appears likely to decrease the ability of Mediterranean forest to recover (Fernandez-Manso et al., 2016). Increases in fire severity may interact with the selection pressure of drought to further alter Mediterranean ecosystems. Moreover, habitat conditions such as soil water availability, nutrient content and soil seed bank as well as occurrence of alien species may also affect the post-disturbance of local species diversity and community composition (Doblas-Miranda et al., 2015;van Blerk et al., 2021a,b;Keeley and Brennan, 2012). Assessing the combined effects of these multiple dimensions of change is a major challenge, but our results highlight the important role of drought stress during the recovery process in leading to long-term alterations in species diversity, composition and abundance dynamics. Through their destabilising role within ecosystems, disturbances bring plant communities to a state in which they are most sensitive to climate change. The changes for PC1 score were analysed by the interaction of treatment (control and drought) and year using the gls function in nlme package. Plot was selected as the random factor in the best model. The mean score of each treatment × year combination was calculated by groupwiseMean function from rcompanion package. The blue lines show the control, while orange lines show drought. The error bars represented the confidence interval of the plots. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

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.