ENSO amplifies global vegetation resilience variability in a changing climate

Domain: physics
Published: 2025-12-04
Editorial outcome: accept

Abstract
A thorough understanding of vegetation resilience to climate variability is critical for sustaining ecosystem functions and terrestrial carbon sinks. Although the El Niño-Southern Oscillation (ENSO) is a key driver of global extreme weather events and vegetation dynamics, its impacts on vegetation resilience remain unclear. Here we estimate global present-day (1981–2018) and future (2015–2100) vegetation resilience using a lag-1 autocorrelation analysis of global leaf area index (LAI) time series and investigate its teleconnection to ENSO. Our findings reveal that ENSO significantly affects vegetation resilience across 53% of the global vegetated area. Within these regions, 15% are linked primarily to large-scale atmospheric synchrony with ENSO, 51% are mainly shaped by ENSO-driven local climate anomalies, and the remaining 34% are influenced by both processes. Future projections suggest that the area impacted via ENSO-driven climate anomalies may expand by 7-10%, with Eastern Siberia and northern North America newly affected. Our study provides a coherent global assessment of vegetation resilience sensitivity to ENSO, identifies teleconnected hotspots and potential influential pathways, and informs targeted restoration and climate-adaptive ecosystem governance under climate change.

Introduction
Vegetation resilience refers to the capacity of vegetation to maintain or recover its structure and function following disturbances1. However, the intensification of global climate change poses a significant threat to this resilience. Previous studies have predominantly focused on the relationship between resilience and temperature and precipitation patterns2–4, often neglecting the remote impacts of global climate oscillations such as the El Niño-Southern Oscillation (ENSO). As the largest natural perturbation to the global climate on an interannual timescale5, ENSO significantly affects rainfall and temperature dynamics, thereby reshaping the distribution of meteorological hazards such as droughts and floods6,7, as well as secondary hazards such as wildfires8. These disruptions affect core ecological functions, including species composition, nutrient cycling, and energy flow, ultimately undermining vegetation resilience9. Moreover, with global warming, ENSO variability is expected to increase10, potentially amplifying its influence on plant growth and ecosystem stability worldwide11,12. Understanding the impact of ENSO on vegetation resilience is crucial for guiding nature-based solutions—such as reforestation, biodiversity protection, and climate-smart land management—that contribute to the UN Sustainable Development Goals13.

ENSO is a coupled ocean–atmosphere phenomenon involving periodic fluctuations in sea temperature and atmospheric pressure in the tropical Pacific14. It typically develops during the boreal summer and decays in the following spring15,16, with peak intensity usually occurring between December and February17. Through large-scale atmospheric circulation anomalies, ENSO alters regional patterns of temperature, precipitation, and solar radiation18,19, thereby disrupting vegetation biophysical processes20. Specifically, the alternating warm (El Niño) and cold (La Niña) phases of ENSO exert complex and spatially heterogeneous vegetation responses12. El Niño typically raises temperature and solar radiation in tropical and subtropical regions, yet its impact on precipitation varies, causing wetter conditions in some areas and drought in others, potentially exacerbating or alleviating vegetation water stress21–23. For instance, El Niño-driven shifts in the Intertropical Convergence Zone modify regional precipitation patterns, decreasing precipitation in West and northern Central Africa, while increasing it in East and southern Central Africa24, thus creating a contrasting north-south gradient in vegetation dynamics. Conversely, La Niña generally exhibits the opposite climatic signals, particularly over tropical land areas where temperature and precipitation patterns respond in antiphase, thereby partially offsetting the effects of El Niño25,26. However, this compensatory mechanism is modulated by variation in the location of equatorial warming and by nonlinear ocean–atmosphere interactions, leading to spatial asymmetries in the timing and intensity of vegetation response to ENSO phases27. Consequently, some regions, such as the Amazon rainforest, are highly sensitive to El Niño28,29, while others, like Australian shrublands, experience alternating impacts from both phases30.

More critically, high-intensity ENSO events can lead to irreversible declines in vegetation resilience. A notable example is the 2015–2016 extreme El Niño event, during which humid forests in Africa and the Americas experienced sustained reductions in aboveground carbon storage. Even after climate conditions normalized in 2017, recovery remained slower than in dryland ecosystems31,32. This contrast highlights that the effects of ENSO on vegetation are not solely determined by event intensity but are also influenced by ecosystem type, regional climatic conditions, and water availability33. Nevertheless, the magnitude and spatial distribution of ENSO’s effects on vegetation resilience remain poorly understood. This knowledge gap is particularly concerning given projections that ENSO events will undergo significant shifts in frequency, intensity, and impact patterns under a warming climate34,35, potentially exhibiting unprecedented disturbance characteristics that further complicate global vegetation resilience dynamics. Thus, a systematic assessment of the multi-scale impacts of diverse ENSO phases is urgently needed, with a particular focus on temporal and spatial variability across ecosystems and climate regions, as well as sustainability and evolution in future scenarios.

The critical slowing down theory (CSD) provides a framework for detecting ecosystem instability by suggesting that declining recovery rates precede state transitions, which can be identified by statistical signals like increased variance and higher first-order autocorrelation36–38. In this article, we applied this framework to quantify global-scale changes in vegetation resilience using remotely sensed Leaf Area Index (LAI) data since 1981. We hypothesized that ENSO influences on global vegetation resilience operate via two complementary perspectives: the atmospheric synchronization teleconnection, which captures the large-scale statistical association between ENSO events and resilience anomalies, independently of local climate variables; while the climate-mediated teleconnection, which reflects indirect influence transmitted through ENSO-driven anomalies in temperature, precipitation, and radiation (see Methods and Supplementary Fig. 1). Together, these pathways reflect ENSO’s dual role as a synchronizer of global variability and a modulator of local climate conditions. Additionally, we used simulations from the CMIP6 model to analyze how ENSO’s impacts on vegetation resilience may evolve in the future. We present a coherent global-scale evaluation of how vegetation resilience responds to ENSO, pinpointing teleconnected hotspots and revealing key pathways of influence. These findings improve our understanding of the teleconnection between vegetation growth and extreme climate, supporting nature-based solutions and informing sustainable development efforts.

Results
Historical characteristics of vegetation resilience dynamics
We quantified vegetation resilience at the pixel level using lag-1 autocorrelation (AC) derived from 5-year moving windows of GLASS LAI data (1981–2018), and assessed its long-term trend (see Methods). This global assessment indicates that 41.56% of vegetated areas have experienced a statistically significant decline in resilience over the past 40 years, with considerable spatial heterogeneity (Fig. 1a). The most significant declines in resilience are observed in northern North America, the Amazon, northern Europe, eastern Australia, and northern and Southeast Asia, which exhibit the greatest increase in AC (Fig. 1a). Compared with the trend in variance, another CSD indicator, 78.4% of areas exhibit consistent resilience trends in both AC and variance estimates (see Methods, Fig. 1b, c).

Spatial patterns of historical changes in global vegetation resilience contrast with the dynamics of global vegetation greening. Regions such as the Amazon, south-central Africa, Southeast Asia, and western Oceania have experienced a widespread increase in greenness—measured by rising LAI—over the last four decades (see Methods, Fig. 1d), yet these areas have concurrently witnessed significant declines in resilience (Fig. 1a). In addition, over 40% of greening vegetation in global open shrublands, grasslands, evergreen broadleaf forests, and deciduous broadleaf forests exhibited a decline in resilience (Fig. 1e and Supplementary Fig. 2). These findings highlight the simultaneous occurrence of accelerated vegetation greening and declining resilience worldwide.

Controlling factors of historical changes in global vegetation resilience
We used a random forest (RF) regression framework at the global scale to attribute temporal changes in vegetation resilience (AC) to ENSO variability and other environmental factors. This model focused on in-situ climate variability and ENSO as the key explanatory variables, incorporating other static in-situ factors such as topography and climate seasonality as constraints to reduce background interference (see Methods). The analysis reveals that climate seasonality exerts the strongest influence on shaping interannual variation, followed by climatic variability and topography (Fig. 2a; model R2 = 0.79). Among the controlling variables, temperature and radiation seasonality exert the strongest climatic control on resilience, and ecosystems experiencing low intra-annual variability are especially prone to resilience loss. Elevation emerges as the key terrain factor, with partial dependence analyses showing that high, steep, and well-drained landscapes are particular hotspots of vulnerability (Fig. 2e and f).

Focusing on baseline in-situ climatic factors (excluding ENSO effects), temperature had the strongest influence on resilience, followed by radiation, while precipitation had a relatively minor effect (Fig. 2a). These factors exhibited non-monotonic impacts, with resilience peaking under optimal climatic conditions (approximating climatological means, temperature = 275–280 K; precipitation = 50–100 mm/a; solar radiation = 100–150 W/m2), and declining when conditions deviated from this range (Fig. 2b, and Supplementary Fig. 3). Notably, at a global scale, low precipitation, high temperature, and elevated solar radiation had a stronger impact on vegetation resilience compared to high precipitation, low temperature, and low radiation (Fig. 2b).

Although in-situ factors collectively explain a larger proportion of the resilience variation, ENSO teleconnections remain a critical global-scale forcing, driving resilience anomalies. ENSO primarily influences vegetation resilience through the climate-mediated teleconnection, driven by ENSO-induced climate anomalies, which exert a stronger impact than ENSO’s atmospheric synchronization teleconnection (Fig. 2a). Specifically, reduced vegetation resilience is closely associated with anomalously high temperatures, low precipitation, and enhanced solar radiation during El Niño events, as well as with low-temperature and low-precipitation anomalies under La Niña conditions (Fig. 2d). While these climate-mediated effects dominate overall, the atmospheric synchronization teleconnections of ENSO still matter: vegetation responds more rapidly and strongly to El Niño events than to La Niña during the event year (Fig. 2b). However, the lagged response of vegetation to ENSO is often a result of regional disparities, including distance from the tropical oceans and local topography like elevation and terrain complexity. With increasing lag time, the global impact of extreme La Niña events becomes more pronounced and, by the third year (the lag = 2), surpasses that of El Niño in shaping vegetation resilience dynamics (Fig. 2c). This temporal difference likely stems from the contrasting durations of the two phenomena: La Niña often persists into the second year and re-strengthens in the winter, whereas El Niño tends to terminate shortly after reaching its peak stage39.

An independent analysis targeting intact forest regions (see Methods) indicates that after excluding areas affected by human activities, vegetation responses to ENSO remain broadly consistent with the global-scale patterns (Supplementary Fig. 4). This suggests that while anthropogenic disturbances can influence vegetation responses at local scales, the large-scale climate–vegetation coupling associated with ENSO persists.

ENSO synchronization teleconnection affecting vegetation resilience
We applied a Pearson Chi-square test to evaluate whether vegetation resilience is temporally synchronized with ENSO over a long historical period, after removing the effects of in-situ climate variables (see Methods). The resulting spatial pattern indicates that 26% of vegetated regions exhibited significant relationships between vegetation resilience anomalies and the ENSO events (P < 0.05), with resilience demonstrating higher sensitivity to El Niño than La Niña events (Fig. 3a). Specifically, 6% of vegetated areas showed increased resilience during La Niña events, while 3.3% exhibited decreased resilience. During El Niño events, 9.4% of areas experienced increased resilience and 7.4% experienced decreased resilience (Fig. 3a). The detrimental impacts of El Niño were predominantly concentrated in northern Asia, Europe and America, and the southern Amazon, whereas the negative effects of La Niña were primarily distributed in southern Amazon and Central Asia and Africa (Fig. 3a). Evergreen broadleaf forests, grasslands and open shrublands were identified as the vegetation types most affected by these climatic phenomena (Supplementary Fig. 5). In contrast, the positive impacts of El Niño and La Niña are more spatially dispersed, with affected vegetation resilience distributed across most continents (Fig. 3a). These spatial distributions remained consistent across varying lag times. However, it is noteworthy that the spatial extent of El Niño’s synchronous negative impact decreased with increasing lag time, while the area experiencing negative impacts from La Niña remained relatively stable. Beyond a twelve-month lag, La Niña’s negative impact area exceeds that of El Niño (Supplementary Fig. 6), whereas the global-scale Random Forest analysis indicates that La Niña’s cumulative impact becomes dominant only after three years (Fig. 2c). This apparent difference reflects the distinct metrics captured by the two approaches—spatial extent versus aggregated strength—yet both consistently point to La Niña’s growing dominance with increasing lag time.

ENSO climate-mediated teleconnection affecting vegetation resilience
RF inference (Fig. 2d) suggests that ENSO-induced climate anomalies trigger resilience changes. By overlapping resilience-climate (Supplementary Fig. 7) and ENSO-climate (Supplementary Fig. 8) linkages, we reveal spatially heterogeneous coupling patterns, where ENSO-driven climate anomalies negatively affected 23% of vegetation resilience, predominantly in the Amazon rainforest and Brazilian Plateau region in central South America, the Congo Basin in south-central Africa, and northeastern Australia (Fig. 3h). Conversely, positive effects were observed in 22% of regions, primarily in northern Australia, Sub-Saharan Africa, and the Southeast Plains of South America (Fig. 3h). Among all anomalous events, 59% of the positive impacts were attributed to El Niño and 41% were linked to La Niña, while 63% of the negative impacts were associated with El Niño and 37% with La Niña (Fig. 3f, g). Particularly, the resilience of evergreen broadleaf forests in the Amazon and savannas in the Brazilian Plateau were negatively affected by El Niño-induced anomalies, including decreased precipitation, increased temperature, and increased solar radiation (Fig. 3c–e and Supplementary Figs. 8 and 9). In contrast, La Niña-driven anomalies had considerable positive effects on the resilience of grasslands and open shrublands in northern Australia, although parts of the mountainous regions in northeastern Australia experienced negative resilience impacts from La Niña, driven by reductions in radiation, temperature, and precipitation (Fig. 3c-e and Supplementary Figs. 8 and 9). On a global scale, evergreen broadleaf forests, open shrublands, woody savannas, and savannas experienced the most pronounced negative impacts, whereas grasslands showed resilience linked to ENSO-driven solar radiation and precipitation variability (Fig. 3h and Supplementary Figs. 8 and 9).

Impact of ENSO on the long-term vegetation resilience trend
To analyze the impact of ENSO events on the long-term resilience trend, we compared the historical AC trend between ENSO-affected and non-ENSO-affected regions. ENSO-affected regions, defined as areas experiencing synchronization or climate-mediated impacts, cover 53% of global vegetation. Within these regions, 15% are influenced directly by ENSO’s synchronization, 51% are affected solely through ENSO-induced climate anomalies (for example, abnormal temperature, precipitation, and radiation), and the remaining 34% are subject to both effects. In terms of ENSO phases, 61% of affected regions are linked to El Niño, while 39% are associated with La Niña.

Regardless of considering ENSO as a whole or its individual phases, the historical decline in vegetation resilience was significantly more pronounced in negatively affected regions compared to positively affected and non-ENSO-affected regions (p < 0.05) (Fig. 3b and Supplementary Table 1–3). In contrast, the positive impact of ENSO did not lead to a corresponding increase in resilience. Specifically, no significant difference was observed between positively affected or non-ENSO-affected regions under the overall impact of ENSO or La Niña (Fig. 3b and Supplementary Tables 1 and 3). Moreover, the resilience decline in positively affected regions under the impact of El Niño was higher than that in non-ENSO-affected regions (Fig. 3b and Supplementary Table 2). This suggests that even in regions where ENSO initially enhances vegetation resilience, long-term instability caused by alternating climate stress events, asymmetry in impact intensity, and lagging effects may contribute to the “pseudo-gain” trap of the positive ENSO effect.

Future global vegetation resilience and its response to future changes in ENSO
Leveraging simulated LAI and climate projections from CMIP6, we evaluated future trajectories of vegetation resilience and examined their links to climate change across three emissions scenarios (SSP126, SSP245, and SSP370) (see Methods). Our results indicate that 33-39% of vegetated regions are anticipated to undergo a further decline in resilience across three future climate scenarios, with an average decline trend of 0.12 in these affected areas. Notably, Western Siberia and Northern Europe were identified as emerging areas of declining vegetation resilience in all three scenarios (Fig. 4a). Regarding ENSO teleconnections, projected warming is expected to weaken the synchronization impact of ENSO on vegetation resilience, decreasing from 21.4% under SSP126 to 21.1% under SSP245 and 13.5% under SSP370 (P < 0.05) (Fig. 4b). In parallel, ENSO-driven climate anomalies are projected to strengthen, affecting approximately 51–55% of global vegetation areas and marking a 6-10% increase compared to historical levels (Fig. 4c). These changes in ENSO-resilience patterns highlight newly vulnerable regions. For instance, under three scenarios, extensive regions across Siberia, as well as northern North America under SSP126 and SSP370, are projected to experience adverse impacts from La Niña-associated radiation and temperature anomalies (Fig. 4c and Supplementary Fig. 10). These findings underscore the increasing complexity of ENSO-vegetation interactions in the context of climate change and emphasize the urgent need for targeted research and conservation strategies in these newly identified vulnerable regions.

Methods
Historical vegetated and regional climate data
The historical vegetation dynamic was indicated by the long-term leaf area index (LAI) product of the Global Land Surface Satellite (GLASS) V5 dataset, which was generated from Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite data69. The dataset spans 1981–2018, with a temporal resolution of eight days and a spatial resolution of 0.05°70,71. After outlier (LAI > 10 or LAI < 0) removal, the LAI data were aggregated into monthly means. As a critical structural variable, LAI captures key vegetation processes such as photosynthesis, respiration, and transpiration72. Moreover, LAI is more resistant to saturation than Normalized Difference Vegetation Index (NDVI)72,73 and provides a longer and more consistent temporal record than Vegetation Optical Depth (VOD)74, making it a robust indicator for evaluating vegetation resilience over time. Regional climate data were obtained from the ERA5-land reanalysis dataset, which provides monthly near-surface (2 meters) air temperature (°C), monthly precipitation (m), and monthly surface solar radiation (j/m2). These variables were processed and downloaded via the Google Earth Engine platform, with surface solar radiation converted to W/m². All data were bilinearly aggregated to a 0.5° resolution, as is typically recommended75. ENSO variability was represented by monthly sea surface temperature anomalies (SSTA) in the Niño 3.4 region (5°S–5°N and 170°W–120°W), obtained from the NOAA ERSST V5 dataset76. All variables were analyzed over the period 1981–2018. Historical trajectories of these variables are shown in Supplementary Fig. 11.

To project future changes, we use the simulated climate and vegetation time series from the sixth phase of the Coupled Model Intercomparison Project (CMIP6)77, including monthly SST, LAI, near-surface air temperature, precipitation, and surface solar radiation. These projections were generated under the Scenario Model Intercomparison Project (SSP126, SSP245, SSP370), covering the period 2015 to 2100. The climate factors (precipitation, temperature, radiation, and SST) and LAI were derived as multi-model averages to mitigate potential model-specific biases and provide a more robust, ensemble-based estimate of future climate projections (the models are listed in Supplementary Tables 4–6). All variables were bilinearly aggregated to a horizontal grid of 1° × 1°.

ENSO variability was represented by the Monthly SSTA averaged over the Niño 3.4 region. Monthly SSTAs in future projected scenarios were constructed by removing the monthly climatology (historical monthly mean), and applying a quadratic detrending to remove the long-term warming trend78,79, then standardizing by dividing by historical monthly SSTA standard deviation63. The projected SSTA trajectories under different scenarios are shown in Supplementary Fig. 11.

Land cover and terrain data
We used land cover data from the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover type product (MCD12C1)80. The maps employed the International Geosphere-Biosphere Programme (IGBP) classification scheme, provided at yearly intervals with 0.05° resolution. For this study, we assumed that land cover remains unchanged, and we utilized the global distribution of vegetated land. We used the Global Intact Forest Landscapes (IFL) dataset for the year 200081, which provides a subset of the LAI time series, to assess the resilience of intact forests. Global elevation data in the year 2000 was derived from the Shuttle Radar Topography Mission (SRTM) Version 3 product, which offers a resolution of 30 m82. The slope was calculated based on this elevation data. The Height Above Nearest Drainage (HAND) data, obtained from Google Earth Engine83, were derived from SRTM elevation data and provide a more accurate representation of groundwater accessibility and hydrological buffering capacity84. Finally, all data were aggregated to a spatial resolution of 0.5° and 1°, respectively.

Climate zones data
The Köppen-Geiger World map is used to define the climate zones in our analysis. We aggregate the 31 climate zones into five major zones (tropical, arid, humid temperate, temperate, continental, and polar) (Supplementary Fig. 3).

Vegetation resilience calculation
Assessing vegetation resilience generally requires estimating the time needed to return to equilibrium85, but large-scale direct measurement is challenging. Natural systems experience continuous random disturbances86, and external shocks vary in intensity and frequency across regions87, making experimental assessments of recovery rates difficult. However, theoretical studies suggest that a system’s response to disturbances can be predicted through its internal variability37,88. The dynamics of an ecosystem can be modeled using the Ornstein-Uhlenbeck process89, using the following equation:

dXt = \theta(\mu - Xt)dt + \sigma dWt (1)

Where Xt is the system state variable, \mu is the steady-state mean of the system, \theta is the recovery rate, \sigma is the noise intensity, and dWt is the standard Wiener process. Discretizing this process over the timestep \Delta t, the equation is:

Xt+\Delta t = Xt e^{-\theta\Delta t} + \mu(1 - e^{-\theta\Delta t}) + \sigma \sqrt{\frac{1 - e^{-2\theta\Delta t}}{2\theta}} \epsilon t (2)

If the timestep \Delta t is small, the terms involving \mu and the scaling of noise can be approximated, leading to the model yielding the first-order autoregressive process67:

Xt+1 = \alpha Xt + \sigma\eta (3)

The lag-1 autocorrelation coefficient (\alpha) can be shown:

\alpha = e^{-\theta\Delta t} (4)

The variance (VAR) of the discretized time series can be shown:

VAR(X) = \frac{\sigma^2}{1 - e^{-2\theta\Delta t}} = \frac{\sigma^2}{2\theta} (5)

These equations indicate that as a system approaches a critical threshold, recovery \theta slows, leading to an increase in variance and AC. Thus, these metrics are considered to be early warning signals of ecosystem state shifts36.

In this study, lag-1 AC was used as a resilience indicator, and variance was also calculated to verify the robustness of AC. The lag-1 AC was estimated using a first-order autoregressive model fitted via least squares:

Xt = \alpha Xt+1 + \epsilon (6)

Where \alpha is the AC coefficient. \epsilon is the residual. Xt+1 is the lag-1 time series.

The variance was estimated as the average of the squared deviations from the mean:

VAR(X) = \frac{1}{n-1} \sum_{i=1}^{n} (Xi - \mu)^2 (7)

Where Xi is the ith element in the LAI subset, \mu is the mean of the subset, and n is the number of elements in the subset.

In this study, the GLASS LAI dataset and CMIP6-derived LAI data under three scenarios were used to calculate historical and future vegetation resilience. However, prior to measuring the resilience metrics, it is necessary to remove the impact of long-term trends of LAI time series. STL (Seasonal-Trend decomposition using LOESS) is a widely used time series decomposition method. It gradually extracts trend and seasonal components through an iterative smoothing process and finally obtains the residual90. We employed STL decomposition through the stl() function in the ‘stats’ package in R for each grid. This method robustly separates the LAI time series into three components: a long-term trend, a fixed seasonal cycle trend (by using ‘periodic’ option in s.window parameter), and a residual component91,92. Then, the residual components of the LAI time series were utilized to calculate the lag-1 AC and variance. The ar.ols() and var() function in ‘stats’ package in R were employed to calculate the lag-1 autocorrelation and variance on the sliding windows of 5 years. To investigate the effect of window size on the resilience trends, we also calculated AC with a sliding window of 10 years (Supplementary Fig. 12). Note that only the regions where AC > 0 were involved in this study. For the LAI data preprocessing, we did not fill the data gap because the relationship between the AC of the gap-free and gapped time series has been proven to be very close to being one-to-one67. We also excluded the grids of croplands (where human management may strongly disturb the response of crops to ENSO) and non-vegetation types, so that our analysis focuses on the vegetation of non-crop ecosystems.

