Tree cover change proves stability and instability in tropical ecosystems

1 Stockholm Resilience Centre, Stockholm University, Stockholm, Sweden 2 Bolin Centre for Climate Research, Stockholm University, Stockholm, Sweden 3 Department of Water Management, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Delft, The Netherlands 4 Department of Physical Geography, Faculty of Geosciences, Utrecht University, Utrecht, The Netherlands * Author to whom any correspondence should be addressed Email: chandrakant.singh@su.se


Main
Climate change and deforestation reduces the resilience of rainforest ecosystems 1,2 , and thus compromise their capacity to remain forests despite various perturbations 3,4 . Resilience can be quantified and analysed by constructing a 'stability landscape' (Fig. 1), in which valleys ('basins of attraction') represent 'stable states' and hilltops represent 'unstable states' under transition. Resilience is then measured as the width of the basin of attraction around a stable state, which erodes towards bifurcation points (i.e., a point where stable and unstable states collide, becoming unstable) 1,2 (Fig. 1a). Within a basin of attraction, stabilising feedbacks help the ecosystem retain its structural and functional characteristics against perturbations 5 .
The ecosystem will eventually return to its native stable state ('minimum' of the basin) when perturbations on the system are released (Fig. 1b,c). Beyond a basin of attraction, i.e., trespassing a threshold ('maximum' of the basin), self-amplifying feedbacks will instead shift an ecosystem to an alternative stable state 1,5 . A better understanding of stability and resilience is helpful to evaluate the potential of ecosystem adaptation and systemic risks under future (climatological or non-climatological) modifications to their conditions 6 .
Due to the lack of analysis of dynamics through time series 7,8 , our present understanding about the stability landscape of the terrestrial tropical ecosystems is based on the frequency distribution of tree cover 1,[9][10][11] , essentially making a space-for-time assumption (Fig. 1a).
According to this methodology, frequency distribution determines the size (i.e., width and depth) of the basin of attraction in the conceptual stability landscape, which is then interpreted to be ecosystems stability (deep basin, more stable and vice versa) and resilience (wider basin, more resilient and vice versa) across time 1,12 (Fig. 1a). The availability of longer time series of remote sensing data now allow for a better representation of these ecological states and resilience across time 8,13 .
Here, for the first time, to our knowledge, a time series of tree cover spanning two decades is analysed to investigate rainforest stability and resilience. It is known that the ecosystem's response towards any perturbations should be captured in the transient state (i.e., as tree cover change (ΔTC)) of the ecosystem 14,15 . We hypothesise that the transient state of the ecosystem should resemble the stability landscape found by the space-for-time assumption (Fig. 1a).
Thus, a highly resilient ecosystem will not show considerable ΔTC over time, whereas a lowly resilient ecosystem will. is, originally, based on the frequency distribution of the tree cover (space-for-time assumption 1,[9][10][11] ). This study substitutes 'tree cover frequency' with magnitudes of 'tree cover change over time' for South America and Africa (spatio-temporal) across different classes of precipitation, which we hypothesise should resemble the original landscape. Stable and unstable states (i.e., equilibria) correspond to the valleys (i.e., local minima) and hilltops (i.e., local maxima) in the stability landscapes, respectively. (b,c) Resilience of an individual ecosystem across the stability landscape is represented as the width of the basin of the attraction around a stable state, which declines towards the bifurcation points (i.e., a point where stable and unstable states collide; depicted in a,b). Perturbations push the ecosystem towards the hilltop, whereas the ecosystem returns to its stable state when these perturbations are released.
Our hypothesis suggests a correlation between ΔTC and resilience of the ecosystems.
Previous research overlooks any such correlation and only considers the hydroclimatespecifically mean precipitation ( ) -when quantifying forest resilience 1,16 . Recent insights, however, hint towards the necessity to also incorporate the buffering capacity of the forest ecosystems, an aspect that is often lacking when representing the ecohydrology of tropical ecosystems 13,17,18 . To include this buffering capacity, we use root zone storage capacity (S r ) representing the maximum amount of the subsoil moisture available to ecosystems buffering water shortage during dry periods 17,19 to quantify the resilience of the ecosystem. This aspect acknowledges that ecosystems respond to water-stress by actively investing in their aboveand below-ground structures to maximise their hydrological benefits 17 . Thus, the resulting resilience metric, by also explicitly considering the ecosystems' adaptive and buffering strategies, should be consistent with actual ΔTC.

