Pronounced loss of Amazon rainforest resilience since the early 2000s

1Global Systems Institute, University of Exeter, Exeter, UK. 2Earth System Modelling, School of Engineering & Design, Technical University of Munich, Munich, Germany. 3Potsdam Institute for Climate Impact Research, Potsdam, Germany. ✉e-mail: c.a.boulton@exeter.ac.uk There is widespread concern about the resilience of the Amazon rainforest to land-use change and climate change. The Amazon is recognized as a potential tipping element in the Earth’s climate system1, is a crucible of biodiversity2 and usually acts as a large terrestrial carbon sink3,4. The net ecosystem productivity (carbon uptake flux) of the Amazon has, however, been declining over the last four decades and, during two major droughts in 2005 and 2010, the Amazon temporarily turned into a carbon source, due to increased tree mortality5–7. Several studies have suggested that deforestation8 and anthropogenic global warming9,10, especially in combination, could push the Amazon rainforest past critical thresholds11,12 where positive feedbacks propel abrupt and substantial further forest loss. Two types of positive feedback are particularly important. First, localized fire feedbacks amplify drought and associated forest loss by destroying trees13 and the fire regime itself may ‘tip’ from localized to ‘mega-fires’14. Second, deforestation and forest degradation, whether due to direct human intervention or droughts, reduce evapotranspiration and hence the moisture transported further westward, reducing rainfall and forest viability there15 and establishing a large-scale moisture recycling feedback. Net rainfall reduction may in turn reduce latent heating over the Amazon to the extent that it weakens the low-level circulation of the South American monsoon8. Model projections of future changes in the Amazon rainforest differ widely9,16–18. Early studies showed that the Amazon rainforest may exhibit strong dieback by the end of the twenty-first century9,19. Both pronounced drying in tropical South America and a weak CO2 fertilization effect18 contributed to this result, with dieback also more common under stronger greenhouse gas emission scenarios17. Other studies based on varying general circulation and vegetation model components show a wider range of results20,21. Nevertheless, the forest may be ‘committed’ to dieback despite appearing stable at the end of model runs16. This highlights the importance of measuring the changing dynamic stability of the forest alongside its mean state. Given the uncertainty in model projections, we directly analyse observational data for signs of resilience loss in the Amazon. The mean state of a system is not usually informative of changes in resilience; either one can change whilst the other remains constant16,17. Thus, higher-order statistical characteristics that respond more sensitively to destabilization than the mean need to be considered to quantify resilience. To measure the changing resilience of the Amazon rainforest, we use a stability indicator used to predict the approach of a dynamical system towards a bifurcation-induced critical transition. The predictability arises from the phenomenon of critical slowing down22,23 (CSD): as the currently occupied equilibrium state of a system becomes less stable, it responds more sluggishly to short-term perturbations (for example, weather variability for the Amazon). This loss of resilience, which is itself typically defined24 as the return rate from perturbations, reflects a weakening of negative feedbacks that maintain stability. The behaviour can be detected by an increase in lag-1 autocorrelation (AR(1)) in time series capturing the system dynamics25,26. It may also manifest as an increase in variance over time but variance can also be easily influenced by changing variability of the perturbations driving the system27. Increasing AR(1) has been used to detect CSD before bifurcation-induced state transitions in a number of systems, including but not limited to climate25,28 and ecology29. In particular, CSD has recently been detected in reconstructions of western Greenland ice sheet height changes30 as well as of the variability of the Atlantic Meridional Overturning Circulation31. A caveat, highlighted by analysis of model projections before Amazon dieback27, is that a system should be forced slower than its intrinsic response time scale for CSD to occur (Methods). Hence, the absence of CSD may not rule out the possibility of a forthcoming critical transition. Conversely, increasing AR(1) can sometimes occur for other physical reasons. A space-for-time substitution has previously revealed that tropical forest resilience as measured by mean AR(1) (on a grid cell basis) is lower for less annual rainfall sums32 but changes of Amazon resilience over time have not been investigated so far. We investigate controls on the resilience of the Amazon vegetation system and how its resilience has changed over the last three decades, in terms of a changing AR(1) coefficient as estimated from satellite-derived vegetation data. For comparison, we also investigate changes in variance over time, as a secondary indicator Pronounced loss of Amazon rainforest resilience since the early 2000s