Vegetation resilience and greenness trend
To assess the historical vegetation resilience trends, we calculated the Kendall tau value for the AC and the variance time series using a 5-year sliding window at each grid from 1981 to 2018. This rank-based correlation coefficient is independent of data distribution, allowing for consistent trend comparisons across regions. Using the cor.test() function with the Kendall method from the R package ‘stats’, we estimated resilience trends with a 95% confidence level. The trend directions of variance and AC were largely consistent across most regions (Fig. 1c). Additionally, we calculated the AC trend using a 10-year window, which indicated a spatial pattern qualitatively similar to that of the 5-year window (Fig. 1a and Supplementary Fig. 12). AC trends under future scenarios were also calculated (Fig. 4).

Long-term LAI trend was quantified using Kendall’s tau at p < 0.05 to represent vegetation greenness change (Fig. 1d). We then overlaid the AC and LAI trend maps to assess the spatial coupling between greenness and resilience.

Random forest model
We utilized random forest regression to identify the nonlinear relationship between long-term vegetation resilience dynamics and environmental factors. Random forest regression is a robust machine-learning model that can decouple the interactions and identify key drivers of change93. By integrating spatial features to capture regional differences and time series data to address temporal changes, it effectively analyzes global patterns of vegetation resilience under climate change94. In this research, the response variable was the annual time series of AC, which characterizes vegetation resilience. Predictor variables included a combination of climatic and terrain factors. Climatic predictors consisted of the annual mean near-surface (2 m) air temperature, annual mean of monthly total precipitation, and annual mean of monthly accumulated surface shortwave radiation, as well as the annual ENSO index. The annual ENSO index was represented by the mean sea surface temperature anomaly (SSTA) over the Niño 3.4 region, calculated as the average over December–February (DJF) using the NOAA ERSST v5 dataset to capture the peak phase of ENSO events17. All climatic variables were computed at the pixel level. The regression model incorporated an interaction term between ENSO and climate factors, as well as a lag term. Climate factors from the preceding three years represented the lagged effect, while the interaction terms (temperature × ENSO, precipitation × ENSO, and radiation × ENSO) characterized how resilience sensitivity to climate changes varies with ENSO intensity95,96. DEM, slope, and HAND were included in the model as the terrain predictors to control potential interference with terrain background heterogeneity. In addition, intra-annual variations in temperature and precipitation may influence vegetation resilience by favoring species with distinct adaptive traits. To account for potential differences in resilience due to climatic seasonality, we included multi-year mean seasonality indices (monthly coefficient of variation) of precipitation, radiation, and temperature as covariates in the model. We employed the ‘ranger’ package in R to train the model with 40% of the global vegetated pixels and validated it using another 10% of pixels, both randomly selected after stratified allocation based on vegetation type proportions. The number of decision trees was set as 500, as the model’s performance had reached stability at this level, and the number of splits was determined as the square root of the number of selected predictors. We then plotted partial dependence plots to illustrate the relationships between AC and the predictors and their marginal contribution to the AC trend. To examine whether human activities influence the observed climate–resilience coupling, we used intact forest ecosystems—largely free from direct land-use impacts—as a control group and trained an additional RF model using data exclusively from these areas (Supplementary Fig. 4).

ENSO’s impact on resilience
To map the spatial distribution of ENSO’s synchronization impact on vegetation resilience, we first removed the effects of in-situ climate by fitting, at each grid cell, a multiple linear regression of resilience on precipitation, temperature, and radiation. We then used a Chi-square test, a nonparametric test method suitable for categorical variables97, to examine whether El Niño or La Niña events triggered anomalies in adjusted vegetation resilience. Positive and negative anomalies in vegetation resilience were defined as standardized AC anomalies (Z-scores) exceeding +2 or falling below –2, respectively, corresponding approximately to the upper and lower 5% tails of the distribution. Similarly, the definition of El Niño or La Niña events was defined based on the commonly accepted criterion: a monthly SSTA greater than 1°C or less than –1°C98,99. El Niño and La Niña were tested separately. When both were significant for a cell (FDR-adjusted P < 0.05), we labeled the dominant phase as the one showing the stronger association. To analyze the lagged effect, Chi-square coefficients were calculated between the ENSO event time series and AC anomaly time series in a lag of 3, 6, 9, 12, and 24 months (Supplementary Fig. 6). We also performed a similar analysis for the three scenarios in the future period (2015–2100) (Fig. 4b). This approach helps reveal the spatial imprint of ecosystem responses that are temporally aligned with ENSO rhythms, including delayed effects, beyond what can be explained by in-situ climate variability.

We then assessed the climate-mediated impact of ENSO, representing the portion of ENSO’s influence on resilience that is mediated through in-situ climatic variables (precipitation, temperature, and solar radiation). This was quantified by evaluating how ENSO-driven climate anomalies affected vegetation resilience. Linear regression was used to measure the co-variation between climatic variability and AC dynamics (Supplementary Fig. 7), and Chi-square tests were applied to assess the association between ENSO events and climatic anomalies (P < 0.05, Supplementary Fig. 8). By spatially overlaying these relationships, the impact on resilience was coded as either positive (“1”) or negative (“−1”). For example, if resilience and precipitation were positively correlated, and ENSO events predominantly caused negative precipitation anomalies in the region, the ENSO-precipitation impact on resilience was recorded as “−1”. We summed the climate-mediated effects of key ENSO–climate variables to quantify their overall impact (Fig. 3c–h). This analysis was also conducted for future scenarios (Supplementary Fig. 10).

