# Increased drought effects on the phenology of autumn leaf senescence

Chaoyang Wu, Jie Peng, Philippe Ciais, Josep Peñuelas, Huanjiong Wang, Santiago Beguería, T. Andrew Black, Rachhpal Jassal, Xiaoyang Zhang, Wenping Yuan, et al.

► To cite this version:

Chaoyang Wu, Jie Peng, Philippe Ciais, Josep Peñuelas, Huanjiong Wang, et al.. Increased drought effects on the phenology of autumn leaf senescence. Nature Climate Change, 2022, 10.1038/s41558-022-01464-9. http://hal-03779021

HAL Id: hal-03779021

https://hal.science/hal-03779021v1

Submitted on 7 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Increased drought effects on the phenology of autumn leaf senescence

Chaoyang Wu ID 1,2 ⌜, Jie Peng ID 1,2,3 ⌜, Philippe Ciais ID 4, Josep Peñuelas ID 5,6, Huanjiong Wang ID 1,2, Santiago Beguería 7, T. Andrew Black 8, Rachhpal S. Jassal ID 8, Xiaoyang Zhang 9, Wenping Yuan 10, Eryuan Liang ID 11, Xiaoyue Wang 1,2, Hao Hua ID 1,2, Ronggao Liu ID 12, Weimin Ju ID 13, Yongshuo H. Fu 14 and Quansheng Ge ID 1,2 ⌜

Global warming delays the autumn date of foliar senescence (DFS) in recent decades, with positive implications for growing season length and therefore global carbon storage. However, warming-associated drought, leading to water limitation, may conversely stimulate earlier DFS. Using ground observations since 1940s and 34 years of satellite greenness data (1982–2015) over the Northern Hemisphere (>30°N), we show the increased impact of drought on DFS. Earlier DFS is linked to decreased precipitation under warming and weaker drought resistance associated with various plant functional traits. For example, iso-hydric plants with strict regulation of water status may drop leaves fast during droughts. We derive an improved set of phenol-ogy models based on this influence and project earlier DFS by the end of the century, particularly at high latitudes (>50°N). Our results limit uncertainties in the later end of plant growth with warming, aiding estimation of carbon uptake of terrestrial ecosystems.

Plant phenology—the timing of life events such as leaf unfolding and senescence—has been extensively used in ecosystem models to reconstruct historical carbon (C) uptake changes as well as to predict future variations<fnref n="1" />,². The accuracy of C flux monitoring relies strongly on the appropriate representation of complex mechanisms that regulate the interaction between plant phenology and climate<fnref n="1" />. It has been previously shown that global warming produces earlier spring leaf unfolding over mid-high latitude regions³, imposing negative feedbacks on climate⁴. There is also growing evidence that autumnal phenology (that is, the date of foliar senescence, DFS) plays an essential role in regulating the length of growing season⁵ and therefore the interannual variability of C uptake⁶. The process of leaf senescence is triggered by the accumulation of abscisic acid⁷ and accompanied by degradation of chlorophyll levels and other pigments (for example, β-carotene and lutein) and remobilization of nutrients⁸ and it could also be regulated by environmental changes, such as temperature, radiation and water availability⁹. Despite the general agreement that a higher autumn temperature delays the autumn DFS¹⁰, several studies have revealed that the explanatory power of the effect of temperature on DFS is low and that a warmer autumn may be also associated with an earlier DFS and decreased annual C uptake¹¹. Therefore, a deeper understanding of the mechanisms that drive the interannual variability of the DFS is urgently required.

Less attention has been paid to the role of drought in plant phenology compared to the effects of temperature. The effects

of drought on plant phenology require further investigation because water availability may contribute more than temperature to the interannual variability of phenology, especially in arid and semi-arid ecosystems¹². Furthermore, water shortage induced by drought may lead to leaf damage and cause substantial reductions in ecosystem productivity¹³,¹⁴. Therefore, investigating drought's role in the DFS and its temporal trend may provide important insights into the understanding of future phenological variations with climate change. Thus, we suggest that over the past several decades, the autumn DFS may have exhibited a non-stationary response to drought with global warming and frequent drought events and such effects may shorten the previously predicted prolonged DFS driven by continuous warming at the end of this century. To this end, we detected the changing response of DFS to drought using a moving window method, analysed the possible reasons and improved DFS modelling and prediction. Two sources of DFS observations, three independent water availability indicators and other climatic and plant physiological datasets were used here: (1) 71,478 ground phenological observations at 905 sites in Europe recorded since the 1940s (Supplementary Fig. 1), (2) 34 years of satellite greenness data over the Northern Hemisphere (>30°N) for the period 1982–2015 (the latest third generation of the normalized difference vegetation index, NDVI3g), (3) the ERA5-Land and TerraClimate soil moisture (SM), (4) the standardized precipitation–evapotranspiration index (SPEI) and (5) the TerraClimate vapour pressure deficit (VPD) representing the atmospheric moisture and other gridded

<footnote n="1">¹The Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographical Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing, China. ²University of the Chinese Academy of Sciences, Beijing, China. ³State Key Laboratory of Grassland Agro-Ecosystems, College of Ecology, Lanzhou University, Lanzhou, China. ⁴Laboratoire des Sciences du Climat et de l'Environnement, IPSL-LSCE CEA CNRS UVSQ, Gif sur Yvette, France. ⁵CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona, Spain. ⁶CREAF, Cerdanyola del Valles, Barcelona, Spain. ⁷Estación Experimental de Aula Dei, Consejo Superior de Investigaciones Científicas, Zaragoza, Spain. ⁸University of British Columbia, Vancouver, British Columbia, Canada. ⁹Geospatial Sciences Center of Excellence (GSCE), Department of Geography, South Dakota State University, Brookings, SD, USA. ¹⁰School of Atmospheric Sciences, Center for Monsoon and Environment Research, Sun Yat-Sen University, Guangzhou, China. ¹¹CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing, China. ¹²State Key Laboratory of Resources and Environmental Information System, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing, China. ¹³International Institute for Earth System Science, Nanjing University, Nanjing, China. ¹⁴College of Water Sciences, Beijing Normal University, Beijing, China. ⌜e-mail: wucy@igsnrr.ac.cn; pengjie.18b@igsnrr.ac.cn; geqs@igsnrr.ac.cn</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

