# Contrasting responses of autumn-leaf senescence to daytime and night-time warming

Chaoyang Wu, Xiaoyue Wang, Huanjiong Wang, Philippe Ciais, Josep Peñuelas, Ranga Myneni, Ankur Desai, Christopher Gough, Alemu Gonsamo, Andrew Black, et al.

► To cite this version:

Chaoyang Wu, Xiaoyue Wang, Huanjiong Wang, Philippe Ciais, Josep Peñuelas, et al. . Contrasting responses of autumn-leaf senescence to daytime and night-time warming. Nature Climate Change, 2018, 8 (12), pp.1092-1096. 10.1038/s41558-018-0346-z . hal-02975897

HAL Id: hal-02975897

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

Submitted on 29 Sep 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>
This is the accepted version of the article: Wu, C. et al. *Contrasting responses of autumn-leaf senescence to daytime and night-time warming* in *Nature climate change* (Ed. Nature publishing), vol. 8 (Dec. 2018), p. 1092-1096.

Available at: DOI 10.1038/s41558-018-0346-z

This version is published under a “All rights reserved” license.
1 **Contrasting responses of autumn leaf senescence to daytime**
2 and nighttime warming

3 Chaoyang Wu<sup>1,2*</sup>, Xiaoyue Wang<sup>1,2</sup>, Huanjiong Wang<sup>1,2*</sup>, Philippe Ciais<sup>3</sup>, Josep Peñuelas<sup>4,5</sup>, Ranga B. Myneni<sup>6</sup>, Ankur R. Desai<sup>7</sup>, Christopher M. Gough<sup>8</sup>, Alemu Gonsamo<sup>9</sup>, Andrew T. Black<sup>10</sup>, Rachhpal S. Jassal<sup>10</sup>, Weimin Ju<sup>11</sup>, Wenping Yuan<sup>12</sup>, Yongshuo Fu<sup>13</sup>, Miaogen Shen<sup>14</sup>, Shihua Li<sup>15</sup>, Ronggao Liu<sup>16</sup>, Jing M. Chen<sup>9</sup>, Quansheng Ge<sup>1,2*</sup>

1. The Key Laboratory of Land Surface Pattern and Simulation, Institute of Geographical Sciences and Natural Resources Research, Chinese Academy of Sciences, Beijing 100101, China;

2. University of the Chinese Academy of Sciences, Beijing 100049, China;

3. Laboratoire des Sciences du Climat et de l'Environnement, IPSL-LSCE CEA CNRS UVSQ, 91191, Gif sur Yvette, France;

4. CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, Barcelona 08193, Catalonia, Spain.

5. CREAF, Cerdanyola del Valles, Barcelona 08193, Catalonia, Spain;

6. Department of Earth and Environment, Boston University, 685 Commonwealth Avenue, Boston, MA 02215, USA;

7. Department of Atmospheric and Oceanic Sciences, University of Wisconsin-Madison, Madison, WI, USA;

8. Department of Biology, Virginia Commonwealth University, Richmond, VA 23284-2012, USA;

9. Department of Geography and Planning, University of Toronto, 100 St. George St., Toronto, ON, M5S 3G3, Canada;

10. Faculty of Land and Food Systems, University of British Columbia, Vancouver, BC, Canada;

11. International Institute for Earth System Science, Nanjing University, 22 Hankou Road, Nanjing, 210093, China;

12. School of Atmospheric Sciences, Center for Monsoon and Environment Research, Sun Yat-Sen University, Guangzhou 510275, China;

13. College of Water Sciences, Beijing Normal University, 100875, Beijing, China;

14. CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China;

15. School of Resources and Environment, University of Electronic Science and Technology of China, Chengdu 611731, China;

16. State key Laboratory of Resources and Environmental Information system, Institute of Geographic Sciences and Natural Resources Research, CAS, Beijing 100101, China

*Corresponding authors: C Wu (wucy@igsnrr.ac.cn), H Wang (wanghj@igsnrr.ac.cn) and Q Ge (geqs@igsnrr.ac.cn).
Plant phenology is a sensitive indicator of climate change¹⁻⁴, and plays a significant role in regulating carbon uptake by plants⁵⁻⁷. Previous studies have focused on spring leaf-out by daytime temperature and the onset of snowmelt time⁸⁻⁹, but the drivers controlling leaf senescence date (LSD) in autumn remain largely unknown¹⁰⁻¹². Using long-term ground phenological records (14536 time series since the 1900s) and satellite greenness observations dating back to the 1980s, we show that rising preseason maximum daytime (Tday) and minimum nighttime (Tnight) temperatures had contrasting effects on the timing of autumn LSD in the Northern Hemisphere (>20°N). If higher Tday leads to an earlier or later LSD, an increase in Tnight systematically drives LSD to occur oppositely. Contrasting impacts of daytime and nighttime warming on drought stress may be the underlying mechanism. A new LSD model considering these opposite effects improved autumn phenology modeling, and predicted an overall earlier autumn LSD by the end of this century compared with traditional projections. These results challenge the notion of prolonged growth under higher autumn temperatures, suggesting instead that leaf senescence in the Northern Hemisphere will begin earlier than currently expected, causing a positive climate feedback.

