Higher tree diversity reduces the likelihood of Amazon tipping points

The Amazon forest is influenced by strong feedback loops between its biotic and abiotic components. Local forest loss increases CO₂ emissions, which, in turn, drives climate change, raising temperatures and reducing rainfall, causing further forest loss. Additionally, forest loss disrupts important forest-rainfall cycles, threatening the overall forest stability. These feedbacks make the system vulnerable to tipping points, where parts of the forest could transition to a degraded state. Slower recovery to short-term disturbances, hereafter named reduced stability, is considered an early warning indicator of such tipping points. However, the role of tree species diversity in regulating this vulnerability remains poorly understood, especially across spatial scales. To examine how tree species diversity impacts tipping point likelihoods across multiple spatial scales, we used modelled tree species diversity data at the alpha (local), beta (asynchrony across local communities), and gamma (regional) scales. We quantified tipping likelihood on the same scales using temporal autocorrelation trends in monthly satellite-derived vegetation productivity time series over 2001–2019. Our findings reveal higher tipping likelihoods at the alpha level (25 km²) compared to the gamma level (209 903 km²), indicating that Amazonian tipping points are more likely to occur locally than regionally or basin-wide. We also observe significant but weak positive linear relationships between tree species diversity and stability at both alpha and beta scales. This emphasizes both the importance of biodiversity conservation at multiple spatial scales and the complexity of understanding the stability of the Amazon forest.

1 Introduction

Tropical forests play a crucial role in the global carbon cycle. While historically serving as substantial carbon sinks, there is growing evidence that certain tropical forests, including the Amazon, are transitioning to carbon sources (Bennett et al., 2023). This shift is likely attributable to a combination of local and regional disturbances, such as extreme drought events and deforestation (Brienen et al., 2015). Given that forest loss in the Amazon can affect the forest's internal water cycle, leading to further forest loss (Zemp et al., 2017) and global climate change (Artaxo, 2023), these recurrent disturbances could push the Amazon to a tipping point, leading to a change in its functional state (Flores et al., 2024). The stability of the Amazon forest to such disturbances can be quantified using changes in temporal autocorrelation at lag-1 (TAC) of a variable representing an ecosystem state (Boulton et al., 2022). TAC measures the correlation between successive time points in the time series, with an increase indicating a greater similarity between the current and previous ecosystem states over time (Dakos et al., 2012). According to mathematical theory, when dynamical systems such as the Amazon lose stability and approach a tipping point, they are expected to show a slowing recovery to short-term disturbances (translated as an increase in TAC) – a phenomenon known as critical slowing down (Boulton et al., 2022; Scheffer et al., 2009). For easier interpretation, we will rather use the term “tipping likelihood” hereafter.

The Amazon forest is one of the most biodiverse regions on the planet, with an estimated 16 000 tree species (ter Steege et al., 2020). These species are not equally distributed along the basin due to the observed gradients of environmental conditions including precipitation, seasonality, soil type and fertility, and exposure to disturbance (Luize et al., 2024; ter Steege et al., 2023). The predicted increase in climate change impacts and deforestation across the Amazon could reduce its tree species richness by fifty percent by 2050 (Gomes et al., 2019), but little is known about how these diversity gradients impact tipping likelihoods in the Amazon forest (Flores et al., 2024; Hutchison et al., 2018). This can be partly attributed to the scarcity of field data on multi-scale tropical tree diversity. However, the emergence of modelled global tree species diversity datasets that cover multiple spatial scales can help to address this (Keil and Chase, 2019; Liang et al., 2022a). Furthermore, similar to biodiversity patterns, tipping likelihood in the Amazon forest might vary significantly from local to regional scales (Keil and Chase, 2019), but this spatial dependency has not yet been quantified (Lenton et al., 2022).

The insurance hypothesis suggests that more diverse ecosystems have a larger buffer against disturbances, owing to the varied species responses to environmental fluctuations (Yachi and Loreau, 1999). Consequently, diversity may increase stability by reducing the likelihood of tipping points. Empirical support for the diversity-stability hypothesis has been observed in both grassland and forest ecosystems (Isbell et al., 2015; Liu et al., 2022), although contrasting results have been found as well (Grossiord et al., 2014). The relationship between diversity and ecosystem stability has primarily been investigated at local scales (Gonzalez et al., 2020), but diversity can be defined on multiple scales. Alpha and gamma diversity refer to diversity at local and regional scales, respectively, while beta diversity represents differences in diversity across local communities (Buckley and Jetz, 2008; Keil and Chase, 2019). In the Amazon, beta diversity is shaped by biogeographic gradients, such as variations in soil fertility (ter Steege et al., 2006), as well as historical factors, including the long history of plant domestication by pre-Columbian peoples (Levis et al., 2017).

Using a similar rationale, ecosystem stability can be expanded from local (alpha stability) to regional metacommunities (gamma stability) (Wang and Loreau, 2016). In this context, beta stability is defined as the spatial variation in stability responses across local communities, hereafter referred to as spatial asynchrony. Multiple stability frameworks exist, and we adopt the one by Van Meerbeek et al. (2021). Here, stability encompasses all system properties determining the size, length, and irreversibility of changes following a disturbance. In other frameworks, this is also referred to as non-local stability or resilience (Dakos and Kéfi, 2022). Combining both frameworks, we anticipate that tree diversity and tipping likelihood in the Amazon will relate across multiple scales (Gonzalez et al., 2020; Wang and Loreau, 2016). Based on the insurance hypothesis, we expect higher diversity to be associated with higher stability, or equivalently, a lower likelihood of tipping, at the alpha scale (Liu et al., 2022). On the other hand, there is no defined universal relationship between beta diversity and ecosystem stability. This relationship rather depends on the existing gradients in abiotic heterogeneity, habitat isolation, and species pool richness (Van Der Plas et al., 2023).

