Forests are a major terrestrial carbon sink, but the increasing frequency and intensity of climate-driven disturbances such as droughts, fires and biotic agent outbreaks is threatening carbon uptake and sequestration. Determining how climate-driven disturbances may alter the capacity of forest carbon sinks in a changing climate is crucial. Here we show that the sensitivity of gross primary productivity to subsequent water stress increased significantly after initial drought and fire disturbances in the conterminous United States. Insect outbreak events, however, did not have significant impacts. Hot and dry environments generally exhibited increased sensitivity. Estimated ecosystem productivity and terrestrial carbon uptake decreased markedly with future warming scenarios due to the increased sensitivity to water stress. Our results highlight that intensifying disturbance regimes are likely to further impact forest sustainability and carbon sequestration, increasing potential risks to future terrestrial carbon sinks and climate change mitigation.

Terrestrial vegetation stores 450 Pg carbon (C) and sequesters ∼1.9 PgC yr -1 (ref. 1) (net carbon sink), which counterbalances approximately 20% of global carbon emissions from fossil-fuel burning. As the primary driver of the terrestrial carbon sink, forests have a key role in regulating terrestrial ecosystems and the carbon cycle. Forests in the conterminous United States (CONUS) sequester 173 TgC yr -1 (ref. 2), and offset 9.7% of anthropogenic carbon emissions annually. Forest-based strategies to mitigate climate change, such as reforestation, improved forest management and avoidance of forest loss, have been proposed as potentially impactful 'nature-based climate solutions' alongside dramatic reductions in fossil-fuel emissions 3,4 . However, disturbances that are sensitive to climate, such as droughts, fires and insect outbreaks, decrease forest productivity, increase tree mortality and decrease carbon storage, at least at short timescales [5][6][7][8] . At regional scales and over longer timescales, changes in disturbance regimes (for example, more-frequent and/or more-severe disturbances) impair forest resistance (the capacity of the ecosystem to maintain its state and function 9 ) and increase the risk of decreasing long-term carbon storage, which is crucial for mitigating climate change 10 . The frequency and intensity of prevalent disturbances, particularly droughts, fires and those from biotic agents (for example, insect outbreaks), are projected to increase in response to global warming [11][12][13] and will probably play a pivotal role in future forest carbon sink. Identifying the changes and dynamics of forest gross primary production (GPP) and how forests respond to environmental stressors after disturbances is paramount in systematically managing terrestrial ecosystems and effectively mitigating climate change.

Climate-driven disturbances have both direct and indirect effects on forest GPP. Direct effects often involve a concomitant decrease in GPP during disturbances. For example, the severe heat and drought event in 2003 led to a 30% decrease in ecosystem GPP in Europe 14 , resulting in a strong anomalous net forest carbon source. However, these direct effects are typically short-lived for droughts, with forest GPP rebounding to pre-disturbance levels within a few months to a year 15 . Although, if trees die during fires and insect outbreaks, forest GPP recovery may take multiple years or longer. Indirect effects refer to changes in the sensitivity of forest GPP to climate stressors after disturbances, which can persist for several years. The sensitivity of forest Artic e https://doi.org/10.1038/s41558-024-02022-1 indirect effects of disturbances on forest GPP is currently lacking. It remains unknown whether forest GPP becomes more or less sensitive to water stress after disturbances. Quantifying the long-term changes in the sensitivity of GPP to water stress in response to disturbances is crucial for enhancing our understanding and modelling the impacts of climate change on forest carbon cycling in the twenty-first century.

We aimed to investigate whether the sensitivity of ecosystem GPP to water stress changes after severe droughts, fires and insect outbreaks. We examined the factors driving these changes and assessed their implications for carbon uptake. Leveraging long-term remotely sensed GPP data in CONUS, we performed regression analysis to GPP to water stress, particularly water availability, is a critical measure of response that indicates the 'resistance' of forests to environmental variability and their capacity to sequester carbon. High sensitivity (low resistance) to water stress often signifies a high vulnerability to water deficits and climatic extremes, frequently preceding an increase in forest mortality 16 . The sensitivity of forest productivity, including GPP, tree-ring width, basal area growth and greenness, to drought can be influenced by various factors, such as environmental conditions (soil, topography and climate), stand composition (species and age), plant functional traits (wood density and hydraulic traits) and human management 9,[17][18][19][20][21][22][23][24] . However, a comprehensive understanding of the not applicable. The distribution maps (4 km) for fires and insect outbreaks were aggregated to 20 km for visual display. f-h, The change in sensitivity in climate space (mean annual temperature (MAT) versus mean annual precipitation (MAP) after drought (f), fire (g) and insect outbreak (h); 1 °C × 100 mm grid).