Climate change over the last several decades has modified the dates of plant flowering, leaf emergence, growth stages, and senescence, collectively termed phenology¹³ with substantial ecological and environmental consequences⁴. Both observations and model simulations have found that air temperature has a positive influence on the onset of plant growth in the Northern Hemisphere (NH), e.g., higher spring temperature triggers earlier
62 leaf-out and flowering dates and hence extends the growing season⁸,¹⁴⁻¹⁵. In contrast to
63 those extensive research efforts on spring phenology, autumn phenology, particularly leaf
64 senescence date (LSD), is more challenging to understand, and has not received
65 sufficient attention¹⁶,¹⁷, while also serving as an important indicator of changing foliar
66 physiological properties. Yet, autumn phenology may be as important as spring in
67 regulating the interannual variability of carbon balance⁷.

68
69 LSD has been occurring later in most regions over the last few decades¹⁸, but providing
70 an explanation for this change is difficult⁹. An increase in global temperature is assumed
71 to be a driver of LSD trends¹⁹, but studies indicated that the contribution of temperature to
72 LSD variability is low, especially compared to spring phenology²⁰,²¹. We argue that
73 ignoring the asymmetric effects²² of daytime maximum temperature (Tday) versus
74 nighttime minimum temperature (Tnight) and their differing impacts on LSD, contributes to
75 the reported overall low contribution of temperature to LSD variability. To test this, we
76 used measured and gridded preseason (defined as months from June to LSD) Tday and
77 Tnight in the NH, together with LSD data from three different datasets: (a) long-term
78 phenological observations at ground sites from 14536 time series since the 1900s (Fig.
79 S1), (b) the latest third generation of the Normalized Difference Vegetation Index (NDVI,
80 GIMMS3g.v1) for 1982-2015, and (c) NDVI and enhanced vegetation index (EVI) from the
81 Moderate Resolution Imaging Spectroradiometer (MODIS) products for 2001-2015.

82
83 Preseason forcing had a better predictive strength on LSD than either summer or autumn
climate forcing alone (Fig. S2). Because preseason $T_{day}$ and $T_{night}$ were highly correlated, we used a partial correlation to remove the effects of $T_{night}$ and of precipitation and radiation (similarly for $T_{night}$) to investigate the response of LSD to $T_{day}$. Correlations were classified into four types, (A) $T_{day}+T_{night}^{+}$, (B) $T_{day}-T_{night}^{-}$, (C) $T_{day}+T_{night}^{-}$ and (D) $T_{day}-T_{night}^{+}$, where $T^{+}$, $T^{-}$ represent positive and negative partial correlation coefficient of T with LSD.

Overall, all three datasets suggested that the onset of autumn LSD responded oppositely to $T_{day}$ and $T_{night}$. The proportions of ground sites of Types A and B were significantly lower than those of Type C and D (Fig. 1a). More significant R values for both $T_{day}$ and $T_{night}$ were found within Types C and D, with only two and one records out of 2231 time series having significant R within Type A and B, respectively. These results from ground sites are consistent with those for the two satellite greenness products (Fig. 1b, c). Types C and D together accounted for 83.7 and 80.0% for GIMMS3g and MODIS pixels, respectively. Only 0.8 and 1.5% of the pixels had the same sign of response of LSD to $T_{day}$ or $T_{night}$ (i.e. significant pixels for Types A+B) for GIMMS3g and MODIS, respectively. The GIMMS3g dataset contained different fractions of Types C and D (45.6% vs. 38.1%), but the compositions of Type C and D in GIMMS3g (i.e., contrasting effects of night and day temperatures) became more consistent with the MODIS results when the overlapping periods between the two sensors is considered (Figs. S3-4). More details on the fractions of the four correlation types for different vegetation types are provided in Figs. S5-6.

**Figure 1** Frequency of the partial correlation coefficient (R) between leaf senescence date
(LSD) and daytime (T<sub>day</sub>), nighttime (T<sub>night</sub>) temperatures. (a) Data for 14536 time series of ground sites, (b) the GIMMS3g dataset for 1982-2015, and (c) the MODIS product for 2001-2015. T<sub>+</sub>, T<sub>-</sub> represent positive and negative partial correlation coefficient of T with LSD. Significance was set at P < 0.05.

