Camille Fournier De Lauriere

Katharina Runge

Gabriel Smith

Vasilis Dakos

Sonia Kefi

Thomas Crowther

Miguel Berdugo

Uncovering

Sonia Kéfi

Uncovering global drivers threatening vegetation resilience

Keywords: resistance, sensitivity, autocorrelation, recovery, climate change, degradation, aridity

à 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.

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 [START_REF] Johnstone | Changing disturbance regimes, ecological memory, and forest resilience[END_REF][START_REF] Mcdowell | Pervasive shifts in forest dynamics in a changing world[END_REF]. Quantifying those changes in ecological resilience across the globe remains a major conceptual challenge in ecology [START_REF] Dakos | Ecological resilience: What to measure and how[END_REF]. 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 heterogenous patterns [START_REF] Boulton | Pronounced loss of Amazon rainforest resilience since the early 2000s[END_REF][START_REF] Feng | Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale[END_REF][START_REF] Forzieri | Emerging signals of declining forest resilience under climate change[END_REF][START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF][START_REF] Zhang | Increasing sensitivity of dryland vegetation greenness to precipitation due to rising atmospheric CO2[END_REF]. However, there is little known about the mechanisms that drive these patterns [START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF]. 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 [START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF][START_REF] Timpane-Padgham | A systematic review of ecological attributes that confer resilience to climate change in environmental restoration[END_REF], 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 [START_REF] Boulton | Pronounced loss of Amazon rainforest resilience since the early 2000s[END_REF][START_REF] Feng | Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale[END_REF][START_REF] Forzieri | Emerging signals of declining forest resilience under climate change[END_REF][START_REF] Jakovac | Loss of secondary-forest resilience by land-use intensification in the Amazon[END_REF][START_REF] Zhang | Increasing sensitivity of dryland vegetation greenness to precipitation due to rising atmospheric CO2[END_REF]. For example, while warming is considered the main global driver in altering ecosystem resilience [START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF], its effect on other ecosystem properties has been shown to vary across ecosystems [START_REF] Berner | Satellite observations document trends consistent with a boreal forest biome shift[END_REF][START_REF] De Keersmaecker | A model quantifying global vegetation resistance and resilience to short-term climate anomalies and their relationship with vegetation cover[END_REF][START_REF] Pecl | Biodiversity redistribution under climate change: Impacts on ecosystems and human well-being[END_REF][START_REF] Tomiolo | Species migrations and range shifts: A synthesis of causes and consequences[END_REF]. 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 [START_REF] Rillig | The role of multiple global change factors in driving soil functions and microbial biodiversity[END_REF] and non-linear [START_REF] Berdugo | Prevalence and drivers of abrupt vegetation shifts in global drylands[END_REF][START_REF] García-Palacios | Climate mediates the biodiversityecosystem stability relationship globally[END_REF] 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 [START_REF] Kéfi | Advancing our understanding of ecological stability[END_REF][START_REF] Van Meerbeek | Unifying the concepts of stability and resilience in ecology[END_REF], the context dependency of resilience change, the interactions between environmental drivers, and the expected nonlinearity 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 [START_REF] Domínguez-García | Unveiling dimensions of stability in complex ecological networks[END_REF][START_REF] Kéfi | Early warning signals also precede non-catastrophic transitions[END_REF] and empirical studies [START_REF] Eby | Alternative stable states and spatial indicators of critical slowing down along a spatial gradient in a savanna ecosystem[END_REF][START_REF] Majumder | Inferring critical thresholds of ecosystem transitions from spatial data[END_REF] 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 [START_REF] Capdevila | Reconciling resilience across ecological systems, species and subdisciplines[END_REF][START_REF] Hodgson | What do you mean, 'resilient'?[END_REF][START_REF] Ingrisch | Towards a Comparable Quantification of Resilience[END_REF][START_REF] Weinans | Evaluating the performance of multivariate indicators of resilience loss[END_REF].