https://doi.org/10.1038/s41558-024-02022-1

understand the response of plant productivity to variations in water stress, as indicated by widely used drought indices such as the Palmer drought severity index (PDSI) 25 and the standardized precipitationevapotranspiration index (SPEI) 26 . We calculated the sensitivity of forest GPP to water stress (referred to as 'drought sensitivity') and compared it before and after disturbances. Machine-learning models, specifically random forest regression, were used to uncover the drivers and potential mechanisms underlying changes in drought sensitivity. We sought to answer the following research questions. (1) How does drought sensitivity change after severe disturbances across CONUS? (2) How do changes in drought sensitivity vary across different land-cover and ecosystem types? (3) What are the major factors influencing changes in drought sensitivity? (4) How might the observed changes in drought sensitivity affect vegetation carbon uptake under future warming scenarios?

The drought sensitivity across CONUS changed significantly after severe droughts and fires. We illustrated how to calculate the change in GPP drought sensitivity with schematics (Fig. 1a,b). The sensitivity increased significantly after severe droughts and fires (Fig. 1c,d and Supplementary Table 1), where the means of the changes in sensitivity (∆k) were 3.80 ± 0.95 gC m -2 (mean ± s.e.m.; P = 0.0001, generalized least squares (GLS)) and 3.83 ± 0.73 gC m -2 (P = 0, GLS), respectively. Most pixels (59.12%) indicated increased sensitivity after severe droughts, but some pixels in eastern and northwestern CONUS manifested decreased sensitivity (Fig. 1c). Fewer pixels were available for analysing the effects of fires, but the increase in sensitivity after fires was still significant (Fig. 1d), with 58.64% of the available pixels indicating increased sensitivity. The sensitivity, however, did not change significantly after insect outbreaks, -0.79 ± 1.23 gC m -2 (P = 0.52, GLS) (Fig. 1e). Forests in the northwest had decreased sensitivity after insect outbreaks, whereas those in the Rocky Mountains had increased sensitivity. The results were similar when using SPEI to represent water stress (Extended Data Fig. 1a-c), where the sensitivity increased significantly across CONUS after severe droughts and fires, at 5.68 ± 2.06 gC m -2 (P = 0.0058, GLS) and 3.95 ± 1.73 gC m -2 (P = 0.023, GLS), respectively, and decreased significantly after insect outbreaks, -6.22 ± 2.73 gC m -2 (P = 0.022, GLS). We note as well that the patterns were robust when considering only pixels with significant GPP-PDSI relationships (Extended Data Fig. 2 and Supplementary Table 2). In summary, disturbances clearly altered GPP drought sensitivity, but the directions of the change in sensitivity diverged among disturbances and regions. The sensitivity increased significantly (∆k = 6.21 ± 1.06 gC m -2 , P = 0, GLS; Supplementary Table 3) in hot and dry regions (for example, temperature (T) > 10 °C and precipitation (P) < 1,000 mm) (Fig. 1f) and did not change in cold and wet regions (for example, T < 10 °C and P > 2,000 mm) after severe droughts. The sensitivity increased significantly in hot and dry regions after fires (∆k = 5.33 ± 0.72 gC m -2 , P = 0, GLS); there were almost no wet regions (only five pixels) (Fig. 1g) because fire was concentrated in dry regions. After insect outbreaks, the sensitivity did not change in hot and dry regions and decreased significantly in cold and wet regions (∆k = -28.49 ± 6.80 gC m -2 , P = 0.0002, GLS) (Fig. 1h). The intercept of the GPP-PDSI model decreased significantly after disturbances, where the means of the changes in the intercept (∆b) were all significantly lower than zero: -7.91 ± 3.51 gC m -2 (P = 0.024, GLS), -7.73 ± 0.90 gC m -2 (P = 0, GLS) and -7.80 ± 1.94 gC m -2 (P = 0.0001, GLS) (Extended Data Fig. 3).