<figure><img src="image_8.png" /></figure>

**Fig. 1** | Temporal changes in the sensitivity of the autumnal foliar senescence date (DFS) to drought. **a,c,e** Site DFS observations to SM (**a**), SPEI (**c**) and VPD (**e**). **b,d,f**, Satellite-derived DFS to SM (**b**), SPEI (**d**) and VPD (**f**) with the significant $T_{sen+}$ and $T_{sen-}$ groups. NS, non-significant trend. The horizontal line is the median, the cross is the mean and the upper and lower edges are the 75th and 25th percentiles, respectively. Red line indicates the linear regression of mean sensitivity.

variables (temperature, precipitation, water use efficiency (WUE), isohydricity, root length and so on; Methods).

## Results

To obtain an objective and successive temporal variation in drought effects, we used a moving window approach instead of splitting

the overall time series into two parts artificially. On the basis of the broadly covered satellite NDVI3g data and ERA5-Land SM at the first soil layer, we found that using a 15-year window size provided an optimum compromise between accuracy and sample size (Methods; Supplementary Fig. 2). For both ground observations and satellite-derived DFS data, within each moving window, we
<figure><img src="image_9.png" /><figcaption>**Fig. 2 | Underlying mechanisms for enhanced drought effects.** a–f, Mean annual temperature (a), total annual precipitation (b), isohydricity (c), maximum root length (d), VOD (e) and canopy height (f). g, h, The correlations between the DFS sensitivity to SM with WUE (g) and dew point temperature (T<sub>dew</sub>) (h) at the T<sub>sen</sub> + and T<sub>sen</sub> − areas. P, N and NS represent positive, negative and not significant correlations, respectively. **, Significant difference between the two groups using the analysis of variance at P < 0.05.</figcaption></figure>

determined the sensitivity of DFS to drought as the absolute regres-
sion coefficient of the preseason SM in a nonlinear regression that
considered effects from both spring phenology¹⁵ and spring–sum-
mer temperature¹¹ (Methods). The preseason length was defined
as the time duration when the mean SM over this period had the
maximum correlation with DFS (Methods; Supplementary Fig. 3).
Therefore, these sensitivity (dDFS/dSM) values varied across the
time series of moving windows and the slope of the ordinary least
square regression of these sensitivities against time was expressed
as the temporal trend in sensitivity, T<sub>sen</sub>. We defined an increased

impact of drought on DFS (T<sub>sen</sub>+) as those sensitivity values signifi-
cantly increased from the first moving window (for example, 1982-
1996) to the last one (for example, 2001–2015) (P < 0.05). Likewise,
we defined a decreased impact of drought on DFS (T<sub>sen</sub>−) as those
sensitivity values significantly decreased (Methods; Supplementary
Fig. 4). We also detected the time-varying impacts of both SPEI and
VPD on DFS using the same moving window method.

During the period 1958–2015 in Europe, we found that the
sensitivity of DFS to SM showed a significantly increasing trend
(R²=0.20, P=0.00) (Fig. 1a). The same analysis using 34 years of
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 3** | Accuracy comparisons with the enhanced drought effects. The four models were CDD, DM, SIAM and DMT, while their improved versions comprise CDDS, DMS, SIAMS and DMTS, respectively. **a–d**, The evaluation metrics were the proportion of observations showing significant correlations (*P* < 0.05) (**a**), Pearson *R* (**b**), RMSE (**c**) and KGE (**d**). The number shows the mean value and the error bar indicates half of the standard deviation.</figcaption></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

satellite greenness data and SM data provided comparable results with 34.0% of the territory showing $T_{sen+}$ versus 23.0% $T_{sen-}$ (Fig. 1b). The overall sensitivity increased from 2.95±3.09 to 3.33±3.34 d mm⁻¹ for the first and last moving window with an $R^2$ of 0.90 (*P*=0.00). Results for SM at different depths were also provided (Supplementary Fig. 5). In addition to SM, we also calculated the sensitivity of DFS to both SPEI and VPD. For site-level DFS data, using SPEI and VPD also presented increased trends with $R^2$ of 0.58 (*P*=0.00) and 0.39 (*P*=0.00), respectively (Fig. 1c,e). The analyses from satellite-observed DFS were in line with that from SM and proportions for $T_{sen+}$ and $T_{sen-}$ were 33.1% and 21.5% for SPEI (Fig. 1d). These proportions differed slightly for VPD (30.0% and 21.4%) (Fig. 1f). Detailed results based on the first-layer SM indicator for plant functional types are shown in Supplementary Fig. 6. In addition, the results of SM, SPEI and VPD showed overall agreements in their spatial distributions (Supplementary Fig. 7). We also performed analyses using different sizes of moving window and the results were consistent (Supplementary Fig. 8). Overall, two independent phenological datasets and three drought indicators suggest an increased impact of drought on DFS for the past several decades in the Northern Hemisphere.

Given the larger geographic coverage and better representation of the diverse climate conditions and plant functional types in the satellite dataset and considering the dominant role of SM as a dryness stress¹⁶, we used these two datasets to explore the underlying mechanisms of the increased impact of drought on DFS. We grouped regions with enhanced or weakened DFS sensitivity, $T_{sen+}$ and $T_{sen-}$, respectively and found that both experienced comparable

warming of 0.4 °C decade⁻¹ over the past 34 years (Fig. 2a). However, the precipitation trends were significantly different in areas with enhanced, compared with weakened, DFS sensitivity. Areas that exhibited $T_{sen+}$ had almost 13% mean lower precipitation increase (Fig. 2b). This shortage of moisture supply caused different strategies for balancing photosynthesis and water loss because the isohydricity (the degree to which plants regulate their water status and a lower value indicating stricter regulation) of the $T_{sen+}$ was significantly lower (9.6%) than that of the $T_{sen-}$ regions (Fig. 2c). The mean maximum root length, vegetation optical depth (VOD) and canopy height were also found to be lower for $T_{sen+}$ than that for $T_{sen-}$ regions. The substantially higher WUE in the $T_{sen+}$ group (Fig. 2g) may also serve as an expression of increasing water sensitivity of ecosystem productivity as a result of comparable heat supply but less water availability. In addition, we found that the adaption of vegetation to drought may also help to explain the higher $T_{sen+}$ because this group of vegetation tended to increase the dew point temperature ($T_{dew}$) by much earlier leaf drop to avoid high volume water loss (Fig. 2h).

