Asymmetric temperature effect on leaf senescence and its control on ecosystem productivity Open Access

Abstract

Widespread autumn cooling occurred in the northern hemisphere (NH) during the period 2004–2018, primarily due to the strengthening of the Pacific Decadal Oscillation and Siberian High. Yet, while there has been considerable focus on the warming impacts, the effects of natural cooling on autumn leaf senescence and plant productivity have been largely overlooked. This gap in knowledge hinders our understanding of how vegetation adapts and acclimates to complex climate change. In this study, we utilize over 36,000 in situ phenological time series from 11,138 European sites dating back to the 1950s, and 30 years of satellite greenness data (1989–2018), to demonstrate that leaf senescence dates (LSD) in northern forests responded more strongly to warming than to cooling in autumn. Specifically, a 1 °C increase in temperature caused 7.5 ± 0.2 days' delay in LSD, whereas a 1 °C decrease led to an advance of LSD with 3.3 ± 0.1 days (P < 0.001). This asymmetry in temperature effects on LSD is attributed to greater preoverwintering plant-resource acquisition requirements, lower frost risk, and greater water availability under warming than cooling conditions. These differential LSD responses highlight the nonlinear impact of temperature on autumn plant productivity, which current process-oriented models fail to accurately capture. Our findings emphasize the need to account for the asymmetric effects of warming and cooling on leaf senescence in model projections and in understanding vegetation–climate feedback mechanisms.

Introduction

The timing of autumn leaf senescence plays a critical role in the regulation of carbon (C) and nutrient cycling in terrestrial ecosystems (1–5). Understanding the changes in leaf senescence dates (LSD) under climate change is therefore essential for the accurate prediction of future dynamics and feedbacks of global biogeochemical cycles (6). Existing studies have identified varied responses of LSD to climate warming, leading to divergent patterns of LSD (delayed, advanced, or stable) (7–9). Compared with spring phenology, a detailed understanding of autumn phenological dynamics remains lacking (10, 11), where the underlying mechanisms of complex LSD responses to temperature change are unclear (12).

Variability in LSD is driven by multiple climate and biotic factors, including temperature, precipitation, drought, photoperiod, wind speed, spring leaf-out date, and growing-season C assimilation (8, 13–20). Autumn temperatures exert direct effects on leaf physiological status, such as chlorophyll levels, photosynthesis, pigment degradation, and nutrient remobilization, ultimately triggering leaf senescence (21–23). Autumn temperatures may also exert indirect impacts on LSD through altering soil moisture and vapor pressure deficit (VPD) (24–26).

Although understanding of ecological effects of regional cooling under on-going global warming has increasingly gained attention (26–29), patterns of LSD responses to autumn warming and cooling scenarios remain a subject of ongoing debate that limits understanding of climate change impacts and maintains bias in forecasts of vegetation–climate interactions (30–32). For example, empirical field data show a stronger LSD response in European deciduous tree species to warmer than cooler autumn temperatures (30), whereas data from temperature manipulation experiments show larger LSD responses to cooling than to warming and a lack of difference between the two conditions (31, 32). Recent analyses indicate a substantial cooling trend in northern hemisphere (NH) autumn temperatures over the period 2004–2018 that is associated with the strengthening of Pacific Decadal Oscillation and Siberian High (26, 33); thus, long-term, large-scale data provide the opportunity to test for natural temperature-driven regulation and biophysical and biochemical drivers of autumn LSD and ecosystem C uptake.

To this end, we aim to quantify responses of autumn LSD and C update to natural warming and cooling and test for underlying mechanisms using ground and remotely sensed data, comprising a phenological time series of >36,000 records collected from 11,138 sites across Europe between 1951 and 2016 and estimates derived from normalized difference vegetation index (NDVI) data obtained between 1989 and 2018 covering the NH autumn cooling period at middle and high latitudes (>30°N). Specifically, we ask: (i) Does LSD in northern forests respond asymmetrically to natural autumn warming and cooling? (ii) What are the underlying mechanisms for these potential asymmetric patterns? (iii) Do asymmetrical temperature effects on LSD influence plant productivity responses?

LSD responses to natural autumn warming and cooling

