<figure><img src="image_1.png" /></figure>

HAL
open science

# Uncovering global drivers threatening vegetation resilience

Camille Fournier de Lauriere, Katharina Runge, Gabriel Smith, Vasilis Dakos, Sonia Kefi, Thomas Crowther, Miguel Berdugo

<figure><img src="image_2.png" /></figure>

## To cite this version:

Camille Fournier de Lauriere, Katharina Runge, Gabriel Smith, Vasilis Dakos, Sonia Kefi, et al. Uncovering global drivers threatening vegetation resilience. 2023. [hal-04311563](https://hal.umontpellier.fr/hal-04311563)

HAL Id: hal-04311563

https://hal.umontpellier.fr/hal-04311563

Preprint submitted on 28 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Uncovering global drivers threatening vegetation resilience

Camille Fournier de Lauriere*†¹, Katharina Runge†¹, Gabriel Smith¹, Vasilis Dakos², Sonia Kéfi²,
Thomas Crowther¹, and Miguel Berdugo³

¹Institute of Integrative Biology, ETH Zurich (Swiss Federal Institute of Technology), Zurich, 8092 Switzerland

²ISEM, CNRS, University of Montpellier, IRD, EPHE, Montpellier, France

³Institut de Biologia Evolutiva (CSIC-UPF), Barcelona 08003, Spain

†CFL and KR contributed equally to this work.

*Correspondence should be addressed to CFL at camille.fournierdel@gmail.com
Abstract:

1) Context: The resilience of the Earth's vegetation is changing heterogeneously, making it a challenge to unveil what causes these resilience changes. Understanding the driving forces of these changes can help us make informed management decisions to protect and restore ecosystems. Here, we address this gap by identifying the drivers that have caused the resilience of ecosystems to change during the last two decades.

2) Methods: We globally measured two complementary aspects of resilience, namely sensitivity and autocorrelation, which are respectively associated with the resisting and recovering ability of ecosystems. We used a machine learning approach to identify the main environmental, climatic, and anthropogenic drivers of changes in resilience between two periods (the period 2000-2010 vs. that of 2010-2020).

3) Results: We found that in 26% of regions worldwide, vegetation exhibits signs of resilience loss. Moreover, ecosystem's properties (aridity, elevation, anthropization) affect the way vegetation resilience has changed over time. When controlling for these properties, different biomes (forest, grasslands, and savannas) will exhibit similar responses to changes in climate conditions. Regions experiencing intense warming (>0.2°C/decade) have shown a major loss in vegetation resilience. Decreasing productivity is associated with reduced resilience and interacts with warming, exacerbating resilience loss of less productive lands (potentially showing signs of degradation).

4) Conclusions: Warming and degradation appear as major drivers of losses in vegetation resilience across vegetation types. These results raise concerns about the persistence of ecosystems under continued climate change and expected intensification of human activities which, highlights the importance of maintaining the resilience of ecosystems under changing environmental conditions.

Keywords: resistance, sensitivity, autocorrelation, recovery, climate change, degradation, aridity.
## Introduction

The resilience of ecosystems underpins their ability to absorb environmental fluctuations, as well as their ability to recover from disturbances. Therefore, resilience is critical for ensuring constant provision of key ecosystem services for humanity such as productivity or climate regulation. Alarmingly, climate and land use changes are threatening the ability of ecosystems to withstand disturbance across the globe (Johnstone et al., 2016; McDowell et al., 2020). Quantifying those changes in ecological resilience across the globe remains a major conceptual challenge in ecology (Dakos & Kéfi, 2022). Up to now, a wealth of studies, using a diversity of metrics measuring resilience, reveal increases and decreases in ecosystem resilience across the globe, showing very heterogeneous patterns (Boulton et al., 2022; Feng et al., 2021; Forzieri et al., 2022; Smith et al., 2022; Y. Zhang et al., 2022). However, there is little known about the mechanisms that drive these patterns (Smith et al., 2022). Shedding light on what creates these patterns is critical if we are to develop a unified, mechanistic understanding of the threats to ecological systems across the globe (Smith et al., 2022; Timpane-Padgham et al., 2017), as well as to develop strategies to respond and adapt to it (IPCC, 2022).

Although recent studies have shown that increases in temperature, changes in water-availability, and intensified human activities emerge as the dominant drivers of resilience shifts across the globe, the effect of these drivers are highly context dependent (Boulton et al., 2022; Feng et al., 2021; Forzieri et al., 2022; Jakovac et al., 2015; Y. Zhang et al., 2022). For example, while warming is considered the main global driver in altering ecosystem resilience (Smith et al., 2022), its effect on other ecosystem properties has been shown to vary across ecosystems (Berner & Goetz, 2022; De Keersmaecker et al., 2015; Pecl et al., 2017; Tomiolo & Ward, 2018). In addition, as it happens with other ecosystem functional attributes, different stressors such as land degradation and climate change may interact with each other, exerting their effects in a synergistic (Rillig et al., 2019) and non-linear (Berdugo, Gaitán, et al., 2022; García-Palacios et al., 2018) manner. Thus, the effects of warming on resilience changes might be highly complex, which may contribute to an increased heterogeneity in resilience patterns and make it challenging to distinguish between effects of multiple drivers.

Overall, the diversity of resilience metrics (Kéfi et al., 2019; Van Meerbeek et al., 2021), the context dependency of resilience change, the interactions between environmental drivers, and the expected non-linearity in resilience response, create great uncertainty about what determines resilience changes in different regions across the planet. Resolving these uncertainties requires a standardized global analysis of the interactive effects (both linear and non-linear) of the environmental drivers of resilience changes across different vegetation types, worldwide.

Here, we explored the main environmental, climatic, and anthropogenic drivers of changes in vegetation resilience across forest, grassland, and savanna ecosystems worldwide. We quantified both resistance (estimated as a system's response to a given disturbance) and recovery (the time a system takes to return
to its original state). Both theoretical (Domínguez-García et al., 2019; Kéfi et al., 2013) and empirical studies (Eby et al., 2017; Majumder et al., 2019) have suggested that both these key aspects of resilience tend to decline as ecosystems approach points at which they become vulnerable to shifts in function or structure (Capdevila et al., 2021; Hodgson et al., 2015; Ingrisch & Bahn, 2018; Weinans et al., 2021). We measured time series autocorrelation as a proxy of recovery (Forzieri et al., 2022; Smith et al., 2022) and sensitivity as a proxy of resistance (Seddon et al., 2016).

To do this, we used time series vegetation data of the MODIS Enhanced Vegetation Index (EVI) and climate data from 2000 to 2020. The difference in the two metrics was calculated between the second half of this time frame (2010-2020), relative to the first half (2000-2010) in order to identify which regions are experiencing shifts in resilience. We then applied a machine learning-based model followed by interpretable machine learning approaches (SHAP values) to: i) identify key environmental, climatic, and anthropogenic variables that correlate with changes in sensitivity and autocorrelation across different ecosystems; ii) elucidate whether resilience losses are related to climate warming in a global manner or depending on the type of ecosystem investigated; iii) uncover warming interactions with habitat degradation on driving resilience losses.

# Methods

## Monitored sites

We uniformly sampled the Earth's terrestrial surface, with one point every 0.3° (approximately 30 km apart). Then, points with more than 20% of missing data in the Enhanced Vegetation Index (EVI) time series were removed, and 126,145 points remained. Finally, points identified as “Croplands”, “Urban lands”, “Wetlands”, “Permanent snow and ice” and “Water Bodies” were removed either because their closeness to human activity makes EVI dynamics more driven by human arbitrary management, or because water, snow-, and ice-covered regions are not suited for EVI measures (Huete, 2016; Q. Zhang et al., 2020). We were left with a total of 106,454 points, with 17,931 points in evergreen forests, 4,386 in deciduous forests, 8,352 in mixed forests, 10,384 in shrublands, 15,313 in woody savannas, 17,710 in open savannas, 31,204 in grasslands and 1,183 in barren areas. We focused the driver analysis on forests, savannas and grasslands, as 3 principal vegetation types, because covariate data product coverage for shrubland situated in high latitudes was poor (Abatzoglou et al., 2018).

## Time-series data

We extracted the data from different spatial sources using the sites described previously. Monthly time series from 2000 to 2020 were used. EVI and pixel reliability was extracted from the MOD13C2 product (Didan & Huete, 2015), downloaded from the LAADS database, at a 0.05° (approximately 5 km) resolution. Cloud-fraction was extracted from NASA Earth Observations, as a monthly composite of
the MOD06 product (MODIS Atmosphere Science Team, 2015), at a 0.1° (approximately 10 km) resolution. Monthly Minimum and Maximum temperature, as well as Actual Evapotranspiration (AET, derived using a one-dimensional soil water balance model), and Potential Evapotranspiration (PET, obtained through ASCE Penman-Monteith method), were extracted from the TerraClimate database (Abatzoglou et al., 2018), at a 0.05° resolution. The TerraClimate database is a monthly interpolation from WorldClim and the Climatic Research Unit, with data available up to the end of December 2020. EVI pixels were included in the analysis when the quality was “good”, or “marginal” (citation). Our temperature variable was obtained by averaging each monthly minimum and maximum temperature. Water-availability was estimated as the monthly ratio of AET over PET (Seddon et al., 2016). Cloud-fraction is a monthly proportion of cloud-coverage for each pixel, and was our cloudiness index. Our study focuses on the period where all data were available, which was from March 2000 to December 2020 (both included). Data were then sampled at each monitored site.

## Processing time series for sensitivity index calculation

To estimate sensitivity, we used the methodology developed by Seddon et al., (2016), with modifications to increase robustness to outliers and to ensure that sensitivity was comparable across our study period. We generated detrended time series, where seasonal and long-term temporal trends were removed. To do this, we first removed seasonal effects by subtracting from each observation the mean value of its month during the monitored period (2000-2010 and 2010-2020). Then, to remove any signal of longer-term time trends, we linearly regressed the deseasonalized time series against time. We extracted the residuals from this linear regression model, and then normalised them. This produced a dataset of monthly anomaly scores (Z-scores), which measures how extreme an observation is given its position in the time series. We performed this procedure separately for every monitored location and variable (EVI, temperature, water availability and cloudiness).

Each time series at each site was split into two parts, one time series from 2000 to 2010, and one time series from 2010 to 2020, and the following analysis (aimed at calculating sensitivity and autocorrelation) was performed on each dataset, to obtain two measures per site to be compared (one per decade of available data). Months with mean EVI < 0.1 and mean temperature < 0°C were considered dormant, and were excluded from the analysis (see Seddon et al., 2016).

The sensitivity index for each climate variable was calculated using eq. (1), where $var(EVI)$ and $var(clim)$ are the corrected variance for EVI and a given climate variable (temperature, water availability and cloudiness), and $\beta_{clim}$ is a weight that measures how closely related are anomalies in EVI and in that climate variable.