The satellite greenness products also allowed us to evaluate spatial patterns of LSD changes in response to variations in T<sub>day</sub> and T<sub>night</sub> (Fig. 2). For the GIMMS3g data, higher T<sub>day</sub> was associated with a delayed LSD for 10.7% of the pixels (mostly boreal regions) and with an earlier LSD for 7.5% of the pixels (central North America, borders of Eurasia and central China). T<sub>night</sub> had evident opposite influences on LSD than T<sub>day</sub>. The patterns of opposite effects from T<sub>day</sub> and T<sub>night</sub> on LSD were highly spatially consistent in all regions where T<sub>day</sub> and T<sub>night</sub> were significantly correlated with LSD. Similar results were obtained with MODIS observations (Fig. 2b, d). LSD for approximately 20% all pixels was significantly correlated with T<sub>day</sub>, of which 60.1 and 39.9% were negatively and positively correlated, respectively. The area where LSD was positively correlated with T<sub>night</sub> was larger (9.4%) than the area with negative correlations (6.5%).

Vegetation grouped into Köppen-Geiger zones showed contrasting patterns between the effects of T<sub>day</sub> and T<sub>night</sub> on LSD (Fig. 2e, f). Type D was more widely distributed, while Type C was more common for continental climates. Monsoon-influenced but not extremely cold regions and mild climates also had higher proportions of Type C. Grouping these correlation types by vegetation types lead to similar results (Fig. 2g, h). In theory,
we would expect to find Type C more in wet vegetation types and Type D in dry types. The
real world seems to show the same thing but still there could be many locations that do
not neatly fall into that continuum and suggests additional mechanisms may work,
probably adaption.

Figure 2 Spatial distributions of the partial correlation coefficient (R) between leaf senescence date (LSD) and daytime (T<sub>day</sub>), nighttime (T<sub>night</sub>) temperatures. R<sup>+</sup>, R<sup>-</sup> represent positive and negative partial correlation coefficient of T with LSD. (a) LSD vs. T<sub>day</sub> for GIMMS3g, (b) LSD vs. T<sub>day</sub> for MODIS, (c) LSD vs. T<sub>night</sub> for GIMMS3g, and (d) LSD vs. T<sub>night</sub> for MODIS. (e), (f), represent distributions of correlation types in Köppen-Geiger climatic classification using GIMMS3g and MODIS, respectively. (g) and (h) represent distributions of correlation types for vegetation types (see Methods) using GIMMS3g and MODIS, respectively. Significance was set at *P* < 0.05.

Our results suggest that ecological trade-offs, particularly those driven by regional
differences in water stress, may underlie the contrasting relationships between LSD and
T<sub>day</sub> and T<sub>night</sub>. Type C was mostly found in humid regions where water is a less limiting
factor for plant growth. In these cases, a higher T<sub>day</sub>, in the likely absence of severe water
stress, benefits photosynthesis while elevated T<sub>night</sub> increases nighttime leaf respiration.

Explanations for the prevalence of Type D relations in dry regions are more complicated.
The Standardized Precipitation Evapotranspiration Index (SPEI)<sup>23</sup>, an indicator of drought
stress, accounted for the contrasting effects of increases in T<sub>day</sub> and T<sub>night</sub> on LSD for Type D (Fig. 3). Partial correlation data indicate that increased T<sub>day</sub> is negatively correlated with SPEI (Fig. 3a), a stronger sensitivity to drought in dry regions that negatively affects plant growth and consequently leads to an earlier LSD. In contrast, we found that an increase in T<sub>night</sub> is associated with a higher SPEI, that is, wetter conditions, and arguably reduced water stress, which could extend the duration of photosynthesis and lead to delayed LSD (Fig. 3b). The latter mechanism is consistent with the generally positive partial correlation values between evapotranspiration (ET) and T<sub>night</sub>, that is, more soil moisture being available for ET in the late season, and sustaining delayed LSD (Fig. 3f), and with studies showing that water stress accelerates leaf drop in dry ecosystems more so than in humid ecosystems<sup>24</sup>. The responses of radiation to T<sub>day</sub> and T<sub>night</sub> may also be viewed as a further evidence for the linkage between leaf senescence and plant water status to support the contrasting patterns (Fig. S7), given that a higher T<sub>day</sub> was associated with stronger radiation and potentially a higher chance of water stress. These findings suggest that dry regions, in which Type D dominates, may be especially vulnerable to earlier onset of LSD if climate change reduces local precipitation and increases daytime evaporation with rising T<sub>day</sub>.
maximize carbon gain and minimize water losses<sup>25-27</sup>. The ecosystem-scale responses of phenology reflects the scaled responses of ecological dynamics of multiple individual species gaining or losing a competitive advantage in a changing climate, or presenting an induced advantage as a result of land-use change and planting<sup>17,26</sup>. Recent results suggest that the magnitude of phenological change to effects by shifts in plant species composition may be similar as that by climate change<sup>27</sup>, and the autumn phenology may thus change accordingly.