At the species level, ridge and multiple linear regression analyses of in situ long-term records from the Pan European Phenology (PEP) database (34) (see Materials and methods) for four representative temperate tree species (Fig. 1A) indicated a tendency toward asymmetric LSD responses to natural autumn warming and cooling, evidenced by different changes in LSD caused by 1 °C increase. Apart from Aesculus hippocastanum L., LSD's sensitivities to warming were overall greater than that to cooling in autumn for other three species, i.e. Betula pendula Roth, Fagus sylvatica L., and Quercus robur L., using ridge regression (Fig. 1B) and multiple linear regression analyses (Fig. 1C).

For the biome-level remote-sensing analyses, we applied a natural autumn cooling period (2004–2018) to identify warming and cooling sites (see Materials and methods) and test for LSD response patterns (Fig. 2). Results showed that LSD was delayed at warming sites and became earlier at cooling sites for evergreen needleleaf, deciduous broadleaf, and mixed forests (Fig. 2B), and the greater overall biome ridge sensitivity of LSD to autumn warming (0.64 ± 0.008, mean ± SE) than to cooling (0.39 ± 0.007) (P < 0.001) reflected consistent asymmetric LSD responses for evergreen needleleaf (0.63 ± 0.008 vs. 0.57 ± 0.009; P < 0.001), deciduous broadleaf (0.68 ± 0.04 vs. 0.09 ± 0.02; P < 0.001), and mixed forests (0.71 ± 0.02 vs. 0.45 ± 0.008; P < 0.001) (Fig. 2D). Similarly, multiple linear regression analyses showed asymmetric patterns of biome LSD responses to warming and cooling (all forests: 7.5 days/°C ± 0.2 vs. 3.3 ± 0.1, P < 0.001; evergreen needleleaf: 8.2 ± 0.2 vs. 4.8 ± 0.2, P < 0.001; deciduous broadleaf: 4.4 ± 0.3 vs. 2.3 ± 0.2, P < 0.001; mixed forests: 4.6 ± 0.3 vs. 3.3 ± 0.1, P < 0.001) (Fig. S1). Although we compared the LSD responses in each biome, the warming and cooling sites were spatially separated from each other, which might introduce uncertainty. Then, we conducted a temporal analysis to complement the spatial analyses, based on grid cells in which there were shifts in warming and cooling between the period 1989–2003 and 2004–2018 (see Materials and methods; Fig. 2E), and similarly found greater sensitivity of LSD responses to autumn warming than to cooling among the forest biomes (Fig. 2F and Fig. S2).

In line with a previous study of field observations of warm and cold autumn conditions identified from a long-term average temperature threshold (30), our analyses of ground and remotely sensed data showed larger tree species and forest biome LSD responses to autumn warming than to cooling, indicating nonlinear temperature control of leaf senescence; thus, cooling effects on LSD should not be assumed to be the inverse of those to warming (35). In contrast to our findings, temperature manipulation experiments show larger or null LSD responses to autumn cooling than to warming (31, 32), possibly due to confounding effects of unmeasured and/or uncontrollable environmental factors that lead to an underestimation of phenological responses in warming experiments (12, 36, 37).

Mechanisms of asymmetric LSD responses to shifts in autumn temperature

We hypothesize that the asymmetric forest biome and species LSD responses to temperature shifts may be associated with differences in climatic forcing effects, such as cold degree days (CDD), adaptation to frost risk, and/or soil water availability (38, 39, 13). To test these hypotheses, we compared LSD responses to changes in CDD, preseason temperature variability and frost risk, and water stress, based on the Standardized Precipitation-Evapotranspiration Index (SPEI) at the autumn warming and cooling sites from 2004 to 2018 (see Materials and methods).

The first hypothesis accounts for changes in LSD response to exponential decreases in CDD with rising temperatures under warming and cooling conditions (Fig. S3), where sensitivities of LSD to CDD were larger under autumn warming than cooling, regardless of base temperature (5, 10, or 15 °C) for the calculation of CDD (P < 0.001) (Fig. 3A). Comparison of the relations between LSD and CDD under autumn warming and cooling and contrasting CDD base temperatures showed that LSD responses to CDD were larger under warming (P < 0.001) (Fig. 3B and Fig. S4), while random slope model estimates showed a larger LSD response to CDD under warming (P < 0.001), when effects of geographical location (latitude and longitude) were controlled (Fig. S5). Before leaf senescence, trees must absorb ample carbohydrate and nutrient stores to ensure tissue maturation, support overwintering and facilitate budburst in the subsequent spring (5, 30) and under favorable, autumn warming conditions, trees may delay leaf senescence to enhance C uptake (5, 30). Thus, our analyses indicate that under a given shift in CDD caused by warming or cooling, delays in leaf senescence caused by warming may be larger than advances caused by cooling, due to the asymmetry of LSD responses to changes in CDD.

