Autumn canopy senescence has slowed down with global warming since the 1980s in the Northern Hemisphere

Climate change strongly impact vegetation phenology, with considerable potential to alter land-atmosphere carbon dioxide exchange and terrestrial carbon cycle. In contrast to well-studied spring leaf-out, the timing and magnitude of autumn senescence remains poorly understood. Here, we use monthly decreases in Normalized Difference Vegetation Index satellite retrievals and their trends to surrogate the speed of autumn senescence during 1982 – 2018 in the Northern Hemisphere (>30°N). We ﬁ nd that climate warming accelerated senescence in July, but this in ﬂ uence usually reversed in later summer and early autumn. Interestingly, summer greening causes canopy senescence to appear later compared to an advancing trend after eliminating the greening effect. This ﬁ nding suggests that summer canopy greening may counteract the intrinsic changes in autumnal leaf senescence. Our analysis of autumn vegetation behavior provides reliable guidance for developing and para-meterizing land surface models that contain an interactive dynamic vegetation module for placement in coupled Earth System Models.

A utumnal senescence is a critical component of the seasonal evolution of terrestrial ecosystems, which regulates the length of the growing season and terrestrial carbon uptake 1,2 .Yet, such event and its response to changing climate is not well represented in most Earth System Models (ESMs) 3 .Impact of climate change on the trends of autumn senescence as well as its role in terrestrial carbon, water, nutrient cycling, fitness, and distribution of tree species have been reported over the past decades 1,2,[4][5][6][7][8][9] .Many of them focused on the timing of specific phenophase, e.g., leaf coloring and abscission 1 .However, the process of autumn senescence consists of a series of complex metabolic events such as macromolecule degradation, decreased photosynthesis, and most importantly, nutrient reabsorption [10][11][12] .These events are co-dominated by plant internal conditions (e.g., metabolic adjustments, genetic expressions), near-surface meteorological factors (e.g., temperature, solar radiation and wind), water availability and their interactions [12][13][14] .Despite the complex climatic responses of autumn phenology, it remains unclear, however, which stage of the senescence process contributes more to changes in the timing of end of the growing season (EOS).Moreover, previous studies pay too much attention to individual phenophases (e.g., EOS), which may conceal other components of the senescence process, such as the variations in the environmental regulations during the multiple stages of the autumn senescence process 15 .These gaps in knowledge would leave the reliable modelling of autumn senescence in doubt 16 , and finally hamper the ability of ESMs to simulate vegetation dynamics under future climatic scenarios 17,18 .
Here, we provide a timely assessment of autumn senescence speed through the monthly decreases in NDVI (ΔNDVI, see schematic representation in Fig. 1), which is also critical for the extraction of EOS.We first investigate the regulations of climatic factors (temperature, soil moisture, insolation i.e., downward shortwave radiation, and windspeed) on the speed of autumn senescence (i.e., the interannual trends in ΔNDVI).Given that such behavior may be confounded by rising CO 2 -induced summer canopy greening (i.e., uptrends in annual maximum NDVI, NDVI max ), we additionally calculated the relative decrease of NDVI (δNDVI = ΔNDVI/NDVI max ) as a surrogate for ΔNDVI (see definitions in Supplementary Table 1).Then, we quantified the effects of increasing NDVI max (i.e., canopy greening) on ΔNDVI, and consequently, the extracted EOS.Together with this analysis, we discussed the implications of our study on benchmarking phenological modules for future generations of Earth System Models (ESMs).