Figure 3 Partial correlation coefficient (R) between the Standardized Precipitation Evapotranspiration Index (SPEI), evapotranspiration (ET), and daytime (T<sub>day</sub>), nighttime (T<sub>night</sub>) temperatures. (a) SPEI vs. T<sub>day</sub> for GIMMS3g, (b) SPEI vs. T<sub>night</sub> for GIMMS3g, (c) SPEI vs. T<sub>day</sub> for MODIS, (d) SPEI vs. T<sub>night</sub> for MODIS, (e) ET vs. T<sub>day</sub> for MODIS, and (f) ET vs. T<sub>night</sub> for MODIS. Significance was set at *P* < 0.05.

We tested the implications of the observation analysis on future trends in autumn LSD by developing a weighted day-night-temperature growing-degree-day (DN<sub>GDD</sub>) algorithm incorporating these opposite changes in LSD to T<sub>day</sub> and T<sub>night</sub> (see Methods). The new model substantially improved LSD modeling (in terms of R (Figs. S8-10), RMSE (Figs. S11-13) and percentage of significant pixels (Figs. S14-15)) compared with the currently used threshold or GDD methods both for the overall dataset and for vegetation types.

Spatial patterns of improvements using MODIS and GIMMS3g were also investigated
194 (Figs. S16-17). The results from MODIS and the ground sites (Fig. S18) were more
195 consistent with the new model, and the accuracy of the threshold method was much lower,
196 so we used the coefficients from the MODIS data to predict LSD variability by the end of
197 this century using traditional GDD and DN<sub>GDD</sub> algorithms under two Representative
198 Concentration Pathways (RCP) scenarios (RCP 4.5 and RCP 8.5) (Fig. S19, and Fig. 4).

199
200 LSD from the DN<sub>GDD</sub> method was overall earlier than conventional predictions across
201 Köppen-Geiger climatic classification types. Globally, LSD was earlier for about 68% of
202 the terrestrial biosphere under RCP 4.5 and for about 70% under RCP 8.5. LSD was
203 mostly later for central North America, western Russia, and southwestern China. Most
204 vegetation types showed earlier LSD estimates under two RCP scenarios while the
205 temperate grasslands were expected to have later senescence dates.

206
207
208 Figure 4 Leaf senescence date (LSD) differences from the weighted
209 day-night-temperature growing-degree-day (DN<sub>GDD</sub>) and traditional GDD
210 (LSD<sub>DNGDD</sub>-LSD<sub>GDD</sub>) models under two RCP scenarios. (a), (b), (c) represent
211 LSD<sub>DNGDD</sub>-LSD<sub>GDD</sub> under RCP 4.5, and averages of differences for the Köppen-Geiger
212 climatic classification, and vegetation types, respectively. (d), (e), (f) represent
213 LSD<sub>DNGDD</sub>-LSD<sub>GDD</sub> under RCP 8.5, and averages of differences for the Köppen-Geiger
214 climatic classification, and vegetation types, respectively.

215
216 Climatic variability, particularly temperature, has driven phenological changes over the
217 last several decades but has been challenging to predict. Our ability to predict autumn
218 LSD is particularly limited. We are the first to report, using14536 ground time series and
219 more than 30 years of remotely sensed observations, the opposite responses of LSD to
220 daytime and nighttime warming, providing a new perspective to account for the previous
221 low estimation accuracy of autumn LSD when relying solely on mean temperature. A
222 model based on mean temperature cannot correctly predict LSD changes, because LSD
223 responds to Tday and Tnight in opposite directions. Our results also provide a perspective to
224 account for the carry-over effects between spring and autumn phenology, i.e. the start and
225 end of a growing season always move in the same direction²⁸. An earlier start of a season
226 is mainly triggered by higher spring temperatures, with increased growth depleting soil
227 water²⁹, which is then associated with autumn drought, inducing a reduction in growth and
228 consequently leading to an earlier end to the growing season³⁰.

229
230 Our improvement in modeling autumn phenology is a strong and convincing evidence for
231 the value of incorporating daytime and nighttime temperatures in terrestrial models, rather
232 than mean temperature alone. The application of this model projects an overall earlier
233 than currently expected start of autumn senescence in the NH by the end of this century,
234 particularly in dry regions. The earlier data of autumn senescence may be a potentially
235 unrecognized positive feedback to climate change and consequently a weakening in the
236 capability of terrestrial carbon uptake.

237
**Acknowledgments:** This work was funded by the Strategic Priority Research Program of the Chinese Academy of Sciences (XDA19040103), International Cooperation and Exchange Programs of National Science Foundation of China (Sino-German, 41761134082), National Natural Science Foundation of China (41522109) and the key Research Program of Frontier Sciences, CAS (QYZDB-SSW-DQC011). J.P. and P.C. were funded by European Research Council Synergy grant ERC-SyG-2013-610028 IMBALANCE-P. A.R.D. acknowledges support from the Ned P Smith Professorship of Climatology, University of Wisconsin-Madison.