In this research, we use changes in TAC of a satellite-derived proxy of canopy productivity as an indicator of tipping likelihood, which we combine with modelled tree species richness across multiple spatial scales for the Amazon forest. Our objective is to understand the impact of tree species diversity on the likelihood of Amazon tipping points, and quantify how this relationship changes across scales. By examining the interplay between diversity and tipping likelihood across varying scales, our research aims to advance our understanding of the stability of the Amazon rainforest and shed light on the importance of biodiversity conservation in maintaining ecosystem stability in this critical region.

2 Materials and Methods

2.1 Biodiversity across spatial scales

We used published global maps of modelled tree species diversity at both local and regional scales to investigate their impact on tipping likelihoods (Keil and Chase, 2019). They integrated data on tree species richness from 1336 forest plots and 282 countries and other administrative units. Additionally, they included 11 predictors: area, latitude and longitude, biogeographical realm, location on mainland or island, elevation, mean gross primary productivity, mean annual temperature, mean isothermality, seasonality, and precipitation in the driest quarter. All climatic input variables had a spatial resolution of approximately 1 km². Using generalized additive models, they predicted species richness in artificially generated 1 ha plots (alpha diversity plots) and in hexagons spanning 209 903 km² (gamma diversity hexagons) distributed uniformly across the global mainland. Subsequently, 1 ha plots lacking at least one environmental variable and hexagons with less than 50 % of mainland area were excluded, resulting in 9761 alpha diversity plots and 620 gamma diversity hexagons globally. Their model explained more than 90 % of the deviance of the data and showed high predictive power when validated with external data. Additionally, they mapped the ratio of gamma and alpha diversity, denoted as beta diversity, at the 1 ha plot scale. Alpha and gamma diversity represent the predicted overall local and regional tree species richness, respectively. Beta diversity indicates the difference in the diversity patterns between both scales, with higher values indicating that a community is relatively less diverse than its surrounding region.

2.2 Stability across spatial scales

In this study, we define high Amazon forest stability as a low likelihood of reaching a tipping point, using time series of satellite data across the Amazon (as explained below). The delineation of the Amazon rainforest was adopted from Olson et al. (2001). As a proxy for the canopy vegetation productivity of tropical forests, we used the enhanced vegetation index (EVI) (Van Passel et al., 2022). EVI measures the canopy greenness, rather than the woody growth response. Therefore, using EVI as a proxy for ecosystem stability assumes that reductions in stem growth coincide with canopy browning (Janssen et al., 2021). We extracted monthly EVI images from 2001 to 2019 from the daily Moderate Resolution Imaging Spectrometer (MODIS) MCD43C4 product with a spatial resolution of 0.05° (Schaaf and Wang, 2015). EVI is calculated using Eq. (1):

$$\begin{array}{}\text{(1)}& \text{EVI}=\mathrm{2.5}\cdot {\displaystyle \frac{(\text{NIR}-\text{RED})}{(\text{NIR}+\mathrm{6}\cdot \text{RED}-\mathrm{7.5}\cdot \text{BLUE}+\mathrm{1})}},\end{array}$$

where NIR, RED and BLUE are the surface reflectance values of the near-infrared, red and blue MODIS bands, respectively. Pixels with a low Bidirectional Reflectance Distribution Function (BRDF) inversion quality (50 % or more fill values) or classified as outliers were excluded from the time series. Additionally, pixels with over 10 % missing values in their time series were masked. We minimized the potential impact of anthropogenic factors on the stability of the Amazon rainforest by excluding areas that were categorized as non-native forest landscapes (Potapov et al., 2008), with a tree cover below 60 % (DiMiceli et al., 2015), and those subjected to burning from 2001 to 2019 (Giglio et al., 2015).

To assess the regional (gamma) stability response, we first calculated the mean EVI time series for all non-masked 0.05° pixels within the gamma diversity hexagons outlined by Keil and Chase (2019). Subsequently, for each hexagon, the mean EVI time series were detrended and deseasonalized using STL decomposition (seasonal and trend decomposition using Loess; Cleveland et al., 1990), resulting in EVI remainder time series. The calculations for the long-term trend and the seasonality used time windows of 19 and 13 months, respectively (Van Passel et al., 2024a). The decomposed EVI remainder was then used to calculate the lag-1 autocorrelation (one-month lag) with a moving window length of five years. The use of longer time lags could highlight recurring seasonal or annual patterns between vegetation and climatic drivers (Hofhansl et al., 2014). However, we used lag-1 temporal autocorrelation rather than longer time lags to be in line with previous critical slowing down research (Boulton et al., 2022; Forzieri et al., 2022; Runge et al., 2025). The trend in TAC was determined as the slope of the linear regression of the TAC time series from 2001 to 2019, and was used to quantify regional tipping likelihood. For ease of interpretation (i.e., higher values denote higher stability), we used the negative value of the slope in TAC as the gamma stability value. Negative gamma stability values thus indicate an increase in TAC over the 20 year period, indicating decreased stability, and vice versa.