Our second hypothesis relates to tree adaptations to frost damage (39–41), as indicated by studies that show spatial variation in phenology responses to temperature changes, with smaller responses observed at sites experiencing higher local temperature variability (40, 42). As expected, we found larger variations in temperature at cooling sites than sites that experienced warming, corresponding to lower LSD responses (Fig. 3C). Furthermore, the standardized frost risk index we developed to explore differences in frost risk between autumn warming and cooling conditions (See Materials and methods) showed greater risks of frost at cooling sites (Fig. 3D) that lead to conservative response strategies of leaf sensing to avoid frost damage, as indicated by decreased LSD responses under increasing risks of frost (Fig. 3E).

The third hypothesis concerns the potential effects of water stress on leaf senescence. We found that water stress levels were higher at the cooling sites than at the warming sites (Fig. 3F). This observation aligns with findings that water stress can significantly reduce leaf stomatal conductance and photosynthesis and enhance chlorophyll degradation (40, 43–45), potentially lowering the temperature sensitivity of leaf senescence at cooler locations. However, warming and more widespread droughts in the future (46–48) may increase the impact of droughts associated with warming on LSD (13). Consequently, this could result in a diminished LSD response to warming, akin to outcomes observed in temperature manipulation experiments that involve significant temperature variations (e.g. changes > 2 °C in warming and cooling scenarios) (32). This suggests a complex interplay between temperature, water availability, and their combined effects on leaf senescence, underlining the need for comprehensive studies to understand these dynamics.

Links between autumn temperature shifts, LSD responses, and plant productivity

Sensitivities of plant growth and productivity to temperature are highly variable across space and time and are regulated by a range of factors, such as vegetation type, prevailing climate, and temperature variability (49, 50). However, the lack of a mechanistic understanding of the drivers of the high levels of heterogeneity in sensitivities of plant productivity to temperature poses a challenge to the robustness of process-oriented models (50, 51). Indeed, we found the high levels of spatial heterogeneity in autumn productivity sensitivities to temperature across remotely sensed grid cells, as indicated by the climate space of mean annual precipitation (MAP) and mean annual temperature (MAT) (Fig. S6). Hence, this study analyzed the phenological temperature sensitivity of plant productivity and explored the contribution of asymmetric LSD responses to asymmetry in plant productivity under autumn warming and cooling. Four productivity datasets, comprising a kernel NDVI (kNDVI) and three gross primary productivity (GPP) datasets (Global Land Surface Satellite GPP, GLASS-GPP; two-leaf light use efficiency modeled GPP, TLLUE-GPP; near-infrared reflectance of vegetation-based GPP, NIRv-GPP), all lines of evidences confirmed a larger response of autumn productivity to increasing autumn temperatures than under cooling (Fig. 4A). Sensitivities of autumn productivity to temperature were positively correlated with LSD sensitivities to temperature, where LSD sensitivities accounted for 46% of spatial variability in autumn productivity sensitivities to temperature (Fig. 4B), while partial least-squares structural equation modeling (PLS-SEM) and random forest modeling analyses indicated that LSD sensitivity to temperature was the most important driver of autumn productivity sensitivity to temperature (Fig. 4C and Fig. S7). These results indicate that asymmetric LSD responses to autumn warming and cooling may lead to corresponding asymmetric responses in autumn productivity. If we do not account for the asymmetric productivity response to a 1 °C autumn temperature change (i.e. 1 °C warming or cooling), the biases in autumn plant productivity across the three studied forest biomes (i.e. evergreen needleleaf, deciduous broadleaf, and mixed forests) in the NH would amount to nearly 0.031 ± 0.007 Pg C yr⁻¹ (Fig. S8; see Materials and methods), which is half the minimum value within the range of carbon loss from deforestation in the Amazon forest (0.06 to 0.21 Pg C yr⁻¹) (52). As the magnitude of warming and cooling increases (e.g. 2 °C), the biases would increase further, indicating that asymmetric LSD and GPP responses to autumn warming and cooling significantly impact global terrestrial carbon budget.

