Reduced resilience of terrestrial ecosystems locally is not reflected on a global scale

Global climate change likely alters the structure and function of vegetation and the stability of terrestrial ecosystems. It is therefore important to assess the factors controlling ecosystem resilience from local to global scales. Here we assess terrestrial vegetation resilience over the past 35 years using early warning indicators calculated from normalized difference vegetation index data. On a local scale we find that climate change reduced the resilience of ecosystems in 64.5% of the global terrestrial vegetated area. Temperature had a greater influence on vegetation resilience than precipitation, while climate mean state had a greater influence than climate variability. However, there is no evidence for decreased ecological resilience on larger scales. Instead, climate warming increased spatial asynchrony of vegetation which buffered the global-scale impacts on resilience. We suggest that the response of terrestrial ecosystem resilience to global climate change is scale-dependent and influenced by spatial asynchrony on the global scale. Climate change negatively influences terrestrial ecosystem resilience at a local scale but also enhances spatial asynchrony which helps to stabilize ecosystems at a global scale, according to statistical analyses of monthly satellite vegetation data.

G lobal changes are likely to alter the structure and function of Earth ecosystems and may induce critical transitions of systems from one stable regime to another 1,2 . Such nonlinear transitions have been identified in diverse ecosystems, such as lakes, marine areas, and terrestrial forests (e.g., savanna and tropical forest as alternative biome states) [3][4][5] . On a larger scale, it is suggested that the planet Earth systems in the Anthropocene may be tipped into an irreversible "Hothouse Earth" state by small perturbations when resilience decreases 6,7 . Assessing the ecological resilience of Earth systems in the current context of global change is thus one of the major challenges of human society, as it is important for developing adaptive strategies and identifying regions with protection priority.
There are two different definitions of resilience. One widely accepted concept of resilience is the degree to which an ecosystem can tolerate changing driver variables and/or disturbance without shifting a quantitatively different state 8,9 . Another definition of resilience is also called engineering resilience, which is the rate of return of the system to its stable state or stationary distribution after disturbance 10 . The two definitions reflect two aspects of ecosystem stability, one focuses on maintaining the existence of function whereas the other focuses on maintaining recovery rate of function 11 . The critical difference of these two different views lies in assumptions regarding whether alternative stable state exist. Engineering resilience assumed that only one stable state exists such that resilience can be directly measured as the recovery rate after external perturbations. For instance, Li et al. evaluate the engineering resilience changes of gymnosperm to extreme drought in the past century using global tree ring records 12 . However, natural ecosystems may have multiple stable states within a specific range of external conditions 3,5 . Therefore, all uses of resilience in this study relate to the first definition.
Direct and precise measurement of ecological resilience is difficult as it requires to estimate the potential drivers and disturbance regimes to move a system across thresholds and boundaries separating alternative domains. In recent years, ecologists proposed statistical early warning indicators (EWIs) based on the phenomenon of "critical slowing down" (ecosystems are more sensitive to perturbations and show slow recovery from perturbations when a catastrophic regime shift is approached) to indirectly depict the ecological behavior, which is characterized as increasing memory, increasing variability, and flickering [13][14][15][16][17][18] . This method generally quantified as the correlation between the resilience indicator and time, with strong increasing trends indicating a loss of resilience. As different individual metric may have different predictive ability, combination of several indicators into a composite indicator could be more reliable inferring regime shifts than a single leading indicator 5,17,19,20 .
The factors controlling the resilience of a terrestrial ecosystem have been shown to be associated with climate change, soil type, number and types of disturbances, and species characteristics 21,22 . Climate change strongly influences the structure and function of global terrestrial ecosystems 23,24 ; for example, elevating atmospheric CO 2 concentration may increase the likelihood of a vegetation shift from a C 4 -to a C 3 -dominated state 25,26 . However, how terrestrial ecosystem resilience changes in the current context of climate change is still not well understood. In addition, previous studies have separately stressed the response of ecosystem functioning and resilience to climate change based on climate variability 24,[27][28][29][30] or changes in mean climate state 31 , with the relative importance of driving forces (e.g., temperature vs. precipitation) and influencing forms (e.g., mean climate state vs. variability) receiving much less attention.
Previous studies have mapped the global distribution of local resilience of terrestrial ecosystems to climate change 24,31,32 . However, changes in ecosystem resilience at the local scale may be inconsistent with those at the global/regional scale because asynchronous or synchronous variations in terrestrial ecosystems among regions may inhibit or enlarge the variation at the global scale 33 . Increased spatial synchrony suggests that the abundance of geographically disjunct populations tends to fluctuate in a coherent way, leading to increased variability in ecosystem functioning at the regional scale [34][35][36] . The key concern is that if critical transitions of these sensitive regions occur simultaneously, the whole ecosystem may trap into an alternative state, leading to a large amount of ecological or economic losses 37 . Therefore, tracking the historical changes in terrestrial ecosystem resilience from local to global scales can help to better predict the synchronized losses of ecosystem services in the future.
Normalized difference vegetation index (NDVI) is commonly used as a measure of gross primary production of terrestrial ecosystem at large scales [38][39][40] , which enables us to study the resilience of global vegetation. In this study, we calculated the composite EWI (CEWI) using monthly NDVI data at local (i.e., each pixel) and global (i.e., across the world or within a biome), with increasing value of CEWI indicating decreased ecological resilience of terrestrial ecosystems. Moreover, we explored the relative importance of climatic factors (temperature vs. precipitation) and influencing forms (climate mean state vs. climate variability) on vegetation resilience in different biomes. In addition, temporal changes of spatial synchrony on vegetation resilience across the world had received much less attention, despite the importance of synchronized change in vegetation abundance enlarged the threat of population extinction 41,42 . We thus evaluated how climate changes influence spatial synchrony in the past three decades to clarify the potential discrepant results of vegetation resilience from local to global scales.

