Non-symmetric responses of leaf onset date to natural warming and cooling in northern ecosystems

Abstract The northern hemisphere has experienced regional cooling, especially during the global warming hiatus (1998–2012) due to ocean energy redistribution. However, the lack of studies about the natural cooling effects hampers our understanding of vegetation responses to climate change. Using 15,125 ground phenological time series at 3,620 sites since the 1950s and 31-year satellite greenness observations (1982–2012) covering the warming hiatus period, we show a stronger response of leaf onset date (LOD) to natural cooling than to warming, i.e. the delay of LOD caused by 1°C cooling is larger than the advance of LOD with 1°C warming. This might be because cooling leads to larger chilling accumulation and heating requirements for leaf onset, but this non-symmetric LOD response is partially offset by warming-related drying. Moreover, spring greening magnitude, in terms of satellite-based greenness and productivity, is more sensitive to LOD changes in the warming area than in the cooling. These results highlight the importance of considering non-symmetric responses of spring greening to warming and cooling when predicting vegetation-climate feedbacks.


Introduction
Spring leaf onset date (LOD) has advanced in recent decades in northern mid to high latitudes (>30°N) under global warming (1)(2)(3)(4)(5)(6)(7).This advance is highly sensitive to temperature changes, extends the growing season length, and accordingly increases the carbon uptake of terrestrial ecosystems (1,8,9).Understanding the responses of LOD to temperature changes in terms of sign and magnitude is therefore crucial for assessing the influence of climate change on terrestrial ecosystems and its feedback to climate (10,11).Unlike warming effects, however, most existing studies ignore the impact of cooling anomalies on LOD, which may cause biased predictions of vegetation-climate feedbacks (12).
In northern regions, winter and spring temperatures are generally considered the principal drivers of spring LOD.Trees need to accumulate enough winter chilling to end the endodormancy phase and enough spring warming to break the ecodormancy phase, further triggering plant leaf onset (13)(14)(15)(16).The earlier emergence of spring leaves has been associated with warmer temperatures because of easily reaching heating demand (14,17).However, the global mean temperature has not always shown a steady increase, with the global warming hiatus observed between 1998 and 2012 possibly due to an energy redistribution within the oceans (18)(19)(20)(21).Until now, we only know the relative responses of LOD to warming and cooling for some species, e.g.tree saplings and grass, in field experiments (12,22,23).Warming of 1°C in winter/spring led to an advance of 8.8 days in budburst dates of Fagus sylvatica L., whereas 1°C cooling delayed it by 10.9 days (22).Two manipulative experiments in the Tibetan Plateau showed non-significant differences in sensitivities to warming and cooling for grass leaf-out (12,23).Both field-and ecosystemscale analyses have mainly focused on advancing effects of natural warming on LOD, influenced by photoperiod (24,25), precipitation amount (26,27) and frequency (28), and soil water (29) and nutrient availability (30).Besides direct physiological effects, temperature changes regulate ecosystem composition and function, soil moisture, snowmelt, and permafrost degradation in high altitudes and latitudes, which may affect plant spring growth (31)(32)(33)(34)(35).There is limited evidence of non-symmetric or symmetric responses of LOD to natural warming and cooling at the species to ecosystem scales, especially in mature woody biomes (23).The impacts of spring greening timing (i.e.LOD) on spring greening magnitudes (i.e.spring greenness and productivity) during the warming hiatus are also unclear.Therefore, we ask two questions: (i) does spring LOD respond to warming and cooling symmetrically in northern biomes?and (ii) what are the physical and physiological mechanisms related to non-symmetric or symmetric patterns?To this end, we investigated the responses of LOD to natural warming and cooling, using gridded meteorological data (temperature, precipitation, and cloud cover) together with LOD data from two independent data sets: (i) 15,125 ground phenological time series at 3,620 sites across Europe since the 1950s and (ii) normalized difference vegetation index (NDVI) data, for northern mid to high latitudes (>30°N) from 1982 to 2012.