We also derived an improved set of phenology models based on the increased influence of drought on DFS and compared them with the semi-empirical cooling degree days (CDD) model, the process-based model referred to as the Delpierre et al. ¹⁷ model (DM), the spring-influenced autumn phenology model (SIAM) and the DM modified by spring-summer temperature (DMT) (Methods). We found that the DFS could be more accurately predicted when including the temporally changing drought-DFS relationships in the models in terms of the proportions of significant
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

**Fig. 4** | DFS difference predicted by the SIAMS and SIAM (SIAMS minus SIAM) under SSP 1-2.6 and SSP 5-8.5 scenarios. **a-d**, The overall spatial distribution of the mean difference over the Northern Hemisphere (>30°N) (**a,b**) and mean differences and half of the standard deviation among main vegetation types (**c,d**). In **a** and **b**, the letter N indicates an earlier DFS predicted with the SIAMS compared to the SIAM, while P indicates the opposite situation. In **c** and **d**, the coloured number above each bar indicates the mean difference of predicted DFS using two models.

predictions, the correlation coefficients (R), the root mean square error (RMSE) and the Kling-Gupta efficiency (KGE) (Fig. 3). For instance, only 16.9% of the predicted DFS time series were significantly correlated with observations using the CDD model and this percentage increased to 27.4% when drought effects were included. Similarly, the SIAM model generated 30.5% significant estimates but it increased to 47.9% with the new algorithm. The results using R, RMSE and KGE were highly consistent, revealing the value of considering drought effects on DFS for the modelling of phenology. Detailed spatial results for these indicators are provided in Supplementary Figs. 9–12.

Because the SIAM and SIAMS showed the highest accuracy, they were further used to predict future DFS under two climate scenarios, the Shared Socioeconomic Pathways (SSPs) 1-2.6 and 5-8.5. By using projected future climate change under both scenarios during the last two decades of this century (Supplementary Fig. 13), we found that, for most vegetation types, the DFS predicted using

the improved algorithms was earlier in the year than that predicted with the original models (Fig. 4). Under the SSP 1-2.6 scenario, with lower climate forcing, the DFS predicted by these two models showed larger differences for savannas and grasslands, with the occurrence of DFS being earlier by -3.6 ± 17.8 d (mean ± 0.5 s.d.) and -3.4 ± 16.5 d using the SIAMS (Fig. 4c). A higher emission scenario (SSP 5-8.5) would enlarge these differences with a much earlier predicted DFS and the largest difference was found for savannas (-6.3 ± 18.0 d), followed by woody savannas (-4.1 ± 17.9 d) and grasslands (-4.9 ± 15.9 d) (Fig. 4d). Spatial distributions for both scenarios showed that earlier DFS were mostly concentrated at relatively high latitudes (Supplementary Fig. 14).

## Discussion

Temporal responses of plant phenology to global warming are characterized by nonlinear patterns expressed as attenuated effects of temperature on C fluxes in spring¹⁸ and a weakening sensitivity
of spring phenology to temperature¹⁹. We found evidence of the enhanced impact of drought on DFS that was probably attributable to changes in physiological traits associated with relative water shortage. While there were spatial disagreements in the regions of enhanced drought sensitivity among different indicators, probably due to the different perspectives of drought characterization¹³,¹⁶, our results suggest that the relative importance of limiting factors on autumn DFS can vary and reveal that ecosystems have complex adjustments and adaptations to the continuously changing climate.

Plants evolve systematic strategies to maximize resources such as water, CO₂ and nitrogen and adapt to changes in the availability of these resources²⁰. The relative importance of these resources is often determined by their quantity compared to the smallest requirement for plant growth²⁰,²¹. A typical evolutionary response of plants to a limiting resource is often characterized by increasing their sensitivity to resource availability. For example, water availability is the limiting factor in arid and semi-arid ecosystems and plants are more sensitive to changes in water availability than to changes in other factors (for example, the convergence of a maximum sensitivity of growth under severe drought among several functional types of plants²²). Similarly, other variables may become limiting for growth (for example, temperature for boreal ecosystems and solar radiation for tropical ecosystems) when water becomes more available and changes in growth consequently become larger with per unit change of these factors²⁰,²². In our study, both *T*<sub>sen</sub>+ and *T*<sub>sen</sub>− regions experienced comparable temperature rise (Fig. 2a), while the precipitation increase for the *T*<sub>sen</sub>+ area was 13% less than that for *T*<sub>sen</sub>− area (Fig. 2b). Thus, the significantly less precipitation supply could be an important trigger for the enhanced drought sensitivity of DFS.

Plants have evolved different levels of stomatal regulation of water status, from the strictest regulation (isohydry) to less regulation labelled as anisohydry and thus may develop contrasting phenological strategies during a drought²³,²⁴. Here, we found that those plants in our group with enhanced drought sensitivity are more isohydric and maintain a relatively stable and high water potential regardless of drought conditions. These isohydric plants probably develop a conservative phenological strategy by losing leaves under drought conditions to increase root:shoot allocation ratio for root growth²⁴. These processes need the efficient shut down of stomata through abscisic acid release, which is an important trigger of fast leaf drop⁷. In comparison, anisohydric species allow leaf water potential to decrease as soil water declines with drought and are thus relatively drought tolerant²³. Therefore, anisohydric plants are more inclined to continue photosynthesis because of less sensitivity of leaf shedding at the beginning of drought conditions. This consequently caused weakened phenological responses of DFS to drought. A much shallower root length²⁵, the lower VOD²⁶ and canopy height²⁷ may also suggest a higher sensitivity to drought because of, for example, rapid depletion of plant water storage (indicated by VOD and canopy height) under drought conditions and the less potential in absorbing water from deep soil for *T*<sub>sen</sub>+ group. Changes in WUE also serve as an important reason for variations in drought sensitivity. Photosynthesis involves the tradeoff between carbon uptake and the loss of water and WUE is a useful proxy of this tradeoff effect, given that increased WUE has been closely linked to drought conditions²⁸. Our results also showed a much higher proportion of increased WUE for the *T*<sub>sen</sub>+ group and more regions with decreased WUE for *T*<sub>sen</sub>− group (Fig. 2g). The higher proportion of increased *T*<sub>dew</sub> for *T*<sub>sen</sub>+ types adds weight to the adaption of phenological changes of plants because of reduced water loss by a faster leaf drop during drought²⁹.