Statistical significance test of regional AC trend
Similar to the climate-mediated impacts, we also quantified the positive and negative atmospheric synchronization impacts as “+1” and “−1”, respectively, and then summed both synchronization (Fig. 3a) and climate-mediated (Fig. 3h) impacts to derive the total impact. Based on these total impacts, we classified the global vegetated areas into three categories for each ENSO type (El Niño, La Niña, and overall ENSO): (i) ENSO-positively affected regions, where the sum of synchronization and climate-mediated effects is greater than zero; (ii) ENSO-negatively affected regions, where the sum is less than zero; and (iii) no ENSO impact regions, where the sum equals zero. Pairwise comparisons were then conducted between these categories using Student’s t-test and the Mann-Whitney U test11 (Fig. 3b). In each test, 5000 samples were randomly selected from each category, and both tests were repeated 10 times to obtain multiple P-values. If the majority of P-values for a given comparison were significant (P < 0.05), the Benjamini-Hochberg correction was applied to control for Type I errors (false positives) due to multiple comparisons100,101. Results were considered statistically significant when both the raw and adjusted P-values were below 0.05. The statistical test results are presented in Supplementary Tables 1, 2, and 3.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Major comments:\nThe authors investigate the direct and indirect impacts of ENSO on past and future global vegetation resilience using\nsatellite observations and model projections. Although the research is very relevant, the methods sections is currently\nincomplete for full reproducibility of the results. Furthermore, the discussion does not include enough actual discussion of\nthe results. There is no comparison of the results with current literature, nor implications of what these results mean.\nMinor comments:\nIntroduction:\nThe introduction could use some more background information on what ENSO is.\nLine 73: Critical slowing down is more often used than critical slowdown.\nLine 83: Remove the extra word: ‘These Our findings …’.\nMethods:\nThere is no reference for the GLASS LAI dataset included, nor what kind of data GLASS is (combination of AVHRR and\nMODIS).\nFor reproducibility, the authors should include more information about cloud masking/outlier definition/data masking.\nLine 310: The dataset has been moved to a new webpage according to the link, so better to update the link.\nThere is no explanation about the use of LAI rather than VOD/NDVI/… which have been used more often in global resilience\nstudies.\nLine 321-323: ‘Other climate factors were derived exclusively from the CESM2 model, as its historical simulations (1981–\n2014) showed a strong agreement with observation (EAR5)…’ I don’t understand this, why the use of another model? Does\nCMIP6 not provide these variables? I think EAR5 should be ERA5? Add which climate factors are referred to here.\nLine 333: Each MODIS product has a data citation, so replace the link by the actual citation.\nLine 337: The assumption of no land cover change from 2000 to 2100 seems quite unrealistic, so I would at least add some\nexplanation or a comment that this is a strong simplification here.\nLine 412: Why are you using yearly mean AC time series if you previously used a 5-year sliding window on monthly values?\nLine 413: Which ENSO indices are referred to here? How were they calculated? Why only December-February?\nLine 419: Is one year enough to capture all climate legacy effects? Is it necessary to include impact of previous years?\nLine 420: I would expect very strong correlations between ENSO and these climatic variables, is that not problematic for the\nmodel?\nLine 425: Were these pixels randomly chosen or did you take the percentage of different land cover types into account?\nLine 436: How robust is this threshold of -1 or 1 for negative or positive anomalies? Do the results change a lot when\nadapting this threshold?\nResults:\nThe text is not consistent in its terminology: sometimes 1-lag AC and sometimes lag-1 AC, better to be consistent.\nFigure 1e: It would be more consistent to use capital letters for delta_ac and delta_lai in the figure and the caption.\nFigure 2f: Spelling mistake in x axis label. It is also not clear how to interpret the values on the x axis: e.g. what does an\nENSO-temperature value of -200 mean?\nFigure 3: Map a is too small to distinguish between the 7 different colours. Why is (b) in the caption listed after (c-h)? Why\nare there boxes with coloured stripes included in each of the graphs in (b)?\nLines 195-201 and figure 3b: I’m confused with the explanation on which differences are significant and which ones are not,\nbecause all three graphs in figure 3b look the same to me.\nLine 2019: How was the projected vegetation resilience calculated from the CMIP6 data?\nFigure 4b: the maps are too small to distinguish between the 7 different colours.\nDiscussion:\nThe background information about ENSO in the second paragraph belongs in the introduction, rather than in the discussion.\nThe discussion does not include any discussions or conclusions about the results, only that the research did not include\nanthropogenic effects and that there was a lot of uncertainty in the data.", "reviewer_2": "Review of the manuscript: “ENSO amplifies global vegetation resilience variability in a changing climate”\nGeneral Assessment\nThis manuscript presents an innovative global analysis of vegetation resilience using 1-lag autocorrelation (AC) of Leaf Area\nIndex (LAI) to assess the influence of ENSO events (El Niño and La Niña). The authors also examine how these effects\ninteract with average climate conditions. The approach is original and the global scope is commendable. However, the\nmanuscript would benefit from greater clarity in the Results section, particularly in explaining the analytical steps, and a\ndeeper discussion of the ecological mechanisms underlying regional differences in resilience. This manuscript offers a\npromising and novel analysis of ENSO impacts on vegetation resilience. To fully realize its potential, the paper needs\nclearer methodological descriptions, more detailed results interpretation, and a stronger discussion of underlying\nmechanisms and ecological context. Including a stratified analysis of intact vs. disturbed vegetation would also significantly\nstrengthen the study’s conclusions. Below, I provide detailed suggestions to enhance the clarity and impact of the\nmanuscript.\nComments\nClarify the analytical framework throughout the Results section:\nThe Results section of the manuscript currently lacks sufficient context on how the analyses were performed. Please clearly\nstate the type of analysis in each subsection. This will help readers follow the sequence and logic of the findings. Please\ngive more context in the results, including the type of analysis that you performed for each section of the paper. This will help\nthe readers follow what you have done. For example:\nFigure 1a shows the spatial distribution of 1-lag autocorrelation (AC) trend. Blue values are negative and red values are\npositive. What does it mean? You mentioned in the caption “An increase in AC or variance reflects a decline in vegetation\nresilience”, so I assume that red is a decline in vegetation resilience?\nLine 91 - “The most significant declines in resilience 92 are notably observed in northern North America, the Amazon,\nnorthern Europe, eastern Australia, 93 and northern and Southeast Asia ”. What are these results based on? How did you\ncompare these areas to get to the conclusion that these areas are particularly affected?\nFigure 1c- “(c) Trend agreement between two indicators. ”. Please state what indicators these are.\nLine 97 - “Regions such as the Amazon, south-central Africa, Southeast Asia, and western Oceania have experienced a\nwidespread increase in greenness—measured by rising LAI— over the last four decades, yet these areas have concurrently\nwitnessed significant declines in resilience ”. Please also explain how you got these results in the text.\nLine 219 - “Our results indicate” also really requires explaining how you did the analysis.\nIt is unclear what form of temperature you used (e.g., mean annual temperature, maximum temperature). Please specify this,\nand similarly for precipitation and radiation.\nHave you considered climate seasonality (e.g., coefficient of variation or seasonality index of precipitation/temperature)?\nSeasonality could influence vegetation resilience by selecting for species with different adaptive strategies.\nFigure 2a highlights slope as an important predictor, but this result is not discussed in the main text. Please include a\nparagraph interpreting this result—how might slope influence resilience (e.g., via soil depth, drainage, erosion)?\nConsider including terrain-based indicators such as Topographic Position Index (TPI) or Height Above Nearest Drainage\n(HAND), which can better reflect groundwater access and hydrological buffering. This could provide a mechanistic link\nbetween terrain, water availability, and resilience. For reference, see: https://www.nature.com/articles/s41586-024-07568-w\nAddress the role of disturbance more directly:\nOne concern is the inclusion of both disturbed and undisturbed vegetation in the same analysis. These systems may have\nfundamentally different dynamics, especially in their response to ENSO events. I recommend conducting a separate\nanalysis for intact forests, using available datasets such as the GLAD Intact Forest Landscapes map:\nhttps://glad.umd.edu/dataset/intact-forest/overview . This could reveal whether observed resilience declines are driven more\nby land-use history or by climate sensitivity, and it would strengthen the paper’s conclusions.\nExpand the discussion of mechanisms:\nThe discussion would benefit from a more thorough treatment of why some regions are more sensitive to El Niño/La Niña\nthan others. Consider integrating more discussion on potential drivers such as species composition and functional traits, soil\nwater storage capacity, historical exposure to climate variability and human impacts and fragmentation\nBest,\nMatheus Nunes"}, "author_response": "Responses for Manuscript: ENSO amplifies global vegetation\nresilience variability in a changing climate\nWe greatly appreciate the reviewers’ valuable and constructive suggestions. We have carefully\nrevised our manuscript based on their comments. We hope that the revisions adequately address any\nconcerns. Below, you will find our item-by-item responses to the comments. The following different\nfonts are used to facilitate the review.\n Contents underlined are the original comments of the Reviewers.\n Contents in normal Times New Roman font are our response to the Reviewers.\n Contents in blue and italic font are part of the revised manuscript, bold text indicates newly\nadded content in direct response to reviewer comments.\n Green text highlights where the revisions have been made in the manuscript.\n Red-labeled figures are supplementary illustrations included in the response.\nContent\nResponses to Reviewer #1 ............................................................................................................... 3\nComment 1 ................................................................................................................................ 3\nComment 2 ................................................................................................................................ 5\nComment 3 ................................................................................................................................ 6\nComment 4 ................................................................................................................................ 6\nComment 5 ................................................................................................................................ 6\nComment 6 ................................................................................................................................ 7\nComment 7 ................................................................................................................................ 7\nComment 8 ................................................................................................................................ 8\nComment 9 .............................................................................................................................. 11\nComment 10 ............................................................................................................................ 11\nComment 11 ............................................................................................................................ 13\nComment 12 ............................................................................................................................ 15\nComment 13 ............................................................................................................................ 15\nComment 14 ............................................................................................................................ 17\nComment 15 ............................................................................................................................ 18\nComment 16 ............................................................................................................................ 18\nComment 17 ............................................................................................................................ 20\nComment 18 ............................................................................................................................ 20\nComment 19 ............................................................................................................................ 22\nComment 20 ............................................................................................................................ 23\nComment 21 ............................................................................................................................ 25\nComment 23 ............................................................................................................................ 27\nComment 24 ............................................................................................................................ 27\nComment 25 ............................................................................................................................ 28\nResponses to Reviewer #2 ............................................................................................................. 30\nComment 1 .............................................................................................................................. 30\nComment 2 .............................................................................................................................. 32\nComment 3 .............................................................................................................................. 33\nComment 4 .............................................................................................................................. 33\nComment 5 .............................................................................................................................. 34\nComment 6 .............................................................................................................................. 34\nComment 7 .............................................................................................................................. 35\nComment 8 .............................................................................................................................. 37\nComment 9 .............................................................................................................................. 38\nComment 10 ............................................................................................................................ 40\nComment 11 ............................................................................................................................ 42\nThank you very much for your thoughtful comments. We are grateful for your recognition of our\nstudy, and we appreciate your constructive critiques regarding the methodological limitations and\nthe practical implications of the discussion. To address your concerns in detail, we have organized\nour responses by specific points. Below, we provide a point-by-point reply, and we hope our\nrevisions have sufficiently addressed the issues raised.\nComment 1\nThe authors investigate the direct and indirect impacts of ENSO on past and future global vegetation\nresilience using satellite observations and model projections. Although the research is very relevant,\nthe methods sections is currently incomplete for full reproducibility of the results. Furthermore, the\ndiscussion does not include enough actual discussion of the results. There is no comparison of the\nresults with current literature, nor implications of what these results mean.\nReply: We have thoroughly revised and supplemented the Methods section to ensure full\nreproducibility of the results. Specifically, we have:\n1. Provided more detailed descriptions of all datasets used, including the GLASS LAI product\n(Comment 5), the updated SSTA data source (Comment 6), and MODIS landcover data\n(Comment 9). We also explained the rationale for selecting LAI as the vegetation indicator\n(Comment 7).\n2. Clarified and refined all data processing steps, including the calculation of future climate\nvariables from multi-model CMIP6 means (Comment 8) and the selection of sample points for\nthe RF model (Comment 15).\n3. Conducted sensitivity analyses to test the robustness of our results, including: (i) comparing\nresilience trends when using LUH2-projected 2050 landcover versus a unified MODIS\nlandcover mask for both historical and future periods (Comment 10), and (ii) assessing the\ninfluence of different thresholds for defining ENSO-related anomalies and their spatial\ndistributions (Comment 16). Both tests confirmed that our main conclusions remain unchanged.\n4. Added justifications for analytical choices, such as using annual mean AC for RF analysis\n(Comment 11), defining and computing the annual ENSO index (Comment 12), and\nelaborating on the RF model’s ability to disentangle correlated variables (Comment 14).\n5. Expanded the scope of the analysis by incorporating longer lag effects (Comment 13),\nstandardizing the terminology for lag-1 autocorrelation (Comment 17), and further enhancing\nmethodological transparency.\nIn addition, we have included a new Data Availability statement with direct download links to all\ndatasets, and uploaded part of the calculation code to Figshare to further facilitate reproducibility.\nwe have substantially expanded this section to provide a more in-depth interpretation of our results.\nWe have added a more explicit comparison of our findings with previous studies, thereby situating\nour results within the broader literature. We also emphasize the broader implications of our work,\nhighlighting the importance of embedding ENSO variability into ecosystem management\nframeworks, with special attention given to tropical and subtropical regions that are especially\nsensitive to climate anomalies.\nThe revised paper now reads:\n(Pages 9-10 Lines 276-327, Discussion)\nThis research provides a comprehensive analysis of the climatic drivers of vegetation resilience,\nuniquely incorporating the remote impacts of global climate oscillations such as ENSO. Using a\nRF model, we disentangled the effects of local climate variables (temperature, precipitation, and\nsolar radiation) from those of ENSO and their interactions, revealing how ENSO shapes\nresilience patterns across global, regional, and biome-specific scales. Projections further indicate\nevolving impacts of ENSO under future climate scenarios. One key finding is that 51% of global\nvegetation over the past four decades has been influenced by either direct ENSO effects or ENSO–\nclimate interactions. These impacts are both asymmetric and delayed, with El Niño and La Niña\nevents exerting opposing yet destabilizing effects on vegetation resilience, contributing to long-\nterm declines in resilience, most notably in the Amazon, Southeast Asia, Central Africa, southern\nNorth America, and Australia (Fig. 3).\nThese regional hotspots of ENSO sensitivity can be attributed to a complex interplay of climatic\ndrivers, ecological constraints, and human activities. ENSO-driven anomalies in heat and\nmoisture show pronounced spatial variation, exposing ecosystems to different types and\nintensities of extreme events (e.g., heatwaves, droughts, floods) that impair photosynthesis, water\nregulation, and root function1,2. These effects are most severe in ecological hotspots such as\ntropical forest, savanna and grassland. In addition, our model highlights topographic features—\nparticularly high elevation, steep slopes, and well drained regions—as key predictors of low\nresilience and strong ENSO responses (Fig. 2f). These landscapes often support vegetation with\nshallow roots, limited soil depth, and poor water and nutrient retention, resulting in reduced\nbuffering capacity3. In parallel, climatic seasonality drives divergent vegetation adaptive\nstrategies, which in turn influence regional sensitivity to ENSO. While low-latitude regions with\nweak seasonality exhibit strong resilience shifts globally (Fig. 2e), ecosystems in areas with high\nintra-annual precipitation variability—such as northern Asian shrublands, Siberian tundra,\neastern Amazon rainforest, and Sahelian grasslands—also show heightened sensitivity compared\nto adjacent regions with similar rainfall and vegetation types (Fig. 3a). This supports prior\nfindings1,4,5 that stable water availability enhances post-disturbance recovery.\nEqually important, anthropogenic disturbances exacerbate these patterns6. Highly fragmented\ntropical and temperate forests, such as those in the Amazon, Central Africa, and parts of North\nAmerica and Southeast Asia, exhibit heightened sensitivity to ENSO (Fig. 3a), likely due to\nintensified edge effects that increase local temperatures and reduce humidity, amplifying climate\nstressors7-9. Intensive agriculture drives land degradation, biodiversity loss, and hydrological\ndisruption, further destabilizing adjacent ecosystems10,11. For instance, we found that mixed\nforests along the Eurasian agro-ecological transition zone show more negative responses to La\nNiña when located near croplands (Fig. 3a), suggesting landscape configuration modulates\nENSO sensitivity—a mechanism that merits further empirical investigation.\nOur findings reinforce previous evidence that vegetation in climatically unstable environments—\ncharacterized by high variability in light, temperature, and moisture—tends to exhibit lower\nresilience. The identified hotspots of ENSO sensitivity are broadly consistent with earlier studies\non vegetation–climate responses12. For example, we found that vegetation resilience is notably\nlower in the southeastern Amazon, consistent with earlier studies that associate this decline with\nrecurring droughts and floods13,14. Similarly, reduced resilience along the Sahelian margins\nsupports earlier findings on the role of increasing water limitation15. ENSO-related effects\nobserved in Southeast Asia, Australia, Central Africa, and Central America also reflect patterns\ndocumented in prior research16-18. As climate change intensifies, such instability is expected to\nincrease, threatening ecosystem stability across large regions4. It is therefore critical to integrate\nENSO forecasts into ecosystem restoration and climate-informed conservation planning,\nespecially for tropical and subtropical forests. These ecosystems, given their high biodiversity,\ncarbon storage potential, and sensitivity to both climate and human disturbance, should be\nprioritized in global strategies for resilience and adaptation. In contrast, some northern forests,\nsuch as mixed forests in Northern Europe, show increasing resilience and minimal ENSO\nsensitivity (Fig. 1a and Fig. 3a), highlighting the need for regionally differentiated adaptation\nstrategies.\nComment 2\nThe introduction could use some more background information on what ENSO is.\nReply: We have revised the Introduction to include a concise definition of ENSO as a coupled\nocean–atmosphere system characterized by periodic fluctuations in sea surface temperature and\natmospheric pressure in the tropical Pacific. We have also added background on the distinct climatic\neffects of its two phases—El Niño and La Niña—on regional temperature, precipitation, and solar\nradiation, as well as the resulting spatial heterogeneity in vegetation responses. These revisions aim\nto provide readers with sufficient context to understand the relevance of ENSO to vegetation\nresilience.\nThe revised paper now reads:\n(Page 2 Lines 46-65, Introduction)\nENSO is a coupled ocean–atmosphere phenomenon involving periodic fluctuations in sea\ntemperature and atmospheric pressure in the tropical Pacific19. It typically develops during the\nboreal summer and decays in the following spring20,21, with peak intensity usually occurring\nbetween December and February22. Through large-scale atmospheric circulation anomalies,\nENSO alters regional patterns of temperature, precipitation, and solar radiation23,24, thereby\ndisrupting vegetation biophysical processes. ENSO25. Specifically, the alternating warm (El Niño)\nand cold (La Niña) phases of ENSO exert complex and spatially heterogeneous vegetation\nresponses26. El Niño typically raises temperature and solar radiation in tropical and subtropical\nregions, yet its impact on precipitation varies, causing wetter conditions in some areas and\ndrought in others, potentially exacerbating or alleviating vegetation water stress27-29. For instance,\npatterns, decreasing precipitation in West and northern Central Africa, while increasing it in East\nand southern Central Africa30, thus creating a contrasting north-south gradient in vegetation\ndynamics. Conversely, La Niña generally exhibits the opposite climatic signals, particularly over\ntropical land areas where temperature and precipitation patterns respond in antiphase, thereby\npartially offsetting the effects of El Niño31,32. However, this compensatory mechanism is modulated\nby variation in the location of equatorial warming and by nonlinear ocean–atmosphere interactions,\nleading to spatial asymmetries in the timing and intensity of vegetation response to ENSO phases33.\nConsequently, some regions, such as the Amazon rainforest, are highly sensitive to El Niño34,35,\nwhile others, like Australian shrublands, experience alternating impacts from both phases36.\nComment 3\nLine 73: Critical slowing down is more often used than critical slowdown.\nReply: We have replaced “critical slowdown” with the more commonly used term “critical slowing\ndown” in line 80-82 and checked the whole manuscript for consistent usage.\nThe revised paper now reads:\n(Page 2 Lines 80-82, Introduction)\nThe critical slowing down theory (CSD) provides a framework for detecting ecosystem instability\nby suggesting that declining recovery rates precede state transitions, which can be identified by\nstatistical signals like increased variance and higher first-order autocorrelation.\nComment 4\nLine 83: Remove the extra word: ‘These Our findings …’.\nReply: We have removed the extra word “Our” in the sentence.\nThe revised paper now reads:\n(Page 3 Lines 90-92, Introduction)\nThese findings shed new light on the teleconnection of vegetation growth and extreme climate,\nsupporting nature-based solutions and informing sustainable development efforts.\nComment 5\nThere is no reference for the GLASS LAI dataset included, nor what kind of data GLASS is\n(combination of AVHRR and MODIS). For reproducibility, the authors should include more\ninformation about cloud masking/outlier definition/data masking.\nReply: We have now included a more detailed description in the Methods section. Specifically, we\nused the GLASS LAI Version 5 dataset (1981–2018), which is derived from AVHRR data. We\nselected this version due to its longer temporal coverage compared to the more recent Version 6\n(2000–2020).\nWe also added details on data processing as follows:\n(1) Outlier definition: Pixels with LAI values less than 0 or greater than 10 were considered outliers.\n(2) Cloud masking: The GLASS LAI dataset has undergone preprocessing to remove cloud\ncontamination and has no gaps or missing values. Due to its high spatial and temporal continuity,\nno additional cloud masking was applied.\n(3) Data masking: We used the MODIS land cover product to exclude non-vegetated and cropland\nareas, as described in the Vegetation resilience calculation section of the Methods (Page 12, Lines\n442-444, Methods).\nLastly, we added the official download link and citation information for the GLASS LAI dataset.\nThe revised paper now reads:\n(Page 10, Lines 340-345, Methods)\nThe historical vegetation dynamic was indicated by the long-term leaf area index (LAI) product\nof the Global Land Surface Satellite (GLASS) V5 dataset, which was generated from Advanced\nVery High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging\nSpectroradiometer (MODIS) satellite data (http://www.glass.umd.edu/Download.html). The\ndataset spans 1981–2018, with a temporal resolution of eight days and a spatial resolution of\n0.05°37,38. After outlier (LAI >10 or LAI <0) removal, the LAI data was aggregated into monthly\nmeans.\nComment 6\nLine 310: The dataset has been moved to a new webpage according to the link, so better to update\nthe link.\nReply: We have updated the link in the revised manuscript.\nThe revised paper now reads:\n(Page 10, Lines 354-356, Methods)\nENSO variability was represented by monthly sea surface temperature anomalies (SSTA) in the Niño\n3.4 region (5°S – 5°N and 170°W – 120°W), obtained from NOAA ERSSTV5 dataset\n(https://psl.noaa.gov/data/timeseries/month/).\nComment 7\nThere is no explanation about the use of LAI rather than VOD/NDVI/… which have been used\nmore often in global resilience studies.\nis less susceptible to saturation in high-biomass regions such as tropical forests, making it more\nsuitable for capturing vegetation structural dynamics. Although VOD provides complementary\ninformation on vegetation water content, available VOD datasets typically have coarser spatial\nresolution and shorter time spans (e.g., AMSR-E: 2002–2011; SMOS: 2010–present; SMAP: 2015–\npresent). While the fused VODCA dataset extends back to 1987, its construction from multiple\nsensors can introduce inter-sensor inconsistencies. Recent studies have shown that such merging\nmay bias critical slowing down indicators39, potentially affecting resilience detection. In contrast,\nthe LAI dataset we used provides a spatially continuous, gap-free record since the early 1980s,\noffering both high resolution and temporal consistency essential for long-term resilience assessment.\nMoreover, LAI directly reflects canopy structure and leaf biomass, making it more mechanistically\nlinked to ecological functions such as photosynthesis and transpiration. These advantages make LAI\nwell-suited for capturing vegetation recovery signals and interannual variability across diverse\nbiomes. We have added this explanation in the Methods section and cited relevant literature\naccordingly.\nThe revised paper now reads:\n(Page 10, Lines 345-349, Methods)\nAs a critical structural variable, LAI captures key vegetation processes such as photosynthesis,\nrespiration, and transpiration40. Moreover, compared to other indices such as Normalized\nDifference Vegetation Index (NDVI) and Vegetation Optical Depth (VOD), LAI is more\nsaturation resistant and provides longer temporal coverage, making it a robust indicator for\nevaluating vegetation resilience over time.\nComment 8\nLine 321-323: ‘Other climate factors were derived exclusively from the CESM2 model, as its\nhistorical simulations (1981–2014) showed a strong agreement with observation (EAR5)…’ I\ndon’t understand this, why the use of another model? Does CMIP6 not provide these variables?\nI think EAR5 should be ERA5? Add which climate factors are referred to here.\nReply: “EAR5” should indeed be “ERA5.”\nIn the original manuscript, “other climate factors” referred to precipitation, temperature, and\nradiation. We initially selected the CESM2 model (a CMIP6 model) to represent these variables\nbecause its historical simulations (1981–2014) showed relatively high temporal consistency with\nERA5 observational data across our study region. Based on this agreement, we used CESM2 outputs\nto characterize future changes in these variables.\nHowever, following your suggestion, we recognize that using a single model could introduce\n8\nsystematic bias. To ensure methodological consistency and reduce model-specific uncertainties, we\nhave now replaced the single-model approach with a multi-model ensemble mean, following the\nsame strategy used for other climate predictors. The revised methods and results sections have been\nupdated accordingly.\nIn line with these revisions, Fig. 4, Extended Fig. 9, and Supplementary Table 6 have also been\nupdated.\nThe revised paper now reads:\n(Page 8, Lines 259-261, Results)\nIn contrast, the influence of ENSO-climate interactions is anticipated to strengthen, with\napproximately 51-55% of global vegetation areas experiencing indirect impacts, marking a 6-10%\nincrease compared to historical levels (Fig. 4c)\n(Page 8, Lines 268-274, Results)\nFig. 4 Projected vegetation resilience trend and ENSO’s direct and interactive impact on\nvegetation resilience under future scenarios (2015–2100). (a) Resilience trend under SSP 126, SSP\n245, and SSP 370, respectively. (b) Direct impact of ENSO on vegetation resilience in the three\nscenarios. (c) Impact of ENSO-climate interaction on vegetation resilience in three scenarios. Plus\n(+) and minus (-) signs indicate positive and negative impacts, respectively, with numbers indicating\nimpact intensity, calculated as the sum of individual ENSO-climate interaction impacts.\n(Pages 10-11, Lines 365-368, Methods)\nThe climate factors (precipitation, temperature, radiation and SST) and LAI were derived as\n9\nmulti-model averages to mitigate potential model-specific biases and provide a more robust,\nensemble-based estimate of future climate projections (the models are listed in Supplementary\nTables 4-6).\n(Page 24, Lines 803-811, Extended Data)\nExtended Data Fig. 9. Distribution of the predicted impacts of ENSO-climate interactions on\nglobal resilience under the SSP126 (a), SSP245 (b), and SSP370 (c). In each SSP scenario, the\nfirst two figures show the total impact of three climate factor anomalies—temperature, solar\nradiation, and precipitation)—caused by El Niño and La Niña events on resilience, respectively. The\nlast three figures show the individual impacts of each climate anomalies induced by ENSO on\n10\nresilience. The plus (+) and minus (-) signs in the figure indicate positive and negative impacts on\nresilience, respectively, while different colors represent impact intensity. In the total impact figures,\nintensity is defined by the number of individual ENSO-climate anomalies influencing resilience.\n(Page 25, Lines 812-815, Extended Data)\nExtended Data Fig. 10. Predicted future trajectories of AC and annual mean climatic factors for\nglobal vegetation regions (2015-2100) under three CMIP6 scenarios. The scenarios include:\nSSP126, SSP 245, and SSP370.\nComment 9\nLine 333: Each MODIS product has a data citation, so replace the link by the actual citation.\nReply: In accordance with the data citation guidelines, we have replaced the web link with the\nformal citation for the MODIS MCD12C1 land cover product in the revised manuscript.\nThe revised paper now reads:\n(Page 11, Lines 377-378, Methods)\nWe used landcover data from Moderate Resolution Imaging Spectroradiometer (MODIS) land cover\ntype product (MCD12C1) (https://doi.org/10.5067/MODIS/MCD12C1.061).\nComment 10\nLine 337: The assumption of no land cover change from 2000 to 2100 seems quite unrealistic, so I\nwould at least add some explanation or a comment that this is a strong simplification here.\nReply: In this study, a consistent land cover dataset was used across the historical and future periods.\nFor the future scenarios, land cover was only applied as a vegetation mask rather than for vegetation\ntype-specific analyses. This simplification was necessary because available future land-use\n11\nprojections that consider climate change—such as LUH2 (Land-Use Harmonization 2)—adopt\nrelatively coarse vegetation classifications (e.g., forest, grassland, pasture) and do not explicitly\ninclude non-vegetated categories such as water bodies, barren land, or glaciers, which are present\nin current observational land cover products.\nNevertheless, we acknowledge this as a strong simplification. To evaluate its potential impact, we\nrefer to a recent study41 that reconstructed land cover classifications for the SSP245 scenario in 2050\nusing the MAgPIE model based on LUH2 data. The resulting land cover map (Fig. 1 and 2) was\nused to redefine the vegetation mask and recalculate the spatial pattern of projected vegetation\nresilience. The results showed that the large-scale spatial trends were generally consistent with those\nbased on the year 2000 vegetation mask, supporting the robustness of our main findings despite the\nland cover simplification.\nFig. 1. Simulated global vegetation cover type for the year 2050 under SSP245. The Grey areas\nindicate non-vegetation cover that is removed from our analysis a priori.\nmasked using two land-cover datasets. Resilience trends are calculated for each SSP scenario and\nmasked using either (a) MODIS land-cover in 2000 or (b) simulated future land-cover in 2050. Red\nareas indicate significant increases in autocorrelation (AC) or variance, suggesting reduced\necosystem resilience. Trends are calculated using Kendall’s tau.\nComment 11\nLine 412: Why are you using yearly mean AC time series if you previously used a 5-year sliding\nwindow on monthly values?\nReply: As described in the Methods section, the 5-year sliding window is a necessary step in\nconstructing the monthly autocorrelation (AC) time series. Within each window, we remove both\nthe long-term trend and the seasonal trend of vegetation to obtain a residual time series, from which\nthe autocorrelation coefficient is calculated. This approach ensures that the resilience indicator is\nnot confounded by long-term change or seasonality.\nDue to computational constraints associated with pixel-wise Random Forest (RF) modeling at\nthe global scale, we averaged the monthly AC values within each year to derive an annual AC time\nseries for modeling. In response to the reviewer’s comment, we conducted supplementary analyses.\nSpecifically, we constructed a monthly RF model using monthly AC values as response variable and\nthe corresponding monthly climatic variables as predictors. This update was made possible by\naccess to high-performance computing resources and the efficient implementation of the ‘ranger’ R\npackage.\nThe monthly RF model achieves an R² of 0.81 and results show that the importance rankings\nof predictor variables derived from the monthly and annual RF models are largely consistent,\nsuggesting that the temporal resolution of input variables has a limited effect on identifying the\ndominant drivers of long-term vegetation resilience. Furthermore, using the monthly RF model,\nwe explored lagged vegetation–climate responses by constructing partial dependence plots for ENSO\nwith 3-, 6-, 9-, 12-, 18-,24-, and 36- month lags (Fig. 3). These analyses reveal asymmetric and\ndelayed vegetation responses to El Niño and La Niña events, consistent with patterns observed in\nthe annual RF model. Specifically, vegetation tends to respond more rapidly and strongly to El Niño,\nwhile the global impacts of extreme La Niña events intensify at longer lags.\nHowever, we note that in the monthly RF models involving numerous lagged predictors, the\nmodel tends to prioritize variables with lower collinearity. This may lead to higher importance being\nassigned to more temporally distant predictors, thereby introducing ambiguity in interpreting the\nresults. To avoid such confounding effects and to simplify interpretation, we ultimately used the\nannually aggregated variables in the main RF analysis presented in the manuscript.\nFig. 3. Influence of explanatory variables on global vegetation resilience from 1981 to 2018\nbased on monthly variables. (a) Variables’ importance scores from the Random Forest (RF) model.\n(b-f) Partial dependence plots (PDPs) show the influence of key predictors from each variable group\non AC: (b) climatic variable, including ENSO, precipitation, temperature, and radiation; (c) Lag effects\nof ENSO at 3-, 6-, 9-, 12-, 18-, 24-, and 36- month intervals; (d) ENSO–climate interactions.\nPositive and negative interaction terms indicate climate anomalies associated with El Niño and La\nNiña phases, respectively; (d) climatic seasonality, represented by intra-annual coefficient of\nvariation (CV) of precipitation, temperature, and radiation (e.g., temperature_seas indicates\ntemperature seasonality); and (f) topographic factors, including DEM, slope, and Height Above\nNearest Drainage (Hand).\nLine 413: Which ENSO indices are referred to here? How were they calculated? Why only\nDecember-February?\nReply: ENSO was quantified using the December–February (DJF) average sea surface temperature\nanomaly (SSTA) over the Niño 3.4 region, derived from the NOAA Extended Reconstructed Sea\nSurface Temperature (ERSST) v5 dataset. We selected the DJF period because it typically\ncorresponds to the mature phase of ENSO events, when the ocean–atmosphere coupling is strongest\nand the teleconnected climate impacts are most pronounced. This approach is widely used in ENSO-\nrelated ecological and climate studies to capture peak ENSO influence. We have clarified the\nexplanation of the predictor variables used in the Random Forest model. These changes have been\nincorporated into the Methods section of the revised manuscript.\nThe revised paper now reads:\n(Page 13, Lines 465-473, Methods)\nIn this research, the response variable was the annual time series of AC, which characterizes\nvegetation resilience. Predictor variables included a combination of climatic and terrain factors.\nClimatic predictors consisted of the annual mean near-surface (2 m) air temperature, annual mean\nof monthly total precipitation, and annual mean of monthly accumulated surface shortwave\nradiation, as well as the annual ENSO index. The annual ENSO index was represented by the\nmean sea surface temperature anomaly (SSTA) over the Niño 3.4 region, calculated as the\naverage over December–February (DJF) using the NOAA ERSST v5 dataset to capture the peak\nphase of ENSO events22. All climatic variables were computed at the pixel level.\nComment 13\nLine 419: Is one year enough to capture all climate legacy effects? Is it necessary to include impact\nof previous years?\nReply: In our initial model design, we focused on short-term lag effects (within one year), as these\nare commonly associated with the delayed vegetation responses to ENSO events via atmospheric\nand oceanic teleconnection processes, which typically operate within 12 months42. Based on this\nrationale, we originally used a one-year lag window and did not incorporate longer-term climate\nlegacies.\nIn response to the reviewer’s suggestion, we have revised our modeling framework to\nincorporate climatic variables from the previous three years, allowing us to better assess whether\nENSO-induced legacy effects influence vegetation resilience. This revision was applied in both the\nannual mean model and the monthly model (Fig 3).\nyear climatic variables remain dominant, variables with a two-year lag contributed significantly\namong the lagged predictors, whereas the importance of the third-year lag was substantially lower.\nThis suggests that multi-year legacy effects exist but tend to weaken over time.\nIn the random forest model built with monthly-resolution variables, the relative importance of\nlagged variables appeared to increase with longer time lags. We interpret this with caution, as it may\nreflect a model preference for predictors with lower collinearity rather than a genuine ecological\nsignal. Indeed, longer-lag variables tend to be less correlated with current-year conditions and may\nbe favored by the model during predictor selection. We have included a note of this potential\nconfounding effect in the revised manuscript.\nOverall, the inclusion of multi-year lags improves the robustness of our analysis and provides\na more comprehensive view of both immediate and delayed vegetation responses to ENSO. These\nmethodological updates and the associated results have now been incorporated into the revised\nMethods and Results sections.\nThe revised paper now reads:\n(Pages 4-5, Lines 151-153, Results)\nWith increasing lag time, the global influence of extreme La Niña events becomes more\npronounced and, by the third year (the lag = 2), surpasses that of El Niño in shaping vegetation\nresilience dynamics (Fig. 2c). This temporal difference likely stems from the contrasting\ndurations of the two phenomena: La Niña often persists into the second year and re-strengthen\nin the winter, whereas El Niño tends to terminates shortly after reaching its peak stage.\n(Page 5, Lines 161-171, Results)\nVariables’ importance scores from the Random Forest (RF) model. (b-f) Partial dependence plots\n(PDPs) show the influence of key predictors from each variable group on AC: (b) climatic variables,\nincluding ENSO, precipitation, temperature, and radiation; (c) Lag effects of ENSO at 1-, 2-, and\n3-year intervals; (d) ENSO–climate interactions. The interaction terms are represented by the\nproduct of ENSO and climate variables. Positive and negative extremes represent high anomaly\nvalues associated with El Niño and La Niña events, respectively, while values near zero indicate\nlow anomalies; (e) climatic seasonality, represented by intra-annual coefficient of variation (CV) of\nprecipitation, temperature, and radiation (e.g., temperature_seas indicates temperature\nseasonality); and (f) topographic factors, including DEM, slope, and Height Above Nearest\nDrainage (Hand).\n(Page 13, Lines 474-477, Methods)\nThe climate factors in the preceding three year were employed to represent the lagged effect, while\nthe interaction terms (temperature × ENSO, precipitation × ENSO, and radiation × ENSO) captured\nhow resilience sensitivity to climate changes varies with ENSO intensity.\nComment 14\nLine 420: I would expect very strong correlations between ENSO and these climatic variables, is\nthat not problematic for the model?\ntemperature, precipitation, and radiation. However, we employed a Random Forest (RF) model,\nwhich is a non-parametric ensemble learning method that is robust to multicollinearity among\npredictors. Unlike traditional regression models, RF does not require independence between\nvariables and is capable of decoupling nonlinear and interacting effects, allowing it to assess the\nrelative importance of each variable even when correlations exist.\nWe have also carefully examined variable importance rankings and partial dependence plots to\nensure interpretability and to avoid misleading conclusions due to collinearity.\nThe revised paper now reads:\n(Page 13, Lines 462-463, Methods)\nRandom forest regression is a robust machine-learning model that can decouple the interactions\nand identify key drivers of change.\nComment 15\nLine 425: Were these pixels randomly chosen or did you take the percentage of different land cover\ntypes into account?\nReply: In the revised manuscript, we have clarified our sampling strategy. We accounted for land\ncover type proportions when selecting both training and validation samples. Specifically, we first\nstratified the global vegetated pixels according to their land cover types and then randomly selected\nsamples from each group in proportion to their global distribution. This approach ensured that the\ntraining and validation datasets were representative of the full diversity of global vegetation types.\nWe have now added a description of this sampling procedure to the Methods section.\nThe revised paper now reads:\n(Page 13, Lines 483-485, Methods)\nWe employed the ‘ranger’ package in R to train the model with 40% of the global vegetated pixels\nand validated it using another 10% of pixels, both randomly selected after stratified allocation\nbased on vegetation type proportions.\nComment 16\nLine 436: How robust is this threshold of -1 or 1 for negative or positive anomalies? Do the results\nchange a lot when adapting this threshold?\nReply: In our study, we used Z scores to extract standardized anomalies of autocorrelation (AC),\n18\nand initially set the threshold at ±2 (previously misreported as ±1) to define significant positive and\nnegative AC anomalies. This threshold is widely adopted in environmental anomaly detection, as it\napproximately captures the top and bottom 5% of values under a normal distribution. It effectively\nfilters out weak fluctuations while retaining areas with substantial responses. Likewise, we used sea\nsurface temperature anomalies (SSTA) of ±1 (previously misreported as ±0.5) as the threshold to\nidentify ENSO events. These corrections have been incorporated in the revised manuscript.\nIn light of your suggestion, we further conducted a sensitivity analysis to assess the robustness of\nthese thresholds. Specifically, we tested alternative AC anomaly thresholds (e.g., ±1, ±1.5, ±2, ±2.5,\n±3) and ENSO intensities (SSTA ±0.5 and ±1.5). We then compared the spatial distribution and pixel\narea statistics of the affected regions (Fig. 4). The results show that although the total area of\nteleconnection-affected regions decreases with more stringent thresholds, the overall spatial patterns\nremain relatively stable. These findings support the appropriateness and robustness of using ±2 for\nAC anomalies and ±1 for SSTA anomalies in the main analysis.\nautocorrelation (AC) anomalies, as determined by chi-square tests.\nThe revised paper now reads:\n(Page 14, Lines 498-501, Methods)\nPositive and negative anomalies in vegetation resilience were defined as standardized AC anomalies\n(Z-scores) exceeding +2 or falling below –2, respectively, corresponding approximately to the upper\nand lower 5% tails of the distribution. Similarly, the definition of El Niño/La Niña events was defined\nbased on the commonly accepted criterion: a monthly SSTA greater than 1 °C or less than -1 °C43,44.\nComment 17\nThe text is not consistent in its terminology: sometimes 1-lag AC and sometimes lag-1 AC, better\nto be consistent.\nReply: We have corrected all the “1-lag AC” to “lag-1 AC” throughout our revised manuscript.\nThe revised paper now reads:\n(Page 1, Lines 20-22, Abstract)\nHere we estimate global present-day (1981–2018) and future (2015–2100) vegetation resilience\nusing a lag-1 autocorrelation analysis of global leaf area index (LAI) time series and investigate its\nteleconnection to ENSO.\n(Page 3, Lines 98-100, Results)\nThis global assessment indicates that 41.56% of vegetated areas have experienced a statistically\nsignificant decline in resilience over the past 40 years, with considerable spatial heterogeneity (Fig.\n1a).\n(Page 4, Lines 133-134, Results)\nFig. 1. Global patterns of historical vegetation resilience changes. (a) lag-1 autocorrelation (AC)\ntrend and (b) variance trend derived from GLASS LAI data from 1981 to 2018.\nComment 18\nFigure 1e: It would be more consistent to use capital letters for delta_ac and delta_lai in the figure\nand the caption.\nconsistency with the rest of the manuscript.\nThe revised paper now reads:\n(Page 4, Lines 132-141, Results)\nFig. 1. Global patterns of historical vegetation resilience changes. (a) lag-1 autocorrelation (AC)\ntrend and (b) variance trend derived from GLASS LAI data from 1981 to 2018. Trends are quantified\nusing the Kendall’s tau statistic. Red regions indicate significant increases in AC or variance, which\nare indicative of declining vegetation resilience. (c) Trend agreement between AC and variance. (d)\nLAI time series trend (1981–2018). (e) Proportion of autocorrelation increase (ΔAC>0) or\ndecrease (ΔAC<0) in greening (ΔLAI>0) and browning (ΔLAI<0) among different vegetation\ntypes. Vegetation types: ENF, Evergreen Needleleaf Forests; EBF, Evergreen Broadleaf Forests;\nDNF, Deciduous Needleleaf Forests; DBF, Deciduous Broadleaf Forests; MF, Mixed Forests;\nCS, Closed Shrublands; OS, Open Shrublands; WS, Woody Savannas.\nFigure 2f: Spelling mistake in x axis label. It is also not clear how to interpret the values on the x\naxis: e.g. what does an ENSO-temperature value of -200 mean?\nReply: We have corrected the spelling error in the x-axis label of Fig. 2f (now revised as Fig. 2d),\nchanging “ENSO-precitation” to “ENSO-precipitation. Regarding the interpretation of the x-axis\nvalues, the variables “ENSO-temperature,” “ENSO-precipitation,” and “ENSO-radiation” represent\ninteraction terms between the ENSO index and local climate variables—specifically, the product of\nthe ENSO index with local temperature (K), precipitation (mm), and solar radiation (W/m²),\nrespectively. As all local climate variables are strictly positive, these interaction terms reflect the\nmodulation of climate anomalies by ENSO phase. Positive values indicate strong anomalies\ntypically associated with El Niño events, whereas negative values are linked to La Niña conditions.\nValues close to zero suggest weak or no ENSO-related anomalies. The x-axis values are\nunstandardized, and thus their magnitudes depend on the scales of the original input variables.\nThe revised paper now reads:\n(Page 4, Lines 143-147, Results)\nENSO primarily influences vegetation resilience through its interactions with climate variables,\nwhich exert a stronger impact than ENSO’s direct effects (Fig 2a). Among the interaction terms,\nreduced vegetation resilience is closely associated with anomalously high temperatures, low\nprecipitation, and enhanced solar radiation during El Niño events, as well as with high-radiation\nand low-precipitation anomalies under La Niña conditions (Fig. 2d).\n(Page 5, Lines 161-171, Results)\nVariables’ importance scores from the Random Forest (RF) model. (b-f) Partial dependence plots\n(PDPs) show the influence of key predictors from each variable group on AC: (b) climatic variables,\nincluding ENSO, precipitation, temperature, and radiation; (c) Lag effects of ENSO at 1-, 2-, and\n3-year intervals; (d) ENSO–climate interactions. The interaction terms are represented by the\nproduct of ENSO and climate variables. Positive and negative extremes represent high anomaly\nvalues associated with El Niño and La Niña events, respectively, while values near zero indicate\nlow anomalies; (e) climatic seasonality, represented by intra-annual coefficient of variation (CV)\nof precipitation, temperature, and radiation (e.g., temperature_seas indicates temperature\nseasonality); and (f) topographic factors, including DEM, slope, and Height Above Nearest\nDrainage (Hand).\nComment 20\nFigure 3: Map a is too small to distinguish between the 7 different colours. Why is (b) in the caption\nlisted after (c-h)? Why are there boxes with coloured stripes included in each of the graphs in (b)?\nReply: We have revised the figure 3a by simplifying the color scheme: only the dominant\nteleconnection types are now displayed — specifically, positive and negative impacts of El Niño\nand La Niña. This modification improves the readability of the map and enhances the visual contrast\n23\namong categories. To improve clarity, we have also revised the caption of Figure 3 to follow the\nlogical visual order (a → b → c–h). Additionally, we updated the panel title for Figure 3b to explain\nthe purpose of the colored inset boxes. Specifically, the inset in the top-right corner magnifies the\ndifference between regional average AC trends (represented by the dashed lines), which would\notherwise be difficult to distinguish due to their subtle differences. The asterisks embedded in the\nsubplots indicate whether the differences between positively affected, negatively affected, and\nunaffected areas are statistically significant, as determined by both the Whitney U-test and Student’s\nt-test. The significance levels are denoted as p < 0.05 (*), p < 0.01 (**), and p < 0.001 (***). A\ncorresponding clarification has also been added to the revised title of Figure 3b.\nThe revised paper now reads:\n(Pages 7-8, Lines 236-247, Results)\nFig 3. ENSO impact on global vegetation resilience. (a) Direct impact of ENSO on resilience. (b)\nRegional differences in AC trends between affected and unaffected regions during 1981–2018. The\ndashed line represents the regional average. The inset box in the top right corner magnifies the\ndifference in regional averages. Asterisks indicate statistical significance between groups based\non both the Whitney U-test and Student’s t-test: p < 0.05(*); p < 0.01(**); p < 0.001(***). (c-h)\nImpacts of ENSO-climate interaction on vegetation resilience: (c-e) Impacts of ENSO-induced\nanomalies in temperature, solar radiation, and precipitation on vegetation resilience, respectively.\n(f) and (g) Combined impact of these climate anomalies during El Niño and La Niña on vegetation\nresilience, respectively; (h) Overall impact of ENSO-induced climate anomalies on vegetation\nresilience. Plus (+) and minus (-) signs indicate positive and negative impacts, respectively, with\nnumbers indicating impact intensity, calculated as the sum of individual ENSO-climate anomaly\nimpacts.\nLines 195-201 and figure 3b: I’m confused with the explanation on which differences are significant\nand which ones are not, because all three graphs in figure 3b look the same to me.\nReply: In Figure 3b, differences in AC trends among regions were assessed using both the Wilcoxon\nrank-sum test (Mann–Whitney U-test) and Student’s t-test. We acknowledge that the original figure\ndid not clearly depict these differences. To improve clarity, we have added an inset in the top-right\ncorner of each subplot that magnifies the regional average AC trends (dashed lines). Statistical\nsignificance indicators (*, **, ***) have been included to denote differences between regional pairs,\ncorresponding to p < 0.05, p < 0.01, and p < 0.001, respectively. An explanatory note regarding this\nhas also been added to the revised panel title of Figure 3b.\nFor example, in the ENSO impact area panel, the AC trend in positively affected regions (green box)\ndoes not significantly differ from that in unaffected regions (blue box), while both differ\nsignificantly from negatively affected regions (red box), as indicated by ***. Notably, negatively\naffected regions exhibit the highest mean AC trend, implying a greater likelihood of resilience\ndecline, which supports the original statement that “the historical decline in resilience trend was\nsignificantly more pronounced in negatively affected regions compared to positively affected and\nnon-teleconnected regions (p < 0.05).”\nThe revised paper now reads:\n(Pages 7-8, Lines 236-247, Results)\nFig 3. ENSO impact on global vegetation resilience. (a) Direct impact of ENSO on resilience. (b)\ndashed line represents the regional average. The inset box in the top right corner magnifies the\ndifference in regional averages. Asterisks indicate statistical significance between groups based\non both the Whitney U-test and Student’s t-test: p < 0.05(*); p < 0.01(**); p < 0.001(***). (c-h)\nImpacts of ENSO-climate interaction on vegetation resilience: (c-e) Impacts of ENSO-induced\nanomalies in temperature, solar radiation, and precipitation on vegetation resilience, respectively.\n(f) and (g) Combined impact of these climate anomalies during El Niño and La Niña on vegetation\nresilience, respectively; (h) Overall impact of ENSO-induced climate anomalies on vegetation\nresilience. Plus (+) and minus (-) signs indicate positive and negative impacts, respectively, with\nnumbers indicating impact intensity, calculated as the sum of individual ENSO-climate anomaly\nimpacts.\nComment 22\nLine 2019: How was the projected vegetation resilience calculated from the CMIP6 data?\nReply: Thank you for the request for further clarification. This study utilized the GLASS LAI\ndataset alongside LAI data derived from three CMIP6 scenarios to estimate both historical and\nfuture vegetation resilience. The CMIP6 LAI data were processed consistently with the historical\ndata by first applying a LOESS decomposition to extract residuals series from the LAI time series,\nfollowed by calculation of the lag-1 autocorrelation coefficient within a 5-year sliding window.\nVegetation pixels were selected based on land use data to estimate vegetation resilience. These\ndetails have been added to the Methods section under Vegetation resilience calculation.\nThe revised paper now reads:\n(Page 12, Lines 427-444, Methods)\nIn this study, the GLASS LAI dataset and CMIP6-derived LAI data under three scenarios were\nused to calculate historical and future vegetation resilience. However, prior to measuring the\nresilience metrics, it is necessary to remove the impact of long-term trends of LAI time series. STL\n(Seasonal-Trend decomposition using LOESS) is a widely used time series decomposition method.\nIt gradually extracts trend and seasonal components through an iterative smoothing process and\nfinally obtains the residual45. We employed STL decomposition through stl() function in ‘stats’\npackage in R for each grid. This method robustly separate the LAI timeseries into three components:\na long-term trend, a fixed seasonal cycle trend (by using ‘periodic’ option in s.window parameter),\nand a residual component46,47. Then, the residual components of the LAI timeseries were utilized to\ncalculate the lag-1 AC and variance. The ar.ols() and var() function in ‘stats’ package in R were\nemployed to calculate the lag-1 autocorrelation and variance on the sliding windows of 5 years. To\ninvestigate the effect of window size on the resilience trends, we also calculated AC with a sliding\nwindow of 10 years (Extended Data Fig. 11). Note that only the regions where AC > 0 were involved\nin this study. For the LAI data preprocessing, we did not fill the data gap because the relationship\nbetween the AC of the gap-free and gapped time series has been proven to be very close to being\none-to-one 48. We also excluded the grids of croplands (where human management may strongly\ndisturb the response of crops to ENSO) and non-vegetation types, so that our analysis focuses on\nthe vegetation of non-crop ecosystems.\nFigure 4b: the maps are too small to distinguish between the 7 different colours.\nReply: Figure 4b illustrates the direct impact of ENSO on vegetation resilience under the three\nscenarios. We have revised the figure by simplifying the color scheme: only the major impact types\nare now shown — namely, positive and negative impacts of El Niño and La Niña. This improves\nthe clarity of the map and enhances the visual distinction between the categories.\nThe revised paper now reads:\n(Page 8, Lines 268-274, Results)\nFig. 4 Projected vegetation resilience trend and ENSO’s direct and interactive impact on\nvegetation resilience under future scenarios (2015–2100). (a) Resilience trend under SSP 126, SSP\n245, and SSP 370, respectively. (b) Direct impact of ENSO on vegetation resilience in the three\nscenarios. (c) Impact of ENSO-climate interaction on vegetation resilience in three scenarios. Plus\n(+) and minus (-) signs indicate positive and negative impacts, respectively, with numbers indicating\nimpact intensity, calculated as the sum of individual ENSO-climate interaction impacts.\nComment 24\nThe background information about ENSO in the second paragraph belongs in the introduction,\nrather than in the discussion.\nReply: In the revised manuscript, we have removed the background information on ENSO from the\nhave been selectively integrated into the Introduction, where they provide appropriate context for\nthe study.\nThe revised paper now reads:\n(Page 2, Lines 52-63, Introduction)\nEl Niño typically raises temperature and solar radiation in tropical and subtropical regions, yet its\nimpact on precipitation varies, causing wetter conditions in some areas and drought in others,\npotentially exacerbating or alleviating vegetation water stress27-29. For instance, El Niño-driven\nshifts in the Intertropical Convergence Zone (ITCZ) modify regional precipitation patterns,\ndecreasing precipitation in West and northern Central Africa, while increasing it in East and\nsouthern Central Africa30, thus creating a contrasting north-south gradient in vegetation\ndynamics. Conversely, La Niña generally exhibits the opposite climatic signals, particularly over\ntropical land areas where temperature and precipitation patterns respond in antiphase, thereby\npartially offsetting the effects of El Niño31,32. However, this compensatory mechanism is modulated\nby variation in the location of equatorial warming and by nonlinear ocean–atmosphere interactions,\nleading to spatial asymmetries in the timing and intensity of vegetation response to ENSO phases.\nComment 25\nThe discussion does not include any discussions or conclusions about the results, only that the\nresearch did not include anthropogenic effects and that there was a lot of uncertainty in the data.\nReply: In the revised manuscript, we have substantially expanded this section to provide a more in-\ndepth interpretation of our results. We also drew clearer comparisons with previous studies, thereby\nsituating our findings within the broader literature. To enhance mechanistic understanding, we\nincluded further explanation of the spatial heterogeneity in ENSO sensitivity. Moreover, we\nemphasized the broader implications of our work, highlighting the importance of incorporating\nENSO variability into ecosystem management frameworks—particularly for tropical and\nsubtropical regions that are especially vulnerable to climate anomalies.\nThe revised paper now reads:\n(Pages 9-10 Lines 276-327, Discussion)\nThis research provides a comprehensive analysis of the climatic drivers of vegetation resilience,\nuniquely incorporating the remote impacts of global climate oscillations such as ENSO. Using a\nRF model, we disentangled the effects of local climate variables (temperature, precipitation, and\nsolar radiation) from those of ENSO and their interactions, revealing how ENSO shapes\nresilience patterns across global, regional, and biome-specific scales. Projections further indicate\nevolving impacts of ENSO under future climate scenarios. One key finding is that 51% of global\nvegetation over the past four decades has been influenced by either direct ENSO effects or ENSO–\n28\nclimate interactions. These impacts are both asymmetric and delayed, with El Niño and La Niña\nevents exerting opposing yet destabilizing effects on vegetation resilience, contributing to long-\nterm declines in resilience, most notably in the Amazon, Southeast Asia, Central Africa, southern\nNorth America, and Australia (Fig. 3).\nThese regional hotspots of ENSO sensitivity can be attributed to a complex interplay of climatic\ndrivers, ecological constraints, and human activities. ENSO-driven anomalies in heat and\nmoisture show pronounced spatial variation, exposing ecosystems to different types and\nintensities of extreme events (e.g., heatwaves, droughts, floods) that impair photosynthesis, water\nregulation, and root function1,2. These effects are most severe in ecological hotspots such as\ntropical forest, savanna and grassland. In addition, our model highlights topographic features—\nparticularly high elevation, steep slopes, and well drained regions—as key predictors of low\nresilience and strong ENSO responses (Fig. 2f). These landscapes often support vegetation with\nshallow roots, limited soil depth, and poor water and nutrient retention, resulting in reduced\nbuffering capacity3. In parallel, climatic seasonality drives divergent vegetation adaptive\nstrategies, which in turn influence regional sensitivity to ENSO. While low-latitude regions with\nweak seasonality exhibit strong resilience shifts globally (Fig. 2e), ecosystems in areas with high\nintra-annual precipitation variability—such as northern Asian shrublands, Siberian tundra,\neastern Amazon rainforest, and Sahelian grasslands—also show heightened sensitivity compared\nto adjacent regions with similar rainfall and vegetation types (Fig. 3a). This supports prior\nfindings1,4,5 that stable water availability enhances post-disturbance recovery.\nEqually important, anthropogenic disturbances exacerbate these patterns6. Highly fragmented\ntropical and temperate forests, such as those in the Amazon, Central Africa, and parts of North\nAmerica and Southeast Asia, exhibit heightened sensitivity to ENSO (Fig. 3a), likely due to\nintensified edge effects that increase local temperatures and reduce humidity, amplifying climate\nstressors7-9. Intensive agriculture drives land degradation, biodiversity loss, and hydrological\ndisruption, further destabilizing adjacent ecosystems10,11. For instance, we found that mixed\nforests along the Eurasian agro-ecological transition zone show more negative responses to La\nNiña when located near croplands (Fig. 3a), suggesting landscape configuration modulates\nENSO sensitivity—a mechanism that merits further empirical investigation.\nOur findings reinforce previous evidence that vegetation in climatically unstable environments—\ncharacterized by high variability in light, temperature, and moisture—tends to exhibit lower\nresilience. The identified hotspots of ENSO sensitivity are broadly consistent with earlier studies\non vegetation–climate responses12. For example, we found that vegetation resilience is notably\nlower in the southeastern Amazon, consistent with earlier studies that associate this decline with\nrecurring droughts and floods13,14. Similarly, reduced resilience along the Sahelian margins\nsupports earlier findings on the role of increasing water limitation15. ENSO-related effects\nobserved in Southeast Asia, Australia, Central Africa, and Central America also reflect patterns\ndocumented in prior research16-18. As climate change intensifies, such instability is expected to\nincrease, threatening ecosystem stability across large regions4. It is therefore critical to integrate\nENSO forecasts into ecosystem restoration and climate-informed conservation planning,\nespecially for tropical and subtropical forests. These ecosystems, given their high biodiversity,\ncarbon storage potential, and sensitivity to both climate and human disturbance, should be\nprioritized in global strategies for resilience and adaptation. In contrast, some northern forests,\nsuch as mixed forests in Northern Europe, show increasing resilience and minimal ENSO\nsensitivity (Fig. 1a and Fig. 3a), highlighting the need for regionally differentiated adaptation\nstrategies.\nWe sincerely appreciate your thoughtful and constructive comments. We are grateful for your\nrecognition of our study as well as for your valuable suggestions to improve the experimental design\nand deepen the discussion. To address your concerns in detail, we have organized our responses\npoint by point. We hope the revisions we have made fully address the issues you raised.\nComment 1\nClarify the analytical framework throughout the Results section:\nThe Results section of the manuscript currently lacks sufficient context on how the analyses were\nperformed. Please clearly state the type of analysis in each subsection. This will help readers follow\nthe sequence and logic of the findings. Please give more context in the results, including the type of\nanalysis that you performed for each section of the paper. This will help the readers follow what you\nhave done. For example:\nFigure 1a shows the spatial distribution of 1-lag autocorrelation (AC) trend. Blue values are negative\nand red values are positive. What does it mean? You mentioned in the caption “An increase in AC\nor variance reflects a decline in vegetation resilience”, so I assume that red is a decline in vegetation\nresilience?\nReply: We have added brief descriptions of the analytical approaches within each subsection of the\nResults to clarify the methodological framework and enhance the logical flow. This additional\ncontext guides readers through the sequence of analyses and the interpretation of key findings. In\naddition, regarding the reviewer’s question on Figure 1a, we have added clarification in the\nIntroduction explicitly stating the relationship between increases in AC or variance and declining\nresilience. Moreover, we now explicitly indicate in the figure caption that regions with increasing\nAC (shown in red) correspond to a loss of resilience.\nThe revised paper now reads:\n(Page 2, Lines 80-82, Introduction)\nThe critical slowing down theory (CSD) provides a framework for detecting ecosystem instability\nby suggesting that declining recovery rates precede state transitions, which can be identified by\nstatistical signals like increased variance and higher first-order autocorrelation.\n(Page 3, Lines 96-100, Results)\nWe quantified vegetation resilience at the pixel level using lag-1 autocorrelation (AC) derived\nfrom 5-year moving windows of GLASS LAI data (1981–2018), and assessed its long-term trend\n(Methods). This global assessment indicates that 41.56% of vegetated areas have experienced a\nstatistically significant decline in resilience over the past 40 years, with considerable spatial\nheterogeneity (Fig. 1a).\n30\n(Page 4, Lines 132-141, Results)\nFig. 1. Global patterns of historical vegetation resilience changes. (a) lag-1 autocorrelation (AC)\ntrend and (b) variance trend derived from GLASS LAI data from 1981 to 2018. Trends are quantified\nusing the Kendall’s tau statistic. Red regions indicate significant increases in AC or variance,\nwhich are indicative of declining vegetation resilience. (c) Trend agreement between AC and\nvariance. (d) LAI time series trend (1981–2018). (e) Proportion of autocorrelation increase (Δ\nAC>0) or decrease (ΔAC<0) in greening (ΔLAI>0) and browning (ΔLAI<0) among different\nvegetation types. Vegetation types: ENF, Evergreen Needleleaf Forests; EBF, Evergreen Broadleaf\nForests; DNF, Deciduous Needleleaf Forests; DBF, Deciduous Broadleaf Forests; MF, Mixed\nForests; CS, Closed Shrublands; OS, Open Shrublands; WS, Woody Savannas.\n(Page 3, Lines 115-116, Results)\nWe used a random forest regression framework at the global scale to attribute temporal changes\nin vegetation resilience (AC) to ENSO variability and other environmental factors (Methods).\n(Page 6, Lines 175-179, Results)\nTo reveal the global spatial pattern of ENSO’s direct impact on vegetation resilience, we used the\nPearson Chi-square test to assess their teleconnection over a long historical period (Methods).\nsignificant relationships between vegetation resilience anomalies and the ENSO events (P<0.01),\nwith resilience demonstrating higher sensitivity to El Niño than La Niña events (Fig. 3a).\n(Page 6, Lines 198-200, Results)\nRF inference (Fig. 2c-e) suggests that ENSO-induced climate anomalies indirectly trigger resilience\nchanges. By overlapping resilience-climate (Extended Data Fig. 6) and ENSO-climate (Extended\nData Fig. 7) linkages, we reveal spatially heterogeneous coupling patterns,\n(Page 8, Lines 249-251, Results)\nLeveraging simulated LAI and climate projections from CMIP6, we evaluated future trajectories\nof vegetation resilience and examined their links to climate change across three emission\nscenarios (SSP126, SSP245, and SSP370) (Methods). Our results indicate that 33-39% of\nvegetated regions are anticipated to undergo a further decline in resilience across three future\nclimate scenarios, with an average decline trend of 0.12 in these affected areas.\nComment 2\nLine 91 - “The most significant declines in resilience 92 are notably observed in northern North\nAmerica, the Amazon, northern Europe, eastern Australia, 93 and northern and Southeast Asia ”.\nWhat are these results based on? How did you compare these areas to get to the conclusion that\nthese areas are particularly affected?\nReply: The identification of regions with the most significant declines in resilience is based on the\nspatial distribution of trends in lag-1 autocorrelation (AC) from 1981 to 2018 (Fig. 1a). We used\nKendall’s tau to assess AC trends for each pixel, and regions with the strongest positive trends (i.e.,\ngreatest increases in AC) were interpreted as experiencing the most pronounced resilience loss. We\nrevised the text to clarify that this conclusion is drawn from the spatial patterns of statistically\nsignificant positive AC trends, which indicate declining ecosystem resilience.\nThe revised paper now reads:\n(Page 3 Lines 100-102, Results)\nThe most significant declines in resilience are observed in northern North America, the Amazon,\nnorthern Europe, eastern Australia, and northern and Southeast Asia, which exhibit the greatest\nincrease in AC (Fig. 1a).\nFigure 1c- “(c) Trend agreement between two indicators. ”. Please state what indicators these are.\nReply: We have revised the figure title to “Trend agreement between AC and variance” to clearly\nindicate the two indicators being compared. This comparison was included to assess the consistency\nbetween different resilience metrics and to ensure the robustness of the AC-based indicator.\nThe revised paper now reads:\n(Page 4 Lines 133-141, Results)\nFig. 1. Global patterns of historical vegetation resilience changes. (a) lag-1 autocorrelation (AC)\ntrend and (b) variance trend derived from GLASS LAI data from 1981 to 2018. Trends are quantified\nusing the Kendall’s tau statistic. Red regions indicate significant increases in AC or variance, which\nare indicative of declining vegetation resilience. (c) Trend agreement between AC and variance.\n(d) LAI time series trend (1981–2018). (e) Proportion of autocorrelation increase (ΔAC>0) or\ndecrease (ΔAC<0) in greening (ΔLAI>0) and browning (ΔLAI<0) among different vegetation\ntypes. Vegetation types: ENF, Evergreen Needleleaf Forests; EBF, Evergreen Broadleaf Forests;\nDNF, Deciduous Needleleaf Forests; DBF, Deciduous Broadleaf Forests; MF, Mixed Forests;\nCS, Closed Shrublands; OS, Open Shrublands; WS, Woody Savannas.\nComment 4\nLine 97 - “Regions such as the Amazon, south-central Africa, Southeast Asia, and western Oceania\nhave experienced a widespread increase in greenness—measured by rising LAI— over the last four\ndecades, yet these areas have concurrently witnessed significant declines in resilience ”. Please\nalso explain how you got these results in the text.\nReply: We quantified vegetation greenness using the LAI time series and assessed its long-term\ntrend using the non-parametric Kendall trend test (shown in Fig. 1d). Similarly, vegetation resilience\nwas assessed based on the trend of lag-1 autocorrelation (AC) over time (Fig. 1a). By spatially\noverlaying the two trend maps, we identified regions where increasing greenness coincided with\ndeclining resilience, such as the Amazon, south-central Africa, Southeast Asia, and western Oceania.\nWe have revised the manuscript to clarify this analytical procedure in Methods, and we have also\nadded a supplementary figure to illustrate the spatial correspondence between greenness and\nresilience trends.\nThe revised paper now reads:\n(Page 13, Lines 456-458, Methods)\nLong-term LAI trend was quantified using Kendall’s tau at p < 0.05 to represent vegetation\n33\ngreenness change (Fig. 1d). We then overlaid the AC and LAI trend maps to assess the spatial\ncoupling between greenness and resilience.\nComment 5\nLine 219 - “Our results indicate” also really requires explaining how you did the analysis.\nReply: We have clarified the methodology used to derive this result. Specifically, the projected\ntrends of AC under future climate scenarios were calculated using Kendall’s tau coefficient,\nconsistent with our approach to historical trend analysis. This information has now been explicitly\nreferenced in the revised Methods section titled “Vegetation Resilience and Greenness Trend”.\nThe revised paper now reads:\n(Page 8, Lines 249-253, Results)\nLeveraging simulated LAI and climate projections from CMIP6, we evaluated future trajectories\nof vegetation resilience and examined their links to climate change across three emission\nscenarios (SSP126, SSP245, and SSP370) (Methods). Our results indicate that 33-39% of\nvegetated regions are anticipated to undergo a further decline in resilience across three future\nclimate scenarios, with an average decline trend of 0.12 in these affected areas.\n(Page 13, Lines 447-455, Methods)\nTo assess the historical vegetation resilience trends, we calculated the Kendall tau value for the AC\nand variance time series using a 5-year sliding window at each grid from 1981 to 2018. This rank-\nbased correlation coefficient is independent of data distribution, allowing for consistent trend\ncomparisons across regions. Using the cor.test() function with the Kendall method from the R\npackage ‘stats’, we estimated resilience trends with a 95% confidence level. The trend directions of\nvariance and AC were largely consistent across most regions (Fig. 1c). Additionally, we calculated\nthe AC trend using a 10-year window, which indicated a spatial pattern qualitatively similar to that\nof the 5-year window (Fig. 1a and Extended Data Fig. 11). AC trends under future scenarios were\nalso calculated (Fig. 4).\nComment 6\nIt is unclear what form of temperature you used (e.g., mean annual temperature, maximum\ntemperature). Please specify this, and similarly for precipitation and radiation.\nReply: We have now specified the exact forms of the climate variables used in our analysis.\nSpecifically, we used annual mean near-surface (2 m) air temperature, annual mean of monthly total\nprecipitation, and annual mean of monthly accumulated surface shortwave radiation. This\ninformation has been added to the revised Methods section for clarity.\n(Page 13, Lines 467-470, Methods)\nClimatic predictors consisted of the annual mean near-surface (2 m) air temperature, annual\nmean of monthly total precipitation, and annual mean of monthly accumulated surface shortwave\nradiation, as well as the annual ENSO index.\nComment 7\nHave you considered climate seasonality (e.g., coefficient of variation or seasonality index of\nprecipitation/temperature)? Seasonality could influence vegetation resilience by selecting for\nspecies with different adaptive strategies.\nReply: We agree that intra-annual variation in temperature and precipitation may influence\nvegetation resilience by selecting for species with different adaptive traits. In this study, our primary\nobjective was to isolate the effect of interannual climate variability—particularly ENSO-related\nanomalies—on vegetation resilience. However, we recognize that spatial differences in climate\nseasonality may shape vegetation composition and function, which could in turn confound the\nobserved ENSO–resilience relationship. To account for this potential background influence, we\ncalculated seasonality indices (coefficient of variation across months) for precipitation, radiation,\nand temperature, and used their multi-year means to represent regional seasonality amplitude. These\nvariables were included to control for the effects of climate-driven biome distribution. The additions\nare now described in the Methods section under Random forest model.\nThe revised paper now reads:\n(Page 3, Lines 116-121, Results)\nThe analysis reveals that local controlling variables—namely topography and climate\nseasonality—play a more dominant role in shaping interannual variation in vegetation resilience\nthan ENSO teleconnections (Fig 2a; model R2 = 0.79). In particular, seasonal climate regimes\ninfluence resilience by selecting species with distinct adaptive traits. Ecosystems subject to lower\nintra-annual climate variability are more prone to resilience loss (Fig. 2e).\n(Page 5, Lines 161-171, Results)\nVariables’ importance scores from the Random Forest (RF) model. (b-f) Partial dependence plots\n(PDPs) show the influence of key predictors from each variable group on AC: (b) climatic variables,\nincluding ENSO, precipitation, temperature, and radiation; (c) Lag effects of ENSO at 1-, 2-, and\n3-year intervals; (d) ENSO–climate interactions. The interaction terms are represented by the\nproduct of ENSO and climate variables. Positive and negative extremes represent high anomaly\nvalues associated with El Niño and La Niña events, respectively, while values near zero indicate\nlow anomalies; (e) climatic seasonality, represented by intra-annual coefficient of variation (CV)\nof precipitation, temperature, and radiation (e.g., temperature_seas indicates temperature\nseasonality); and (f) topographic factors, including DEM, slope, and Height Above Nearest\nDrainage (Hand).\n(Page 13, Lines 479-483, Methods)\nIn addition, intra-annual variations in temperature and precipitation may influence vegetation\nresilience by favoring species with distinct adaptive traits. To account for potential differences in\nresilience due to climatic seasonality, we included multi-year mean seasonality indices (monthly\ncoefficient of variation) of precipitation, radiation, and temperature as covariates in the model.\nFigure 2a highlights slope as an important predictor, but this result is not discussed in the main text.\nPlease include a paragraph interpreting this result—how might slope influence resilience (e.g., via\nsoil depth, drainage, erosion)?\nReply: We agree that slope and elevation are ecologically meaningful terrain metrics that may\nindirectly modulate vegetation resilience. In our results, topographic variables were indeed\nidentified as important predictors in the random forest model (Fig. 2a), and we observed that regions\ncharacterized by high elevation, steep slopes, or well-drained landscapes are especially vulnerable\nto resilience decline (Fig. 2f).\nHowever, as this study focuses on climate-driven mechanisms and relies on relatively coarse-\nresolution data, the primary role of topographic variables (slope, elevation, and HAND) in our\nanalysis was to control for background terrain heterogeneity and to enhance model predictive\nperformance, rather than to explore specific terrain-related ecological processes such as erosion,\nhydrological redistribution, or soil development. We have now clarified this rationale in the revised\nmanuscript to avoid potential misunderstanding.\nWe also agree that future work at finer spatial scales, using more detailed hydrological and\ngeomorphological data, would be valuable for investigating these terrain–ecology interactions in\ngreater depth.\nThe revised paper now reads:\n(Page 3, Lines 116-123, Results)\nThe analysis reveals that local controlling variables—namely topography and climate\nseasonality—play a more dominant role in shaping interannual variation in vegetation resilience\nthan ENSO teleconnections (Fig 2a; model R2 = 0.79). In particular, seasonal climate regimes\ninfluence resilience by selecting species with distinct adaptive traits. Ecosystems subject to lower\nintra-annual climate variability are more prone to resilience loss (Fig. 2e). Terrain constraints\nfurther highlight that regions characterized by high elevation, steep slopes or well-drained\nlandscapes are especially vulnerable to declines in vegetation resilience (Fig. 2d).\n(Page 9, Lines 292-295, Discussion)\nIn addition, our model highlights topographic features—particularly high elevation, steep slopes,\nand well drained regions—as key predictors of low resilience and strong ENSO responses (Fig.\n2f). These landscapes often support vegetation with shallow roots, limited soil depth, and poor\nwater and nutrient retention, resulting in reduced buffering capacity3. (Page 13, Lines 477-479, Methods)\nDEM, slope, and HAND were included in the model as the terrain predictors to control potential\ninterference with terrain background heterogeneity.\nConsider including terrain-based indicators such as Topographic Position Index (TPI) or Height\nAbove Nearest Drainage (HAND), which can better reflect groundwater access and hydrological\nbuffering. This could provide a mechanistic link between terrain, water availability, and resilience.\nFor reference, see: https://www.nature.com/articles/s41586-024-07568-w\nReply：We have incorporated the Height Above Nearest Drainage (HAND) as a terrain-based\nindicator to better capture groundwater accessibility and hydrological buffering capacity. The\nHAND data have been added to the Methods section under Landcover and terrain data and Random\nForest model. Corresponding modifications have also been made in the Results section to reflect its\ninfluence. This addition helps to strengthen the mechanistic link between terrain, water availability,\nand ecosystem resilience, as recommended. The suggested reference\n(https://www.nature.com/articles/s41586-024-07568-w) has also been cited accordingly.\nThe revised paper now reads:\n(Page 3, Lines 116-123, Results)\nThe analysis reveals that local controlling variables—namely topography and climate\nseasonality—play a more dominant role in shaping interannual variation in vegetation resilience\nthan ENSO teleconnections (Fig 2a; model R2 = 0.79). In particular, seasonal climate regimes\ninfluence resilience by selecting species with distinct adaptive traits. Ecosystems subject to lower\nintra-annual climate variability are more prone to resilience loss (Fig. 2e). Terrain constraints\nfurther highlight that regions characterized by high elevation, steep slopes or well-drained\nlandscapes are especially vulnerable to declines in vegetation resilience (Fig. 2d).\n38\n(Page 5, Lines 161-171, Results)\nFig. 2. Influence of explanatory variables on global vegetation resilience from 1981 to 2018. (a)\nVariables’ importance scores from the Random Forest (RF) model. (b-f) Partial dependence plots\n(PDPs) show the influence of key predictors from each variable group on AC: (b) climatic variables,\nincluding ENSO, precipitation, temperature, and radiation; (c) Lag effects of ENSO at 1-, 2-, and\n3-year intervals; (d) ENSO–climate interactions. The interaction terms are represented by the\nproduct of ENSO and climate variables. Positive and negative extremes represent high anomaly\nvalues associated with El Niño and La Niña events, respectively, while values near zero indicate\nlow anomalies; (e) climatic seasonality, represented by intra-annual coefficient of variation (CV)\nof precipitation, temperature, and radiation (e.g., temperature_seas indicates temperature\nseasonality); and (f) topographic factors, including DEM, slope, and Height Above Nearest\nDrainage (Hand).\n(Page 11, Lines 388-391, Methods)\nThe Height Above Nearest Drainage (HAND) data, obtained from Google Earth Engine\n(https://gee-community-catalog.org/projects/hand/), were derived from SRTM elevation data and\nprovide a more accurate representation of groundwater accessibility and hydrological buffering\ncapacity49.\n39\n(Page 13, Lines 477-479, Methods)\nDEM, slope, and HAND were included in the model as the terrain predictors to control potential\ninterference with terrain background heterogeneity.\nComment 10\nAddress the role of disturbance more directly:\nOne concern is the inclusion of both disturbed and undisturbed vegetation in the same analysis.\nThese systems may have fundamentally different dynamics, especially in their response to ENSO\nevents. I recommend conducting a separate analysis for intact forests, using available datasets such\nas the GLAD Intact Forest Landscapes map: https://glad.umd.edu/dataset/intact-forest/overview .\nThis could reveal whether observed resilience declines are driven more by land-use history or by\nclimate sensitivity, and it would strengthen the paper’s conclusions.\nReply:\nWe fully agree that controlling for potential human influences is essential for accurately assessing\necosystem resilience to climate variability. In response, we conducted additional analyses focused\nexclusively on intact forest ecosystems, using the GLAD Intact Forest Landscapes dataset (as now\ndescribed in the Methods and Extended Data). While land-use history was not explicitly included\nas a covariate in our modeling framework, these ecosystems are minimally affected by land\nmanagement practices and thus offer a clearer signal of climate-driven resilience dynamic.\nAs reported in the Results Fig. 2 and Extended Data Fig. 3, the vegetation responses to ENSO in\nthese undisturbed areas are broadly consistent with the global-scale patterns. This indicates that\nalthough anthropogenic disturbances may modulate local resilience responses, they do not\nfundamentally alter the underlying climate–vegetation coupling associated with ENSO. These\nfindings support the robustness of our main conclusions and underscore the dominant role of\nclimatic drivers over land-use history in shaping large-scale vegetation resilience.\nThe revised paper now reads:\n(Page 5, Lines 156-160, Results)\nAn independent analysis targeting intact forest regions (Methods) indicates that after excluding\nareas affected by human activities, vegetation responses to ENSO remain broadly consistent with\nthe global-scale patterns (Extended Data Fig. 3), which imply that anthropogenic disturbances\nmay modulate but do not fundamentally alter the underlying climate–vegetation coupling driven\nby ENSO.\n(Page 11, Lines 382-384, Methods)\nWe used GLAD Intact Forest Landscapes dataset for the year 2000\n(https://intactforests.org/data.ifl.html), which provide a global map of forest without significant\ndisturbance.\n(Pages 13-14, Lines 489-491, Methods)\nTo examine whether human activities influence the observed climate–resilience coupling, we used\nintact forest ecosystems—largely free from direct land-use impacts—as a control group and\ntrained an additional RF model using data exclusively from these areas (Extended Data Fig. 3).\n(Page 20, Lines 754-765, Extended Data)\nExtended Data. 3. Influence of explanatory variables on intact forest vegetation resilience from\n1981 to 2018. (a) Variables’ importance scores from the Random Forest (RF) model. (b-f) Partial\ndependence plots (PDPs) show the influence of key predictors from each variable group on AC: (b)\nclimatic variables, including ENSO, precipitation, temperature, and radiation; (c) Lag effects of\nENSO at 1-, 2-, and 3-year intervals; (d) ENSO–climate interactions. The interaction terms are\nrepresented by the product of ENSO and climate variables. Positive and negative extremes represent\nhigh anomaly values associated with El Niño and La Niña events, respectively, while values near\nzero indicate low anomalies; (e) climatic seasonality, represented by intra-annual coefficient of\nvariation (CV) of precipitation, temperature, and radiation (e.g., temperature_seas indicates\ntemperature seasonality); and (f) topographic factors, including DEM, slope, and Height Above\nNearest Drainage (Hand).\nExpand the discussion of mechanisms:\nThe discussion would benefit from a more thorough treatment of why some regions are more\nsensitive to El Niño/La Niña than others. Consider integrating more discussion on potential drivers\nsuch as species composition and functional traits, soil water storage capacity, historical exposure to\nclimate variability and human impacts and fragmentation\nReply: In the revised version, we have expanded the discussion to explore possible mechanisms\nunderlying regional differences in vegetation sensitivity to ENSO. Specifically, we now highlight\nseveral potential drivers, including (1) varying exposure to extreme climatic anomalies, (2)\ndifferences in topography-related soil and water storage capacity, (3) climate seasonality that shapes\nregion-specific plant adaptation strategies, and (4) anthropogenic disturbances such as\ndeforestation-induced fragmentation and pollution from nearby croplands. These additions aim to\nprovide a more comprehensive understanding of the processes shaping the observed spatial\nheterogeneity in vegetation response.\nThe revised paper now reads:\n(Pages 9-10, Lines 276-327, Discussion)\nThis research provides a comprehensive analysis of the climatic drivers of vegetation resilience,\nuniquely incorporating the remote impacts of global climate oscillations such as ENSO. Using a\nRF model, we disentangled the effects of local climate variables (temperature, precipitation, and\nsolar radiation) from those of ENSO and their interactions, revealing how ENSO shapes\nresilience patterns across global, regional, and biome-specific scales. Projections further indicate\nevolving impacts of ENSO under future climate scenarios. One key finding is that 51% of global\nvegetation over the past four decades has been influenced by either direct ENSO effects or ENSO–\nclimate interactions. These impacts are both asymmetric and delayed, with El Niño and La Niña\nevents exerting opposing yet destabilizing effects on vegetation resilience, contributing to long-\nterm declines in resilience, most notably in the Amazon, Southeast Asia, Central Africa, southern\nNorth America, and Australia (Fig. 3).\nThese regional hotspots of ENSO sensitivity can be attributed to a complex interplay of climatic\ndrivers, ecological constraints, and human activities. ENSO-driven anomalies in heat and\nmoisture show pronounced spatial variation, exposing ecosystems to different types and\nintensities of extreme events (e.g., heatwaves, droughts, floods) that impair photosynthesis, water\nregulation, and root function1,2. These effects are most severe in ecological hotspots such as\ntropical forest, savanna and grassland. In addition, our model highlights topographic features—\nparticularly high elevation, steep slopes, and well drained regions—as key predictors of low\nresilience and strong ENSO responses (Fig. 2f). These landscapes often support vegetation with\nshallow roots, limited soil depth, and poor water and nutrient retention, resulting in reduced\nbuffering capacity3. In parallel, climatic seasonality drives divergent vegetation adaptive\nstrategies, which in turn influence regional sensitivity to ENSO. While low-latitude regions with\nweak seasonality exhibit strong resilience shifts globally (Fig. 2e), ecosystems in areas with high\nintra-annual precipitation variability—such as northern Asian shrublands, Siberian tundra,\neastern Amazon rainforest, and Sahelian grasslands—also show heightened sensitivity compared\n42\nto adjacent regions with similar rainfall and vegetation types (Fig. 3a). This supports prior\nfindings1,4,5 that stable water availability enhances post-disturbance recovery.\nEqually important, anthropogenic disturbances exacerbate these patterns6. Highly fragmented\ntropical and temperate forests, such as those in the Amazon, Central Africa, and parts of North\nAmerica and Southeast Asia, exhibit heightened sensitivity to ENSO (Fig. 3a), likely due to\nintensified edge effects that increase local temperatures and reduce humidity, amplifying climate\nstressors7-9. Intensive agriculture drives land degradation, biodiversity loss, and hydrological\ndisruption, further destabilizing adjacent ecosystems10,11. For instance, we found that mixed\nforests along the Eurasian agro-ecological transition zone show more negative responses to La\nNiña when located near croplands (Fig. 3a), suggesting landscape configuration modulates\nENSO sensitivity—a mechanism that merits further empirical investigation.\nOur findings reinforce previous evidence that vegetation in climatically unstable environments—\ncharacterized by high variability in light, temperature, and moisture—tends to exhibit lower\nresilience. The identified hotspots of ENSO sensitivity are broadly consistent with earlier studies\non vegetation–climate responses12. For example, we found that vegetation resilience is notably\nlower in the southeastern Amazon, consistent with earlier studies that associate this decline with\nrecurring droughts and floods13,14. Similarly, reduced resilience along the Sahelian margins\nsupports earlier findings on the role of increasing water limitation15. ENSO-related effects\nobserved in Southeast Asia, Australia, Central Africa, and Central America also reflect patterns\ndocumented in prior research16-18. As climate change intensifies, such instability is expected to\nincrease, threatening ecosystem stability across large regions4. It is therefore critical to integrate\nENSO forecasts into ecosystem restoration and climate-informed conservation planning,\nespecially for tropical and subtropical forests. These ecosystems, given their high biodiversity,\ncarbon storage potential, and sensitivity to both climate and human disturbance, should be\nprioritized in global strategies for resilience and adaptation. In contrast, some northern forests,\nsuch as mixed forests in Northern Europe, show increasing resilience and minimal ENSO\nsensitivity (Fig. 1a and Fig. 3a), highlighting the need for regionally differentiated adaptation\nstrategies.\nplants survive while others succumb to drought? New Phytologist 178, 719-739,\ndoi:https://doi.org/10.1111/j.1469-8137.2008.02436.x (2008).\nto drought: Exploring resilience and resistance. CATENA 243, 108156,\ndoi:https://doi.org/10.1016/j.catena.2024.108156 (2024).\nTopographic, soil, and climate drivers of drought sensitivity in forests and shrublands of the\nPacific Northwest, USA. Sci Rep 10, 18486, doi:10.1038/s41598-020-75273-5 (2020).\nNat. Commun. 14, 498, doi:10.1038/s41467-023-36207-7 (2023).\nof the Amazon rainforest. Proceedings of the National Academy of Sciences 106, 20610-20615,\ndoi:doi:10.1073/pnas.0804619106 (2009).\nestimated. Nature 631, 570-576, doi:10.1038/s41586-024-07629-0 (2024).\ncanopies. Nat. Commun. 13, 917, doi:10.1038/s41467-022-28490-7 (2022).\nin forests. Nature Ecology & Evolution, doi:10.1038/s41559-025-02776-7 (2025).\neabo5003, doi:doi:10.1126/science.abo5003 (2023).\nPampa. Nat. Sustain. 6, 662-670, doi:10.1038/s41893-023-01074-z (2023).\ninsect biodiversity worldwide. Nature 605, 97-102, doi:10.1038/s41586-022-04644-x (2022).\nterrestrial ecosystems to climate variability. Nature 531, 229-232, doi:10.1038/nature16986\n(2016).\nof the National Academy of Sciences 114, 4442-4446, doi:10.1073/pnas.1617988114 (2017).\nmore variable rainfall. Nature Geoscience 12, 174-179, doi:10.1038/s41561-019-0312-z (2019).\nGlob. Change Biol. 28, 571-587, doi:https://doi.org/10.1111/gcb.15939 (2022).\nforests. Ecohydrology 12, e2056, doi:https://doi.org/10.1002/eco.2056 (2019).\nthe IOD and the mechanisms for impacts on Australian rainfall. Journal of Climate 24, 3910-\n3923 (2011).\nGlobal Ecology and Biogeography 29, 2039-2051, doi:https://doi.org/10.1111/geb.13174\n(2020).\ndoi:10.1093/nsr/nwy104 (2018).\nCentral Pacific El Niño. Geophys. Res. Lett. 48, e2020GL092367,\ndoi:https://doi.org/10.1029/2020GL092367 (2021).\nClimatology and Variations during the ENSO Developing and Decaying Phases. Journal of\nClimate 21, 6304-6320, doi:https://doi.org/10.1175/2008JCLI2331.1 (2008).\nEarth & Environment 2, 628-644, doi:10.1038/s43017-021-00199-z (2021).\nWIREs Climate Change 15, e861, doi:https://doi.org/10.1002/wcc.861 (2024).\nAtmosphere Interaction, Geophys. Monogr 147, 21-48 (2004).\ndestabilize plant species coexistence via changes in plant–soil feedback. Nature Ecology &\nwarming. npj Climate and Atmospheric Science 6, 8, doi:10.1038/s41612-023-00341-2 (2023).\ntropical carbon recovery. Glob. Change Biol. 30, e17423, doi:https://doi.org/10.1111/gcb.17423\n(2024).\nglobal tropics. Nat. Commun. 12, 5785, doi:10.1038/s41467-021-26048-7 (2021).\nshortwave radiation forcing over the tropical Pacific. Atmos. Chem. Phys. 8, 5565-5577,\ndoi:10.5194/acp-8-5565-2008 (2008).\nAfrican tropics. Journal of Geophysical Research: Biogeosciences 118, 1347-1358,\ndoi:https://doi.org/10.1002/jgrg.20115 (2013).\ndoi:10.1038/s41586-018-0252-6 (2018).\nteleconnections. Journal of Climate 10, 1769-1786 (1997).\n(2020).\nanomaly. Nat. Clim. Chang. 13, 967-974, doi:10.1038/s41558-023-01776-4 (2023).\nwet and dry hydrological extremes reveals onset of a local drought-induced tipping point. Glob.\nChange Biol. 29, 6077-6092, doi:https://doi.org/10.1111/gcb.16933 (2023).\nSystems in Eastern Australia. Geophys. Res. Lett. 50, e2023GL104814,\ndoi:https://doi.org/10.1029/2023GL104814 (2023).\nDerived From MODIS and AVHRR Surface Reflectance. IEEE Transactions on Geoscience\nand Remote Sensing 54, 5301-5318, doi:10.1109/TGRS.2016.2560522 (2016).\nproducts. Agric. For. Meteorol. 246, 218-230,\ndoi:https://doi.org/10.1016/j.agrformet.2017.06.016 (2017).\noccurrence. Proceedings of the National Academy of Sciences 121, e2316924121,\ndoi:doi:10.1073/pnas.2316924121 (2024).\nIndex (LAI): Methods, Products, Validation, and Applications. Reviews of Geophysics 57, 739-\n799, doi:https://doi.org/10.1029/2018RG000608 (2019).\n41 von Jeetze, P. J. et al. Projected landscape-scale repercussions of global action for climate and\nbiodiversity protection. Nat. Commun. 14, 2515, doi:10.1038/s41467-023-38043-1 (2023).\nusing a new indicator. Ecol. Indic. 129, 107872,\ndoi:https://doi.org/10.1016/j.ecolind.2021.107872 (2021).\nunder greenhouse warming. npj Climate and Atmospheric Science 7, 81, doi:10.1038/s41612-\n024-00627-z (2024).\nNat. Commun. 5, 3712, doi:10.1038/ncomms4712 (2014).\nsince the early 2000s. Nat. Clim. Chang. 12, 271-278, doi:10.1038/s41558-022-01287-8 (2022).\nLoess Plateau. Communications Earth & Environment 4, 336, doi:10.1038/s43247-023-01000-\n3 (2023).\nNature Ecology & Evolution 7, 1799-1808, doi:10.1038/s41559-023-02194-7 (2023).\nNature 631, 111-117, doi:10.1038/s41586-024-07568-w (2024).\n45\nresilience variability in a changing climate\nWe greatly appreciate the reviewers’ valuable and constructive suggestions. We have carefully\nrevised our manuscript based on their comments. We hope that the revisions adequately address any\nconcerns. Below, you will find our item-by-item responses to the comments. The following different\nfonts are used to facilitate the review.\n Contents underlined are the original comments of the Reviewers.\n Contents in normal Times New Roman font are our response to the Reviewers.\n Contents in blue and italic font are part of the revised manuscript; bold text indicates newly\nadded content in direct response to reviewer comments.\n Green text highlights where the revisions have been made in the manuscript.\n Red-labeled figures are supplementary illustrations included in the response.\nComment 1\nAbstract/Discussion: I still miss the implications/importance of your findings. Why should the\nscientific community care about these results.\nReply: A thorough understanding of vegetation resilience to climate variability is critical for\nsustaining ecosystem functions and terrestrial carbon sinks. Although previous studies have\ninvestigated the trends and spatial variation in global vegetation resilience, as well as its underlying\nenvironmental drivers, there is still a critical knowledge gap regarding the impacts of El Niño-Southern\nOscillation (ENSO) on global vegetation resilience. This gap hinders accurate projections of future\necosystem function, stability, and evolution under climate change. Our study addresses this by\nconducting a comprehensive global assessment of vegetation resilience sensitivity to ENSO,\nidentifying teleconnected hotspots and potential influential pathways. These findings will enable\ntargeted restoration and climate-informed ecosystem management, making them relevant to\nresearchers in climate science, atmospheric science, terrestrial ecology, and Earth system science.\nIn the revised Abstract and Introduction, we added a concluding sentence to state that our study\nprovides a coherent global assessment of vegetation resilience sensitivity to ENSO and identifies\nENSO-teleconnected hotspots and pathways, highlighting their relevance for restoration and climate-\ninformed governance. In the Discussion, we emphasize that the principal contribution of our work is\nto demonstrate ENSO’s role as a synchronizing agent of resilience loss in climatically unstable\necosystems. We also underscore the practical implications of our hotspot maps for integrating ENSO\nforecasts into ecological restoration, conservation planning, and regionally differentiated adaptation\nstrategies. These additions clarify both the scientific importance and policy relevance of our study.\nThe revised paper now reads:\n(Page 1 Lines 27-30, Abstract)\nOur study provides a coherent global assessment of vegetation resilience sensitivity to ENSO,\nidentifies teleconnected hotspots and potential influential pathways, and informs targeted\nrestoration and climate-adaptive ecosystem governance under climate change.\n(Page 3 Lines 93-97, Introduction)\nWe present a coherent global-scale evaluation of how vegetation resilience responds to ENSO,\npinpointing teleconnected hotspots and revealing key pathways of influence. These findings shed\nnew light on the teleconnection of vegetation growth and extreme climate, supporting nature-based\nsolutions and informing sustainable development efforts.\n(Pages 10-11 Lines 346-359, Discussion)\ncoherent global assessment of vegetation resilience sensitivity to ENSO. The identified ENSO\nsensitivity hotspots—such as the southeastern Amazon, Sahelian margins, Southeast Asia,\nAustralia, Central Africa, and Central America—broadly align with earlier studies linking reduced\nresilience to recurrent droughts, floods, and water limitations1-7, thereby validating our approach.\nMore importantly, the synthesis reveals that ENSO acts as a synchronizing agent, imposing a\ncommon rhythm of resilience loss across climatically unstable ecosystems worldwide. As climate\nchange intensifies, this global perspective becomes critical8. Our findings provide actionable\ninsights for identifying vulnerability hotspots and demonstrate the importance of integrating ENSO\nforecasts into ecosystem restoration and climate-informed conservation planning, with tropical and\nsubtropical forests prioritized for their biodiversity, carbon storage, and climatic sensitivity. In\ncontrast, some boreal forests, such as mixed forests in Northern Europe, show higher resilience and\nminimal ENSO sensitivity (Fig. 1a and Fig. 3a), indicating that ENSO need not be a central concern\nin these regions. These insights point toward the need for regionally differentiated strategies.\nComment 2\nLine 105: typing error, S should be removed: ‘Spatial patterns of historical changes in global\nvegetation resilience contrastS …’\nReply: We have corrected the typo by changing “contrastS” to “contrast.”\nThe revised paper now reads:\n(Page 3 Lines 110-111, Results)\nSpatial patterns of historical changes in global vegetation resilience contrast with the dynamics of\nglobal vegetation greening.\nComment 3\nLines 149-151: I don’t understand this sentence: ‘However, due to spatial heterogeneity in distance\nfrom oceanic forcing and topographic complexity, vegetation often exhibits lagged responses to\nENSO.’ What is spatial heterogeneity in distance from oceanic forcing and topographic complexity?\nReply: We accept our original wording was unclear. What we meant is that the distance of different\nregions from the tropical Pacific (the source region of ENSO signals) and their topographic features\n(e.g., elevation and terrain complexity) influence the propagation of ENSO impacts. These differences\noften lead to lagged vegetation responses to ENSO. We have revised the sentence for clarity in the\nrevised manuscript.\n(Page 5 Lines 159-161, Results)\nHowever, the lagged response of vegetation to ENSO is often a result of regional disparities,\nincluding distance from the tropical oceans and local topography like elevation and terrain\ncomplexity.\nComment 4\nAt line 152, you say that the global influence of extreme LN events surpasses that of ELN events by\nthe third year in shaping vegetation resilience patterns, but at line 193, you say that ‘beyond a nine-\nmonth lag, LN’s negative impact area already exceeds that of EN’. Why is it nine months in the one\ncase and three years in the other? And why is one in months and the other in years?\nReply: The apparent difference between the “nine-month” and “three-year” results is not a\ncontradiction but reflects the different purposes, scales, and resolutions of the two analytical\napproaches. To clarify, we have revised the text in the Results and Discussion sections to make this\ndistinction explicit.\nThe “nine-month” result (Line 193, now corrected to “twelve-month”) comes from a pixel-scale\ncorrelation analysis conducted at monthly resolution. This method is designed to detect the onset\ntiming and spatial extent of ENSO impacts. It shows that the spatial area negatively affected by La\nNiña surpasses that of El Niño within the first nine months, because El Niño’s direct footprint\ndiminishes with increasing lag, whereas La Niña’s footprint remains relatively stable.\nThe “three-year” result (Line 152) is derived from our global-scale Random Forest (RF)\nanalysis, based on annual data. Unlike correlation analysis, the RF framework evaluates the integrated,\ncumulative strength of ENSO’s influence on vegetation resilience. It indicates that La Niña’s overall\nimpact overtakes that of El Niño only by the third year (36-month lag). To ensure robustness, we\nadditionally conducted a complementary RF analysis using monthly data in our previous revision\n(Figure 1). The monthly RF model confirmed that the dominance of La Niña becomes evident at the\n36-month lag and is not yet observed at 24 months.\nIn summary, the two approaches are complementary: the pixel-scale monthly correlation analysis\nidentifies when and where La Niña begins to exceed El Niño in terms of spatial coverage, while the\nglobal-scale RF analysis captures when La Niña’s cumulative influence becomes stronger. Both\nconverge on the conclusion that with increasing lag time, La Niña’s influence intensifies and\neventually surpasses that of El Niño.\nbased on monthly variables. (a) Variables’ importance scores from the Random Forest (RF) model.\n(b-f) Partial dependence plots (PDPs) show the influence of key predictors from each variable group\non AC: (b) climatic variable, including ENSO, precipitation, temperature, and radiation; (c) Lag effects\nof ENSO at 3-, 6-, 9-, 12-, 18-, 24-, and 36- month intervals; (d) ENSO–climate interactions. Positive\nand negative interaction terms indicate climate anomalies associated with El Niño and La Niña phases,\nrespectively; (d) climatic seasonality, represented by intra-annual coefficient of variation (CV) of\nprecipitation, temperature, and radiation (e.g., temperature_seas indicates temperature seasonality);\nand (f) topographic factors, including DEM, Slope, and Height Above Nearest Drainage (Hand).\nThe revised paper now reads:\n(Page 5 Lines 161-163, Results)\nWith increasing lag time, the global impact of extreme La Niña events becomes more pronounced and,\nby the third year (the lag = 2), surpasses that of El Niño in shaping vegetation resilience dynamics\n(Fig. 2c).\n(Page 6 Lines 199-207, Results)\nThese spatial distributions remained consistent across varying lag times. However, it is noteworthy\nthat the spatial extent of El Niño's synchronous negative impact decreased with increasing lag time,\n5\nwhile the area experiencing negative impacts from La Niña remained relatively stable. Beyond a\ntwelve-month lag, La Niña's negative impact area exceeds that of El Niño (Extended Data Fig. 6),\nwhereas the global-scale Random Forest analysis indicates that La Niña’s cumulative impact\nbecomes dominant only after three years (Fig. 2c). This apparent difference reflects the distinct\nmetrics captured by the two approaches—spatial extent versus aggregated strength—yet both\nconsistently point to La Niña’s growing dominance with increasing lag time.\nComment 5\nLines 202: sometimes Amazon rainforest is with a capital R and sometimes with a lowercase, best to\nbe consistent.\nReply: We have standardized the terminology throughout the manuscript and now consistently use\n“Amazon rainforest” with a lowercase “r.”\nThe revised paper now reads:\n(Page 6 Lines 212-214, Results)\nwe reveal spatially heterogeneous coupling patterns, where ENSO-driven climate anomalies\nnegatively affected 23% of vegetation resilience, predominantly in the Amazon rainforest and\nBrazilian Plateau region in central South America,\nComment 6\nLine 221: what are non-teleconnected regions? If these are non-ENSO-affected regions, maybe it’s\nmore clear to just say that?\nReply: The term “non-teleconnected regions” in the original manuscript was intended to refer to\nregions not significantly influenced by ENSO. To improve clarity, we have revised the text throughout\nthe manuscript to use the phrase “non-ENSO-affected regions” instead of “non-teleconnected regions.”\nThis change ensures that the description of ENSO influence is more precise and easily understood by\nreaders.\nThe revised paper now reads:\n(Page 7 Lines 232-233, Results)\nTo analyze the impact of ENSO events on the long-term resilience trend, we compared the historical\nAC trend between ENSO-affected and non-ENSO-affected regions.\n6\n(Page 7 Lines 239-247, Results)\nRegardless of considering ENSO as a whole or its individual phases, the historical decline in\nvegetation resilience was significantly more pronounced in negatively affected regions compared to\npositively affected or non-ENSO-affected regions (p < 0.05) (Fig. 3b and Supplementary Table 1-3).\nIn contrast, the positive impact of ENSO did not lead to a corresponding increase in resilience.\nSpecifically, no significant difference was observed between positively affected and non-ENSO-\naffected regions under the overall impact of ENSO or La Niña (Fig. 3b and Supplementary Table 1\nand 3). Moreover, the resilience decline in positively affected regions under the impact of El Niño was\nhigher than that in non-ENSO-affected regions (Fig. 3b and Supplementary Table 2).\nComment 7\nLines 347-349: Can you provide a reference for this statement?\nReply: We have now added the appropriate reference to support this statement in the revised\nmanuscript.\nThe revised paper now reads:\n(Page 11 Lines 385-388, Methods)\nMoreover, LAI is more resistant to saturation than the Normalized Difference Vegetation Index\n(NDVI)9,10 and provides a longer and more consistent temporal record than Vegetation Optical\nDepth (VOD)11, making it a robust indicator for evaluating vegetation resilience over time.\nEvolution 7, 1790-1798, doi:10.1038/s41559-023-02187-6 (2023).\nMethods, Products, Validation, and Applications. Reviews of Geophysics 57, 739-799,\ndoi:https://doi.org/10.1029/2018RG000608 (2019).\nProceedings of the National Academy of Sciences 121, e2316924121,\ndoi:doi:10.1073/pnas.2316924121 (2024).\nComment 8\nLines 377:378: MODIS data has a citation that you can use instead of just including the DOI, and so\ndo the other datasets in this paragraph.\nReply: We have updated all dataset references in this paragraph to include full citations in accordance\n7\nwith Nature Communications style, replacing the previous DOI-only entries. All dataset references\n(e.g., GLASS LAI, NOAA SSTA, MODIS landcover, CMIP6, GLAD Intact Forest, DEM, and HAND)\nhave been updated to full citations.\nThe revised paper now reads:\n(Page 11 Lines 379-383, Methods)\nThe historical vegetation dynamic was indicated by the long-term leaf area index (LAI) product of the\nGlobal Land Surface Satellite (GLASS) V5 dataset, which was generated from Advanced Very High\nResolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS)\nsatellite data12.\n(Pages 11-12 Lines 393-395, Methods)\nENSO variability was represented by monthly sea surface temperature anomalies (SSTA) in the Niño\n3.4 region (5°S – 5°N and 170°W – 120°W), obtained from the NOAA ERSSTV5 dataset13.\n(Page 12 Lines 399-401, Methods)\nTo project future changes, we use the simulated climate and vegetation time series from the sixth phase\nof the Coupled Model Intercomparison Project (CMIP6)14, including monthly SST, LAI, near-surface\nair temperature, precipitation, and surface solar radiation.\n(Page 12 Lines 415-415, Methods)\nWe used landcover data from the Moderate Resolution Imaging Spectroradiometer (MODIS) land\ncover type product (MCD12C1)15.\n(Page 12 Lines 418-420, Methods)\nWe used the GLAD Intact Forest Landscapes dataset for the year 200016, which provides a global map\nof forests without significant disturbance.\n(Page 12 Lines 421-422, Methods)\nGlobal elevation data in the 2000 year was derived from the Shuttle Radar Topography Mission (SRTM)\nVersion 3 product, which offers a resolution of 30 m17.\n(Page 12 Lines 423-425, Methods)\nThe Height Above Nearest Drainage (HAND) data, obtained from Google Earth Engine18, were\nderived from SRTM elevation data and provide a more accurate representation of groundwater\naccessibility and hydrological buffering capacity.\nComment 9\nLines 433: typing error, S should be added: ‘This method robustly separateS …’\nReply: We have corrected a typing error by changing “separate” to “separates”.\n(Page 13 Lines 468-470, Methods)\nThis method robustly separates the LAI timeseries into three components: a long-term trend, a fixed\nseasonal cycle trend (by using ‘periodic’ option in s.window parameter), and a residual component.\nResponses to Reviewer #2\nComment 1\nThis manuscript addresses an important and timely topic by examining how ENSO influences global\nvegetation resilience, but despite improvements since the first round of review, the paper still requires\nsubstantial clarification and restructuring to ensure conceptual transparency, methodological rigor, and\nconsistency in the interpretation of results.\nReply: We sincerely thank the reviewer for their comments to help us further strengthen our paper.\n(1) To improve conceptual clarity, we now explicitly describe ENSO’s influence on vegetation\nresilience through two distinct pathways: (i) atmospheric synchronization teleconnection, the\nperiodic impacts of large-scale statistical association between ENSO events and resilience\nanomalies, analyzed independently of local climate variables; and (ii) climate-mediated\nteleconnection, defined as ENSO-driven anomalies in precipitation, temperature, and radiation\nthat subsequently alter vegetation resilience. These two pathways collectively reflect ENSO’s dual\nrole as a synchronizer of global variability and a modulator of local climate conditions. To further\nclarify this distinction, we have also added a schematic diagram (Extended Data Fig. 1) in the\nMethods section. (see Comments 5 and 9).\n(2) To improve methodological rigor, we carefully reviewed the entire manuscript and made the\nfollowing revisions. First, when describing the role of climatic factors, we replaced the term local\nwith the more precise term in-situ to avoid ambiguity regarding spatial scale (see Comment 8).\nSecond, we made minor revisions and added more details to the calculation of synchronization\nimpacts, such as removing the in-situ climate anomaly before calculation to obtain the pure ENSO\nsynchronization impact and the rules for determining the dominant event (see Comments 5 and\n9). Third, we added a discussion of the uncertainties associated with the critical slowing down\n(CSD) approach and potential strategies for improvement, thereby enhancing the transparency of\nresilience estimation (see Comment 4). Collectively, these revisions strengthen the\nmethodological clarity and reproducibility of our study without altering the underlying analyses.\n(3) To improve consistency in the interpretation of results, we thoroughly reviewed the Results and\nDiscussion sections to ensure that the interpretation of findings is consistent across the text,\nfigures,\n9\nand tables. Specifically, we clarified the relative contributions of synchronization and climate-\nmediated effects of ENSO to vegetation resilience by categorizing ENSO-affected regions into\nthree groups: Regions where only the synchronization effect (linked primarily to large-scale\nteleconnection with ENSO), Regions where only the climate-mediated effect (ENSO-driven\nclimate anomalies), and Regions where both synchronization and climate-mediated effects were\nstatistically significant (see Comment 2). In addition, we strengthened the global analysis based\non the Random Forest models by providing clearer rankings of variable importance, thereby\nrefining the identification of dominant controlling factors and their influence characteristics (see\nComment 6). Furthermore, we expanded the Discussion to provide a more detailed description of\nthe mechanisms by which ENSO influences vegetation resilience, highlighting both the climate-\nmediated pathways and two potential pathways (see Comment 11). Together, these revisions\nensure a consistent and transparent presentation of our findings throughout the manuscript.\nComment 2\nAbstract: I would recommend being more explicit here: where 68% of vegetation resilience anomalies\nare directly attributed to ENSO, while the remaining 25–32% result from ENSO–climate interaction.\nAlso in abstract: What does this ENSO–climate interaction mean? Being more explicit here gives\nmore opportunity to entice the readers.\nReply: In the revised manuscript, we first clarified the definitions of atmospheric synchronization\nteleconnection (previously termed the “direct effect”) and Climate-mediated teleconnection\n(previously referred to as the “ENSO–climate interactions”) and then specified their relative\ncontributions in the Abstract by breaking this down further. We now specify that ENSO affects 53%\nof the global vegetated area. Within these regions, 15% are linked to the large-scale atmospheric\ncirculation changes associated with ENSO (atmospheric synchronization teleconnection), 51% are\naffected mainly through ENSO-driven anomalies in local temperature, precipitation, and solar\nradiation (Climate-mediated pathway), and the remaining 34% are influenced by both processes.\nThe revised paper now reads:\n(Page 1 Lines 22-25, Abstract)\nOur findings reveal that ENSO significantly affects vegetation resilience across 53% of the global\nvegetated area. Within these regions, 15% are linked primarily to large-scale atmospheric\nsynchrony with ENSO, 51% are mainly shaped by ENSO-driven local climate anomalies, and the\nremaining 34% are influenced by both processes.\n(Pages 2-3 Lines 85-91, Introduction)\nperspectives: The atmospheric synchronization teleconnection, which captures the large-scale\nstatistical association between ENSO events and resilience anomalies, independently of local\nclimate variables; while the climate-mediated teleconnection, which reflects indirect influence\ntransmitted through ENSO-driven anomalies in temperature, precipitation, and radiation (see\nMethods and Extended Data Fig. 1). Together, these pathways reflect ENSO’s dual role as a\nsynchronizer of global variability and a modulator of local climate conditions.\n(Page 7 Lines 235-238, Results)\nWithin these regions, 15% are influenced directly by ENSO’s synchronization, 51% are affected\nsolely through ENSO-induced climate anomalies (for example, abnormal temperature,\nprecipitation, and radiation), and the remaining 34% are subject to both effects. In terms of ENSO\nphases, 61% of affected regions are linked to El Niño, while 39% are associated with La Niña.\nComment 3\nLine 51. Remove the word ENSO.\nReply: We have removed the redundant word “ENSO” from line 51 as suggested.\nThe revised paper now reads:\n(Page 2 Lines 50-52, Introduction)\nThrough large-scale atmospheric circulation anomalies, ENSO alters regional patterns of\ntemperature, precipitation, and solar radiation19,20, thereby disrupting vegetation biophysical\nprocesses21.\nComment 4\nLine 82. You are proposing that you can estimate resilience based on the critical slowing down theory.\nWould you be able to discuss the opportunities and challenges in estimating resilience with this theory\nafter you introduce the term critical slowing down theory? It would make your results more transparent\nand demonstrate that uncertainties could arise from your estimates. In the Discussion you could also\ndevelop these uncertainties and how to minimize them.\nReply: As suggested, we have developed an in-depth discussion of the opportunities and challenges\nin estimating resilience using critical slowing down theory, uncertainties, and ways to minimize them.\nThese additions improve the transparency of our resilience estimates.\n(Page 11 Lines 366-375, Discussion)\nMethodologically, critical slowing down (CSD) indicators—such as increasing variance or\nautocorrelation—provide a theoretical framework for detecting system stability changes without\nabrupt disturbances22. However, they capture only theoretical resilience shifts rather than actual\nrecovery rates or critical transitions23. Moreover, the application of CSD indicators requires long\ntime series, yet large-scale observational records often combine data from multiple sensors with\ndifferent processing steps, introducing biases in trends and significance estimates 24,25. Noise and\nnon-stationarity can further inflate signals, causing false positives. More conservative significance\ntesting (for example, surrogate time-series methods)26, and cross-validation across multiple datasets\nand indicators would strengthen the robustness of resilience estimates.\nComment 5\nLine 85. Can you give an example of the direct effects? You mentioned “which may reflect\ncomponents of ENSO’s impact not fully mediated by local climate variables,” but what are these\ncomponents?\nReply: Thank you for pointing out this ambiguity. Your question has helped us refine the conceptual\nclarity of the direct effect.\nOur distinction between “direct” and “indirect” effects is not meant to imply that these processes are\nphysically separable, but rather provides an operational framing to distinguish complementary\ninfluences. To avoid ambiguity and enhance conceptual clarity, in the revised manuscript, we define\nclimate-mediated teleconnection27,28 (previously referred to as the “indirect effect”) represents the\ninfluence of ENSO transmitted through ENSO-driven anomalies in local temperature, precipitation,\nor solar radiation—an impact pathway with explicit logical inference. The atmospheric\nsynchronization teleconnection29-31 (previously termed the “direct effect”) refers to the statistical\nassociation between ENSO indices and vegetation resilience analyzed independently of local climate\nvariables. This pathway reflects ecosystem responses that are temporally aligned with ENSO rhythms,\nincluding lagged effects, and likely arise from large-scale atmospheric teleconnections (e.g.,\ncirculation or humidity shifts). It captures phase-coherent vegetation dynamics with ENSO,\neven without concurrent local climate anomalies.\nTo the best of our knowledge, no study has directly and quantitatively explored the full spectrum of\nENSO's influencing mechanisms on vegetation resilience using a systematic partitioning framework\nlike ours. Existing literature has primarily focused on the drivers of resilience change (e.g., land use\nchange, disturbance history, water availability) or the ENSO impact on ecosystem functions (e.g.,\n12\nproductivity). Previous research indicates ENSO-induced large-scale circulation anomalies can trigger\nregional extremes such as droughts, floods, and heatwaves, or shift hydroclimatic regimes. These\nevents alter fundamental ecosystem controls—light, temperature, and water availability—that regulate\nphotosynthesis, respiration, and productivity, thereby leading to observable changes in vegetation\nresilience. In addition, we speculate that ENSO may also affect vegetation resilience through\nadditional mechanisms. For example, ENSO-related circulation changes might influence long-range\naerosol transport and deposition, thereby altering radiation balance and nutrient inputs. Likewise,\nENSO-driven fluctuations in commodity prices could, in principle, indirectly shape land-use practices\nand management decisions. While these mechanisms remain conjectural, they illustrate potential\navenues through which ENSO’s influence on resilience could extend beyond the direct modulation of\nlocal climate anomalies.\nWe emphasize that while our dataset does not allow us to directly confirm these specific non-climatic\nmechanisms, we highlight them in the Discussion as plausible pathways and emphasize that further\nresearch is critically needed to test their contribution. Our primary contribution is the spatial\npartitioning of the total ENSO influence based on the dominance of the established climatic versus the\nunmediated macro-scale influences.\nThe revised paper now reads:\n(Pages 2-3 Lines 85-91, Introduction)\nWe hypothesized that ENSO influences on global vegetation resilience operate via two complementary\nperspectives: The atmospheric synchronization teleconnection, which captures the large-scale\nstatistical association between ENSO events and resilience anomalies, independently of local\nclimate variables; while the climate-mediated teleconnection, which reflects indirect influence\ntransmitted through ENSO-driven anomalies in temperature, precipitation, and radiation (see\nMethods and Extended Data Fig. 1). Together, these pathways reflect ENSO’s dual role as a\nsynchronizer of global variability and a modulator of local climate conditions.\n(Page 10 Lines 302-320, Discussion)\nENSO events exert asymmetric teleconnection effects on vegetation resilience, with El Niño–\nassociated regions exhibiting more pronounced declines (Fig. 2). Specifically, over one-third of\nvegetated areas experience significant resilience changes during El Niño events, compared to only\nabout 20% during La Niña (Fig. 3). Mechanistically, these effects are primarily mediated by\nregional climate: ENSO-driven anomalies in temperature, precipitation, and radiation trigger\nextremes—heatwaves, droughts, floods—can impair photosynthesis, water regulation, and root\nfunction32,33, ultimately slow canopy recovery. Our Random Forest analysis confirms that such\nanomalies, particularly those linked to El Niño (high temperature and radiation, water stress) and\nLa Niña (cold and dry conditions), are significantly associated with resilience decline globally (Fig.\n2). Even intact, high-biomass forests are vulnerable: our targeted analysis shows that resilience loss\nin these systems under ENSO-related stresses (Extended Data Fig. 4), consistent with previous\nfindings that tall canopies and deep-root systems depend heavily on stable water supply and become\nvulnerable under persistent ENSO pressure34. In addition to climate-mediated effects, ENSO may\ninfluence ecosystems through other pathways. For example, ENSO-driven circulation anomalies\ncan intensify savanna burning27 and facilitate long-range smoke transport to adjacent forests35,36,\nreducing incoming solar radiation and limiting canopy productivity37. ENSO-modulated shifts in\ncommodity prices (for example, palm oil)38,39, may also accelerate tropical deforestation and\nfragmentation40, further weakening forest resilience. These hypothesized pathways remain\nspeculative and warrant further investigation.\n(Page 15 Lines 529-542, Methods)\nTo map the spatial distribution of ENSO’s synchronization impact on vegetation resilience, we first\nremoved the effects of in-situ climate by fitting, at each grid cell, a multiple linear regression of\nresilience on precipitation, temperature, and radiation. We then used a Chi-square test, a\nnonparametric test method suitable for categorical variables41, to examine whether El Niño or La\nNiña events triggered anomalies in adjusted vegetation resilience. Positive and negative anomalies\nin vegetation resilience were defined as standardized AC anomalies (Z-scores) exceeding +2 or falling\nbelow –2, respectively, corresponding approximately to the upper and lower 5% tails of the\ndistribution. Similarly, the definition of El Niño or La Niña events was defined based on the commonly\naccepted criterion: a monthly SSTA greater than 1 °C or less than -1 °C42,43. El Niño and La Niña\nwere tested separately. When both were significant for a cell (FDR-adjusted P<0.05), we labeled the\ndominant phase as the one showing the stronger association. To analyze the lagged effect, Chi-square\ncoefficients were calculated between the ENSO event time series and AC anomaly time series in a lag\nof 3, 6, 9, 12, and 24 months (Extended Data Fig. 6). We also performed a similar analysis for the\nthree scenarios in the future period (2015–2100) (Fig. 4b). This approach helps reveal the spatial\nimprint of ecosystem responses that are temporally aligned with ENSO rhythms, including delayed\neffects, beyond what can be explained by in-situ climate variability.\nComment 6\nLine 117. “The analysis reveals that local controlling variables—namely topography and climate\nseasonality—play a more dominant role in shaping interannual variation in vegetation resilience than\nENSO teleconnections (Fig. 2a; model R² = 0.79).” Did you mention topography and climate\nseasonality because they were the most important group of variables controlling vegetation resilience?\nI suggest outlining the most important variables here (say, five) in a more organized fashion.\nReply: In the original version, we included two groups of controlling variables—topography and\nclimate seasonality—because they were designed as static background constraints rather than the main\nexplanatory variables of interest. However, our results indicated that these two groups collectively\nexert substantial influence on vegetation resilience, which is why they were highlighted in the text\n(Fig. 2a).\nIn the revised manuscript, we clarified this design choice and, in response to the reviewer’s suggestion,\nprovided a ranked assessment of the variable groups. Local static constraints related to climate\nseasonality have the greatest importance, followed by in-situ dynamic climate factors, topography, and\n14\nfinally ENSO teleconnections. We added a separate paragraph describing the key controlling variables\nwithin these groups. Specifically, temperature and radiation seasonality emerged as the dominant\ncomponents of climate seasonality, while elevation proved to be the primary topographic factor. We\nalso elaborated on the distinct effects associated with these features, thereby making the role of\ncontrolling variables more explicit and organized.\nThe revised paper now reads:\n(Page 3 Lines 120-130, Results)\nWe used a random forest (RF) regression framework at the global scale to attribute temporal\nchanges in vegetation resilience (AC) to ENSO variability and other environmental factors. This\nmodel focused on in-situ climate variability and ENSO as the key explanatory variables,\nincorporating other static in-situ factors such as topography and climate seasonality as constraints\nto reduce background interference (see Methods). The analysis reveals that climate seasonality\nexerts the strongest influence on shaping interannual variation, followed by climatic variability and\ntopography (Fig 2a; model R2 = 0.79). Among the controlling variables, temperature and radiation\nseasonality exert the strongest climatic control on resilience, and ecosystems experiencing low intra-\nannual variability are especially prone to resilience loss. Elevation emerges as the key terrain factor,\nwith partial dependence analyses showing that high, steep, and well-drained landscapes are\nparticular hotspots of vulnerability (Fig. 2d and e).\nComment 7\nLine 119. Unless you want to have Results and Discussion together, this sentence - “In particular,\nseasonal climate regimes influence resilience by selecting species with distinct adaptive traits.” - is\nmore like Discussion, because your study does not involve trait analysis.\nReply: We agree that the original sentence could be interpreted as a discussion point rather than a\ndirect result, since our study does not explicitly analyze plant traits. To address this, we revised the\nResults section to focus solely on the statistical attribution.\nThe revised paper now reads:\n(Page 3 Lines 126-130, Results)\nAmong the controlling variables, temperature and radiation seasonality exert the strongest climatic\ncontrol on resilience, and ecosystems experiencing low intra-annual variability are especially prone\nto resilience loss. Elevation emerges as the key terrain factor, with partial dependence analyses\nshowing that high, steep, and well-drained landscapes are particular hotspots of vulnerability (Fig.\n2d and e).\nComment 8\nLine 124. “Among local climatic factors.” What does local climatic factors mean? And then later\n15\nin line 129 you mention these climatic variables at a global scale. Scale is very important, and results\ncan vary substantially based on scale. I suggest you be very careful and specific throughout the paper\nregarding the use of local and global scales.\nReply: Thank you for highlighting the ambiguity of “local climatic factors” and emphasizing the\nimportance of clearly defining scale—this feedback has helped us refine the clarity and precision of\nour terminology.\nTo clarify, we recognize that the term “local” may have been confusing, as our analysis is global in\nscope. In the revised manuscript, we have either removed “local” or replaced it with “in-situ” to more\nprecisely refer to climate variables measured at the grid-cell scale. These changes ensure a clear\ndistinction between direct climatic controls and remote ENSO teleconnection effects throughout the\nmanuscript.\nWithin this framework, we further distinguish between: baseline in-situ climatic factors (excluding\nENSO effects), representing grid-scale climatic fluctuations unrelated to ENSO (e.g., seasonal cycles,\nlocal processes); and ENSO-driven in-situ climate anomalies, representing the portion of those\nvariables that covaries with ENSO indices. This clarification ensures consistent terminology\nthroughout the manuscript: “in-situ” always denotes grid-cell scale climate data, while modifiers\n(“baseline” vs. “ENSO-driven”) specify whether ENSO influence is included. This two-part\ndistinction enables us to separate ENSO’s teleconnection signal from its modulation of local climate\nconditions.\nThe revised paper now reads:\n(Page 1 Lines 35-36, Introduction)\nPrevious studies have predominantly focused on the relationship between resilience and temperature\nand precipitation patterns.\n(Pages 2-3 Lines 85-91, Introduction)\nWe hypothesized that ENSO influences on global vegetation resilience operate via two complementary\nperspectives: The atmospheric synchronization teleconnection, which captures the large-scale\nstatistical association between ENSO events and resilience anomalies, independently of local climate\nvariables; while the climate-mediated teleconnection, which reflects indirect influence\ntransmitted through ENSO-driven anomalies in temperature, precipitation, and radiation.\n(Page 3 Lines 121-124, Results)\nThis model focused on in-situ climate variability and ENSO as the key explanatory variables,\nincorporating other static in-situ factors such as topography and climate seasonality as constraints\nto reduce background interference (see Methods).\n16\n(Page 3 Lines 131-133, Results)\nFocusing on baseline in-situ climatic factors (excluding ENSO effects), temperature had the\nstrongest influence on resilience, followed by radiation, while precipitation had a relatively minor\neffect (Fig. 2a).\n(Page 4 Lines 150-151, Results)\nAlthough in-situ factors collectively explain a larger proportion of the resilience variation, ENSO\nteleconnections remain a critical global-scale forcing, driving resilience anomalies.\n(Page 9 Lines 293-294, Discussion)\nUsing an RF model, we disentangled the effects of in-situ climate variables (temperature, precipitation,\nand solar radiation) from those of ENSO and their interactions, revealing how ENSO shapes resilience\npatterns across global, regional, and biome-specific scales.\n(Page 10 Lines 337-341, Discussion)\nHighly fragmented tropical and temperate forests, such as those in the Amazon, Central Africa, and\nparts of North America and Southeast Asia, exhibit heightened sensitivity to ENSO (Fig. 3a), likely\ndue to intensified edge effects that increase temperatures and reduce humidity, amplifying climate\nstressors44-46.\n(Page 15 Lines 546-548, Methods)\nWe then assessed the climate-mediated impact of ENSO, representing the portion of ENSO's influence\non resilience that is mediated through in-situ climatic variables (precipitation, temperature, and solar\nradiation).\nComment 9\nLine 143. “ENSO primarily influences vegetation resilience through its interactions with climate\nvariables, which exert a stronger impact than ENSO’s direct effects.” My comment is related to a\nprevious comment on the direct versus indirect effects of ENSO. Can you be more explicit about what\nyou mean by ENSO’s interactions with climatic variables versus the direct effects? Are the anomalies\nin precipitation, temperature, and radiation not part of the direct effects of ENSO? I cannot\nunderstand the rest of the paragraph without better understanding this. Even in the Methods where you say\n“ENSO variability and resilience over time, which we refer to as the ‘direct impact’,” this is still\nvague. What is ENSO variability? Is it based on temperature, precipitation, radiation?\nReply: We agree that ENSO inevitably influences vegetation through its modulation of local climate\nanomalies. Our distinction between “direct” and “indirect” effects is therefore not meant to imply that\nthese processes are physically separable, but rather provides an operational framing to distinguish\ncomplementary influences.\nENSO influences vegetation resilience in the revised manuscript. Climate-mediated teleconnection27,28\n(previously referred to as the “indirect effect”) represents the influence of ENSO transmitted through\nENSO-driven anomalies in local temperature, precipitation, or solar radiation—an impact pathway\nwith explicit logical inference. The atmospheric synchronization teleconnection29-31 (previously\ntermed the “direct effect”) refers to the statistical association between ENSO indices and vegetation\nresilience analyzed independently of local climate variables. This pathway reflects ecosystem\nresponses that are temporally aligned with ENSO rhythms, including lagged effects, and likely arise\nfrom large-scale atmospheric teleconnections (e.g., circulation or humidity shifts). It captures phase-\ncoherent vegetation dynamics with ENSO, even without concurrent local climate anomalies. Together,\nthis framing allows us to disentangle ENSO’s role as both a synchronizer of global variability and a\nmodulator of local climate conditions.\nTo quantify these two perspectives, we combined: (i) a global analysis using a random forest\nmodel, where ENSO indices capture synchronization and interaction terms (ENSO × climate) capture\nclimate-mediated effects; and (ii) a spatially explicit analysis, where synchronization effects are\ndetected as ENSO–resilience correlations excluding local climate variables, and climate-mediated\neffects are identified as overlapping regions of significant ENSO–climate and climate–resilience\nrelationships.\n(Note: To enhance methodological rigor and operational clarity, we enhanced the calculation of the\nsynchronization effect across spatial scales. We first removed the influence of in-situ climate\nvariability on resilience via multiple regression to obtain residual resilience; we then applied chi-\nsquare contingency tests to El Niño and La Niña separately. When both phases were significant, we\nassigned dominance to the phase with the larger absolute effect size, thereby making both the\npreprocessing and the dominance rule explicit. These updates led to small numerical changes but did\nnot alter the overall conclusions.)\nFor clarity, we have also replaced “interactions with climate variables” with “climate-mediated\nteleconnection, driven by ENSO-induced climate anomalies” and “ENSO variability” with “ENSO\nevents” throughout the manuscript.\nThe revised paper now reads:\n(Page 1 Lines 22-25, Abstract)\nOur findings reveal that ENSO significantly affects vegetation resilience across 53% of the global\nvegetated area. Within these regions, 15% are linked primarily to large-scale atmospheric\nsynchrony with ENSO, 51% are mainly shaped by ENSO-driven local climate anomalies, and the\nremaining 34% are influenced by both processes.\n(Pages 2-3 Lines 85-91, Introduction)\nperspectives: The atmospheric synchronization teleconnection, capture reflects the large-scale\nstatistical association between ENSO events and resilience anomalies, independently of local\nclimate variables; while the climate-mediated teleconnection, which reflects indirect influence\ntransmitted through ENSO-driven anomalies in temperature, precipitation, and radiation (see\nMethods and Extended Data Fig. 1). Together, these pathways reflect ENSO’s dual role as a\nsynchronizer of global variability and a modulator of local climate conditions.\n(Pages 4-5 Lines 151-159, Results)\nENSO primarily influences vegetation resilience through the climate-mediated teleconnection,\ndriven by ENSO-induced climate anomalies, which exert a stronger impact than ENSO’s\natmospheric synchronization teleconnection (Fig. 2a). Specifically, reduced vegetation resilience is\nclosely associated with anomalously high temperatures, low precipitation, and enhanced solar\nradiation during El Niño events, as well as with low-temperature and low-precipitation anomalies\nunder La Niña conditions (Fig. 2d). While these climate-mediated effects dominate overall, the\natmospheric synchronization teleconnections of ENSO still matter: vegetation responds more rapidly\nand strongly to El Niño events than to La Niña during the event year (Fig. 2b).\n(Page 5 Lines 174-179, Results)\nFig. 2. Influence of explanatory variables on global vegetation resilience from 1981 to 2018. …(d)\nClimate-mediated teleconnection effects, calculated as the product of ENSO and in-situ climate\nvariables.\n(Page 6 Lines 185-197, Results)\nENSO synchronization teleconnection affecting vegetation resilience\nWe applied a Pearson Chi-square test to evaluate whether vegetation resilience is temporally\nsynchronized with ENSO over a long historical period, after removing the effects of in-situ climate\nvariables (see Methods). The resulting spatial pattern indicates that 26% of vegetated regions\nexhibited statistically significant relationships between vegetation resilience anomalies and the ENSO\nevents (P<0.05), with resilience demonstrating higher sensitivity to El Niño than La Niña events (Fig.\n3a). Specifically, 6% of vegetated areas showed increased resilience during La Niña events, while 3.3%\nexhibited decreased resilience. During El Niño events, 9.4% of areas experienced increased resilience\nand 7.4% experienced decreased resilience (Fig. 3a). The detrimental impacts of El Niño were\npredominantly concentrated in northern Asia, Europe and America, and the southern Amazon,\nwhereas the negative effects of La Niña were primarily distributed in southern Amazon and Central\nAsia and Africa (Fig. 3a). Evergreen broadleaf forests, grasslands and open shrublands were\nidentified as the vegetation types most affected by these climatic phenomena (Extended Data Fig. 5).\n(Page 6 Lines 200-205, Results)\nHowever, it is noteworthy that the spatial extent of El Niño's synchronous negative impact decreased\nwith increasing lag time, while the area experiencing negative impacts from La Niña remained\nrelatively stable. Beyond a twelve-month lag, La Niña's negative impact area exceeds that of El Niño\n(Extended Data Fig. 6), whereas the global-scale Random Forest analysis indicates that La Niña’s\n19\ncumulative impact becomes dominant only after three years (Fig. 2c). This apparent difference reflects\nthe distinct metrics captured by the two approaches—spatial extent versus aggregated strength—yet both\nconsistently point to La Niña’s growing dominance with increasing lag time.\n(Page 6 Line 209-214, Results)\nENSO climate-mediated teleconnection affecting vegetation resilience\nRF inference (Fig. 2c-e) suggests that ENSO-induced climate anomalies trigger resilience changes.\nBy overlapping resilience-climate (Extended Data Fig. 7) and ENSO-climate (Extended Data Fig. 8)\nlinkages, we reveal spatially heterogeneous coupling patterns, where ENSO-driven climate\nanomalies negatively affected 23% of vegetation resilience, predominantly in the Amazon rainforest\nand Brazilian Plateau region in central South America,\n(Page 7 Lines 233-238, Results)\nENSO-affected regions, defined as areas experiencing synchronization or climate-mediated impacts,\ncover 53% of global vegetation. Within these regions, 15% are influenced directly by ENSO’s\nsynchronization, 51% are affected solely through ENSO-induced climate anomalies (for example,\nabnormal temperature, precipitation, and radiation), and the remaining 34% are subject to both\neffects. In terms of ENSO phases, 61% of affected regions are linked to El Niño, while 39% are\nassociated with La Niña.\n(Page 8 Lines 250-258, Results)\nFig 3. ENSO impact on global vegetation resilience. (a) Synchronization impact of ENSO on\nresilience. … (c-h) Climate-mediated teleconnection effects on vegetation resilience: (c-e) Impacts\nof ENSO-induced anomalies in temperature, solar radiation, and precipitation on vegetation resilience,\n20\nrespectively.\n(Page 8 Lines 271-275, Results)\nRegarding ENSO teleconnections, projected warming is expected to weaken the synchronization\nimpact of ENSO on vegetation resilience, decreasing from 21.4% under SSP126 to 21.1% under\nSSP245 and 13.5% under SSP370 (P<0.05) (Fig. 4b). In parallel, ENSO-driven climate anomalies\nis projected to strengthen, affecting approximately 51-55% of global vegetation areas and marking\na 6-10% increase compared to historical levels (Fig. 4c).\n(Page 9 Lines 282-289, Results)\nFig. 4 Projected vegetation resilience trend and ENSO’s synchronization and climate-mediated\nteleconnection impacts on vegetation resilience under future scenarios (2015–2100). (a) Resilience\ntrend under SSP 126, SSP 245, and SSP 370, respectively. (b-c) ENSO impacts on vegetation\nresilience under the three scenarios, including (b) atmospheric synchronization effects and (c)\nclimate-mediated effects transmitted through ENSO-driven anomalies. Plus (+) and minus (-) signs\nindicate positive and negative impacts, respectively, with numbers indicating impact intensity,\ncalculated as the sum of individual grid-level impacts.\n(Page 9 Lines 296-298, Discussion)\nOne key finding is that 53% of global vegetation over the past four decades has been influenced by\nENSO, either through synchronization or climate-mediated effects.\n(Page 10 Lines 302-320, Discussion)\nENSO events exert asymmetric teleconnection effects on vegetation resilience, with El Niño–\nassociated regions exhibiting more pronounced declines (Fig. 2). Specifically, over one-third of\nvegetated areas experience significant resilience changes during El Niño events, compared to only\nabout 20% during La Niña (Fig. 3). Mechanistically, these effects are primarily mediated by\nregional climate: ENSO-driven anomalies in temperature, precipitation, and radiation trigger\nextremes—heatwaves, droughts, floods—can impair photosynthesis, water regulation, and root\nfunction32,33, ultimately slow canopy recovery. Our Random Forest analysis confirms that such\nanomalies, particularly those linked to El Niño (high temperature and radiation, water stress) and\nLa Niña (cold and dry conditions), are significantly associated with resilience decline globally (Fig.\n2). Even intact, high-biomass forests are vulnerable: our targeted analysis shows that resilience loss\nin these systems under ENSO-related stresses (Extended Data Fig. 4), consistent with previous\nfindings that tall canopies and deep-root systems depend heavily on stable water supply and become\nvulnerable under persistent ENSO pressure34. In addition to climate-mediated effects, ENSO may\ninfluence ecosystems through other pathways. For example, ENSO-driven circulation anomalies\ncan intensify savanna burning27 and facilitate long-range smoke transport to adjacent forests35,36,\nreducing incoming solar radiation and limiting canopy productivity37. ENSO-modulated shifts in\ncommodity prices (for example, palm oil)38,39, may also accelerate tropical deforestation and\nfragmentation40, further weakening forest resilience. These hypothesized pathways remain\nspeculative and warrant further investigation.\n(Page 10 Lines 321-323, Discussion)\nAlthough ENSO affects vegetation resilience through these synchronization or climate-mediated\nmechanisms, ecosystem responses vary due to interactions between ENSO’s asymmetric\nteleconnections, ecosystem traits, and human activities.\n(Page 14 Lines 508-512, Methods)\nThe regression model incorporated an interaction term between ENSO and climate factors, as well\nas a lag term. Climate factors from the preceding three years represented the lagged effect, while\nthe interaction terms (temperature × ENSO, precipitation × ENSO, and radiation × ENSO) characterized\nhow resilience sensitivity to climate changes varies with ENSO intensity.\n(Page 15 Lines 529-544, Methods)\nTo map the spatial distribution of ENSO’s synchronization impact on vegetation resilience, we first\nremoved the effects of in-situ climate by fitting, at each grid cell, a multiple linear regression of\nresilience on precipitation, temperature, and radiation. We then used a Chi-square test, a\nnonparametric test method suitable for categorical variables41, to examine whether El Niño or La\nNiña events triggered anomalies in adjusted vegetation resilience. Positive and negative anomalies\nin vegetation resilience were defined as standardized AC anomalies (Z-scores) exceeding +2 or falling\nbelow –2, respectively, corresponding approximately to the upper and lower 5% tails of the\ndistribution. Similarly, the definition of El Niño or La Niña events was defined based on the commonly\naccepted criterion: a monthly SSTA greater than 1 °C or less than -1 °C42,43. El Niño and La Niña\nwere tested separately. When both were significant for a cell (FDR-adjusted P<0.05), we labeled the\ndominant phase as the one showing the stronger association. To analyze the lagged effect, Chi-square\ncoefficients were calculated between the ENSO event time series and AC anomaly time series in a lag\nof 3, 6, 9, 12, and 24 months (Extended Data Fig. 6). We also performed a similar analysis for the\nthree scenarios in the future period (2015–2100) (Fig. 4b). This approach helps reveal the spatial\nimprint of ecosystem responses that are temporally aligned with ENSO rhythms, including delayed\neffects, beyond what can be explained by in-situ climate variability.\n(Page 15 Lines 546-556, Methods)\nWe then assessed the climate-mediated impact of ENSO, representing the portion of ENSO's\ninfluence on resilience that is mediated through in-situ climatic variables (precipitation,\ntemperature, and solar radiation). This was quantified by evaluating how ENSO-driven climate\nanomalies affected vegetation resilience. Linear regression was used to measure the co-variation\nbetween climatic variability and AC dynamics (Extended Data Fig. 7), and Chi-square tests were\napplied to assess the association between ENSO events and climatic anomalies (P < 0.05, Extended\nData Fig. 8). By spatially overlaying these relationships, the impact on resilience was coded as either\npositive (“1”) or negative (“-1”). For example, if resilience and precipitation were positively\ncorrelated, and ENSO events predominantly caused negative precipitation anomalies in the region,\nthe ENSO-precipitation impact on resilience was recorded as “-1”. We summed the climate-mediated\neffects of key ENSO–climate variables to quantify their overall impact (Fig. 3c–h).\n(Pages 15-16 Lines 560-566, Methods)\nSimilar to the climate-mediated impacts, we also quantified the positive and negative atmospheric\nsynchronization impacts as “+1” and “-1”, respectively, and then summed both synchronization\n(Fig. 3a) and climate-mediated (Fig. 3h) impacts to derive the total impact. Based on these total\nimpacts, we classified the global vegetated areas into three categories for each ENSO type (El Niño,\nLa Niña, and overall ENSO): (i) ENSO-positively affected regions, where the sum of synchronization\nand climate-mediated effects is greater than zero;\n(Page 21 Lines 832-834, Methods)\nExtended Data Fig. 1. Two pathways by which ENSO affects vegetation resilience: atmospheric\nsynchronization and in-situ climate mediation.\n(Page 22 Lines 844-849, Extended Data)\nExtended Data. 4. Influence of explanatory variables on intact forest vegetation resilience from 1981\nto 2018. …(d) Climate-mediated teleconnection effects, calculated as the product of ENSO and\nclimate variables.\n(Page 23 Line 856-858, Extended data)\nExtended Data Fig. 5. Area affected by ENSO's atmospheric synchronization across different\nvegetation cover types.\n(Page 23 Line 864-867, Extended data)\nExtended Data Fig. 6. ENSO’s atmospheric synchronization teleconnection and vegetation\nresilience at different time lags. (a-f) Spatial map of significant synchronization teleconnections for\nlag ranging from 0 to 24 months.\n(Page 25 Lines 881-885, Extended Data)\nExtended Data Fig. 9. Areas of different vegetation types affected by ENSO through climate-mediated\nteleconnections. Blue and orange represent positive and negative impacts, respectively. (a) represents\nthe total area affected by ENSO-driven climate anomalies, while (b) and (c) show the impact areas\nassociated with El Niño- and La Niña-induced climate anomalies, respectively.\n(Page 26 Lines 895-896, Extended Data)\nExtended Data Fig. 10. Distribution of the predicted impacts of ENSO-driven climate anomalies on\nglobal resilience under SSP126 (a), SSP245 (b), and SSP370 (c).\nComment 10\nLine 158. “which imply that anthropogenic disturbances may modulate but do not fundamentally\nalter the underlying climate–vegetation coupling driven by ENSO.” Just excluding the areas does\nnot leave room for this interpretation. There is a lot of field evidence based on tree and trait\nmeasurements showing different effects of ENSO events on intact versus disturbed forests.\nReply: In fact, field studies have documented differential ENSO or climate extreme effects between\nintact and disturbed forests. For example, previous studies found that mortality rates in disturbed\nforests remained elevated for years following extreme drought events, compared to less-disturbed\nforests47,48. However, these findings also demonstrate that intact forests are not immune to ENSO-\ndriven anomalies. Both forest types exhibit sensitivity to ENSO-related stressors, reinforcing the\nbroad-scale nature of the ENSO–vegetation resilience linkage observed in our study.\nWe have revised the sentence to avoid overstating the implications. We intended to highlight that the\nlarge-scale ENSO–vegetation teleconnection remains evident even when excluding regions affected\nby human activities, but we agree that this does not imply anthropogenic disturbances have no\nfundamental effects. The revised text now states that anthropogenic disturbances can influence local\nresponses, while the large-scale ENSO–vegetation coupling remains detectable.\nThe revised paper now reads:\n(Page 5 Lines 169-171, Results)\nThis suggests that while anthropogenic disturbances can influence vegetation responses at local\nscales, the large-scale climate–vegetation coupling associated with ENSO persists.\nComment 11\nThe Discussion is still very limited and did not attempt to discuss the mechanisms by which forest\nresilience is affected by ENSO.\nReply: We have revised the Discussion to explicitly address the mechanisms by which ENSO affects\nforest resilience. First, we describe the climate-mediated pathway, in which ENSO-driven anomalies\nin temperature, precipitation, and radiation impair vegetation functioning and delay recovery. We\nsupport this with new evidence from our Random Forest analysis, which links El Niño–associated\nhigh temperatures and radiation, as well as La Niña–related cold and dry conditions, to reduced\nresilience globally. We also highlight that even high-biomass forests are vulnerable, as demonstrated\nby our intact forest analysis. Second, we outline two plausible synchronization teleconnection\npathways beyond climate mediation: an atmospheric route (e.g., ENSO-driven savanna burning and\nlong-range smoke transport to adjacent forests), and a market-mediated route (e.g., ENSO-induced\ncommodity price shifts accelerating deforestation). This addition provides a clearer mechanistic\nexplanation, as requested.\nThe revised paper now reads:\n(Page 10 Lines 302-320, Discussion)\nENSO events exert asymmetric teleconnection effects on vegetation resilience, with El Niño–\nassociated regions exhibiting more pronounced declines (Fig. 2). Specifically, over one-third of\nvegetated areas experience significant resilience changes during El Niño events, compared to only\nabout 20% during La Niña (Fig. 3). Mechanistically, these effects are primarily mediated by\nregional climate: ENSO-driven anomalies in temperature, precipitation, and radiation trigger\nextremes—heatwaves, droughts, floods—can impair photosynthesis, water regulation, and root\nfunction32,33, ultimately slow canopy recovery. Our Random Forest analysis confirms that such\nanomalies, particularly those linked to El Niño (high temperature and radiation, water stress) and\nLa Niña (cold and dry conditions), are significantly associated with resilience decline globally (Fig.\n2). Even intact, high-biomass forests are vulnerable: our targeted analysis shows that resilience loss\nin these systems under ENSO-related stresses (Extended Data Fig. 4), consistent with previous\nfindings that tall canopies and deep-root systems depend heavily on stable water supply and become\nvulnerable under persistent ENSO pressure34. In addition to climate-mediated effects, ENSO may\ninfluence ecosystems through other pathways. For example, ENSO-driven circulation anomalies\ncan intensify savanna burning27 and facilitate long-range smoke transport to adjacent forests35,36,\nreducing incoming solar radiation and limiting canopy productivity37. ENSO-modulated shifts in\ncommodity prices (for example, palm oil)38,39, may also accelerate tropical deforestation and\nfragmentation40, further weakening forest resilience. These hypothesized pathways remain\nspeculative and warrant further investigation.\nvariable rainfall. Nature Geoscience 12, 174-179, doi:10.1038/s41561-019-0312-z (2019).\nterrestrial ecosystems to climate variability. Nature 531, 229-232, doi:10.1038/nature16986 (2016).\nof the National Academy of Sciences 114, 4442-4446, doi:10.1073/pnas.1617988114 (2017).\nGlob. Change Biol. 28, 571-587, doi:https://doi.org/10.1111/gcb.15939 (2022).\nforests. Ecohydrology 12, e2056, doi:https://doi.org/10.1002/eco.2056 (2019).\nIOD and the mechanisms for impacts on Australian rainfall. Journal of Climate 24, 3910-3923\n(2011).\nGlobal Ecology and Biogeography 29, 2039-2051, doi:https://doi.org/10.1111/geb.13174 (2020).\nCommun. 14, 498, doi:10.1038/s41467-023-36207-7 (2023).\nEvolution 7, 1790-1798, doi:10.1038/s41559-023-02187-6 (2023).\n(LAI): Methods, Products, Validation, and Applications. Reviews of Geophysics 57, 739-799,\ndoi:https://doi.org/10.1029/2018RG000608 (2019).\nProceedings of the National Academy of Sciences 121, e2316924121,\ndoi:doi:10.1073/pnas.2316924121 (2024).\nstudies. International Journal of Digital Earth 6, 5-33 (2013).\nMatthew J. Menne, Thomas M. Smith, Russell S. Vose, and Huai-Min Zhang. NOAA Extended\nReconstructed Sea Surface Temperature (ERSST), Version 5. NOAA National Centers for\nEnvironmental Information, doi:doi:10.7289/V5T72FNM (2017).\nexperimental design and organization. Geoscientific Model Development 9, 1937-1958 (2016).\nCMG V061. NASA Land Processes Distributed Active Archive Center,\ndoi:https://doi.org/10.5067/MODIS/MCD12C1.061 (2022).\n2000 to 2013. Science advances 3, e1600821 (2017).\nDistributed Active Archive Center,\ndoi:https://doi.org/10.5067/MEASURES/SRTM/SRTMGL1.003 (2013).\nGeneral Assembly (2016).\nWIREs Climate Change 15, e861, doi:https://doi.org/10.1002/wcc.861 (2024).\nAtmosphere Interaction, Geophys. Monogr 147, 21-48 (2004).\ndestabilize plant species coexistence via changes in plant–soil feedback. Nature Ecology &\nEvolution 6, 546-554, doi:10.1038/s41559-022-01700-7 (2022).\nResilience: Inferring the Chance of a Critical Transition. Annual Review of Ecology, Evolution,\nand Systematics 46, 145-167, doi:10.1146/annurev-ecolsys-112414-054242 (2015).\nOverturning Circulation. Nat. Clim. Chang. 11, 680-688 (2021).\n27\nSyst. Dynam. 14, 173-183, doi:10.5194/esd-14-173-2023 (2023).\nNature Ecology & Evolution 7, 1799-1808, doi:10.1038/s41559-023-02194-7 (2023).\nindicators of observation-based fingerprints of the Atlantic Overturning Circulation. Nat. Commun.\n14, 8344, doi:10.1038/s41467-023-44046-9 (2023).\ndoi:10.1038/s41467-023-36052-8 (2023).\nRes. Lett. 51, e2024GL108176, doi:https://doi.org/10.1029/2024GL108176 (2024).\nand leaf fall in a mixed temperate forest. Sci. Total Environ. 850, 157751,\ndoi:https://doi.org/10.1016/j.scitotenv.2022.157751 (2022).\natmospheric formaldehyde. Glob. Change Biol. 28, 1809-1822,\ndoi:https://doi.org/10.1111/gcb.15880 (2022).\nFor. Meteorol. 345, 109830, doi:https://doi.org/10.1016/j.agrformet.2023.109830 (2024).\nplants survive while others succumb to drought? New Phytologist 178, 719-739,\ndoi:https://doi.org/10.1111/j.1469-8137.2008.02436.x (2008).\nto drought: Exploring resilience and resistance. CATENA 243, 108156,\ndoi:https://doi.org/10.1016/j.catena.2024.108156 (2024).\nclimate change risks for forests. Forest Ecology and Management 259, 660-684,\ndoi:https://doi.org/10.1016/j.foreco.2009.09.001 (2010).\nCape Verde and the Amazon rainforest. Geophys. Res. Lett. 36,\ndoi:https://doi.org/10.1029/2009GL037923 (2009).\nCommunications Earth & Environment 4, 154, doi:10.1038/s43247-023-00795-5 (2023).\n5413, doi:10.1038/s41467-018-07921-4 (2018).\nPredictability. American Journal of Agricultural Economics 100, 239-263,\ndoi:https://doi.org/10.1093/ajae/aax060 (2018).\n380, 1064-1069, doi:doi:10.1126/science.adf2983 (2023).\nefforts to reduce deforestation in Brazil. Nat. Commun. 13, 5476, doi:10.1038/s41467-022-33213-\nz (2022).\nVegetation in the West Pacific Region. Geophys. Res. Lett. 49, e2021GL096666,\ndoi:https://doi.org/10.1029/2021GL096666 (2022).\nunder greenhouse warming. npj Climate and Atmospheric Science 7, 81, doi:10.1038/s41612-024-\n00627-z (2024).\nCommun. 5, 3712, doi:10.1038/ncomms4712 (2014).\nNat. Commun. 13, 917, doi:10.1038/s41467-022-28490-7 (2022).\nforests. Nature Ecology & Evolution, doi:10.1038/s41559-025-02776-7 (2025).\ndoi:doi:10.1126/science.abo5003 (2023).\nAmazonian forests. Proceedings of the National Academy of Sciences 118, e2019377118,\ndoi:doi:10.1073/pnas.2019377118 (2021).\nNat. Clim. Chang. 13, 967-974, doi:10.1038/s41558-023-01776-4 (2023).\n28\nresilience variability in a changing climate\nWe greatly appreciate the reviewers’ valuable and constructive suggestions. We have carefully revised\nour manuscript based on their comments. We hope that the revisions adequately address any concerns.\nBelow, you will find our item-by-item responses to the comments. The following different fonts are\nused to facilitate the review.\n Contents underlined are the original comments of the Reviewers.\n Contents in normal Times New Roman font are our response to the Reviewers.\n Contents in blue and italic font are part of the revised manuscript; bold text indicates newly\nadded content in direct response to reviewer comments.\n Green text highlights where the revisions have been made in the manuscript.\nResponses to Reviewer #1\nComment 1. The revisions have made the manuscript easier to understand and show the relevance of\nthe research clearly. I only have some minor comments left:\nReply: We sincerely thank the reviewer for the positive feedback and constructive comments. We are\nglad that the revisions have improved the clarity and readability of the manuscript. We have addressed\nthe remaining minor comments in detail below.\nComment 2. - In Fig 3, the figure caption mentions the inset box in the top right corner of b, but the\ninset boxes have been removed in this version, not sure if that was supposed to happen, because it is\nnot mentioned in the revisions?\nReply: We thank the reviewer for noticing this oversight. The inset boxes were unintentionally\nremoved during figure revision, and we have now restored them in Fig. 3b to match the description in\nthe caption.\nThe revised paper now reads:\n(Pages 8 Lines 250-262, Results)\nFig 3. ENSO impact on global vegetation resilience. (a) Synchronization impact of ENSO on resilience.\n(b) Regional differences in AC trends between affected and unaffected regions during 1981–2018. The\ndashed line represents the regional average. The inset box in the top right corner magnifies the\ndifference in regional averages. Asterisks indicate statistical significance between groups based on\nboth the Whitney U-test and Student’s t-test: p < 0.05(*); p < 0.01(**); p < 0.001(***). (c-h) Climate-\nmediated teleconnection effects on vegetation resilience: (c-e) Impacts of ENSO-induced\nanomalies in temperature, solar radiation, and precipitation on vegetation resilience, respectively. (f)\nand (g) Combined impact of these climate anomalies during El Niño and La Niña on vegetation\nresilience, respectively; (h) Overall impact of ENSO-induced climate anomalies on vegetation\nresilience. Plus (+) and minus (-) signs indicate positive and negative impacts, respectively, with\nnumbers indicating impact intensity, calculated as the sum of individual ENSO-climate anomaly\nimpacts.\nComment 3. - Line 273: Change ‘is’ to ‘are’: ‘In parallel, ENSO-driven climate anomalies ARE\nprojected to strengthen …’\nReply: We thank the reviewer for catching this grammatical oversight. The sentence has been\ncorrected accordingly in the revised manuscript.\nThe revised paper now reads:\n(Pages 8 Lines 273-274, Results)\nIn parallel, ENSO-driven climate anomalies are projected to strengthen, affecting approximately 51-\n55% of global vegetation areas and marking a 6-10% increase compared to historical levels (Fig. 4c).\nComment 4. - Lines 306-308: This sentence is not grammatically correct I think: ‘ENSO-driven\nanomalies in temperature, precipitation, and radiation trigger extremes— heatwaves, droughts,\nfloods—can impair photosynthesis, water regulation, and root function ultimately slow canopy\nrecovery.\n- Lines 311-314: Same for this sentence: ‘our targeted analysis shows that resilience loss in these\nsystems under ENSO-related stresses (Extended Data Fig. 4), consistent with previous findings that\ntall canopies and deep-root systems depend heavily on stable water supply and become vulnerable\nunder persistent ENSO pressure.’\nReply: We appreciate the reviewer’s careful attention to these grammatical issues. Both sentences\nhave been revised for clarity and grammatical correctness.\nThe revised paper now reads:\n(Pages 10 Lines 306-307, Results)\nENSO-driven anomalies in temperature, precipitation, and radiation trigger extremes—such as\nheatwaves, droughts, floods—that can impair photosynthesis, water regulation, and root function1,2,\nultimately slow canopy recovery.\n(Pages 10 Lines 311-314, Results)\nour targeted analysis shows that resilience loss occurs in these systems under ENSO-related stresses\n(Extended Data Fig. 4), consistent with previous findings that tall canopies and deep-root systems\ndepend heavily on stable water supply and become vulnerable under persistent ENSO pressure3.\nResponses to Reviewer #2\nComment 1. The authors have addressed all my comments carefully. I feel the paper is much stronger\nand clear after revisions. One very minor comment:\nReply: We thank the reviewer for the encouraging comments and are pleased to hear that the revisions\nhave strengthened and clarified the manuscript. We have carefully addressed the remaining minor\ncomment as detailed below.\nComment 2.Line 136 - The term drought may not be appropriate here, as it implies a deviation from\nlong-term climatological norms. What you likely mean is low precipitation. Please revise for accuracy.\nReply: We agree that “low precipitation” is a more accurate description than “drought,” and have\nrevised the sentence accordingly in the Results section.\nThe revised paper now reads:\n(Pages 4 Lines 136-137, Results)\nNotably, at a global scale, low precipitation, high temperature, and elevated solar radiation had a\nstronger impact on vegetation resilience compared to high precipitation, low temperature, and low\nradiation (Fig. 2b).\nReference\nplants survive while others succumb to drought? New Phytologist 178, 719-739,\ndoi:https://doi.org/10.1111/j.1469-8137.2008.02436.x (2008).\nto drought: Exploring resilience and resistance. CATENA 243, 108156,\ndoi:https://doi.org/10.1016/j.catena.2024.108156 (2024).\nclimate change risks for forests. Forest Ecology and Management 259, 660-684,\ndoi:https://doi.org/10.1016/j.foreco.2009.09.001 (2010)."}]

Decision letter