T here is widespread concern about the resilience of the Amazon rainforest to land-use change and climate change. The Amazon is recognized as a potential tipping element in the Earth's climate system 1 , is a crucible of biodiversity 2 and usually acts as a large terrestrial carbon sink 3,4 . The net ecosystem productivity (carbon uptake flux) of the Amazon has, however, been declining over the last four decades and, during two major droughts in 2005 and 2010, the Amazon temporarily turned into a carbon source, due to increased tree mortality [5][6][7] . Several studies have suggested that deforestation 8 and anthropogenic global warming 9,10 , especially in combination, could push the Amazon rainforest past critical thresholds 11,12 where positive feedbacks propel abrupt and substantial further forest loss. Two types of positive feedback are particularly important. First, localized fire feedbacks amplify drought and associated forest loss by destroying trees 13 and the fire regime itself may 'tip' from localized to 'mega-fires' 14 . Second, deforestation and forest degradation, whether due to direct human intervention or droughts, reduce evapotranspiration and hence the moisture transported further westward, reducing rainfall and forest viability there 15 and establishing a large-scale moisture recycling feedback. Net rainfall reduction may in turn reduce latent heating over the Amazon to the extent that it weakens the low-level circulation of the South American monsoon 8 . Model projections of future changes in the Amazon rainforest differ widely 9,[16][17][18] . Early studies showed that the Amazon rainforest may exhibit strong dieback by the end of the twenty-first century 9,19 . Both pronounced drying in tropical South America and a weak CO 2 fertilization effect 18 contributed to this result, with dieback also more common under stronger greenhouse gas emission scenarios 17 . Other studies based on varying general circulation and vegetation model components show a wider range of results 20,21 . Nevertheless, the forest may be 'committed' to dieback despite appearing stable at the end of model runs 16 . This highlights the importance of measuring the changing dynamic stability of the forest alongside its mean state. Given the uncertainty in model projections, we directly analyse observational data for signs of resilience loss in the Amazon.
The mean state of a system is not usually informative of changes in resilience; either one can change whilst the other remains constant 16,17 . Thus, higher-order statistical characteristics that respond more sensitively to destabilization than the mean need to be considered to quantify resilience. To measure the changing resilience of the Amazon rainforest, we use a stability indicator used to predict the approach of a dynamical system towards a bifurcation-induced critical transition. The predictability arises from the phenomenon of critical slowing down 22,23 (CSD): as the currently occupied equilibrium state of a system becomes less stable, it responds more sluggishly to short-term perturbations (for example, weather variability for the Amazon). This loss of resilience, which is itself typically defined 24 as the return rate from perturbations, reflects a weakening of negative feedbacks that maintain stability. The behaviour can be detected by an increase in lag-1 autocorrelation (AR(1)) in time series capturing the system dynamics 25,26 . It may also manifest as an increase in variance over time but variance can also be easily influenced by changing variability of the perturbations driving the system 27 . Increasing AR(1) has been used to detect CSD before bifurcation-induced state transitions in a number of systems, including but not limited to climate 25,28 and ecology 29 . In particular, CSD has recently been detected in reconstructions of western Greenland ice sheet height changes 30 as well as of the variability of the Atlantic Meridional Overturning Circulation 31 . A caveat, highlighted by analysis of model projections before Amazon dieback 27 , is that a system should be forced slower than its intrinsic response time scale for CSD to occur (Methods). Hence, the absence of CSD may not rule out the possibility of a forthcoming critical transition. Conversely, increasing AR(1) can sometimes occur for other physical reasons. A space-for-time substitution has previously revealed that tropical forest resilience as measured by mean AR(1) (on a grid cell basis) is lower for less annual rainfall sums 32 but changes of Amazon resilience over time have not been investigated so far.
We investigate controls on the resilience of the Amazon vegetation system and how its resilience has changed over the last three decades, in terms of a changing AR(1) coefficient as estimated from satellite-derived vegetation data. For comparison, we also investigate changes in variance over time, as a secondary indicator The resilience of the Amazon rainforest to climate and land-use change is crucial for biodiversity, regional climate and the global carbon cycle. Deforestation and climate change, via increasing dry-season length and drought frequency, may already have pushed the Amazon close to a critical threshold of rainforest dieback. Here, we quantify changes of Amazon resilience by applying established indicators (for example, measuring lag-1 autocorrelation) to remotely sensed vegetation data with a focus on vegetation optical depth (1991-2016). We find that more than three-quarters of the Amazon rainforest has been losing resilience since the early 2000s, consistent with the approach to a critical transition. Resilience is being lost faster in regions with less rainfall and in parts of the rainforest that are closer to human activity. We provide direct empirical evidence that the Amazon rainforest is losing resilience, risking dieback with profound implications for biodiversity, carbon storage and climate change at a global scale.
for CSD. The main dataset we use is from the Vegetation Optical Depth Climate Archive (VODCA) 33 but we also analyse the NOAA Advanced Very-High-Resolution Radiometer's (AVHRR) Normalized Difference Vegetation Index (NDVI) 34 for comparison. Vegetation optical depth (VOD) has been previously used to estimate changes in vegetation biomass 35 , whereas NDVI is more commonly used to measure the greenness of vegetation (that is, photosynthetic activity 36 ), which can saturate at dense grass cover. VOD, a microwave-derived product, does not saturate and remains sensitive to changes also at high biomass density 35 . We use the longer, Ku-band product from VODCA, which has a resolution of 0.25° × 0.25°, from which we create a monthly dataset by taking the mean values and for direct comparison we rescale the NDVI data to the same resolution. For a check of robustness, we also use the C-band product from VODCA, which spans a shorter time period. We focus on two stressors of the Amazon that may cause resilience changes-precipitation decline and human influence.