$$Sensitivity_{clim} = \left( \frac{var(EVI)}{var(clim)} * \beta_{clim} \right) (1)$$
To obtain the corrected variance, only photosynthetically active months were considered. Instead of using directly the variance of EVI and climate, we followed the approach of Seddon et al., by using the residuals from a quadratic model between variance and mean, fitted for each variable, to prevent any influence of the known relationship between mean-variance (Supplementary Figure 1) in the metric. Because we are interested on comparing two time periods, we derived the relationship var-mean on the two combined periods. By doing so, we control for trends of change in the mean NDVI due to greening or browning, which are apparent in the study period and might yield significantly different relationship in 2000-2010 respect to 2010-2020. Besides, the relationship derived is meant to correct for a universal var-mean relationship existing across time and space. These residuals were rescaled from 0 to 100 (based on the range of the residuals calculated for the two periods, to keep them comparable), to ensure they were positive.

To estimate $\beta_{clim}$ for each climate variable, we quantify how an anomaly in the EVI time series is observed when there is an anomaly in each climate time series, taking into account the autocorrelation of the system (by including the 1-month lagged EVI time series in the regression; equivalent to perform an autoregressive model). To do that, because we expect our 4 variables (3 climatic and autocorrelation) to be non-independent, we performed an independent Principal Component Analysis (PCA) regression (Jolliffe, 2002), containing the monthly Z scores for each variable across years, only for the active months and each site. To minimise the effect of outliers, a robust linear regression (instead of a regular linear model used in Seddon et al.) was used between the EVI Z-scores and the 4 first PCA axis. To get a monthly estimation of the effect of our initial Z variables (climate and lag) on the EVI Z scores, significant coefficients extracted from the robust regression were multiplied by the loading scores and summed for each corresponding variable. These monthly coefficients were averaged and the absolute value of this average results in weights, ($\beta_{clim}$ in the equations).

As a result of eq. (1), we then obtained a sensitivity for each climate variable (temperature, water availability, and cloudiness). From now on, the term sensitivity will be used to describe the sum of these three climate sensitivities, resulting in a total sensitivity (Vegetation Sensitivity Index, as in Seddon et al., 2016). We obtained some values of betas that are too large (either negative or positive). This usually occurs when PCA regressions have outliers and they indicate that the regression is ill-conditioned. To ensure that these values were not biasing our results, we conducted for each period an outlier removal procedure by filtering out points with a sum of the betas in the 2.5% upper and lower quantiles in each period (and therefore the sensitivity and autocorrelation of these sites is not calculated).

Although the autocorrelation can be inferred directly from the weights of the memory effect in the PCA regression process, we calculated the autocorrelation as AR1 separately, using a simple correlation test between the EVI Z scores and the 1-month lagged EVI Z scores for the active season. We did so in
order to use more directly comparable methodology with other studies (Boulton et al., 2022; Forzieri et al., 2022; Smith et al., 2022). In any case, the memory effect obtained by Seddon et al.'s method and the AR1 were positively correlated (ρ = 0.57 and 0.61, for 2000-2010 and 2010-2020 periods, respectively).

For each site, the changes in sensitivity and autocorrelation were calculated by subtracting values calculated for 2000-2010 to that calculated for 2010-2020, splitting our time series in half, for both metrics. Therefore, because both metrics are negatively related to resilience, positive changes in this difference between periods indicate a decrease in resilience.

As the change in sensitivity was very heterogeneous, we performed a test to assess significant changes between the two periods. To obtain multiple sensitivity assessments for each period, we used a resampling method to reconstruct time series of EVI, temperature, water availability, and cloud-cover. We randomly selected years for each decade, using uniform sampling with replacement, allowing us to create 100 different time-series, representative of the relationships between EVI and climate conditions, for both periods. We calculated the sensitivity scores in both periods for these 100 iterations, and performed a Student's test, identifying sites with significant changes in sensitivity. 22% and 21% of points showed a significant increase and decrease in sensitivity respectively. We found that significant differences were exhibiting patterns similar to raw sensitivity change, making us confident in our measure of change in sensitivity (Supplementary Figure 2 and 3). Such a process is not suitable to test AR1, as the time series should not be segmented. Nevertheless, our AR1 patterns reflect trends found in (Forzieri et al., 2022; Smith et al., 2022).

As we uncovered relationships between change in climatic averages over the period and change in resilience metrics, we used change in temperature from WorldClim in the model, instead of the initial time-series for the other variables. We further tested for statistical artifacts between the resilience metrics and trends in each component of the time series. To do this, we simulated time series of each variable, with a uniform random distribution for EVI, water availability, and cloud-cover, and a linear positive trend in temperature. All values were simulated within bounds comparable to those of the initial dataset. We found no relationship between trend in each climate variable, productivity, and the changes in resilience (Supplementary Figure 4).

To create the maps, monitored points (separated 0.3° from each other) were aggregated into 0.5° (approximately 50 km) pixels, with each pixel value computed as the mean of the points contained in the pixel, to reduce noise in the data visualisation and improve readability.
Estimating drivers of the change in resilience using machine-learning (XGBooster)

First, we compared gains/losses in resilience between different biomes and forest types (Dinerstein et al., 2017). To identify if each biome tends to gain or lose resilience, we tested the significance of the change in sensitivity and autocorrelation through a Wilcoxon-Mann-Whitney test, as the normality conditions were not respected. To compare changes in resilience over time of different biomes, we performed pairwise permutation tests of both changes in sensitivity and autocorrelation between biomes with 100 permutations. To determine what is driving changes in sensitivity and autocorrelation within or across each biome, we used a machine learning based model: XGBooster, from the "xgboost" package (Chen & Guestrin, 2016). The model confronts change in resilience with variables (later referred to as drivers if significant in the model), and uses a machine-learning algorithm to derive a tree classification model that results in a prediction of the value for the change in resilience. We aim to use our models for inference rather than for prediction (Addicott et al., 2022; Berdugo, Gaitán, et al., 2022), since the resilience metrics can be inferred directly through remote sensing products. We thus focus here on identifying the statistical effects attributable to variables of interest chosen a priori based on prior literature (Supplementary Table 1) instead of maximising predictive power. The data used for the model as well as detailed methods about the model are described below.

a. Data

A number of static and dynamic variables were considered as potential drivers of resilience changes. The static drivers are detailed in Table 1 and are related to site-specific characteristics of a given pixel, either historical climate (averaged for the period 1970-2000), vegetation structure, topography, soil, or human influence. We included drivers that have been previously associated with an intrinsic mechanism of ecosystem resilience (vegetation structure, topography, aridity), or can be an external variable affecting these mechanisms (human influence, change in climate). The specific rationale for including each static variable, including references that support this inclusion are depicted in Supplementary Table 1. The dynamic variables used were changes in average temperature, water availability, cloudiness, EVI, and vegetation structure during the studied period (2000-2010 vs 2010-2020), as we expect climate, productivity (here, remote-sensed greenness) and land-cover change, to impact resilience. A decrease in productivity can also be a valuable indicator of land degradation (Burrell et al., 2020; Eckert et al., 2015), however, it is important to recognize that this relationship is complex and context-dependent (Easdale et al., 2018; Yengoh et al., 2016). The cloudiness, water availability and EVI changes were computed from the time series used to estimate sensitivity. Temperature change was calculated from the WorldClim database (Fick & Hijmans, 2017). We removed dormant months when calculating the overall average for water availability, cloudiness and EVI, based on the same criteria from the sensitivity analysis for each site (decadal monthly average of temperature < 0°C and EVI < 0.1). The changes in each of these variables were calculated as the difference between the 2010-2020 (2010-2018 for temperature from WorldClim) decadal mean and that of 2000-2010 derived from
monthly means. Changes in vegetation structure were estimated as changes in land cover between 2001 and 2020 using the MODIS land cover database (Friedl & Sulla-Menashe, 2019). Land cover types were re-classified into three vegetation structures: Forests, Savannas and Grasslands, and each change from one to another was used as a potential driver in the model.

To remove very highly correlated variables, only one of every highly correlated pair of variable (i.e. |ρ| > 0.7) was kept in the final model. Mean Annual Temperature was correlated to Soil Carbon Content, and Mean Annual Precipitation was correlated to Aridity Index. Temperature has been shown to impact resilience (Feng et al., 2021) and Aridity has been widely used to characterise various ecosystems (Berdugo et al., 2020; Tijerín-Triviño et al., 2022; UNEP, 1992). For these reasons, we chose to keep Mean Annual Temperature and Aridity in the model.

## b. Machine-learning model parametrization

The functioning of our model is conditioned by hyper-parameters, which determine the way in which the decision tree is constructed. Minimum child weight and Minimum splitting loss both change the degree to which the algorithm can create new nodes. Learning rate, which is how much the algorithm weights each driver, and Maximum depth, changing the size of the tree, for these two, high value can cause overfitting. Subsample size is a ratio of the training data used by the model, sampled once every tree generation, and can prevent overfitting. To estimate the best values for these parameters, a Bayesian optimization was used, from the "ParBayesianOptimization" package (Snoek et al., 2012) to minimise RMSE on the validation dataset using a cross validation procedure, thus preventing overfitting. We used the "gbtree" booster and a squared loss regression as a parameter of the learning task, because they are suited to our continuous variables. After tuning, values for each parameter are shown in Supplementary Table 2. 70% of our data was used to train the model, and 30% for validation. Then, the drivers with the biggest importance were identified, i.e., the proportion of additional variance explained by a given driver to the model, compared to a model without this driver (Supplementary Figure 5).

To obtain the effects of the drivers on the sensitivity/autocorrelation changes, we used SHapley Additive exPlanations (SHAP) values. SHAP values are based on information theory and are used to interpret machine-learning outcomes (Lundberg & Lee, 2017). Specifically, these values reveal the contribution of the value of a given driver to a given specific prediction. The sum of all SHAP values (one per predictor used) for a data point results in the final prediction, in a process similar to partial dependence regressions. This means that SHAP values of a given predictor are the effects of that predictor alone, all other predictors being controlled. Thus, plotting SHAP values of a given predictor against the values of that same driver allows us to characterise the effect of that driver on the data (i.e., interpreting the effect) and to highlight potential thresholds and interactions between drivers (Figures 2 to 4). To assess the robustness of our approach, different subsets of the initial data were used to run the
model, with an optimization of the parameters for each subset. We didn't observe any major difference
observed on accuracy or when plotting SHAP values (Supplementary Figure 6). Finally, to add
transparency, an interactive interface to explore relationships between SHAP values and variables
included in the regressive model, as well as interactions, is available on this link:
"https://cfournierdel.shinyapps.io/SHAPValuesFournier/".

To test if observed patterns were homogeneous across vegetation structures, the regressive model was
also run for 3 separate vegetation types: forests, savannas and grasslands, based on the classification of
the MODIS land cover database in 2001.

