Ecosystems in China have become more sensitive to changes in water demand since 2001

Abstract

Changes in heat and moisture significantly co-alter ecosystem functioning. However, knowledge on dynamics of ecosystem responses to climate change is limited. Here, we quantify long-term ecosystem sensitivity based on weighted ratios of vegetation productivity variability and multiple climate variables from satellite observations, greater values of which indicate more yields per hydrothermal condition change. Our results show ecosystem sensitivity exhibits large spatial variability and increases with the aridity index. A positive temporal trend of ecosystem sensitivity is found in 61.28% of the study area from 2001 to 2021, which is largely attributed to declining vapor pressure deficit and constrained by solar radiation. Moreover, carbon dioxide plays a dual role; which in moderation promotes fertilization effects, whereas in excess may suppress vegetation growth by triggering droughts. Our findings highlight moisture stress between land and atmosphere is one of the key prerequisites for ecosystem stability, offsetting part of the negative effects of heat.

Introduction

Global-scale warming has led to distinct changes in regional climate patterns, including increasing temperature, increasing intra- and inter-annual precipitation variability, dry spells and extreme climate events¹. These unprecedented climate changes are projected to alter vegetation-atmosphere feedback and ecosystem structure (e.g., spatial configuration of species)², thus threatening ecosystem functioning^3,4. Consequently, characterizing terrestrial ecosystem responses to climate change is critical to support a broad suite of Sustainable Development Goals (SDGs), including taking urgent action to combat climate change, sustainably managing forests, halting and reversing land degradation, and halting biodiversity loss^2,5.

Quantifying aggregated biophysical responses of vegetation to climate change is key for understanding how ecosystem functions adapt to global warming^6,7,8. Terrestrial vegetation serves as an important link between the land and atmosphere⁹ through the coupled water-energy-carbon cycle¹⁰. Progress has recently been made in quantifying the vegetation response to climate stresses, but the focus has been on a single meteorological parameter^11,12,13,14. Increased annual rainfall in tropical savannas has been confirmed to change the ecosystem structure by favoring woody species over herbaceous vegetation¹⁵. Studies also observed high ecosystem sensitivity to soil moisture in many semi-arid and arid regions⁵ and reported that the sensitivity to precipitation in drylands was closely related to CO₂^16,17. However, separate consideration of climate metrics can impede our full understanding of ecosystem functions under interrelated factors because of the mutual dependence of vegetation growth on water and energy availability¹⁸, as well as flexible physio-ecological regimes¹⁹. Therefore, there is an urgent need to generate ecosystem responses to compound climate changes to cope with big-picture adaptive strategies^20,21.

Many statistical methods have been used to model the mean ecosystem sensitivity to periodical climatic conditions at spatial or space-for-time substitution patterns²². However, there is no clear strategic and mechanic consensus of temporal patterns of ecosystem sensitivity²³. Recently, an approach based on the autoregression model has identified areas sensitive to climate variability and quantified the relative importance of climate variables in driving local productivity patterns^11,22. The state-of-the-art earth observations enable better monitoring of ecosystem changes²⁴. The enhanced vegetation index (EVI) from Moderate Resolution Imaging Spectroradiometer (MODIS) sensors improves sensitivity to densely vegetated areas by correcting for canopy background noise and atmospheric conditions²⁵. In addition, advances in machine learning applications with increasing availability of earth observations are facilitating the study of ecosystem feedback mechanisms without any premise of explicit functional relationships, e.g., extreme sklearn random forest (RF) models and eXtreme gradient boosting (XGBoost) models. These optimized eco-climate indicators and methodologies allow us to more precisely establish and attribute ecosystem responses to environments, yet they have not been addressed directly from observations.