Results and Discussion
Changes in the speed of autumn senescence.Across the Northern Hemisphere (>30°N, NH, thereafter), divergent trends in ΔNDVI can be discerned among the senescence months, with larger increasing trend found in July (1.8 × 10 −4 yr −1 ) and September (2.1 × 10 −4 yr −1 ), whereas decreasing trends were found for August (−3.7 × 10 −5 yr −1 ) and October (−1.7 × 10 −4 yr −1 ; Fig. 2a).There is also an increasing trend in ΔNDVI for the entire senescence period (July to October, JASO) (Fig. 2b), although noting this may be expected due to rising NDVI max under planetary greening (Fig. 2c).Dramatic regional variations, however, existed in ΔNDVI trends across the NH (Fig. 2d-i).Specifically, we found ΔNDVI increased in northwestern North America in July, which dominated the acceleration of autumnal senescence during the senescence period (Fig. 2d, h).In contrast, negative ΔNDVI trends were mainly distributed in northeastern North America and eastern Eurasia during August, which contributed to a slowing of the speed of autumn senescence (Fig. 2e, h).In September, contrasting changes in ΔNDVI were observed between the high and middle latitudes of Eurasia (Fig. 2f), resulting in opposite trends in the speed of autumn senescence (e.g., acceleration in middle Eurasia but deceleration in higher latitudes) (Fig. 2h).We also found that the speed of autumn senescence slowed down in the west and east of North America in October, when vegetation was already in dormant phase across most high latitudes (Fig. 2g).
Overall, ΔNDVI over the entire senescence period increased in 53.1% of the NH (significant in 17.9%, Fig. 2h), implying an overall acceleration of autumn senescence would occur in these areas.Although acceleration of senescence was found in more than half of the NH, the remaining NDVI (NDVI left in the month of EOS, see definition in Fig. 1 and Supplementary Table 1) still exhibited positive trends in 69.4% of the NH (Fig. 2i).These findings are anticipated due mainly to the apparent summer greening (i.e., a prevalent increase in NDVI max , Supplementary Fig. 1).Indeed, when averaged across NH, NDVI max increased significantly at a rate of 8.4 × 10 −4 yr −1 (p < 0.001) during the senescence period (Fig. 2c).Meanwhile, ΔNDVI over JASO general increased, but its magnitude (2.1 × 10 −4 yr −1 , p > 0.05) is weaker than that of NDVI max .As a result, the remaining NDVI had an overall positive trend (6.3 × 10 −4 yr −1 , p < 0.001).Our analysis, therefore, reveals that changes in the speed of autumn senescence (ΔNDVI) and the peak greenness (NDVI max ) can affect vegetation greenness at the later stage of the growing season.
Climatic regulation on the speed of autumn senescence.To investigate the climatic regulations on the speed of autumn senescence, we conducted partial correlation analyses to isolate the influence of insolation, wind speed, soil moisture, and temperature (Supplementary Fig. 2).The factor with the strongest partial correlation with ΔNDVI was defined as the dominant driver of autumn senescence (see Methods).Again, the dominant factors that accelerate or slow down the speed of autumn senescence were heterogeneous and discernible in their patterns over the NH (Supplementary Fig. 3).
Our noted strong spatial variations in ΔNDVI trends may be explained by the divergent biological responses of plants to In this study, autumn senescence starts at Peak Of Season (POS) when NDVI reaches its annual maximum (NDVI max ), and ends at the End Of Season (EOS).The monthly decrease in NDVI (ΔNDVI) is estimated as the NDVI difference between the current and the following months (in positive values), and the remaining NDVI is the NDVI left in the EOS month (in positive values).The sum of monthly ΔNDVI and the remaining NDVI is therefore equal to NDVI max .Detailed definitions of these items can be found in Supplementary Table 1.
climate change across different climatic regions.In fact, climatic regulations on ΔNDVI trends (higher values indicate a greater senescence speed, Fig. 1) were mainly associated with local hydrothermal conditions (Fig. 3).At the beginning of the autumn senescence period (that is, July), ΔNDVI correlated positively with temperature, but negatively with soil moisture in colder (mean annual temperature, MAT < 5 °C) and relatively drier (mean annual precipitation, MAP < 500 mm) regions (Fig. 3a, f).This finding indicated a possible acceleration of senescence caused by warming, because climatic warming may in turn lead to severe soil water deficits.Meanwhile, insolation was identified as a major promotor of the speed of autumn senescence in 33.2% of NH (Fig. 3k).Interestingly, the overall temperature effects on ΔNDVI shifted from positive (July) to negative (August, September, and October) as autumn senescence progresses.The negative temperature dominance on ΔNDVI trends over the NH increased from ~31% (August and September, Fig. 3l, m) to 46.9% (October, Fig. 3n).These results imply that a warming climate could accelerate autumn senescence in earlier stages of the senescence period, probably due to the detrimental summer extreme temperatures.On the other hand, it might significantly slow down the speed of senescence in its later stages by prolonging the exposure of vegetation to favorable summer climate.
Unlike such reversal of temperature regulation on the speed of autumn senescence, the effects of insolation, soil moisture and windspeed remained almost consistent along the senescence period.To be specific, insolation was identified as the positive dominant factor of ΔNDVI, particularly over the warmer regions (MAT > −5 °C) and the lower latitudes (Fig. 3).Because the increasing insolation might raise the risk of leaf photo-oxidative damage, which could promote the progress of leaf senescence 19 .As for soil moisture and wind speed, their influence is only apparent in a few sporadic areas.
While for the remaining NDVI (higher values indicate slower senescence speed, Fig. 1), it was generally positively related to temperature and insolation, but negatively related to wind speed (Fig. 3e, j, o).The positive temperature effect was prevalently observed over the NH (Supplementary Figs.2x, 3f), which confirmed again that warming climate contribute to slow down the speed of senescence in late autumn.Followed by insolation, its positive and dominant effects on the remaining NDVI were mainly distributed in North Siberia and North America (with MAT < 0 °C) (Supplementary Figs.2e, 3e), indicating the negative light regulation on leaf senescence was stronger than that of warming for these regions 20 .We also found that the wind speed exerted negative controls on the remaining NDVI in 34.3% of NH (Fig. 3e, o) through its role in promoting leaf abscission and intensifying drought stress by increasing evapotranspiration 7 .The influence of soil moisture expressed a contrasting regional difference (Supplementary Fig. 2q).More specifically, in terms of inland Eurasia and North America, the remaining NDVI was positively regulated by soil moisture (Supplementary Fig. 2q), indicating that increases in soil moisture can slow down the speed of autumn senescence in these water-limited ecosystems.However, in high-latitude Siberia (>50°N), soil moisture was negatively correlated with the remaining NDVI (Supplementary Fig. 2q), mainly due to the anoxia conditions under high soil moisture 21 and limited nutrient availability in these regions of permafrost 6 .
The effect of summer greening on the extraction of autumn phenology.Parallel to climate change, the observed summer greening (that is, an increase of 8.4 × 10 −4 yr −1 in NDVI max , p < 0.001, Fig. 2c) could modulate the extracted autumn phenology and its temporal trend.We hypothesize that, accompanied with summer greening, the increasing mature leaves in the peak season would contribute to an increasing trend in ΔNDVI, thus creating an illusion of faster senescence (Fig. 4a).Such effect could be greatly reduced by replacing monthly senescence (ΔNDVI) with the relative decrease in NDVI (δNDVI = ΔNDVI/ NDVI max ).Summer greening, in a similar way, would contribute to delayed EOS calculated from fixed NDVI thresholds (EOS a ) than that from fixed percentage of NDVI max (EOS p ), implying the influence of changing NDVI max in EOS detection (Fig. 4b, also see definitions in Supplementary Table .

1).
Notably, trends in δNDVI generally presented similar spatial patterns as that of ΔNDVI in each of the senescence months (Supplementary Fig. 4a-e).After separating the contribution of trends in NDVI max and δNDVI to trends in ΔNDVI (Supplementary Fig. 5a-e, see Methods), we found that the increase in NDVI max accounted for ~40% of ΔNDVI trends in NDVI for parts of Siberia and central Eurasia in August and September (Supplementary Fig. 5b, c).Despite this, the overall contribution of NDVI max to ΔNDVI trends was marginal for other regions and periods (e.g., only 1.7, 5.9, 5.0 and 2.6% for Jul, Aug, Sep and Oct, respectively).Such influence was confirmed by the relatively small portion of significant correlations between trends in NDVI max and ΔNDVI for these months (Supplementary Fig. 6a-d).Besides, similar spatial patterns of dominant climatic factors were also observed in δNDVI, further indicating their tight connection (Supplementary Fig. 3).
We also found that the influence of increasing NDVI max on trends of the remaining NDVI was apparent (36.7% on average, Supplementary Fig. 5f), which was in line with the widespread significant correlations between trends in NDVI max and the remaining NDVI over the NH (Supplementary Fig. 6f).This leads to higher fractions of positive trend areas for the remaining NDVI (69.4% with 30.1% significant, Fig. 2i) than its percentage-based counterpart (55.3% with 14.7% significant, Supplementary Fig. 4f).Obviously, the impact of conspicuous summer greening on the remaining NDVI may have cascaded effects on the extraction of EOS.
We therefore further compared the temporal trends of EOS a and EOS p (definitions in Fig. 4b, also see definitions in Supplementary Table. 1) during the period 1982-2018.Across the NH, we found areas occupied with delayed EOS a (60.2%, significant in 22.2%, Fig. 5a) exceeds that of EOS p (51.1%, significant in 15.1%, Fig. 5b).As a result, an overall delayed trend of EOS a (0.6 days dec −1 ) but advanced EOS p (−0.4 days dec −1 ) trend was presented in NH (Fig. 5a-c).This divergence summed up to 1.0 days dec −1 NDVI max -induced delay in EOS a .Given that such an advance in autumn senescence was supported by a metaanalysis based on field observations 22 , we reckon that EOS p may be a better indicator of autumnal leaf phenology, which captures the physiological changes of plants.EOS a , carrying the additional effect of NDVI max , my instead, serve as a better indicator for canopy greenness browning.Our results thus infer that the increase in the number of mature leaves in the canopy under summer greening may counteract the intrinsic changes in autumnal leaf senescence.
Among the studied plant functional types (PFTs), we also observed clear mismatch between trends in EOS a and EOS p (Fig. 5d).In line with the NH average (a combination of all PFTs), we found the opposite contributions of NDVI max (delay) and EOS p (advance) to EOS a trends in grasslands, savannas and shrublands in boreal and arctic regions.However, the delaying effect of summer greening could hardly compensate for the large advancement in EOS p , as these ecosystems generally locate in resource-limited areas (Supplementary Fig. 7).Growth enhancement (greening) during the earlier season may prematurely consume water and nitrogen, which would in turn speed up the progress of vegetation senescence and cause earlier end of growing season [23][24][25] .Indeed, the largest advancement in EOS p was observed in grasslands of central Eurasia, which were characterized with severe water limitations and significant vegetation greening (Supplementary Fig. 1).
The observed delay in EOS p for the woody species could be ascribed to the higher resource availability in their relatively wetter and fertile habitants, and their ability to accumulate large nutrient and carbohydrate reserves in autumn 26,27 .The increase in NDVI max contribute to later EOS a than EOS p , albeit varied across these PFTs (Fig. 5d).We quantified the contribution of increasing NDVI max to the delayed EOS a through the residue proportion of trends after subtracting EOS p from EOS a (i.e., the proportion of orange bars within the transparent purple bars, last section of methods).As a result, the increase in NDVI max contributed to more than half of the delay in EOS a for woody savannas (69.2%) and temperate shrublands (88.5%).However, such influence is minor in forests, contributing to 20.3% and 9.5% in mixed forests and deciduous broadleaf forests, respectively.On the contrary, the decrease in NDVI max in deciduous needleleaf forests advanced EOS a by 29.3%, comparing to EOS p .
Remarkably, NDVI max -induced divergence between EOS a and EOS p trends exerted strong controls on their responses to climate change (Supplementary Figs. 2, 8).Despite the observed clear positive temperature dominance for EOS a and EOS p (Supplementary Fig. 3f, r), the difference in the sensitivities of EOS a and EOS p to temperature was examined (e.g., EOS a (γ a ) and EOS p (γ p ), see Methods).The partial least squares regressions analysis suggested that both EOS a and EOS p had positive temperature sensitivities in most areas, while negative sensitivities were concentrated in xeric grasslands of central Eurasia and USA (Supplementary Fig. 9a, b).Note that EOS a presented larger temperature sensitivities than EOS p in 69.5% of northern areas, e.g., γ a (2.1 days °C−1 ) was more than 50% higher than γ p (1.3 days °C−1 ) when averaged over the NH (Supplementary Fig. 9c), due mainly to the additional NDVI max greening trend carried by EOS a .Indeed, the prominent differences between γ a and γ p were discovered in grasslands and temperate shrublands (Supplementary Fig. 9d), where EOS a trends were predominantly contributed by changes in NDVI max (Fig. 5d).
Implications for model benchmarking.The mismatch between phenology metrics based on absolute (ΔNDVI, EOS a ) and relative (δNDVI, EOS p ) changes in NDVI indicates the considerable effect of summer greening on the extraction of autumn phenology.In other word, the increase in the number of mature leaves in canopy under summer greening might confound the actual changes in the timing of leaf senescence.Such an effect is most clearly shown by the mismatch between EOS a and EOS p .EOS a , derived from fixed thresholds of NDVI, is a direct reflection of drop in canopy greenness, with a delayed trend in EOS a indicating the prolonged canopy browning process.EOS p , derived from fixed percentages of NDVI max , instead, better captures senescence of leaves under the overall changes in canopy greenness.As shown in our results, the widely observed global greening (Supplementary Fig. 1) in this study and previous studies 28,29 , has, indeed, counteracted the intrinsic advance of leaf senescence by an apparent delay in canopy browning.Despite that these two metrics may represent phenology at different ecological levels (e.g., leaf vs. canopy), leaf phenology is more closely related to the time when plants enter winter dormancy and is a key regulator of vegetation dynamics and terrestrial carbon sequestration 2 .The absence of eliminating the summer greening effect in the extraction of autumn phenology (i.e., threshold-based methods: EOS a ) may introduce bias to our understanding of the leaf senescence process and consequently, the vegetation-climate interactions and the terrestrial carbon cycle.Since the timing of leaf senescence is crucial for reconstructing historical plant carbon uptake and predicting its future changes in Earth system models, it is essential to update phenology modules embedded in state-of-the-art land surface schemes using accurate representations of these processes.This will enable a better simulation of the terrestrial carbon cycle 30,31 .
However, leaf senescence so far is one of the most difficult processes to be parameterized in terrestrial ecosystem models, which is restricted by our incomplete understandings of its complex mechanisms.To our knowledge, DGVMs coupled with land surface models generally determine the onset of plant dormancy by tracking the thresholds of modelled or prescribed leaf area index 30 .For example, in the ORCHIDEE model, it is defined as the date when the modelled leaf area index (LAI) drops below 0.2 32 .Evidence from our study suggests that fixed thresholds, however, may not accurately reproduce the timings of autumn leaf phenology when the summer greening effect is not thoroughly considered.Other models, for example, the Ecosystem Demography model version 2 (ED2), used the prescribed satellite-derived EOS estimates to parameterize the date of leaf drop during the year 33 .In this case, an accurate estimation of EOS, especially taking the effect of summer greening into consideration, is of great importance for accurate model parameterizations.Moreover, a large body of models simulates the timing of leaf shed by a function of air temperature and (or) day length [34][35][36] .Although these models directly simulate leaf phenology, caution is still needed for benchmarking models against satellite-derived results, since phenology based on absolute changes in vegetation greenness may be not sufficient to characterize the trajectories of leaf senescence.
We also found that EOS a , carrying the additional greening trend at POS, generally exhibits higher temperature sensitivities than EOS p .This discrepancy was in support of findings in Zhang et al. (2020), who reported that the response of greenness-based phenology to climate change was greater than that of photosynthesis-based phenology.Therefore, if the apparent sensitivity of satellite-based leaf phenology (with the absence of a summer greening effect in the extraction of phenology) is used to parameterize leaf behavior in DGVMs, the future carbon sink potential of terrestrial ecosystems might be overestimated.
Our study provides a novel angle on the speed of autumn senescence and its response to climate change over the NH.We find that warming climate accelerates senescence in July, but this influence usually reversed in later summer and early autumn.Insolation generally accelerates senescence across the entire senescence period.Our analysis also verifies a summer greening (increase in NDVI max )-induced delay in EOS, as canopy greening may counteract the intrinsic changes in leaf phenology, and lead to higher temperature sensitivities for EOS carrying the additionally greening trend.Global warming and increasing atmospheric CO 2 have been reported to continually contribute to earth greening, especially in NH 28,29 .Our analysis thus allows us to reconsider the linkage between changes in apparent canopy greenness and changes in intrinsic leaf physiology, which could provide an important benchmark for those developing DGVMs for use in ESMs.