The modelled alpha and beta diversity plots from Keil and Chase (2019) are spatially distributed point values at regular intervals of approximately 82 km. To integrate this data with the pixel-based 0.05° satellite time series, we established a buffer zone surrounding the diversity plots to determine which satellite pixels corresponded to specific diversity values. To assess their robustness, multiple buffer zones of varying sizes were employed. Specifically, buffer zones with radii of 8.5 km (size of 225 km²), 14.1 km (625 km²), and 19.7 km (1225 km²) were utilized to match the area of three by three, five by five, and seven by seven MODIS 0.05° EVI pixels, respectively. Within each buffer zone, all pixels were assigned the same alpha and beta diversity values. In the analysis using the smallest buffer size of 225 km², a total of 1866 EVI pixels within 497 buffer zones were considered. These numbers expanded to 5130 and 9904 EVI pixels in 551 and 591 buffer zones when using the larger buffer sizes of 625 and 1125 km², respectively (Table A1).

Using the same methodology as for the gamma scale, alpha stability was derived as the negative value of the TAC trend in the EVI time series of the pixels included in a buffer zone. Spatial asynchrony (or beta stability) indicates the spatial variability in stability responses within a region. To estimate it, we followed the definition of beta diversity by Keil and Chase (2019), and calculated it on the local scale for all pixels within a buffer zone using Eq. (2):

$$\begin{array}{}\text{(2)}& \text{Spatial asynchrony}={\displaystyle \frac{\text{Gamma stability}+\mathrm{0.1}}{\text{Alpha stability}+\mathrm{0.1}}},\end{array}$$

Given that both gamma and alpha stability values varied between approximately −0.01 and 0.01, we added 0.1 to both stability metrics to ensure that both the numerator and denominator are always positive. A spatial asynchrony value above one indicates lower local stability compared to the stability of the surrounding region, while a value below one indicates a relatively higher local stability compared to the regional stability. This transformation of adding a constant to both numerator and denominator is necessary to avoid having a positive gamma and negative alpha stability value or a negative gamma and positive alpha stability resulting in the same spatial asynchrony value. However, to quantify whether this transformation shifts the ratio away from representing the true proportional relationship between regional and local stability, we also split our dataset into those with strictly positive or negative alpha stability values. For each category, we also calculated spatial asynchrony as the ratio of gamma and alpha stability without adding a constant, and these results are presented in Appendix A.

Apart from increasing TAC, another potential indicator of tipping points is an increase in the variance of the system's state variable (Ditlevsen and Johnsen, 2010; Scheffer et al., 2009). To quantify the trend in variance, we calculated the standard deviation (SD) of the remainder of the decomposed EVI time series with a moving window of five years. We calculated alpha and gamma stability using only the 0.05° pixels that showed a consistent positive or negative trend for both TAC and SD, following the theory that an increase in TAC cannot serve as evidence of increased tipping likelihoods if there is no corresponding increase in variance (Ditlevsen and Johnsen, 2010). We repeated all analyses below using only the subset of pixels with consistent trends to provide a more nuanced understanding of the stability dynamics of the Amazon forest (see Appendix A).

2.3 Environmental heterogeneity

To account for the diverse conditions within the Amazon forest, we incorporated different environmental variables known to impact how tropical forests respond to perturbations. As topography modulates the embolism resistance of tropical tree species (Mattos et al., 2023), we included elevation (from NASA JPL, 2013). Additionally, information on the soil sand and clay content from Poggio et al. (2021) was integrated, as soil texture influences the water infiltration speed and water retention capacity of the soil (Van Passel et al., 2022). Lastly, we included seasonality to describe the intra-annual precipitation variability in the Amazon forest, derived from monthly TerraClimate precipitation data from 1980 to 2019 (Abatzoglou et al., 2018) and calculated using the Seasonality Index (SI; Walsh and Lawler, 1981). Low SI values reflect evenly distributed rainfall throughout the year, while high values indicate that all rainfall is concentrated in a single month. We took the mean value of each environmental variable per 0.05° pixel and 209 903 km² hexagon to include in the following analyses.

2.4 Extreme drought occurrences

Extreme drought occurrences are prevalent perturbations in tropical forests (IPCC, 2021). Therefore, we also included variables characterizing the drought history of each pixel. Extreme drought occurrences were identified spatially and temporally using cumulative water deficit (CWD) anomalies (Van Passel et al., 2022). We used precipitation time series from TerraClimate, which merges the high-spatial resolution climatological normals from the WorldClim dataset with time-varying data with a coarser spatial resolution from the Climatic Research Unit time series data version 4.0 (CRU TS4.0) and the Japanese 55 year Reanalysis (JRA-55) (Abatzoglou et al., 2018). CWD was calculated for all pixels within the Amazon using monthly TerraClimate precipitation (P) time series from 1980 to 2019 and a fixed evapotranspiration (E) of 100 mm per month, using the following rule (Aragão et al., 2007):

$$\begin{array}{c}{\displaystyle }{\text{If CWD}}_{n-\mathrm{1}}-E+P_n<\mathrm{0};\\ {\displaystyle }{\text{Then CWD}}_n={\text{CWD}}_{n-\mathrm{1}}-E+P_n;\\ {\displaystyle }{\text{Else CWD}}_n=\mathrm{0}\end{array}$$