Terrestrial ecosystems in China have a strong potential for carbon sequestration and are critical for the world’s carbon pool²⁶. However, the dynamics of ecosystems functioning as current carbon sinks under global warming in the future are uncertain, especially in relation to vegetation responses^27,28. The objectives of our study are to provide evidence of changes and drivers of the overall ecosystem sensitivity to land‒atmosphere coupling, hereafter referred to as ecosystem sensitivity, ESI. We first calculated ESI according to Seddon’s approach²² at the pixel level. Specifically, by sliding moving windows of years, we derived the long-term series of ecosystem (productivity) sensitivity to climate changes (here, larger ESI indicates more yields per hydrothermal condition change). ESI contained separate measures of vegetation sensitivity to water (SI_wua)²⁹, solar radiation (represented by cloudiness, SI_cld) and temperature (SI_tmp) by integrating various satellite-based products from 2001 to 2021³⁰. Next, we examined the spatiotemporal changes in ESI based on a bootstrap algorithm across different climate zones and vegetation types. Finally, we explored the relative importance of a set of biotic and abiotic factors in controlling ESI changes in China using the machine learning and the Shapley Additive exPlanations (SHAP) model.

Results

Ecosystem sensitivity to climate change

Normalized values of ESI can be compared in space and time and thereby provide the quantitative magnitude of vegetation changes per unit of climate change. The distributions of annual mean ESI showed significant variations over different ecoclimatic zones (Fig. 1a). Among the five forest biomes, the highest sensitivity values were identified in southern tropical monsoon rainforests, followed by northeastern temperate mixed coniferous forests, cold temperate coniferous forests and subtropical broad-leaved evergreen forests, and the lowest sensitivity was observed in warm temperate deciduous broad-leaved forests (Fig. 1b). Overall, non-forest biomes (mainly related to northern herbaceous vegetation) had relatively lower sensitivity (temperate deserts, grasses, and shrubs) than forests.

Fig. 1: Spatial patterns of multiyear average ecosystem sensitivity (ESI).

a Distributions of mean ESI. b ESI in various biomes (Supplementary Table 1). c ESI changes along the aridity index. TDS temperate desert zone, TGR temperate grassland zone, CCF cold temperate coniferous forest zone, MCF temperate mixed coniferous forest zone, QTP Qinghai–Tibet Plateau alpine vegetation zone, DBF warm temperate deciduous broad-leaved forest zone, EBF subtropical broad-leaved evergreen forest zone, and TRF tropical monsoon rainforest and rainforest zone. Lines in (b) at the top and bottom represent the maximum and minimum values, respectively. The extents of the colored boxes indicate the 25th and 75th percentiles. SI_wua, SI_cld, and SI_tmp refer to vegetation sensitivity to water availability, cloud and temperature, respectively (Supplementary Table 2). Shades indicate the 95% confidence interval calculated in each aridity bin through bootstrapping (N = 3000). The orange, yellow, green, and blue ranges represent the arid (0.05 ≤ AI < 0.2), semi-arid (0.2 ≤ AI < 0.5), dry sub-humid (0.5 ≤ AI < 0.65) and humid regions (AI ≥ 0.65), respectively.

Due to complex hydrothermal combinations in space, vegetation responses to meteorological factors in China are also expected to vary with aridity gradients. In drylands, ecosystem sensitivity to the availability of water, solar radiant energy, and temperature all monotonically increased as the aridity index increased. Nevertheless, in humid regions, we found that vegetation sensitivity to cloudiness and temperature was prone to less fluctuation (Fig. 1c), as well as a decline in ecosystem sensitivity to available water because of the sufficient supply of precipitation and soil moisture (Fig. 1c and Supplementary Fig. 1). Our results suggest that dryness affects ecosystem sensitivity to different climatic elements in a nonlinear manner (Fig. 1c).

Changes in sensitivity

We hypothesized that vegetation may gradually become more adaptable to climate change over time or become more responsive to moisture or heat. Integrating the changes at the pixel level, approximately three-fifths of the study area showed significant positive trends (Kendall test, p < 0.05) in ESI during the past two decades (Fig. 2a), which was further evidenced by robust uptrends when sliding the regression windows from 7 to 13 years and mean-aggregating the ESI in space, despite with different relative value sizes (Fig. 2a). ESI trends displayed obvious heterogeneity in the distribution and direction (Fig. 2b), with the humid and dry subhumid regions having the highest proportions of positive trends and semi-arid and arid regions sharing slightly lower proportions. Considering drylands and humid regions exhibited uniformly wider coverage in sensitivity to temperature and water than to solar radiation (Fig. 2c), we inferred an enhanced dominant effect of the water-energy balance on ecosystem dynamics in the future (Supplementary Fig. 2).