Tree cover change in relation to stability equilibria
Our spatio-temporal analysis consistently shows low ΔTC for ecosystems at both high (>75%) and low (<10%) tree cover, whereas high ΔTC is observed for ecosystems at intermediate (30-60%) tree cover (Fig. 2a,c). A low ΔTC for both high and low tree cover ecosystems can be the result of either a minimal perturbation on the ecosystem over the last two decades (2000-2019), or a robust adaptive mechanism that is able to offset the experienced perturbations without considerable change in the ecosystem structure 17 , which we, therefore, perceive as 'stable'. Conversely, a high ΔTC at intermediate tree cover (Fig.   2a,c) implies that the ecosystems in these ranges have been potentially influenced by either strong perturbations 20 (e.g., deforestation) causing significant changes to their ecosystem structure, or the adaptive mechanism has modified the ecosystem structure to efficiently utilise available resources 17 (e.g., tree mortality under water stress to make more moisture available for the rest of the ecosystem, or tree growth under the influence of wetter climate 21 ), which could result in them undergoing the observed regime shift 1,12 . Thus, we consider such ecosystems with relatively high structure changes (i.e., ΔTC) as 'unstable'. These spatio-temporal patterns against different levels (Fig. 2a,c), therefore, further strengthen the presence of stability and instability, which previous studies observed using a space-for-time assumption 1,[9][10][11] , can also manifest as actual ΔTC over time across the broader tree cover structures.
A closer look at the alternative stable states (i.e., stable-low and -high tree cover bins representing a series of numerical ranges highlighted in dark brown and green, respectively, in Fig. 2a,c and spatially highlighted in Fig. 2b,d) reveals certain dissimilarities across the precipitation classes. Stable-high tree cover bins decrease gradually with decreasing ( Fig.   2a,c), therefore, implying the inability of the forest ecosystems to maintain their dense structural characteristics under drier conditions 17 , which makes them undergo a self-propagating shift to a savannah-like open-canopy structure 1,22 . Reversely, stable-low tree cover bins decrease with increasing (Fig. 2a,c). Here, increase in wetter conditions in the ecosystem helps suppress fire and prevents fire-driven seedling mortality 22 , and drives more soil water storage under a wetter climate 23 , thus promoting forest growth and colonisation 1,24 .
Nevertheless, the shifting potential, in both these cases, generally manifests itself as a relatively high ΔTC within the stable extent (e.g., relatively high ΔTC for < 985 mm yr -1 for South America, and < 1,468 mm yr -1 for Africa at a tree cover > 70% in Fig. 2), with some exceptions ( Supplementary Fig. 2). . The total samples are equally divided into four classes. Here, each individual bin within each class corresponds to 2,500 samples. From all these bins, the one with the least ΔTC is considered stable, and the one with the most ΔTC is unstable (see methods). These stable bins at low (dark brown) and high (dark green) tree cover in (a,c) are spatially plotted in (b,d). The unstable bins on either side of ΔTC = 0 correspond to tree cover loss (red) or gain (blue). The relative surface area in (a,c) represents the portion of total sample surface area (separately for South America and Africa) on either side of the ΔTC = 0. The tree cover extent of stable and unstable bins in (a,c) are spatially plotted in (b,d). The white regions in (b,d) correspond to excluded land cover ( Supplementary Fig. 1d).
Interestingly, we also observe that for most of the classes, the extent of the unstable bins (i.e., ranges highlighted in red and blue in Fig. 2a,c) is almost similar for both tree cover loss and gain segments. This equal potential for both tree cover loss and gain at intermediate tree cover was already hypothesised in a space-for-time based approach 1 (Fig. 1a), and is confirmed by observable evidence at field scale (e.g., forest loss under increasingly drier conditions 25 or increasing forest growth under El Niño-southern oscillation influenced wet conditions 21 ). Our spatio-temporal approach provides empirical evidence to this ΔTC potential at continental scales, as well as proof that the ecosystem change leading to a regime shift is indeed intensified at intermediate tree cover 1 (Fig. 2). This change in tree cover structure across different levels, thus, seem to follow our spatio-temporal hypothesis ( Fig.   1a and 2).
But why can forest ecosystems maintain stability at different levels and how does that relate to ΔTC (Fig. 2)? The results from our S r analysis show that forest ecosystems maintain their tree cover structure at decreasing by increasing investment in their subsoil structure 17 ( Fig. 3). When going from high to low tree cover, we clearly observe a steep increase in S r with decreasing within the stability extent of tree cover from 85 % to 75% (Fig. 3a) in South America. For Africa, although only a small portion of the forest is in this comparatively low S r stable state (Fig. 3b), we still observe a steep increase in S r near the stable-high tree cover state. Here, the least ΔTC within this stability extent of stable-high tree cover ecosystems' reveals that forests respond to the change in by investing in their S r ( Fig.   2 and 3). This S r investment in reality is the vertical and lateral growth of roots which allows for more subsoil moisture storage, thereby assisting the forest ecosystems' to maintain their (stable) dense tree cover structure even under hydroclimatic stresses 17 . However, this balancing feedback of S r investment to keep the ecosystems in a stable-high tree cover state starts to change as we move to the intermediate tree cover.
At (unstable) intermediate tree cover, we find ΔTC to gradually increase and maximise around 40-50% tree cover (Fig. 2). We also find that the steep increase in S r gradually maximises around the 70-60% tree cover and remains unchanged between 60-30% tree cover ( Fig. 3), thus suggesting causation between maximum S r investment and changes to ecosystem structure 17 . When analysing the changes to the forest ecosystems' structure against varying levels of drought and fire stress at local scale (Supplementary Fig. 3