Where n represents each month in the time series. Months with negative CWD values represent periods when monthly precipitation is insufficient to offset both the evaporation occurring that month and any precipitation shortfall carried over from the previous month, indicating a water deficit. The CWD dataset was then used to calculate the mean and standard deviation of CWD per month. Standardized anomalies were determined per 0.05° pixel by subtracting the monthly mean from the pixel value and dividing the result by the monthly standard deviation. Pixels with CWD anomalies below −1.96 had significantly larger water deficits than average (with p<0.05). Extreme drought periods were defined as starting with at least two months of significantly dry CWD anomalies and ending when the anomaly became positive. Since a fixed evapotranspiration value was used, the CWD anomalies in this study represent meteorological droughts rather than hydrological droughts. For each drought event, the intensity was calculated as the absolute value of the minimal CWD anomaly value during the drought period, while the duration was defined as the number of months that the drought lasted. The drought history of each pixel was then characterized using drought frequency (i.e., the number of separate drought events between 2001 and 2019), along with the average drought intensity and duration of all occurred droughts.

2.5 Diversity-stability relationship across spatial scales

All analyses were performed using the R statistical environment (R version 4.2.3; R Core Team, 2024). Using the R package piecewiseSEM (Lefcheck, 2016), piecewise structural equation modelling was used to explore the relationships between environmental factors, drought-related variables, diversity, and stability across various spatial scales, considering each buffer size. We first constructed a conceptual model, where we hypothesized that higher local diversity increases ecosystem stability and therefore lowers tipping likelihood, due to its buffering effect against disturbances (Fig. 1) (Liu et al., 2022; Wang and Loreau, 2016). Furthermore, we anticipated that the spatial variability in stability responses between local communities (i.e., spatial asynchrony) would be driven by the spatial variability in diversity between them (i.e., beta diversity). Lastly, we expected spatial asynchrony to exert a more pronounced positive impact on gamma stability than local stability. This expectation stems from the high spatial heterogeneity and species turnover characterizing the Amazon forest (Keil and Chase, 2019; Qiao et al., 2022).

Figure 1 Conceptual model with hypotheses for all causal relationships in the structural equation model. The dashed boxes represent stability across spatial scales, while the dotted boxes represent the diversity components. The different hypotheses are: (1) A region will be more stable if it contains more stable communities; (2) A region will be more stable due to higher spatial variability in stability responses between communities; (3) Spatial variability in stability responses will be driven by the spatial variability in diversity; (4) A region will be more diverse due to higher spatial variability in diversity between pixels; (5) Soil, climate variability, and drought-related variables will impact local stability and spatial asynchrony; (6) A region will be more diverse if it contains more diverse communities; and (7) A community will be more stable if it is more diverse due to the insurance hypothesis.

The structural equation model incorporated three linear mixed-effects models with an exponential spatial autocorrelation structure to model the diversity-stability relationship on the alpha, beta, and gamma scales, using the R package nlme (Pinheiro et al., 2021). The models met all the assumptions of a linear model, that is, linearity of the data, normality of the residuals, homogeneity of residuals variance, and independence of the error terms. Considering that all alpha stability values within a given buffer zone shared the same alpha diversity value, “buffer zone ID” was included as a random effect in the alpha and gamma stability models. Additionally, for the beta-scale model, where beta diversity and spatial asynchrony were computed using a single gamma value per region, both “buffer zone ID” and “region ID” were included as nested random effects. We also included a generalized mixed-effects model for the diversity relationship across scales in the structural equation model, with “region ID” as a random effect and a Poisson distribution for the gamma diversity data, using the R package MASS (Venables and Ripley, 2002). We included multiple climatic and soil variables known to influence the stability of the Amazon to disturbances, along with variables describing the extreme drought history of each pixel (see above), as an example of a recurring disturbance within the Amazon forest. However, seasonality was not included directly in the alpha stability model due to its high variance inflation factor when combined with alpha diversity. Therefore, we included seasonality as a covariate of alpha stability instead. To test the robustness of the modelled diversity data, we also repeated the alpha-scale diversity-stability analysis using two other publicly available tree species diversity datasets (Liang et al., 2022a; ter Steege et al., 2023), with the results included in Appendix A.

3 Results

When including only the satellite pixels overlapping with the regularly distributed alpha diversity plots, we observed negative stability values, indicating increased tipping likelihood, for 40 % of the included EVI pixels at the alpha scale and for 20 % of the regions at the gamma scale (Figs. 2 and 3). Additionally, more than 70 % of the EVI pixels exhibited a spatial asynchrony value above one, indicating relatively higher tipping likelihoods in the local forest communities than in the larger region where they are situated.

Figure 2 Stability (i.e., the negative TAC trend) and diversity (i.e., modelled tree species richness) for the Amazon region across spatial scales using the smallest buffer size of 225 km². (a, b) Alpha stability and diversity; (c, d) spatial asynchrony and beta diversity; and (e, f) gamma stability and diversity. In (a) and (c), the individual pixels are shown bigger than their real size for increased visibility. The original maps are shown in Fig. A1.