We measured time series autocorrelation as a proxy of recovery [START_REF] Forzieri | Emerging signals of declining forest resilience under climate change[END_REF][START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF] and sensitivity as a proxy of resistance [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF].

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)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020), relative to the first half (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(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 [START_REF] Huete | Vegetation's responses to climate variability[END_REF][START_REF] Zhang | 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[END_REF]. 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 [START_REF] Abatzoglou | TerraClimate, a highresolution global dataset of monthly climate and climatic water balance from 1958-2015[END_REF].

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 [START_REF] Didan | MOD13C2 MODIS/Terra Vegetation Indices Monthly L3 Global 0.05Deg CMG V006[END_REF], 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 [START_REF] Abatzoglou | TerraClimate, a highresolution global dataset of monthly climate and climatic water balance from 1958-2015[END_REF], 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 [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF]. Cloudfraction 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 [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF], 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 [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF].

The sensitivity index for each climate variable was calculated using eq. (1), where 𝑣𝑎𝑟(𝐸𝑉𝐼) and 𝑣𝑎𝑟(𝑐𝑙𝑖𝑚) are the corrected variance for EVI and a given climate variable (temperature, water availability and cloudiness), and 𝛽 𝑐𝑙𝑖𝑚 is a weight that measures how closely related are anomalies in EVI and in that climate variable.

𝑆𝑒𝑛𝑠𝑖𝑡𝑖𝑣𝑖𝑡𝑦 𝑐𝑙𝑖𝑚 = ( 𝑣𝑎𝑟(𝐸𝑉𝐼) 𝑣𝑎𝑟(𝑐𝑙𝑖𝑚) * 𝛽 𝑐𝑙𝑖𝑚 ) (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 𝛽 𝑐𝑙𝑖𝑚 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 [START_REF] Jolliffe | Principal Component Analysis[END_REF], 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, (𝛽 𝑐𝑙𝑖𝑚 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 [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF]. 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 [START_REF] Boulton | Pronounced loss of Amazon rainforest resilience since the early 2000s[END_REF][START_REF] Forzieri | Emerging signals of declining forest resilience under climate change[END_REF][START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF]. 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 and3). 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 [START_REF] Forzieri | Emerging signals of declining forest resilience under climate change[END_REF][START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF].

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 [START_REF] Dinerstein | An Ecoregion-Based Approach to Protecting Half the Terrestrial Realm[END_REF]. 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 [START_REF] Chen | XGBoost: A Scalable Tree Boosting System[END_REF]. 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 [START_REF] Addicott | Toward an improved understanding of causation in the ecological sciences[END_REF][START_REF] Berdugo | Prevalence and drivers of abrupt vegetation shifts in global drylands[END_REF], 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 [START_REF] Burrell | Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification[END_REF][START_REF] Eckert | Trend analysis of MODIS NDVI time series for detecting land degradation and regeneration in Mongolia[END_REF], however, it is important to recognize that this relationship is complex and context-dependent [START_REF] Easdale | Phases or regimes? Revisiting NDVI trends as proxies for land degradation[END_REF][START_REF] Yengoh | Limits to the Use of NDVI in Land Degradation Assessment[END_REF]. 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 [START_REF] Fick | WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas[END_REF]. 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 [START_REF] Friedl | MCD12Q1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V006[END_REF]. 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 [START_REF] Feng | Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale[END_REF] and Aridity has been widely used to characterise various ecosystems [START_REF] Berdugo | Global ecosystem thresholds driven by aridity[END_REF][START_REF] Tijerín-Triviño | Identifying Forest Structural Types along an Aridity Gradient in Peninsular Spain: Integrating Low-Density LiDAR, Forest Inventory, and Aridity Index[END_REF]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 [START_REF] Snoek | Practical Bayesian Optimization of Machine Learning Algorithms[END_REF] 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 [START_REF] Lundberg | A Unified Approach to Interpreting Model Predictions[END_REF]. 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 and8. 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" [START_REF] Fong | chngpt: Threshold regression model estimation and inference[END_REF] 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 [START_REF] Berdugo | Prevalence and drivers of abrupt vegetation shifts in global drylands[END_REF][START_REF] Mahalanobis | On the Generalized Distance in Statistics[END_REF]. 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 and3). 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. 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.

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 2 =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, cloudcover, 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). andd) 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. 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.

Discussion

Changes in resilience show heterogenous patterns across the globe, in accordance to previous work [START_REF] Forzieri | Emerging signals of declining forest resilience under climate change[END_REF][START_REF] Smith | Empirical evidence for recent global shifts in vegetation resilience[END_REF], 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 [START_REF] Smith | Reliability of vegetation resilience estimates depends on biomass density[END_REF]. 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 [START_REF] Dale | 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[END_REF][START_REF] De Groot | Climate change impacts on future boreal fire regimes[END_REF][START_REF] Lenihan | Response of vegetation distribution, ecosystem productivity, and fire to climate change scenarios for California[END_REF] sometimes coupled with human activity [START_REF] Gou | Intra-annual relationship between precipitation and forest disturbance in the African rainforest[END_REF][START_REF] Rito | Precipitation mediates the effect of human disturbance on the Brazilian Caatinga vegetation[END_REF]. 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 [START_REF] Beckage | Fire Feedbacks with Vegetation and Alternative Stable States[END_REF][START_REF] Johnstone | Changing disturbance regimes, ecological memory, and forest resilience[END_REF][START_REF] Scheffer | Early-warning signals for critical transitions[END_REF]. 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 [START_REF] Zhao | Plant Responses to Heat Stress: Physiology, Transcription, Noncoding RNAs, and Epigenetics[END_REF], leading to higher sensitivity [START_REF] Yi | Tree mortality in a warming world: Causes, patterns, and implications[END_REF] and slower recruitment [START_REF] Hatfield | Temperature extremes: Effect on plant growth and development[END_REF][START_REF] Lippmann | Development of Wild and Cultivated Plants under Global Warming Conditions[END_REF]. 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 [START_REF] Nolan | Limits to post-fire vegetation recovery under climate change[END_REF].

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 [START_REF] Berner | Satellite observations document trends consistent with a boreal forest biome shift[END_REF][START_REF] Boulton | Pronounced loss of Amazon rainforest resilience since the early 2000s[END_REF][START_REF] Sebald | Mixing tree species at different spatial scales: The effect of alpha, beta and gamma diversity on disturbance impacts under climate change[END_REF][START_REF] Spasojevic | Scaling up the diversity-resilience relationship with trait databases and remote sensing data: The recovery of productivity after wildfire[END_REF]. However, we provide the first global empirical evidence showing a clear effect of losing EVI, or productivity (potentially linked to land degradation, [START_REF] Burrell | Anthropogenic climate change has driven over 5 million km2 of drylands towards desertification[END_REF][START_REF] Eckert | Trend analysis of MODIS NDVI time series for detecting land degradation and regeneration in Mongolia[END_REF] 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 [START_REF] Thompson | Forest Resilience, Biodiversity, and Climate Change. Secretariat of the Convention on Biological Diversity[END_REF]. In other words, when human activities lead to the degradation of ecosystems, we do not only lose the productivity of those systems [START_REF] Barlow | Anthropogenic disturbance in tropical forests can double biodiversity loss from deforestation[END_REF][START_REF] Foley | Global Consequences of Land Use[END_REF][START_REF] Heinrichs | Habitat degradation and loss as key drivers of regional population extinction[END_REF], 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 [START_REF] Rillig | The role of multiple global change factors in driving soil functions and microbial biodiversity[END_REF]. On the other hand, rising CO2 concentration and temperature may increase vegetation productivity [START_REF] Dusenge | Plant carbon metabolism and climate change: Elevated CO2 and temperature impacts on photosynthesis, photorespiration and respiration[END_REF][START_REF] Gustafson | Do rising temperatures always increase forest productivity? Interacting effects of temperature, precipitation, cloudiness and soil texture on tree species growth and competition[END_REF], 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 [START_REF] Bentz | Climate Change and Bark Beetles of the Western United States and Canada: Direct and Indirect Effects[END_REF][START_REF] Halofsky | Changing wildfire, changing forests: The effects of climate change on fire regimes and vegetation in the Pacific Northwest, USA[END_REF][START_REF] Lenihan | Response of vegetation distribution, ecosystem productivity, and fire to climate change scenarios for California[END_REF]. 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 [START_REF] Choat | Triggers of tree mortality under drought[END_REF] 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 [START_REF] Kariuki | Assessing Interactions between Agriculture, Livestock Grazing and Wildlife Conservation Land Uses: A Historical Example from East Africa[END_REF], inducing fire regimes or invasions [START_REF] Archibald | Managing the human component of fire regimes: Lessons from Africa[END_REF][START_REF] Gorgone-Barbosa | Disturbance as a factor in breaking dormancy and enhancing invasiveness of African grasses in a Neotropical Savanna[END_REF], but could also be misidentified degraded forests (exhibiting low percent tree cover) by the MODIS land classification [START_REF] Friedl | MCD12Q1 MODIS/Terra+Aqua Land Cover Type Yearly L3 Global 500m SIN Grid V006[END_REF]. As such, we emphasise that unique sitespecific environmental and biological properties must always be considered when assessing the resilience of a system [START_REF] Dakos | Ecological resilience: What to measure and how[END_REF].

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 [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF][START_REF] Wu | Assessing Terrestrial Ecosystem Resilience using Satellite Leaf Area Index[END_REF]. 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 [START_REF] Bachofen | Phenotypic plasticity versus ecotypic differentiation under recurrent summer drought in two droughttolerant pine species[END_REF][START_REF] Richter | Phenotypic plasticity facilitates resistance to climate change in a highly variable environment[END_REF]. 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 [START_REF] Berdugo | Species-specific adaptations determine how aridity and biotic interactions drive the assembly of dryland plant communities[END_REF][START_REF] Berdugo | Global ecosystem thresholds driven by aridity[END_REF]. Vegetation adapted to aridity levels drier than 0.5 is very specialised at the organism level, probably including strategies to evade environmental stress [START_REF] Berdugo | Ecological mechanisms underlying aridity thresholds in global drylands[END_REF]. 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 [START_REF] Smith | Reliability of vegetation resilience estimates depends on biomass density[END_REF] 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.

Figure 1 :

 1 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.

Figure 2 :

 2 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

Figure 3 :

 3 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).

Figure 4 :

 4 Figure 4: Change in productivity interacts with change in temperature to influence resilience. 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.

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.

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/"

Conflict of Interest Statement:

The authors do not declare any conflict of interested.

Supplementary:

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 [START_REF] Seddon | Sensitivity of global terrestrial ecosystems to climate variability[END_REF]. 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. 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 timeseries. 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.

Supplementary Figure

Supplementary