**Contributions:** C.W., H.W. and Q.G. designed the research. C.W. wrote the first draft of the paper. J.P. and C.P. extensively revised the writing. H.W. performed the site model simulations. X.W. performed remote sensing model simulations. All the authors contributed to writing the paper.

**Correspondence and requests for materials should be addressed C.W.**

**Competing interests:** The authors declare no competing financial interests.

## References

1. Fu, Y. H., et al (2015). Declining global warming effects on the phenology of spring leaf unfolding. Nature, 526, 104-107.

2. Xia, J., et al (2015). Joint control of terrestrial gross primary productivity by plant
260 phenology and physiology. Proceedings of the National Academy of Sciences, 112,
261 2788-2793.

262 3. Buitenwerf, R., et al (2015). Three decades of multi-dimensional change in global leaf
263 phenology. Nature Climate Change, 5, 364.

264 4. Peñuelas, J., & Filella, I. (2009). Phenology feedbacks on climate change. Science, 324,
265 887-888.

266 5. Richardson, A. D., et al (2013). Climate change, phenology, and phenological control of
267 vegetation feedbacks to the climate system. Agricultural and Forest Meteorology, 169,
268 156-173.

269 6. Keenan, T.F., et al. (2014) Net carbon uptake has increased through warming-induced
270 changes in temperate forest phenology. Nature Climate Change, 4, 598-604.

271 7. Wu, C., et al (2013). Interannual variability of net ecosystem productivity in forests is
272 explained by carbon flux phenology in autumn. Global Ecology and Biogeography, 22,
273 994-1006.

274 8. Piao, S., et al (2015). Leaf onset in the northern hemisphere triggered by daytime
275 temperature. Nature communications, 6, doi:10.1038/ncomms7911.

276 9. Pulliainen, J., et al. (2017). Early snowmelt significantly enhances boreal springtime
277 carbon uptake. Proceedings of the National Academy of Sciences, 114,11081-11086.

278 10. Liu, Q., et al. (2016). Delayed autumn phenology in the Northern Hemisphere is
279 related to change in both climate and spring phenology. Global change biology, 22,
280 3702-3711.

281 11. Keenan, T.F., & Richardson, A. D. (2015). The timing of autumn senescence is
affected by the timing of spring phenology: implications for predictive models. Global Change Biology, 21, 2634-2641.

12. Gill, A.L., et al (2015). Changes in autumn senescence in northern hemisphere deciduous trees: a meta-analysis of autumn phenology studies. Annals of botany, 116, 875-888.

13. Myneni, R.B., et al (1997). Increased plant growth in the northern high latitudes from 1981 to 1991. Nature, 386, 698-702.

14. Suni, T., et al (2003). Air temperature triggers the recovery of evergreen boreal forest photosynthesis in spring. Global Change Biology, 9, 1410-1426.

15. Richardson, A.D., et al (2012). Terrestrial biosphere models need better representation of vegetation phenology: results from the North American Carbon Program Site Synthesis. Global Change Biology, 18, 566-584.

16. Gallinat, A.S., et al (2015). Autumn, the neglected season in climate change research. Trends in Ecology & Evolution, 30, 169-176.

17. Walther, G.R., et al (2002). Ecological responses to recent climate change. Nature, 416, 389-395.

18. Zhu, W., et al (2012). Extension of the growing season due to delayed autumn over mid and high latitudes in North America during 1982-2006. Global Ecology and Biogeography, 21, 260-271.

19. Garonna, I., et al (2014). Strong contribution of autumn phenology to changes in satellite-derived growing season length estimates across Europe (1982-2011). Global Change Biology, 20, 3457-3470.
20. Yang, Y., et al (2015). Changes in autumn vegetation dormancy onset date and the climate controls across temperate ecosystems in China from 1982 to 2010. Global change biology, 21, 652-665.

21. Jeong, S. J., et al (2011). Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982-2008. Global Change Biology, 17, 2385-2399.

22. Peng, S., et al (2013). Asymmetric effects of daytime and night-time warming on Northern Hemisphere vegetation. Nature, 501, 88-92.

23. Beguería, S., et al (2014). Standardized precipitation evapotranspiration index (SPEI) revisited: parameter fitting, evapotranspiration models, tools, datasets and drought monitoring. International Journal of Climatology, 34, 3001-3023.

24. Meng, T.T, et al (2009). Plant morphometric traits and climate gradients in northern China: a meta-analysis using quadrat and flora data. Annals of Botany, 104, 1217-1229.

25. Prasad, V.K., et al. (2007). Spatial patterns of vegetation phenology metrics and related climatic controls of eight contrasting forest types in India - analysis from remote sensing datasets. Theoretical and Applied Climatology, 89:95-107