Figure 3 Linear diversity-stability relationships on the (a) alpha, (b) beta, and (c) gamma scale, for the smallest buffer zone of 225 km². The lines in (a) and (b) show significant positive relationships while considering random effects and spatial autocorrelation (p<0.05), while (c) does not have a significant relationship. The histograms on the top and the right side show the distribution of the diversity and stability values, respectively, following the colour scale of Fig. 2. The black lines in the right histograms indicate zero in (a) and (c), and 1 in (c), with the percentage of values above and below.

Consistent with our hypotheses, we found significant but weak positive linear relationships between diversity and stability on both the alpha scale (marginal R²=0.01, conditional R²=0.18) and the beta scale (marginal R²=0.02, conditional R²=0.41) (Figs. 3 and 4, Table A1). However, we did not find a significant relationship on the gamma scale when accounting for spatial autocorrelation (Fig. 3). These significant relationships on the alpha and beta scale persisted when using larger buffer sizes, and when restricting the analysis to only pixels with consistent TAC and SD trends (Tables A1 and A2). We also found a significant positive relationship between alpha stability and diversity for the globally modelled tree species diversity dataset from Liang et al. (2022a), but no significant relationship with the Amazon-based diversity estimates (Fig. A2).

Figure 4 Representation of the structural equation model involving the causal relationships between stability, diversity, environment, and drought history across spatial scales, using the smallest buffer size of 225 km². Model Fisher's C=41.4 (p=0.08). The R²m and R²c values next to the response variables represent the marginal and conditional R² values, respectively. The single-headed straight arrows represent causal pathways, and the double-headed curved arrows represent the covarying variables. The numbers on the arrows represent the significant effect sizes of the standardized path coefficients.

Moreover, we found significant positive effects of both alpha stability and spatial asynchrony on gamma stability, with spatial asynchrony having a slightly larger effect size (Fig. 4). Notably, local environmental conditions and drought history variables did not exert significant effects on alpha stability. In contrast, precipitation seasonality, elevation, and drought frequency did significantly impact spatial asynchrony. When incorporating larger buffer sizes or when restricting the analysis to the pixels with consistent TAC and SD trends, most of these significant relationships remained robust (Figs. A3 and A4). Furthermore, higher soil sand content was found to significantly increase alpha stability in many of these models. When splitting the dataset into negative and positive alpha stability values to avoid the use of the constant in Eq. (2), the larger effect size of spatial asynchrony on gamma stability compared to that of alpha stability remained robust, as did the impact of precipitation seasonality on spatial asynchrony (Fig. A5). We also found positive impacts of alpha and beta diversity on alpha stability and spatial asynchrony, respectively, when only including the positive alpha stability values, although the impact of beta diversity became insignificant for the negative alpha stability values.

4 Discussion

We integrated changes in satellite-derived TAC of canopy productivity with modelled tree species richness to understand the impact of tree species diversity on the likelihood of Amazonian tipping points. By incorporating stability and diversity data across multiple spatial scales, we showed that forest areas characterized by higher alpha tree species diversity exhibit lower tipping likelihoods. More diverse forest communities thus tend to be further away from a potential critical threshold to a degraded state. Additionally, communities with greater dissimilarity in tree species composition (i.e., higher beta diversity) contribute to increased heterogeneity in stability responses across the region, confirming the positive diversity-stability relationship in the Amazon at both alpha and beta scales. This highlights the critical role of biodiversity conservation in upholding ecosystem stability within the Amazon forest. However, the low R² values of the positive relationships also show the complexity of understanding the drivers of ecosystem stability in the Amazon forest.

4.1 Likelihood of Amazon forest tipping points

Our analysis reveals that approximately one-third of the analyzed pixels within the Amazon forest exhibited increased tipping likelihoods over the 20 year study period, which decreases to one-fifth of the Amazon at the regional scale. This aligns closely with the findings of Van Passel et al. (2024a) for the entire Amazon, indicating that the subset of pixels investigated in this study can be considered representative of the broader Amazon rainforest. Additionally, more than 70 % of the EVI pixels displayed a spatial asynchrony value above one, indicating relatively higher tipping likelihoods locally than regionally. This disparity can possibly be attributed to significant variations in environmental conditions, including topography, microclimate and soil types, at the local scale (Ismaeel et al., 2024; Mattos et al., 2023; Zuquim et al., 2023). Such variations likely contribute to greater variability, leading to a higher likelihood of local tipping points. In contrast, the analysis at the regional scale may smooth out these local differences, making broader regional patterns of the region more apparent and resulting in less pronounced regional trends. This observed pattern of increased stability from the local to the regional levels parallels findings in research on drought resistance across spatial scales in the Amazon (Janssen et al., 2020). It underscores the importance of considering various spatial scales when investigating ecosystem stability. Furthermore, this pattern also suggests that the Amazon is more likely to experience local transitions to a degraded state than to reach a regional or system-wide critical threshold (Franco et al., 2025; Rammig et al., 2010).

4.2 Diversity-stability relationship across spatial scales

The significant diversity-stability relationships on the alpha and beta scales highlight the role of tree species diversity and its spatial heterogeneity in protecting the Amazon forest from reaching a tipping point. However, the low explanatory power of these relationships could potentially be due to the high functional redundancy in Amazonian forests. At small scales, additional species can still make a difference in terms of niche complementarity, whereas, at larger spatial scales, they may lead to functional redundancy (Poorter et al., 2015). Addressing the low R² values and the unexplained variability of the TAC trends will require mechanistic modelling to disentangle the underlying causal processes.