To improve the robustness of the patterns uncovered with SHAP values from the continuous regressive
model, we ran a binary classifier using the points identified as losing or gaining sensitivity. To do that,
we used a XGBOOST model with a binary logistic regression, and optimized our hyper-parameters
following the same Bayesian process described previously. SHAP value plots for this model are
displayed in Figure 2, and Supplementary Figure 7 and 8. This model reveals that significant changes
in sensitivity are driven by similar processes as continuous changes in sensitivity. Using the continuous
regressor allows for a detailed analysis of the scales of changes in both sensitivity and autocorrelation.

Thresholds for SHAP Values

Visual inspection of SHAP value plots suggested the possible presence of aridity and change in
temperature thresholds (i.e., points in the gradient where there is a sudden change in either slope or
intercept of the relationship between the driver and its associated SHAP value). To verify this
quantitatively, we used R package "chngpt" (Fong et al., 2017) to identify points at which the
relationship between SHAP values and drivers changed abruptly. To reduce the effect of outliers when
identifying these thresholds, points were weighted based on Mahalanobis distance to the mean of the
overall distribution for aridity or temperature (Berdugo, Gaitán, et al., 2022; Mahalanobis, 1936). We
ran this analysis on SHAP values for the whole dataset, as well as specifically for forests, savannas and
grasslands. The threshold was only quantified when the model including a change point had a lower
AIC than a regular linear model.

Results

Geographical patterns of resilience change

We explored the responses of both aspects of resilience (resilience and recovery) examine where there
were similarities and differences between the two. Across all sites, 27% exhibit increased vegetation
resilience between the two periods analysed (i.e. they showed a lower sensitivity and autocorrelation in
the second period compared to the first one). In contrast, another 26% of sites have been losing
resilience in terms of both properties. However, 47% of study sites exhibit opposing trends for resistance
and recovery metrics, resulting in an unclear pattern of resilience change. Therefore, slightly more than
half the sites have been found to either lose or gain both sensitivity and autocorrelation, but when
considering all points, the correlation between the changes in the two metrics is low (0.14, Pearson
Correlation test). When looking respectively at forests, grasslands, and savannas, 27%, 26%, and 28%
are losing resilience and 26%, 29%, and 31% are gaining resilience.

Spatially, the changes in resilience are not homogeneous, but instead they form mosaics of hot and cold
spots of resilience gain and loss (Figure 1 and Supplementary Figures 2 and 3). Some regions show a
rise in both autocorrelation and sensitivity, such as the Australian deserts, the *Caatinga* region in South
America, and parts of the boreal region. Meanwhile, Eastern Europe, parts of Asia, Africa and North
America are gaining resistance and are recovering faster. In the boreal region, we observe a very high
variability for both indices. Deciduous forests have gotten more sensitive, with tropical forests mostly
gaining sensitivity. In contrast, evergreen and sclerophyllous forests have become less sensitive
(Supplementary Figure 9). This pattern is not observed for changes in autocorrelation.

<figure><img src="image_4.png" /><figcaption>**Figure 1:** Consensus on different metrics of resilience changes in the last two decades. Global map showing the direction of change in resilience (sensitivity and autocorrelation as AR1) between 2000-2010 and 2010-2020. Autocorrelation was calculated on consecutive months during the growing season. Yellow values translate an increase in autocorrelation and in sensitivity, related to a decrease in resilience. Blue and green pixels are decreasing one property of resilience but increase on the other. Most pixels are either purple or yellow, which means that both properties of resilience are usually changing in the same direction. Urban, croplands, wetlands, barren and snow-dominated pixels are mapped in grey. Pixel resolution: 0.5° ≈ 55km.</figcaption></figure>
## Drivers of resilience change

We used a machine-learning model to assess the drivers of resilience change, and to compare their relative impacts. In general, the overall precision of our models was not very high (Pearson's R²=0.24 and 0.38 for change in sensitivity and change in autocorrelation, when looking at the relationship between prediction and actual value). However, we aim to use our models for inference rather than for prediction (see methods *Estimating drivers of the change in resilience using machine-learning (XGBooster)*). Despite the model's low predictive power, we can still infer relationships between resilience metrics and predictive variables, chosen *a priori* based on prior literature (Supplementary Table 1).

Overall, mean annual temperature and changes in temperature are the most important drivers of sensitivity change and, along with change in EVI, are the most important for autocorrelation change (Supplementary Figure 5). We also highlight the importance of change in water-availability, cloud-cover, and EVI for both metrics (Supplementary Figure 5). Elevation, aridity and human modification are also key variables in explaining the observed patterns. Surprisingly, intact landscapes and forest integrity, which are related to protected areas, are not explaining changes in resilience, potentially because our measure of human modification, a very important driver in the model, is already capturing the effect of land-protection. Finally, we find that variables driving resilience changes had different importance across different vegetation structures (Supplementary Figure 10). Mean temperature is the most important driver for both sensitivity and autocorrelation across vegetation types, but human modification shows higher importance than climate on savannas, compared to forests and grasslands. In a similar way, changes in our average climate properties have higher importance on forests than on grasslands and savannas.

In order to characterise the relationship between the change in resilience and each driver, we used SHapley Additive exPlanations (SHAP) values, which measure the contribution of the value of a given driver to the model's prediction (analogous to partial dependence effect plots). By plotting these SHAP values against the values of their respective driver, we can visualize these relationships and easily identify thresholds, and interactions between drivers (Figures 2 to 4).

We find that higher altitude is related to increasing resilience (Supplementary Figure 11). Moreover, ecosystems at different aridity levels reveal different responses in resilience change (Figure 2). Specifically, the more arid systems tend to exhibit an increase in autocorrelation over time, thus becoming slower to recover from disturbances, while they show an increase in resistance (i.e., less sensitivity. Furthermore, we see that ecosystems gain most resistance around the 0.5 aridity value, but do not gain more resistance after that threshold (Figure 2).
<figure><img src="image_5.png" /></figure>

**Figure 2:** **Extreme aridity influences resilience changes.** Contribution (SHAP values being the predicted change in a prediction, controlling for all other predictors) of aridity to the probability of being classified as significantly a) increasing or b) decreasing in Sensitivity. c) and d) display SHAP values contributing to the estimation of the continuous change in sensitivity and autocorrelation, respectively. SHAP values for a) and b) are extracted from the binary classifier models, while SHAP values for c) and d) are extracted from the continuous regressor models. Here, we calculate aridity as 1-Aridity Index, so that values between 0 and 1 are water-limited. A more arid system will have decreasingly smaller probability of increasing its sensitivity, meaning more arid systems do not lose their resistance. Very wet systems tend to decrease in sensitivity, meaning they exhibit greater resistance.

Temperature also affects resilience in a non-linear way (Figure 3 and Supplementary Figure 4), as ecosystems experiencing warming show gradual increasing autocorrelation (meaning a decreasing ability to recover), and a rise greater than 0.2°C results in increased sensitivity (losing the ability to
resist disturbances). The 0.2°C threshold for a rise in sensitivity is similar when looking at different vegetation structures, but different vegetation structures will exhibit a different response for autocorrelation. Indeed, we observe an increased autocorrelation for grasslands and savannas experiencing a temperature decline of more than 0.2°C, which is not as prevalent for forests. A sharp rise in temperature can induce a decrease in autocorrelation (increase in recovering capacity), with a threshold around the 0.5°C increase for savannas and 0.6°C for forests.

<figure><img src="image_6.png" /></figure>

**Figure 3: Change in temperature influences resilience changes.** Contribution (SHAP values being the predicted change in a prediction, controlling for all other predictors) of average change in temperature to the change in resilience properties. Predicted values for both change in sensitivity (a) and autocorrelation (b) models, run for each of 3 vegetation structures (forest, grassland and savanna). We see a subsequent increase in both metrics after a 0.2 and 0.4°C increase in the last decade (meaning reduced resilience).

Ecosystems that have become more productive tend to have increased resilience, and vice versa, as we have a negative relationship between change in EVI and change in sensitivity and autocorrelation (Figure 4). Changes in temperature exhibit a positive interaction with the latter relationship, indicating that places displaying an increase of more than 0.2°C will lose, or gain more resilience, for a respective loss or gain in EVI. These relationships are constant along vegetation structures, except for
autocorrelation in grasslands, which does not show an interaction between rise in temperature and change in EVI. Nevertheless, there is a general pattern for change in sensitivity across different vegetation structures, with a greater than 0.2°C increase resulting in a gradual loss in the ability to resist, with further warming leading to subsequent resilience loss.

<figure><img src="image_7.png" /><figcaption>**Figure 4: Change in productivity interacts with change in temperature to influence resilience.**</figcaption></figure>

Contribution (SHAP values being the predicted change in a prediction, controlling for all other predictors) of average change in productivity (here, estimated through EVI), interacting with the average change in temperature, to the change in resilience properties. Predicted values for both change in sensitivity (a) and autocorrelation (b) models, red meaning gaining more temperature. There is a negative relationship between change in both metrics and EVI. Linear models fit lines show steeper slopes for changes in temperature above 0.2°C. R² for the fitted lines are 0.6 and 0.55 for points experiencing a respective warming greater than 0.2°C and less than 0.2°C.

## Discussion

Changes in resilience show heterogenous patterns across the globe, in accordance to previous work (Forzieri et al., 2022; Smith et al., 2022), with highest heterogeneity in the boreal region. Autocorrelation and sensitivity are often changing in the same direction (Figure 1 and Supplementary Figure 2), yet are not highly correlated, supporting the inclusion of both metrics when identifying drivers of resilience (Smith & Boers, 2023b). We show that these patterns are created by complex interactions between changes in climatic conditions and location-specific drivers, with subtle differences for different vegetation types. We also highlight the non-linearity of the impact of these drivers on resilience, with temperature rising more than 0.2°C resulting in subsequent loss in resilience, interacting with decreases in productivity, leading to a further decline in resilience.
Our results suggest that global resilience change is primarily driven by changes in climate. Indeed, temperature, water availability, and cloud-cover, are rapidly shifting, with direct impacts on both resistance (sensitivity) and recovery (autocorrelation). Climate change influences drought and fire dynamics (Dale et al., 2001; de Groot et al., 2013; Lenihan et al., 2008) sometimes coupled with human activity (Gou et al., 2022; Rito et al., 2017). Since the disturbance regime can determine which vegetation structure is occurring at a given place, changes in this regime can be a major driver of vegetation shift resulting in land-cover change, thus affecting resilience (Beckage et al., 2009; Johnstone et al., 2016; Scheffer et al., 2009). Nevertheless, temperature seems to be the primary force resulting in losses of resilience, with a notable nonlinear response of ecosystems to global warming (Figure 3 and Supplementary Figure 7). Indeed, our results suggest that an increase greater than 0.2°C will result in the loss of ecosystem resilience, regardless of the mechanisms, which should be investigated in further studies. Aside from more frequent and intense disturbances, we hypothesize warmer temperatures could rise physiological stress in the vegetation (Zhao et al., 2020), leading to higher sensitivity (Yi et al., 2022) and slower recruitment (Hatfield & Prueger, 2015; Lippmann et al., 2019). Finally, our study quantified the effect of a rapid increase in temperature, which may involve different mechanisms, but can be symptomatic of the effects of a longer-term increase (Nolan et al., 2021).