Fig. 2: Changes in ecosystem sensitivity index (ESI).

a Distribution of ESI trends. We changed the windows from 7 years to 13 years to test the robustness of ESI trends to the choice of window length. b ESI trends for different biomes. The trend was calculated from the Kendall coefficient between the ESI and time series. Positive values represent increasing trends, and negative values represent decreasing trends. A total of 61.28% of the areas of interest showed a significant increasing trend (p < 0.05). The red points refer to the mean values. c Proportions of positive trends in vegetation sensitivity (Kendall test, p < 0.05) across drylands and humid regions. Abbreviations of biomes and sensitivity please see Fig. 1.

Attribution of sensitivity trends

We performed the SHAP algorithm in the RF and XGBoost models over undisturbed areas to identify the relative importance of multiple factors contributing to observed trends in ESI. After excluding predictors that contained similar information, the high forecast accuracy (over 84%) and low mean squared error (less than 0.1) shown in Supplementary Table 3 for both models were supposed to provide compelling evidence that ESI trends could be almost entirely explained by an optimal combination of meteorological, hydrological, and anthropogenic predictors³¹. More specifically, rankings of the relative importance of multiple drivers in the two models were almost the same (Fig. 3), further signifying the reliability of the attribution results. The outputs demonstrated that the observed spatial patterns of ESI trends were strongly correlated with changes in (1) vapor pressure deficit (defined as the difference between saturated water vapor pressure and actual water vapor pressure), (2) inter-annual solar radiation and (3) atmospheric CO₂ concentration. Precipitation was the primary water supply that regulated the vegetation response to climate change (Supplementary Fig. 3).

Fig. 3: Relative importance of multiple factors controlling trends in ESI based on SHAP values.

a Output of the random forest (RF) model. b Output of the eXtreme Gradient Boosting (XGBoost) model. Abbreviations of influencing factors please see Methods. Changes of variables refer to the Kendall’s tau coefficient.

To obtain an in-depth understanding of how ecosystems cope with interconnected climate change and human activities, we further calculated ESI trends under various interactions of drivers. We grouped the results from Fig. 3 with respect to the identified main controls and discovered that positive ESI trends were strongest in regions with the most substantial decrease in vapor pressure deficit and the greatest amount of solar radiation. Carbon dioxide changes played a dual role in ESI variability (Fig. 4a), with moderate increasing trends benefiting ESI and overdose having a negative effect. Moreover, downward trends of the points in Fig. 4b illustrated that ESI trends were also restricted by the energy supply system (such as downwelling solar radiation and evapotranspiration) at the time when atmospheric water demand was relatively low. Our results indicated that factors could exhibit synergistic or antagonistic interactions. Overall, we suggest that the major increasing ESI trends under diverse perturbations imply that declining water pressure and appropriately increasing photosynthesis may surpass the negative effect of declining solar energy.

Fig. 4: ESI trends grouped by vapor pressure deficit (VPD) trends, solar radiation (SOLAR) trends and carbon dioxide trends.

a Colors indicate mean values of trends in ESI, integrated by changes in VPD and carbon dioxide; the numbers of grid cells in each group are represented by the size of the “+”. b Similar to (a) but grouped by VPD and SOLAR trends.

Discussion

We applied comprehensive ESI metrics and principal regression analysis to identify sensitive ecological hotspots from a more realistic perspective. This is an improved method to quantify ecosystem functions that simultaneously considers the total absolute effect of multiple climatic elements on primary production after reducing covariance.