The observed positive relationship between local diversity and stability is consistent with the insurance hypothesis. Higher species diversity promotes diversity of responses to various local disturbances, mitigating the overall impact on local communities, and thereby, lowering the chances of the forest to transition into a degraded state. This positive relationship is in line with prior research investigating stability of (sub)tropical forests, where stability was quantified as either a low tipping likelihood, resistance to drought events, or temporal stability of productivity (Hutchison et al., 2018; Liu et al., 2022; Ouyang et al., 2021).

The significant positive relationship between beta diversity and spatial asynchrony indicates that communities with greater dissimilarity in composition contribute to increased heterogeneity in stability responses across the region, a phenomenon referred to as the spatial insurance effect (Wang and Loreau, 2016). While there has been less research on this relationship compared to the alpha scale, especially in tropical forests, empirical evidence in temperate biomes has also demonstrated the existence of the spatial insurance effect (Liang et al., 2022b; Qiao et al., 2022). Beta diversity is especially important in driving ecosystem functioning in abiotically heterogeneous landscapes (Van Der Plas et al., 2023). The Amazon forest is a clear example of such a landscape, with its highly diverse soils shaped by a complex geological history (Hoorn et al., 2010; Tuomisto et al., 2019), and further influenced by variations in precipitation seasonality and local topography that strongly influence local hydrological conditions (Costa et al., 2023; Mattos et al., 2023).

The absence of a significant causal relationship between diversity and stability on the regional scale aligns with previous findings in temperate ecosystems (Hautier et al., 2020; Qiao et al., 2022), but contradicts other findings (Liang et al., 2022b). This lack of significance could be due to the small sample size (Fig. 3c), which could hide any possible real effect due to low statistical power. At the same time, in a meta-analysis of observed drought responses in the Amazon, regional-scale drought responses of productivity were also reported to be of smaller magnitude and significance compared to those on the local scale (Janssen et al., 2020). At the regional scale, aggregating local dynamics often leads to a smoothing of temporal variability, as asynchronous responses among local communities can offset each other. This weakening of the temporal signal, combined with high functional redundancy, may buffer regional-scale ecosystem functioning and thus weaken the apparent strength of diversity–stability relationships. As a result, even substantial changes at the local scale may translate into smaller, less detectable effects when averaged over large areas. This, in combination with the relatively low R² values observed between diversity and stability on both the alpha and beta scales, emphasizes that positive relationships between species richness and stability tend to be weak in tropical forests and diminish with increasing spatial scales (Ouyang et al., 2021; Qiao et al., 2022).

4.3 Environmental and climatic drivers

Our study also revealed significant effects of various environmental and drought-related variables on stability at both the alpha and beta scales. Specifically, more sandy soils were associated with lower local tipping likelihoods, following earlier findings demonstrating higher drought resistance on more sandy soils in the Amazon forest (Van Passel et al., 2022). While this result may seem surprising, given that higher sand content typically correlates with lower nutrient availability and higher aluminium saturation (Laurance et al., 1999), it might be attributed to the higher hydraulic margin of tree species in low-resource environments (Oliveira et al., 2021).

Similarly, spatial asynchrony was significantly influenced by multiple variables. Firstly, it was lower in more seasonal forests within the Amazon, indicating that more seasonal forest communities exhibit relatively lower tipping likelihoods compared to the larger region where they are located. Conversely, forest communities located at higher elevations or experiencing more drought events in the 20 year period showed higher spatial asynchrony. Existing literature suggests negative impacts of all three variables – seasonality, elevation, and drought frequency – on local tropical forest stability (Liu et al., 2022; Mattos et al., 2023; Ouyang et al., 2021; Van Passel et al., 2022). However, our data reveals interesting nuances in how these impacts on the likelihood of Amazonian tipping points change when increasing the spatial scale (Fig. A6). Notably, the impact of higher seasonality becomes more negative (although non-significant) when increasing the spatial scale, which could be attributed to the gradual changes in seasonality over large spatial extents. In contrast, the (non-significant) negative effects of elevation and drought frequency on the local scale, both variables that vary more strongly between local forest communities, become more positive on the regional scale. It is also important to consider that EVI, as a proxy for forest productivity, is not necessarily fully coupled with environmental conditions such as drought or seasonality. For instance, some tree species may increase their productivity under drier conditions due to phenology or responses to strong drought events (Janssen et al., 2021), which could potentially weaken the direct link between these environmental drivers and stability metrics. Furthermore, as the choice of precipitation dataset influences the variability in droughts stress across the Amazon (Papastefanou et al., 2022), using another dataset may yield different causal relationships. TerraClimate showed slightly better validation performance than CRU when compared to weather station data, although its long-term precipitation trends are inherited from CRU (Abatzoglou et al., 2018). A comparison of multiple precipitation datasets over the Amazon found CRU-based drought estimates to be broadly consistent with other products (Papastefanou et al., 2022). Nevertheless, reliance on a single dataset can still lead to regional over- or underestimation of drought occurrences.

These results highlight that the likelihood of Amazonian tipping points is impacted differently by various environmental variables, depending on the spatial scale considered. Most research on the Amazon's response to climate change has focussed either on plot-level studies (Brienen et al., 2015; Esquivel-Muelbert et al., 2019) or coarse-scale satellite time series (Boulton et al., 2022; Van Passel et al., 2024a), often overlooking this scale-dependent aspect of the forest's stability.