Changes in productivity, measured by change in greenness of the system over time is also positively related to resilience (Figure 4). The link between reduced diversity (functional and species), reduced greenness, human activities, and resilience loss has already been documented (Berner & Goetz, 2022; Boulton et al., 2022; Sebald et al., 2021; Spasojevic et al., 2016). However, we provide the first global empirical evidence showing a clear effect of losing EVI, or productivity (potentially linked to land degradation, Burrell et al., 2020; Eckert et al., 2015) on resilience, even when controlling for a number of site-specific drivers. This supports the idea that degradation of a system would also reduce its ability to withstand potential disturbances (Thompson et al., 2009). In other words, when human activities lead to the degradation of ecosystems, we do not only lose the productivity of those systems (Barlow et al., 2016; Foley et al., 2005; Heinrichs et al., 2016), we also degrade the ability of that ecosystem to resist and recover from further stressors and disturbances.

Most interestingly, the relationship between EVI changes and resilience metrics is even stronger when temperature is changing beyond 0.2°C (Figure 4). This interaction indicates that global warming could positively interact with ecosystem degradation, exacerbating resilience loss of degraded lands. It also adds to the body of literature showing that combined effects of multiple global change drivers may result in unexpected detrimental effects for ecosystem functioning and stability (Rillig et al., 2019). On the other hand, rising CO₂ concentration and temperature may increase vegetation productivity (Dusenge et al., 2019; Gustafson et al., 2017), which our results relate to a gain in resilience.
Conversely, rising temperature and more biomass have been suggested to lead to greater drought, fire, and parasite outbreak severity, implying a drastic change in vegetation structure (Bentz et al., 2010; Halofsky et al., 2020; Lenihan et al., 2008). For that reason, uncertainty remains surrounding the persistence and resilience of re-greening ecosystems that are also getting warmer.

Drivers of resilience change may vary among the different vegetation types we examined (Figure 3, Supplementary Figures 10, 12 and 13). Forests show more sensitivity to climate, perhaps related to increased droughts (Choat et al., 2018) and savannas are more responsive to changes in human-driven pressures. Indeed, savannas might be more subject to changes in human activities such as grazing (Kariuki et al., 2021), inducing fire regimes or invasions (Archibald, 2016; Gorgone-Barbosa et al., 2016), but could also be misidentified degraded forests (exhibiting low percent tree cover) by the MODIS land classification (Friedl & Sulla-Menashe, 2019). As such, we emphasise that unique site-specific environmental and biological properties must always be considered when assessing the resilience of a system (Dakos & Kéfi, 2022).

In that regard, we find that ecosystems at different aridity thresholds might have different responses concerning resilience change (Figure 2). Indeed, arid systems are on average more resistant and recover slower than non-arid regions (Seddon et al., 2016; Wu & Liang, 2020). We show that these arid systems exhibit increasing resistance but are becoming slower to recover from disturbances in the recent years. This might be related to their ability to show more plasticity in response to drought episodes, thus displaying remarkable ability to resist water-stress and disturbances (Bachofen et al., 2021; Richter et al., 2012). Climate change thus imposes a scenario in which the adaptive resistance to a disturbed environment on most arid lands (i.e., changes in weather fluctuations) becomes more obvious. Yet, at the same time, it emphasises their poor ability to recover when their resistance limit is crossed. Moreover, we see that systems gain most resistance around a 0.5 aridity value, but do not gain more resistance after that threshold. This threshold has been shown to be related to a decrease in productivity, with drastic changes in soil properties, vegetation structures, and species composition associated with adaptations to extreme aridity (Berdugo et al., 2017, 2020). Vegetation adapted to aridity levels drier than 0.5 is very specialised at the organism level, probably including strategies to evade environmental stress (Berdugo, Vidiella, et al., 2022). Our findings support that the specialization of organisms might make them not responsive to climate change, which has a significant impact on the dynamics of the ecosystem as a whole.

In this study, we uncovered various mechanisms of how the resilience of ecosystems changes over a short period of time. However, there are several limitations to our approach that should be considered when interpreting our results. While our global study provides valuable insights into how our earth system works, it is important to note that local knowledge such as disturbance dynamics, community
composition, and the social context is always needed to make informed conservation and restoration decisions. Furthermore, decreased productivity is not always related to land degradation. There may be other forces at play that can affect productivity, and further research is needed to fully understand these relationships. Additionally, it is important to remember that correlation does not imply causation. Our results should be interpreted with this in mind and serve to point towards possible mechanisms of ecosystem resilience change. Finally, the performance of our models is moderate, which may be due to heterogeneity in the data (Smith & Boers, 2023b) or unknown drivers that have not been proposed. However, despite these sources of uncertainty, our approach allowed us to see through the statistical noise to identify consistent underlying trends. We hope that these trends encourage future experimental research to understand the mechanisms underpinning these global patterns.

## Conclusion

Our study provides insights into the global patterns of changes in vegetation resilience, showing the importance of site-specific conditions, as well as global climatic drivers, across vegetation types. Most importantly, we find that warming is decreasing the resilience of ecosystems, in synergy with shifts in greenness, leading to further loss of resilience in degraded habitats that also get warmer. We hypothesize these signals to be related to the higher frequency and intensity of disturbances but call further research to identify the mechanisms underlying our findings. Nevertheless, our findings reveal that even ecosystems seemingly unresponsive to climate change may be losing their ability to resist and recover from disturbances. For these ecosystems, future disturbances could trigger abrupt, irreversible shifts. Thus, nuanced, system-specific insights are crucial for implementing effective protection and restoration strategies. This targeted approach is essential for preserving ecosystem persistence and sustaining the ecosystem natural functions and services upon which humanity relies.
**Acknowledgements:**

GRS acknowledges funding from DOB Ecology and from the Marc R. Benioff Revocable Trust, which, in collaboration with the World Economic Forum, made this work possible. M.B. is supported by a Ramón y Cajal grant from the Spanish Ministry of Science (RYC2021-031797-I). We wish to thank Tom Lauber for his assistance with the analysis.

**Conflict of Interest Statement:**

The authors do not declare any conflict of interested.

**Data Availability:**

The data and code that support the findings of this study are not yet made openly available at this stage of publication. SHAP values plots from the models can be visualized using "https://cfournierdel.shinyapps.io/SHAPValuesFournier/"
**References:**

Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A., & Hegewisch, K. C. (2018). TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. *Scientific Data*, 5(1), Article 1. https://doi.org/10.1038/sdata.2017.191

Addicott, E. T., Fenichel, E. P., Bradford, M. A., Pinsky, M. L., & Wood, S. A. (2022). Toward an improved understanding of causation in the ecological sciences. *Frontiers in Ecology and the Environment*, 20(8), 474–480. https://doi.org/10.1002/fee.2530

Amatulli, G., Domisch, S., Tuanmu, M.-N., Parmentier, B., Ranipeta, A., Malczyk, J., & Jetz, W. (2018). A suite of global, cross-scale topographic variables for environmental and biodiversity modeling. *Scientific Data*, 5(1), Article 1. https://doi.org/10.1038/sdata.2018.40

Archibald, S. (2016). Managing the human component of fire regimes: Lessons from Africa. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 371(1696), 20150346.
https://doi.org/10.1098/rstb.2015.0346

Bachofen, C., Perret-Gentil, A., Wohlgemuth, T., Vollenweider, P., & Moser, B. (2021). Phenotypic plasticity versus ecotypic differentiation under recurrent summer drought in two drought-tolerant pine species. *Journal of Ecology*, 109(11), 3861–3876. https://doi.org/10.1111/1365-2745.13762

Baldeck, C. A., Harms, K. E., Yavitt, J. B., John, R., Turner, B. L., Valencia, R., Navarrete, H., Davies, S. J., Chuyong, G. B., Kenfack, D., Thomas, D. W., Madawala, S., Gunatilleke, N., Gunatilleke, S., Bunyavejchewin, S., Kiratiprayoon, S., Yaacob, A., Supardi, M. N. N., & Dalling, J. W. (2013). Soil resources and topography shape local tree community structure in tropical forests. *Proceedings of the Royal Society B: Biological Sciences*, 280(1753), 20122532.
https://doi.org/10.1098/rspb.2012.2532

Barlow, J., Lennox, G. D., Ferreira, J., Berenguer, E., Lees, A. C., Nally, R. M., Thomson, J. R., Ferraz, S. F. de B., Louzada, J., Oliveira, V. H. F., Parry, L., Ribeiro de Castro Solar, R., Vieira, I. C. G., Aragão, L. E. O. C., Begotti, R. A., Braga, R. F., Cardoso, T. M., de Oliveira, R. C., Souza Jr, C. M., ... Gardner, T. A. (2016). Anthropogenic disturbance in tropical forests can double biodiversity loss from deforestation. *Nature*, 535(7610), Article 7610.
https://doi.org/10.1038/nature18326
Beckage, B., Ellingwood, C., & University of Vermont. (2009). Fire Feedbacks with Vegetation and Alternative Stable States. *Complex Systems*, 18(1), 159–173.
https://doi.org/10.25088/ComplexSystems.18.1.159

Bentz, B. J., Régnière, J., Fettig, C. J., Hansen, E. M., Hayes, J. L., Hicke, J. A., Kelsey, R. G., Negrón, J. F., & Seybold, S. J. (2010). Climate Change and Bark Beetles of the Western United States and Canada: Direct and Indirect Effects. *BioScience*, 60(8), 602–613.
https://doi.org/10.1525/bio.2010.60.8.6

Berdugo, M., Delgado-Baquerizo, M., Soliveres, S., Hernández-Clemente, R., Zhao, Y., Gaitán, J. J., Gross, N., Saiz, H., Maire, V., Lehmann, A., Rillig, M. C., Solé, R. V., & Maestre, F. T. (2020). Global ecosystem thresholds driven by aridity. *Science*, 367(6479), 787–790.
https://doi.org/10.1126/science.aay5958

Berdugo, M., Gaitán, J. J., Delgado-Baquerizo, M., Crowther, T. W., & Dakos, V. (2022). Prevalence and drivers of abrupt vegetation shifts in global drylands. *Proceedings of the National Academy of Sciences*, 119(43), e2123393119. https://doi.org/10.1073/pnas.2123393119

Berdugo, M., Maestre, F. T., Kéfi, S., Gross, N., Bagousse-Pinguet, Y. L., & Soliveres, S. (2017). *Species-specific adaptations determine how aridity and biotic interactions drive the assembly of dryland plant communities* (p. 147181). bioRxiv. https://doi.org/10.1101/147181

Berdugo, M., Vidiella, B., Solé, R. V., & Maestre, F. T. (2022). Ecological mechanisms underlying aridity thresholds in global drylands. *Functional Ecology*, 36(1), 4–23.
https://doi.org/10.1111/1365-2435.13962

Berner, L. T., & Goetz, S. J. (2022). Satellite observations document trends consistent with a boreal forest biome shift. *Global Change Biology*, 28(10), 3275–3292.
https://doi.org/10.1111/gcb.16121