26. Peñuelas J., et al. (2013). Evidence of current impact of climate change on life: a walk from genes to the biosphere. Global Change Biology, 19, 2303-2338.

27. Wolf, A.A., et al. (2017). Flowering phenology shifts in response to biodiversity loss. Proceedings of the National Academy of Sciences, 114, 3463-3468.

28. Fu, Y. S., et al (2014). Variation in leaf flushing date influences autumnal senescence and next year's flushing date in two temperate tree species. Proceedings of the National
326 Academy of Sciences, 111, 7355-7360.

327 29. Wolf, S., et al (2016). Warm spring reduced carbon cycle impact of the 2012 US summer drought. Proceedings of the National Academy of Sciences, 113, 5880-5885.

329 30. Peñuelas, J., et al (2017). Shifting from a fertilization-dominated to a warming-dominated period. Nature Ecology & Evolution, 1, 1438-1445.

331

332 Methods

333 1. Phenological observation data

334 We used observations of leaf senescence date (LSD) from three independent phenological datasets.

336 1) The Pan European Phenological Database (PEP725; http://www.pep725.eu/), an open-access database with long term plant phenological observations (since 1868) from 19608 sites and 78 species across 25 European countries.

339 2) The China Phenological Observation Network (CPON), with data since 1963 for >100 species at 42 sites across China.

341 3) LSD data for two tree species (*Acer palmatum* and *Ginkgo biloba*) at 54 meteorological stations in South Korea for 1989-2007³¹.

343 The definitions of LSD notably differ among the datasets. LSD for the PEP725, CPON, and Korean datasets is defined as the date when 50, 90, and 20% of the tree leaves, respectively, change color from green to red or yellow. We removed outliers using the methods³² to exclude potential biases and inadequate degrees of freedom and focused on time series with at least 15 years of records for 1900-2015. We thus analyzed 14536 LSD
348 time series for 24 species (Table S1).

349

350 2. LSD derived from satellite data

351 LSD in the Northern Hemisphere was determined using two satellite-derived
352 vegetation indices, the Normalized Difference Vegetation Index (NDVI) and the Enhanced
353 Vegetation Index (EVI)³³. Both NDVI and EVI are direct indicators of vegetation growth and
354 have been widely applied for investigating vegetation phenology³⁴. We used two datasets
355 to reduce the uncertainties caused by a single data source: GIMMS NDVI third-generation
356 (NDVI3g) data derived from the Advanced Very High Resolution Radiometer (AVHRR)
357 and NDVI and EVI derived from the Moderate Resolution Imaging Spectroradiometer
358 (MODIS). The GIMMS NDVI3g data have a spatial resolution of 1/12°, a half-month
359 interval, and a 34-year time span (1982-2015). The MODIS 16-day composite product
360 MOD13C1 (Collection 6) includes both NDVI and EVI with a 0.05° resolution for
361 2001-2015.

362 We eliminated the impacts of areas with sparse vegetation on the results by first
363 excluding pixels with annual NDVI <0.1 or annual EVI <0.08³⁵. A Savitzky-Golay filter was
364 then used to smooth the NDVI (EVI) time series³⁶. We then estimated LSD using two
365 methods.

366 The first method was a dynamic-threshold approach, which uses an annually defined
367 threshold for each pixel based on the NDVI ratio:

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

369 where NDVI is the daily NDVI and NDVI<sub>max</sub> and NDVI<sub>min</sub> are the annual daily maximum
and minimum of NDVIs. The NDVI<sub>ratio</sub> ranges from 0 to 1. LSD is determined when NDVI<sub>ratio</sub> decreases to 0.5 in autumn<sup>37,38</sup>.

The second method was based on a series of piecewise logistic functions. The NDVI time series were first divided into two sections by the maximum daily NDVI in each year, and a double logistic function was applied to fit each section<sup>39</sup>:

$$y(t) = a_1 + (a_2 - a_7t) \left[ \frac{1}{1 + e^{(a_3 - t) / a_4}} - \frac{1}{1 + e^{(a_5 - t) / a_6}} \right] \quad (2)$$

LSD was then defined as the time when the curvature changing rate reached its last local maximum value.

For GIMMS3g data, we calculated LSD using NDVI from both the dynamic-threshold approach and the piecewise logistic function method. Since MODIS sensor provides EVI, we further used EVI-based logistic function method to derive LSD. To sum up, for GIMMS3g data, average LSD from threshold and logistic function method were used, and for MODIS, an additional LSD from EVI-based logistic function method was used (not for MODIS NDVI data).

