ARTICLE https://doi.org/10.1038/s43247-023-00835-0 OPEN Autumn canopy senescence has slowed down with global warming since the 1980s in the Northern Hemisphere Yichen Zhang 1, Songbai Hong1, Qiang Liu 2 ✉, Chris Huntingford 3, Josep Peñuelas 4,5, Sergio Rossi6, Ranga B. Myneni7 & Shilong Piao 1,8 ✉ Climate change strongly impact vegetation phenology, with considerable potential to alter land-atmosphere carbon dioxide exchange and terrestrial carbon cycle. In contrast to well- 1234567890():,; studied spring leaf-out, the timing and magnitude of autumn senescence remains poorly understood. Here, we use monthly decreases in Normalized Difference Vegetation Index satellite retrievals and their trends to surrogate the speed of autumn senescence during 1982–2018 in the Northern Hemisphere (>30°N). We find that climate warming accelerated senescence in July, but this influence usually reversed in later summer and early autumn. Interestingly, summer greening causes canopy senescence to appear later compared to an advancing trend after eliminating the greening effect. This finding suggests that summer canopy greening may counteract the intrinsic changes in autumnal leaf senescence. Our analysis of autumn vegetation behavior provides reliable guidance for developing and para- meterizing land surface models that contain an interactive dynamic vegetation module for placement in coupled Earth System Models. 1 Sino-French Institute for Earth System Science, College of Urban and Environmental Sciences, Peking University, Beijing 100871, China. 2 PLECO (Plants and Ecosystems), Department of Biology, University of Antwerp, Universiteitsplein 1, 2610 Wilrijk, Belgium. 3 UK Centre for Ecology and Hydrology, Wallingford, Oxfordshire OX10 8BB, UK. 4 CREAF, Cerdanyola del Valles, Barcelona 08193 Catalonia, Spain. 5 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona 08193 Catalonia, Spain. 6 Laboratoire sur les écosystèmes terrestres boréaux Département des Sciences Fondamentales, Université du Québec à Chicoutimi, Chicoutimi, Quebec G7H 2B1, Canada. 7 Department of Earth and Environment, Boston University, 685 Commonwealth Avenue, Boston, MA US- 02215, USA. 8 State Key Laboratory of Tibetan Plateau Earth System, Resources and Environment, Institute of Tibetan Plateau Research, Chinese Academy of Sciences, Beijing 100085, China. ✉email: qiang.liu@uantwerpen.be; slpiao@pku.edu.cn COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv 1 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 A utumnal senescence is a critical component of the seaso- (see definitions in Supplementary Table 1). Then, we quantified nal evolution of terrestrial ecosystems, which regulates the the effects of increasing NDVImax (i.e., canopy greening) on length of the growing season and terrestrial carbon ΔNDVI, and consequently, the extracted EOS. Together with this uptake1,2. Yet, such event and its response to changing climate is analysis, we discussed the implications of our study on bench- not well represented in most Earth System Models (ESMs)3. marking phenological modules for future generations of Earth Impact of climate change on the trends of autumn senescence as System Models (ESMs). well as its role in terrestrial carbon, water, nutrient cycling, fit- ness, and distribution of tree species have been reported over the Results and Discussion past decades1,2,4–9. Many of them focused on the timing of spe- Changes in the speed of autumn senescence. Across the cific phenophase, e.g., leaf coloring and abscission1. However, the Northern Hemisphere (>30°N, NH, thereafter), divergent trends process of autumn senescence consists of a series of complex in ΔNDVI can be discerned among the senescence months, with metabolic events such as macromolecule degradation, decreased larger increasing trend found in July (1.8 × 10−4 yr−1) and Sep- photosynthesis, and most importantly, nutrient re- tember (2.1 × 10−4 yr−1), whereas decreasing trends were found absorption10–12. These events are co-dominated by plant inter- for August (−3.7 × 10−5 yr−1) and October (−1.7 × 10−4 yr−1; nal conditions (e.g., metabolic adjustments, genetic expressions), Fig. 2a). There is also an increasing trend in ΔNDVI for the entire near-surface meteorological factors (e.g., temperature, solar senescence period (July to October, JASO) (Fig. 2b), although radiation and wind), water availability and their interactions12–14. noting this may be expected due to rising NDVImax under pla- Despite the complex climatic responses of autumn phenology, it netary greening (Fig. 2c). Dramatic regional variations, however, remains unclear, however, which stage of the senescence process existed in ΔNDVI trends across the NH (Fig. 2d–i). Specifically, contributes more to changes in the timing of end of the growing we found ΔNDVI increased in northwestern North America in season (EOS). Moreover, previous studies pay too much attention July, which dominated the acceleration of autumnal senescence to individual phenophases (e.g., EOS), which may conceal other during the senescence period (Fig. 2d, h). In contrast, negative components of the senescence process, such as the variations in ΔNDVI trends were mainly distributed in northeastern North the environmental regulations during the multiple stages of the America and eastern Eurasia during August, which contributed to autumn senescence process15. These gaps in knowledge would a slowing of the speed of autumn senescence (Fig. 2e, h). In leave the reliable modelling of autumn senescence in doubt16, and September, contrasting changes in ΔNDVI were observed finally hamper the ability of ESMs to simulate vegetation between the high and middle latitudes of Eurasia (Fig. 2f), dynamics under future climatic scenarios17,18. resulting in opposite trends in the speed of autumn senescence Here, we provide a timely assessment of autumn senescence (e.g., acceleration in middle Eurasia but deceleration in higher speed through the monthly decreases in NDVI (ΔNDVI, see latitudes) (Fig. 2h). We also found that the speed of autumn schematic representation in Fig. 1), which is also critical for the senescence slowed down in the west and east of North America in extraction of EOS. We first investigate the regulations of climatic October, when vegetation was already in dormant phase across factors (temperature, soil moisture, insolation i.e., downward most high latitudes (Fig. 2g). shortwave radiation, and windspeed) on the speed of autumn Overall, ΔNDVI over the entire senescence period increased in senescence (i.e., the interannual trends in ΔNDVI). Given that 53.1% of the NH (significant in 17.9%, Fig. 2h), implying an such behavior may be confounded by rising CO2-induced sum- overall acceleration of autumn senescence would occur in these mer canopy greening (i.e., uptrends in annual maximum NDVI, areas. Although acceleration of senescence was found in more NDVImax), we additionally calculated the relative decrease of than half of the NH, the remaining NDVI (NDVI left in the NDVI (δNDVI = ΔNDVI/NDVImax) as a surrogate for ΔNDVI month of EOS, see definition in Fig. 1 and Supplementary Table 1) still exhibited positive trends in 69.4% of the NH (Fig. 2i). These findings are anticipated due mainly to the apparent summer greening (i.e., a prevalent increase in NDVImax, Supplementary Fig. 1). Indeed, when averaged across NH, NDVImax increased significantly at a rate of 8.4 × 10−4 yr−1 (p < 0.001) during the senescence period (Fig. 2c). Meanwhile, ΔNDVI over JASO general increased, but its magnitude (2.1 × 10−4 yr−1, p > 0.05) is weaker than that of NDVImax. As a result, the remaining NDVI had an overall positive trend (6.3 × 10−4 yr−1, p < 0.001). Our analysis, therefore, reveals that changes in the speed of autumn senescence (ΔNDVI) and the peak greenness (NDVImax) can affect vegetation greenness at the later stage of the growing season. Climatic regulation on the speed of autumn senescence. To investigate the climatic regulations on the speed of autumn senescence, we conducted partial correlation analyses to isolate Fig. 1 Schematic representation of NDVI variations during the autumn the influence of insolation, wind speed, soil moisture, and tem- senescence process. In this study, autumn senescence starts at Peak Of perature (Supplementary Fig. 2). The factor with the strongest Season (POS) when NDVI reaches its annual maximum (NDVImax), and partial correlation with ΔNDVI was defined as the dominant ends at the End Of Season (EOS). The monthly decrease in NDVI (ΔNDVI) driver of autumn senescence (see Methods). Again, the dominant is estimated as the NDVI difference between the current and the following factors that accelerate or slow down the speed of autumn months (in positive values), and the remaining NDVI is the NDVI left in the senescence were heterogeneous and discernible in their patterns EOS month (in positive values). The sum of monthly ΔNDVI and the over the NH (Supplementary Fig. 3). remaining NDVI is therefore equal to NDVImax. Detailed definitions of these Our noted strong spatial variations in ΔNDVI trends may be items can be found in Supplementary Table 1. explained by the divergent biological responses of plants to 2 COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 ARTICLE Fig. 2 Trends in ΔNDVI over individual months and the entire senescence period from 1982 to 2018. a Temporal trends of ΔNDVI for senescence months in the Northern Hemisphere (NH) (recalling that a positive value of ΔNDVI corresponds to a decrease in NDVI). b Temporal trends of ΔNDVI over the senescence period (from July to October, JASO) and remaining NDVI at the end of the autumn senescence. c Temporal trends of the annual maximum NDVI (NDVImax). The bars inserted in a indicate trends in ΔNDVI for each month, the whole senescence period, and the remaining NDVI. d–i Spatial patterns of ΔNDVI trends during 1982–2018 for senescence months (d–g), the entire senescence period (h), and remaining NDVI (i). Black dots in d–i indicate significant trends at 0.05 level, while the inserted stacked bars indicate the frequency distributions of ΔNDVI trends. climate change across different climatic regions. In fact, climatic extreme temperatures. On the other hand, it might significantly regulations on ΔNDVI trends (higher values indicate a greater slow down the speed of senescence in its later stages by senescence speed, Fig. 1) were mainly associated with local prolonging the exposure of vegetation to favorable summer hydrothermal conditions (Fig. 3). At the beginning of the autumn climate. senescence period (that is, July), ΔNDVI correlated positively Unlike such reversal of temperature regulation on the speed of with temperature, but negatively with soil moisture in colder autumn senescence, the effects of insolation, soil moisture and (mean annual temperature, MAT < 5 °C) and relatively drier windspeed remained almost consistent along the senescence (mean annual precipitation, MAP < 500 mm) regions (Fig. 3a, f). period. To be specific, insolation was identified as the positive This finding indicated a possible acceleration of senescence dominant factor of ΔNDVI, particularly over the warmer regions caused by warming, because climatic warming may in turn lead to (MAT > −5 °C) and the lower latitudes (Fig. 3). Because the severe soil water deficits. Meanwhile, insolation was identified as increasing insolation might raise the risk of leaf photo-oxidative a major promotor of the speed of autumn senescence in 33.2% of damage, which could promote the progress of leaf senescence19. NH (Fig. 3k). Interestingly, the overall temperature effects on As for soil moisture and wind speed, their influence is only ΔNDVI shifted from positive (July) to negative (August, apparent in a few sporadic areas. September, and October) as autumn senescence progresses. The While for the remaining NDVI (higher values indicate slower negative temperature dominance on ΔNDVI trends over the NH senescence speed, Fig. 1), it was generally positively related to increased from ~31% (August and September, Fig. 3l, m) to temperature and insolation, but negatively related to wind speed 46.9% (October, Fig. 3n). These results imply that a warming (Fig. 3e, j, o). The positive temperature effect was prevalently climate could accelerate autumn senescence in earlier stages of observed over the NH (Supplementary Figs. 2x, 3f), which the senescence period, probably due to the detrimental summer confirmed again that warming climate contribute to slow down COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv 3 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 Fig. 3 The dominant climatic factors of the trends in ΔNDVI and the remaining NDVI. a–e The distribution of dominant climatic factors was negatively correlated with ΔNDVI for the senescence months (a–d) and the remaining NDVI (e) in the climate space. f–j The distribution of dominant climatic factors was positively correlated with ΔNDVI for the senescence months (f–i) and the remaining NDVI (j) in the climate space. Each climate bin was defined by 0.5 °C intervals of mean annual temperature (MAT) and 20 mm intervals of mean annual precipitation (MAP). The factor that appeared most frequently for each bin was presented in these figures. k–o The percentages of pixels occupied by dominant climatic factors, which positively (left bars, add to 100%) and negatively (right bars, add to 100%) correlated with ΔNDVI for senescence months (k–n) and the remaining NDVI (o). Note that the positive correlation with ΔNDVI or the negative correlation with the remaining NDVI indicates that this factor can accelerate the speed of autumn senescence, and vice versa. the speed of senescence in late autumn. Followed by insolation, its (Fig. 3e, o) through its role in promoting leaf abscission and positive and dominant effects on the remaining NDVI were intensifying drought stress by increasing evapotranspiration7. The mainly distributed in North Siberia and North America (with influence of soil moisture expressed a contrasting regional MAT < 0 °C) (Supplementary Figs. 2e, 3e), indicating the negative difference (Supplementary Fig. 2q). More specifically, in terms light regulation on leaf senescence was stronger than that of of inland Eurasia and North America, the remaining NDVI was warming for these regions20. We also found that the wind speed positively regulated by soil moisture (Supplementary Fig. 2q), exerted negative controls on the remaining NDVI in 34.3% of NH indicating that increases in soil moisture can slow down the speed 4 COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 ARTICLE Fig. 4 Schematic representation of the influence of greening on the speed of monthly senescence and the detection of EOS. a For each senescence month, monthly senescence is proxied by the decrease in NDVI (ΔNDVI) or percentage of decrease in NDVI relative to NDVImax (δNDVI = ΔNDVI/ NDVImax). Peak of season (POS) greening would lead to an apparent increase in ΔNDVI even if δNDVI stays invariant. b In correspondence to ΔNDVI and δNDVI, end of season can be determined by fixed NDVI thresholds (EOSa) or fixed percentage of NDVImax (EOSp). If NDVImax were kept unchanged, these two methods would deliver the same EOS, whereas increases in NDVImax would results to a later EOSa than EOSp. of autumn senescence in these water-limited ecosystems. How- remaining NDVI may have cascaded effects on the extraction of ever, in high-latitude Siberia (>50°N), soil moisture was EOS. negatively correlated with the remaining NDVI (Supplementary We therefore further compared the temporal trends of EOSa Fig. 2q), mainly due to the anoxia conditions under high soil and EOSp (definitions in Fig. 4b, also see definitions in moisture21 and limited nutrient availability in these regions of Supplementary Table. 1) during the period 1982–2018. Across permafrost6. the NH, we found areas occupied with delayed EOSa (60.2%, significant in 22.2%, Fig. 5a) exceeds that of EOSp (51.1%, significant in 15.1%, Fig. 5b). As a result, an overall delayed trend The effect of summer greening on the extraction of autumn of EOSa (0.6 days dec−1) but advanced EOSp (−0.4 days dec−1) phenology. Parallel to climate change, the observed summer trend was presented in NH (Fig. 5a–c). This divergence summed greening (that is, an increase of 8.4 × 10−4 yr−1 in NDVImax, up to 1.0 days dec−1 NDVImax-induced delay in EOSa. Given that p < 0.001, Fig. 2c) could modulate the extracted autumn phe- such an advance in autumn senescence was supported by a meta- nology and its temporal trend. We hypothesize that, accompanied analysis based on field observations22, we reckon that EOSp may with summer greening, the increasing mature leaves in the peak be a better indicator of autumnal leaf phenology, which captures season would contribute to an increasing trend in ΔNDVI, thus the physiological changes of plants. EOSa, carrying the additional creating an illusion of faster senescence (Fig. 4a). Such effect effect of NDVImax, my instead, serve as a better indicator for could be greatly reduced by replacing monthly senescence canopy greenness browning. Our results thus infer that the (ΔNDVI) with the relative decrease in NDVI (δNDVI = ΔNDVI/ increase in the number of mature leaves in the canopy under NDVImax). Summer greening, in a similar way, would contribute summer greening may counteract the intrinsic changes in to delayed EOS calculated from fixed NDVI thresholds (EOSa) autumnal leaf senescence. than that from fixed percentage of NDVImax (EOSp), implying the Among the studied plant functional types (PFTs), we also influence of changing NDVImax in EOS detection (Fig. 4b, also see observed clear mismatch between trends in EOSa and EOSp definitions in Supplementary Table. 1). (Fig. 5d). In line with the NH average (a combination of all Notably, trends in δNDVI generally presented similar spatial PFTs), we found the opposite contributions of NDVImax (delay) patterns as that of ΔNDVI in each of the senescence months and EOSp (advance) to EOSa trends in grasslands, savannas and (Supplementary Fig. 4a–e). After separating the contribution of shrublands in boreal and arctic regions. However, the delaying trends in NDVImax and δNDVI to trends in ΔNDVI (Supple- effect of summer greening could hardly compensate for the large mentary Fig. 5a–e, see Methods), we found that the increase in advancement in EOSp, as these ecosystems generally locate in NDVImax accounted for ~40% of ΔNDVI trends in NDVI for resource-limited areas (Supplementary Fig. 7). Growth enhance- parts of Siberia and central Eurasia in August and September ment (greening) during the earlier season may prematurely (Supplementary Fig. 5b, c). Despite this, the overall contribution consume water and nitrogen, which would in turn speed up the of NDVImax to ΔNDVI trends was marginal for other regions and progress of vegetation senescence and cause earlier end of periods (e.g., only 1.7, 5.9, 5.0 and 2.6% for Jul, Aug, Sep and Oct, growing season23–25. Indeed, the largest advancement in EOSp respectively). Such influence was confirmed by the relatively small was observed in grasslands of central Eurasia, which were portion of significant correlations between trends in NDVImax characterized with severe water limitations and significant and ΔNDVI for these months (Supplementary Fig. 6a–d). vegetation greening (Supplementary Fig. 1). Besides, similar spatial patterns of dominant climatic factors The observed delay in EOSp for the woody species could be were also observed in δNDVI, further indicating their tight ascribed to the higher resource availability in their relatively connection (Supplementary Fig. 3). wetter and fertile habitants, and their ability to accumulate large We also found that the influence of increasing NDVImax on nutrient and carbohydrate reserves in autumn26,27. The increase trends of the remaining NDVI was apparent (36.7% on average, in NDVImax contribute to later EOSa than EOSp, albeit varied Supplementary Fig. 5f), which was in line with the widespread across these PFTs (Fig. 5d). We quantified the contribution of significant correlations between trends in NDVImax and the increasing NDVImax to the delayed EOSa through the residue remaining NDVI over the NH (Supplementary Fig. 6f). This leads proportion of trends after subtracting EOSp from EOSa (i.e., the to higher fractions of positive trend areas for the remaining NDVI proportion of orange bars within the transparent purple bars, last (69.4% with 30.1% significant, Fig. 2i) than its percentage-based section of methods). As a result, the increase in NDVImax counterpart (55.3% with 14.7% significant, Supplementary Fig. 4f). contributed to more than half of the delay in EOSa for woody Obviously, the impact of conspicuous summer greening on the savannas (69.2%) and temperate shrublands (88.5%). However, COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv 5 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 Fig. 5 Mismatch between end-of-season (EOS) determined by fixed thresholds of NDVI (EOSa) and fixed percentages of NDVImax (EOSp). a, b Spatial patterns of trends in EOSa and EOSp in 1982–2018. Black dots indicate a significant trend at the 0.05 level. The stacked inserted bars indicate the frequency distributions of trends in EOSa and EOSp. c Difference between trends in EOSa and EOSp (trends in EOSa minus trends in EOSp). d Contribution of trends in EOSp (green bars) and NDVImax (i.e., EOSa-EOSp, orange bars) to trends in EOSa (purple bars) for different vegetation types (see Methods). MF mixed forests, DBF deciduous broadleaf forests, DNF deciduous needleleaf forests. such influence is minor in forests, contributing to 20.3% and 9.5% carried by EOSa. Indeed, the prominent differences between γa in mixed forests and deciduous broadleaf forests, respectively. On and γp were discovered in grasslands and temperate shrublands the contrary, the decrease in NDVImax in deciduous needleleaf (Supplementary Fig. 9d), where EOSa trends were predominantly forests advanced EOSa by 29.3%, comparing to EOSp. contributed by changes in NDVImax (Fig. 5d). Remarkably, NDVImax-induced divergence between EOSa and EOSp trends exerted strong controls on their responses to climate change (Supplementary Figs. 2, 8). Despite the observed clear Implications for model benchmarking. The mismatch between positive temperature dominance for EOSa and EOSp (Supple- phenology metrics based on absolute (ΔNDVI, EOSa) and relative mentary Fig. 3f, r), the difference in the sensitivities of EOSa and (δNDVI, EOSp) changes in NDVI indicates the considerable effect of EOSp to temperature was examined (e.g., EOSa (γa) and EOSp summer greening on the extraction of autumn phenology. In other (γp), see Methods). The partial least squares regressions analysis word, the increase in the number of mature leaves in canopy under suggested that both EOSa and EOSp had positive temperature summer greening might confound the actual changes in the timing of sensitivities in most areas, while negative sensitivities were leaf senescence. Such an effect is most clearly shown by the mismatch concentrated in xeric grasslands of central Eurasia and USA between EOSa and EOSp. EOSa, derived from fixed thresholds of (Supplementary Fig. 9a, b). Note that EOSa presented larger NDVI, is a direct reflection of drop in canopy greenness, with a temperature sensitivities than EOSp in 69.5% of northern areas, delayed trend in EOSa indicating the prolonged canopy browning e.g., γa (2.1 days °C−1) was more than 50% higher than γp process. EOSp, derived from fixed percentages of NDVImax, instead, (1.3 days °C−1) when averaged over the NH (Supplementary better captures senescence of leaves under the overall changes in Fig. 9c), due mainly to the additional NDVImax greening trend canopy greenness. As shown in our results, the widely observed 6 COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 ARTICLE global greening (Supplementary Fig. 1) in this study and previous Methods studies28,29, has, indeed, counteracted the intrinsic advance of leaf Datasets. In this study, NDVI during the 1982–2018 period was extracted from senescence by an apparent delay in canopy browning. Despite that the third-generation product of NASA’s Global Inventory Modeling and Mapping Studies (GIMMS) group37. Various issues, such as orbital drift, calibration, viewing these two metrics may represent phenology at different ecological geometry and volcanic aerosols, have been corrected in this product38. This dataset levels (e.g., leaf vs. canopy), leaf phenology is more closely related to consists of fortnightly NDVI observations at a spatial resolution of 0.083° (~8 km). the time when plants enter winter dormancy and is a key regulator of The monthly NDVI data was estimated through the maximum composition of vegetation dynamics and terrestrial carbon sequestration2. The each month’s two records. In this analysis, we screened out pixels dominated by cropland and evergreen forests (inferred from the MODIS Landcover Product, absence of eliminating the summer greening effect in the extraction MCD12C1, IGBP classification)39, because their seasonal cycles are unclear. of autumn phenology (i.e., threshold-based methods: EOSa) may Consistent with the NDVI time span, the monthly mean 2-m surface temperature introduce bias to our understanding of the leaf senescence process was derived from CRU.TS.4.05 with 0.5° spatial resolution40. Monthly insolation and consequently, the vegetation-climate interactions and the ter- was extracted and summed from CRU-JRA v2.2 with 0.5° spatial and 6-h temporal resolution41. The monthly surface soil moisture was acquired from the C3S dataset restrial carbon cycle. Since the timing of leaf senescence is crucial for provided by the European Center for Medium-Range Weather Forecasts reconstructing historical plant carbon uptake and predicting its (ECMWF) with a 0.25° spatial resolution, which was then interpolated to future changes in Earth system models, it is essential to update 0.5° × 0.5°. The monthly wind speed was obtained from the ERA5 reanalysis data at phenology modules embedded in state-of-the-art land surface a spatial resolution of 0.25° 42 and was interpolated to 0.5° × 0.5°. schemes using accurate representations of these processes. This will enable a better simulation of the terrestrial carbon cycle30,31. Phenology extraction methods. Our study was conducted across north of 30°N However, leaf senescence so far is one of the most difficult (excluding the subtropical regions), where plants have clear seasonal variations and are capable to extract phenological metrics. The Autumn senescence period is processes to be parameterized in terrestrial ecosystem models, defined as the time window between POS (Supplementary Fig. 10) and EOS which is restricted by our incomplete understandings of its (Supplementary Fig. 11). The decrease in NDVI (ΔNDVI) during the entire complex mechanisms. To our knowledge, DGVMs coupled with senescence period was defined as the difference in NDVI between POS month and land surface models generally determine the onset of plant EOS month. Similarly, the monthly ΔNDVI was the NDVI difference between the current and the following month. The remaining NDVI i.e., NDVI at EOS month, dormancy by tracking the thresholds of modelled or prescribed refers to the remaining greenness in the final stage of autumn senescence. For the leaf area index30. For example, in the ORCHIDEE model, it is minor regions with later POS (e.g., later than July, only 6.2% of the study area), the defined as the date when the modelled leaf area index (LAI) drops calculation of ΔNDVI also starts from POS month, NDVI data prior to the POS below 0.232. Evidence from our study suggests that fixed month was excluded from the analysis. Besides monthly senescence, we further obtained threshold-based EOS (EOSa). thresholds, however, may not accurately reproduce the timings To avoid the influence of potential snow cover on EOS extraction, we first used the of autumn leaf phenology when the summer greening effect is not spline function to interpolate the monthly temperature to a daily basis, and then thoroughly considered. Other models, for example, the Ecosystem obtained the boundary of the thermal growing season by a sequence of five days Demography model version 2 (ED2), used the prescribed below 0 °C. Before EOS extraction, NDVI pixels outside this thermal growing satellite-derived EOS estimates to parameterize the date of leaf season were replaced by that of the temporally nearest snow-free date. We than determined EOSa according to the following steps (Polyfit-Maximum method43; drop during the year33. In this case, an accurate estimation of Supplementary Fig. 12): (1) calculate the averaged NDVI curves during the past 37- EOS, especially taking the effect of summer greening into year (1982–2018), (2) estimate the rates of changes in the average NDVI curves as: consideration, is of great importance for accurate model NDVIratio(t)=[NDVI(t + 1)-NDVI(t)]/NDVI(t), (3) determine the time T when parameterizations. Moreover, a large body of models simulates the minimum NDVIratio occurred, (4) use the NDVI at time T (i.e., NDVI(T)) as the threshold for EOSa extraction, (5) apply the threshold to determine EOSa for the timing of leaf shed by a function of air temperature and (or) each year. The Polyfit-maximum method have been successfully used in many of day length34–36. Although these models directly simulate leaf previous phenological studies44–46, and its outcomes are comparable with other phenology, caution is still needed for benchmarking models commonly used methods6. The priorly determined NDVI threshold also gives us against satellite-derived results, since phenology based on the chance to correspondently obtain a percentage-based threshold, which is then used for further analysis (see detection of EOSp in next paragraph). absolute changes in vegetation greenness may be not sufficient The extraction of ΔNDVI and EOSa is probably confounded by global greening, to characterize the trajectories of leaf senescence. especially in the context of the prominent increase in NDVImax (i.e., the effect of We also found that EOSa, carrying the additional greening summer canopy greening, Fig. 4 and Supplementary Fig. 1). To reduce such an trend at POS, generally exhibits higher temperature sensitivities impact, we employed another method based on the percentage of NDVImax. During each month and the entire senescence period, we calculated the percentage than EOSp. This discrepancy was in support of findings in Zhang of decrease in NDVI compared to the peak season NDVI as δNDVI = ΔNDVI/ et al. (2020), who reported that the response of greenness-based NDVImax. Similarly, the percentage of remaining NDVI relative to NDVImax was phenology to climate change was greater than that of calculated as the division between the remaining NDVI and NDVImax. The photosynthesis-based phenology. Therefore, if the apparent percentage-based EOS (EOSp) was extracted from almost the same steps as EOSa, sensitivity of satellite-based leaf phenology (with the absence of except that the percentage of NDVI(T) relative to NDVImax (NDVI(T)/NDVImax) were used as thresholds to detect yearly EOSP (step 4, Supplementary Fig. 12). The a summer greening effect in the extraction of phenology) is used differences in ΔNDVI and δNDVI, and EOSa and EOSp gives us chance to to parameterize leaf behavior in DGVMs, the future carbon sink investigate the effect of NDVImax on autumn phenology detection (see last section potential of terrestrial ecosystems might be overestimated. of methods). Our study provides a novel angle on the speed of autumn senescence and its response to climate change over the NH. We find Climatic regulations on the process of autumn senescence. Benefiting from that warming climate accelerates senescence in July, but this previous studies, we examined the effects of four climatic factors i.e., insolation, influence usually reversed in later summer and early autumn. wind speed, soil moisture, and temperature on the leaf senescence process6,7,14,47. Note that we did not take precipitation into consideration, due mainly to the Insolation generally accelerates senescence across the entire following aspects: (1) Soil moisture is extensively correlated with precipitation in senescence period. Our analysis also verifies a summer greening the study area (Supplementary Fig. 13) and (2) soil moisture is more closely and (increase in NDVImax)-induced delay in EOS, as canopy greening directly related to plant water uptake than precipitation5,48. The lagged effect of may counteract the intrinsic changes in leaf phenology, and lead to climate on vegetation dynamics was also considered for each of the phenological matrices (ΔNDVI, δNDVI, EOSa and EOSp) following previous studies6,8,15. higher temperature sensitivities for EOS carrying the additionally Firstly, simple correlations were used to determine the preseason length of each greening trend. Global warming and increasing atmospheric CO2 climatic factors, which was defined as the period when the mean value of factors have been reported to continually contribute to earth greening, having the largest absolute simple correlation with phenological metrics (Supple- especially in NH28,29. Our analysis thus allows us to reconsider the mentary Figs. 14, 15). The preseason length ranged from the month of phenology metrices to its preceding 3 months, with a monthly time step. The preseason was linkage between changes in apparent canopy greenness and changes defined for each of the four phenological metrics and climatic variables. We then in intrinsic leaf physiology, which could provide an important performed partial correlation analysis between each of the phenological metrices benchmark for those developing DGVMs for use in ESMs. and climatic factors during the preseason (Supplementary Figs. 2, 8). It enables us COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv 7 ARTICLE COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 to explore the relationship between plant phenology and single climatic factor 4. Chen, L. et al. Leaf senescence exhibits stronger climatic responses during while eliminating the effects of the remaining factors. The climatic factor with the warm than during cold autumns. Nat Clim Chang 10, 777–780 (2020). largest positive (lowest negative) partial correlation coefficient was defined as the 5. Wang, X., Wu, C., Liu, Y., Peñuelas, J. & Peng, J. Earlier leaf senescence dates positive (negative) dominant factor (Supplementary Fig. 3). are constrained by soil moisture. Glob Change Biol 29, 1557–1573 (2022). 6. Liu, Q. et al. Delayed autumn phenology in the Northern Hemisphere is Temperature sensitivity of EOS. Our partial correlation analysis suggested dis- related to change in both climate and spring phenology. Glob Change Biol 22, tinct temperature dominances on EOSa and EOSp, while similar patterns were not 3702–3711 (2016). clear for ΔNDVI and δNDVI. We then performed a partial least squares regression 7. Wu, C. et al. Widespread decline in winds delayed autumn foliar senescence between EOSa (also EOSp) and the four climatic factors, and assigned the coeffi- over high latitudes. Proceedings of the National Academy of Sciences 118, cient of temperature as the temperature sensitivity of EOSa (γa) and for EOSp (γb). https://doi.org/10.1073/pnas.2015821118 (2021). The coefficient was commonly used as a proxy of changes in phenological dates per 8. Wu, C. et al. Increased drought effects on the phenology of autumn leaf unit increase in temperature (i.e., temperature sensitivity)49–51. senescence. Nat Clim Chang 12, 943–949 (2022). 9. Yu, H., Zhou, G., Lv, X., He, Q. & Zhou, M. Environmental factors rather than productivity drive autumn leaf senescence: evidence from a grassland in situ Separation of the effects of NDVImax in ΔNDVI and EOSa. In this study, ΔNDVI simulation experiment. Agric For Meteorol 327, https://doi.org/10.1016/j. for each month of senescence is also expressed as: agrformet.2022.109221 (2022). ΔNDVI ¼ δNDVI ´ NDVImax ð1Þ 10. Estiarte, M. & Peñuelas, J. Alteration of the phenology of leaf senescence and This formula inherently indicates that ΔNDVI is jointly determined by δNDVI fall in winter deciduous species by climate change: effects on nutrient and NDVImax. Trends of ΔNDVI during 1982–2018, therefore, can be separated proficiency. Glob Change Biol 21, 1005–1017 (2015). into the trends of δNDVI and NDVImax based on the multiplication derivation rule 11. Keskitalo, J., Bergquist, G., Gardeström, P. & Jansson, S. A Cellular Timetable ((uv)’=u’v + v’u, u and v are variables). The following approximation was used: of Autumn Senescence. Plant Physiol 139, 1635–1648 (2005). 12. Lim, P. O., Kim, H. J. & Gil Nam, H. Leaf Senescence. Annu Rev Plant Biol 58, tr½ΔNDVI  tr½δNDVI ´ mean½NDVImax þ tr½NDVImax ´ mean½δNDVI ð2Þ 115–136 (2007). where tr[ΔNDVI], tr[δNDVI] and tr[NDVImax] are trends of ΔNDVI, δNDVI and 13. Deslauriers, A., Rossi, S. & Sevanto, S. Metabolic memory in the phenological NDVImax during the period 1982–2018, respectively. mean[NDVImax] and events of plants: looking beyond climatic factors. Tree Physiol 39, 1272–1276 mean[δNDVI] are the multi-year averages of NDVImax and δNDVI. tr[δNDVI] x (2019). mean[NDVImax] and tr[NDVImax] x mean[δNDVI] thus denote the contribution 14. Fracheboud, Y. et al. The Control of Autumn Senescence in European Aspen. of δNDVI trends and NDVImax trends to ΔNDVI trends, respectively. To test the Plant Physiol 149, 1982–1991 (2009). reliability of such approximation, we additionally checked whether the sum of the 15. Piao, S., Wang, J., Li, X., Xu, H. & Zhang, Y. Spatio‐temporal changes in the terms on the right-hand side of Eq. 2 (i.e., the partitioning approach) is equivalent speed of canopy development and senescence in temperate China. Glob to the left-hand side. The feasibility of such partitioning is ensured by a perfect 1:1 Change Biol 28, 7366–7375 (2022). distribution of values on both side of Eq. 2 (Supplementary Fig. 16). 16. Vitasse, Y., Bresson, C. C., Kremer, A., Michalet, R. & Delzon, S. Quantifying As for the EOS, once the thresholds of NDVI(T) and NDVI(T)/NDVImax were phenological plasticity to temperature in two temperate tree species. Funct defined from the multi-year averaged NDVI curve (i.e., step 1 to step 4; Ecol 24, 1211–1218 (2010). Supplementary Fig. 12) the determination of EOSa and EOSp for each year (step 5, 17. Anav, A. et al. Evaluation of Land Surface Models in Reproducing Satellite based on the prescribed thresholds) is less likely to be influenced by the peak times or Derived Leaf Area Index over the High-Latitude Northern Hemisphere. Part decrease speed of NDVI. For example, if NDVImax keeps constant (e.g., 0.8), II: Earth System Models. Remote Sens 5, 3637–3661 (2013). NDVI(T) for EOSa is 0.2 and the NDVI(T)/NDVImax for EOSp is 25%. No matter 18. Murray-Tortarolo, G. et al. Evaluation of Land Surface Models in Reproducing how the peak time and decrease speed changes, the time when NDVI value decrease Satellite-Derived LAI over the High-Latitude Northern Hemisphere. Part I: to 0.2 always coincides with the time when a relative decrease to 25% occurred Uncoupled DGVMs. Remote Sens 5, 4819–4838 (2013). (Supplementary Fig. 17). Therefore, we can conclude that if NDVImax stays the same, 19. Wu, Z. et al. Atmospheric brightening counteracts warming‐induced delays in EOSa and EOSp should be equal. Difference in EOSa and EOSp can thus be used to autumn phenology of temperate trees in Europe. Glob Ecol Biogeogr 30, indicate asynchronization of EOS resulted from the increasing NDVImax. 2477–2487 (2021). 20. Zhang, Y., Commane, R., Zhou, S., Williams, A. P. & Gentine, P. Light limitation regulates the response of autumn terrestrial carbon uptake to Data availability warming. Nat Clim Chang 10, 739–743 (2020). All observational data sets that we used are publicly available. The GIMMS3g NDVI data 21. Mainiero, R. & Kazda, M. Effects of Carex rostrata on soil oxygen in relation sets covering the period 1982–2015 are publicly available at https://data.tpdc.ac.cn/en/ to soil moisture. Plant Soil 270, 311–320 (2005). data/9775f2b4-7370-4e5e-a537-3482c9a83d88/, datafiles covering the period 2016–2018 22. Liu, H. et al. Phenological mismatches between above- and belowground plant were directly requested from J. Tucker, the creator of GIMMS NDVI3g dataset. The responses to climate warming. Nat Clim Chang 12, 97–102 (2021). CRU.TS.4.05 mean 2-m surface temperature data sets are available at https://data.ceda.ac. 23. Leuzinger, S., Zotz, G., Asshoff, R. & Korner, C. Responses of deciduous forest uk/badc/cru/data/cru_ts/cru_ts_4.05. The CRU.JRA v2.2 6-h insolation data sets are trees to severe drought in Central Europe. Tree Physiol 25, 641–650 (2005). available at https://data.ceda.ac.uk/badc/cru/data/cru_jra/cru_jra_2.2. The C3S soil 24. Lian, X. et al. Summer soil drying exacerbated by earlier spring greening of moisture data sets are available at https://cds.climate.copernicus.eu/cdsapp#!/dataset/ northern vegetation. Sci Adv 6, https://doi.org/10.1126/sciadv.aax0255 (2020). satellite-soil-moisture?tab=overview. The ERA5 reanalysis wind speed data sets are 25. Zhang, Y., Keenan, T. F. & Zhou, S. Exacerbated drought impacts on global available at https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-single- ecosystems due to structural overshoot. Nat Ecol Evol 5, 1490–1498 (2021). levels-monthly-means?tab=overview. 26. Dietze, M. C. et al. Nonstructural Carbon in Woody Plants. Annu Rev Plant Biol 65, 667–687 (2014). 27. Fu, Y. S. H. et al. Variation in leaf flushing date influences autumnal senescence and next year’s flushing date in two temperate tree species. Code availability Proceedings of the National Academy of Sciences 111, 7355–7360 (2014). All the analyses and figures are made using MATLAB R2020a and Power Point 2019, and 28. Piao, S. et al. Characteristics, drivers and feedbacks of global greening. Nat Rev the codes are available from the corresponding author. Earth Environ 1, 14–27 (2019). 29. Zhu, Z. et al. Greening of the Earth and its drivers. Nat Clim Chang 6, 791–795 (2016). Received: 24 November 2022; Accepted: 5 May 2023; 30. Richardson, A. D. et al. Terrestrial biosphere models need better representation of vegetation phenology: results from the North American Carbon Program Site Synthesis. Glob Change Biol 18, 566–584 (2012). 31. Keenan, T. F. et al. Net carbon uptake has increased through warming- induced changes in temperate forest phenology. Nat Clim Chang 4, 598–604 (2014). References 32. Krinner, G. et al. A dynamic global vegetation model for studies of the coupled 1. Piao, S. et al. Plant phenology and global climate change: Current progresses atmosphere-biosphere system. Glob Biogeochem Cycle 19, https://doi.org/10. and challenges. Glob Change Biol 25, 1922–1940 (2019). 1029/2003GB002199 (2005). 2. Piao, S. et al. Net carbon dioxide losses of northern ecosystems in response to 33. Medvigy, D., Wofsy, S. C., Munger, J. W., Hollinger, D. Y. & Moorcroft, P. R. autumn warming. Nature 451, 49–52 (2008). Mechanistic scaling of ecosystem function and dynamics in space and time: 3. Gallinat, A. S., Primack, R. B. & Wagner, D. L. Autumn, the neglected season Ecosystem Demography model version 2. J Geophys Res 114, https://doi.org/ in climate change research. Trends Ecol Evol 30, 169–176 (2015). 10.1029/2008JG000812 (2009). 8 COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-023-00835-0 ARTICLE 34. Caffarra, A., Donnelly, A. & Chuine, I. Modelling the timing of Betula Acknowledgements pubescens budburst. II. Integrating complex effects of photoperiod into This study was supported by the National Natural Science Foundation of China process-based models. Clim Res 46, 159–170 (2011). (41988101). 35. Dufrêne, E. et al. Modelling carbon and water cycles in a beech forest. Ecol Model 185, 407–436 (2005). 36. Jeong, S. J. & Medvigy, D. Macroscale prediction of autumn leaf coloration Author contributions S.L.P., Y.C.Z., S.B.H. and Q.L. designed the research; Y.C.Z. performed analysis; Y.C.Z. throughout the continental United States. Glob Ecol Biogeogr 23, 1245–1254 wrote the first draft of the manuscript. C.H., J.P., S.R. and R.B.M. contributed to the (2014). interpretation of the results and to the text. 37. Tucker, C. J. et al. An extended AVHRR 8‐km NDVI dataset compatible with MODIS and SPOT vegetation NDVI data. Int J Remote Sens 26, 4485–4498 (2010). 38. Kaufmann, R. K. et al. Effect of orbital drift and sensor changes on the time Competing interests series of AVHRR vegetation index data. IEEE Trans Geosci Remote Sensing 38, The authors declare no competing interests. 2584–2597 (2000). 39. Friedl, M. A. et al. MODIS Collection 5 global land cover: Algorithm refinements and characterization of new datasets. Remote Sens Environ 114, 168–182 (2010). Additional information 40. Harris, I., Osborn, T. J., Jones, P. & Lister, D. Version 4 of the CRU TS Supplementary information The online version contains supplementary material monthly high-resolution gridded multivariate climate dataset. Sci Data 7, available at https://doi.org/10.1038/s43247-023-00835-0. https://doi.org/10.1038/s41597-020-0453-3 (2020). 41. Kobayashi, S. et al. The JRA-55 Reanalysis: General Specifications and Basic Correspondence and requests for materials should be addressed to Qiang Liu or Shilong Characteristics. Journal of the Meteorological Society of Japan Ser II 93, 5–48 (2015). Piao. 42. Bell, B. et al. The ERA5 global reanalysis: Preliminary extension to 1950. Q J R Meteorol Soc 147, 4186–4227 (2021). Peer review information Communications Earth & Environment thanks the anonymous 43. Piao, S., Fang, J., Zhou, L., Ciais, P. & Zhu, B. Variations in satellite-derived reviewers for their contribution to the peer review of this work. Primary Handling Editor: phenology in China’s temperate vegetation. Glob Change Biol 12, 672–685 (2006). Aliénor Lavergne. 44. Jeong, S. J., Ho, C. H., Gim, H. J. & Brown, M. E. Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere Reprints and permission information is available at http://www.nature.com/reprints for the period 1982-2008. Glob Change Biol 17, 2385–2399 (2011). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 45. White, M. A. et al. Intercomparison, interpretation, and assessment of spring published maps and institutional affiliations. phenology in North America estimated from remote sensing for 1982-2006. Glob Change Biol 15, 2335–2359 (2009). 46. Cong, N. et al. Changes in satellite-derived spring vegetation green-up date and its linkage to climate in China from 1982 to 2010: a multimethod analysis. Open Access This article is licensed under a Creative Commons Glob Change Biol 19, 881–891 (2013). Attribution 4.0 International License, which permits use, sharing, 47. Gill, A. L. et al. Changes in autumn senescence in northern hemisphere adaptation, distribution and reproduction in any medium or format, as long as you give deciduous trees: a meta-analysis of autumn phenology studies. Ann Bot 116, appropriate credit to the original author(s) and the source, provide a link to the Creative 875–888 (2015). Commons license, and indicate if changes were made. The images or other third party 48. Liu, L. et al. Soil moisture dominates dryness stress on ecosystem production material in this article are included in the article’s Creative Commons license, unless globally. Nat Commun 11, https://doi.org/10.1038/s41467-020-18631-1 (2020). indicated otherwise in a credit line to the material. If material is not included in the 49. Fu, Y. H. et al. Declining global warming effects on the phenology of spring article’s Creative Commons license and your intended use is not permitted by statutory leaf unfolding. Nature 526, 104–107 (2015). regulation or exceeds the permitted use, you will need to obtain permission directly from 50. Menzel, A. et al. European phenological response to climate change matches the copyright holder. To view a copy of this license, visit http://creativecommons.org/ the warming pattern. Glob Change Biol 12, 1969–1976 (2006). licenses/by/4.0/. 51. Yang, Y., Guan, H., Shen, M., Liang, W. & Jiang, L. Changes in autumn vegetation dormancy onset date and the climate controls across temperate ecosystems in China from 1982 to 2010. Glob Change Biol 21, 652–665 (2015). © The Author(s) 2023 COMMUNICATIONS EARTH & ENVIRONMENT | (2023)4:173 | https://doi.org/10.1038/s43247-023-00835-0 | www.nature.com/commsenv 9