Prevalence and drivers of abrupt vegetation shifts in global drylands

Prevalence and drivers of abrupt vegetation shifts in global drylands Miguel Berdugo a,b,c,1,2 , Juan J. Gait an d,e,f , Manuel Delgado-Baquerizo g,h , Thomas W. Crowther a , and Vasilis Dakos i,2 Edited by Pablo Marquet, Pontificia Universidad Catolica de Chile, Santiago, Chile; received January 2, 2022; accepted September 9, 2022 The constant provision of plant productivity is integral to supporting the liability of ecosystems and human wellbeing in global drylands. Drylands are paradigmatic examples of systems prone to experiencing abrupt changes in their functioning. Indeed, space-fortime substitution approaches suggest that abrupt changes in plant productivity are widespread, but this evidence is less clear using observational time series or experimental data at a large scale. Studying the prevalence and, most importantly, the unknown drivers of abrupt (rather than gradual) dynamical patterns in drylands may help to unveil hotspots of current and future dynamical instabilities in drylands. Using a 20-y global satellite-derived temporal assessment of dryland Normalized Difference Vegetation Index (NDVI), we show that 50% of all dryland ecosystems exhibiting gains or losses of NDVI are characterized by abrupt positive/negative temporal dynamics. We further show that abrupt changes are more common among negative than positive NDVI trends and can be found in global regions suffering recent droughts, particularly around critical aridity thresholds. Positive abrupt dynamics are found most in ecosystems with low seasonal variability or high aridity. Our work unveils the high importance of climate variability on triggering abrupt shifts in vegetation and it provides missing evidence of increasing abruptness in systems intensively managed by humans, with low soil organic carbon contents, or around specificaridity thresholds. These results highlight that abrupt changes in dryland dynamics are very common, especially for productivity losses, pinpoint global hotspots of dryland vulnerability, and identify drivers that could be targeted for effective dryland management. abrupt shifts jproductivity dynamics jdryland ecology jremote sensing Human activities are exerting an unprecedented level of stress on terrestrial ecosystems worldwide. Drylands, covering over half of the terrestrial surface (1), are among the most vulnerable terrestrial biomes in the face of desertification; that is a permanent decline in ecosystem productivity associated with climate change and land degradation (2). Space-for-time substitution studies suggested that declines in dryland productivity occur in an abrupt manner along aridity gradients (3) and with land use change (4). However, whether and why temporal changes in dryland productivity also occur in an abrupt fashion is largely unknown, limiting our capacity to predict whether human environmental pressures on dryland ecosystems will result in dysfunctional ecosystems, and compromise human livelihood, in a drier and more degraded world. Drastic and abrupt changes in ecosystem productivity that are maintained through time indicate the occurrence of a profound and long-lasting change in an ecosystem’s structural or functional feature (5). In drylands, abrupt shifts have been predominantly studied theoretically or using space-for-time substitutions (6), but empirical evidence for shifts over time of dryland vegetation is more scarce at global scales. Indeed, several studies on the shape of vegetation dynamics have occurred in the last years using breakpoint detection methods in time series of indices related to productivity (such as the Normalized Difference Vegetation Index, NDVI). Mostly focused on changes in the slope (the trend) of productivity trajectories (7–9), these studies have proven of high interest for identifying turning points, nuancing our understanding of long-term trends of productivity in drylands especially focused on rainfall efficiency (10, 11). Yet, research documenting the occurrence of abrupt shifts, defined as a level shift in the intercept in productivity that is persistent once it occurs (i.e., sudden step changes in the intercept, rather than the slope, that change abruptly NDVI values and do not revert immediately) remains scarce at a global scale. Because an accumulation of such abrupt shifts can be considered an indication of unstable ecosystems potentially shifting between alternative states (12), studying the occurrence of hotspots of abrupt shifts and identifying their drivers is an important step to understand under which conditions an ecosystem is most vulnerable. Indeed, identifying which environmental factors determine the occurrence of dryland abrupt shifts (e.g., climate, soil characteristics, or human activities) is key to understanding the mechanisms for the emergence of regime shifts, and thus for developing prevention Significance After decades of analyzing productivity patterns worldwide, we still ignore the prevalence and drivers of abrupt (rather than gradual) dynamical patterns of productivity in drylands. Here, we unveil that the global prevalence of hotspots of abrupt shifts of productivity in drylands is very high, especially for negative abrupt changes, and we relate abrupt dynamical patterns with climatic, topographical, edaphic, and human factors. Our work contributes to evaluating dryland sensitivity to abrupt losses in productivity in the coming years. Author affiliations: a Institute of Integrative Biology, Department of Environment Systems Science, ETH Z€ urich, 8092 Z€ urich, Switzerland; b Institut de Biologia Evolutiva (CSIC-UPF), Passeig Maritim de la Barceloneta 37, 8003 Barcelona, Spain; c ICREA-Complex Systems labDepartment of Biology, UPF-PRBB, Doctor Aiguader 88Barcelona, 8003 Spain; d Instituto Nacional de Tecnolog ıa Agropecuaria (INTA), Instituto de Suelos-CNIA, 1686 Buenos Aires, Argentina; e Departamento de Tecnolog ıa, Universidad Nacional de Luj an, 6700 Luj an, Buenos Aires, Argentina; f Consejo Nacional de Investigaciones de Argentina (CONICET), Buenos Aires, Argentina; g Unidad asociada CSICUPO (BioFun), Universidad Pablo de Olavide, Sevilla, 41013 Spain; h Laboratorio de Biodiversidad y Funcionamiento EcosistemicoInstituto de Recursos Naturales y Agrobiologia de Sevilla, Av. Reina Mercedes 10Sevilla, 41012, Spain; and i Institut des Sciences de l’Evolution de Montpellier (ISEM), CNRS, University of Montpellier, EPHE, IRD, 34090 Montpellier, France Author contributions: M.B. and V.D. designed research; M.B. and V.D. performed research; J.J.G., M.D.-B., and T.W.C. contributed new reagents/analytic tools; M.B., J.J.G., and V.D. analyzed data; and M.B., M.D.-B., T.W.C., and V.D. wrote the paper. The authors declare no competing interest. This article is a PNAS Direct Submission. Copyright © 2022 the Author(s). Published by PNAS. This article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND). 1 To whom correspondence may be addressed. Email: [email protected]. 2 M.B. and V.D. contributed equally to this work. This article contains supporting information online at http://www.pnas.org/lookup/suppl/doi:10.1073/pnas. 2123393119/-/DCSupplemental. Published October 17, 2022. PNAS 2022 Vol. 119 No. 43 e2123393119 https://doi.org/10.1073/pnas.2123393119 1of10 RESEARCH ARTICLE | ECOLOGY Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. efforts against dryland abrupt desertification. Among the several complex mechanisms that have been described as generators of abrupt dynamics, positive feedbacks (changes that reinforce themselves due to interactions between different ecosystem attributes) are suggested to be common in dryland ecosystems (13). Previous empirical studies have suggested that dryland productivity should respond abruptly more often in places where aridity is close to a specific threshold [∼0.8, where abrupt decreases of plant cover are described through spatial gradients (3)], where soils are less fertile [due to the link between soil hydrology and fertility in drylands, (14)], or where human activities such as grazing are intense [due to the amplification of vegetation destruction through plant-plant positive feedbacks, (4)]. Yet, despite the strong conceptual basis for such expectations, the drivers of abrupt shifts in drylands have never been investigated using temporal series of vegetation productivity at global scales and remain largely undetermined. Here, we analyze the temporal dynamics in ecosystem NDVI for more than 40,000 natural dryland ecosystems worldwide, to determine whether abrupt changes are common in drylands, and to identify which environmental factors explain the incidence of abrupt temporal dynamics in plant productivity. Specifically, we compiled a 20-y satellite-derived temporal assessment of ecosystem productivity using NDVI from the Moderate Resolution Imaging Spectroradiometer (MODIS) satellite. NDVI is used as a surrogate of greenness or plant fractional cover (15), which is tightly linked to ecosystem productivity (16), notably in drylands (17), where it also relates with overall ecosystem functioning (18) and is a common proxy used to determine ecosystem degradation and desertification (19). We first classified temporal trajectories of NDVI in three classes attending to their overall trend (neutral, positive, and negative). Neutral dynamics represent average lack of changes in NDVI over time and can include cancelling positive and negative changes in NDVI. We then quantified the percentage of ecosystems following linear, curvilinear, or abrupt patterns within positive and negative trends. By doing this, we aim to classify the most common types of ecosystem dynamics in drylands worldwide. We argue that, rather than what has been applied in previous studies about turning points, identification of abrupt and persistent changes in NDVI may allow further insights on the prevalence of abruptness in temporal trends of vegetation, a key feature linked to ecosystem instability. We then conducted an attribution analysis of NDVI time series to investigate potential environmental drivers of NDVI (i.e., linear relationships between NDVI time series and temperature, water balance, and CO 2 time series). Finally, we investigated with machine learning models the associationbetweentheoccurrenceofabruptpositiveandnegative changes with multiple site-specific environmental factors including historical climatic conditions (e.g., aridity, precipitation variability), descriptors of climate change (e.g., trends in precipitation, number of drought years), environmental characteristics (e.g., soil properties, topography), and human impact (e.g., human footprint, ruminants density). Results and Discussion Our analyses revealed that almost half of all ecosystems in global drylands have been subject to (positive and negative) changes in NDVI over the last 20 y (44% of all studied ecosystems, Fig. 1B). In the rest of the studied ecosystems, we found neutral NDVI trajectories suggesting that many ecosystems in drylands have maintained an overall constant supply of plant productivity over the last two decades. Lack of sustained changes in temporal productivity in global drylands are probably associated with the typically reported slow responses of dryland vegetation to climatic temporal fluctuations. The occurrence of neutral, positive, or negative dynamics was associated with the trend in water availability observed for each ecosystem during the study period (SI Appendix, Fig. S1), although higher abundance of positive rather than negative changes was probably also indicating a positive CO 2 fertilization effect (20). This is indeed suggested by a higher influence of CO 2 dynamics on NDVI for positive trends in respect to the other trajectory types (SI Appendix,Figs.S2andS3). Our analyses further show that, when significant temporal changes in NDVI are reported, these dynamics mostly follow abrupt trajectories (Fig. 1B), suggesting that nonlinear drastic changes in vegetation have occurred in many drylands worldwide over the last two decades and they are becoming increasingly more common in recent years (SI Appendix, Fig. S4). These results from long-term time-series observations are in partial agreement with previous studies on shifting trends and abrupt turning points, revealing the prevalence of abrupt shifts in global drylands (7–11). However, it should be noted that our methodology (more focused on abrupt level shifts than on trend changes) and the timespan in which we focused do not allow a complete comparison with other studies (SI Appendix, Fig. S5). The timespan of our analyses is limited by data availability to 20 y (we restricted our analyses to MODIS data, see Materials and Methods). As such, as the abrupt changes identified here are generally maintained for at least 5 y or more, we cannot rule out the possibility that they represent transient changes that could be reversed in the near future. Abrupt changes were found in both positive and negative temporal ecosystem trajectories. However, our analyses revealed that the prevalence in abrupt changes is greater in the case of negative dynamics (i.e., in 2 out of 3 ecosystems, Fig. 1B). Abrupt losses of dryland productivity have been attributed to aridity increase, overgrazing, or soil degradation (3, 4, 6); here we demonstrate that the prevalence of such behavior in global drylands is dominant when negative changes occur. On the other hand, we also found a high prevalence of abrupt positive changes (Fig. 1B). Positive abrupt changes in vegetation or productivity are anticipated as a result of changes in land use such as land abandonment (10, 21). However, the interaction of different factors may be also responsible for the high proportion found. For instance, the CO 2 fertilization may interact with higher water availability in “good”years yielding abrupt rather than gradual increases in productivity. Such interactions, however, are not emerging clearly from the attribution analyses (SI Appendix, Figs. S2 and S3), which shows more importance of water availability in positive and negative abrupt changes compared to the rest, but also a higher amount of unexplained variance. This higher amount of unexplained variance for abrupt changes compared to the other types suggests that abrupt changes are intrinsic and they may be better explained by site-specific characteristics or by the occurrence of specific events rather than overall dynamics of environmental drivers. By mapping the global distribution of abrupt negative and positive dynamics (Fig. 1Aand SI Appendix, Figs. S6 and S7), we identified some important hotspots showing abrupt shifts in NDVI. For example, semiarid forests have gone through important drastic reductions in NDVI over the last two decades. Some of them are historically consistent, e.g., in South America, in the Caatinga region, (Fig. 1A), where droughts in the 2012–2015 period have been reported to drastically reduce NDVI (22). Further analyses show that abrupt negative types are associated with several changes in land structure such as desertification (transformation of vegetated land into deserts), forest loss (transformations of forests into open ecosystems or 2of10 https://doi.org/10.1073/pnas.2123393119 pnas.org Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. into savannahs), and savannah degradation into open ecosystems (SI Appendix, Fig. S8). In contrast, we found regions from central Asia experiencing abrupt regreening during the 21st century (Fig. 1A). Such regreening has been previously associated with land abandonment (21) and, more recently, to afforestation practices carried out in central China (23) but they could also be related to contractions of the Gobi Desert in recent years (24). Indeed, we find that changes from desert systems into vegetated land, from open ecosystems into savannahs, and from savannahs into forest are associated with positive abrupt shifts of NDVI (SI Appendix, Fig. S8). Moreover, abrupt dynamical types (especially negative abrupt changes) negatively associate with unchanged landscapes, strongly suggesting that is unlikely to experience abrupt changes of NDVI and still maintain the same ecosystem structure. This result is even more intriguing if taking into account that the landscape units used in these analyses have been previously suggested to be alternative stable states of terrestrial ecosystems (25). Further analyses should be conducted in situ to disentangle fine-scale plant community changes that can also trigger abrupt shifts in NDVI. Assessing the environmental factors associated with abrupt shifts in plant productivity is essential to minimize the negative impacts and maximize the positive influence of these drastic dynamics in the maintenance of human wellbeing and ecosystems. To understand if there is any underlying pattern driving the probability to experience abrupt shifts in NDVI, we used a machine learning Random Forest (RF) algorithm to identify which site-specificclimatic, soil, human, and topographic factors were associated with hotspots of abrupt positive and negative shapes (vs. curvilinear and linear types of the same trend). Our models exhibited moderate to high accuracy on classifying either positive abrupt shifts Curvilinear Abrupt (50.62%) (67.19%) (15.42%) (17.38%) (83.08%) (16.92%) (35.92%) Linear Abrupt Curvilinear Linear Linear Curvilinear Neutral Negative Positive (36%) (56%) (8%) B 2005 2010 2015 2005 2010 2015 2005 2010 2015 2005 2010 2015 2005 2010 2015 2005 2010 2015 2005 2010 2015 2005 2010 2015 IV DN IV D N IVDN I V DN I VD N IVD N NDVI NDVI A Abrupt Abrupt Curvilinear Curvilinear Curvilinear Linear Linear Linear Negative non abrupt (curvilinear or linear) Negative abrupt Positive abrupt Positive non abrupt Neutral (curvilinear or linear) (13.46%) Fig. 1. Abrupt changes in ecosystem productivity dynamics govern global drylands. Map on classified annual trajectories of NDVI based on their trend and shape estimated for 41,830 dryland sites across the globe monitored between 2000–2019 (A). Panel Bshows a pie plot presenting the percentage of neutral (yellow), positive (blue), and negative (red) NDVI trajectory trends. The ring around the pie plot indicates the percentage of linear, quadratic (curvilinear), and abrupt trajectory shapes for each trajectory trend. Panels around the pie show an example time series of each classified trajectory shape (y axis: standardized NDVI units; x axis: years, 2000–2019, dashed line indicates type of fitted trajectory). PNAS 2022 Vol. 119 No. 43 e2123393119 https://doi.org/10.1073/pnas.2123393119 3of10 Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. (62.7%) or negative abrupt shifts (70.4%, SI Appendix,TableS1). Although the model accuracy could have improved had we information on unknown drivers, our analysis succeeded in identifying the most important environmental factors explaining abrupt positive and negative changes in NDVI in global drylands (Fig. 2). Some of these environmental factors are important for both positive and negative abrupt changes. We show, for instance, a high relevance of climate variability for vegetation dynamics (Fig. 3A), as it is often argued in drylands (26, 27). However, by showing that higher rainfall interannual variability is associated with dryland vulnerability to abrupt shifts (both for positive and negative changes), our results suggest that the increasing climate variability such as that expected with climate change may enhance the frequency and severity of abrupt vegetation shifts. This is something only hypothesized by few model studies (26), usually disregarded on studies about thresholds (6), and with high importance for certain phenomena in certain regions such as in South America due to variations in the strength of El Ni~no-Southern oscillation (28). Rainfall seasonality and temperature seasonality were also important predictors of positive and negative abrupt changes in NDVI, but their effects were opposite between negative and positive changes (Fig. 3 Band C). Ecosystems with greater rainfall seasonality were associated with negative abrupt shifts, while the opposite pattern was found for positive abrupt shifts (Fig. 3B). On the contrary, ecosystems with more constant climates (less seasonality in rainfall and temperature) supported abrupt positive increases in NDVI (Fig. 3 Band C), suggesting that these ecosystems are more efficient in capitalizing rainfall increases (29). This may have important implications to understand the impacts of climatic seasonality on ecosystem productivity (30), e.g., monsoonal climates in India (31). Factors associated with climatic changes in water balance such as trends in precipitation and number of drought years in the study period were important drivers of negative abrupt changes, indicating clearly that droughts are important triggers of negative abrupt changes in NDVI (Fig. 2 and SI Appendix,Fig.S9). Although of lower statistical importance for both positive and negative abrupt changes (Fig. 2), sites with high soil organic carbon exhibited less frequent abrupt negative changes, and from a certain amount of soil organic carbon the occurrence of positive abrupt shifts sharply increased (Fig. 3G). This could be due to the role of soil carbon on modulating soil water availability feedbacks in drylands (14). Regarding aridity, we identified a very contrasting behavior between positive and negative abrupt changes. The probability of experiencing negative abrupt shifts declined above an aridity level ∼0.8—the transition from semiarid to arid ecosystems— and peaked just before this value. In contrast, the probability of a positive abrupt shift increased above the same value (∼0.8, Fig. 3D). This means that the aridity threshold of ∼0.8 is separating contrastingly different dynamical zones. In sites with aridity around 0.8, ecosystems suffering droughts or sustained losses of water incomes (that is, exhibiting negative trends in NDVI) are highly vulnerable to experiencing abrupt negative loss in NDVI, whereas in sites with higher than 0.8 aridity values, ecosystems experiencing wetting trends (that is, positive trends in NDVI) are prone to exhibit abrupt greening. This is consistent with the existence of a marked aridity threshold at values ∼0.8 affecting vegetation development as evidenced in previous studies using spatial gradients (3). This result also suggests that sites with aridity around 0.8 will be highly variable; able to show rapid decreases or increases in NDVI depending on interannual climatic variations. Ecosystems supporting these levels of aridity will be particularly challenging to manage. In fact, our results can help to explain abrupt changes in dryland ecosystems in the past. For example, we could explain why rapid desertification that occurred in the Sahel in the 1970s [a zone with aridity values around 0.8, whose drought in the 70s actually gave birth to the term of desertification, (32)] recovered so fast after wet years during the 1980s (33). Evidently, we cannot conclude that crossing aridity values of 0.8 a given year would lead to abrupt shifts in productivity, which should be studied by examining the relationship between dynamical changes in aridity and productivity taking into account both how strong is the change in aridity and its duration. However, our results do show an accumulation of negative abrupt shifts in NDVI before this value and an accumulation of positive abrupt change after, which matches expectations and expands the implications of the existence of aridity thresholds into vegetation dynamics. We also found environmental factors specifically relevant for explaining positive abrupt changes in NDVI. For example, when mean temperature of the driest quarter was below 0 (freezing conditions), the probability of abrupt positive changes decreased, whereas it increased when the mean temperature was Population density Depth to bedrock Ruminants density Sand content Soil organic carbon Elevation Solar radiation Productivity (NDVI) Aridity Depth of water table Trend in Precipitation Rainfall of coldest quarter Mean temperature of driest quarter Human footprint Diurnal temperature range Precipitation variability Mean Temperature of wettest quarter Rainfall seasonality Trend in PDSI Years of drought Temperature seasonality 0.00 0.05 0.10 0.15 Importance Negative abrupt changes Population density Depth to bedrock Years of drought Sand content Elevation Soil organic carbon Human footprint Depth of water table Trend in Precipitation Aridity Rainfall of coldest quarter Diurnal temperature range Trend in PDSI Productivity (NDVI) Temperature seasonality Mean Temperature of wettest quarter Rainfall seasonality Ruminants density Precipitation variability Solar radiation Mean temperature of driest quarter 0.000 0.025 0.050 0.075 Importance Positive abrupt changes AB Fig. 2. Drivers of abrupt changes have contrasting importance depending on the trend of NDVI analyzed. Ranked importance for negative (A) and positive (B) abrupt changes. Importance is measured as gain and represents the improvement in accuracy of the RF model when a factor is included. The model fitted is a RF algorithm confronting abruptness (1 if abrupt, 0 otherwise) with 21 uncorrelated environmental factors summarizing climatic, soil, topographical, and human-related drivers. Total number of instances used to train/validate the model is 1955/839 for negative abrupt changes and 9135/3915 for positive abrupt changes. 4of10 https://doi.org/10.1073/pnas.2123393119 pnas.org Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. above 0 (warm conditions, SI Appendix, Fig. S9). Also worth mentioning is the higher importance of human-related factors (human footprint and ruminants density) in positive than negative abrupt changes (Fig. 2), with lower probability of abrupt positive shifts when these two factors were above a certain level (Fig. 3 Eand F). This result has important consequences as it corroborates that human pressure may impede rapid recovery of vegetation in drylands (4, 34) and it can be used as a management guide for the limits of human pressure under which restoration of vegetation will be most successful. Moreover, both high human footprint and ruminants density also indicate higher probability of abrupt losses of NDVI in areas more influenced by humans. Overall, the drivers of abrupt changes identified in our study indicate that multiple climatic conditions (e.g., temperature and rainfall) need to be taken into account simultaneously to explain when and why abrupt changes in ecosystem productivity occur. Our results also match expectations for aridity thresholds and highlight known dryland ecological (e.g., soil organic carbon) and socioeconomic (e.g., human footprint and grazing) feedbacks as drivers of abrupt shifts. One limitation on assessing the drivers of abrupt shifts is that the spatial scale used (around 1 km) does not account for fine-scale process, such as ecohydrological feedbacks or specific human disturbances, which can be more important at higher resolution, even if the identified relationships suggest promising results in this direction. Another limitation is that mechanistic knowledge cannot be inferred from the described relationships between drivers and the probability of abrupt shifts. Future studies should be conducted in situ to prospect more mechanistic knowledge of these drivers at the proper scales. Our study provides an atlas of potential sites that might be of interest to conduct such in situ observations. Ecosystem productivity has been relatively constant in half of all global drylands over the last 20 y. Yet, our results show that in the other half of dryland ecosystems plant productivity is changing in a predominantly abrupt rather than gradual way. We provide evidence for the hypothesized links between abrupt negative shifts and critical aridity thresholds, soil fertility, and grazing. Climatic variability and seasonality are the most consistent environmental factors explaining abrupt changes in NDVI. In contrast, the high incidence of abrupt positive changes found in low seasonality drylands suggests that restoration efforts could be particularly efficient to enhance plant productivity when climate variability is weak. Such context dependency is essential to anticipate responses to degradation processes and facilitate restoration efforts when managing vegetation in drylands. In this study, we created a global atlas of abrupt shifts in short-to-midterm NDVI dynamics, which helps to identify global hotspots of unstable dynamics pointing into abrupt degradation, but also an opportunity of productivity recovery in global drylands. This knowledge is integral to providing tools to improve the management of global drylands and to maintain or restore vegetation productivity and its associated ecosystem services in the largest biome on Earth. Materials and Methods Site Selection. We initially selected 59,000 sites of 30 arc-seconds (about 1km 2 ) interspersed 30 arc-minutes from each other across global drylands, which are defined as the zones of the Earth with aridity index (ratio of annual precipitation to annual potential evapotranspiration) lower than 0.65 (35). The aridity index was obtained from ref. 36 as the average aridity index for the period 1970–2000 with a spatial resolution of 30 arc-seconds. Because the purpose of our study is to estimate the drivers of abrupt shifts, we discarded areas where vegetation is primarily controlled by humans (such as urban or cultivated lands), where abrupt shifts can be a consequence of arbitrary human decisions (e.g., addition of fertilizers, rotatory crops, change in crop type, etc.), whose incidence cannot be controlled at a global scale with the variables we manage. Those areas, plus water bodies, were masked out using FAO’slanduse and land-cover maps (37) (spatial resolution of 30 arc-seconds). Finally, we also discarded sites classified as hyper-arid (aridity index <0.05), as those are sites usually lacking any vegetation and, thus, their NDVI estimates (see next section) are highly inaccurate (10). After this two-step filtering, 43,080 sites of the total 59,000 (71%) remained for the analysis. NDVI Trajectories. We gathered information on the NDVI as a surrogate of plant productivity. NDVI is an index that combines the spectral reflectance measurements acquired in the red and near-infrared bands and provides a global measure of the “greenness”of vegetation across the Earth’s landscapes (15, 38). It has been successfully linked with vegetation productivity in several studies Fig. 3. Effects of environmental drivers of abrupt shifts depending on the trend of productivity analyzed. Effects of seven factors (A-G: precipitation variability, rainfall seasonality, temperature seasonality, aridity, human footprint, ruminants density and soil organic carbon) for negative (red) and positive (blue) abrupt changes. The effect is expressed in SHAP values, which measure the impact of each factor on the model output (probability of abrupt change). SHAP valuesare derived for a given factor value in a process analogous to partial dependence plots. CV: coefficient of variation, SD: standard deviation. Rest of legend as Fig. 2. PNAS 2022 Vol. 119 No. 43 e2123393119 https://doi.org/10.1073/pnas.2123393119 5of10 Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. because it is directly related to the photosynthetically active radiation absorbed by plant canopies (39, 40). NDVI data were extracted from the MOD13Q1 product from the MODIS sensor (41), which provides 23 images per year (every 16 d to biweekly) with a pixel size of 250 m by 250 m in the period between 2000 and 2020 (https://modis.ornl.gov/data.html). MODIS data are geometrically and atmospherically corrected and include a quality index (QA) of data reliability based on the environmental conditions in which the data were recorded. The QA ranges from 0 to 1 (good quality data) to 4 (raw or absent data), with 2 being covered by snow and 3 being covered by clouds (42). To match the spatial scale of NDVI resolution (250 ×250 m) with that of our study sites (30 arc-seconds), for each site and each biweekly measurement we extracted all pixels up to 1 km 2 (16 pixels per time point) and averaged them across the 16 pixels for obtaining a spatial composite of the site per time point. Then the 23 biweekly time points of each year were averaged to remove the influence of seasonal variation, obtaining yearly NDVI averages. When averaging (either across space or time), we only used data that had a QA of 0, 1, or 2. Pixels with QA =2 (i.e., snow cover) were assumed to have NDVI =0 (as freezing temperature implies negligible photosynthesis). For the average across pixels, if more than 7 pixels were missing (out of the 16 possible), we assigned to this biweekly time point a missing value. Similarly, annual averages of NDVI were assigned a missing value when more than 10 (out of the 23 possible time points) were missing. In the end, we produced a time series of annual average NDVI that corresponded to a 20-y trajectory for each site. Out of the total 43,080 trajectories, we only used 41,830 (97%) sites for which we had at least 17 y of NDVI data (from the 20 possible years). MODIS data are only available since 2000. Although other products are available for longer periods such as the Advanced Very High Resolution Radiometer (AVHRR; which has data since 1981), we focused our study on the last 20 y of MODIS for several reasons. First, MODIS sensors provide higher spatial resolution than other products with larger timespans (AVHRR has resolution ∼9km and MODIS ∼250 m). This is important because, when rescaling to 1 km 2 grids (to match scales of the maps we used in the identification of drivers), we count with 16 pixels per location with values of NDVI in the MODIS data ensuring less noise and less missing values. This is critical for fitting the models to be used in the classification of trajectory types (see Trajectory Shapes section), which are highly sensitive to outliers (43). Second, other sensors such as LANDSAT may yield larger timespans at the expense of joining data from two sensors (LANDSAT 5 and LANDSAT 7), but, although LANDSAT data are corrected to be comparable, this correction might not be perfect in all cases inducing systematic mismatches in the time series that may affect the trajectory type classification procedure. Third, 20 y of data provided by MODIS is thought to be sufficient to capture vegetation changes dynamically (44). Fourth, by focusing on the last 20 years of productivity we are capturing the period that is mostly influenced by humans in terms of climate and land use change (45). It is also worth mentioning that previous studies on NDVI trends usually correct NDVI trends by precipitation in each year (46). This is a way of removing the fluctuations in precipitation and focusing only on trajectories of sustained degradation or recovery of productivity based on rainfall use efficiency for analyzing abrupt shifts (47). We, however, did not correct NDVI trends in the classification process but rather studied the effects of precipitation changes in a separate analysis of the drivers of trajectory types. We did so because first, we wanted to understand changes in productivity as a metric of ecosystem state. In dynamical studies, mixing the signals of the main driver (precipitation) and that of the state variable (productivity) usually confounds rather than clarifies the existence of a given shift (12). For instance, if the NDVI of a system decreases because rainfall decreases so that it changes into a low NDVI state, correcting by rainfall (e.g., by using rainfall use efficiency) could mask this change by yielding equal rainfall use efficiency when the system actually becomes less productive (12). Correcting by precipitation is important in trends analysis especially if we want to assess long-term trends in productivity (48) or when we focus on the functional efficiency of the system (21). However, here we are not interested in synchrony between precipitation and NDVI, nor on functional efficiency, but rather on changes in the productivity state of the system. Second, and most importantly, the approach of removing effects of precipitation from the NDVI signal relies on the assumption that the relationship between NDVI and precipitation is linear and known, yet, recent studies showed that sensitivity of productivity to climatic conditions (which in drylands basically measures the effects of precipitation on productivity) decreases with increasing aridity (3, 47, 49). This, together with the fact that increases in rainfall fluctuations are expected in the most arid lands (3), would render very noisy signals and, most importantly, this noise would be aridity dependent (higher in the most arid places than in less arid ones). Trajectory Shapes. To classify trajectory shapes, for each site we modeled the relationship between NDVI and year using different statistical models, each assuming a particular equation with an associated shape. In particular, we wanted to know whether the trajectories were linear (i.e., monotonic trends or no trends), curvilinear (with an acceleration or deceleration that makes the trajectory nonlinear), or abrupt (with a discontinuous trend tracking one particular NDVI value up to a year in which the value changed suddenly and then maintained until the end of the trajectory). These behaviors represent underlying mechanisms that inform about the dynamics of the ecosystem functioning with important consequences for the ecosystem stability (12) and monitoring procedures (50). To perform this classification, we used three statistical models, each associated with one of these behaviors: linear models (associated with linear shapes), quadratic models (associated with curvilinear shapes and able to identify smooth trend/slope changes), and step models (associated with abrupt shapes). All models fitted NDVI values vs. time (year). We fitted linear and quadratic models using generic lm functions in R (able to account for slope changes) and step models using the R package chngpt (51). We also fitted a linear model without a slope term to check the fitting to a no-trend typology, which is linear in shape but with only one parameter (intercept) and no slope. Prior to fitting, we normalized the NDVI data by subtracting the mean of the site and dividing by the SD of the site. This allowed us to compare coefficients across models. To select the best fitted model for each trajectory, we compared the Akaike Information Criteria (AIC) values of each fit. AIC is based on the log-likelihood of a given fit corrected by the complexity of the model. A lower value indicates a model that fits better with the data, but AIC differences lower than 2 units usually are taken as similarly good candidates (52). In those cases, we chose the simplest model based on the number of parameters of each model. We considered no-trend models the simplest, followed by linear, step, and quadratic. Although step and quadratic models have the same number of parameters (three), step regressions are considered simpler because the mathematical operations used only involve summation and not multiplication. To account for potential uncertainty in classifying trajectories due to their noisy nature and relative short length, we did not use a single regression of each model to the complete trajectory of a site. Instead, we bootstrapped each regression 100 times and compared between models for each bootstrapped iteration. In each bootstrap, we selected time points with a probability that is weighted by the inverse of the Mahalanobis distance of the point to the center of the distribution of year and NDVI values. This means that, in the bootstrap, we selected less frequently influential points (53). We registered the number of times that each model was selected as best-fit out of the 100 bootstraps, to identify the best fitted shape for each trajectory and from the shape that scored the highest percentage between all types. We used this percentage as a measure of uncertainty for best-fit shape of each trajectory (hereafter confidence of trajectory classification). Because the step regressions are sensitive to outliers (especially at the beginning and end of a timeseries) and harder to compare with other models, we used three additional criteria for testing the accuracy in classifying trajectories as “step-change.”First, we neglected step fits where the change point occurred in the three first or three last years of the trajectory. This ensures that abrupt changes are not detected because of anomalous points in the extremes. Also, this means that abrupt changes detected are remaining at least 4 y after the change, pointing to certain stability of the change detected. Second, for the trajectories that were classified as “step-change,”we registered the change point position (i.e., the year in which the trend is detected to change abruptly) and calculated the mean and SD of these change points across the 100 bootstrap iterations of each site. Trajectories in which step models were better fitted but presented high variability in the change point position are likely not correctly classified as “step-change.”To decide the value of change point SD that is critically influencing anomalous steps, we related the certainty in the bootstrap selection and the SD of change points in all 6of10 https://doi.org/10.1073/pnas.2123393119 pnas.org Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. sites. We found that both parameters were related: for trajectories in which SD of change point was lower than 0.68 there was a negative strong correlation between SD of change point and certainty in the bootstrap selection, whereas higher SDs in change point showed similar low certainty in this selection procedure (SI Appendix, Fig. S10). Thus, we only classified as “stepchange”those trajectories exhibiting lower SD of change points than 0.68. Third, we only considered abrupt changes that had at least 70% confidence in the bootstrap classification. This latter threshold was selected because it marks a segmented pattern between confidence in the model classification and the R-squared of the model residuals (SI Appendix, Fig. S11). This means that confidence of trajectory classification higher than 0.7 is associated with higher accuracy on NDVI fitting through time. Finally, we classified trajectories as “step-change”only when NDVI before and after the change point were significantly different according to a Wald test. We did so to ensure that the change produced after the change point was strong enough, matching the definition of a regime shift (shifts occurring between two states that are significantly different in functioning or structure). If a trajectory shape that was initially classified as “step-change”did not fulfill all these criteria, its type was reassigned to the shape of the next best-fit detected in the bootstrap (according to AIC). In other words, after this filtering procedure, the trajectory was assigned as linear-shape if it better fitted to linear models or no trend models or curvilinear-shape if it better fitted to quadratic models. It is important to note that the fulfilling of these criteria also allowed us to minimize errors produced due to multiple testing in our methods (note that after trajectory classification the results of these analyses are introduced into another one). These errors induced due to multiple testing are important to control and usually are handled by introducing more conservative approaches aimed at increasing robustness of metrics (54). Our approach allowed us to avoid influential points and errors in step classification minimizing these errors, but noise-induced errors and error propagation are still possible. Trajectory Trends. Each trajectory shape (linear, curvilinear, or abrupt) was also classified according to its overall trend as positive (increasing trend), negative (decreasing trend), or neutral (no trend). We did this based on the parameters of the bootstrapped regressions. Linear trajectories were classified as positive or negative according to the sign of the slope parameter, or neutral if the best-fit model was the no-trend model (a model with only an intercept). Abrupt shapes were classified as positive or negative according to the sign of the step parameter. Curvilinear shapes were classified according to the combination of aand bparameters of their corresponding quadratic regression (in a model of the form: y∼ax 2 +bx+c). It must be noted that the standardization of NDVI and year values in the previous section allows comparison in the same terms (i.e., aand bparameters have the same units across sites). Also note that the collection of parameters for each trajectory is performed by bootstrapping, thus allowing comparison of parameters with an associated confidence interval within a given trajectory. Attribution Analysis. We conducted an attribution analysis aimed at measuring the dependence of NDVI fluctuations in the time series of a given site on known climatic drivers that also fluctuate through time. Given the limited number of observations we had, we decided to include only 3 climatic factors in the attribution analysis, namely temperature, water availability (palmer drought severity index), and CO 2 . Usually, similar approaches have also included radiation as an important factor in this type of analysis (49, 55). We discarded it because radiation is arguably less limiting in dryland ecosystems (49) and also because its interannual variation usually covaries with that of temperature and water availability (49). Instead, we found more interesting testing for the interaction between CO 2 and water availability. This interaction is important in the contextofclimatechangebecauseCO 2 fertilization is argued to counteract the effect of droughts (48, 56) and there is an increasing debate about which of these factors may be more important (57, 58). We performed this attribution analysis as is shown in Jiao et al. (55). We downloaded information of climatic attributes from TERRACLIMATE (59), a platform for monthly simulated climatic data commonly used in these types of analyses (55). CO 2 concentration was taken from Mauna Loa (Scripps CO 2 program). For each site, we basically related the climatic attributes with annual NDVI from our time series using a linear model with interacting terms between CO 2 and water availability, and the temperature as a single term. Resulting regressions were then examined with the package relaimpo in R (60), in order to extract from them the explained variance per attribute (55). Then we grouped these explained variances according to trajectory types identified and conducted an analysis of variance (with Honest Significance Difference post hoc test) to distinguish homogeneous groups. Landscape Change Analysis. Because part of the abrupt behaviors identified can be a consequence of a profound structural change in the landscape (e.g., a change from forest into grassland), we conducted an additional analysis to associate changes through time of landscape classes with trajectory types. To do so, we used the MCD12Q2 MODIS product, which assigns the landscape class into one of the most common vegetation types and that have yearly information since 2000 (61). This product has temporal resolution of 1 year, meaning that the landscape is reassigned each year. We took the most common landscape class from the first 5 y of the study (2000–2005) and from the last five (2015–2019). We confronted all possible changes and reassigned them into more understandable categories meant to maximize representation of the vegetation classes that have been thought to be alternative states in ecosystems (25, 62, 63). Forests, Savannahs, and open ecosystems (grasslands and shrublands) were aggregated, aiming at representing key well-known alternative states of dryland ecosystems (25, 64, 65). We additionally summarized changes from grasslands into shrublands and vice versa, aimed at capturing shrub encroachment [another plausible alternative state shift (62)]. Finally, we also made distinctions of changes between barren soils (deserts) and other vegetated lands to account for desertification processes (63). We performed a χ 2 test between the different landscape changes possible and the trajectory types. We took the standardized residuals of the χ 2 test and plot them in a matrix map. Absolute values of residuals lower than the critical χ 2 value (3.46 for pairwise comparisons) were considered to be not significant and were set to 0. Environmental Factors. To assess the dependence of trajectory shapes and trends on the environmental conditions of a given site, we selected several variables summarizing site-specific conditions including: soil nutrient status, topographical attributes, human influence, and a set of comprehensive climate attributes related to climate change within the time series, historical climate legacies, continentality, intraand interannual climate variability, and the coupling between precipitation and temperature seasonality. For soil status, we downloaded from global maps (66) the following information on soil characteristics: sand content, nitrogen content, and organic carbon content. Soil nutrient content (carbon and nitrogen contents) is an important limiting factor of plant physiology, and it has been experimentally linked to photosynthetic rates of plants (67, 68). Soil sand content and soil carbon are tightly linked to soil water redistribution (notably infiltration), which reflects effective precipitation that is a major driver of NDVI (69, 70). We also introduced depth to water table (from ref. 71) and depth to bedrock (from ref. 66) to control for under-developed soils, which may induce special dynamics of soil water (increasing water runoff and showing very low water retention) thus affecting productivity. For topographic variables, we extracted elevation from ALOS-DSM digital elevation models (72). We did so because differences in photosynthetic responses due to air pressure depend on altitude (73, 74). For human influence, we extracted information of the human footprint (from ref. 75), ruminants density (from ref. 76), and population density (from ref. 77). These three factors have been widely used as an indirect standardized measurement of human impact in ecosystems (78, 79). In particular, human footprint is a metric measuring the density of humans, the distance to cities, roads or infrastructures, and the overall change induced in land of the zone (75). We used ruminants density as proxy of grazing intensity, which decreases photosynthetic production due to predation. Climate change or climate variation during the study period was assessed by downloading annual precipitation from the CHIRPS database (80) for the same period as NDVI and the annual Palmer Drought Severity Index (PDSI), which calculates precipitation anomalies with respect to evaporation (81) (extracted from ref. 59). Although both indices are proxies of water availability, precipitation is an indicator of water inputs, whereas PDSI is an indicator of droughts and water PNAS 2022 Vol. 119 No. 43 e2123393119 https://doi.org/10.1073/pnas.2123393119 7of10 Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. balance. We calculated for each site the trend of precipitation and PDSI using the Kendall tau correlation between precipitation/PDSI and year in those sites where we could retrieve the two indices for at least 14 y. Also, we calculated the number of years in the study period in which a drought has been declared (PDSI values lower than 3, see ref. 81). This allows accounting for sporadic water shortages apart from the particular trend in PDSI. We also gathered information on historical (averaged for the period 1980–2000) climate for each site. These climatic characteristics refer to the overall climate typology of a site. We extracted climatic characteristics from WorldClim (82), which include a wide range of variables related to average climate (i.e., precipitation averages), climate intra-annual variability (i.e., temperature seasonality and precipitation seasonality), and continentality (i.e., mean diurnal range of temperature), as well as coordination between temperature and precipitation within a year (i.e., mean temperature of wettest and driest quarter, precipitation of hottest and coldest quarter). To these variables, we added aridity, measured as 1-aridity index (described in Site selection). The aridity index is widely used as a good surrogate of net water availability (83, 84), being the main abiotic factor driving dryland structure and functioning (85). This metric is also commonly used to identify both the current (1, 2) and future (86, 87) extension of global drylands and in management/policy activities related with desertification and land use management (2). Also, we extracted variables related to interannual variability of precipitation, measured as the coefficient of variation of interannual precipitation for the period 1990–2000, which summarizes uncertainty in rainfall across years; an intrinsic characteristic defining important dryland features related to their ability to reach stable equilibrium (27, 88). Finally, we also included solar radiation (from ref. 89), which is important for photosynthesis as it provides the main source of energy for carboxylation processes. We used the averaged NDVI for the study period, as well as the total cover of vegetation (extracted from ref. 90) as factors that summarize the overall productivity of a site. These factors are important control variables that ensure that the detection of abrupt shifts is not related to the overall productivity in a system (e.g., abrupt shifts would be less common if productivity is already low just because they are harder to identify). We also included changes in landscape class as obtained in previous section as dummy variables in our models. We did so to account for plausible land use changes. From all these site-specific variables we excluded those that had a correlation higher than 0.7 to avoid problems of multicolineality. This means that we excluded soil nitrogen and sand content (which were highly correlated with soil organic carbon), total cover (correlated with NDVI and aridity), and rainfall (correlated with NDVI). All the variables extracted are provided with a resolution of 1 km 2 except for the PDSI and the precipitation (both ∼5km 2 ) where we compromised the spatial accuracy in exchange of temporal resolution (note that both indices are used to derive climatic trends). Machine Learning Model. We analyzed which environmental factors affected trajectory trends and trajectory shapes. To do this, we fitted a RF model to the trendtypeswefoundintheclassification using as covariates climatic, topographical, and soil variables (see Environmental Factors). RF is an algorithm that builds different classification trees using different sets of predictors and bootstrapped sites and formulates a consensus model that best fits the data (91). We used xgboost algorithms in R to perform the RF algorithm (92). For all models fitted, we used 70% of the data to train the RF keeping the rest for cross validation. For positive and negative trend types, we fitted two separate RF models using as the response variable a binomial variable indicating whether the shape of the trajectory was abrupt or not (i.e., grouping linear and curvilinear shapes and confronting them to abrupt shapes) to assess the relative importance of environmental factors on predicting the abruptness within each trend type. We used gbtree as a booster option. Because of the nature of our response variable, we used the “error of classification”as an evaluation metric and defined the objective function as binary logistic. Apart from these fixed parameters, we also tuned the next ones: eta (learning rate), maximum depth, gamma, minimum child weight, percentage of subsample per tree, number of folds, percentage of feature per tree, and number of trees. All these parameters are important when deciding how much the model should learn from the data to perform optimally without overfitting (a complete description of each one is done in SI Appendix, Table S2). The hyperparameter tuning was performed using a Bayesian optimization procedure using parBayesianOpt package in R and the resulting parameters are reported in SI Appendix,TableS2. In this procedure we also used 10-fold cross validation using xgb.cv function to avoid overfitting (error of test set is used instead of train set to decide the best performance, see 93). To interpret the RF outcome, we used SHAP values (Shapley Additive explanation 94). SHAP values are a way of interpreting machine learning models based on information theory (94, 95). In a nutshell, the SHAP value is derived from a regression tree model for a given feature and prediction. Its value is the effect of the predictor of interest in the model output for a given prediction and, thus it is provided in the same units than the response variable. SHAP values are actually homologous to evaluating the expression β 1 *x 1 in a regular multiple regression (y∼β 1 *x 1 +β 2 *x 2 +…+β n *x n ). This means that, essentially, a given predicted value of the model is the summation of all SHAP values obtained from the model given the values of predictors (96). By plotting the values of predictor vs. the associated SHAP values, we obtain a response curve analogous to the effects of that predictor over the response variable (after controlling for the effect of all others). SHAP values are widely used in machine learning (97), economics (94), security (95), and ecology (98). We extracted SHAP values with the package fastshap (99) and built-in functions of the package xgboost in R. The results obtained by this model were not influenced by spatial autocorrelation. Results of semivariograms on the residuals of these models showed a characteristic autocorrelation scale of 20–50 km for negative/positive models (SI Appendix,Fig.S12). We repeated the models taken this into account by weighting points proportionally to the number of points closer than this scale. By doing so, we underscore the influence of aggregated points, correcting the SHAP values, and importance of variables. Results were very similar to the ones reported in the main manuscript (SI Appendix, Figs. S13 and S14). Data, Materials, and Software Availability. Data for the analysis are freely available from NASA MODIS. Ancillary data used in the machine learning models are also freely available in their respective repositories. A table of sites used and a data frame containing the resulting classification of trajectories is available in figshare (100). The code used in the analyses is also available in figshare as an R-markdown describing step by step what is done up to the consecution of final figures of the manuscript. [Codes and processed data] data have been deposited in [figshare] (10.6084/ m9.figshare.16546278) (101). ACKNOWLEDGMENTS. The authors want to thank Fernando T. Maestre and Santiago Soliveres for useful comments and discussions in previous stages of the manuscript. M.B. acknowledges funding from Spanish Ministry of Science (Juan de la Cierva Fellowship FJCI-2018-036520-I). M.D.-B. acknowledges a Ram on y Cajal grant (RYC2018-025483-I) from Spanish Ministry of Science. T.W.C. acknowledges financial support from DOB ecology. 1. R. Pr av alie, Drylands extent and environmental issues. A global approach. Earth Sci. Rev. 161, 259–278 (2016). 2. M. Cherlet et al., Part Introduction, in World Atlas of Desertification, M. Cherlet et al., Eds. (Publication Office of the European Union, Luxembourg, ed. 3, 2018), pp. 11–19 3. M. Berdugo et al., Global Ecosystem thresholds driven by aridity. Science (80-) 367, 787–790 (2020). 4. B. T. Bestelmeyer et al., Desertification, land use, and the transformation of global drylands. Front. Ecol. Environ. 13,28–36 (2015). 5. C. Folke et al., Regime shifts, resilience, and biodiversity in ecosystem management. Annu. Rev. Ecol. Evol. Syst. 35, 557–581 (2004). 6. T. Filatova, J. G. Polhill, S. van Ewijk, Regime shifts in coupled socio-environmental systems: Review of modelling challenges and approaches. Environ. Model. Softw. 75, 333–347 (2016). 7. R. De Jong, J. Verbesselt, A. Zeileis, M. E. Schaepman, Shifts in global vegetation activity trends. Remote Sens. 5, 1117–1133 (2013). 8. R. de Jong, J. Verbesselt, M. E. Schaepman, S. de Bruin, Trend changes in global greening and browning: Contribution of short-term trends to longer-term change. Glob. Change Biol. 18, 642–655 (2012). 9. S. Horion et al., Mapping European ecosystem change types in response to land-use change, extreme climate events, and land degradation. Land Degrad. Dev. 30, 951–963 (2019). 10. P. N. Bernardino et al., Global-scale characterization of turning points in arid and semi-arid ecosystem functioning. Glob. Ecol. Biogeogr. 29, 1230–1245 (2020). 11. C. Abel et al., The human–environment nexus and vegetation–rainfall sensitivity in tropical drylands. Nat. Sustain. 4,25–32 (2020). 8of10 https://doi.org/10.1073/pnas.2123393119 pnas.org Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161. 12. Z. Ratajczak et al., Abrupt change in ecological systems: Inference and diagnosis. Trends Ecol. Evol. 33, 513–526 (2018). 13. S. K efi, M. Holmgren, M. Scheffer, When can positive interactions cause alternative stable states in ecosystems? Funct. Ecol. 30,88–97 (2016). 14. L. Turnbull et al., Understanding the role of ecohydrological feedbacks in ecosystem state change in drylands. Ecohydrology 5, 174–183 (2012). 15. T. N. Carlson, D. A. Ripley, On the relation between NDVI, fractional vegetation cover, and leaf area index. Remote Sens. Environ. 62, 241–252 (1997). 16. J. Wang, P. M. Rich, K. P. Price, W. D. Kettle, Relations between NDVI and tree productivity in the central Great Plains. Int. J. Remote Sens. 25, 3127–3138 (2004). 17. G. T. Yengoh, D. Dent, L. Olsson, A. E. Tengberg, C. J. Tucker, “Applications of NDVI for land degradation assessment BT”in Use of the Normalized Difference Vegetation Index (NDVI) to Assess Land Degradation at Multiple Scales: Current Status, Future Trends, and Practical Considerations, G. T. Yengoh, D. Dent, L. Olsson, A. E. Tengberg, C. J. Tucker III, Eds. (Springer International Publishing, Cham, 2016), pp. 17–25. 18. J. J. Gait an et al., Evaluating the performance of multiple remote sensing indices to predict the spatial variability of ecosystem structure and functioning in Patagonian steppes. Ecol. Indic. 34, 181–191 (2013). 19. A. L. Burrell, J. P. Evans, M. G. De Kauwe, Anthropogenic climate change has driven over 5 million km 2 of drylands towards desertification. Nat. Commun. 11, 3853 (2020). 20. X. Lu, L. Wang, M. F. McCabe, Elevated CO 2 as a driver of global dryland greening. Sci. Rep. 6, 20716 (2016). 21. S. Horion et al., Revealing turning points in ecosystem functioning over the Northern Eurasian agricultural frontier. Glob. Change Biol. 22, 2801–2817 (2016). 22. V. Paiva Alcoforado Rebello, A. Getirana, O. C. Rotunno Filho, V. Lakshmi, Spatiotemporal vegetation response to extreme droughts in eastern Brazil. Remote Sens Appl Soc Environ 18, 100294 (2020). 23. C. Chen et al., China and India lead in greening of the world through land-use management. Nat. Sustain. 2, 122–129 (2019). 24. T. Sternberg, H. Rueff, N. Middleton, Contraction of the Gobi desert, 2000–2012. Remote Sens. 7, 1346–1358 (2015). 25. J. G. Pausas, W. J. Bond, Alternative biome states in terrestrial ecosystems. Trends Plant Sci. 25, 250–263 (2020). 26. S.Schwinning,O.E.Sala,M.E.Loik,J.R.Ehleringer,Thresholds,memory,andseasonality: Understanding pulse dynamics in arid/semi-arid ecosystems. Oecologia 141,191–193 (2004). 27. S. Sullivan, R. Rohde, On non-equilibrium in arid and semi-arid grazing systems. J. Biogeogr. 29, 1595–1618 (2002). 28. H. F. Diaz, M. P. Hoerling, J. K. Eischeid, ENSO variability, teleconnections and climate change. Int. J. Climatol. 21, 1845–1862 (2001). 29. V. Guttal, C. Jayaprakash, Self-organization and productivity in semi-arid ecosystems: Implications of seasonality in rainfall. J. Theor. Biol. 248, 490–500 (2007). 30. T. M. Lenton, et al., Tipping elements in the Earth’s climate system. Proc. Natl. Acad. Sci. 105, 1786–1793 (2008). 31. A. Zhisheng et al., Global monsoon dynamics and climate change. Annu. Rev. Earth Planet. Sci. 43,29–77 (2015). 32. A. R. E. Sinclair, J. M. Fryxell, The Sahel of Africa: Ecology of a disaster. Can. J. Zool. 63, 987–994 (1985). 33. U. Helld en, Desertification: Time for an assessment? Ambio 20, 372–383 (1991). 34. S. Ravi, D. D. Breshears, T. E. Huxman, P. D’Odorico, Land degradation in drylands: Interactions among hydrologic–aeolian erosion and vegetation dynamics. Geomorphology 116, 236–245 (2010). 35. N. J. Middleton, D. S. Thomas, World atlas of desertification (United Nations Environment Programme/Edward Arnold, London, 1992). www.worldcat.org/isbn/0340555122. 36. R. J. Zomer, A. Trabucco, D. A. Bossio, L. V. Verchot, Climate change mitigation: A spatial analysis of global land suitability for clean development mechanism afforestation and reforestation. Agric. Ecosyst. Environ. 126,67–80 (2008). 37. J. Latham, R. Cumani, I. Rosati, M. Bloise, (2014) FAO Global Land Cover (GLC-SHARE) BetaRelease 1.0 Database, Land and Water Division (Rome). https://data.apps.fao.org/map/catalog/ srv/eng/catalog.search#/metadata/ba4526fd-cdbf-4028-a1bd-5a559c4bff38. 38. N. Pettorelli et al., Using the satellite-derived NDVI to assess ecological responses to environmental change. Trends Ecol. Evol. 20, 503–510 (2005). 39. C. J. Tucker, C. Vanpraet, E. Boerwinkel, A. Gaston, Satellite remote sensing of total dry matter production in the Senegalese Sahel. Remote Sens. Environ. 13, 461–474 (1983). 40. PRINCE SD, Satellite remote sensing of primary production: Comparison of results for Sahelian grasslands 1981-1988. Int. J. Remote Sens. 12, 1301–1311 (1991). 41. K. Didan, (2015) MOD13Q1 MODIS/Terra Vegetation Indices 16-Day L3 Global 250m SIN Grid V006 [Dataset]. NASA EOSDIS LP DAAC. doi.org/10.5067/MODIS/MOD13Q1.006. Accessed 20 May 2020. 42. C. O. Justice et al., An overview of MODIS Land data processing and product status. Remote Sens. Environ. 83,3–15 (2002). 43. P. Fearnhead, G. Rigaill, Changepoint detection in the presence of outliers. J. Am. Stat. Assoc. 114, 169–183 (2017). 44. S. Cusser, J. Helms 4th, C. A. Bahlai, N. M. Haddad, How long do population level field experiments need to be? Utilising data from the 40-year-old LTER network. Ecol. Lett. 24, 1103–1111 (2021). 45. A. H. Cooper, T. J. Brown, S. J. Price, J. R. Ford, C. N. Waters, Humans are the most significant global geomorphological driving force of the 21st century. Anthr Rev 5, 222–229 (2018). 46. D. M. Browning, J. J. Maynard, J. W. Karl, D. C. Peters, Breaks in MODIS time series portend vegetation change: Verification using long-term data in an arid grassland ecosystem. Ecol. Appl. 27, 1677–1693 (2017). 47. A. M. Ukkola et al., Annual precipitation explains variability in dryland vegetation greenness globally but not locally. Glob. Chang. Biol. 27, 4367–4380 (2021). 48. S. Piao et al., Characteristics, drivers and feedbacks of global greening. Nat. Rev. Earth Environ. 1,14–27 (2020). 49. A. W. R. Seddon, M. Macias-Fauria, P. R. Long, D. Benz, K. J. Willis, Sensitivity of global terrestrial ecosystems to climate variability. Nature 531, 229–232 (2016). 50. K. N. Suding, R. J. Hobbs, Threshold models in restoration and conservation: A developing framework. Trends Ecol. Evol. 24, 271–279 (2009). 51. Y. Fong, Y. Huang, P. B. Gilbert, S. R. Permar, chngpt: Threshold regression model estimation and inference. BMC Bioinformatics 18, 454 (2017). 52. H. Akaike, A new look at the statistical model identification. IEEE Trans. Automat. Contr. 19, 716–723 (1974). 53. X.-Q.Liu,F.Gao,Y.Wu,Z.Zhao,Detectingoutliersandinfluential points: An indirect classical Mahalanobis distance-based method. J. Stat. Comput. Simul. 88,1–21 (2018). 54. J. Cort es et al., Where are global vegetation greening and browning trends significant? Geophys. Res. Lett. 48, e2020GL091496 (2021). 55. W. Jiao et al., Observed increasing water constraint on vegetation growth over the last three decades. Nat. Commun. 12, 3777 (2021). 56. A. Koutavas, CO2 fertilization and enhanced drought resistance in Greek firs from Cephalonia Island, Greece. Glob. Change Biol. 19, 529–539 (2013). 57. C. D. Allen, D. D. Breshears, N. G. McDowell, On underestimation of global vulnerability to tree mortality and forest die-off from hotter drought in the Anthropocene. Ecosphere 6, art129 (2015). 58. J. Pe~ nuelas et al., Shifting from a fertilization-dominated to a warming-dominated period. Nat. Ecol. Evol. 1, 1438–1445 (2017). 59. J. T. Abatzoglou, S. Z. Dobrowski, S. A. Parks, K. C. Hegewisch, TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958-2015. Sci. Data 5, 170191 (2018). 60. U. Gr€ omping, Relative importance for linear regression in R: The package relaimpo. J. Stat. Softw. 17,1–27 (2006). 61. M. Friedl, J. Gray, D. Sulla-Menashe, MCD12Q2 MODIS/Terra+Aqua Land Cover Dynamics Yearly L3 Global 500m SIN Grid V006 (NASA EOSDIS Land Processes DAAC, 2019) [Data set]. 62. P. D’Odorico, G. S. Okin, B. T. Bestelmeyer, A synthetic review of feedbacks and drivers of shrub encroachment in arid grasslands. Ecohydrology 5, 520–530 (2012). 63. P. D’Odorico, A. Bhattachan, K. F. Davis, S. Ravi, C. W. Runyan, Global desertification: Drivers and feedbacks. Adv. Water Resour. 51, 326–344 (2013). 64. M. Hirota, M. Holmgren, E. H. Van Nes, M. Scheffer, Global resilience of tropical forest and savanna to critical transitions. Science (80-) 334, 232–235 (2011). 65. J. Verbesselt et al., Remotely sensed resilience of tropical forests. Nat. Clim. Chang. 6, 1028 (2016). 66. T. Hengl et al., SoilGrids250m: Global gridded soil information based on machine learning. PLoS One 12, e0169748 (2017). 67. C. Terrer, S. Vicca, B. A. Hungate, R. P. Phillips, I. C. Prentice, Mycorrhizal association as a primary control of the CO 2 fertilization effect. Science (80-) 353,72–74 (2016). 68. C. Terrer et al., Ecosystem responses to elevated CO 2 governed by plant-soil interactions and the cost of nitrogen acquisition. New Phytol. 217, 507–522 (2018). 69. A. F. Feldman et al., Moisture pulse-reserve in the soil-plant continuum observed across biomes. Nat. Plants 4, 1026–1033 (2018). 70. J. Mart ınez-Fern andez, A. Gonz alez-Zamora, L. Almendra-Mart ın, Soil moisture memory and soil properties: An analysis with the stored precipitation fraction. J. Hydrol. (Amst.) 593, 125622 (2021). 71. Y. Fan, H. Li, G. Miguez-Macho, Global patterns of groundwater table depth. Science (80-) 339, 940–943 (2013). 72. T. Tadono et al., Precise global DEM generation by ALOS PRISM. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci. 2, 71 (2014). 73. J. Gale, Plants and altitude–revisited. Ann. Bot. 94, 199 (2004). 74. N. Kumar, S. Kumar, S. K. Vats, P. S. Ahuja, Effect of altitude on the primary products of photosynthesis and the associated enzymes in barley and wheat. Photosynth. Res. 88,63–71 (2006). 75. O. Venter et al., Sixteen years of change in the global terrestrial human footprint and implications for biodiversity conservation. Nat. Commun. 7, 12558 (2016). 76. M. Gilbert et al., Global distribution data for cattle, buffaloes, horses, sheep, goats, pigs, chickens and ducks in 2010. Sci. Data 5, 180227 (2018). 77. A. J. Florczyk et al. (2019) GHSL data package 2019 (Luxembourg). https://data.europa.eu/doi/10. 2760/290498. Accessed 13 January 2021. 78. J. Sanderman, T. Hengl, G. J. Fiske, Soil carbon debt of 12,000 years of human land use. Proc. Natl. Acad. Sci. 114, 9575–9580 (2017). 79. M. Delgado-Baquerizo et al., Human impacts and aridity differentially alter soil N availability in drylands worldwide. Glob. Ecol. Biogeogr. 25,36–45 (2016). 80. C. Funk et al., The climate hazards infrared precipitation with stations–a new environmental record for monitoring extremes. Sci. Data 2, 150066 (2015). 81. W. C. Palmer, Meteorological Drought (US Department of Commerce, Weather Bureau, 1965). 82. R. J. Hijmans et al., Very high resolution interpolated climate surfaces for global land areas. Int. J. Climatol. 25, 1965–1978 (2005). 83. M. Delgado-Baquerizo et al., Decoupling of soil nutrient cycles as a function of aridity in global drylands. Nature 502, 672–676 (2013). 84. C. Wang et al., Aridity threshold in controlling ecosystem nitrogen cycling in arid and semi-arid grasslands. Nat. Commun. 5, 4799 (2014). 85. F. T. Maestre et al., Structure and functioning of dryland ecosystems in a changing world. Annu. Rev. Ecol. Evol. Syst. 47, 215–237 (2016). 86. J. Huang, H. Yu, X. Guan, G. Wang, R. Guo, Accelerated dryland expansion under climate change. Nat. Clim. Chang. 6, 166–171 (2015). 87. J. Yao et al., Accelerated dryland expansion regulates future variability in dryland gross primary production. Nat. Commun. 11, 1665 (2020). 88. H. von Wehrden, J. Hanspach, P. Kaczensky, J. Fischer, K. Wesche, Global assessment of the non-equilibrium concept in rangelands. Ecol. Appl. 22, 393–399 (2012). 89. D. N. Karger et al., Climatologies at high resolution for the earth’s land surface areas. Sci. Data 4, 170122 (2017). 90. C. Dimiceli et al.,MOD44B MODIS/Terra Vegetation Continuous Fields Yearly L3 Global 250m SIN Grid V006 (NASA EOSDIS Land Processes DAAC, 2015). 91. T. K. Ho, “Random decision forests”in Proceedings of 3rd International Conference on Document Analysis and Recognition (1995), August 15, 1995, Montreal, Canada. Publisher: IEEE Computer society. vol.1, pp. 278–282. PNAS 2022 Vol. 119 No. 43 e2123393119 https://doi.org/10.1073/pnas.2123393119 9of10 Downloaded from https://www.pnas.org by CSIC UNIDAD DE RECURSOS DE INFORMAC on November 3, 2022 from IP address 161.111.158.161.