At high latitudes (or elevations), snow cover is important for regional climate and arrives early in autumn and potentially masking evergreen vegetation. However, we suggested that using a Savitzky-Golay filter could solve the noise from a "sudden" change in the time series of NDVI due to snow<sup>36</sup>. In particular, a study showed that snowfall had little influence on determining EOS in western Arctic Russia<sup>41</sup>. For high elevations, our previous analysis on Tibetan Plateau showed that for more than 98% of regions snow happened later than LSD<sup>42</sup>.
392 **3. Climatic data**

393 We used the CRU-TS 4.00 dataset with a spatial resolution of 0.5°×0.5° for
394 1901-2015. We extracted monthly data for Tmax, Tmin, Tmean, precipitation, and cloud cover
395 from this dataset for analyzing LSD from in-situ observations and the two remote sensing
396 data. We modeled past and future LSD by temperature by acquiring daily gridded data for
397 maximum and minimum temperature with a spatial resolution of 0.5° from NOAA Earth
398 System Research Laboratory Physical Sciences Division for 1982-2015. We used daily
399 Tmax and Tmin simulated by the CCSM 4 model under two climatic scenarios (RCP4.5 and
400 RCP8.5) for future climatic data (2081-2100). These data were from an open-access
401 database of the Coupled Model Intercomparison Project Phase 5.

402

403 **4. Analyses**

404 We used partial correlation analyses to determine the responses of LSD to Tday and
405 Tnight. The reason is that directly correlating LSD to Tday would give misleading results
406 since Tnight is a confounding variable that is numerically related to both LSD and Tday,
407 violating independence of variables in multiple correlation tests. Thus, using the partial
408 correlation between LSD and Tday will measure the degree of association with the effect of
409 a set of controlling random variables removed (e.g., Tnight, precipitation, radiation), given
410 that these factors have shown strong influences on LSD¹⁰, ²⁰. Since Echer and Silva (2014)
411 demonstrated that clouds are the main atmospheric factor modulating the surface
412 incidence of solar radiation⁴³, cloud cover data was used to model the effect of radiation
413 on LSD, as similarly conducted in previous analyses⁸. Correlations were classified into
four types, (A) T<sub>day</sub>+T<sub>night</sub>+, (B) T<sub>day</sub>-T<sub>night</sub>-, (C) T<sub>day</sub>+T<sub>night</sub>- and (D) T<sub>day</sub>-T<sub>night</sub>+, where T+, T-
represent positive and negative partial correlation coefficient of T with LSD. An R of at
least 0.514 for MODIS is required for the significance test (p = 0.05), but this value
decreases to 0.339 for the longer GIMMS3g data. These analyses were investigated for
both Köppen-Geiger climatic classifications and vegetation types (Table S2<sup>40</sup>). Crops
were excluded because their signal may result from changes in cropping or harvest cycles,
rather than from climate change. Furthermore, since at low latitudes plant phenology of
tropical and subtropical areas responds to other factors than temperature, regions with
latitudes < 20° N were also excluded.

Current phenology algorithms in most terrestrial-biosphere models are solely based
on temperatures in the preceding months<sup>44-45</sup>. We determined the length of the preseason
whose average T<sub>day</sub> had the most influence on LSD by calculating the partial correlation
coefficients between LSD and mean T<sub>day</sub> during 0, 1, 2, ... *n* months prior to LSD,
controlling for corresponding mean T<sub>night</sub>, total precipitation, and radiation. The maximum
range (*n*) of the preseason is generally from June to the multiyear mean date of LSD (see
Fig. S20 for example). The partial correlation coefficients with the highest absolute values
were then used in the following analysis. We obtained the relationship between LSD and
T<sub>night</sub> the same way but replacing T<sub>day</sub> with T<sub>night</sub>. This analytical procedure was applied for
observed LSD from ground sites and derived LSD from the MODIS and GIMMS NDVI<sup>3g</sup>
data.

434

5. Models for predicting LSD
Our results indicated that LSD responded oppositely to T<sub>day</sub> and T<sub>night</sub>, so we developed a weighted day-night-temperature growing-degree-day (DN<sub>GDD</sub>) algorithm from observations to model LSD, and compared the algorithm with currently used threshold and GDD models based on T<sub>mean</sub><sup>46</sup>.

The threshold model was the simplest method. We calculated average T<sub>mean</sub> for five days before LSD in each year and used the multiyear mean value as the threshold to model LSD. If T<sub>mean</sub> was lower than the threshold for five consecutive days from 1st July, the last date was considered the LSD.

GDD was calculated as:

$$GDD(d) = \max(T_b - T_{mean}(d), 0) \qquad (3)$$

$$GDD_{threshold} = \sum_{d=d_0}^{LCD} GDD(d) \qquad (4)$$

where T<sub>b</sub> is the base temperature set to 15, 20, 25, and 30 °C, T<sub>mean</sub>(d) is the mean daily temperature, and d<sub>0</sub> is the date on which the calculation begins (1st July in this study). LSD is the observed or derived date of leaf coloring in each year. The date when GDD(d) exceeded the multiyear average GDD threshold was defined as LSD.