Besnard, S., Koirala, S., Santoro, M., Weber, U., Nelson, J., Gütter, J., Herault, B., Kassi, J., N'Guessan, A., Neigh, C., Poulter, B., Zhang, T., & Carvalhais, N. (2021). Mapping global forest age from forest inventories, biomass and climate data. *Earth System Science Data*, 13(10), 4881–4896.
https://doi.org/10.5194/essd-13-4881-2021

Blauhut, V., Stahl, K., Stagge, J. H., Tallaksen, L. M., De Stefano, L., & Vogt, J. (2016). Estimating drought risk across Europe from reported drought impacts, drought indices, and vulnerability
factors. *Hydrology and Earth System Sciences*, 20(7), 2779-2800. https://doi.org/10.5194/hess-20-2779-2016

Boulton, C. A., Lenton, T. M., & Boers, N. (2022). Pronounced loss of Amazon rainforest resilience since the early 2000s. *Nature Climate Change*, 12(3), 1–8. https://doi.org/10.1038/s41558-022-01287-8

Burrell, A. L., Evans, J. P., & De Kauwe, M. G. (2020). Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification. *Nature Communications*, 11(1), Article 1. https://doi.org/10.1038/s41467-020-17710-7

Capdevila, P., Stott, I., Oliveras Menor, I., Stouffer, D. B., Raimundo, R. L. G., White, H., Barbour, M., & Salguero-Gómez, R. (2021). Reconciling resilience across ecological systems, species and subdisciplines. *Journal of Ecology*, 109(9), 3102–3113. https://doi.org/10.1111/1365-2745.13775

Chen, T., & Guestrin, C. (2016). XGBoost: A Scalable Tree Boosting System. *Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining*, 785–794. https://doi.org/10.1145/2939672.2939785

Choat, B., Brodribb, T. J., Brodersen, C. R., Duursma, R. A., López, R., & Medlyn, B. E. (2018). Triggers of tree mortality under drought. *Nature*, 558(7711), Article 7711. https://doi.org/10.1038/s41586-018-0240-x

Crowther, T. W., Glick, H. B., Covey, K. R., Bettigole, C., Maynard, D. S., Thomas, S. M., Smith, J. R., Hintler, G., Duguid, M. C., Amatulli, G., Tuanmu, M.-N., Jetz, W., Salas, C., Stam, C., Piotto, D., Tavani, R., Green, S., Bruce, G., Williams, S. J., ... Bradford, M. A. (2015). Mapping tree density at a global scale. *Nature*, 525(7568), Article 7568. https://doi.org/10.1038/nature14967

Dakos, V., & Kéfi, S. (2022). Ecological resilience: What to measure and how. *Environmental Research Letters*, 17(4), 043003. https://doi.org/10.1088/1748-9326/ac5767

Dale, V. H., Joyce, L. A., McNulty, S., Neilson, R. P., Ayres, M. P., Flannigan, M. D., Hanson, P. J., Irland, L. C., Lugo, A. E., Peterson, C. J., Simberloff, D., Swanson, F. J., Stocks, B. J., & Wotton, B. M. (2001). Climate Change and Forest Disturbances: Climate change can affect forests by altering the frequency, intensity, duration, and timing of fire, drought, introduced
species, insect and pathogen outbreaks, hurricanes, windstorms, ice storms, or landslides.

*BioScience*, 51(9), 723-734. https://doi.org/10.1641/0006-3568(2001)051[0723:CCAFD]2.0.CO;2

de Groot, W. J., Flannigan, M. D., & Cantin, A. (2013). Climate change impacts on future boreal fire regimes. *Forest Ecology and Management*, 294, 35-44.
https://doi.org/10.1016/j.foreco.2012.09.027

De Keersmaecker, W., Lhermitte, S., Tits, L., Honnay, O., Somers, B., & Coppin, P. (2015). A model quantifying global vegetation resistance and resilience to short-term climate anomalies and their relationship with vegetation cover. *Global Ecology and Biogeography*, 24(5), 539-548.
https://doi.org/10.1111/geb.12279

Didan, K., & Huete, A. (2015). *MOD13C2 MODIS/Terra Vegetation Indices Monthly L3 Global 0.05Deg CMG V006* [dataset]. NASA EOSDIS Land Processes DAAC.
https://doi.org/10.5067/MODIS/MOD13C2.006

Dinerstein, E., Olson, D., Joshi, A., Vynne, C., Burgess, N. D., Wikramanayake, E., Hahn, N., Palminteri, S., Hedao, P., Noss, R., Hansen, M., Locke, H., Ellis, E. C., Jones, B., Barber, C. V., Hayes, R., Kormos, C., Martin, V., Crist, E., ... Saleem, M. (2017). An Ecoregion-Based Approach to Protecting Half the Terrestrial Realm. *BioScience*, 67(6), 534-545.
https://doi.org/10.1093/biosci/bix014

Domínguez-García, V., Dakos, V., & Kéfi, S. (2019). Unveiling dimensions of stability in complex ecological networks. *Proceedings of the National Academy of Sciences*, 116(51), 25714-25720.
https://doi.org/10.1073/pnas.1904470116

Dusenge, M. E., Duarte, A. G., & Way, D. A. (2019). Plant carbon metabolism and climate change: Elevated CO2 and temperature impacts on photosynthesis, photorespiration and respiration.
*New Phytologist*, 221(1), 32-49. https://doi.org/10.1111/nph.15283

Easdale, M. H., Bruzzone, O., Mapfumo, P., & Tittonell, P. (2018). Phases or regimes? Revisiting NDVI trends as proxies for land degradation. *Land Degradation & Development*, 29(3), 433-445. https://doi.org/10.1002/ldr.2871
Eby, S., Agrawal, A., Majumder, S., Dobson, A. P., & Guttal, V. (2017). Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem. *Global Ecology and Biogeography*, 26(6), 638–649. https://doi.org/10.1111/geb.12570

Eckert, S., Hüsler, F., Liniger, H., & Hodel, E. (2015). Trend analysis of MODIS NDVI time series for detecting land degradation and regeneration in Mongolia. *Journal of Arid Environments*, 113, 16–28. https://doi.org/10.1016/j.jaridenv.2014.09.001

ESA Land Cover CCI project team; & Defourny, P. (2016). *Dataset Record: ESA Land Cover Climate Change Initiative (Land_CoverCCI): Land Surface Seasonality Products*. https://catalogue.ceda.ac.uk/uuid/7c114fc6e2884c1f9ca107e7a502fdbf

Feng, Y., Su, H., Tang, Z., Wang, S., Zhao, X., Zhang, H., Ji, C., Zhu, J., Xie, P., & Fang, J. (2021). Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale. *Communications Earth & Environment*, 2(1), 1–11. https://doi.org/10.1038/s43247-021-00163-1

Fick, S. E., & Hijmans, R. J. (2017). WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. *International Journal of Climatology*, 37(12), 4302–4315. https://doi.org/10.1002/joc.5086

Foley, J. A., DeFries, R., Asner, G. P., Barford, C., Bonan, G., Carpenter, S. R., Chapin, F. S., Coe, M. T., Daily, G. C., Gibbs, H. K., Helkowski, J. H., Holloway, T., Howard, E. A., Kucharik, C. J., Monfreda, C., Patz, J. A., Prentice, I. C., Ramankutty, N., & Snyder, P. K. (2005). Global Consequences of Land Use. *Science*, 309(5734), 570–574. https://doi.org/10.1126/science.1111772

Fong, Y., Huang, Y., Gilbert, P. B., & Permar, S. R. (2017). chngpt: Threshold regression model estimation and inference. *BMC Bioinformatics*, 18(1), 454. https://doi.org/10.1186/s12859-017-1863-x

Food and Agriculture Organization of the United Nations. (2005). *Global Forest Resources Assessment Update 2005*. FAO.

Forzieri, G., Dakos, V., McDowell, N. G., Ramdane, A., & Cescatti, A. (2022). Emerging signals of declining forest resilience under climate change. *Nature*, 608(7923), Article 7923. https://doi.org/10.1038/s41586-022-04959-9
Friedl, M., & Sulla-Menashe, D. (2019). *MCD12Q1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V006* [dataset]. https://doi.org/10.5067/MODIS/MCD12Q1.006

Ganteaume, A., Camia, A., Jappiot, M., San-Miguel-Ayanz, J., Long-Fournel, M., & Lampin, C. (2013). A Review of the Main Driving Factors of Forest Fire Ignition Over Europe. *Environmental Management*, 51(3), 651–662. https://doi.org/10.1007/s00267-012-9961-z

García-Palacios, P., Gross, N., Gaitán, J., & Maestre, F. T. (2018). Climate mediates the biodiversity–ecosystem stability relationship globally. *Proceedings of the National Academy of Sciences*, 115(33), 8400–8405. https://doi.org/10.1073/pnas.1800425115

Génin, A., Majumder, S., Sankaran, S., Schneider, F. D., Danet, A., Berdugo, M., Guttal, V., & Kéfi, S. (2018). Spatially heterogeneous stressors can alter the performance of indicators of regime shifts. *Ecological Indicators*, 94, 520–533. https://doi.org/10.1016/j.ecolind.2017.10.071

Ghosh, T., L. Powell, R., D. Elvidge, C., E. Baugh, K., C. Sutton, P., & Anderson, S. (2010). Shedding Light on the Global Distribution of Economic Activity. *The Open Geography Journal*, 3(1). https://benthamopen.com/ABSTRACT/TOGEOGJ-3-147

Gibbs, H. K., & Ruesch, A. (2008). *New IPCC Tier-1 Global Biomass Carbon Map for the Year 2000* [dataset]. Environmental System Science Data Infrastructure for a Virtual Ecosystem. https://doi.org/10.15485/1463800

Gorgone-Barbosa, E., Pivello, V. R., Baeza, M. J., & Fidelis, A. (2016). Disturbance as a factor in breaking dormancy and enhancing invasiveness of African grasses in a Neotropical Savanna. *Acta Botanica Brasilica*, 30, 131–137. https://doi.org/10.1590/0102-33062015abb0317

Gou, Y., Balling, J., Sy, V. D., Herold, M., Keersmaecker, W. D., Slagter, B., Mullissa, A., Shang, X., & Reiche, J. (2022). Intra-annual relationship between precipitation and forest disturbance in the African rainforest. *Environmental Research Letters*, 17(4), 044044. https://doi.org/10.1088/1748-9326/ac5ca0

Grantham, H. S., Duncan, A., Evans, T. D., Jones, K. R., Beyer, H. L., Schuster, R., Walston, J., Ray, J. C., Robinson, J. G., Callow, M., Clements, T., Costa, H. M., DeGemmis, A., Elsen, P. R., Ervin, J., Franco, P., Goldman, E., Goetz, S., Hansen, A., ... Watson, J. E. M. (2020). Anthropogenic modification of forests means only 40% of remaining forests have high
ecosystem integrity. *Nature Communications*, 11(1), Article 1. https://doi.org/10.1038/s41467-020-19493-3