We observed that woody vegetation was more sensitive than herbaceous vegetation to climate change (Fig. 1), contrary to the findings of sensitivity of primary production to precipitation across the United States¹⁶. We acknowledge that the literature is messy regarding sensitivity metrics, in which a popular calculation is the coefficient from regression, linear mixed-effects models and other slopes. Here, we are interested in the variance ratios of vegetation productivity variability and each of the climate variables, but it would not make much sense to determine the variance ratios of a location for a climate variable that was not very important in driving productivity patterns; thus, we scaled the ratios by the coefficients. This study practiced another definition of sensitivity and may provide complementarity for the precise assessment of the relative changes in ecosystems to climate variability. Our results also revealed that forests of different shapes and physiologies respond differently to environmental elements. Needle-leaved forests in northern China were most sensitive to variations in water availability, while broad-leaved forests in mild temperate and tropical zones shared the greatest sensitivity to cloudiness (Supplementary Fig. 4), which was consistent with recent satellite observations of forest-cloud cover interactions³². Compared to clear-sky conditions, moderate radiation conditions with a certain amount of increased cloud cover in Supplementary Fig. 5 led to a decrease in total solar radiation (Supplementary Fig. 6). The increase in scattered radiation was most beneficial to the net carbon uptake of forest ecosystems, especially for the forests in southern China during the last two decades (Supplementary Fig. 6)³³. In contrast, angiosperm plants with flatter, broader leaves exhibited high evapotranspiration, which supplied abundant water vapor for cloud formation and maintained large-scale water circulation (Supplementary Fig. 6)³⁴. The low albedo and high roughness of broad-leaved forests promote the partitioning of more solar energy into turbulent heat flows, which heightens the turbulent mixing and convective instability in the boundary layer^32,35. The present research determined that cloudiness in the last decade has promoted the growth of subtropical broadleaf evergreen forests in the Northern Hemisphere^36,37,38. Our work has obtained a more comprehensive picture of the compound climate effects on ecosystem productivity, enabling a better understanding of the underlying drivers and potential mechanisms by encompassing water and energy variables. Of note, our identified possible sensitive hotspots are instructive for ecological protection and carbon emission reduction strategies worldwide.

We discovered a robustly increasing ecosystem sensitivity to water, cloudiness, and temperature changes since 2001 (Fig. 2), which was mainly attributed to water demand (i.e., vapor pressure deficit) and modulated by CO₂ concentration and solar radiation changes (Supplementary Fig. 3). We indicated that climatic factors greatly influenced vegetation sensitivity in three main ways: by determining the water requirements, by regulating photosynthesis and by supplying solar energy. Recognizing these dominant processes greatly helps safeguard sensitive regions under water limitations in the changing environment.

The increasing sensitivity to water availability in China facilitated vegetation greenness (Supplementary Fig. 2), which cooccurred in global vegetation resilience linked to water availability and variability³⁹. Sensitivity to cloudiness increased in most forests and decreased in grasslands (Supplementary Fig. 2). The temperature increased in 44.20% of China and decreased in 31.89% of China (Supplementary Fig. 6). Our findings emphasized the amplified, either degradative or facilitated, compound effects of climate change on greenness over the last two decades. To better anticipate carbon stocks in terrestrial ecosystems, dynamic vegetation modeling should consider changes in nonlinear feedback.

VPD substantially controlled the increasing trend of ESI by offsetting the detrimental effects of the excessive CO₂ concentration trend and solar radiation trend. The declining VPD improved plant stomatal conductance and photosynthetic efficiency, leading to gains in their biomes (Supplementary Fig. 6)⁴⁰. Numerous forests with adaptable stomatal closure strategies were present where there displayed low vapor pressure deficit. The woody vegetation of wetter habitats was more inclined to open their stomata when atmospheric CO₂ levels declined and humidity increased⁴¹, which resulted in extensive exposure and a high sensitivity to environmental changes (Fig. 4a). In pivotal biographical zones, such as China’s Loess Plateau and the Qinling Mountains, the change in vegetation sensitivity was highly susceptible to the SPEI and VPD⁴². The growth rate of sensitivity increased with the slowing rate of VPD, implying that there will be a factorial amplification of the moisture effect on vegetation response mechanisms in the future, in accordance with the recent observed increase in the vegetation water constraint of vegetation growth^7,43. Our results offer empirical evidence of an accelerated ecosystem productivity adaptation regime under global warming and urge specific water allocation and management efforts.