and 4; see
Supplementary Methods), we observe that unstable state forests -in comparison to stable-high tree cover ecosystems -have often maximised their S r investment and show deterioration to a savannah-like state. These deteriorations are not sudden but gradual over time, thus, suggesting that there exists a certain maximum investment potential beyond which the shift from forest to a savannah state becomes eminent 17 , which manifests itself as relatively high ΔTC over time for the unstable forest ecosystems (Fig. 2, Supplementary Fig.   3 and 4). Considering S r along with , therefore, has allowed us to evaluate the invisible buffering responses of forest ecosystems specifically catered towards efficiently optimising the available water resources and modifying their tree cover structure, and thus is able to manifest the shifts between the transient (stable and unstable) states as different magnitudes of ΔTC.

Resilience of the rainforest
Resilience based on logistic regression predicts the probability of the occurrence of a forest ecosystem (tree cover > 50%) as a function of both and S r for respective continents (Supplementary Table 1  The +S r -based resilience metric shows that the resilience of a large portion of the Congo rainforest is higher than previously presumed (based on only) 1,16 , whereas the resilience of Amazon rainforests shows minor differences ( Fig. 4 and Supplementary Fig. 5). Due to the unique evolutionary history of their respective ecology and climatology 26 , high wet-season precipitation has allowed for Amazonian rainforests species to have a larger subsoil storage (i.e., S r ) to buffer the water-deficit than the Congo rainforests 23,27 . The grass species in Congo rainforests, on the other hand, have evolved to be highly water-efficient 28 , which reduces the competitiveness between trees and grasses for moisture uptake 17 , thereby increasing the resilience of the overall rainforest ecosystem, even with low S r , against water-deficit.
Including S r in our resilience metric, therefore, has allowed us to capture this grass species-induced drought coping strategy in Congo rainforests, which otherwise is hard to detect with just . Nevertheless, the resilience of both rainforest ecosystems are both declining due to increasing regional climatic risks 29 , and combined feedbacks from local deforestation and human-induced fires 3,4 .
Validation with actual ΔTC shows that the +S r -based resilience estimates perform better than only the -based resilience (Supplementary Fig. 6). This performance based on ΔTC further strengthens our original hypothesis that the lower the ΔTC, the more resilient the ecosystem and vice versa ( Fig. 1a and Supplementary Fig. 6). Although is an important variable defining the broad influence of the water cycle of the ecosystem, considering S r as well includes local-scale ecosystem adaptation 17 , which is essential for forests to buffer and withstand hydroclimatic changes, and is thus able to better represent the resilience of the rainforest ecosystems (Supplementary Table 1). This better representation of the ecosystem resilience can, therefore, play a crucial role for management and conservation efforts 30 . Here, a value of one implies a forest ecosystem with highest resilience, and zero implies a forest ecosystem with lowest resilience. Comparing the two resilience metrics, we observed that by considering only (in resilience calculation), the resilience estimates show considerable differences for the Amazon and Congo rainforests (exact difference in Supplementary Fig. 5). Regions with tree cover ≤ 50% and human-influenced land use (see methods) are masked.
We conclude that our observation-based spatio-temporal approach analyses ΔTC over the last two decades and provides empirical evidence of alternative stable and unstable states in the tropical terrestrial ecosystem of South America and Africa. We observe low ΔTC for the ecosystems at >75% and <10% tree cover, which we define as stable ecosystems. In contrast, high ΔTC manifests itself at intermediate tree cover of 30-60%, which we deem as unstable as ecosystems in these ranges are undergoing regime shifts. The tree cover ranges of stability and instability, thus, resembles the stability landscape of the previous space-for-time substitution based approach. Analysing spatio-temporal patterns with S r suggests a trade-off between stability and ΔTC, which reveals that forest ecosystems maintain stability by investing in their subsoil resources. However, maximising this investment potential leads to forest ecosystems undergoing regime shifts -manifesting as ΔTC -towards oncoming perturbations, which otherwise can not be easily explained by only -based space-for-time approach. Only by modifying the existing, commonly used, -based resilience metric with an extended +S r metric, we consider both the influence of hydroclimate (i.e., ) and the adaptive capacity of the ecosystem (i.e., S r ) in defining the resilience of the rainforest ecosystems. Comparing the two resilience metrics, we find that the previous metric underestimates the resilience for a large part of Congo rainforests. Furthermore, the +S r resilience metric shows better consistency with actual ΔTC, thus strengthening its performance over the -based metric. Overall, this study accounts for the ecosystems temporal and adaptation dynamics which are becoming increasingly important to assess the

Study area
This paper focuses on the tropical ecosystems of South America and Africa, but the whole study area is slightly larger: 12°N-50°S for South America and 20°N-35°S for Africa. We have used a global administrative database from the Food and Agriculture Organisation (FAO; http://www.fao.org/geonetwork/) to define geographical boundaries for each country and do not have any political intentions behind our research.

Data
We used remotely-sensed gauge-corrected precipitation and evaporation data for our analysis. While i and ii were obtained at 0.5°, iii was obtained at 0.083°spatial resolution. All three evaporation datasets were obtained at a monthly timescale for the years 2001-2012. We downscaled these datasets from monthly to daily timescale using the daily estimate of the ERA5 36 evaporation at 0.25° spatial resolution.
The above-ground structure of the ecosystem was analysed using the remotely-sensed Ultimately, all the mentioned above datasets were spatially interpolated to 250 m using bilinear interpolation, except for the land-use dataset, which was interpolated using nearest-neighbour interpolation.

Spatio-temporal approach for determining ecosystem states
To evaluate these stable/unstable states, a sample size (n) of 1,000,000 pixels each -from both continents -from all the 250 m 250 m pixels was chosen, and analysed separately for × South America and Africa. This sample was used to determine the tree cover change (ΔTC) in the ecosystem structure in the last two decades as follows:  Fig. 2a and c), with each class representing 25% of the total land area. We further separated each of these classes into samples of tree cover gain (i.e., ΔTC ≥ 0) and tree cover loss (ΔTC < 0) (Fig. 2a and c).
After classifying, we binned the samples into separate bins sorted by mean tree cover (i.e., ; Fig. 2a and c), such that each bin represented an equal area (i.e., 2500 sampled 2000−2019 pixels = 156.25 km 2 ). Lastly, to relate stable and unstable states with the ecosystem's structural change, the 13.4% of the bins with the highest change (i.e., highest ΔTC median ) from all the classes combined were categorised as unstable). Furthermore, 38.2% of the bins with lowest change (i.e., lowest ΔTC median ) were classified as stable. The justification behind selecting the stable and unstable regions was based on the area under the distribution curve ( Supplementary Fig. 7).
Finally, these stable and unstable regions, which were analysed separately for tree cover gain and loss pixels at each class were plotted spatially at 250 m resolution (Fig. 2). For example, our sample analysis suggests that the unstable region for tree cover loss in South America at class of 0-985 mm yr -1 falls approximately between 40-60% ( Fig.   2000−2019   2a). This will be spatially plotted in reality (population) for all the pixels falling between 40-60% at of 0-985 mm yr -1 for the pixels where ΔTC < 0.

Root zone storage capacity
For our analysis, we have considered root zone storage capacity (S r ; derived from daily precipitation and evaporation data) to represent the intrinsic capacity of the ecosystem to absorb and adapt to water-stress conditions (defined here as a deficit in soil water availability inhibiting plant growth). S r is the maximum amount of available subsurface moisture that vegetation can store and utilize through their roots for transpiration during dry periods (i.e., periods in which evaporation is greater than precipitation, irrespective of the seasons) 17,19,38 .
Plants can increase their S r by expanding their roots in the soil laterally as well as vertically.
We adopted the mass-balance approach by ref. 17 to derive S r from precipitation and evaporation estimates (Supplementary Method-1 in Supplementary Information). The underlying assumption of this approach is that ecosystems would not invest in expanding their storage capacity more than necessary to bridge the water-deficits it experiences 19 .

Forest resilience and validation
We adapted ref. 1 Supplementary Information). The logistic regression predicts the probability of forest (tree cover > 50%) as a function of the independent variable. Ref. 1 had only considered as the independent variable. However, the new resilience metric proposed in this study also considered S r as an independent variable representing the drought buffer capacity of the forest ecosystems. Here, we experimented with and S r independently and its combination, and chose the best performing model to represent the ecosystem state (Supplementary Table 1). We modified the S r values for all the regions with tree cover < 30%