Gustafson, E. J., Miranda, B. R., De Bruijn, A. M. G., Sturtevant, B. R., & Kubiske, M. E. (2017). Do rising temperatures always increase forest productivity? Interacting effects of temperature, precipitation, cloudiness and soil texture on tree species growth and competition. *Environmental Modelling & Software*, 97, 171-183.
https://doi.org/10.1016/j.envsoft.2017.08.001

Halofsky, J. E., Peterson, D. L., & Harvey, B. J. (2020). Changing wildfire, changing forests: The effects of climate change on fire regimes and vegetation in the Pacific Northwest, USA. *Fire Ecology*, 16(1), 4. https://doi.org/10.1186/s42408-019-0062-8

Hatfield, J. L., & Prueger, J. H. (2015). Temperature extremes: Effect on plant growth and development. *Weather and Climate Extremes*, 10, 4-10. https://doi.org/10.1016/j.wace.2015.08.001

Heinrichs, J. A., Bender, D. J., & Schumaker, N. H. (2016). Habitat degradation and loss as key drivers of regional population extinction. *Ecological Modelling*, 335, 64-73.
https://doi.org/10.1016/j.ecolmodel.2016.05.009

Hengl, T., Jesus, J. M. de, Heuvelink, G. B. M., Gonzalez, M. R., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., & Kempen, B. (2017). SoilGrids250m: Global gridded soil information based on machine learning. *PLOS ONE*, 12(2), e0169748. https://doi.org/10.1371/journal.pone.0169748

Hodgson, D., McDonald, J. L., & Hosken, D. J. (2015). What do you mean, 'resilient'? *Trends in Ecology & Evolution*, 30(9), 503-506. https://doi.org/10.1016/j.tree.2015.06.010

Huang, M., Piao, S., Ciais, P., Peñuelas, J., Wang, X., Keenan, T. F., Peng, S., Berry, J. A., Wang, K., Mao, J., Alkama, R., Cescatti, A., Cuntz, M., De Deurwaerder, H., Gao, M., He, Y., Liu, Y., Luo, Y., Myneni, R. B., ... Janssens, I. A. (2019). Air temperature optima of vegetation productivity across global biomes. *Nature Ecology & Evolution*, 3(5), Article 5.
https://doi.org/10.1038/s41559-019-0838-x

Huete, A. (2016). Vegetation's responses to climate variability. *Nature*, 531(7593), 181-182.
https://doi.org/10.1038/nature17301
Ingrisch, J., & Bahn, M. (2018). Towards a Comparable Quantification of Resilience. *Trends in Ecology & Evolution, 33*(4), 251–259. https://doi.org/10.1016/j.tree.2018.01.013

IPCC. (2022). *Climate Change 2022: Impacts, Adaptation, and Vulnerability*. Contribution of Working Group II to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. https://www.ipcc.ch/report/ar6/wg2/downloads/report/IPCC_AR6_WGII_FinalDraft_FullReport.pdf

IUCN. (2016, August 11). *World Database on Protected Areas*. IUCN. https://www.iucn.org/theme/protected-areas/our-work/world-database-protected-areas

Jakovac, C. C., Peña-Claros, M., Kuyper, T. W., & Bongers, F. (2015). Loss of secondary-forest resilience by land-use intensification in the Amazon. *Journal of Ecology, 103*(1), 67–77. https://doi.org/10.1111/1365-2745.12298

Johnstone, J. F., Allen, C. D., Franklin, J. F., Frelich, L. E., Harvey, B. J., Higuera, P. E., Mack, M. C., Meentemeyer, R. K., Metz, M. R., Perry, G. L., Schoennagel, T., & Turner, M. G. (2016). Changing disturbance regimes, ecological memory, and forest resilience. *Frontiers in Ecology and the Environment, 14*(7), 369–378. https://doi.org/10.1002/fee.1311

Jolliffe, I. T. (2002). *Principal Component Analysis*. Springer-Verlag. https://doi.org/10.1007/b98835

Karger, D. N., Conrad, O., Böhner, J., Kawohl, T., Kreft, H., Soria-Auza, R. W., Zimmermann, N. E., Linder, H. P., & Kessler, M. (2017). Climatologies at high resolution for the Earth’s land surface areas. *Scientific Data, 4*(1), 1–20. https://doi.org/10.1038/sdata.2017.122

Kariuki, R. W., Western, D., Willcock, S., & Marchant, R. (2021). Assessing Interactions between Agriculture, Livestock Grazing and Wildlife Conservation Land Uses: A Historical Example from East Africa. *Land, 10*(1), Article 1. https://doi.org/10.3390/land10010046

Kéfi, S., Dakos, V., Scheffer, M., Van Nes, E. H., & Rietkerk, M. (2013). Early warning signals also precede non-catastrophic transitions. *Oikos, 122*(5), 641–648. https://doi.org/10.1111/j.1600-0706.2012.20838.x

Kéfi, S., Domínguez-García, V., Donohue, I., Fontaine, C., Thébault, E., & Dakos, V. (2019). Advancing our understanding of ecological stability. *Ecology Letters, 22*(9), 1349–1356. https://doi.org/10.1111/ele.13340
Kennedy, C. M., Oakleaf, J. R., Theobald, D. M., Baruch-Mordo, S., & Kiesecker, J. (2019). Managing the middle: A shift in conservation priorities based on the global human modification gradient. *Global Change Biology*, 25(3), 811–826. https://doi.org/10.1111/gcb.14549

Koutsias, N., Arianoutsou, M., Kallimanis, A. S., Mallinis, G., Halley, J. M., & Dimopoulos, P. (2012). Where did the fires burn in Peloponnisos, Greece the summer of 2007? Evidence for a synergy of fuel and weather. *Agricultural and Forest Meteorology*, 156, 41–53. https://doi.org/10.1016/j.agrformet.2011.12.006

Lenihan, J. M., Bachelet, D., Neilson, R. P., & Drapek, R. (2008). Response of vegetation distribution, ecosystem productivity, and fire to climate change scenarios for California. *Climatic Change*, 87(1), 215–230. https://doi.org/10.1007/s10584-007-9362-0

Lippmann, R., Babben, S., Menger, A., Delker, C., & Quint, M. (2019). Development of Wild and Cultivated Plants under Global Warming Conditions. *Current Biology*, 29(24), R1326–R1338. https://doi.org/10.1016/j.cub.2019.10.016

Loudermilk, E. L., O’Brien, J. J., Goodrick, S. L., Linn, R. R., Skowronski, N. S., & Hiers, J. K. (2022). Vegetation’s influence on fire behavior goes beyond just being fuel. *Fire Ecology*, 18(1), 9. https://doi.org/10.1186/s42408-022-00132-9

Lundberg, S., & Lee, S.-I. (2017). A Unified Approach to Interpreting Model Predictions. *arXiv:1705.07874 [Cs, Stat]*. http://arxiv.org/abs/1705.07874

Mahalanobis, P. C. (1936). On the Generalized Distance in Statistics. *Proceedings of the National Institute of Science of India*, 2, 49–55.

Majumder, S., Tamma, K., Ramaswamy, S., & Guttal, V. (2019). Inferring critical thresholds of ecosystem transitions from spatial data. *Ecology*, 100(7), e02722. https://doi.org/10.1002/ecy.2722

McDowell, N. G., Allen, C. D., Anderson-Teixeira, K., Aukema, B. H., Bond-Lamberty, B., Chini, L., Clark, J. S., Dietze, M., Grossiord, C., Hanbury-Brown, A., Hurtt, G. C., Jackson, R. B., Johnson, D. J., Kueppers, L., Lichstein, J. W., Ogle, K., Poulter, B., Pugh, T. A. M., Seidl, R., ... Xu, C. (2020). Pervasive shifts in forest dynamics in a changing world. *Science*, 368(6494), eaaz9463. https://doi.org/10.1126/science.aaz9463
MODIS Atmosphere Science Team. (2015). *MOD06_L2 MODIS/Terra Clouds 5-Min L2 Swath 1km and 5km* [dataset]. NASA Level 1 and Atmosphere Archive and Distribution System. https://doi.org/10.5067/MODIS/MOD06_L2.006

Nolan, R. H., Collins, L., Leigh, A., Ooi, M. K. J., Curran, T. J., Fairman, T. A., Resco de Dios, V., & Bradstock, R. (2021). Limits to post-fire vegetation recovery under climate change. *Plant, Cell & Environment, 44*(11), 3471-3489. https://doi.org/10.1111/pce.14176

Pecl, G. T., Araújo, M. B., Bell, J. D., Blanchard, J., Bonebrake, T. C., Chen, I.-C., Clark, T. D., Colwell, R. K., Danielsen, F., Evengård, B., Falconi, L., Ferrier, S., Frusher, S., Garcia, R. A., Griffis, R. B., Hobday, A. J., Janion-Scheepers, C., Jarzyna, M. A., Jennings, S., ... Williams, S. E. (2017). Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being. *Science, 355*(6332), eaai9214. https://doi.org/10.1126/science.aai9214

Potapov, P., Hansen, M. C., Laestadius, L., Turubanova, S., Yaroshenko, A., Thies, C., Smith, W., Zhuravleva, I., Komarova, A., Minnemeyer, S., & Esipova, E. (2017). The last frontiers of wilderness: Tracking loss of intact forest landscapes from 2000 to 2013. *Science Advances, 3*(1), e1600821. https://doi.org/10.1126/sciadv.1600821

Richter, S., Kipfer, T., Wohlgemuth, T., Calderón Guerrero, C., Ghazoul, J., & Moser, B. (2012). Phenotypic plasticity facilitates resistance to climate change in a highly variable environment. *Oecologia, 169*(1), 269-279. https://doi.org/10.1007/s00442-011-2191-x

Rillig, M. C., Ryo, M., Lehmann, A., Aguilar-Trigueros, C. A., Buchert, S., Wulf, A., Iwasaki, A., Roy, J., & Yang, G. (2019). The role of multiple global change factors in driving soil functions and microbial biodiversity. *Science, 366*(6467), 886-890. https://doi.org/10.1126/science.aay2832

Rito, K. F., Arroyo-Rodríguez, V., Queiroz, R. T., Leal, I. R., & Tabarelli, M. (2017). Precipitation mediates the effect of human disturbance on the Brazilian Caatinga vegetation. *Journal of Ecology, 105*(3), 828-838. https://doi.org/10.1111/1365-2745.12712

Scheffer, M., Bascompte, J., Brock, W. A., Brovkin, V., Carpenter, S. R., Dakos, V., Held, H., van Nes, E. H., Rietkerk, M., & Sugihara, G. (2009). Early-warning signals for critical transitions. *Nature, 461*(7260), Article 7260. https://doi.org/10.1038/nature08227

Sebald, J., Thrippleton, T., Rammer, W., Bugmann, H., & Seidl, R. (2021). Mixing tree species at different spatial scales: The effect of alpha, beta and gamma diversity on disturbance impacts
under climate change. *Journal of Applied Ecology*, 58(8), 1749–1763.
https://doi.org/10.1111/1365-2664.13912