Responses of LOD to natural warming and cooling
We used long-term in situ LOD observations of seven European dominant tree species (5), derived from the PEP725 database that provides the longest and most comprehensive phenological records, to study the LOD responses to warming and cooling (Fig. 1) (see Methods).The results overall indicated the non-symmetric LOD responses to warming and cooling, i.e. five out of seven species (Tilla cordata P < 0.05, Fagus sylvatica P < 0.05, Betula pendula P < 0.05, Alnus glutinosa P < 0.01, and Aesculus hippocastanum P < 0.001) were more sensitive to cooling than to warming (Fig. 1B).Consistent with previous field experiment for Fagus sylvatica (22), our results highlight the stronger responses of LOD to natural cooling than to warming, indicating a nonlinear temperature control of LOD.
To quantify the variations of non-symmetric LOD response, we defined a LOD non-symmetric index calculated by the difference between the LOD sensitivities to warming and cooling; the positive value indicates stronger sensitivity to cooling, while the negative value suggests stronger sensitivity to warming (see Methods).Using the LOD simulated by a growing-degree-day (GDD) algorithm (see Methods), we also found that LOD was more sensitive to cooling than to warming during the warming hiatus (Fig. S3).For future projections (2016-99), we also found a non-systematic response but with different patterns, i.e.LOD becomes more sensitive to warming than to cooling under the scenarios of highest baseline of carbon emissions (RCP8.5)(Fig. S4).