We further analyzed sensitivities of autumn productivity to temperature under warming and cooling conditions using state-of-the-art Dynamic Global Vegetation Models (DGVMs) (see Materials and methods) and found that 11 out of 16 models showed larger responses of autumn productivity to warming compared with cooling (Fig. S9). The sensitivity comparison of the combined results from 16 models reveals a stronger, though not statistically significant, response to warming than to cooling (P = 0.13; Fig. 4D). These findings suggest that current state-of-the-art models partially but limitedly capture the observed differences in responses of autumn productivity to temperature changes under warming and cooling. One potential reason could be related to the complexity of phenological control and its climatic responses. Incorporating a detailed characterization of the asymmetric LSD responses may enhance the representation of asymmetric productivity responses to autumn warming and cooling in process-oriented models, thereby improving the accuracy of seasonal vegetation dynamics predictions.

Conclusions

The phenological responses of trees to changing temperature will strongly regulate terrestrial carbon storage under climate change. By exploring the plant responses to warming vs. cooling, our study allows us to test for the presence of acclimation, which will underpin the extent and duration of phenological shifts under climate warming. Both in situ observations and satellite data suggest that the timing of LSD in northern forest biomes responds asymmetrically to natural variations in autumn temperatures, showing a more pronounced response to warming than to cooling. This disparity in LSD responses may be attributed to factors such as the greater sensitivity of LSD to CDD, lower risks of frost, and greater water availability under warming than cooling. Such asymmetric LSD reactions may, in turn, lead to unequal changes in autumn plant productivity in response to warming and cooling, a complexity that current DGVMs struggle to accurately represent. This new understanding is crucial for refining vegetation and climate models, thereby improving the accuracy of future carbon cycle and climate projections.

Materials and methods

In situ LSD observation

