Changes and influencing factors of ecosystem resilience in China

The multifunctionality and sustainability of ecosystems are strongly dependent on their ability to withstand and recover from disturbances—that is, ecosystem resilience (ER). However, the dynamics and attributes of ER remain largely unknown, especially in China, where climatic and anthropogenic pressures are high. In this study, we evaluated spatiotemporal patterns of ER in China from 2001 to 2020 using solar-induced chlorophyll fluorescence. We estimated the relative independent importance of climate change, CO2, and anthropogenic factors on changes in ER signals. The results showed that more than half of the ecosystems in the study area have experienced ER gain followed by ER loss during the past two decades. Before breakpoints (BPs), climate change explained 58.29% of the ER change associated with increasing precipitation. After BPs, 65.10% of the ER change was most affected by CO2, and drought from rising temperature further deteriorated ER loss. We highlight that relationships between changes in ER and climate are spatially heterogeneous and suggest increased negative radiative effects of CO2, associated with global warming, on ecosystem stability due to the saturated canopy photosynthesis. These findings have crucial implications for future climate change mitigation, carbon peak, and carbon neutrality targets.


Introduction
The stability of terrestrial ecosystems plays a fundamental role in maintaining the biogeochemical and hydrological cycles as well as the long-term provision of essential goods and services for human life. However, ecosystems are increasingly disturbed by climate change and numerous human activities, with uncertain consequences (Boulton et al 2022, Wu et al 2022). For example, frequent extreme heat events in areas such as the Amazon rainforest have increased the risk of widespread forest mortality (Ciemer et al 2019, Hubau et al 2020. Warming may also promote vegetation growth in the alpine regions of the Northern Hemisphere (Wang et al 2020a). The ecosystem responses to multivariate interference are dynamic, nonlinear, and largely unknown.
International environmental policies must address this gap to achieve carbon neutrality because ecosystem changes control global biodiversity and carbon pools (Hubau et al 2020, Gatti et al 2021. Ecosystem resilience (ER) is a strong indicator of ecosystem stability and refers to the ability of a system to persist and recover in the face of perturbations (Dakos et al 2015, Scheffer et al 2015. The generic phenomenon of critical slowing down (CSD) has attracted considerable attentions. As a system gradually loses resilience and becomes increasingly unstable, it recovers more slowly from short-term perturbations (Scheffer et al 2009, Liu et al 2019. This slowness causes the fluctuation of the system state to be relatively close, with an increasing correlation to the previous state at a one-time step (the lag-1 auto-correlation [AR(1)] in the time series) (Dakos and Bascompte 2014, van der Bolt et al 2018, Liu et al 2019 and wider variability of natural fluctuations estimated in terms of the variance. Regional studies have measured ER, but present limitations, including that (i) most are based on mean climatic state and cannot reveal the ER variation; (ii) practical application of CSD has been confined for rainforests and other vegetation types, and aridity gradients outside humid tropical areas have rarely been addressed; and (iii) corresponding attributes of ER dynamics are lacking, as well as climatic, anthropogenic, and CO 2 elements, which can influence productivity and biomes (Smith et al 2022). Thus, the regional assessment of ER is a requisite task.
Solar-induced chlorophyll fluorescence (SIF) is an optical signal that is emitted from vegetation when sunlight activates photosynthesis, which is linked mechanistically to gross primary productivity (GPP) and has shown a near-linear relationship with GPP at the ecosystem scale when measured by satellite (Li and Xiao 2019, Wei et al 2022). It is superior to conventional vegetation indices for monitoring complex photosynthesis across diverse terrestrial ecosystems, as it can separate vegetation and soil backgrounds and respond more rapidly to environmental stress (Moya et al 2004, Jiao et al 2019. Sensors such as the Orbiting Carbon Observatory-2 (OCO-2) have been used to retrieve SIF (Li and Xiao 2019, Ryu et al 2019). However, there have been insufficient studies using SIF to investigate ER to date.
Terrestrial ecosystems in China are among the most important carbon sinks in the Northern Hemisphere (Wang et al 2020b). China is under high resource pressure, with 6% of the total freshwater supporting nearly 20% of the global population (Hu et al 2022). Since the 21st century, severe climatic and increased anthropogenic disturbances have further affected ecosystem structure and function in complex ways and differing magnitudes. Here, we attempt to evaluate the dynamics of ER in China and effects of climate change, carbon dioxide concentration and human activities on ER changes, focusing on three key issues: how ER changes over time and space; whether there are differences in ER along aridity and biome gradients; and how natural-social factors influence ER. By answering these questions, we aim to identify underlying processes that encourage ecosystem growth rather than using visible greenness or vegetation cover as an indicator, for more accurate prediction of ecological risk areas.

Study area
To reduce the effects of classification errors, only grids with the same dominant land-cover classes in China from 2001 to 2020 were used in our land cover-specific analyses (Forzieri et al 2017). We removed any pixels that represented areas subjected to human land use, such as cropland, urban and crop/natural vegetation mosaic cover, and unvegetated areas, including deserts, permanent snow, and ice (figure S1). We mapped 161 234 grid cells representing nine land cover classes, including broadleaved and coniferous tree cover, shrubland, grassland, and mosaic forests, as our study objects for calculating ER (table S1). Afterwards, we averaged the regional values to represent the dominant land cover attributes.

Land cover
We used the Land Cover Classification System developed by the United Nations Food and Agriculture Organization (table S2) (https://cds.climate.copernicus.eu/cdsapp#!/dataset/ satellite-land-cover?tab=overview). Moreover, for more accurate identification of unchanged vegetation, we combined the 30 m annual land cover data set, and Landsat-derived annual China land cover data set (CLCD) from 1990 to 2019 (Yang and Huang 2021). It is the latest fine-resolution land cover data set produced for China using substantial Landsat images via a post-processing method incorporating spatial-temporal filtering and logical reasoning to further improve the spatial-temporal consistency of CLCD (https://zenodo.org/record/ 4417810). The overall accuracy of CLCD outperforms that of MCD12Q1, ESACCI_LC, FROM_GLC, and GlobeLand30. Afterwards constructing the datasets, we averaged the regional values of ER to represent the dominant land cover attributes. Based on the annual vegetation types and their geographical distributions from multiple datasets (www.resdc.cn), we divided China into eight biome zones (figure S2).

ER indicators
To remove the linear trend and seasonality, we initially adopted seasonal trend decomposition using Loess (STL) to split the vegetation datasets into the seasonal, trend, and residual components for each grid cell (Boulton et al 2022). The different span (in lags) of the Loess window (t. window) for trend extraction was tested (figure S5) and chosen as 19 months in the analysis. The residual component was then used to measure AR(1) and variance as robust indicators of the CSD (Dakos et al 2008). We set a sliding window length of 60 months to create a time series for AR(1), with 181 records for each pixel in the main text. We then tested the results at different window lengths to verify their robustness to changing windows (figures S5-S7), as demonstrated in previous studies (Boulton et al 2022). To ensure reliability, the coefficient of variation (hereinto, variance) was calculated as a secondary indicator (Dakos et al 2012) where α 1 is the autoregressive coefficient, and ∈ t is Gaussian white noise process. V is the variance of variable Q t and µ is the mean of variable Q t . An increase in AR(1) indicates that the system state has become unstable. The loss of resilience manifests as an increase in the variance over time. All analyses were performed using MATLAB 2021a and R 4.0.5.

Breakpoint (BP) identification and trends in ER
For the AR(1) monthly time series and variance per pixel, we detected and analyzed breaks for additive seasonal trends (BFAST) without selecting a reference period, setting a threshold, or defining a change trajectory (Verbesselt et al 2010). This has also been widely used as an alarm system to provide information on when and where changes have occurred using the R packages' bfast' , 'tseries' , and 'rlang' . The shortest segment was set to 0.25 of the total data length to limit the minimum number of observations in each segment (Li et al 2022). The most significant BPs were tested using the Chow test (P < 0.05) (Verbesselt et al 2012).
We calculated the Kendall rank correlation coefficient (Kendall τ ) to describe the temporal trends of AR(1) and variance (Liang et al 2021). A positive Kendall τ value implies that time series are always increasing, negative implies decreasing, and 0 indicates there is no overall trend (Pirnia et al 2018). We tested statistical significance based on the original data and a set of 1000 surrogate simulated time series (Boers 2018).

Principal component regression (PCR)
In addition to precipitation (ppt), temperature (tmp) and solar radiation (srad) as climatic factors that affect vegetation dynamics and ecosystem function (Seddon et  Therefore, this study investigated the relative contributions of climatic factors, carbon dioxide concentration and anthropogenic factors to annual changes in ER by the relationships between AR(1) and multiple factors. Considering there may exist multicollinearity among these factors, we used the PCR method to separate contributions (Li et al 2022). All explanatory variables were tested for no multicollinearity by variance inflation factor values (VIF < 5). PCR is an effective way to isolate the relative independent importance of natural and anthropogenic factors on AR(1) for each pixel (Seddon et al 2016). Akaike information criterion was applied to further avoid overfitting (Dannenberg et al 2019). To eliminate seasonal cyclical pseudo-correlations, we abstracted trend components of these time series using STL (Song et al 2022). To eliminate multicollinearity, we converted the trend series for the six factors into independent principal components: where PC1, PC2, PC3, PC4, PC5, and PC6 represent the principal components of climatic factors and CO 2 , respectively. Before entering these into the model, we transformed all trend series into z-score anomalies: ppt ′ , tmp ′ srad ′ , vpd ′ , sm ′ , and CO 2 ′ (Zhang et al 2021b).
To isolate the relative importance of climate variables and CO 2 on ER, we performed multiple linear regressions between the principles and AR(1) (equation (4)): where RES denotes the residual. For the principal components with significant relationship with AR(1) (P < 0.1), we multiplied the loading scores for each variable (α i , i = 1, 2, 3, 4, 5, 6) using PCR coefficients (k ij , j = 1, 2, 3, 4, 5, 6) and summed the scores. Therefore, ∑ 6 i = 1 α i k ij was used to estimate the relative importance of each variable in driving changes in AR(1). We rescaled the environmental weights to values between 0 and 1 to compare their relative contributions (Lü et al 2015). The linear trend of the residuals was defined as the anthropogenic factor where β 0 is the intercept and the significant coefficient β 1 (P < 0.05) reflects the anthropogenic effect (Li et al 2022).

Spatial pattern of ER in China
The loss of ER can be detected by an increase in AR(1) as a less stable system with a lower return rate from perturbations. Thus, AR(1) at different locations provides the relative size of local ER. The distribution of ER in China was characterized as being high (low AR (1)) in the south and low in the north ( figure 1(a)). For variance, the estimate is influenced by noise at the pixel level and exists theoretical uncertainty from stochastic dynamics (Smith et al 2022).
Thus, the pattern of variance showed different spatial heterogeneity in ER with AR(1) (figure S3), which also has been found in previous research (Boulton et al 2022). The annual mean AR(1) decreased with increasing AI ( figure 1(b)). When integrating AR(1) according to land cover, the tropical monsoon rainforest, tropical rainforest, and temperate mixed coniferous forest zones displayed a lower AR(1) than other biome types ( figure 1(c)).

Reversal of ER changes from gain to loss
The ER in China experienced different amplitudes and directions of changes between 2001 and 2020 ( figure S4). The shifts for changes in ER represented by AR(1) and variance were both in 2013-2014, accounting for approximately 29.04% and 26.97% of the regions, respectively (table S3). Earlier BPs of AR(1) time series occurred in northeast Inner Mongolia, bordering Greater Khingan Prefecture and southeast China, suggesting locations with higher sensitivity to external stressors (figures 2(a) and (d)). Approximately 66.54% of the study area experienced AR(1) trend reversals, of which 50.21% increased before BPs and decreased after BPs, and 16.33% showed the opposite trend (figures 2(b) and (c)). Similarly, we detected BPs of the variance and observed an increasing trend in 58.86% of the grid cells before BPs and a decreasing trend in 66.92% after BPs (figures 2(e) and (f)). The trends of resilience signals from sliding t. windows of STL and window lengths of autocorrelation all showed great robustness (figures S5-S7). In addition, AR(1) calculated from different vegetation datasets had different values but the same reversal trends and shifted in 2010-2015  (figure S8). The consistent shifts of ER from AR(1) and variance suggest the potential of variance as a supplementary for the resilience change, as long-term trends can eliminate the noise effects on a certain location. Thus, we focused on AR(1) changes in the latter ER analysis.
We further explored the spatiotemporal evolution of ER within each biome across the aridity gradient. The ER conversion from gain to loss had the largest proportion (52.23%) in HM. Most of the warm temperate deciduous broad-leaved (78.42%) and subtropical broad-leaved evergreen (53.21%) forests experienced a reversal trend (figure 3). These dissimilar patterns along the AI and in different vegetation regions suggest that atmosphere-vegetation interactions may be one of the causes for ER contrasting changes.

Attribution of changes towards ER
The relative contributions of multiple factors to ER changes were quantified for each grid cell using an empirical approach (that is, the weights related to different climatic variables derived from principal component regression, see section 2.3.2, PCR) (figure 4). Before BPs, climate change explained 58.29% of AR(1) trends, followed by CO 2 (39.79%) and human activities (0.23%); after BPs, CO 2 had the strongest effect on AR(1) changes (65.10%), followed by climate change (33.86%) and human activities (0.07%). Before BPs, positive anthropogenic impacts were mainly observed in the southeast Three-North region, including southeast Inner Mongolia and northwest Heilongjiang Provinces (35.02%, AR(1) decreased at a rate of 0.17 × 10 −3 yr −1 , suggesting an increase in ER). After BPs, climatic factors associated with increasing water demand from substantial increases in PET had a negative effect of −1.01 × 10 −3 yr −1 across central China (57.82%). With increasing temperature and largescale drought, rising CO 2 displayed the effect on more regional resilience loss, where the CO 2 fertilization effect on plant growth may have already reached saturation or been offset (56.34%, AR(1) increased at a rate of 11.51 × 10 −3 yr −1 , suggesting a decrease in ER). Anthropogenic impacts weakened with the maximum average negative effect in 17.36% of humid areas (AR(1) increased at a rate of 0.04 × 10 −3 yr −1 ).

Discussion
Multiyear averaged ER in China displayed spatial heterogeneity in the stability of ecosystem recovery from local disturbance. The spatial patterns suggested that the tropical rainforests had greater ER than cold temperate coniferous forests (figure 1). Compared with the global resilience tested by enhanced vegetation index and vegetation optical depth, ordered ER by land cover types was consistent, from high to low, rainforest, broadleaf, temperate mixed coniferous, cold temperate coniferous, and temperate grassland Xia 2019, Smith et al 2022). However, different from the ER change shift towards a global decline since the early 2000s, the overall BPs of ER signals time series in China were detected around 2014 (Smith et al 2022). The aggregated time series of ER signals in China exhibited spatial heterogeneity and were more similar with global AR(1) than variance, which also shifted around 2014, especially for evergreen vegetation (Smith et al 2022).
ER changes have been significantly affected by interactions of multiple climatic factors. Sensitivity analysis revealed that before BPs, ER was most sensitive to changes in soil moisture and vapor pressure deficit ( figure S9(a)). The accelerated hydrological cycle with increasing precipitation positively affected ER response to drought in northern arid areas of China. Notably, precipitation does not always have a positive effect on all terrestrial ecosystems. For example, excessive precipitation negatively affected tropical monsoon rainforests due to haze and clouds (Verbesselt et al 2016). After BPs, ER was most sensitive to changes in temperature ( figure S9(b)), when weaker southeast trade winds led to a warm, dry winter in northern China, related to the PDO (Gao et al 2020, Chen and Li 2022) (figure S10). The temperature across the entire region increased rapidly, contributing to a corresponding upward trend in PET (Oren et al 1999). Regional case studies have demonstrated that PET contributed more to long-term aridity trends than the precipitation . Furthermore, temperature affected ecosystems with large uncertainty. Recent intensifying global warming has benefited forests in north-eastern China as the increased heat enhanced physiological activities, such as vegetation photosynthesis and water recycling. We found that there were areas with significantly decreasing trends of soil moisture even though the overall national trends were increasing (figure S14(d)), for example, on the Loess Plateau, which had lower resilience than other semi-arid areas (Zhu et al 2019) (figure 1(a)). The rising temperature increased PET (figure S14(a)), which further decreased the soil moisture when the precipitation supplement was insufficient (Xu et al 2022). This decrease in soil moisture reduced evaporative cooling, which increased soil temperature and atmospheric moisture demand (Miralles et al 2014). Moreover, soil supplied necessary nutrients to vegetation for recovery from disturbances. These positive temperature-associated couplings combined with multiple water and biological cycles resulted in high atmospheric moisture demand, which made an active impact to the earlier BPs and widespread reduction in ER after BPs (figure 2). Together, the present spatiotemporal results demonstrated that the relationships between ER and climate conditions varied among different regions and were also influenced by ocean-atmospheric circulations (figure S10) (Kim et al 2019, Wang et al 2022). Our comprehensive analysis provides important findings that help understand the response of ecosystem functions to terrestrial-atmospheric variables.
Besides limited by energy and water, the spatially heterogeneous BPs were also impacted by human activities. Human activities, particularly the largescale environmental protection programs in China (figure S11), have improved ecosystem functions and enhanced ecosystem services (Lu et al 2018, Rashid et al 2020), such as soil erosion control and water retention (Chen et al 2021b). Before BPs, ecological protection such as the 4th phase of the Three-North Shelter Forest Programme has benefited ER gains, as well as the increase in leaf area index (in 64.24% of China) and SIF (71.51%) (figures S12 and S13). However, the effects gradually turned negative predominantly in the Loess Plateau region. Increasing precipitation and temperature in China in recent years have promoted temporary vegetation growth, but this was at the cost of soil moisture and was unsustainable (Feng et al 2016).
The effect of rising atmospheric CO 2 concentrations on the ER change is complicated, as is the coordination of fertilization and radiative effects of CO 2 (Wang et al 2020c, Winkler et al 2021. CO 2 fertilization involves two pathways: (i) vegetation lowers stomatal conductance and stomatal density when the atmosphere is filled with accumulated CO 2 (Fatichi et al 2016). Consequently, decreasing water loss by decreasing transpiration can improve water use efficiency, which further improves ER (Adams et al 2020). (ii) Increased CO 2 can stimulate carbon assimilation, thereby increasing the plant biomass (Luo et al 2004). However, the radiative effect of CO 2 has both benefits and disadvantages by influencing climatic change: (i) warming can prolong the growing season of temperature-limited plants and promote ER (Walker et al 2021). (ii) Long-term high-intensity aridity can disrupt the ecosystem (Seleiman et al 2021, Seidl et al 2017). This may depend on plant species characteristics, hydraulic diversity, and local environment (González de Andrés et al 2018, Zhang et al 2022b). The study found that a two-sided effect of CO 2 on ER has occurred in China (figures 4(c) and (d)) and suggested more focus on the radiative effect as canopy photosynthesis may saturate at high CO 2 concentrations under the ongoing global warming.

Conclusion
We inferred that the ecosystems in China have experienced resilience gain followed by resilience loss since 2001 from CSD indicators, although vegetation is continued greening. BPs in ER signals change were predominantly concentrated between 2013 and 2014. Before BPs, climate change greatly impacted ER, whereas after BPs, CO 2 had a more considerable effect on ER. We suggest that more attention be given to increasing negative radiative effects of CO 2 on ecosystem stability. By simultaneously considering climatic, anthropogenic and carbon fertilization effect, a fuller explanation of regional changes in ecosystem functions can be obtained. Such work can provide early warnings for hotspots at risk of ecosystem degradation and help design effective vegetation-based mitigation strategies for the goals of 'carbon peak' and 'carbon neutrality' .

Data availability statements
All data that support the findings of this study are included within the article (and any supplementary files). All data used in this study is publicly