In addition, the ESI trend was greatly influenced by the continued increase in CO₂. The positive ESI trend in Fig. 4a indicates that an appropriate increase in atmospheric carbon dioxide can increase vegetation sensitivity by increasing photosynthesis and increasing leaf area⁴⁴. However, we also found that a Kendall trend of CO₂ greater than or equal to 0.09 may reduce vegetation sensitivity, implying a declining CO₂ fertilization effect, perhaps by reducing stomatal conductance and resulting in warming, eventually altering precipitation regimes, more extreme weather events and so on (Fig. 4a)⁴⁵. Efficient photosynthesis may be saturated if the CO₂ concentration exceeds a certain threshold, and long-term warming can even disrupt the stability of ecosystems. Some recent work has indeed declared that the global CO₂ fertilization impact decreased between 1982 and 2015, and CO₂-induced vegetation browning is projected to occur in future decades^46,47.

In this study, we observed a pronounced increase in vegetation sensitivity to climate change in China (Fig. 2), as well as overall greening since 2001 (Supplementary Fig. 5a). However, there may inevitably be some biases in this study because we incorporated multiple satellite-based and reanalysis outputs. We have endeavored to eliminate errors by calculating all the sensitivity indicators from MODIS sensors and comparing hydrometeorological data from different station observations and remote datasets. Moreover, we tested the increasing sensitivity trend by using different sliding time windows (from 7 to 13 years) to perform the principal regression analysis and produced similar results (Fig. 2a). To categorize sensitivity trends, we used boosting and bagging algorithms together. The high accuracy and extremely low mean squared error indicated the strong predictive power of machine learning models in capturing the dynamics of vegetation sensitivity (Supplementary Table 3). Moreover, ecosystem sensitivity has broader response characteristics in addition to precipitation, temperature and solar radiation, such as ecosystem in the face of pests and diseases⁴⁸, or under extreme climate events⁴⁹. We should also apply system dynamic sciences and causal diagnosis methods to digest the mechanisms of changes in ecosystem sensitivity. Further in-depth work is still urgently needed to quantitatively portray the direction and magnitude of the ecosystem response to each climate factor in the light of the trade-offs between water-energy limitations and ecosystem functionings.

Conclusions

In summary, this study showed a broad increase in ESI over the last two decades in China, which was primarily driven by decreasing water demand (i.e., vapor pressure deficit) and modulated by CO₂ and solar radiation changes. These findings demonstrated that woody plants with greater latent heat display higher sensitivity to cloudiness than herbaceous plants, and needle-leaved forests have the highest ESI due to their amplified response to the available water. Our results were derived via explainable machine learning, which goes beyond purely correlation-based analyses by effectively isolating the influence of diverse climatic and anthropogenic factors. Overall, the detected increase in ESI reflects an enhanced ecosystem productivity response regime to climate variations under global warming. By identifying regions of strong and increasing ESI, our study highlights hotspots where water constraints can have severe impacts on vegetation growth and thus alter carbon-climate feedback.

Methods

Study area

China comprises a rich variety of terrestrial ecosystems from deserts and grasslands to forests, which share allocated water-heat resources, thereby retaining the spatial structure and functioning of ecosystems⁵⁰. Chinese ecosystems are facing increased threats from natural and anthropogenic disturbances. In this study, human land use, including cropland areas, areas with urban and crop/rainfed/natural vegetation mosaic, and unvegetated areas, such as deserts, permanent snow, and ice cover were masked. To reduce the impacts of classification error and land cover change on ecosystem sensitivity, only stable pixels—defined as areas with no change in the dominant land cover class during 2001–2021 were considered. We also excluded areas with flooded tree cover, fresh, saline or brackish water (Supplementary Fig. 7). According to the Land Cover Classification System (LCCS) developed by the United Nations (UN) Food and Agriculture Organization (FAO) (Supplementary Table 1)⁵¹, we categorized the study area into eight dominant eco-climatic zones: temperate desert, temperate grassland, cold temperate coniferous forest, temperate mixed coniferous forest, Qinghai–Tibet Plateau alpine vegetation, warm temperate deciduous broad-leaved forest, subtropical broad-leaved evergreen forest, tropical monsoon rainforest and tropical rainforest.