We used ground LSD observations, the date of 50% autumnal coloring of tree leaves (BBCH code 94), across Europe from the long-term plant phenological observation database, the PEP Project (PEP725, http://www.pep725.eu/) (34). We collected all available phenological time series (>36,000) from 11,138 study sites since the 1950s for four dominant tree species, comprising A. hippocastanum L., B. pendula Roth, F. sylvatica L., and Q. robur L., which had sufficient autumn warming and cooling samples (Fig. 1A). The median absolute deviation (MAD) method was applied to identify and remove erroneous data points in PEP725 records (30). For a time series of LSD_t1, LSD_t2, …, LSD_t at each site, the MAD is calculated as follows:

The data point exceeding 2.5 times the MAD was deemed an outlier and consequently eliminated from the LSD time series before analysis.

Remotely sensed greenness-derived LSD

We used the bi-weekly PKU Global Inventory Modeling and Mapping Studies (GIMMS) NDVI, at 1/12° spatial resolution, to estimate LSD between 1989 and 2018. The PKU GIMMS NDVI was produced utilizing a biome-specific back-propagation neural network algorithm that using GIMMS NDVI3g product and Landsat NDVI samples, and its temporal coverage was extended to 2022 following consolidation with the Moderate-Resolution Imaging Spectroradiometer (MODIS) NDVI using the random forest method. This new NDVI product efficiently eliminated the effects of sensor degradation and satellite orbital drift and showed high overall accuracy assessed by Landsat NDVI samples. To mitigate snow effects, we replaced all compromised NDVI values with the winter (December–February) mean of snow-free NDVI values across all years and then the NDVI time series was reconstructed using a Savitzky–Golay filter to eliminate abnormal values (53, 54); sparse vegetation cover was removed by eliminating grid cells with mean annual NDVI value < 0.1.

We estimated LSD using the double-logistic function and dynamic-threshold approaches, which were widely used methods for estimating phenological dates from remotely sensed data (13, 15, 28). A double-logistic function was fitted to the NDVI time series, followed by calculation of the second-order derivative for the fitted curve (55, 56), and LSD was determined as the date corresponding to the second local maximum in the second half year:

where y(t) is the NDVI value at the day of year (DOY), t; a is the background NDVI value; and, b–e are the double-logistic function parameters.

For the dynamic-threshold approach, we computed NDVI ratio (⁠$\mathrm{Rati}{\mathrm{o}}_{\mathrm{day}}$⁠) annually for each NDVI time series as follows (57):

where $\mathrm{NDV}{\mathrm{I}}_{\min }$ and $\mathrm{NDV}{\mathrm{I}}_{\max }$ are the minimum and maximum NDVI of each year, respectively, and $\mathrm{NDV}{\mathrm{I}}_{\mathrm{day}}$ is the daily NDVI. LSD is determined as the time (DOY) when $\mathrm{Rati}{\mathrm{o}}_{\mathrm{day}}$ decreased to 0.5 (57). The average LSD of the two approaches was calculated as the ultimate LSD.

Plant productivity data

We used two types of plant productivity indicator, comprising the newly developed kNDVI and GPP; kNDVI is a widely used indicator of GPP (58) and it was calculated based on a simplified operational index version, which was expressed as $\mathrm{kNDVI}=\tanh (\mathrm{NDV}{\mathrm{I}}^2)$⁠, using the monthly MODIS MOD13C2 v061 product (59). We used three independent GPP products to minimize uncertainty stemming from a single product: the eight-day GLASS GPP (60), the monthly TLLUE modeled GPP (61); and the monthly NIRv based GPP (62). All kNDVI/GPP data were resampled to 1/12° to match the spatial resolution of the PKU GIMMS NDVI.

Process-oriented modeled GPP

We employed modeled GPP data from 16 state-of-the-art DGVMs that participated in the TRENDY (Trends and drivers of the regional scale sources and sinks of carbon dioxide) project version 11, comprising CABLE-POP, CLM5.0, CLASSIC, DLEM, IBIS, ISBA-CTRIP, ISAM, JSBACH, LPJ-GUESS, LPX-Bern, JULES, ORCHIDEE, OCN, SDGVM, VISIT-NIES, and YIBs (63). We used GPP data from the S3 simulations that are based on all time-varying forcings of CO₂, climate, nitrogen deposition, and land use.

Climate data

We used daily European gridded observational (E-OBS v27.0e) climate data at 0.1° spatial resolution (temperature, precipitation, and radiation) produced by ECA & D (European Climate Assessment & Dataset) project (64) in species-scale analysis, with PEP725 data to calculate temperature changes and determine temperature-relevant optimal preseason length (see Analyses), and daily Multi-Source Weather (MSWX) climate data at 0.1° spatial resolution (65) were used in the biome-scale analysis, with satellite-derived LSD; VPD was calculated using “RHtoVPD” function of “plantecophys” package in R language, based on MSWX climate data (monthly temperature, atmospheric pressure, and relative humidity) (66). We obtained the monthly SPEI dataset (three-month scale) from Consejo Superior de Investigaciones Científicas to calculate soil water availability changes (67).

Cold degree day

We computed cold day degrees to assess climate forcing under warming and cooling (38):

where $\mathrm{CD}{\mathrm{D}}_d$ represents the CDD for date d and $T_{\mathrm{daily}}$ is the daily mean temperature; base temperature (⁠$T_{\mathrm{base}}$⁠) was set 5, 10, and 15 °C (13) and $\mathrm{CD}{\mathrm{D}}_d=0$ when $T_{\mathrm{daily}}$ > $T_{\mathrm{base}}$⁠. We computed the accumulated CDD during the period 1st July (⁠$d_0$⁠) to LSD (15), using the daily MSWX temperature dataset.

Analyses

We analyzed ground and remotely sensed based data at forest tree species and biome scales, respectively. For ground-based observations, we utilized PEP725 LSD and E-OBS climate data to identify autumn warming and cooling samples for each tree species, where we selected the highest coefficient from a partial correlation analysis (P < 0.05) to determine the optimal preseason period, based on the timing of the greatest impact of temperature on LSD, whilst accounting for the effects of radiation and precipitation; temperature effects on LSD were analyzed for 8 to 120 d before mean LSD, with a step of 8 d. Then, we calculated LSD responses to temperature changes during the optimal preseason using ridge regression and multiple linear regression analyses, while controlling for effects of precipitation and radiation. For ridge regression analysis, we calculated and normalized (0–1) anomalies for each variable, prior to estimation of sensitivities. We used linear least-squares regression, utilizing year as an independent variable, to calculate trends in autumn temperature (September to November) at P < 0.05. Given the diverse temporal spans of ground observations, we utilized a 15-year moving window to detect autumn warming and cooling periods within each LSD time series for individual tree species; then, we calculated species average LSD responses to autumn warming and cooling, based on site mean values, and differences in response to warming and cooling were analyzed using Student's t test method at P < 0.05. Process framework for identifying warming and cooling samples can be found in Fig. S10.

For analysis of remotely sensed data, we applied the NH autumn cooling period (2004–2018) to detect autumn warming and cooling grid cells for forest biome types obtained from the MODIS MCD12C1 IGBP land cover product. Due to a lack of warming grid cells for evergreen broadleaf and deciduous needleleaf forests, analyses of LSD responses were restricted to evergreen needleleaf, deciduous broadleaf, and mixed forests. Differences in mean biome LSD responses between warming and cooling were analyzed utilizing Student's t test method at P < 0.05. We analyzed differences in LSD responses (Student's t test at P < 0.05) to shifts in warming and cooling between the periods 1989–2003 and 2004–2018, based on grid cells in which warming was followed by cooling or vice versa. We found no evidence for effects of day length on LSD responses to shifts in autumn temperature (Fig. 2D and Figs. S1 and S11).