Results
Map of high-resolution ecosystem resilience. Overall, the CEWI increased in 64.5% of terrestrial vegetated grids (Fig. 1a). The spatial pattern of Kendall's τ varied considerably with latitude and elevation, with the fastest increase in the high latitudes of the Northern Hemisphere (Fig. 1b). Specifically, ecosystems in North eastern America, Eastern Siberia, Europe, and Australia experienced a stronger increase in CEWI.
The results of zonal statistics showed that the medians of τ for all biomes (or across the world) were all significantly larger than zero (Fig. 1c, p < 0.001), indicating decreasing resilience at the pixel level during the study period. Further, the medians of Kendall's τ were significantly different among biomes (p < 0.001). Tundra had the largest value, followed by deserts and xeric shrubland, temperate grassland, and montane grasslands.
Relative importance of climatic drivers and influencing forms on ecosystem resilience. We explored the relative importance of temperature (i.e., mean temperature and coefficient of variation (CV) of temperature) and precipitation (i.e., mean precipitation and CV of precipitation) in CEWI changes at the local scale (Fig. 2). For different biomes, temperature was a more important explanatory variable than precipitation (p < 0.001) except boreal forest and temperate grassland (Fig. 2e, p < 0.001 for the statistic test with the null hypothesis that "the relative importance of temperature is equal to precipitation"). The temperature changes contribute the most variation of CEWI in temperate broadleaf and mixed forests, followed by tundra, montane grasslands, and tropical and subtropical moist broadleaf forests (Fig. 2e). The change in climatic mean state (i.e., mean temperature and mean precipitation) rather than variability (CV of temperature and CV of precipitation) was the main cause of decreased resilience in most parts of the world (Fig. 3). Tropical and subtropical moist broadleaf forests, montane grasslands, and deserts and xeric shrubland were the most sensitive biomes to changes in mean climate state (Fig. 3e).
Our results showed that changes in the mean climate state may be a key driver of decreased terrestrial resilience at pixel scales, with more than 50% of the world's areas showing the same direction of climate changes (changes in mean temperature and precipitation) and trends of CEWI (located in the first or third quadrant) (Fig. 4a, c). Specifically, 55.8% of terrestrial regions showed simultaneously a positive trend for mean temperature and a positive trend of CEWI (Fig. 4a).
Ecosystem resilience and spatial synchrony at the global scale. Using the averaged NDVI across the world or biomes, we carried out a similar analysis at the global scale. The global CEWI of terrestrial ecosystems decreased significantly (τ = −0.137, p < 0.001), indicating that there was no evidence of reduced resilience of global terrestrial ecosystems in the past 35 years (Fig. 5). This could be also evidenced from Fig. 4a, c, which showed that climate warming (or increasing precipitation) lead to negative trends of CEWI at the global scale (red points in the fourth quadrant). For different biomes, similar results were also observed for tropical and subtropical moist broadleaf forests and montane grasslands (both p < 0.001), while temperate broadleaf and mixed forests, temperate grassland and tundra had positive trends (Fig. 5).
Variation in spatial asynchrony of vegetation productivity provides a possible explanation for the opposite EWI trends when scaling up from local to global scales (Fig. 6). Global-scale CVs (CV G ) of the world significantly decreased, while local-scale CVs (CV P ) increased in the past decades, which is consistent with the contrasting trends of CEWI in local and global scales. We found spatial asynchrony of NDVI significantly increased in the past decades, which implies that vegetation tends to fluctuate in a more inconsistent way through time. For different biomes, we found spatial asynchrony increased for tropical and subtropical moist broadleaf forests, temperate grasslands, and montane grasslands, while temperate broadleaf and mixed forests and tundra showed decreased spatial asynchrony.
We found spatial asynchrony of temperature decreased in the past decades and was positively correlated with asynchrony, indicating that increasing spatial correlations of climate may synchronize the vegetation fluctuation. However, results of standardized multiple linear regression analysis showed that changes in the mean temperature was the dominant driving factor that increased spatial asynchrony of vegetation (Table 1).