The mechanisms under non-symmetric LOD responses
Exploring physiological mechanisms under LOD responses is challenging.Here, we hypothesized that the non-symmetric LOD responses to warming and cooling might be related to changes in chilling accumulation (CA, the amount of chilling during endodormancy), heat requirement (HR, the accumulated forcing temperature required for leaf onset), and water availability (16,(36)(37)(38).To test these hypotheses, we determined the changes in CA, HR, and water stress using a drought index (the Standardized Precipitation Evapotranspiration Index, SPEI) in warming and cooling areas during the warming hiatus (see Methods).First, we confirmed a dual role of temperature in controlling LOD variations with an exponential decay-like relationship between chilling days and forcing degrees (16) (Fig. 3A).Grouping grid cells into warming and cooling indicated that warming reduced CA and HR while cooling increased the two variables (Fig. 3B, C; Figs.S5-S6).Trees with more CA in the phase of endodormancy might need more HR to break ecodormancy for reactivating growth (13,39).Expectedly, cooling grids showed more changes in HR caused than warming grids both in terms of trees (forests) and low vegetation (shrublands, savannas, woody savannas, and grasslands) (Fig. 3C; Figs.S6-S7).Non-symmetric changes in CA and HR caused by warming and cooling might follow the non-symmetric LOD responses.
On the other hand, we found that warming and cooling are associated with soil water availability changes (i.e.SPEI trend), further affecting LOD when controlling the effects of precipitation and radiation (Fig. S8).Decreased SPEI generally accompanies abundant sunshine in the warming areas, and these processes together lead to earlier LOD (40-43) (Fig. S8).In the current climate, the severity of preseason drying may not reach a turning point that could cause a delaying effect on LOD (40).Before the turning point, the elevated preseason temperature and radiation in drought may advance LOD (3,40,44).In contrast, cooling benefited maintaining soil water availability (Fig. S8A), offsetting the advancing effect caused by drought stress (40,41) and leading to delayed LOD (Fig. S8B).We found stronger effects of warming on water availability than effects of cooling (Fig. 3D).Considering the opposite and non-symmetric effects on soil water availability, the non-symmetric LOD responses to warming and cooling might be partially offset.For future projections (2016-99), we found a reversion of warming and cooling effect sizes, that is LOD will be more sensitive to warming than to cooling under RCP8.5 (Fig. S4B).Apart from the projection uncertainty caused by models and datasets, we proposed two potential reasons.First, future climate change may alter current non-symmetric patterns of chilling accumulation and heating requirements under warming and cooling, especially with temperature increases and precipitation variations (3,28).Second, warming-related drying stress might adjust climatic responses leading to a warming-dominant control on spring plant growth (45)(46)(47).

Connections among temperature change, LOD, and spring greening magnitude
Spring (from March to May) accumulated gross primary productivity (GPP) and mean NDVI were used as proxies of spring greening magnitude during the warming hiatus.Spring greening magnitudes were negatively correlated with spring greening timing (i.e.LOD) (Fig. 4), which suggested that earlier LOD would increase plant carbon uptake in spring (40,48,49).However, these relationships were significantly different in the warming and cooling areas by using covariance analysis (P < 0.001) (50).We also used the random slope model to present the relationship between spring GPP/NDVI and LOD in the temporal scale when controlling for latitudes and longitudes of grid cells as random factors (Fig. S9).To reduce the uncertainty brought by the fixed period used (i.e. from March to May), we calculated the spring greening magnitude with accumulated GPP and mean NDVI during the period from LOD to the maturity (i.e. the date corresponding to the maximum NDVI in the GIMMS NDVI3g time series) and obtained the similar results with fixed period used (Fig. S10).As confirmed by these independent lines of evidence, spring GPP/NDVI was more sensitive to LOD in the warming areas than in the cooling (Fig. 4; Figs.S9-S10), suggesting that the increase in spring plant productivity caused by 1-day LOD advance by warming was overall greater than the decrease by 1-day LOD delay by cooling.These results call for caution concerning model-based climatic responses of GPP and aid in understanding vegetation-climate feedbacks.

Conclusions
Using both ground records and satellite observations, we found non-symmetric LOD responses to natural warming and cooling, i.e. the LOD of northern biomes exhibited stronger responses to cooling than to warming.The underlying mechanism might be The four GPP/NDVI datasets all showed that the differences in the slopes between warming and cooling conditions were significant (P < 0.001) by using covariance analysis.Significance code for differences: ***P < 0.001.

He et al. | 5
associated with stronger variations of CA and HR by cooling, which could be partially offset by warming-associated drying.Moreover, the spring greening magnitude was more sensitive to LOD changes in the warming areas than in the cooling.Our findings provide a new conceptual framework of LOD responses to climate change, which is enlightening for model improvements and projections.

In situ LOD observation
We applied the in situ LOD observations from the Pan European Phenology Project (PEP725), which is an open database that contains long-term plant phenological observations from 25 European countries (http://www.pep725.eu/)(51).The date of the first visible foliar stalk for tree species (BBCH code 11) was used.All available records (15,125) from 1951 to 2018 were collected from 3,620 sites for seven European dominant tree species (5), i.e.Aesculus hippocastanum, Alnus glutinosa, Betula pendula, Fagus sylvatica, Fraxinus excelsior, Quercus robur, and Tilia cordata (Fig. 1A).

Satellite greenness-based LOD
We used Global Inventory Modeling and Mapping Studies (GIMMS) NDVI3g v1 data to derive LOD between 1982 and 2012 (52).The NDVI3g v1 data are derived from optical surface reflectance measurements taken by a series of NOAA-AVHRR satellites.
Corrections for intersensor calibration, orbital drifts, and stratospheric aerosols from volcanic eruptions have made it the most consistent long-term satellite vegetation dataset currently available (52,53).To remove snow effects, we replaced all contaminated NDVI with the mean of snow-free NDVI values from all years in winter (December-February) (54).A modified Savitzky-Golay filter was then used to eliminate abnormal values and reconstruct the NDVI time series (55).Furthermore, we eliminated sparse vegetation by removing grids with a mean annual NDVI value of less than 0.1.We applied two methods, i.e. the dynamic threshold approach and the double-logistic function, to estimate LOD to minimize the uncertainty from a single method.The two methods show similar results, so we calculated the average LOD from the two methods as the final LOD.
In the first method, we calculated NDVI ratios annually for each pixel as follows: where NDVI day is daily NDVI and NDVI min and NDVI max are the minimum and maximum NDVI of each year, respectively.A threshold ratio of 0.5 was used to determine LOD.
In the second method, we fitted the NDVI time series with a double-logistic function and then calculated the second-order derivative of the fitted curve.LOD was defined as the time when the rate of change in curvature reached its first local maximum in spring.
where t is time in days and y(t) is the NDVI value at time t. a is the initial background NDVI value, and b − e are parameters of this function.

Spring greening magnitude
We used mean NDVI and accumulated GPP from March to May as spring greening magnitude.The GIMMS NDVI3g v1 dataset was used to calculate spring mean NDVI.We used GPP data from three independent sources, i.e. the monthly satellite-based nearinfrared reflectance (NIRv) GPP, the daily light response function (LRF) GPP, and the 8-day two-leaf light use efficiency model (TL-LUE) GPP datasets.The NIRv GPP dataset has good performance at capturing seasonal and interannual variations of terrestrial GPP at a global scale (56).LRF GPP was estimated by an ecosystem-level physiological method using an asymptotic light response function between incoming photosynthetically active radiation (PAR) and GPP, which well represents the response observed at high spatiotemporal resolutions (57).TL-LUE GPP dataset distinguished GPP of sunlit and shaded leaves, suitable for studying the variations in seasonal cycles of GPP over many years (58).

Climatic data
In satellite-derived analysis, we used Multi-Source Weather (MSWX) temperature, precipitation, and radiation data with daily temporal resolution and 0.1° spatial resolution (59) in partial correlation analysis to determine optimal preseason length and the site-specific period before LOD with the highest absolute partialcorrelation coefficient (see Analyses).European gridded observational (E-OBS v23.1e) daily climate data with a spatial resolution of 0.1° were used in the phenological observation analysis at the species scale with PEP725 data.This dataset was provided by the ECA&D (European Climate Assessment & Dataset) project (60).
The daily temperature of the GFDL-ESM2M model in ISIMIP2b (Inter-Sectoral Impact Model Intercomparison Project 2b simulation round) was used to predict LOD from 2016 to 2099 under future scenarios (RCP4.5 and RCP8.5), and surface downwelling shortwave radiation and precipitation data were used to calculate the LOD responses to warming and cooling.

Chilling and forcing models
We used 10 chilling models to measure the number of chilling days (i.e.CA) and 8 GDD models for estimating the heating requirement (HR) for the LOD.Chilling models 1-8, 11, and 12 and GDD models 1-8 in (16) were used in our study.

Model for predicting LOD
To predict future LOD, we used a two-phase parallel model (PM) with optimal parameters (61).In contrast to the one-phase models (e.g.growing degree days [GDD] and spring warming [SW] models) that focus only on forcing accumulation, PM assumes that the forcing accumulation cannot begin until a critical threshold (C crit ) of the chilling state (S c , daily sum of chilling rates) is reached (62).The first phase of PM is chilling accumulation.A triangle function (Eq. 3) was used to describe the daily rate of chilling (R c ) (63), and S c began to accumulate after September 1 st of the preceding year (t c ) (Eq. 4): where R c is the daily rate of chilling.T is the daily mean temperature.T a , T b , and T c are three model parameters.S c is the daily sum of chilling rates and begins to accumulate after September 1 st of the preceding year (t c ).
The second phase of PM is forcing accumulation, and the day that the state of force (S f ) achieved its critical value (F crit ) was used to determine the modeled LOD.
where R f is the daily rate of forcing and starts from January 1 st of the current year (t 0 ).T d is a temperature threshold to establish the requirement for beginning the accumulation of forcing (Eq.6) and fulfilling C crit .A f , alpha, beta, F crit , and C crit are model parameters.K is an adjustment factor to ensure that the accumulation of forcing occurs after the chilling state (C crit ) is fulfilled.K min is another model parameter that determines the minimum potential of an unchilled bud to respond to the forcing temperature (63).Finally, the date when S f exceeds F crit is regarded as the LOD (Eq.8).PM parameters, including A f , alpha, beta, F crit , C crit , T a T b , T c , T d , and K min , were calibrated optimally by implementing the particle swarm optimization (PSO) algorithm (SPSO-2011) at each pixel, based on 31 years of satellite-derived LOD (1982-2012) and gridded air temperature data with daily scale.The set of optimal parameters was employed when the RMSE value between the observed and modeled LOD was the lowest.With the optimal parameters, we used PM to predict future LOD under scenarios RCP4.5 and RCP8.5.It should be noted that the uncertainty of future climatic projections might, to some degree, undermine the robustness of future LOD and its responses to warming and cooling.

Analyses
We performed ground-and satellite-based analyses at the species and biome scales, respectively.We used seven dominant tree species with long-term in situ observations of LOD in Europe from the PEP725 database.The biome types were obtained based on MCD12C1 land cover product (64).
To identify warming and cooling periods/grid cells, PEP725 phenology and E-OBS climate data were applied for ground-based analysis.The satellite-derived LOD and MSWX climate data were employed for the biome-based analysis.The relevant periods for preseason temperature impacts on phenology vary among biomes, species, and locations (37).To determine the optimal preseason during which average temperature had the largest influence on phenology, we computed the partial correlation coefficients between average temperature and LOD, controlling the effects of precipitation and radiation, from 0 to 6 months before the mean LOD with a step of 8 days (65).The optimal preseason length was the period with the highest absolute partial correlation coefficient.We then calculated the temperature trend during the optimal preseason length by linear least-squares regression analysis with year as the independent variable.This study used statistical significance at a 0.05 level for partial correlation and trend analysis.Due to the different time lengths of PEP725 records, we used a 15-year moving window to obtain warming and cooling periods in the time series of each station.In the satellite-based analysis, we applied the warming hiatus period (1998-2012) to identify warming and cooling grid cells.In addition, we identified grid cells with warming during 1982-98 and cooling during 1998-2012 for temporal analysis.For future projections, we determined the temperature-relevant preseason and computed the temperature changes during preseason for all grids within a moving window of 15 years and then identified the warming and cooling grids in each moving window from 2016 to 2099 based on ISIMIP2b climatic datasets.Finally, we calculated the non-symmetric index within a moving window of 15 years (Fig. S4).
Due to nonlinear temperature responses, the LOD responses to temperature changes were computed by log-log regression to avoid potential statistical artifacts using the linear method (66,67).In the ground-based analysis, we calculated the average value of the LOD responses to warming and cooling in each station and then the average value of the LOD responses to warming and cooling in each species.In the satellite-based analysis, we calculated the average value of the LOD responses to warming and cooling for each biome at the spatial (warming vs. cooling during 1998-2012) and temporal (warming during 1982-98 vs. cooling during 1998-2012) scales.Finally, the non-symmetric LOD response to warming and cooling was determined by Student's t-test method (at least P < 0.05).
We compared the changes in the number of chilling days, accumulated forcing degrees, and a water indicator at the warming and cooling grids during 1998-2012 to explore the mechanisms under the non-symmetric/symmetric LOD responses to warming and cooling.To calculate the water indicator, we employed monthly SPEI data at a spatial resolution of 0.5° from the SPEI base v. 2.7 at Consejo Superior de Investigaciones Científicas (CSIC) (68).The SPEI data consisted of multiscale monthly SPEI from 1 to 48 months; we selected the 3-month SPEI to capture the short-term water deficit (69).We calculated trends of chilling days, forcing degrees, and SPEI in the relevant periods for preseason temperature by the linear least-squares regression method.Then, we used Student's t-test method to check whether non-symmetric patterns exist.
We applied GPP and NDVI data to compute spring greening magnitude between March and May from 1998 to 2012.The regressions between spring greening magnitude and LOD in the warming and cooling conditions were created by the least-squares linear regression method.Then, we tested the statistical significance of the difference in the slopes of GPP/NDVI-LOD regressions between warming and cooling conditions by covariance analysis based on a procedure in (50).Finally, we applied a random slope model ("lme4" package in R4.2.0) to compare GPP/NDVI-LOD regressions in warming and cooling scenarios at the temporal scale when grid cells' latitudes and longitudes were used as random factors.

Fig. 1 .
Fig. 1.The distribution of PEP725 sites and comparisons of LOD responses to warming and cooling at the species scale.A) The locations of PEP725 sites for long-term in situ LOD observations of seven European dominant tree species since the 1950s.B) LOD responses (log-transfer, see Methods) to warming and cooling at the species scale from PEP725 ground observation data.The warming and cooling samples were obtained by using P < 0.05 for temperature changes and partial correlation analysis.The bar represents the standard error.Student's t-test was used to test the significance of the difference between the warming and cooling conditions.Significance code for differences: ***P < 0.001, **P < 0.01, and *P < 0.05.