We tested for drivers of shape of LSD responses to shifts in temperature, by analyzing sensitivities of LSD to CDD, frost risk, and water availability between autumn warming and cooling conditions during the period 2004–2018. Differences in the sensitivities of LSD to CDD were tested using linear mixed models (30), covariance analysis (68, 69), and random slope models (70). The linear mixed model analyses were based on linear regression analysis of site-based warming and cooling sensitivities of LSD to CDD that were then compared with direction of autumn temperature change as a fixed effect and random intercepts among forest biomes. Covariance analysis comprised a comparison of fitted regression slopes between LSD and CDD under autumn warming and cooling conditions, while a random slope model compared temporal LSD–CDD regression slopes between autumn warming and cooling conditions, with grid cell latitude and longitude as random factors (70).

We defined a frost risk index (0–1) that was standardized across the number of frost days in the studied warming and cooling grid cells, based on daily minimum temperature <0 °C during the preseason period (71, 72), and we calculated mean SPEI (three-month scale) during the preseason period. Linear mixed models were used to test for differences in frost risk and mean SPEI between autumn warming and cooling conditions.

Autumn plant productivity was estimated from mean kNDVI and accumulated GPP from September to November (26, 28) and sensitivities of autumn productivity to autumn temperature were calculated using multiple linear regressions, controlling for effects of radiation and precipitation. PLS–SEM and random forest models were used to examine effects of sensitivity of LSD to temperature S(LSD|T) on sensitivities of autumn productivity (i.e. kNDVI) to temperature S(kNDVI|T) and differences in autumn GPP sensitivities to warming and cooling estimated by 16 DGVMs were analyzed using Student's t test at P < 0.05. We quantified the extent to which the asymmetry GPP responses impact carbon budget, based on the asymmetry responses derived from three GPP datasets (i.e. GLASS-GPP, TLLUE-GPP, and NIRv-GPP) and the coverage of the three studied forest biomes in the NH (i.e. evergreen needleleaf, deciduous broadleaf, and mixed forests) using MODIS MCD12C1 land cover dataset (Fig. 4A and Fig. S8). The bias when not accounting for the asymmetric productivity response was calculated based on the product of the forest area and the difference in GPP responses to 1 °C warming and cooling (Fig. 4A and Fig. S8).

Supplementary Material

Supplementary material is available at PNAS Nexus online.

Funding

This work is supported by the National Natural Science Foundation of China (41921001).

Author Contributions

L.H. and Z.-L.L. designed this study. L.H. carried out the analysis and drafted the initial manuscript. J.W., J.P., C.M.Z., T.W.C., Y.F., W.Z., J.X., Z.-L.L., and X.W. provided critical insights on results interpretation and text editing. All authors contributed to the interpretation of the results and approved the final manuscript.

Data Availability

All study data are public and included in SI appendix. The specific links for data used in this study can be found in SI Appendix, Table S1. The codes used for data analysis are available on Zenodo at https://zenodo.org/doi/10.5281/zenodo.10677097.

References

Author notes