Seddon, A. W. R., Macias-Fauria, M., Long, P. R., Benz, D., & Willis, K. J. (2016). Sensitivity of global terrestrial ecosystems to climate variability. *Nature*, 531(7593), 229–232.
https://doi.org/10.1038/nature16986

Smith, T., & Boers, N. (2023a). Global vegetation resilience linked to water availability and variability. *Nature Communications*, 14(1), Article 1. https://doi.org/10.1038/s41467-023-36207-7

Smith, T., & Boers, N. (2023b). Reliability of vegetation resilience estimates depends on biomass density. *Nature Ecology & Evolution*, 1–10. https://doi.org/10.1038/s41559-023-02194-7

Smith, T., Traxl, D., & Boers, N. (2022). Empirical evidence for recent global shifts in vegetation resilience. *Nature Climate Change*, 1–8. https://doi.org/10.1038/s41558-022-01352-2

Snoek, J., Larochelle, H., & Adams, R. P. (2012). Practical Bayesian Optimization of Machine Learning Algorithms. *arXiv:1206.2944 [Cs, Stat]*. http://arxiv.org/abs/1206.2944

Spasojevic, M. J., Bahlai, C. A., Bradley, B. A., Butterfield, B. J., Tuanmu, M.-N., Sistla, S., Wiederholt, R., & Suding, K. N. (2016). Scaling up the diversity–resilience relationship with trait databases and remote sensing data: The recovery of productivity after wildfire. *Global Change Biology*, 22(4), 1421–1432. https://doi.org/10.1111/gcb.13174

Thompson, I., Mackey, B., McNulty, S., & Mosseler, A. (2009). Forest Resilience, Biodiversity, and Climate Change. *Secretariat of the Convention on Biological Diversity, Montreal. Technical Series No. 43*. 1-67., 43, 1–67.

Tijerín-Triviño, J., Moreno-Fernández, D., Zavala, M. A., Astigarraga, J., & García, M. (2022). Identifying Forest Structural Types along an Aridity Gradient in Peninsular Spain: Integrating Low-Density LiDAR, Forest Inventory, and Aridity Index. *Remote Sensing*, 14(1), 235.
https://doi.org/10.3390/rs14010235

Timpane-Padgham, B. L., Beechie, T., & Klinger, T. (2017). A systematic review of ecological attributes that confer resilience to climate change in environmental restoration. *PLOS ONE*, 12(3), e0173812. https://doi.org/10.1371/journal.pone.0173812
Tomiolo, S., & Ward, D. (2018). Species migrations and range shifts: A synthesis of causes and consequences. *Perspectives in Plant Ecology, Evolution and Systematics, 33*, 62–77.
https://doi.org/10.1016/j.ppees.2018.06.001

Trabucco, A., & Zomer, R. (2019). *Global Aridity Index and Potential Evapotranspiration (ET0) Climate Database v2* [dataset]. figshare. https://doi.org/10.6084/m9.fishshare.7504448.v3

Tuanmu, M.-N., & Jetz, W. (2015). A global, remote sensing-based characterization of terrestrial habitat heterogeneity for biodiversity and ecosystem modelling: Global habitat heterogeneity. *Global Ecology and Biogeography, 24*(11), 1329–1339. https://doi.org/10.1111/geb.12365

UNEP. (1992). *World Atlas of Desertification. First Edition*. United Nations Environment Programme.
https://library.wur.nl/WebQuery/isric/2269522

Van Meerbeek, K., Jucker, T., & Svenning, J.-C. (2021). Unifying the concepts of stability and resilience in ecology. *Journal of Ecology, 109*(9), 3114–3132. https://doi.org/10.1111/1365-2745.13651

Weinans, E., Quax, R., van Nes, E. H., & Leemput, I. A. van de. (2021). Evaluating the performance of multivariate indicators of resilience loss. *Scientific Reports, 11*(1), Article 1. https://doi.org/10.1038/s41598-021-87839-y

Wu, J., & Liang, S. (2020). Assessing Terrestrial Ecosystem Resilience using Satellite Leaf Area Index. *Remote Sensing, 12*(4), Article 4. https://doi.org/10.3390/rs12040595

Yengoh, G. T., Dent, D., Olsson, L., Tengberg, A. E., & Tucker, C. J. (2016). Limits to the Use of NDVI in Land Degradation Assessment. In G. T. Yengoh, D. Dent, L. Olsson, A. E. Tengberg, & C. J. Tucker III (Eds.), *Use of the Normalized Difference Vegetation Index (NDVI) to Assess Land Degradation at Multiple Scales: Current Status, Future Trends, and Practical Considerations* (pp. 27–30). Springer International Publishing. https://doi.org/10.1007/978-3-319-24112-8_4

Yi, C., Hendrey, G., Niu, S., McDowell, N., & Allen, C. D. (2022). Tree mortality in a warming world: Causes, patterns, and implications. *Environmental Research Letters, 17*(3), 030201.
https://doi.org/10.1088/1748-9326/ac507b
Young, D. J. N., Stevens, J. T., Earles, J. M., Moore, J., Ellis, A., Jirka, A. L., & Latimer, A. M. (2017). Long-term climate and competition explain forest mortality patterns under extreme drought. *Ecology Letters*, 20(1), 78–86. https://doi.org/10.1111/ele.12711

Zhang, Q., Yao, T., Huemmrich, K. F., Middleton, E. M., Lyapustin, A., & Wang, Y. (2020). Evaluating impacts of snow, surface water, soil and vegetation on empirical vegetation and snow indices for the Utqiaġvik tundra ecosystem in Alaska with the LVS3 model. *Remote Sensing of Environment*, 240, 111677. https://doi.org/10.1016/j.rse.2020.111677

Zhang, Y., Gentine, P., Luo, X., Lian, X., Liu, Y., Zhou, S., Michalak, A. M., Sun, W., Fisher, J. B., Piao, S., & Keenan, T. F. (2022). Increasing sensitivity of dryland vegetation greenness to precipitation due to rising atmospheric CO2. *Nature Communications*, 13(1), Article 1. https://doi.org/10.1038/s41467-022-32631-3

Zhao, J., Lu, Z., Wang, L., & Jin, B. (2020). Plant Responses to Heat Stress: Physiology, Transcription, Noncoding RNAs, and Epigenetics. *International Journal of Molecular Sciences*, 22(1), 117. https://doi.org/10.3390/ijms22010117
**Supplementary:**

<figure><img src="image_8.png" /><figcaption>**Supplementary Figure 1:** Relationship between variance and mean value for EVI and 3 climate variables (mean temperature, mean water availability, and cloudiness). The left panels (1.) show the residuals for the 2000-2010 period, while the right panels (2.) show the residuals for the 2010-2020 period. We use the residuals from these relationships to calculate sensitivity as in Seddon et al. 2016.</figcaption></figure>
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>**Supplementary Figure 2: Resistance and Recovery changes in the last decade.** Global map of difference in 1) sensitivity to climate (metric opposite to resistance) and 2) EVI-anomaly autocorrelation (as AR1) between 2010-2020 and 2000-2010. Autocorrelation was calculated on consecutive months during the growing season. Red values translate an increase in sensitivity or autocorrelation, related to a decrease in resilience while blue is the opposite. Urban, croplands, wetlands, barren and snow-dominated pixels are mapped in grey. Sensitivity was multiplied by 1.5 for readability purposes. Autocorrelation was multiplied by 100 for readability purposes. Pixel resolution: 0.5° ≈ 55 km at the equator.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

**Supplementary Figure 3: Significant Changes in Sensitivity.** The significance of sensitivity changes was tested through a process similar to a bootstrap. On 100 independent instances, and for each period of interest (2000 to 2010 and 2010 to 2020), months were randomly sampled to re-create time a time-series. Sensitivity was then calculated for each of these re-sampled time-series, to obtain 100 sensitivity scores for each period, and each site of interest. A Student test was then done for each site to determine if sensitivity was higher, or lower, for each period. Here, the t-statistic of the test is plotted if $p < 0.05$ in either of these tests. We reveal 21% of the sites exhibited significant increasing sensitivity, while 22% revealed significant reduced sensitivity, after outlier removal. Pixel resolution: 0.3° ≈ 30 km at the equator.
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

**Supplementary Figure 4:** No link between overall trend in the initial datasets and the processed resilience metrics. To ensure our main results were not due to statistical artifacts and trends in the initial data, we simulated a dataset with random values of EVI and climatic metrics. We simulated the temperature to include a positive trend. Here, we plotted the change in sensitivity (a and b) and in autocorrelation (c and d) resulting from running the analysis on the simulated data and the associated trend in temperature and EVI.
<figure><img src="image_16.png" /></figure>

**Supplementary Figure 5: What is impacting resilience changes?** Importance of each driver in global models, calculated by the proportion of variance lost when that driver is removed from the global models. Land conversion is calculated as the sum of variance explained by each specific type of land conversion, such as forest to grassland. On the left, (a) the graph corresponds to the model describing the change in sensitivity, and on the right (b) the change in autocorrelation.
a)

Sensitivity

1

<figure><img src="image_17.png" /></figure>

2

<figure><img src="image_18.png" /></figure>

3

<figure><img src="image_19.png" /></figure>

4

<figure><img src="image_20.png" /></figure>

Autocorrelation

1

<figure><img src="image_21.png" /></figure>

2

<figure><img src="image_22.png" /></figure>

3

<figure><img src="image_23.png" /></figure>

4

<figure><img src="image_24.png" /></figure>

b)

Sensitivity

1

<figure><img src="image_25.png" /></figure>

2

<figure><img src="image_26.png" /></figure>

3

<figure><img src="image_27.png" /></figure>

4

<figure><img src="image_28.png" /></figure>

Aridity

<figure><img src="image_29.png" /></figure>

SHAP Value

Autocorrelation

1

<figure><img src="image_30.png" /></figure>

2

<figure><img src="image_31.png" /></figure>

3

<figure><img src="image_32.png" /></figure>

4

<figure><img src="image_33.png" /></figure>

Aridity

<figure><img src="image_34.png" /></figure>
<figure><img src="image_35.png" /></figure>

**Supplementary Figure 6: Sensitivity analysis; different sets of optimized parameters don't change the main patterns observed.** Contribution of different variables; (a) Change in Temperature, (b) Aridity, (c) Change in EVI to the change in resilience properties. Here, we calculate aridity as 1-Aridity Index, so that values between 0 and 1 are water-limited. Predicted values for both change in sensitivity (left column) and autocorrelation models, for 4 sets of parameters.
<figure><img src="image_36.png" /><figcaption>**Supplementary Figure 7: Changes in Temperature influence resilience changes (binary classifier models).** Contribution (SHAP values being the predicted change in a prediction, controlling for all other predictors) of change in temperature to the probability of being classified as significantly a) increasing or b) decreasing in Sensitivity. SHAP values are extracted from the binary classifier models.</figcaption></figure>
<figure><img src="image_37.png" /></figure>