Considering the temporal changes in water constraint, we projected earlier DFS by the end of this century using the improved SIAMS compared to its original version. A later DFS was found to contribute to an overall increase in annual C uptake by temperate forests² but may offset C uptake due to extended ecosystem respiration

in higher latitudes⁶. The delayed DFS with warming, as predicted by ecosystem models without including this changing effect of water availability, may overestimate C uptake of terrestrial ecosystems. In particular, a larger proportion of the areas with comparatively earlier DFS was observed at northern high latitudes (>50°N) (Fig. 4 and Supplementary Fig. 14). This presents an urgent need for the re-evaluation of C uptake status because these regions harbour large C stocks and have very high sensitivity to global warming³⁰. Given that future precipitation change is spatially heterogeneous and projected droughts could be more frequent and severe³¹, our analyses of the accelerated impact of drought on DFS over the past decades reveal the unrecognized processes of water deficiency in regulating plant growth with warming.

## Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-022-01464-9.

Received: 29 November 2021; Accepted: 29 July 2022;
Published online: 5 September 2022

## References

1. 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).

2. Keenan, T. F. et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. *Nat. Clim. Change* **4**, 598–604 (2014).

3. Piao, S. L. et al. Leaf onset in the northern hemisphere triggered by daytime temperature. *Nat. Commun.* **6**, 6911 (2015).

4. Penuelas, J., Rutishauser, T. & Filella, I. Phenology feedbacks on climate change. *Science* **324**, 887–888 (2009).

5. Garonna, I. et al. Strong contribution of autumn phenology to changes in satellite-derived growing season length estimates across Europe (1982–2011). *Glob. Change Biol.* **20**, 3457–3470 (2014).

6. Piao, S. L. et al. Net carbon dioxide losses of northern ecosystems in response to autumn warming. *Nature* **451**, 49–52 (2008).

7. Zhao, Y. et al. ABA receptor PYL9 promotes drought resistance and leaf senescence. *Proc. Natl Acad. Sci. USA* **113**, 1949–1954 (2016).

8. Keskitalo, J., Bergquist, G., Gardestrom, P. & Jansson, S. A cellular timetable of autumn senescence. *Plant Physiol.* **139**, 1635–1648 (2005).

9. Liu, Q. et al. Delayed autumn phenology in the Northern Hemisphere is related to change in both climate and spring phenology. *Glob. Change Biol.* **22**, 3702–3711 (2016).

10. Wu, C. Y. et al. Contrasting responses of autumn-leaf senescence to daytime and night-time warming. *Nat. Clim. Change* **8**, 1092–1096 (2018).

11. Zani, D., Crowther, T. W., Mo, L., Renner, S. S. & Zohner, C. M. Increased growing-season productivity drives earlier autumn leaf senescence in temperate trees. *Science* **370**, 1066–1071 (2020).

12. Zhang, Y., Parazoo, N. C., Williams, A. P., Zhou, S. & Gentine, P. Large and projected strengthening moisture limitation on end-of-season photosynthesis. *Proc. Natl Acad. Sci. USA* **117**, 9216–9222 (2020).

13. Grossiord, C. et al. Plant responses to rising vapor pressure deficit. *New Phytol.* **226**, 1550–1566 (2020).

14. Ciais, P. et al. Europe-wide reduction in primary productivity caused by the heat and drought in 2003. *Nature* **437**, 529–533 (2005).

15. Keenan, T. F. & Richardson, A. D. The timing of autumn senescence is affected by the timing of spring phenology: implications for predictive models. *Glob. Change Biol.* **21**, 2634–2641 (2015).

16. Liu, L. B. et al. Soil moisture dominates dryness stress on ecosystem production globally. *Nat. Commun.* **11**, 4892 (2020).

17. Delpierre, N. et al. Modelling interannual and spatial variability of leaf senescence for three deciduous tree species in France. *Agric. For. Meteorol.* **149**, 938–948 (2009).

18. Piao, S. L. et al. Weakening temperature control on the interannual variations of spring carbon uptake across northern lands. *Nat. Clim. Change* **7**, 359–363 (2017).

19. Fu, Y. S. H. et al. Declining global warming effects on the phenology of spring leaf unfolding. *Nature* **526**, 104–107 (2015).
20. Seastedt, T. R. & Knapp, A. K. Consequences of nonequilibrium resource availability across multiple time scales: the transient maxima hypothesis. *Am. Nat.* **141**, 621–633 (1993).

21. Korner, C. Paradigm shift in plant growth control. *Curr. Opin. Plant Biol.* **25**, 107–114 (2015).

22. Huxman, T. E. et al. Convergence across biomes to a common rain-use efficiency. *Nature* **429**, 651–654 (2004).

23. McDowell, N. et al. Mechanisms of plant survival and mortality during drought: why do some plants survive while others succumb to drought. *New Phytol.* **178**, 719–739 (2008).

24. Nolan, R. H. et al. Differences in osmotic adjustment, foliar abscisic acid dynamics, and stomatal regulation between an isohydric and anisohydric woody angiosperm during drought. *Plant Cell Environ.* **40**, 3122–3134 (2017).

25. Fan, Y., Miguez-Macho, G., Jobbágy, E. G., Jackson, R. B. & Otero-Casal, C. Hydrologic regulation of plant rooting depth. *Proc. Natl Acad. Sci. USA* **114**, 10572–10577 (2017).

26. Choat, B. et al. Triggers of tree mortality under drought. *Nature* **558**, 531–539 (2018).

27. Giardina, F. et al. Tall Amazonian forests are less sensitive to precipitation variability. *Nat. Geosci.* **11**, 405–409 (2018).

28. Kannenberg, S. A., Driscoll, A. W., Szejner, P., Anderegg, W. R. L. & Ehleringer, J. R. Rapid increases in shrubland and forest intrinsic water-use efficiency during an ongoing megadrought. *Proc. Natl Acad. Sci. USA* https://doi.org/10.1073/pnas.2118052118 (2021).