Discussion
Spatial pattern of terrestrial ecosystem resilience. We identified regions with decreased resilience of terrestrial ecosystems at high spatial resolutions. Low resilience means that the ecosystem has a lower ability to resist environmental perturbations and a lower recovery speed after disturbance 43 , which increases the collapse risk to an alternative regime. Thus, our work had great implications for the conservation and restoration of sensitive regions. We observed that arctic tundra, Europe, northeast America, Australia, and northeast China were vulnerable regions, which was generally in agreement with the results of multiple studies assessing how current vegetation responds to climate change 24,31,44 . For different ecosystem types, tundra, deserts and xeric shrubland, temperate grasslands, and montane grassland have the lowest resilience, which is generally consistent with previous findings that tundra and mountainous regions are the most vulnerable environments to climate change 45 . Shrubland is often intermediary successional communities, which is sensitive to climate change, especially in the tundra biome 46,47 . Tropical forest is the biome with the highest levels of diversity, and previous study have shown that a more diverse biome is often more productive and more resilient to external perturbations 48 .
Resilience indicators in this study builds on the theory of critical slowing down, which has been widely used to assess the risk approaching the tipping point 5,19,20,[49][50][51] . However, these indicators are not manifested in all cases where regime shifts occur. Thus, the caveats of these indicators should be discussed here. Resilience indicators may sometimes cause a wrong signal (e.g., a regime shift without signals (no alarms) or a signal occurs without experiencing a regime shift (false alarms)). No alarms may be caused by the process or observation error, measuring the wrong variable, infrequent or too short time series data, and interference with spatial processes 52 . False alarms can arise when the environmental shocks have some memory in the temporal evolution. In such cases, EWIs would just reflect changes in the patterns of stochasticity, which may have little relevance to resilience. In addition, EWIs do not always indicate bistability (e.g., ecosystem have abrupt shifts but without hysteresis to environmental changes). Due to these limitations, it is more meaningful to adapt these indicators as an agent of ecosystem resilience rather than as forecasting tools in advance of critical transitions.
Effects of climate change and human pressure on ecosystem resilience. Understanding the controlling factors of global terrestrial ecosystem resilience is a key research priority. The strength of ecological resilience is often related to positive feedbacks between vegetation and abiotic environments 53 . For instance, forest transpiration can enhance soil moisture and local microclimatic precipitation, which facilitate vegetative growth and allow vegetation to recover rapidly. Previous studies showed that~33% of Amazon rainfall originates from vegetation transpiration within the same basin 54 , which can buffer against drought in seasons with reduced precipitation. The positive feedbacks can also be evidenced by the conclusion that the biomass resilience of Neotropical secondary forests strongly depends on water availability 55 . In the context of global warming, high temperature would result in increases in tree mortality in the tropics due to the higher evaporation and increasing aridity 56 . Our results indicated that boreal forest regions were more sensitive to precipitation rather than temperature, which is consistent with previous findings that water availability is becoming increasingly important for boreal forest productivity 57 . In undisturbed boreal regions, an old-age forest may have a decreased probability of fire spread, which stabilizes the community state and reduces its sensitivity to climate change 58 . Under future climate change, frequent periods of low water availability may increase boreal forest decline through effects such as peak summer heatwaves, outbreaks of pests, and increases in fire frequency, which decrease the strength of feedbacks and may lead to transitions from boreal forest to open shrubland or grassland 59 .
Besides climate change, human-caused land cover change may also exert a major threat to ecosystem resilience (e.g., biodiversity loss) 60 . For instance, the land use of deforestation can damage the positive feedback between the forest canopy cover and soil moisture, leading to hotter and drier microclimates than those of natural habitats 61 . Future conservation should focus on regions with increased exposure from climate change and human activity, as well as regions with decreased resilience.
Enhanced global scale resilience and spatial asynchrony of terrestrial ecosystems. Spatial scaling of resilience or stability is key to understanding the sensitivity of ecosystems to habitat destruction. Spatial asynchrony is the major factor driving the scaling of stability from local communities to regional ecosystems 62 . In our data, the temporal trends of spatial asynchrony explain the contrasting patterns of stability between local and global scales. Specifically, at larger scales vegetation resilience may increase through time if spatial asynchrony increased even when local resilience is low (Fig. 6). The spatial asynchrony of primary production among regions has increased in recent decades (Fig. 6), and such compensatory dynamics among different areas provide insurance effects and decrease the variability of the terrestrial ecosystem at the global scale. Such spatial insurance effects were suggested to stabilize global crop production in a recent study 63 .
The strength of spatial asynchrony may result from β diversity and environmental heterogeneity, which is related to dispersal, species invasions, consumer-resource interactions and climatic drivers 35,64,65 . Previous studies had suggested that ecosystem resilience might increase with spatial scales due to scale-mediated effects of diversity 66 . β diversity has considered to be an important driver of spatial asynchrony 33,67 , which suggested that different species compositions show asynchronous responses to spatially correlated environments. Moreover, global warming are large-scale disturbances and thus often increased the spatial correlation in the environment (i.e., decrease the asynchrony of environment change; Supplementary Fig. 1), which can decrease the environment heterogeneity, synchronize the vegetation dynamics, and thereby decrease the global vegetation stability. Our results showed that spatial asynchrony of climate drivers exhibited a positive relationship with that of vegetation productivity (Table 1), consistent with recent findings in pests and birds 35,36 . Moreover, we found that climate warming is the dominant driving factor rather than environment heterogeneity, which increases the spatial asynchrony of vegetation variability directly (Table 1). This result indicates that ecosystems in separate regions may have asynchronous responses to climate warming. In short, our analysis suggests that climate changes have contrasting effects on ecological resilience at local versus global scales, by decreasing local-scale resilience but increasing global-scale resilience through enhancing vegetation asynchrony among geographically disjunct areas.
To conclude, in this study, we identified regions and biomes with decreased terrestrial ecosystem resilience at high spatial and temporal resolutions using satellite data. We found that terrestrial ecosystem resilience decreased locally but not globally. Changes in temperature and climate mean state had a higher explanation than precipitation and climate variability to composite resilience indicator, respectively, which can help to promote management practices related to vegetation regeneration and restoration. Increased spatial asynchrony of ecosystem variance in recent decades is probably the key to understanding stabilized terrestrial ecosystems at the global scale. Furthermore, our results showed that climate warming and variability had a positive influence on the spatial asynchrony of vegetation dynamics, which offset the decreased ecological resilience at the local scale. The remaining challenge is to predict how resilience will change in response to future anthropogenic pressures (e.g., frequent global transportation facilitates seed dispersal and species invasion), which may lead to homogenization of species composition, decreasing spatial asynchrony, and lower global terrestrial ecosystem resilience. In this study, although composite early warning metric showed an increasing trend in most regions, we did not intend to provide an exact threshold that provided information on how close or when the biosphere might fall into a planetary scale shift. More research should focus on threshold-related positive feedback loops and possible synergistic effects among multiple anthropogenic pressures and biological responses to make the Earth ecosystems more resilient.