methodology for determining resilience using logistic regressions (Supplementary Methods-3 in
to 99 th percentile of each continent's S r . This is because we assume that forest ecosystems will maximise their capacity (i.e., maximise S r ) before transitioning to a savannah 17 . Lastly, we validate the resilience estimates of and +S r combination for both tree cover loss and gain samples separately against observed ΔTC to assess the performance of our metric using spearman rank correlation. A high positive or negative spearman correlation would indicate a high strength and consistency between the resilience estimates and ΔTC.

Tree cover change proves stability and instability in tropical ecosystems
Where T is the drought return period (i.e., 20 years in this study), is the mean annual accumulated deficit , for the years 2001-2012, is the standard deviation of the sample. Also, is the reduced mean and S n is σ −1 the reduced standard deviation, which for n = 11 years equal to 0.4996 and 0.9676, respectively 39 .

Supplementary Method-2: Perturbation trends
Ecosystems lose their structural integrity due to both climatological and non-climatological factors (i.e., human-influenced). In this study, we only focus on environmental changes that are detectable using remotely sensed datasets. To analyse this, we considered two variables (i) fire, and (ii) drought severity. Fire modifications to the ecosystems were analysed using a global time series of burned areas (named FireCCI51; km 2 ) derived from Moderate Resolution Imaging Spectroradiometer (MODIS). The data was procured for the years 2001-2019 at a resolution of 250 m x 250 m 41 . ESA Globcover dataset was used to remove pixels with human land use and non-terrestrial land cover from FireCCI51 dataset.
To analyse the drought severity, we used standardised precipitation evaporation index (SPEI) 42 . SPEI is the modified extension of the standardised precipitation index. It is derived from climate-based datasets and has been widely used to determine the influence of droughts with varying magnitude and duration on the natural and human-influenced systems 6 . SPEI integrates both rainfall (i.e., supply) and potential evaporation (i.e., demand) to capture climate trends (positive for wet and negative for dry climate). For our analysis, we experimented with different SPEI potential evaporation equations (Penman-Monteith 43 and Thornthwaite 44 ). However, we did not observe any significant differences between the general SPEI trends. For our study, we used present the SPEI-12 (i.e. 12-month rolling average SPEI) data derived from the Climatic Research Unit (CRU) dataset, where potential evaporation was based on the Penman-Monteith equation. The SPEI-12 data was directly procured for the years 2000-2018 at 0.5°latitude 0.5°longitude resolution. SPEI-12 algorithm aggregates precipitation and × potential evaporation over a period of 12 months, and is thus convenient for evaluating annual (consistent with tree cover dataset) implications of droughts.

Supplementary Tables
Supplementary Table 1 | Performance of logistic regression models considered for determining the resilience of the tropical ecosystem. The models are evaluated using (shown in increasing order of) the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC). The model with least AIC and BIC is used in this study (Fig. 4). Supplementary Table 2 | Parameters of the logistic regression that predicts the probability of tree cover as a function of and S r . Tree cover > 50% is considered as a high tree cover ecosystem, and tree cover ≤ 50% is considered as a low tree cover ecosystem. All parameters are statistically significant (p-value < 0.05).