Our DN<sub>GDD</sub> model improved upon the original GDD model and was calculated by:

$$GDD(d) = k \cdot \max(T_b - T_{day}(d), 0) + (1-k) \cdot \max(T_b - T_{night}(d), 0) \qquad (5)$$

where T<sub>day</sub>(d) is the daily maximum temperature, T<sub>night</sub>(d) is the daily minimum temperature, and k is the weighting factor. When 0<k<1, the effects of T<sub>day</sub> and T<sub>night</sub> on LSD are consistently positive; When k>1 or k<0, the effects of T<sub>day</sub> and T<sub>night</sub> on LSD are opposite. In order to determine the value of k, we first calculated the ratio of R<sub>day</sub> and R<sub>night</sub> for each station or pixel, and found 99.9% of the ratio values were between -10 and 10 for
both ground and satellite data (Figure S21). In other words, the level of $T_{day}(T_{night})$ effect could be 1 to 10 times than the level of $T_{night}(T_{day})$ effect (note that $T_{day}$ represents $T_{day}$ with the effects of $T_{night}$ removed). Therefore, the values of *k* ranged from -1 to 2 (see Table S3). In addition, when *k* tends to infinity, the effects of $T_{day}$ and $T_{night}$ on LSD are opposite with same level.

We evaluated the accuracy and obtained the most appropriate parameters of the models by calculating the correlation coefficient (R) and the root mean square error (RMSE) between modeled and observed LSD. $T_b$ and *k* with the lowest RMSE were considered the most appropriate values for each site or pixel.

## Data Availability

The data that support the findings of this study are available from the corresponding author upon request.

## References

31. Park, C.K., et al., (2017). Spatial and temporal changes in leaf coloring date of Acer palmatum and Ginkgo biloba in response to temperature increases in South Korea. PLoS ONE, 12.

32. Schaber, J., & Badeck, F.W. (2002). Evaluation of methods for the combination of phenological time series and outlier detection. Tree Physiology, 22, 973-982.

33. Huete, A., et al., (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote sensing of environment, 83, 195-213.
34. Zhang, X., et al. (2003). Monitoring vegetation phenology using MODIS. Remote sensing of environment, 84, 471-475.

35. Shen, M., et al., (2014). Increasing altitudinal gradient of spring vegetation phenology during the last decade on the Qinghai-Tibetan Plateau. Agricultural and Forest Meteorology, 189, 71-80.

36. Chen, J., et al., (2004). A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky-Golay filter. Remote sensing of environment, 91, 332-344.

37. Piao, S., et al., (2011). Changes in satellite-derived vegetation growth trend in temperate and boreal Eurasia from 1982 to 2006. Global Change Biology, 17, 3228-3239.

38. White, M.A., et al. (1997). A continental phenology model for monitoring vegetation responses to interannual climatic variability. Global Biogeochemical Cycles, 11, 217-234.

39. Elmore, A.J., et al., (2012). Landscape controls on the timing of spring, autumn, and growing season length in mid-Atlantic forests. Global Change Biology, 18, 656-674.

40. Olson, D. M., et al. (2001). Terrestrial ecoregions of the world: a new map of life on earth. Bioscience, 51, 933-938.

41. Zeng, H., & Jia, G. (2013). Impacts of snow cover on vegetation phenology in the arctic from satellite data. Advances in Atmospheric Sciences, 30, 1421-1432.

42. Wang, X. et al. 2018. Snow cover phenology affects alpine vegetation growth
502 dynamics on the Tibetan Plateau: Satellite observed evidence, impacts of biomes,
503 and climate drivers. Agricultural and Forest Meteorology, 256, 61-74.

504 43. Silva, A. A., & Echer, M. P. D. S. (2013). Ground-based measurements of local cloud
505 cover. Meteorology & Atmospheric Physics, 120, 201-212.

506 44. Delpierre, N., et al. (2009). Modelling interannual and spatial variability of leaf
507 senescence for three deciduous tree species in france. Agricultural & Forest
508 Meteorology, 149, 938-948.

509 45. Richardson, A.D., et al., (2012). Terrestrial biosphere models need better
510 representation of vegetation phenology: results from the North American Carbon
511 Program Site Synthesis. Global Change Biology, 18, 566-584.

512 46. Krinner, G., et al., (2005). A dynamic global vegetation model for studies of the
513 coupled atmosphere-biosphere system. Global Biogeochemical Cycles, 19,
514 doi:10.1029/2003GB002199.
<figure><img src="image_2.png" /></figure>

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

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

Partial correlation coefficient

-1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 non-significant
<figure><img src="image_6.png" /></figure>

Partial correlation coefficient

-1.0 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1.0 non-significant

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