4.4 Drivers of gamma stability

Consistent with our hypothesis (Fig. 1), our findings indicate that spatial asynchrony exerts a more pronounced positive impact on gamma stability than local stability. This greater relative effect of spatial asynchrony aligns with prior research investigating stability across a large temperate forest ecosystem (Qiao et al., 2022) and can partly be attributed to the high spatial heterogeneity and species turnover characterizing the Amazon forest (Keil and Chase, 2019).

Our results imply that higher alpha and beta diversity exert indirect stabilizing effects on the regional scale by enhancing alpha stability and spatial asynchrony, respectively. The predominant influence of spatial asynchrony, and thus beta diversity, on gamma stability resonates with studies conducted in temperate forests. These studies similarly highlighted the stronger effects of beta diversity compared to alpha diversity on both ecosystem stability and multifunctionality (Sebald et al., 2021; Van Der Plas et al., 2016). The significance of beta diversity could imply that the high spatial turnover in species composition across the Amazon is closely linked to a high spatial turnover in functional diversity. This, in turn, would lead to local changes in ecosystem functioning that can scale up to large-scale changes in the provision of multiple ecosystem functions, with stability being just one aspect (Mori et al., 2018; Van Der Plas et al., 2016).

4.5 Limitations and perspectives

The tree species richness data used in this study were obtained from published global model predictions rather than direct plot measurements. This approach was necessary due to the absence of spatially comprehensive tree species data covering multiple spatial scales across the Amazon. However, it is important to note that the accuracy of the diversity-stability relationship discussed here is dependent on the accuracy of the diversity predictions. Although the diversity model accounted for more than 90 % of the deviance of the data, exhibited low relative uncertainty in the Amazon forest compared to other global biomes, and performed well in the Amazon when validated with external data (Supplementary Fig. 3 in Keil and Chase, 2019), caution should still be exercised in interpreting these results. The significant positive relationship between alpha stability and the globally modelled tree species diversity from Liang et al. (2022a) supports our findings, but no significant relationship was found with the Amazon-based estimates from ter Steege et al. (2023) (Fig. A2). Although the former dataset demonstrated greater explanatory power for South America compared to that of ter Steege et al. (2023), this further shows the complexity of understanding the drivers of ecosystem stability in the Amazon forest.

Furthermore, we anticipate that the stabilizing effects of diversity in the Amazon forest across spatial scales would be more pronounced when incorporating functional instead of species diversity metrics. Functional diversity, particularly related to hydraulic traits, has been identified as a crucial factor in explaining the variation of drought responses in tropical forests (Barros et al., 2019). Heterogeneity in the vegetation structure has also been identified as a key driver of stability in temperate forest ecosystems (Qiao et al., 2023), but this heterogeneity is more effectively captured through functional than taxonomic diversity. Prior research has noted a gradual shift of Amazon tree communities towards more dry-affiliated species in response to climate change, although delayed due to the extended generation times of tropical trees (Esquivel-Muelbert et al., 2019). These shifts in compositional dynamics might enhance the stability of the Amazon forest against future droughts, but they could also lead to reduced carbon stocks and other ecosystem functions (Rius et al., 2023), and reduced biodiversity, given the dominance of wet-affiliated tree species among Amazonian trees (Esquivel-Muelbert et al., 2017). Hence, considering the temporal dynamics of diversity would offer valuable insights into the diversity-stability relationship in the Amazon. However, data on tree species and functional diversity within the Amazon remain limited, especially on multiple spatial and temporal scales (Keil and Chase, 2019). Recently, remotely-sensed spectral asynchrony, which captures the spatial heterogeneity of species' functional responses across distinct pixels, has been proposed as a readily monitorable metric for assessing the impacts of species diversity in seasonally dry tropical forests (Mazzochini et al., 2024), offering an alternative to ground-based plant diversity data.

Lastly, the use of 20 year EVI-based trends in TAC as a proxy for forest stability has its limitations. Our assumption that the phenology in the Amazon largely explains the remotely sensed EVI patterns was based on previous research (Anderegg et al., 2019; Wu et al., 2022; Wu et al., 2018). However, because EVI does not capture forest structure, changes in EVI TAC do not necessarily reflect increased tree mortality but rather shifts in canopy characteristics. Incorporating additional satellite products that capture other forest characteristics, such as solar-induced chlorophyll fluorescence (SIF) or vegetation optical depth (VOD), would improve the robustness of the results. We also restricted our analysis to regions without fires or land-use change, ensuring that differences in disturbance history between pixels were primarily linked to drought occurrences, which we explicitly accounted for. Nonetheless, since both deforestation and degradation can affect tropical forests up to 100 km away (Araujo et al., 2023; Butt et al., 2023), it is unlikely that this masking completely removed their influence. Finally, while the 20 year timescale might be too short to capture long-term stability trends in the Amazon, the inclusion of three widespread drought years within this period enabled us to quantify how diversity influences medium-term stability responses.

To gain deeper insights into the complex relationships between diversity and tipping likelihoods in the Amazon, there is a need for the expansion and integration of more extensive diversity datasets. Nevertheless, our findings highlight the importance of conserving tree species diversity, rather than just forest cover, in maintaining ecosystem stability within the Amazon.