**Supplementary Figure 8: Changes in greenness and in temperature interact to influence resilience changes (binary classifier models).** Contribution (SHAP values being the predicted change in a prediction, controlling for all other predictors) of change in EVI (greenness) to the probability of being classified as significantly a) increasing or b) decreasing in Sensitivity. SHAP values are extracted from the binary classifier models.
<figure><img src="image_38.png" /></figure>

**Supplementary Figure 9: How is the resilience of different forest types changing?** Average change in resilience (the labels a) and b) respectively characterize the change in sensitivity and autocorrelation), for different forest types, according to the Resolve classification. Letters indicate a significant difference between forest types (Pairwise permutation tests, with 100 permutations), and a star indicates if the change in resilience is significantly different from 0 (Wilcoxon test, p < 0.05).
<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

**Supplementary Figure 10: What is impacting resilience changes for different vegetation structures?** Importance of each driver, calculated as the proportion of additional variance explained by a given driver to a model. We display the importance for each model focusing on a specific vegetation structure. Changes from one vegetation type to another were included in the model but not displayed in that Figure, as this concerns very few sites, resulting in low importance. The proportion of explained variance were rescaled so that the total without changes in vegetation structure was equal to 1. The label a) corresponds to the model describing the change in sensitivity, and b) the change in autocorrelation.
<figure><img src="image_41.png" /><figcaption>**Supplementary Figure 11: Elevation influences resilience changes.** Contribution (SHAP values being the predicted change in a prediction, controlling for all other predictors) of elevation to the change in resilience properties. Predicted values for both change in sensitivity (a) and autocorrelation (b) models. Autocorrelation and sensitivity are decreasing with greater elevation.</figcaption></figure>
<figure><img src="image_42.png" /></figure>

**Supplementary Figure 12: Change in productivity influences resilience changes in interaction with change in temperature across vegetation structures.** Contribution (SHAP values) of average change in productivity, interacting with the average change in temperature, to the change in resilience properties. Predicted values for both change in sensitivity and autocorrelation models, red meaning gaining more temperature. There is a negative relationship between change in both metrics and EVI. Linear models fit lines show steeper slopes for changes in temperature above 0.2°C. R² for the fitted lines are 0.6 and 0.55 for points experiencing a respective warming greater than 0.2°C and less than 0.2°C.
<figure><img src="image_43.png" /><figcaption>**Supplementary Figure 13: Aridity influences resilience changes across vegetation structures.**</figcaption></figure>

Contribution of aridity to the change in resilience properties. Here, we calculate aridity as 1-Aridity Index, so that values between 0 and 1 are water-limited. Predicted values for both change in sensitivity and autocorrelation models. Autocorrelation is increasing with greater aridity while the trend for change in sensitivity is different for values of aridity above 0 (water limited ecosystems). Since different ecosystem types show different aridity levels, the trends are not the same as the model including all structures for autocorrelation. While savannas and forests seem to exhibit increasing autocorrelation with increasing aridity, grasslands show a similar response for autocorrelation and sensitivity, showing increasing ability to recover for highly arid regions.
<figure><img src="image_44.png" /><figcaption>**Supplementary Figure 14: Change in temperature relationship with change in sensitivity thresholds.** Threshold detection ran with package “chgpts”, on SHAP values for the global model (a), as well as for each vegetation structure (forest (b), grassland (c), savanna (d)).</figcaption></figure>
<figure><img src="image_45.png" /></figure>

<figure><img src="image_46.png" /></figure>

<figure><img src="image_47.png" /></figure>

<figure><img src="image_48.png" /></figure>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>

<figure><img src="image_51.png" /></figure>

<figure><img src="image_52.png" /></figure>

<figure><img src="image_53.png" /></figure>

<figure><img src="image_54.png" /></figure>

<figure><img src="image_55.png" /></figure>

<figure><img src="image_56.png" /><figcaption>**Supplementary Figure 15: Change in temperature relationship with change in autocorrelation thresholds.** Threshold detection ran with package “chgpts”, on SHAP values for the global model (a), as well as for each vegetation structure (forest (b), grassland (c), savanna (d)).</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>**Supplementary Figure 16: Aridity relationship with change in sensitivity thresholds.** Threshold detection ran with package “chgpts”, on SHAP values for the global model (a), as well as for each vegetation structure (forest (b), grassland (c), savanna (d)).</figcaption></figure>
<figure><img src="image_58.png" /></figure>

<figure><img src="image_59.png" /></figure>

<figure><img src="image_60.png" /></figure>

<figure><img src="image_61.png" /></figure>

<figure><img src="image_62.png" /></figure>

<figure><img src="image_63.png" /></figure>

<figure><img src="image_64.png" /></figure>

<figure><img src="image_65.png" /></figure>

<figure><img src="image_66.png" /></figure>

<figure><img src="image_67.png" /></figure>

<figure><img src="image_68.png" /></figure>

<figure><img src="image_69.png" /><figcaption>**Supplementary Figure 17: Aridity relationship with change in autocorrelation thresholds.**</figcaption></figure>

Threshold detection ran with package "chgpts", on SHAP values for the global model (a), as well as for each vegetation structure (forest (b), grassland (c), savanna (d)).
<figure><img src="image_70.png" /><figcaption>**Supplementary Figure 18: Temperature changes in the last decade.** Global map of difference in temperature between 2010-2020 and 2000-2010, from the WorldClim database. Orange pixels show an increase greater than the first threshold related to resilience loss in our model, while red pixels are exhibiting more than double that increase, which is linked to further resilience loss. Dark-grey pixels are related to vegetated areas that are not experiencing such warming. Urban, croplands, wetlands, barren and snow-dominated pixels are mapped in light-grey. Pixel resolution: 0.5° ≈ 55km.</figcaption></figure>
<table><caption>**Supplementary Table 1:** Final static variables in the Random Forest model, citation and DOI link of data source, as well as the reason for inclusion in the model.</caption><thead><tr><th>Type</th><th>Name</th><th>Details</th><th>Citation</th><th>DOI - link</th><th>Reason for inclusion</th></tr></thead><tbody><tr><td rowspan="4">Human impact</td><td>GDP</td><td>Disaggregated GDP</td><td>(Ghosh et al., 2010)</td><td>http://dx.doi.org/10.2174/1874923201003010147</td><td rowspan="4">Human activities can influence the disturbance dynamics as well as the vegetation structure and composition of a system, all of which can affect resilience (Blauhut et al., 2016; Ganteaume et al., 2013; Loudermilk et al., 2022)</td></tr><tr><td>Human Modification</td><td>Cumulative measure of human influence</td><td>(Kennedy et al., 2019)</td><td>https://doi.org/10.1111/gc.b.14549</td></tr><tr><td>Intact Landscapes</td><td>Amount of area that is either an Intact Forest Landscape or a protected area with high IUCN protection status</td><td>(IUCN, 2016; Potapov et al., 2017)</td><td>www.protectedplanet.net (09.03.2022); https://doi.org/10.1126/sciadv.1600821</td></tr><tr><td>Forest Landscape Integrity</td><td>Degree of forest modification for the beginning of 2019</td><td>(Grantham et al., 2020)</td><td>https://doi.org/10.1038/s41467-020-19493-3</td></tr><tr><td>Climate</td><td>Aridity Index</td><td>Ratio between average annual precipitation and potential evapotranspiration</td><td>(Trabucco & Zomer, 2019)</td><td>https://doi.org/10.6084/m9.figshare.7504448.v3</td><td>Aridity has been shown to drive resilience (Smith & Boers, 2023a)</td></tr></tbody></table>
<table><tr><td></td><td>Mean Annual Temperature</td><td>Mean annual daily mean air temperatures averaged over 1 year and for 1981-2010</td><td>(Karger et al., 2017)</td><td>https://doi.org/10.1038/sdata.2017.122</td><td>Temperature drives vegetation productivity, which is related to resilience (Huang et al., 2019)</td></tr><tr><td rowspan="4">Vegetation structure</td><td>Mean Forest Age</td><td>Mean Forest Age for 10% tree cover cutoff</td><td>(Besnard et al., 2021)</td><td>https://doi.org/10.5194/esd-13-4881-2021</td><td>Can affect vegetation structure, which can influence disturbance dynamics (Loudermilk et al., 2022)</td></tr><tr><td>Tree Density</td><td>Stem Density</td><td>(Crowther et al., 2015)</td><td>https://doi.org/10.1038/nature14967</td><td>High density and biomass can increase mortality to disturbances (Koutsias et al., 2012; Young et al., 2017)</td></tr><tr><td>Biomass</td><td>Global Biomass (IPCC Tier 1 Calculated Values) for 2000</td><td>(Gibbs & Ruesch, 2008)</td><td>https://doi.org/10.15485/1463800</td><td></td></tr><tr><td>Texture Dissimilarity</td><td>Dissimilarity (EVI Heterogeneity); Difference in EVI between adjacent pixels</td><td>(Tuanmu & Jetz, 2015)</td><td>https://doi.org/10.1111/geb.12365</td><td>Spatial heterogeneity has been documented to influence resilience (Génin et al., 2018)</td></tr><tr><td></td><td>Burnt Areas Probability</td><td>Probability of Burnt Areas occurrence, based on observations over the 2001-2019 period</td><td>(ESA Land Cover CCI project team; & Defourny, 2016)</td><td>https://cataloque.ceda.ac.uk/uuid/7c114fc6e2884c1f9ca107e7a502f</td><td>Directly related to disturbance dynamics (FAO, 2005)</td></tr></table>
<table><tr><td rowspan="3">Topography and soil</td><td>Elevation</td><td>Elevation (meters)</td><td>(Amatulli et al., 2018)</td><td>https://doi.org/10.1038/sd ata.2018.40</td><td rowspan="3">Site specific control variables, that might also influence species composition and vegetation types (Baldeck et al., 2013)</td></tr><tr><td>Slope</td><td>Slope</td><td>(Amatulli et al., 2018)</td><td>https://doi.org/10.1038/sd ata.2018.40</td></tr><tr><td>Sand Content</td><td>Sand content (50–2000 micro meter) at 0.00m</td><td>(Hengl et al., 2017)</td><td>https://doi.org/10.1371/jo urnal.pone.0 169748</td></tr></table>
<table><caption>**Supplementary Table 2:** Values of each hyperparameter for both the change in sensitivity and autocorrelation models, after Bayesian optimization.</caption><thead><tr><td>Parameter</td><td>Value for change in Sensitivity</td><td>Value for change in Autocorrelation</td></tr></thead><tbody><tr><td>Minimum child weight</td><td>20</td><td>20</td></tr><tr><td>Minimum splitting loss</td><td>0</td><td>30</td></tr><tr><td>Learning rate</td><td>0.1</td><td>0.1</td></tr><tr><td>Maximum depth</td><td>11</td><td>15</td></tr><tr><td>Subsample size</td><td>1</td><td>1</td></tr></tbody></table>