Materials and methods
Framework. Based on the time-detrended NDVI and ecoregion data, we calculated seven simple EWIs for each pixel and within each biome. We then selected four suitable indicators based on the variance inflation factor (VIF) to generate the CEWI. Climate data was also collected to detect the driving factors of CEWI change at the local scale. Subsequently, the potential difference of CEWI trends were compared when scaling up from local to global scales, which is proved under a research framework of scaling stability (i.e., changes of spatial asynchrony of NDVI in the past decades). All data processing and analysis mentioned in this section were conducted in R with rgdal and earlywarnings packages [68][69][70] . The methodological framework of this study was illustrated in Supplementary Fig. 2.
Datasets. As the longest existing NDVI dataset, Global Inventory Monitoring and Modeling System (GIMMS) NDVI 3 g is an ideal data source for calculating EWIs 71,72 . We sourced the dataset from National Aeronautics and Space Administration (NASA) Ames Research Center (https://ecocast.arc.nasa.gov/data/ pub/gimms/3g.v1/) and denoised it using the maximum value composition method 73 . The following rules were further applied to ensure the accurate calculation of EWIs: (1) the pixels that contained over 40% null values were identified as waters or ice sheets and removed; (2) the pixels located in barren sites were also excluded because their NDVI values were seriously affected by soil 74 ; and (3) we filled the null values using linear interpolation when they were less than 40%. After the above processes, we finally obtained a monthly NDVI dataset consisting of 414 images (from July 1981 to December 2015) and 1,822,400 pixels for each image. To better understand ecological resilience dynamics in various biomes, we focused on eight major biomes selected from the World Wide Fund (WWF) terrestrial ecoregions (Supplementary Fig. 3) 75 , namely Tropical and Subtropical Moist Broadleaf Forests (TrMBF), Temperate Broadleaf and Mixed Forests (TeBF), Boreal Forests/Taiga (BoFT), Tropical and Subtropical Grasslands, Savannas, and Shrublands (TrG), Temperate Grasslands, Savannas, and Shrublands (TeG), Montane Grasslands and Shrublands (MoG), Tundra (Tu), and Deserts and Xeric Shrublands (DXS), which totally accounted for about 92% of the vegetated area. Moreover, the data were resampled to the same spatial resolution as that of the NDVI dataset (8 km or 0.083°) based on the majority rule.
Global monthly temperature and precipitation data were collected as potential drivers of EWI change from the Climatic Research Unit (CRU TS V4.03, https:// crudata.uea.ac.uk/cru/data/hrg/cru_ts_4.03/) 76 . We rescaled the climate data through bilinear interpolation, as is classically done and recommended 77 .
Composite early warning indicator. Seven simple EWIs namely, autocorrelation at first lag (ACF1), autoregressive coefficient of a first-order model (AR1), standard deviation (SD), skewness (SK), kurtosis (KURT), return rate (RR), and density ratio (DR), within a sliding window of monthly NDVI were calculated to evaluate the risk of potential critical transition (Eqs. 1 -7) 68 : where z t is the subset of time-detrended ecological state series (i.e., NDVI in this study) within the current sliding window; μ and σ 2 z are the mean and variance of z t , respectively. n is the size of window; z tþ1 is the first-order lag series; ε t is the residuals obtained by the ordinary least squares method; SP(*) is the spectral density function in the power spectrum analysis, and the density ratio is calculated as the spectral density at low frequency (i.e., SPð0:05Þ) to that at high frequency (i.e., SPð0:5Þ). In these resilience indicators, ACF1, AR1, RR, and DR are early warning signals to capture the rising memory, and SD, SK, and KURT are early warning signals to capture the rising variability and flickering 68 .
The choice of detrending method and window size is vital to EWI calculation. The detrending method is applied to make the time series data satisfy the stationary hypothesis. Since NDVI had been shown to vary linearly in recent years 78,79 , we used a linear model to filter out its interannual trend. To eliminate the influence of seasonal variation in NDVI, we set the window size to an integer multiple of 12 (a year, actually 36 months in this study). However, for four indicators derived from the first-order lag model (i.e., ACF1, AR1, RR, and DR), we set the window size to 37 to ensure that the EWI series was not affected by seasonal fluctuations (Supplementary Fig. 4). After the above processing, we could finally obtain seven simple EWI series of length 379 from an NDVI series.
Although previous studies had proposed multiple EWIs in capturing ecological resilience, we adopted the CEWI to make conclusions more conservative because a single metric could only capture changes either in the rising memory or rising variability 68 . However, we did not use all simple EWIs to calculate the CEWI since some of them were highly correlated. For example, ACF1 and AR1 are equal when z t obeys the standard normal distribution (Eqs. 1 and 2). Hence, we determined suitable indicators using a backward selection, which terminated with the condition that the VIFs for all remaining indicators were less than ten (Supplementary  Table 1). Consequently, four indicators (i.e., ACF1, SD, SK, and KURT) were selected to generate the CEWI. Then, the original EWI value is divided by the longrun standard deviation after subtracting the long-run average (z-score transformation), and the CEWI was calculated as follows (Eq. 8): where CEWI t is the CEWI for time step t; EWI 0 t;i is the z-score transformed ith simple EWI for time step t. N is the number of resilience indicators (the value here is four).
Data analyses. Rank correlation coefficients between the EWI and time, such as Kendall's τ, have been statistically confirmed to be suitable metrics to assess CEWI trends 14 . A more important reason should be emphasized for choosing Kendall's τ instead of other traditional metrics (e.g., slope in the linear model) in this study. The sensitivity of NDVI to vegetation coverage will significantly decrease with the growth of coverage 80 , which may result in a smaller CEWI variation range and thus a smaller slope in dense vegetation areas than in sparse vegetation areas. Therefore, rank correlation coefficients represented by Kendall's τ can compensate for this type of deficiency and make EWI trends comparable across regions. A larger positive values of Kendall's τ indicate a loss of resilience and a greater risk of ecosystem collapse 68 . We thus generate "risk evaluation maps" of terrestrial vegetation based on the value of Kendall's τ of CEWI 20 . Spatial patterns of the trends for each simple EWI was also shown in Supplementary Fig. 5.
We employed variation partitioning and multiple linear regression to evaluate the impact of climate factors on the CEWI. In practice, using the sliding window size of 36, four series, namely, the mean temperature, mean precipitation, CV of temperature, and CV of precipitation, were first created. According to different combinations of these variables ( Supplementary Fig. 6), we then identified the major driving factors from two perspectives using the variation partitioning method: (1) whether temperature or precipitation mainly affected the CEWI (Supplementary Fig. 6a) and (2) whether climatic mean state or climatic variability dominated the variation in the CEWI (Supplementary Fig. 6b). In addition, we defined the relative importance on the basis of the above results to make the difference between two components in contribution more intuitive. For example,   Columns from left to right show the temporal variations of global instability (i.e., the instability across the world or within a biome, CV G ), local instability (i.e., the instability for each pixel, CV L ), and spatial asynchrony (φ). TrMBF tropical and subtropical moist broadleaf forests, TeBF temperate broadleaf and mixed forests, BoFT boreal forests/taiga, TrG tropical and subtropical grasslands, savannas, and shrublands, TeG temperate grassland, savannas, and shrublands, MoG montane grasslands and shrublands, Tu tundra, DXS deserts and xeric shrublands. * , ** , and *** indicate p < 0.1, p < 0.05, and p < 0.01, respectively.
the relative importance (RI T ) of temperature could be derived from Eq. 9: where I T and I P are the variations of CEWI contributed by temperature and precipitation (%), respectively ( Supplementary Fig. 6). Furthermore, multiple linear regression was introduced to determine the direction and intensity of the influence of each variable on the CEWI 81,82 : where X M T , X M P , X CV T , and X CV P are the aforementioned four series, which are the mean of temperature, mean of precipitation, CV of temperature, and CV of precipitation, respectively; β M T , β M P , β CV T , and β CV P are the corresponding regression coefficients; β C is the constant term; and ε is the error term. C i is the derivative of the product term with respect to time and represents the change rate of CEWI, which is directly affected by a certain variable. In this study, we also called C M T (C M P , C CV T , CV CV P ) the contribution of the mean temperature (mean precipitation, CV of temperature, and CV of precipitation, respectively) to the CEWI change 81,82 . Spatial patterns of regression coefficients for each climatic variable were shown in Supplementary Fig. 7.
To explain the potential inconsistent trends of the CEWI at different scales, we calculated the spatial asynchrony (φ) of the NDVI. This indicator comes from the research framework of scaling stability and can reflect the correlation of ecological resilience at different scales in a mathematical sound way 33,63 . In the context of current study, the calculation processes of φ can be expressed as Eqs. 12-14.
where CV t;G , CV t;L , and φ t are the global instability (i.e., the instability across the world or within a biome), local instability (i.e., the instability for each pixel), and spatial asynchrony for the time step t, respectively. X t;i;j is the element located in the ith row and jth column of a covariance matrix X t , which is calculated by the time-detrended NDVI series of all pixels within a biome. Precisely, X t;i;j ¼ ðz t;i À μ i Þðz t;j À μ j Þ, where z t is the subset of time-detrended NDVI series within a sliding window; and μ is the mean of z t ; μ G is the non-time-detrended mean of global NDVI; and n is the number of pixels within a biome. In view of the different usages of the covariance matrix, we expect CV t;G to have a similar trend with the CEWI derived from the averaged NDVI across the world or within a biome, and expect CV t;L to have a similar trend as that of CEWI for each pixel.
To detect the main driving factors of spatial asynchrony of vegetation, standardized multiple linear regression was used. The potential driving factors include mean temperature, CV of temperature, spatial asynchrony of temperature, mean precipitation, CV of precipitation, and spatial asynchrony of precipitation.

Code availability
Code for drawing all main figures can be found at https://github.com/fengyuhao1030/ Resilience_Terrestrial_Ecosystem.
Received: 4 September 2020; Accepted: 8 April 2021; Table 1 Standardized regression coefficients of driving factor analysis of spatial asynchrony. In each regression analysis, the dependent variable is the asynchrony derived from NDVI series, and the independent variables include the mean, CV, and asynchrony of temperature (TEM) and precipitation (PRE).