5 Conclusions

By incorporating stability and diversity data across multiple spatial scales, we found that Amazonian tipping points are more likely to occur locally than regionally or basin-wide. The weak but significant positive relationships between diversity and stability at both alpha and beta scales highlight both the importance of biodiversity conservation and the complexity of understanding and predicting tropical forest stability. In the face of climate change and more frequent extreme drought occurrences across the Amazon, preserving high tree species diversity at multiple spatial scales could act as a valuable buffer against the Amazon transitioning from a carbon sink to a carbon source.

Appendix A

Figure A1 Stability on the alpha and beta scale using the smallest buffer size of 225 km², with the pixels shown at the correct size. (a, b) Alpha stability; and (c, d) spatial asynchrony. (b) and (d) show a zoomed-in detail of the Amazon-scale maps.

Figure A2 Comparison with other publicly available tree species richness estimates across the Amazon forest. Liang et al. (2022a) refers to a globally modelled tree species diversity product with a spatial resolution of 0.025°, based on 1.3 million sample plots (Liang et al., 2022). Their model explained 95 % of the variance in tree species richness in South America. Ter Steege et al. refers to a modelled tree species diversity map of the Amazon with a spatial resolution of 0.1°, based on 2046 plots (ter Steege et al., 2023). Their model explained 71 % of tree species richness in the Amazon. (a) Density comparison of the different tree species diversity estimates in the Amazon. Keil and Chase refers to the main alpha diversity product used in the manuscript. (b) Significant positive linear relationship between alpha diversity from Liang et al. (2022) and alpha stability. (c) Non-significant relationship between alpha diversity from ter Steege et al. (2023) and alpha stability. (b) and (c) show the marginal effect plots of the models including environmental and climatic variables.

Figure A3 Representation of the structural equation model involving the causal relationships between stability, diversity, environment and drought history across spatial scales, using the trend in TAC calculated as lag-1 autocorrelation, for the buffer sizes of (a) 625 km²; and (b) 1225 km². (a) Model Fisher's C=31.1 (p=0.15). (b) Model Fisher's C=24.8 (p=0.21). The R²m and R²c values next to the response variables represent the marginal and conditional R² values, respectively. The full blue and red arrows indicate significantly positive and negative coefficients (p<0.05), while the dotted lines indicate non-significant effects. The single-headed straight arrows represent causal pathways, and the double-headed curved arrows represent the covarying variables. The numbers on the arrows represent the significant effect sizes of the standardized path coefficients.

Figure A4 Representation of the structural equation model involving the causal relationships between stability, diversity, environment and drought history across spatial scales, using only the 0.05° EVI pixels with a trend in TAC and in SD that do not show significant opposite effects, for the buffer sizes of (a) 225 km², (b) 625 km², and (c) 1225 km². (a) Model Fisher's C=36.7 (p=0.13). (b) Model Fisher's C=31.4 (p=0.21). (c) Model Fisher's C=36.6 (p=0.08). The R²m and R²c values next to the response variables represent the marginal and conditional R² values, respectively. The single-headed straight arrows represent causal pathways, and the double-headed curved arrows represent the covarying variables. The numbers on the arrows represent the significant effect sizes of the standardized path coefficients.

Figure A5 Representation of the structural equation model involving the causal relationships between stability, diversity, environment and drought history across spatial scales, including only (a) negative alpha stability values, and (b) positive alpha stability values. In (a), higher spatial asynchrony values represent a higher regional than local tipping likelihood, which is opposite from the interpretation in the main manuscript, hence the switched signs of some of the relationships. In (b), higher spatial asynchrony values represent a higher local than regional tipping likelihood, similar to the interpretation in the main manuscript. The R²m and R²c values next to the response variables represent the marginal and conditional R² values, respectively. The single-headed straight arrows represent causal pathways, and the double-headed curved arrows represent the covarying variables. The numbers on the arrows represent the significant effect sizes of the standardized path coefficients.

Figure A6 Linear relationships between stability and (a, b) seasonality; (c, d) elevation; and (e, f) drought frequency, on the (a, c, e) alpha and (b, d, f) gamma scale. Spatial autocorrelation was included using an exponential correlation structure in all the models. Dashed and full lines indicate non-significant and significant linear relationships (with p<0.05), respectively. Elevation, seasonality and drought frequency on the gamma scale were calculated as the mean values of all EVI pixels included in the smallest buffer zone of 225 km².

Table A1 Results of the linear diversity-stability relationships on the alpha and beta scale, using the trend in TAC for all EVI pixels, for the three buffer sizes with the 0.05° pixels. R²m and R²c are the marginal and conditional R² values, respectively. Significant effect sizes are indicated with asterisks (${}^{\ast \ast \ast }$ p<0.001; ${}^{\ast \ast }$ p<0.01; ^∗ p<0.05).

Table A2 Results of the linear diversity-stability relationships on the alpha and beta scale, using only the EVI pixels with a trend in TAC and in SD that do not show significant opposite effects, for the three buffer sizes with the 0.05° pixels. R²m and R²c are the marginal and conditional R² values, respectively. Significant effect sizes are indicated with asterisks (${}^{\ast \ast \ast }$ p<0.001; ${}^{\ast \ast }$ p<0.01; ^∗ p<0.05).