MODIS data

Considering that remote sensing has trade-offs among temporal-spatial-spectral resolution, we retrieved four monthly time series of key climatic and ecosystem variables that covered the period of 2001–2021 from the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) sensors (Supplementary Table 2). We used the monthly MOD13C2 Version 6.1 product, providing EVI values at a 0.05° resolution (https://lpdaac.usgs.gov/products/mod13c2v061/). The EVI is an “optimized” vegetation index with improved sensitivity over high biomass regions and reduced atmospheric aerosol effects⁵². Water availability, defined as the ratio of actual evapotranspiration to potential evapotranspiration (ET/PET), was calculated using Google Earth Engine from the 8-day composite MOD16A2 Version 6.1 dataset at a 500-m pixel resolution (https://lpdaac.usgs.gov/products/mod16a2v061/)²². The original 8-day products were mean-composited to monthly values, and then from monthly to yearly values. To match water availability data with EVI data, we resampled ET/PET to a resolution of 0.05° × 0.05° using the nearest neighbor method. We employed the clear-sky days from the MOD13C2 version 6.1 product to obtain a cloudiness measure (the proportion of cloudy to clear-sky days) (https://lpdaac.usgs.gov/products/mod13c2v061/)⁵³. The temperature was acquired from MOD13C2⁵⁴. These data were also aggregated to a resolution of 0.05°.

Hydrometeorological data

We extracted gridded meteorological and water balance variables during 2001–2021 from the TerraClimate dataset (https://climatedataguide.ucar.edu/climate-data/terraclimate-global-high-resolution-gridded-temperature-precipitation-and-other-water), which is suitable for climate-impact analyses in ecological and hydrological systems. Explanatory variables include precipitation, potential evapotranspiration (Penman‒Monteith)⁵⁵, shortwave solar radiation, vapor pressure deficit and near surface soil moisture (0~10 cm)^56,57. These data were mean-composited from monthly to yearly values and were aggregated to 0.05° using the nearest neighbor resampling.

Vegetation continuous fields (VCF)

The yearly percent tree cover at 250-m grid cells was acquired from the MOD44B MODIS Vegetation Continuous Fields dataset (https://ladsweb.modaps.eosdis.nasa.gov/archive/allData/6/MOD44B/). VCF products give a continuous, quantitative portrayal of land surface cover with improved spatial detail and hence are widely used in environmental modeling and monitoring applications⁵⁸.

Aridity index (AI)

We used meteorological data from the China Meteorological Data Network (https://data.cma.cn/) to calculate the AI, which was defined as the ratio of precipitation to potential evapotranspiration⁵⁹. Temperature and precipitation data were obtained by spatial interpolation of relevant records at meteorological stations from 2001 to 2021 using Anusplin software version 4.4⁶⁰. Potential evapotranspiration was calculated based on the Penman‒Monteith evapotranspiration equation⁶¹. We also compared the calculated AI with the China Meteorological Data Network from the TerraClimate datasets. According to the consistent patterns of AI values, China was grouped into five aridity classes (Supplementary Fig. 8a): (1) hyper-arid (HD, AI < 0.05), (2) arid (AD, 0.05 ≤ AI < 0.2), (3) semi-arid (SD, 0.2 ≤ AI < 0.5), (4) dry sub-humid (SH, 0.5 ≤ AI < 0.65), and (5) humid regions (HM, AI ≥ 0.65). HD regions were excluded from the analysis because of the relatively low vegetation cover and data nonavailability.

Land cover types

Constant vegetation maps were downloaded from the United Nations (UN) Food and Agriculture Organization (FAO)⁵¹. In this work, we focused on tree cover, mosaic tree cover, shrub cover, mosaic herbaceous cover, shrubland cover, and grassland cover to investigate the vegetation response to climate change. Combined with the distribution of climate zones in China from a 1:1,000,000 vegetation map (https://www.resdc.cn), we further divided China into eight biome zones (Supplementary Fig. 8b): temperate desert zone (TDS); temperate grassland zone (TGR); cold temperate coniferous forest zone (CCF); temperate mixed coniferous forest zone (MCF); Qinghai–Tibet Plateau alpine vegetation zone (QTP); warm temperate deciduous broad-leaved forest zone (DBF); subtropical broad-leaved evergreen forest zone (EBF); and tropical monsoon rainforest zone (TRF).

CO₂ emissions data

We extracted monthly anthropogenic CO₂ emissions from the CarbonTracker modeling system developed by the National Oceanic and Atmospheric Administration (NOAA) to keep track of sources (emissions to the atmosphere) and sinks (removal from the atmosphere) of carbon dioxide worldwide (https://gml.noaa.gov/ccgg/carbontracker/).

Population maps

The total number of people per grid cell was collected from the statistics of the census and administrative units of WorldPop (https://www.worldpop.org/) and the Center for International Earth Science Information Network (CIESIN) in ~2010 and then extrapolated to other years. The improved spatial demographic database was prepared for high-quality applied research. We used yearly population density maps to examine the continuous social pressure on ESI.

Ecosystem sensitivity to climate variability

We constructed a flow chart (Supplementary Fig. 9) to illustrate the data preprocessing and sensitivity analysis. ESI in this study was a proxy of the sum of the response magnitude of vegetation to water availability (the ratio of actual evapotranspiration to potential evapotranspiration, WATER), cloudiness (the proportion of cloudy to clear-sky days, SOLAR), and temperature (TMP). To improve data comparability and accuracy, we performed zero-mean normalization of the long-term MODIS datasets⁶² and filtered out growing season months with mean EVI below 0.1 and mean monthly temperature less than 0 °C. To reduce any impact of collinearity between variables, we first implemented a principal component analysis (PCA) on the data from all years, including EVI, WATER, SOLAR, TMP and 1-month lagged EVI (EVI_t−1). We obtained the corresponding principal components (PCs). Then, we calculated the sensitivity proxy according to Seddon²² by regressing PCs and established spatiotemporal patterns by sliding 7-year blocks at the pixel level (2001–2007, 2002–2008, …, 2015–2021). The variance ratios between the ecosystem and each principal climate component were used as proxies for ecosystem sensitivity by Eq. (1):

$${{{{{{\rm{RI}}}}}}}=\log (({{{{{{\rm{climate}}}}}}\_{{{{{\rm{anom}}}}}}}+1)/({{{{{{\rm{EVI}}}}}}\_{{{{{\rm{anom}}}}}}}+1))$$

(1)

where climate_anom and EVI_anom referred to anomalies of variables by subtracting long-term mean values of climatic factors and EVI, respectively. Ecosystem sensitivity index is defined as the log10-transformed ratio of climate variables and EVI (e.g., climate_anom, EVI_anom). Next, we quantified the relative weights of climatic factors from the significant coefficients (p < 0.05) of multiple regression of PCs (principal component regression, PCR) using Eq. (2):

$${{{{{{\rm{EVI}}}}}}}=\alpha \cdot {{{{{{\rm{PC}}}}}}}1+\beta \cdot {{{{{{\rm{PC}}}}}}}2+\gamma \cdot {{{{{{\rm{PC}}}}}}}3$$

(2)

where the slopes of different components ${{{{{\rm{\alpha }}}}}}$, ${{{{{\rm{\beta }}}}}}$, and ${{{{{\rm{\gamma }}}}}}$ are the corresponding climate weights of the three principal components (${{\mbox{PC}}}1$, ${{\mbox{PC}}}2$, and ${{\mbox{PC}}}3$). ESI was further modified by scaling the ratios by climate weights, as evidenced by correctly reproducing the response of the carbon cycle to essential local climate variability using Eq. (3):

$${{{{{{\rm{ESI}}}}}}}=\alpha \cdot {{{{{{\rm{RI}}}}}}\_{{{{{\rm{wua}}}}}}}+\beta \cdot {{{{{{\rm{RI}}}}}}\_{{{{{\rm{cld}}}}}}}+\gamma \cdot {{{{{{\rm{RI}}}}}}\_{{{{{\rm{tmp}}}}}}}$$

(3)

where ${{\mbox{RI\_wua}}}$, ${{\mbox{RI\_cld}}}$ and ${{\mbox{RI\_tmp}}}$ are the original sensitivity index of ecosystem response to water availability, solar radiation and temperature. ESI is the sum of the weighted vegetation sensitivity to water, temperature and solar radiation (SI_tmp, SI_wua, SI_cld). To verify the robustness, we also compared ESI patterns against different time windows (W = 7, 9, 11, 13 years). All the calculation was finished in R and MATLAB.

Changes in sensitivity

We assessed the trends of annual ESI for the two-decade study period using Kendall τ ranging from −1 to 1, which is a ranked correlation coefficient with nonparametric hypothesis testing⁶³. A positive Kendall τ value indicates an enhanced response of the ecosystem to environmental variables, and vice versa. Kendall τ is a comparable and robust index owing to temporal dependence and represents the changing rate of ESI.

Drivers of changes in sensitivity

To estimate the relative importance of various climatic and anthropogenic variables (predictors) potentially driving long-term changes in ESI (targets), we implemented both the bagging and boosting models in machine learning⁵, which can isolate marginal contributions of each predictor on the target variable⁶⁴. The technical scheme in Supplementary Fig. 10 shows the execution process of the two algorithms. The RF model reduces variance by constructing multiple mutually independent learners based on a bootstrap aggregating strategy⁶⁵. The XGBoost model also reduces bias by constructing strong, interlinked evaluators based on a boosting tree strategy⁶⁶. For both models, we used RandomizedSearch to evaluate the results for different combinations of hyperparameters and thus determined the optimal parameters (for RF: number of estimators: 2500, maximum features: 30%, random state: 42; for XGBoost: training size: 30%, learning rate: 0.015, number of estimators: 400, maximum depth of tree: 8, minimum value of loss function: 0, sum of weights of the smallest leaf node samples: 10). In particular, we applied the out-of-bag error (OOB error) to measure the performance of the random forest and directly modeled all the predictors and target variables⁶⁷. For XGBoost, we split the data into testing and training parts. To avoid overfitting, we additionally performed fivefold cross-validation.

Moreover, we quantified the marginal contributions of predictors to ESI trends in each model by employing the SHAP method originating from the coalitional game theory⁶⁸, which has been widely used to reflect the influence of features in each sample in explainable machine learning. All the processes were conducted using the packages of “xgboost”, “RandomForestRegressor”, “shap”, “sklearn.model_selection” and “sklearn.metrics” in Python software (version: 3.8)^69,70.

Climate changes are likely to have profound impacts on the structure and functioning of Earth’s ecosystems. Here we selected climatic factors related to water, heat, and energy to identify the potential drivers of ecosystem sensitivity trends, including incoming shortwave solar radiation (SOLAR), temperature (TMP), precipitation (PRE)⁷¹, vapor pressure deficit (VPD) and soil moisture (SM)⁴⁸, potential evapotranspiration (PET)⁷², and aridity index (AI)⁷³. The elevated CO₂ has been demonstrated to promote vegetation growth and thus affects ecosystem structure and functioning⁷⁴. We also included the vegetation cover fraction (VCF) to provide the information of vegetation structure. Additionally, we used population (POP) as a proxy of human activities.