Methods
Datasets.In this study, NDVI during the 1982-2018 period was extracted from the third-generation product of NASA's Global Inventory Modeling and Mapping Studies (GIMMS) group 37 .Various issues, such as orbital drift, calibration, viewing geometry and volcanic aerosols, have been corrected in this product 38 .This dataset consists of fortnightly NDVI observations at a spatial resolution of 0.083°(~8 km).The monthly NDVI data was estimated through the maximum composition of each month's two records.In this analysis, we screened out pixels dominated by cropland and evergreen forests (inferred from the MODIS Landcover Product, MCD12C1, IGBP classification) 39 , because their seasonal cycles are unclear.Consistent with the NDVI time span, the monthly mean 2-m surface temperature was derived from CRU.TS.4.05 with 0.5°spatial resolution 40 .Monthly insolation was extracted and summed from CRU-JRA v2.2 with 0.5°spatial and 6-h temporal resolution 41 .The monthly surface soil moisture was acquired from the C3S dataset provided by the European Center for Medium-Range Weather Forecasts (ECMWF) with a 0.25°spatial resolution, which was then interpolated to 0.5°× 0.5°.The monthly wind speed was obtained from the ERA5 reanalysis data at a spatial resolution of 0.25°4 2 and was interpolated to 0.5°× 0.5°.
Phenology extraction methods.Our study was conducted across north of 30°N (excluding the subtropical regions), where plants have clear seasonal variations and are capable to extract phenological metrics.The Autumn senescence period is defined as the time window between POS (Supplementary Fig. 10) and EOS (Supplementary Fig. 11).The decrease in NDVI (ΔNDVI) during the entire senescence period was defined as the difference in NDVI between POS month and EOS month.Similarly, the monthly ΔNDVI was the NDVI difference between the current and the following month.The remaining NDVI i.e., NDVI at EOS month, refers to the remaining greenness in the final stage of autumn senescence.For the minor regions with later POS (e.g., later than July, only 6.2% of the study area), the calculation of ΔNDVI also starts from POS month, NDVI data prior to the POS month was excluded from the analysis.
Besides monthly senescence, we further obtained threshold-based EOS (EOS a ).To avoid the influence of potential snow cover on EOS extraction, we first used the spline function to interpolate the monthly temperature to a daily basis, and then obtained the boundary of the thermal growing season by a sequence of five days below 0 °C.Before EOS extraction, NDVI pixels outside this thermal growing season were replaced by that of the temporally nearest snow-free date.We than determined EOS a according to the following steps (Polyfit-Maximum method 43 ; Supplementary Fig. 12): (1) calculate the averaged NDVI curves during the past 37year (1982-2018), (2) estimate the rates of changes in the average NDVI curves as: NDVI ratio (t)=[NDVI(t + 1)-NDVI(t)]/NDVI(t), (3) determine the time T when the minimum NDVI ratio occurred, (4) use the NDVI at time T (i.e., NDVI(T)) as the threshold for EOS a extraction, (5) apply the threshold to determine EOS a for each year.The Polyfit-maximum method have been successfully used in many of previous phenological studies [44][45][46] , and its outcomes are comparable with other commonly used methods 6 .The priorly determined NDVI threshold also gives us the chance to correspondently obtain a percentage-based threshold, which is then used for further analysis (see detection of EOS p in next paragraph).
The extraction of ΔNDVI and EOS a is probably confounded by global greening, especially in the context of the prominent increase in NDVI max (i.e., the effect of summer canopy greening, Fig. 4 and Supplementary Fig. 1).To reduce such an impact, we employed another method based on the percentage of NDVI max .During each month and the entire senescence period, we calculated the percentage of decrease in NDVI compared to the peak season NDVI as δNDVI = ΔNDVI/ NDVI max .Similarly, the percentage of remaining NDVI relative to NDVI max was calculated as the division between the remaining NDVI and NDVI max .The percentage-based EOS (EOS p ) was extracted from almost the same steps as EOS a , except that the percentage of NDVI(T) relative to NDVI max (NDVI(T)/NDVI max ) were used as thresholds to detect yearly EOS P (step 4, Supplementary Fig. 12).The differences in ΔNDVI and δNDVI, and EOS a and EOS p gives us chance to investigate the effect of NDVI max on autumn phenology detection (see last section of methods).
Climatic regulations on the process of autumn senescence.Benefiting from previous studies, we examined the effects of four climatic factors i.e., insolation, wind speed, soil moisture, and temperature on the leaf senescence process 6,7,14,47 .Note that we did not take precipitation into consideration, due mainly to the following aspects: (1) Soil moisture is extensively correlated with precipitation in the study area (Supplementary Fig. 13) and (2) soil moisture is more closely and directly related to plant water uptake than precipitation 5,48 .The lagged effect of climate on vegetation dynamics was also considered for each of the phenological matrices (ΔNDVI, δNDVI, EOS a and EOS p ) following previous studies 6,8,15 .Firstly, simple correlations were used to determine the preseason length of each climatic factors, which was defined as the period when the mean value of factors having the largest absolute simple correlation with phenological metrics (Supplementary Figs.14, 15).The preseason length ranged from the month of phenology metrices to its preceding 3 months, with a monthly time step.The preseason was defined for each of the four phenological metrics and climatic variables.We then performed partial correlation analysis between each of the phenological metrices and climatic factors during the preseason (Supplementary Figs. 2, 8).It enables us to explore the relationship between plant phenology and single climatic factor while eliminating the effects of the remaining factors.The climatic factor with the largest positive (lowest negative) partial correlation coefficient was defined as the positive (negative) dominant factor (Supplementary Fig. 3).
Temperature sensitivity of EOS.Our partial correlation analysis suggested distinct temperature dominances on EOS a and EOS p , while similar patterns were not clear for ΔNDVI and δNDVI.We then performed a partial least squares regression between EOS a (also EOS p ) and the four climatic factors, and assigned the coefficient of temperature as the temperature sensitivity of EOS a (γ a ) and for EOS p (γ b ).The coefficient was commonly used as a proxy of changes in phenological dates per unit increase in temperature (i.e., temperature sensitivity) [49][50][51] .
Separation of the effects of NDVI max in ΔNDVI and EOS a .In this study, ΔNDVI for each month of senescence is also expressed as: This formula inherently indicates that ΔNDVI is jointly determined by δNDVI and NDVI max .Trends of ΔNDVI during 1982-2018, therefore, can be separated into the trends of δNDVI and NDVI max based on the multiplication derivation rule ((uv)'=u'v + v'u, u and v are variables).The following approximation was used: where tr[ΔNDVI], tr[δNDVI] and tr[NDVI max ] are trends of ΔNDVI, δNDVI and NDVI max during the period 1982-2018, respectively.mean[NDVI max ] and mean[δNDVI] are the multi-year averages of NDVI max and δNDVI.tr [δNDVI] x mean[NDVI max ] and tr[NDVI max ] x mean[δNDVI] thus denote the contribution of δNDVI trends and NDVI max trends to ΔNDVI trends, respectively.To test the reliability of such approximation, we additionally checked whether the sum of the terms on the right-hand side of Eq. 2 (i.e., the partitioning approach) is equivalent to the left-hand side.The feasibility of such partitioning is ensured by a perfect 1:1 distribution of values on both side of Eq. 2 (Supplementary Fig. 16).
As for the EOS, once the thresholds of NDVI(T) and NDVI(T)/NDVI max were defined from the multi-year averaged NDVI curve (i.e., step 1 to step 4; Supplementary Fig. 12) the determination of EOS a and EOS p for each year (step 5, based on the prescribed thresholds) is less likely to be influenced by the peak times or decrease speed of NDVI.For example, if NDVI max keeps constant (e.g., 0.8), NDVI(T) for EOS a is 0.2 and the NDVI(T)/NDVI max for EOS p is 25%.No matter how the peak time and decrease speed changes, the time when NDVI value decrease to 0.2 always coincides with the time when a relative decrease to 25% occurred (Supplementary Fig. 17).Therefore, we can conclude that if NDVI max stays the same, EOS a and EOS p should be equal.Difference in EOS a and EOS p can thus be used to indicate asynchronization of EOS resulted from the increasing NDVI max .

Fig. 1
Fig. 1 Schematic representation of NDVI variations during the autumn senescence process.In this study, autumn senescence starts at Peak Of Season (POS) when NDVI reaches its annual maximum (NDVI max ), and ends at the End Of Season (EOS).The monthly decrease in NDVI (ΔNDVI) is estimated as the NDVI difference between the current and the following months (in positive values), and the remaining NDVI is the NDVI left in the EOS month (in positive values).The sum of monthly ΔNDVI and the remaining NDVI is therefore equal to NDVI max .Detailed definitions of these items can be found in Supplementary Table1.

Fig. 2
Fig. 2 Trends in ΔNDVI over individual months and the entire senescence period from 1982 to 2018. a Temporal trends of ΔNDVI for senescence months in the Northern Hemisphere (NH) (recalling that a positive value of ΔNDVI corresponds to a decrease in NDVI).b Temporal trends of ΔNDVI over the senescence period (from July to October, JASO) and remaining NDVI at the end of the autumn senescence.c Temporal trends of the annual maximum NDVI (NDVI max ).The bars inserted in a indicate trends in ΔNDVI for each month, the whole senescence period, and the remaining NDVI.d-i Spatial patterns of ΔNDVI trends during 1982-2018 for senescence months (d-g), the entire senescence period (h), and remaining NDVI (i).Black dots in d-i indicate significant trends at 0.05 level, while the inserted stacked bars indicate the frequency distributions of ΔNDVI trends.

Fig. 3
Fig. 3 The dominant climatic factors of the trends in ΔNDVI and the remaining NDVI.a-e The distribution of dominant climatic factors was negatively correlated with ΔNDVI for the senescence months (a-d) and the remaining NDVI (e) in the climate space.f-j The distribution of dominant climatic factors was positively correlated with ΔNDVI for the senescence months (f-i) and the remaining NDVI (j) in the climate space.Each climate bin was defined by 0.5 °C intervals of mean annual temperature (MAT) and 20 mm intervals of mean annual precipitation (MAP).The factor that appeared most frequently for each bin was presented in these figures.k-o The percentages of pixels occupied by dominant climatic factors, which positively (left bars, add to 100%) and negatively (right bars, add to 100%) correlated with ΔNDVI for senescence months (k-n) and the remaining NDVI (o).Note that the positive correlation with ΔNDVI or the negative correlation with the remaining NDVI indicates that this factor can accelerate the speed of autumn senescence, and vice versa.

Fig. 4
Fig. 4 Schematic representation of the influence of greening on the speed of monthly senescence and the detection of EOS. a For each senescence month, monthly senescence is proxied by the decrease in NDVI (ΔNDVI) or percentage of decrease in NDVI relative to NDVI max (δNDVI = ΔNDVI/ NDVI max ).Peak of season (POS) greening would lead to an apparent increase in ΔNDVI even if δNDVI stays invariant.b In correspondence to ΔNDVI and δNDVI, end of season can be determined by fixed NDVI thresholds (EOS a ) or fixed percentage of NDVI max (EOS p ).If NDVI max were kept unchanged, these two methods would deliver the same EOS, whereas increases in NDVI max would results to a later EOS a than EOS p .

Fig. 5
Fig. 5 Mismatch between end-of-season (EOS) determined by fixed thresholds of NDVI (EOS a ) and fixed percentages of NDVI max (EOS p ). a, b Spatial patterns of trends in EOS a and EOS p in 1982-2018.Black dots indicate a significant trend at the 0.05 level.The stacked inserted bars indicate the frequency distributions of trends in EOS a and EOS p .c Difference between trends in EOS a and EOS p (trends in EOS a minus trends in EOS p ). d Contribution of trends in EOS p (green bars) and NDVI max (i.e., EOS a -EOS p , orange bars) to trends in EOS a (purple bars) for different vegetation types (see Methods).MF mixed forests, DBF deciduous broadleaf forests, DNF deciduous needleleaf forests.