29. Liu, Q. et al. Extension of the growing season increases vegetation exposure to frost. *Nat. Commun.* https://doi.org/10.1038/s41467-017-02690-y (2018).

30. Schuur, E. A. G. et al. Climate change and the permafrost carbon feedback. *Nature* **520**, 171–179 (2015).

31. Samaniego, L. et al. Anthropogenic warming exacerbates European soil moisture droughts. *Nat. Clim. Change* **8**, 421–426 (2018).

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Springer Nature or its licensor holds exclusive rights to this article under a publishing agreement with the author(s) or other rightsholder(s); author self-archiving of the accepted manuscript version of this article is solely governed by the terms of such publishing agreement and applicable law.

© The Author(s), under exclusive licence to Springer Nature Limited 2022
Methods

**Site-level ground LUD and DFS observations.** Ground leaf unfolding date (LUD) and DFS observations across Europe were represented by the date of first visible leaf stalk (BBCH scale³²=11) and date of autumnal leaf colouring (50%) (BBCH=94), respectively, compiled by the Pan European Phenology project (PEP725) (http://www.pep725.eu/). Consecutive records for at least 30 years for a given site and species were selected after choosing a moving window size of 15 years on the basis of the satellite data. As a result, a total of 71,478 site-year records were available for 9 plant species at 905 sites during the period 1946–2015. Due to the data availability, analyses based on SM and VPD started from 1958 to 2015. A more detailed description is presented in Supplementary Table 2.

**Satellite-derived LUD and DFS data.** Satellite LUD and DFS for 1982–2015 in the Northern Hemisphere (>30°N) were determined using the 15-d NDVI time series from the third generation of the Global Inventory Modeling and Mapping Studies (GIMMS) from the Advanced Very High-Resolution Radiometer (AVHRR), that is, the GIMMS NDVI3g dataset. This dataset has been widely used to model phenology over the Northern terrestrial ecosystems. We reconstructed the NDVI time series by reducing the snow and cloud impacts on the NDVI temporal trajectory before extracting LUD and DFS. First, we only focused on the vegetated pixels with clear seasonal variations and excluded those non-vegetated and agricultural lands. The pixels were selected on the basis of the multiple standards developed by ref. ³³. This also helps to solve the concerns on the landcover change issue because pixels showing low NDVI value (<0.1) even in a single year were removed from later analyses. Second, due to the absence of snow information, potential snow cover was eliminated using daily temperature data. For a given pixel, the snow-contaminated phase was identified when the daily mean air temperature was <0°C for at least 5 continuous days. The dormancy period was defined as the snow-contaminated phases during the period from the autumn to the spring of the next year. The snow-influenced NDVs during year *i* were replaced by the average of five maximum NDVs for two consecutive dormant phases (during the second half of year *i*−1, year *i* and the first half of year *i*+1)³⁴. Third, we used the modified Savitzky–Golay method, an effective and widely used tool for minimizing atmospheric contamination, to smooth the NDVI data³⁵. Fourth, we determined LUD and DFS using the following three algorithms: (1) the dynamic-threshold method³⁶, (2) a piecewise logistic function³⁷ and (3) a modified double-logistic function³⁸.

For the dynamic-threshold method, we linearly interpolated the 15-d smoothed NDVI into daily values generating the NDVI<sub>ratio</sub>, the LUD and DFS were set as the dates when the NDVI<sub>ratio</sub> became 0.5 in spring and autumn, respectively:

$$ \text{NDVI}_{\text{ratio}} = (\text{NDVI}_{t} - \text{NDVI}_{\text{min}}) / (\text{NDVI}_{\text{max}} - \text{NDVI}_{\text{min}}) \quad (1) $$

where NDVI<sub>t</sub> denotes the daily NDVI on day of year (DOY) *t*, NDVI<sub>max</sub> is the peak growth value during the summer and NDVI<sub>min</sub> is the minimum value during the year.

For the piecewise logistic function, the smoothed NDVI time series was divided into two sections before and after its annual peak value (a). Both parts were fitted separately into daily curves using the logistic function with different sets of fitting parameters:

$$ \text{NDVI}_{t} = \begin{cases} \frac{c_1}{1+e^{-c_1+t+\beta_1}} + d_1, & t \le \alpha \\ \frac{c_2}{1+e^{-c_2+t+\beta_2}} + d_2, & t > \alpha \end{cases} \quad (2) $$

where parameters *d*<sub>1</sub> and *d*<sub>2</sub> are background NDVI values, *c*<sub>1</sub> and *c*<sub>2</sub> are the differences between the peak and the background values and *a* and *b* are fitting parameters.

For the third method, the smoothed 15-d NDVI was fitted with the modified double-logistic function, which produced NDVI<sub>t</sub>:

$$ \text{NDVI}_{t} = \alpha_{1} + \frac{\alpha_{2}}{1+e^{-\beta_{1}(t-\beta_{1})}} - \frac{\alpha_{3}}{1+e^{-\beta_{2}(t-\beta_{2})}} \quad (3) $$

where *α*<sub>1</sub> is the background NDVI value, *α*<sub>2</sub> is the difference between the peak in summer and the lowest value in spring and *α*<sub>3</sub> is the difference between the peak in summer and the minimum NDVI in late autumn. *∂*<sub>1</sub> and *∂*<sub>2</sub> are the curvature parameters for the ascending and descending phases, and *β*<sub>1</sub> and *β*<sub>2</sub> are the mean date of spring leaf out and autumn senescence, respectively.

The LUD and DFS for the last two logistic functions were subsequently determined as the DOY at the local maximum and minimum of the first derivative of the NDVI curves, respectively. Further analyses were based on the averages of the values obtained by using the three methods.

**Water availability indicators.** We used three individual water availability indicators, including SM, SPEI and VPD to test the changing response of DFS to drought.

For the SM data, three sets of SM data were used to satisfy the needs of different time spans of the site DFS, satellite DFS and projected data. (1) The monthly ERA5-Land SM data (m³ m⁻³) for the period 1982–2015 was derived

from the European Center for Medium-Range Weather Forecast (ECMWF) ERA5-Land reanalysis dataset³⁹, which was produced by incorporating model data together with global observations using the physics laws. It started from 1981 to the present with a grid size of 9 km and a vertical resolution of four soil layers: 0–7 cm, 7–28 cm, 28–100 cm and 100–289 cm. (2) The monthly TerraClimate SM data⁴⁰ (mm) for 1958–2015 with a spatial resolution of 1/24° were obtained to match the longer duration of ground DFS observations. It was synthesized by combining the WorldClim with the Climatic Research Unit Time-series v.4.0 (CRU_TS v.4.0) and the Japanese 55-yr reanalysis using the climatologically aided interpolation. (3) To predict future DFS, the mass of water in the upper soil layer of 0–10 cm (kg m⁻²) for 2081–2100 was obtained from the higher-resolution version of the Max Planck Institute Earth System Model (MPI-ESM1.2-HR) under the scenarios of SSP 1-2.6 and SSP 5-8.5 (ref. ⁴¹; https://esgf-node.llnl.gov/search/cmpip6/).

The SPEI comprises the standardized difference between precipitation and potential evapotranspiration (PET), that is, the climatic water balance, during a fixed period⁴². The SPEI data at any given month can be calculated at multiple timescales or aggregation times, so, the SPEI at the *n*-month time scale (*n*-month SPEI) provided an overall climatic water balance over *n* continuous months. The FAO-56 Penman-Monteith (PM) algorithm is used to calculate PET in SPEI⁴³ and it has been proved that SPEI outperformed other drought indices in detecting drought events, such as the standardized precipitation index (SPI) and the Palmer drought severity index (PDSI)⁴⁴,⁴⁵. Here, we adopted the timescales of 1–12 months to detect the effect of water status during the previous single year⁴⁶. The SPEI data used here were derived from the SPEIbase v.2.5 dataset (https://spei.csic.es/database.html), which was calculated using the recently released 4.03 version of CRU_TS¹⁷. It has a spatial resolution of 0.5°, a monthly time interval and covers the historical period of 1901–2015.

Monthly VPD data (hPa) from 1958 to 2015 were derived from the TerraClimate dataset⁴⁰ for both site and satellite-derived DFS analyses.

**Other climatic and plant physiological datasets.** Daily temperature product for 1982–2015 from the NASA Global Land Data Assimilation System (GLDAS)⁴⁸,⁴⁹ was used to identify the snow-covered pixels when extracting LUD and DFS from the GIMMS NDVI3g time series and was also used to build DFS models. These data have a spatial resolution of 0.25° and a 3-h temporal interval. To predict DFS under future scenarios, we also obtained daily temperature (2081–2100) from the MPI-ESM1.2-HR dataset under both SSP 1-2.6 and 5-8.5 scenarios.

Monthly temperature and precipitation data for 1982–2015 were acquired from the CRU_TS v.4.03 (ref. ⁴⁷). WUE was determined as the ratio of GPP to evapotranspiration. GPP was produced by a revised light use efficiency model with a temporal resolution of 8 d and a pixel size of 0.05° (ref. ⁵⁰). The monthly evapotranspiration data were generated on the basis of a modified process-based land surface evapotranspiration/heat fluxes algorithm, with a resolution of 8 km and spanning from 1982 to 2013³¹. Both data were obtained from Scientific Data (https://www.nature.com/sdata/). The monthly T<sub>dew</sub> data (1982–2015) were derived from the ERA5-Land dataset with a spatial resolution of 0.1°.

We acquired the iso/anisohydray data which were produced on the basis of the Ku-Band backscatter data from QuikSCAT⁵². These have a pixel size of 0.05° and cover the whole globe. Extreme isohyric behaviour is indicated by a value near zero, while a value >1 suggests extreme anisohyric behaviour. The Ku-Band VOD data for 1987–2015 at a pixel size of 0.25° were derived from the VOD Climate Archive (VODCA)⁵³. This dataset was produced using VOD observations from multiple sensors. The maximum root length at a grid of 0.0083° was obtained from the EartH2Observe²⁵ and the canopy height was obtained from the Global 1 km Forest Canopy Height dataset (https://webmap.ornl.gov/ogc/dataset.jsp?dg_id=10023_1).

**Calculation of T<sub>sen</sub>.** To detect the temporal changes in the sensitivity of DFS to drought, we adopted a moving window method. Given the 34 years of satellite observation in our study, the window size was set from 13 years to 23 years. For a fixed window size, we computed the drought sensitivity of DFS using a nonlinear regression method within each moving window. Take SM as an example:

$$ \text{DFS} = a(\text{SM}) + b(\text{LUD}) + c(\text{Tss}) + d(\text{SM} \times \text{LUD}) + e(\text{SM} \times \text{Tss}) + \varepsilon \quad (4) $$

In equation (4), DFS represents the ground- or satellite-derived DFS time series, SM and Tss are mean soil moisture over its preseason length and spring-summer mean temperature time series, respectively, and LUD represents the ground- or satellite-derived leaf unfolding dates. Values *a*, *b*, *c*, *d* and *e* are regression coefficients and *ε* is the residual error. For the preseason length of SM, we defined it as the period when the mean SM over this period had the absolute largest partial correlation with DFS during the period 1982–2015, controlling the effects of LUD and Tss. The preseason length ranged from the month of DFS to its preceding *n* month (n ≤ 12) with a step of one month. Before applying the nonlinear regression and partial correlation, all variables were linearly detrended during the whole study period (1982–2015) to better focus on the interannual fluctuations.

Thus, the sensitivity of DFS to drought was determined as the absolute value of the coefficient *a* and it represented the changes in the DFS in response to per mm change in the preseason SM.
Accordingly, the temporal trend of sensitivity, T<sub>sen</sub>, was extracted from the slope of the ordinary least square regression of the sensitivities for the series of moving windows. The significantly increasing and decreasing trends in T<sub>sen</sub> were denoted by T<sub>sen</sub> + and T<sub>sen</sub> - at P < 0.05, respectively, while NS indicates a non-significant trend.

**Window size selection.** For each window size (13–23 yr), we calculated the percentages of significant partial correlation (P < 0.05) between the detrended DFS and preseason SM, controlling LUD and Tss. We selected a window size of 15 yr to perform subsequent analyses, as it provided robust results while resulting in a sufficiently high number of windows as required for the analysis of T<sub>sen</sub> (Supplementary Fig. 2).

**Calculations of R, RMSE and KGE.** The metrics for the accuracy evaluation between the predicted and observed DFS included Pearson correlation (R), the proportion of observations showing significant R at P < 0.05, RMSE and KGE. KGE has high recognition in hydrological model calibration as it balances the correlation, bias and ratio of variances<sup>24</sup>. A KGE value closer to 1 indicates better model performance.

$$R = \frac{\sum_{i=1}^{n} (p_i - \mu_p)(o_i - \mu_o)}{\sqrt{\sum_{i=1}^{n} (p_i - \mu_p)^2} \sqrt{\sum_{i=1}^{n} (o_i - \mu_o)^2}} \quad (5)$$

$$\text{RMSE} = \sqrt{\frac{\sum_{i=1}^{n} (o_i - p_i)^2}{n}} \quad (6)$$

$$KGE = 1 - \sqrt{(R - 1)^2 - \left(\frac{\mu_p}{\mu_o} - 1\right)^2 - \left(\frac{\sigma_p}{\sigma_o} - 1\right)^2} \quad (7)$$

where p<sub>i</sub> and o<sub>i</sub> represent the predicted and observed DFS at year i, respectively; μ<sub>p</sub> and μ<sub>o</sub> are the mean values of the predicted and observed DFS over n years, respectively and the total number of years is n. σ<sub>p</sub> and σ<sub>o</sub> are their standard deviation values, respectively.

**Phenology models.** The four original DFS models included the CDD<sup>35</sup>, DM<sup>17</sup>, SIAM<sup>15</sup> and DMT<sup>11</sup> and their improved versions by adding SM were referred to as the CDDS, DMS, SIAMS and DMTS, respectively. Furthermore, we also adopted a growing degree day model (GDD) to predict future LUD for SIAM-related DFS projection.

The GDD model is a classic and widely accepted spring phenology model. The temperature threshold was zero and the only parameter was the accumulated temperature. The temperature started to accumulate when the daily temperature rose above zero and it stopped on the day of leaf unfolding.

The CDD model predicted the DFS using the accumulated cold temperature, since this parameter may be an important cue for leaf colouring or fall:

$$\text{CDD}_t = \sum_{t=t_0}^{t_y} \max(T_{\text{th}} - T_{(t)}, 0) \quad (8)$$

$$\text{DFS} = t_y, \text{ if CDD}_t \ge \text{CDD}_{\text{th}} \quad (9)$$

where CDD<sub>t</sub> represents the accumulated cooling degree days from day t<sub>0</sub> to day t and T<sub>(t)</sub> represents the daily mean temperature on day t. The whole process started from day t<sub>0</sub> when the daily temperature dropped below the critical temperature (T<sub>th</sub>) and it ended at day t<sub>y</sub> when the CDD<sub>t</sub> reached the acquired cooling degree (CDD<sub>th</sub>). Finally, the DFS was determined as the DOY t<sub>y</sub>. In our study, the T<sub>th</sub> was set from 0 to 25 °C with an interval of 5 °C.

DM was developed on the basis of the CDD model considering the impact of both temperature and photoperiod on leaf colouring:

$$S_{\text{sen}}(t) = S_{\text{sen}}(t-1) + R_{\text{sen}}(t) \quad (10)$$

$$R_{\text{sen}}(t) = \begin{cases} [T_b - T_{(t)}]^x \times f[P_{(t)}]^y, & \text{if } t \ge \text{DOY}_s \\ 0, & \text{if } t < \text{DOY}_s \end{cases} \quad (11)$$

$$f[P_{(t)}] = \frac{P_{(t)}}{P_s} \quad \text{or} \quad [P_{(t)}] = 1 - \frac{P_{(t)}}{P_s} \quad (12)$$

$$\text{DFS} = t_y, \text{ if } S_{\text{sen}}(t) \ge Y_{\text{crit}} \quad (13)$$

where S<sub>sen</sub>(t) and R<sub>sen</sub>(t) denote the status and rate of leaf colouring on day t, respectively. T<sub>(t)</sub> is the daily mean temperature on day t, T<sub>b</sub> is the base temperature, P<sub>(t)</sub> is the photoperiod length on day t and P<sub>s</sub> is the critical photoperiod. The values of x and y could be 0, 1 or 2 but they cannot be zero at the same time. The process

of leaf colouring began on day DOY<sub>s</sub> when the daily temperature dropped to T<sub>b</sub> and the photoperiod decreased to P<sub>s</sub> and it ended at day t<sub>y</sub> when the state of leaf colouring, S<sub>sen</sub>(t), reached its critical value (Y<sub>crit</sub>). Finally, the DFS was determined as the DOY of t<sub>y</sub>.

SIAM and DMT are modified versions on the basis of DM. For SIAM, the Y<sub>crit</sub> value was linearly related to LUD anomaly (equation (14)); while it was linearly linked to spring-summer temperature for DMT (equation (15)).

$$Y_{\text{crit}} = a + b \times \text{LUD}_a \quad (14)$$

$$Y_{\text{crit}} = a + b \times \text{Tss} \quad (15)$$

LUD<sub>a</sub> in equation (14) is the LUD anomaly (deviation from the multiyear mean LUD) and Tss in equation (15) is the spring-summer mean temperature.

The improved models referred to as CDDS, DMS, SIAMS and DMTS were developed by incorporating the SM data into CDD, DM, SIAM and DMT, respectively. In these new models, the CDD<sub>th</sub> (for CDDS) or Y<sub>crit</sub> (for DMS, SIAMS and DMTS) was nonlinearly correlated with preseason SM. For example, the DMS was formulated by the following equation:

$$Y_{\text{crit}} = a + b \times \text{SM}_{\text{pre}} + c \times \text{SM}_{\text{pre}}^2 \quad (16)$$

where SM<sub>pre</sub> is the preseason SM at the month of DOY<sub>s</sub>, a, b and c are coefficients of the quadratic regression.

The 34 years of satellite-derived DFS time series at each pixel were used to parametrize our models. To reduce uncertainties caused by small sample size, we used fivefold cross validation to determine model parameters. The set of optimal parameters was selected when the KGE value between the modelled and observed DFS was the highest.

## Data availability

All data used in this study are included in the article and the Supplementary Information. The specific link for each dataset can be found in Supplementary Table 1.

## Code availability

All the data analyses and modelling were performed using MATLAB. The codes for T<sub>sen</sub> calculation and the eight phenology models as well as the satellite LUD and DFS data used in our study are available at https://doi.org/10.5281/zenodo.6892387. Other codes and relevant data are available upon request to the corresponding authors.

## References

32. Templ, B. et al. Pan European Phenological database (PEP725): a single point of access for European data. *Int. J. Biometeorol.* **62**, 1109–1113 (2018).

33. Shen, M. et al. Increasing altitudinal gradient of spring vegetation phenology during the last decade on the Qinghai-Tibetan Plateau. *Agric. For. Meteorol.* **189**, 71–80 (2014).

34. Zhang, X. Y. Reconstruction of a complete global time series of daily vegetation index trajectory from long-term AVHRR data. *Remote Sens. Environ.* **156**, 457–472 (2015).

35. Chen, J. et al. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. *Remote Sens. Environ.* **91**, 332–344 (2004).

36. White, M. A. et al. Intercomparison, interpretation, and assessment of spring phenology in North America estimated from remote sensing for 1982–2006. *Glob. Change Biol.* **15**, 2335–2359 (2009).

37. Zhang, X. et al. Monitoring vegetation phenology using MODIS. *Remote Sens. Environ.* **84**, 471–475 (2003).

38. Gonsamo, A., Chen, J. M., Price, D. T., Kurz, W. A. & Wu, C. Y. Land surface phenology from optical satellite measurement and CO<sub>2</sub> eddy covariance technique. *J. Geophys. Res.* **117**, G03032 (2012).

39. Muñoz, S. ERA5-Land Monthly Averaged Data from 1981 to Present (C3S CDS, date accessed:10-8-2021); https://doi.org/10.24381/cds.68d2bb30

40. Abatzoglou, J. T., Dobrowski, S. Z., Parks, S. A. & Hegewisch, K. C. TerraClimate, a high-resolution global dataset of monthly climate and climatic water balance from 1958–2015. *Sci. Data* **5**, 170191 (2018).

41. Müller, W. A. et al. A Higher-resolution version of the Max Planck Institute Earth System Model (MPI-ESM1.2-HR). *J. Adv. Model. Earth Syst.* **10**, 1383–1413 (2018).

42. Vicente-Serrano, S. M. et al. Response of vegetation to drought time-scales across global land biomes. *Proc. Natl Acad. Sci. USA* **110**, 52–57 (2013).

43. Allen, R. G., Smith, M., Pereira, L. S. & Perrier, A. An update for the calculation of reference evapotranspiration. *ICID Bull.* **43**, 64–92 (1994).

44. Gampe, D. et al. Increasing impact of warm droughts on northern ecosystem productivity over recent decades. *Nat. Clim. Change* https://doi.org/10.1038/s41558-021-01112-8 (2021).
45. Sheffield, J., Wood, E. F. & Roderick, M. L. Little change in global drought over the past 60 years. *Nature* **491**, 435–438 (2012).

46. Peng, J., Wu, C. Y., Zhang, X. Y., Wang, X. Y. & Gonsamo, A. Satellite detection of cumulative and lagged effects of drought on autumn leaf senescence over the Northern Hemisphere. *Glob. Change Biol.* **25**, 2174–2188 (2019).

47. Harris, I., Jones, P. D., Osborn, T. J. & Lister, D. H. Updated high-resolution grids of monthly climatic observations—the CRU TS3.10 Dataset. *Int. J. Climatol.* **34**, 623–642 (2014).

48. Beaudoing, H., Rodell, M. & NASA/GSFC/HSL. GLDAS Noah Land Surface Model L4 3 Hourly 0.25 × 0.25 Degree Version 2.0 (GES DISC, 2015); https://doi.org/10.5067/342OHQM9AK6Q

49. Beaudoing, H., Rodell, M. & NASA/GSFC/HSL. GLDAS Noah Land Surface Model L4 3 Hourly 0.25 × 0.25 Degree Version 2.1 (GES DISC, 2016); https://doi.org/10.5067/E7TYRXPJKWOQ

50. Zheng, Y. et al. Improved estimate of global gross primary production for reproducing its long-term variation, 1982–2017. *Earth Syst. Sci. Data* **12**, 2725–2746 (2020).

51. Zhang, K. et al. Vegetation greening and climate change promote multidecadal rises of global land evapotranspiration. *Sci. Rep.* https://doi.org/10.1038/srep15956 (2015).

52. Li, Y. et al. Estimating global ecosystem isohydry/anisohydry using active and passive microwave satellite data. *J. Geophys. Res.* **122**, 3306–3321 (2017).

53. Moesinger, L. et al. The global long-term microwave Vegetation Optical Depth Climate Archive (VODCA). *Earth Syst. Sci. Data* **12**, 177–196 (2020).

54. Gupta, H. V., Kling, H., Yilmaz, K. K. & Martinez, G. F. Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modelling. *J. Hydrol.* **377**, 80–91 (2009).

55. Botta, A., Viovy, N., Ciais, P., Friedlingstein, P. & Monfray, P. A global prognostic scheme of leaf onset using satellite data. *Glob. Change Biol.* **6**, 709–725 (2000).

Acknowledgements

This work was funded by the Strategic Priority Research Programme of the Chinese Academy of Sciences (XDA19040103), the National Natural Science Foundation of China (42125101) and the CAS Interdisciplinary Innovation Team (JCTD-2020-05). J.P. and P.C. were funded by European Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P. J.P. was also financially supported by the Fundación Ramon Areces grant ELEMENTAL-CLIMATE, the Spanish Government grant PID2019-110521GB-I00 and the Catalan Government grant SGR 2017-1005.

Author contributions

C.W., J.P. and Q.G. designed the research. C.W. and J.P. wrote the first draft of the manuscript. J.P. performed data analyses and remote-sensing model simulations. X.W. and H.H. contributed to model simulation. P.C., J.P. and S.B. substantially revised the manuscript with intensive suggestions. H.W., A.B., R.J., X.Z., W.Y., E.L., R.L., W.J. and Y.F. contributed to the writing of the manuscript.

Competing interests

The authors declare no competing interests.

Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41558-022-01464-9.

**Correspondence and requests for materials** should be addressed to Chaoyang Wu, Jie Peng or Quansheng Ge.

**Peer review information** *Nature Climate Change* thanks Ainong Li, Constantin Zohner and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at www.nature.com/reprints.