Coefficient South America Africa
(Intercept) a 0.  Fig. 1a and 2). We observe that the (a) high tree cover ecosystems (not-stable; tree cover > 65% at < 985 mm yr -1 in Fig. 2a) at the foot of the Andes (marked in red) is fed not only by local , but also gets substantial moisture from glacial runoff 46,47 . These ecosystems are, therefore, able to sustain themselves even at low precipitation. (b) On the other hand, low tree cover ecosystems (not-stable; tree cover < 40% at > 1,834 mm yr -1 in Fig. 2a) between (lowland) Colombia and Venezuela (marked in red) suffers from low nutrient soil characteristics, rapid leaching and an extreme precipitation seasonality, thus promoting a savannah ecosystem even at high precipitation 48,49 . Fig. 3 | Influence of climate and fire in influencing the tree cover of the ecosystem at different regions of interest (ROI) for South America. The regression corresponds to the changes in mean annual tree cover values in the respective ROI. The shade around the regression line corresponds to the 95% confidence interval. Here, near zero slope values signify no tree cover change, positive slopes signify a tree cover gain, and negative slopes signify a tree cover loss over time. The 12-month standardised precipitation and evaporation index (SPEI) signifies progress of the wet (positive; blue) and dry (negative; red) climate in the region. The black lines -on the top -evaluates the total burnt area due to the influence of fire in the respective ROI.  Fig. 2). The dots in blue (i.e., ΔTC ≥ 0) and red (i.e., ΔTC < 0) correspond to our ( and S r ) resilience metric, whereas dots in black correspond to -derived resilience estimates. The shaded regions represent the 1 st and 3 rd quantile. The statistical test calculates the spearman rank correlation (Sp. corr.) coefficient with associated p-value.

Supplementary
Supplementary Fig. 7 | Tree cover change (ΔTC) distribution for South America. The portion within the grey and outside the brown lines represents the stable and unstable regions, respectively. Here, P(z) refers to the probability distribution of the curve, μ is the mean (i.e., -0.26), and σ is the standard deviation (i.e., 10.81) of the distribution. ΔTC for Africa shows a similar trend, however, is not shown in this study.