results
We use the Amazon basin as our study region and focus on those grid cells that have a ≥80% evergreen broadleaf (BL) fraction according to the MODIS Land Cover Type product in 2001 37  Figure 1 shows that mean changes in BL fraction in this region correspond well to changes in VOD. Averaged across the Amazon study region we find overall decreasing VOD over 2001-2016, corresponding to the observed decrease in the number of grid cells that have BL ≥ 80% each year (Fig. 1a). Between 2001 and 2016, the BL fraction has changed most prominently in the south-eastern parts of the Amazon basin, along parts of the Amazon River and in some northern parts of the basin (Fig. 1b). Changes in VOD have a similar spatial pattern to changes in BL fraction, with decreases concentrated around the south-eastern edges of the forest (Fig. 1c). NDVI, in contrast, does not agree temporally or spatially with the changes in BL fraction ( Supplementary Fig. 1), with NDVI increases observed in the south-eastern parts of the Amazon where deforestation rates are known to be highest. On the individual grid cell level, changes in BL are strongly correlated with changes in VOD (Pearson's r = 0.556; Supplementary Fig. 2a) compared to changes in NDVI (r = −0.133; Supplementary Fig. 2b). This echoes previous in-situ comparisons between VOD and NDVI 38 . Hence, we focus our analysis on VOD in the following, with results for NDVI in the Supplementary Figures.
We begin our resilience analysis by focusing on the temporal changes of AR(1), computed in sliding windows from the nonlinearly detrended and deseasonalized VOD time series ( Fig. 2 and Methods). We remove forested (BL ≥ 80%) grid cells that have any human land use in them (Methods), resulting in 6,369 grid cells being analysed. The spatial distribution of the AR(1) tendency, measured by the Kendall's rank correlation coefficient τ (Methods) at each grid cell, shows that the AR(1) increases in most of the grid cells comprising the Amazon rainforest (Fig. 2a,b). Likewise, the time series calculated from the mean AR(1) value across our study area each month shows a substantial increase over time, particularly from the early 2000s (Fig. 2c). We observe some stable or decreasing AR(1) values around the tributaries of the Amazon River, suggesting increasing resilience. VOD values can be influenced by open water but this should be minimized by looking only at grid cells where BL ≥ 80%. Previously, floodplain forests near the river, which cover 14% of the basin, have been shown to be much less resilient than the non-flooded forests throughout the Amazon 39 . However, when we compare grid cells in the Amazon basin floodplains and those outside (Methods), we find no significant differences between the resilience signals (Mann-Whitney U-test P = 0.579), except for a smaller Kendall τ value for floodplains from 2003 onwards (0.863 compared to 0.915; Supplementary Fig. 3).
Overall, most (76.2%) of the grid cells show increasing AR(1) values from the early 2000s onwards and hence, loss of resilience ( Fig. 2b), as well as 77.8% of grid cells over the full time period. Using alternative methods of detrending the VOD time series (Methods) yields similar results (Supplementary Figs. 4 and 5), as does varying the window length used to estimate AR(1) (Supplementary Fig. 6). The results are also robust to the choice of the BL fraction threshold used to determine forest, finding increasing AR(1) also when either BL ≥ 90% or BL ≥ 40% is used ( Supplementary Fig. 7). Furthermore, restricting the analysis to those grid points that have BL ≥ 80% through the period 2001-2016, rather than only checking grid cells for BL ≥ 80% in 2001, shows similar signals of resilience loss ( Supplementary Fig. 8).
A predominance of increasing AR(1) trends is also found for the NDVI time series since 2003 ( Supplementary Fig. 9).
To try and determine the causes for the detected resilience loss across the Amazon basin, we explore the relationships between the AR(1) trends and mean annual precipitation (MAP, estimated from the CHIRPS dataset 40 ), as well as between the AR(1) trends and the distance from human land use (Methods). It has previously been suggested that drier forest is less resilient 32 as well as forest nearer human land use 41 . We also include distance from roads alongside land use but this restricts the domain of analysis to Brazil to avoid biases by heterogeneities in road data across different countries (Supplementary Fig. 10; Methods). Figure 3 compares the spatial patterns of the AR(1) trends, MAP and human land use. Relationships with both explanatory variables are discernible, with (less common) decreases in AR(1) (Fig. 3a) being found in regions of high MAP in the north of the region (Fig. 3b) and further away from human land use (Fig. 3c). We find no confounding relationship between MAP and human land use; they are only very weakly correlated with each other (Spearman's ρ = 0.057, P < 0.001). Hence, we can consider them as separate relationships. The computed minimum distances to human land use and roads should be interpreted as upper bounds because for the full region we do not include roads, and for the Brazil distances our dataset will not include non-federal or non-state roads, which also have human activity associated with them. Furthermore, the classification of grid cells that contain human land use at the spatial resolution used in this analysis is likely to only detect large farms and settlements.
To further explore the relationship between MAP and AR(1) trends, we create mean AR(1) time series on a moving MAP band of 500 mm (Methods). These bands show broadly the same behaviour as the region overall (Fig. 4a), with all bands showing a significant decrease in resilience post-2003 (P < 0.001 for all MAP bands). The increase in AR (1) post-2003 appears least pronounced for the highest rainfall band (3,500-4,000 mm). The strength of resilience loss increases as the MAP band decreases below 3,500-4,000 mm (Fig. 4b). For NDVI, the same relationship is also observed ( Supplementary Fig. 11a,b). However, due to a large decrease in NDVI AR(1) pre-2003 across the region, analysing the full NDVI AR(1) time series yield decreasing AR(1) Kendall τ coefficients for the higher MAP bands.
To further explore the relationship between resilience and the distance to human land use, we calculate mean AR(1) time series on 50 km distance bands. Our results show that increases in AR(1) post-2003 are stronger for grid cells closer to human land use (Fig. 5a). Above 200-250 km away from human land use the signal of loss of resilience becomes less pronounced (Fig. 5b). Using the subset of Brazilian grid cells (n = 3,797) to include roads in our measurement of distance to human activity (Supplementary Fig. 10; Methods) generally reduces distances but we find a similar relationship, with decreases in Kendall τ coefficients observed up to 75 km (Fig. 5c,d). We note in both cases that, at longer distances, the number of grid cells used to calculate the AR(1) time series is lower (red lines and right-hand y axis in Fig. 5b,d). The remaining areas also tend to become more separated (the two vertical lines in Fig. 5b,d mark where <100 and <50 grid cells were left for the calculations, respectively). This helps to explain the more variable AR(1) time series at greater distances (for example, the yellow line in Fig. 5a) and the more fluctuating results for Kendall τ coefficients at greater distances (Fig. 5b,d). NDVI time series also show a loss of resilience from the early 2000s, in grid cells that are closer than 200 km from human land use ( Supplementary Fig. 11c,d). We reiterate that the stated minimum distances to human land use should be viewed as upper limits with, for example, selective logging and other intrusions expected to be closer to the forest than the grid cells with major human land use and major roads. Comprehensive robustness tests, using alternative datasets and indicators, are presented in the Methods. In particular, we find overall consistent results when using the variance instead of the AR(1) as measure of resilience.

Discussion
We reiterate that changes in the mean state of a system do not directly relate to changes in resilience. Model studies show that large parts of the Amazon rainforest can be committed to dieback 16 before showing a strong change in mean state. Indeed, from our CSD indicators we infer a marked loss of Amazon rainforest resilience since the early 2000s, in vast areas where the BL fraction has not strongly decreased (compare Figs. 1b and 2a or Supplementary Fig. 8b).
Given that lower baseline MAP (Fig. 4) and greater proximity to human interference (Fig. 5) are both statistically associated with greater loss of resilience, we hypothesize that low MAP and increasing human interference could both be contributing to the large-scale loss of resilience (Fig. 2). What remains to be explained is why these two factors might play such an important role and why the large-scale resilience loss started in the early 2000s.
Previous work 32 has shown that regions with lower MAP have lower absolute forest resilience, conceivably because vegetation is more water stressed and struggles to regulate its internal water content. VOD, being dependent on this water content, would consequently adjust more slowly to perturbations. Our finding that resilience has been lost faster in lower MAP regions, additionally suggests that vegetation in regions with more pronounced aridity stress is at greater risk of losing resilience. Large parts of the study region show decreasing MAP. However, the spatial pattern of MAP change (as measured by the difference between the means for January 1998 to December 2002 and January 2012 to December 2016; Supplementary Fig. 12) is different to that of AR(1) increases (Fig. 2) and we find no spatial correlation between these MAP changes and VOD AR(1) Kendall τ (Spearman's ρ between the spatial field of MAP change and the spatial field of VOD AR(1) Kendall  (1) in or near the drought years. These peak values may appear earlier due to the aR(1) time series being calculated on a moving window, compared to the SST anomalies being a monthly mean. τ equals 0.04). Increases in dry-season length (DSL; Supplementary  Fig. 12) reported in several recent studies [42][43][44][45] , might conceivably be a better explanatory variable but again we find no spatial correlation between DSL change and VOD AR(1) Kendall τ (ρ = −0.08). The lack of spatial correlations for both MAP and DSL could be due to the relatively short period to measure rainfall trends and for DSL due to the discrete nature of DSL values, which are given in units of full months, compared to the continuous τ values.
Despite a lack of spatial correlations, existing understanding and larger-scale aggregate measures suggest that climate variability may be among multiple factors contributing to the observed Amazon resilience loss since the early 2000s (Fig. 6). Notably, the Amazon shows signs of resilience loss during a period with three 'one-in-a-century' droughts 10,[46][47][48] . Sea surface temperature (SST) anomalies in the northern tropical Atlantic Ocean (from the HadISST 49 datatset; Methods) from around 2000 onwards have been generally positive compared to climatology, consistent with a shift of the Atlantic Multidecadal Oscillation (AMO) to its positive phase ( Supplementary Fig. 13), although reductions in anthropogenic northern-hemisphere aerosol cooling may also play a role 10 . These positive northern tropical Atlantic SST anomalies led-via the associated northward shift of the Intertropical Convergence Zone-to drier conditions in the Amazon and, in particular, to two severe droughts in 2005 and 2010 46,48 (Fig. 6a). These two drought events are associated with corresponding peaks in the spatial-mean AR(1) time series, superimposed on the overall positive trend (see arrows in Fig. 6b). These peaks are also found in the separate AR(1) time series in Fig. 4a, appearing about 2.5 yr ahead due to the time series being plotted at the end of the 5-yr sliding windows used to calculate the AR(1) there. Moreover, a third, El Niño-driven drought in 2015/16 is accompanied by an increased overall rate of resilience loss at the very end of the time range for which the VOD data are available. The decrease in AR(1) before the early 2000s may also be linked to internal climate variability; the AMO was in its negative phase ( Supplementary Fig. 13), consistent with negative SST anomalies in the northern tropical Atlantic (Fig. 6) and wetter conditions in the Amazon. The fact that the AR(1) increase since the early 2000s is statistically strongly significant suggests that it is not just due to natural climate variability.
Increasing human land use also appears to be contributing to the observed Amazon resilience loss, with human land-use areas increasing in both reach and intensity ( Fig. 6c and Supplementary  Fig. 14). Notably, the expansion of human land use accelerates after 2010, in an interval that also shows accelerated resilience loss ( Fig. 6b) but less striking northern tropical Atlantic SST anomalies (Fig. 6a). Greater proximity to human land use can increase disturbance factors such as direct removal of trees, construction of roads and fires, conceivably reducing absolute resilience (Fig. 5) and making the forest more prone to resilience loss.
Other factors, including rising atmospheric temperatures in response to anthropogenic greenhouse gas emissions, may additionally have negative effects on Amazon resilience (and are contributing to the warming of northern tropical Atlantic SSTs; Fig. 6a). Furthermore, the rapid change in climate is triggering ecological changes but ecosystems are having difficulties in keeping pace. In particular, the replacement of drought-sensitive tree species by drought-resistant ones is happening slower than changes in (hydro)meteorological conditions 50 , potentially reducing forest resilience further.
In summary, we have revealed empirical evidence that the Amazon rainforest has been losing resilience since the early 2000s, risking dieback with profound implications for biodiversity, carbon storage and climate change at a global scale. We further provided empirical evidence suggesting that overall drier conditions, culminating in three severe drought events, combined with pronounced increases in human land-use activity in the Amazon, probably played a crucial role in the observed resilience loss. The amplified loss of Amazon resilience in areas closer to human land use suggests that reducing deforestation will not just protect the parts of the forest that are directly threatened but also benefit Amazon rainforest resilience over much larger spatial scales.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41558-022-01287-8.

Datasets.
We use the Amazon basin (http://worldmap.harvard.edu/data/ geonode:amapoly_ivb, accessed 28 January 2021) as our region of study. To determine the grid cells that are contained within Brazil for a subset of analysis, we use the 'maps' package in R (v.3.3.0; https://CRAN.R-project.org/package=maps). This is also used in the plotting of country outlines. The main dataset used to determine forest health is from VODCA 33 , of which we use the Ku-band product. These data are available at 0.25° × 0.25° at a monthly resolution from January 1988 to December 2016. We also use NOAA AVHRR NDVI 34 . For precipitation data, we use the CHIRPS dataset 40 downloaded from Google Earth Engine at a monthly resolution. Finally, to determine land cover types, we used the IGBP MODIS land cover dataset MCD12C1 (ref. 37 ). All these datasets are at a higher spatial resolution than the VODCA dataset and thus we downscale them to match the lower resolution. Our SST data comes from HadISST 49 , where we define a North Atlantic region (15-70° W, 5-25° N), for which we take the spatial mean. The mean monthly cycle is then removed to produce anomalies.
For the vegetation datasets that we measure the resilience indicators on (below), we use STL decomposition (seasonal and trend decomposition using Loess) 51 using the stl() function in R. This splits time series in each grid cell into an overall trend, a repeating annual cycle (by using the 'periodic' option for the seasonal window) and a residual component. We use the residual component in our resilience analysis. The first 3 yr of data had large jumps in VOD which were seen when testing other regions of the world as well as in the Amazon region. Hence, we restrict our analysis to the period January 1991 to December 2016.
To test the robustness of the detrending, we also vary the size of the trend window in the stl() function. The results from these alternatively detrended time series are shown in Supplementary Fig. 4. The results are also robust to varying the window used to calculate the seasonal component rather than using 'periodic'; at the strictest plausible value of 13, we still see the same increases in AR(1) (Supplementary Fig. 5).
For the AMO index shown in Supplementary Fig. 13, data come from the Kaplan SST dataset and can be downloaded from https://psl.noaa.gov/data/ timeseries/AMO/.
Grid cell selection. We use the IGBP MODIS land cover dataset at the resolution described above to determine which grid cells to use in our analysis. The dataset is available at an annual resolution from 2001 to 2018 (but we only use the time series up to 2016 to match the time span of our VOD and NDVI datasets). To focus on changes in forest resilience, we use grid cells where the evergreen BL fraction is ≥80% in 2001. Grid cells are treated as human land-use area if the built-up, croplands or vegetation mosaics fraction is >0%. We remove grid cells that have human land use in them from our forest analysis, regardless of if there is ≥80% BL fraction in the grid cell.
We measure the minimum distance between forested Amazon basin grid cells and human land-use grid cells in 2016 (believing this to be the most cautious and least biased way to measure distance) using the latitude and longitude of each grid point and computing the great-circle distance. We use human land-use grid cells over a larger area than the basin, so that we can determine the closest distance to human land use, regardless of whether this human land use lies within the basin. We also measure the minimum distance from human land use or roads in Brazil, where we have reliable data on state and federal roads (https://datacatalog. worldbank.org/dataset/brazil-road-network-federal-and-state-highways). As in the main text, we reiterate that these minimum distances can be viewed as the maximum distance from human land use as our data will not include roads for the full Amazon basin, or non-federal or non-state roads in Brazil that will have human activity associated with them.
To ensure that the pattern of changes in resilience is not a consequence of more settlements being in the southeast of the region, combined with the gradient of rainfall from northwest to southeast typical of the rainforest, we measure the correlation between MAP and the distances from the urban grid cells, which is very weak (Spearman's ρ = 0.109, P < 0.001) and as such we are confident that there are separate processes that causes these relationships.

Resilience indicator AR(1).
We measure our resilience indicator on the residual component of the decomposed vegetation time series. We focus on AR(1), which provides a robust indicator for CSD before bifurcation-induced transitions and has been widely used for this purpose 23,25,32 . We measure it on a sliding window length equal to 5 yr (60 months). The sliding window creates a time series of the AR(1) coefficient in each location. Our results are robust to the sliding window length used, as shown in Supplementary Fig. 6.
From linearization and the analogy to the Ornstein-Uhlenbeck process, it holds approximately that for discrete time steps of width Δt (1 month in the case at hand): where κ is the linear recovery rate. A decreasing recovery rate κ implies that the system's capability to recover from perturbations is progressively lost, corresponding to diminishing stability or resilience of the attained equilibrium state. From the above equation it is clear that the AR(1) increases with decreasing κ. The point at which stability is lost and the system will undergo a critical transition to shift to a new equilibrium state, corresponds to κ = 0 and AR(1) = 1, respectively.
Measuring AR(1) across the whole time series provides information about the characteristic time scales of the two vegetation datasets we use 26 . Inverting κ gives the characteristic time scale of the system; for the VOD, we find 1/κ = 1.240 months, whereas for the NDVI, we find 1/κ = 0.838 months when using the mean AR(1) value across the region. This suggests that, in accordance with our interpretation of the two satellite-derived variables, the NDVI is more sensitive to shorter-term vegetation changes such as leaf greenness, while the VOD's Ku-band is sensitive to longer-term changes such as variability in the thickness of forest stems.
Creation and tendency of AR(1) and variance time series. For analyses where either MAP bands or distance bands are used to create an AR(1) or variance series, we calculate the mean AR(1) or variance value in each month for forested (BL ≥ 80%), non-human land-use Amazon basin grid cells, from which the tendency of this mean series can be calculated. Alternatively, the Kendall τ for each band can be calculated by taking the mean Kendall τ for each individual grid cell that is within the band. Results from the first option are shown in Figs. 4 and 5 and results from the second method in Supplementary Fig. 15 for AR(1).
The tendencies of the CSD indicators are determined in terms of Kendall τ. This is a rank correlation coefficient with one variable taken to be time. Kendall τ values of 1 imply that the time series is always increasing, −1 implies that the time series is always decreasing and 0 indicates that there is no overall trend. Following previous work 25,52,53 , we test the statistical significance of positive tendencies using a test based on phase surrogates that preserve both the variance and the serial correlations of the time series from which the surrogates are constructed. Specifically, we compute the Fourier transform of each time series for which we want to test the significance of Kendall τ, then randomly permute the phases and finally apply the inverse Fourier transform. Since this preserves the power spectral density, it also preserves the autocorrelation function due to the Wiener-Khinchin theorem. For each time series this procedure is repeated 100,000 times to obtain the surrogates. Kendall τ is computed for each surrogate to obtain the null model distribution (corresponding to the assumption of the same variance and autocorrelation but no underlying trend), from which we calculate a P value by calculating the proportion of surrogates that have a higher Kendall τ value.

Robustness tests.
To account for the possibility of human deforestation interfering with the signals we observe (which may not necessarily be detected by the MODIS Land Cover dataset) we also use the Hansen forest loss dataset 54 to determine grid cells to remove in an alternate analysis. The original Hansen dataset is at a 0.00025° resolution, 1,000 times higher than the VOD dataset and as such for each VOD grid cell we measure the percentage of Hansen grid cells that show some forest loss over the time period. Note that this dataset does not specify if the observed loss is natural or caused by human deforestation. Excluding any VOD grid cells that contain more than a conservative 5% of lost forest grid cells according to this dataset and running the analysis in the main paper shows similar results. Supplementary Figs. 16-18 are recreations of Figs. 2, 4 and 5, respectively.
The loss of forest resilience observed as increasing AR(1) in both vegetation indices is supported by another indicator of CSD, namely increasing variance 28of both VOD (Supplementary Fig. 19) and NDVI ( Supplementary Fig. 20). Variance is more strongly affected by changes in the frequency and amplitude of the forcing of a system and as such results could be biased towards individual events. However, we assessed the precipitation time series for changes in variance and found no relationship with the variance signals of VOD and NDVI ( Supplementary Fig. 21). Nevertheless, AR(1) is considered the more robust indicator 55 . As another test of robustness, we partition the grid cells into those that are in floodplains and those that are not ( Supplementary Fig. 3). Floodplain data are part of the NASA Large Scale Biosphere-Atmosphere Experiment (LBA-ECO) 56 . We also calculate the resilience signals for the C-band product of VOD for comparison ( Supplementary Fig. 22). Despite the smaller temporal scale of this product, we still see increases in both AR(1) and variance. To account for a change in the number of satellites used to calculate VOD, for the Ku-band we also recreate the dataset by sampling a single random day per month rather than taking a monthly average, to mimic a constant satellite pass for the whole time period (Supplementary Fig. 23). Although this expectedly affects the absolute values of AR(1) and variance, their relative changes over time remain unaffected. To further test the robustness of our results, we looked for similar signals of resilience change in terms of trends in AR(1) in addition to variance in the precipitation time series used, as a change in the forcing could have an impact on the forest that could mistakenly be interpreted as a vegetation resilience loss. However, as for the variance, there is no clear increase in the AR(1) of precipitation, nor do the spatial patterns of both indicators reveal any relationship between changing precipitation AR(1) and variance and the observed vegetation resilience loss. Hence, we are confident that changes in precipitation forcing are not driving the vegetation AR(1) signals.