The drought sensitivity of forests generally increased after severe disturbances, with some notable differences among forest types. The sensitivity increased after severe droughts for evergreen needleleaf, evergreen broadleaf and deciduous broadleaf forests (Extended Data Fig. 4, Fig. 2a and Supplementary Table 1), with the largest increase in deciduous broadleaf forests (5.47 ± 2.22 gC m -2 ; P = 0.014, GLS). Evergreen and deciduous broadleaf forests were more sensitive to water stress after fires, and evergreen needleleaf forests were less sensitive, at -0.73 ± 0.45 gC m -2 (Fig. 2b). The change in the sensitivity was significant only in deciduous broadleaf forests, at 6.67 ± 3.08 gC m -2 (P = 0.031, GLS). The drought sensitivity of forests did not change after insect outbreaks. The changes in the sensitivity of shrubland and grassland were large but not always significant. The changes in drought sensitivity were highly unlikely to be due to long-term trends because undisturbed regions had different trends in sensitivity during 1982-2018 compared with disturbed regions with the same land-cover type (Extended Data Fig. 5). These results confirmed that the sensitivity of forests was affected by severe disturbances and that the sensitivity tended to increase after disturbances. The results were comparable when using SPEI (Extended Data Fig. 1d-f

https://doi.org/10.1038/s41558-024-02022-1

and deciduous broadleaf forests exhibited significantly increased sensitivity, at 23.50 ± 8.14 gC m -2 (P = 0.0041, GLS) and 16.18 ± 7.04 gC m -2 (P = 0.022, GLS), after severe droughts.

The intercept of the GPP-PDSI model decreased in forests (Fig. 2 and Supplementary Table 1). The decreases were due to lower biomass and foliar area caused by disturbances, leading to decreased forest productivity. The directions of the changes in the intercept were not consistent for grassland and shrubland. It is very intriguing that the intercept of shrubland increased after all three disturbances, which might indicate that the direct effects (GPP decrease) of disturbances on shrubland can recover quickly. However, the intercept change is not related strictly to the absolute GPP change since we used detrended GPP. From the perspective of vegetation structure, shrubs have smaller leaf area, frequent resprouting behaviours and shorter canopy heights than forests, which might benefit the recovery of shrubland productivity after disturbances.

Land.cover Trend.LAI Mean.CO 2 Mean.LAI Trend.CO 2 GPP.Anomaly Mean.PDSI Trend.T Trend.PDSI Trend.SRAD Interaction.TP MAT Trend.SM a Importance (%) MAT (f), trend of SRAD (g) and mean annual PDSI (h). (i-l) Drivers of the change in the sensitivity of GPP to water stress after severe insect outbreak: importance of drivers (i) and random forest partial dependence of the change in sensitivity on the three most important drivers: mean annual SRAD (j), trend of CO 2 (k) and trend of SRAD (l). The solid black line is the average, and the shading shows the range (from minimum to maximum) of the partial dependence from 100 runs of random forest models. LAI, leaf area index.

Several climatic factors play a crucial role in driving the change in drought sensitivity. These factors include the trend of soil moisture, mean annual temperature, the interaction between temperature and reversed precipitation and downward surface short-wave radiation (SRAD). For drought, the random forest model explained 67% of the variation (R 2 = 0.67; Extended Data Fig. 6a) in the change in sensitivity, with the trend of soil moisture emerging as the most important driver (Fig. 3a). Mean annual temperature and the interaction between temperature and reversed precipitation ranked as the second and third most important drivers, respectively. The change in sensitivity increased with decreasing trends of soil moisture (Fig. 3b), indicating that lower soil moisture levels contributed to higher post-disturbance sensitivity. Similarly, the change in sensitivity increased with mean annual temperature (Fig. 3c), suggesting that regions experiencing higher temperatures were more likely to exhibit increased sensitivity after disturbances. The interaction between temperature and reversed precipitation (Fig. 3d) had similar effects to mean annual temperature, with higher values (indicating hot and dry regions) associated with increased drought sensitivity.

For fires and insect outbreaks, the random forest models explained 37% and 28% of the variations in sensitivity change, respectively (Extended Data Fig. 6b,c). Mean annual temperature emerged as the most important driver for fires (Fig. 3e), where regions with higher temperatures were associated with increases in sensitivity after the disturbance (Fig. 3f), particularly when mean annual temperature exceeded approximately 15 °C. Increasing, increasing downward surface short-wave radiation (particularly when the trend of SRAD is greater than zero) and a high PDSI were linked to an increase in sensitivity after fires. In the case of insect outbreaks, mean annual downward surface short-wave radiation emerged as the most important driver (Fig. 3i), with regions experiencing high radiation (for example, >180 W m -2 ) demonstrating increased sensitivity, whereas regions with low radiation exhibited decreased sensitivity after the disturbance (Fig. 3j). This result aligns with the distribution of sensitivity change after insect outbreaks (Fig. 1e), where northwest regions generally displayed decreased sensitivity and low levels of short-wave radiation. (c) and +2 °C warming (d). 'Annual' indicates mean annual total ∆C across CONUS, and the asterisk indicates the mean is significant at the 0.05 level (two-sided, both GLS models and t test). Multiple comparisons are not applicable. Panels c and d are aggregated to 20 km for visual display.

https://doi.org/10.1038/s41558-024-02022-1

Increasing trends in CO 2 and short-wave radiation also contributed to the increase in sensitivity after insect outbreaks (Fig. 3k,l).

The change in drought sensitivity of GPP has notable implications for carbon uptake and loss. We quantified carbon uptake responses to water stress and assessed the recovery time by utilizing pixels with a post-disturbance period of at least 16 years (Methods). It took approximately 5 years for the sensitivity to return to the pre-disturbance level after severe droughts and around 6 years after fires (Extended Data Fig. 7b,c). Owing to the relatively short post-disturbance time (9 years: 2010-2018), the recovery time for insect outbreaks was not calculated. We further estimated the potential change in carbon uptake resulting from the sensitivity change. As a first-order exploration, the change in carbon uptake was calculated as the difference in GPP when considering the changes in sensitivity and intercept (equation ( 2)). For drought, the mean annual change in carbon uptake over the 5 year period following the disturbance was -9.94 ± 2.99 TgC (Fig. 4a) across CONUS, indicating a reduction in carbon absorption. For fires, the mean annual change in carbon uptake over the 6 year period after the disturbance was -1.50 ± 0.03 TgC (Fig. 4c). To assess the effects of future warming, PDSI under the +2 °C warming scenario was used to estimate the change in carbon uptake in response to the sensitivity change. For drought, the change in carbon uptake under the +2 °C warming scenario was estimated to be -11.21 ± 2.42 TgC (Fig. 4b), indicating even lower carbon absorption than historical conditions. For fires, the change in carbon uptake was estimated to be -1.63 ± 0.04 TgC under the warming scenario (Fig. 4d).

Severe disturbances, such as severe droughts, fires and insect outbreaks, have substantial effects on forest ecosystems. They all have the potential to alter forest composition, leading to a shift towards early succession species 27 , and can cause physiological damage to surviving trees. Both effects influence the sensitivity of post-disturbance forest productivity to water availability. This study specifically examines the indirect effects of severe disturbances and highlights the changes in GPP drought sensitivity following these extreme events. We have observed that severe disturbances, such as severe droughts and fires, tend to increase the sensitivity of forest productivity to water availability. This is particularly notable in deciduous broadleaf forests, possibly due to drought legacy effects 28 and disturbance-induced damage, such as embolism and overheating. These factors make trees more susceptible to subsequent water stress. In addition, for deciduous trees, physiological traits such as shallow roots 29 and thin bark probably contribute to increased susceptibility to fire and drought damage. By contrast, for evergreen needleleaf forests, the sensitivity did not change and even decreased after fires and insect outbreaks. This phenomenon may be attributed to decreases in stand density following disturbances. In general, stand density has increased in many evergreen needleleaf forests in the western USA due to historical fire suppression activities 30 . Disturbances can relax overstocked conditions and reduce competition 31 for water. This result suggests that thinning holds the potential to alleviate water stress in certain conifer forests. In addition, gymnosperm-dominated forests, mainly needleleaf forests in the western USA, have shown notable shifts characterized by decreases in P50 (water potential at which 50% of conductivity is lost) and increases in the hydraulic safety margin (that is, the difference between P50 and the minimum water potential experienced) 32 in response to climate-driven mortality, making these forests more drought tolerant. Non-forested ecosystems, such as shrublands and grasslands, also exhibit heightened drought sensitivity following disturbances, especially after severe droughts and fires. This increased sensitivity may be attributed to their proximity to absolute biogeographic and climate thresholds. Shrubs and grasses predominantly thrive in arid regions characterized by high solar radiation and temperature coupled with low water availability. Moreover, dry regions, particularly shrublands and grasslands, reveal stronger correlations between GPP and PDSI (Extended Data Fig. 2b). Climate emerges as a crucial factor influencing sensitivity changes, with hot and dry regions experiencing increased sensitivity and cold and wet regions showing decreased sensitivity (Fig. 1). Previous research 28 also indicates that plants in arid regions exhibit stronger drought legacy effects compared with those in wet regions.

Droughts, fires and insect outbreaks have varied effects on ecosystems. In general, understory species such as grasses and herbs are sensitive to water availability 33 and usually senesce rapidly due to water deficits. During droughts, both the understory and overstory vegetation will be constrained. The absolute GPP values of understory grasses might recover in the next year; however, the GPP sensitivity to water (indirect effects) of the whole ecosystem, including both the understory and overstory vegetation, may not recover as quickly. The sensitivity of grasslands increased significantly (P = 0, GLS) after droughts (Fig. 2a), indicating that the GPP sensitivity of grasses may not necessarily recover at the same pace as the absolute GPP. For overstory trees, drought legacy effects can last for years 28 . A similar situation goes for fires, where both the understory grasses and overstory trees might be burned during fires. The absolute GPP of understory grasses can recover quickly; however, the GPP sensitivity might not (Fig. 2b). For overstory trees, fire damage, such as heat stress (for example, heat emboli) and biomass consumption, can cause a long-lasting sensitivity change. For insect outbreaks, understory grasses are not the target of widespread insects such as bark beetles 7 , and thus the impacts of insects on understory vegetation will be small. The unchanged (or even decreased) sensitivity after insect outbreaks could be due to the decreases in stand density and shifts in composition.

Across CONUS, the increased sensitivity of GPP to water stress leads to substantial decreases in carbon uptake after severe droughts and fires. These decreases are logical as photosynthesis and carbon uptake tend to decline more when sensitivity increases at the same level of water stress, and climate change is likely to bring more-frequent and severe droughts in many regions. Our estimates suggest that the decrease in carbon uptake can persist for approximately 5-6 years after the disturbances. This reduction in carbon uptake hampers the capacity of terrestrial ecosystems to sequester anthropogenic carbon emissions. Furthermore, our findings highlight the accelerated decline in carbon uptake under warming scenarios. Less carbon would be absorbed in response to increased sensitivity under warming scenarios compared with historical conditions. Many regions in CONUS are projected to experience hotter and drier conditions under future warming scenarios 13 , leading to increased water stress on plant growth. As a result, ecosystems' ability to absorb carbon would decline, exacerbated by increased water stress and vegetation sensitivity. Regions such as western North America face high risks of carbon loss and species loss due to climate change, as indicated by global assessments 34 .

Climate change is expected to increase the frequency and severity of disturbances in many regions. From the perspective of climate policy and management, relying on planting more trees to counterbalance carbon emissions from fossil-fuel burning can be challenging, particularly when ecosystem productivity is threatened by disturbances that are climate dependent. Our results highlight pervasive changes in the drought sensitivity of GPP in forests after disturbances and indicate meaningful effects on carbon uptake. Increased sensitivity increases the vulnerability of ecosystems to drought and could lead to substantial decreases in ecosystem carbon uptake under future warming scenarios. These long-term dynamics are important for evaluating the capacity of terrestrial ecosystems as carbon sinks for carbon management, nature-based climate solutions and net-zero pledges.

https://doi.org/10.1038/s41558-024-02022-1

We used GPP to indicate the photosynthetic capacity and productivity of forests. Four state-of-the-art long-term GPP datasets covering CONUS were used: the Numerical Terradynamic Simulation Group (NTSG) Landsat GPP dataset 35 , the Global Land Surface Satellite (GLASS) GPP dataset 36 , the revised Eddy Covariance-Light Use Efficiency (EC-LUE) model-derived GPP dataset 37 and the near-infrared reflectance of vegetation (NIRv)-based GPP dataset 38 . The Numerical Terradynamic Simulation Group Landsat GPP dataset (1986-2021) provided 16 day 30 m GPP data based on Landsat data and climatic variables across CONUS. (See Supplementary Information for details on these GPP products.) The four sets of GPP data were resampled to 4 km and aggregated to the annual level. Anomalies were calculated and detrended for each pixel for each GPP product. The four detrended anomalies were averaged for each pixel for all subsequent analysis to avoid inconsistencies and biases among the GPP products. The following analysis was based on the average anomalies of the four GPP products from 1982 to 2018.

We utilized two widely recognized drought indices, the PDSI 25,39 and the SPEI 26 , to quantify water stress and calculate changes in GPP drought sensitivity. The PDSI is a standardized metric derived from a two-layer soil-water balance model, where negative values indicate dry conditions and positive values indicate wet conditions. We obtained monthly historical PDSI data at a 4 km resolution covering the period between 1982 and 2018 from TerraClimate 40 , which provides climatic variables such as temperature, precipitation, soil moisture, potential evapotranspiration (PET) and downward surface SRAD. The PDSI data from TerraClimate were based on the Penman-Monteith equation-based PET. Monthly PDSI values were averaged to generate annual PDSI data. Severe drought disturbances were defined as annual PDSI values below -3 (refs. 25,41). Other thresholds such as -2 and -4 were also tested, and the results were comparable when using -3 and -4. The -2 threshold usually indicates moderate droughts, which have limited effects on ecosystems. The SPEI is a multi-scalar drought index that captures atmospheric water deficits by considering the difference between precipitation and PET. In our analysis, we used monthly precipitation and PET data from TerraClimate at a 4 km resolution to calculate the monthly 4 km SPEI12 (scale = 12 months) between 1982 and 2018 using the 'SPEI' package in R. SPEI12 was selected as we used annual data in this study. Monthly SPEI12 values were averaged to the annual level, and severe droughts were defined as SPEI12 values below -1.2 (ref. 9). Other thresholds such as -1 and -1.5 were also tried, and the results were comparable when using -1.2 and -1.5. The -1 threshold is usually related to moderate droughts, which have limited effects on ecosystems.

Maps of annual fire severity were obtained from the Monitoring Trends in Burn Severity (MTBS) database, which has provided fire maps for the USA since 1984. The 30 m resolution maps of annual fire severity between 1984 and 2018 were downloaded and aggregated to 4 km. Low-severity pixels were discarded because low-severity fires generally have very small effects on forests. The proportion of burned 30 m pixels (moderate-and high-severity pixels) within each 4 km grid was calculated when aggregating the 30 m pixels to 4 km data. A 4 km grid was treated as burned when the proportion was >10%. Plots from the Forest Inventory and Analysis (FIA) programme were used to identify damage caused by insect outbreaks. Agent codes of 10-19 for the FIA plots indicate that insects caused mortality. The rates of mortality of basal area due to insects were used to indicate insect damage. The mean annual rate of mortality for each plot during 2000-2009 was taken from ref. 11. We aggregated the plots to 4 km and calculated the average rate for each pixel. An average rate >0.3% was used to define insect outbreaks. The species of biotic agents (particularly insects here) were not provided by the FIA data.

Land-cover maps at a resolution of 500 m from the Terra and Aqua combined MODIS (Moderate Resolution Imaging Spectroradiometer) Land Cover Type (MCD12Q1) version 6 from 2001 to 2020 were used to identify different vegetation types. Land-cover type 5 was used, where croplands and non-vegetation (water, urban and barren) classes were removed. Pixels with changes in land cover in 2001-2020 were also removed. We aggregated the 500 m land-cover maps to 4 km to match the GPP data. The land-cover map for 2001 is shown in Extended Data Fig. 4.

The slope of a linear regression (GPP anomaly versus PDSI) was used to represent the sensitivity of photosynthesis to water availability because most pixels (∼70%) in CONUS presented a linear relationship between GPP anomaly and PDSI (Extended Data Fig. 2). Those pixels presenting nonlinear relationships were dominated by grasslands. Thus, although a minority of areas do not have linear GPP-PDSI relationships, our results are robust to including only significant linear relationship pixels (Extended Data Fig. 2d-f and Supplementary Table 2), and a linear model represents the most parsimonious, comparable and understandable model of drought sensitivity; thus, we used it for all analyses. We compared simple linear regression and multiple linear regression when deriving the sensitivity of GPP to water stress, where the sensitivity from the two models was significantly correlated (Extended Data Fig. 8). For droughts, the slope of the linear regression between pre-drought GPP anomalies and PDSI was calculated to indicate pre-drought sensitivity (k 0 ). Similarly, post-drought sensitivity (k 1 ) was calculated using post-drought GPP anomalies and PDSI. The difference between the post-drought and pre-drought slopes was treated as the change in sensitivity (∆k).

where b 0 and b 1 are intercepts. The change in sensitivity derived with SPEI12 was calculated in the same way shown in equation (1), where SPEI12 was used to indicate water stress. We required at least 8 years of data for regression when calculating the slopes. According to the histogram of drought return intervals in Extended Data Fig. 9a, eight was the first break point for drought return intervals, and the disturbance return intervals of most pixels are higher than 8 years. We also tried 6 years (the trough) and 10 years (the second break point) as the minimum for regression, and the corresponding changes in sensitivity, which are comparable to the change in sensitivity (∆k) in Fig. 2a, are shown in Extended Data Fig. 9c,d. A threshold of 8 years strikes a reasonable balance between sample size and stability. A lower threshold, such as 6 years, could make the results of linear regression unstable, whereas a higher threshold would decrease the sample size. For fire (Extended Data Fig. 9b), the return intervals of most pixels were greater than 10 years, indicating that any values lower than or equal to ten were acceptable. Therefore, using eight as the minimum number of years for regression was reasonable. The length of pre-disturbance and post-disturbance data could vary, and 8 years was the minimum we used for regression analysis. The first drought was checked for each pixel. The change in sensitivity was calculated when we had at least eight data points for a regression analysis before and after the first drought. Otherwise, the second drought would be checked. This pixel was discarded if all droughts were not qualified for calculating the change in sensitivity. Continuous droughts were treated as one drought when calculating the change in sensitivity. When there were two or more drought events suitable for calculating ∆k, the first one would be used. The same procedure was applied to fires, producing the change in sensitivity after a fire. For insect outbreaks, we calculated the slope between GPP anomalies and PDSI before 2000 to represent pre-outbreak sensitivity, and the slope after 2009 was treated as the Artic e https://doi.org/10.1038/s41558-024-02022-1 post-outbreak sensitivity. We resorted to this approximation due to the limited availability of standardized FIA plots, which were accessible only from 2000 onwards and repeated every 5 or 10 years. Although the strategy used to calculate the change in sensitivity for insect outbreaks may not be ideal, it was a reasonable approach given the data constraints. To validate our findings, we conducted a sensitivity analysis by separately analysing the four GPP anomalies and the PDSI data. The results (Extended Data Fig. 10) are comparable to those obtained using the average anomalies of the four GPP products. The sensitivity change of each land-cover type based on a single GPP product is presented in Supplementary Table 4.

To account for spatial autocorrelations in the changes in sensitivity (∆k) and intercept (∆b), we used the GLS model 42 and examined their significance. We utilized the gls function from the 'nlme' package in R. The exponential correlation structure, corExp, was selected to assess the significance of the mean ∆k and mean ∆b by fitting GLS models with the formulas '∆k ~ 1' and '∆b ~ 1' (regression with only the constant term), respectively. We experimented with five correlation structures: exponential, sphere, Gaussian, ratio and linear. The exponential correlation structure was always ranked in the top two correlation structures for minimizing Akaike information criterion. So to save time for computation (GLS models are extremely time-consuming when dealing with large samples), we decided to use the exponential correlation structure throughout the manuscript. Owing to computational intensity, we aggregated the ∆k and ∆b maps for drought disturbances using a factor of three (4 km to 12 km) during the calculation. If the P value (two-sided) was below 0.05 based on the GLS, the constant term (mean ∆k or mean ∆b) was considered significant.

We used random forest regression 43 , a widely used machine-learning model, to capture the change in sensitivity (equation ( 1)) and identify the primary drivers responsible for this change. Random forest regression, which utilizes an ensemble of decision trees, is well suited for analysing high-dimensional data, such as large sample sizes with numerous predictors and complex structure. To train the random forest model, we utilized the 'randomForest' package in R. The response variable was ∆k, and the predictors included various variables: T, P, soil moisture, downward surface SRAD, PDSI, leaf area index 44 and CO 2 concentration 45 . For each pixel, we calculated the long-term mean (for example, mean annual temperature) and trend (for example, the trend of temperature, which is the slope of linear regression between T and year) of these variables using data from 1982 to 2018. In addition, we incorporated the interaction between T and P as a predictor, which represented the mean quantiles of T and reversed P. First, we calculated the quantiles (0-1) of T, and high T had high quantiles. Second, we reversed P by multiplying by -1 (-1 × P) and calculated the quantiles of the reversed P. Finally, we averaged the quantiles of T and reversed P. The average quantiles can represent the interaction between T and P, where high quantiles indicate hot and dry conditions. GPP anomalies in the disturbed years and land-cover types were also included as extra predictors. Before training the models, highly correlated predictors were removed. For example, when mean annual precipitation and mean annual soil moisture exhibited a high correlation (|r| > 0.7), the predictor with a lower correlation with ∆k was eliminated. In the random forest models, we used 500 decision trees, and the number of splits was determined as the square root of the number of selected predictors (⌊√n⌋). The leaf area index data were obtained from the half-month 8 km Global Inventory Modeling and Mapping Studies LAI4g product (1982-2020) 44 , which were resampled (bilinear) to 4 km and averaged to the annual level. Similarly, the global monthly 1° historical CO 2 concentration data 45 from 1982 to 2013 were resampled (bilinear) to 4 km and averaged to the annual level. Separate random forest models were developed for each of the three disturbances, as the drivers varied among them. The resulting partial dependence plots generated by the random forest models illustrated the relationships between ∆k and the predictors. Predictor importance was assessed on the basis of the decrease in node purity, where the importance of a predictor was calculated as the decrease in node purity attributed to that predictor divided by the sum of the decrease in node purity for all predictors, multiplied by 100%. We ran the random forest models 100 times and averaged the responses (importance and partial dependence) of the 100 models (Fig. 3).

Recovery time was determined as the duration required for the sensitivity to return to its pre-disturbance level. Because the sensitivity could increase following disturbances, ecosystems may take years to recover their original states. To derive the recovery time, we selected pixels with a post-disturbance time of at least 16 years, ensuring a minimum of two independent 8 year periods to calculate sensitivity. A moving-window strategy, utilizing 8 year intervals, was used to compute post-disturbance sensitivity for each selected pixel. Within each moving window (1-8 years, 2-9 years, 3-10 years and so on), the sensitivity was determined as the slope of linear regression between GPP anomaly and PDSI. We considered ecosystems to have recovered when the post-disturbance sensitivity was equal to or lower than the pre-disturbance sensitivity (k 0 ). The recovery time was defined as the first year of the moving window when the post-disturbance sensitivity was ≤k 0 . The calculation of recovery time is illustrated in the schematic provided in Extended Data Fig. 7a, which exemplifies a recovery time of 5 years. If the sensitivity decreased after disturbances, the recovery time would be identified when the post-disturbance sensitivity was ≥k 0 . The average recovery times for severe droughts and fires were 4.8 and 5.6 years, respectively. However, the recovery time for insect outbreaks could not be determined due to the limited post-disturbance period of 9 years (2010-2018).

The change in carbon uptake was estimated on the basis of the change in GPP. The following equation represents the change in GPP (∆GPP) in response to the change in sensitivity after disturbances:

where ∆GPP indicates the change in productivity, and ∆k and ∆b represent the changes in sensitivity and intercept, respectively (as shown in Fig. 1 and Extended Data Fig. 3). Historical PDSI data after disturbances were used to derive the change in carbon uptake (∆C = ∆GPP). For severe droughts, the closest integer to 4.8 (5 years) was used, whereas for severe fires, the closest integer to 5.6 (6 years) was used. Each pixel was associated with five PDSI values after severe droughts and six PDSI values after severe fires, which were inputted into equation ( 2) to calculate the corresponding ∆C. For example, if a severe drought occurred in the tenth year, PDSI values in the eleventh, twelfth, thirteenth, fourteenth and fifteenth years would be used in equation ( 2), producing five ∆GPP values. The mean annual ∆C per pixel after severe droughts and fires is presented in Fig. 4a,c, indicating a reduction in carbon assimilation through photosynthesis due to disturbances.

In warming scenarios, we utilized 4 km monthly precipitation and PET data under +2 °C warming obtained from ref. 46 to generate PDSI values under the +2 °C scenario with the 'scPDSI' package in R. The precipitation and PET data under the +2 °C scenario corresponded to the nominal years of 1985-2015. Monthly PDSI values under the warming scenario were calculated and averaged to the annual level. The change in carbon uptake was estimated using the annual PDSI values under the warming scenario following the same approach outlined in equation ( 2), with ∆k and ∆b obtained from Fig. 1 and Extended Data Fig. 3, respectively. The mean annual ∆C under the warming scenario is depicted in Fig. 4.