Fig. 2 .
Fig. 2. Warming and cooling grid cells used in this study and comparisons of LOD responses to warming and cooling at the biome scale.A) Locations of warming and cooling grid cells during the warming hiatus (1998-2012).Grid cells with P < 0.05 for temperature changes and partial correlation analysis were retained (see Methods).B) The biome types of grid cells.C) Trends in LOD in warming and cooling areas for biomes from 1998 to 2012.D) LOD responses (log-transfer) to warming and cooling at the biome scale from satellite-based LOD data.The bar represents the standard error.Student's t-test was used to test the significance of the difference between the warming and cooling conditions.Significance code for differences: ***P < 0.001, **P < 0.01, and *P < 0.05.

Fig. 3 .
Fig. 3. Comparisons of changes for chilling accumulation (CA), heat requirement (HR), and water availability indicator (i.e.SPEI) in warming and cooling areas obtained from satellite-based analysis during the warming hiatus.A) Relationship between CA and HR.The green line indicates an exponential decay regression fitted using CA and HR.B, C, and D) present the trends of CA, HR, and SPEI in the warming and cooling areas during 1998-2012, respectively.The bar represents the standard error.Student's t-test was used to test the significance of the difference in the absolute values of trends for CA, HR, and SPEI in the warming and cooling areas.Significance code for differences: ***P < 0.001 and **P < 0.01.

Fig. 4 .
Fig.4.Connections among temperature change, LOD, and spring greening magnitude.The comparison for regressions between spring GPP/NDVI and LOD in the warming and cooling areas using the satellite-based NIRv GPP A), the LRF GPP B), TL-LUE GPP C), and GIMMS NDVI D), respectively.The spring-accumulated GPP and mean NDVI were calculated during the period from March to May.The red and blue kernel density plots represent the density distribution of warming and cooling grids in GPP/NDVI-LOD space, respectively.The four GPP/NDVI datasets all showed that the differences in the slopes between warming and cooling conditions were significant (P < 0.001) by using covariance analysis.Significance code for differences: ***P < 0.001.