Five deciduous broadleaved species of two trees (Fraxinus chinensis and Salix babylonica) and three shrubs (Jasminum nudiflorum, Forsythia suspensa, and Syringa oblata), widely distributed in temperate regions of China, were selected for the study. The sampling site was located in the Olympic Forest Park in northern Beijing (40°01′3.00″N, 116°23′2.98″E). Meteorological data (daily minimum and maximum temperatures and precipitation) for the period 1951–2021 in the study area were obtained from the China Meteorological Data Service Center.¹ The annual mean temperature and annual total precipitation in Beijing are 12.47°C and 586.92 mm (average for 1981–2010), respectively, and a trend for significant warming in the study area over the past 70 years is apparent (Figure 1A).

Twigs of each species were collected from 8 October, 2018 to 26 March, 2019. To assess the seasonal change in dormancy depth, twigs were sampled every 3–4 days, resulting in 42 collections in total. Because the pre-treatment time for measuring cold hardiness was longer than that for estimating dormancy depth, separate samples for determination of cold hardiness were collected every 15–20 days, resulting in 10 collections in total (Figure 1B).

Five individuals of each species were selected for measurement of dormancy depth. One twig (bearing winter buds) approximately 20–30 cm long was collected from each individual on each sampling date. Samples were transported to the laboratory immediately after collection. We placed five twigs from each species in a glass bottle (0.25 L) containing tap water and transferred the bottles to growth chambers (Ningbo Jiangnan Instrument Factory). The growth chambers were set to day/night temperatures of 25°C/15°C and a photoperiod of 14 h/10 h (light/dark) with a flux density of 100 μmol⋅m^–2⋅s^–1. A tray filled with water was placed inside the growth chamber to ensure that the relative humidity was maintained at approximately 70%. To prevent blockage of the xylem in the twigs, the base of the stems was pruned, and the water in the bottles was replaced weekly.

The development of each twig was monitored every 2 days. The first leaf date (FLD) and first flowering date (FFD) were recorded. The FLD was defined as the date when at least one leaf in a twig was fully unfolded, corresponding to BBCH 11. The FFD was defined as the date when the first flower on a twig had fully opened, corresponding to BBCH 60 (Meier, 2001). We only investigated the first spring event of each species. For two species that flowered earlier than leaf unfolding (J. nudiflorum and F. suspensa), the FFD was recorded. For the other three species (F. chinensis, S. babylonica, and S. oblata), the FLD was recorded.

The growing degree day (GDD) requirement of FFD or FLD in the growth chamber was used to represent the dormancy depth (Dantec et al., 2014) or to assess the relationship between chilling and GDD requirement (Laube et al., 2014). The GDD requirement was calculated as the accumulation of temperature greater than a certain threshold. For twigs sampled before 1 January, 2019, only the GDD in the growth chamber was calculated because the ambient temperature was below 5°C, whereas for twigs sampled after this date, the GDD already accumulated outdoors was also included. The equation used to calculate GDD was:

where GDD_R represents the GDD requirement of FLD and FFD; t₁ is the date when GDD starts to accumulate (1 January, 2019); t_s is the sampling date; t_f is the date of onset of FLD or FFD averaged from five twigs; T_o(t) is the daily mean temperature outdoors at date t; T_g(t) is the daily mean temperature inside the growth chamber, which was 20.83°C based on the day/night temperature and photoperiod set in the growth chamber; and T_b is the threshold temperature for GDD, which was set at 5°C based on previous studies (Kramer, 1994; Peaucelle et al., 2019; Wang et al., 2020a). All twigs were observed for 55 days after the sampling date (corresponding to 871°C⋅days, i.e., the average GDD from January to May in Beijing).

To describe the seasonal variation in dormancy depth, the variation in GDD_R (average of five individuals per species) with sampling date was plotted. Based on the shape of the curve, an exponential function was used to fit the changes in GDD_R (after the date of deepest dormancy):

where a, b, and c are parameters, and date is the number of days since 1 October, 2018. The parameters were estimated using the least-squares method.

Previous studies showed that the GDD_R of spring events is associated with chilling accumulation (CA) during winter (Fu et al., 2014; Wang et al., 2020b). Therefore, CA can be used as a bioclimatic index of dormancy depth. CA was calculated with the equation:

where CA is the chilling accumulation. t₀ is the date when CA started to accumulate (1 November, 2018), t_s is the sampling date. T_o(t) is the daily mean temperature at date t.

An exponential function between CA and GDD_R was established for each species:

where a₁, b₁, and c₁ are parameters estimated using the least-squares method.

For each species, 30 twigs bearing dormant buds from 5 individuals were collected. The collected samples were randomly divided into six groups and preserved in plastic bags. The collected samples were first washed with deionized water. Subsequently, five groups were placed in programmable freezers each maintained at a different target temperature, and one group was placed in a 5°C refrigerator (living control). The target temperatures of the five programmable freezers were adjusted according to the changes in cold hardiness (Table 1; Lenz et al., 2013; Vitasse et al., 2014; Vitra et al., 2017). All programmable freezers were cooled from 5°C at a rate of 6°C h^–1 until the target temperature was attained. The samples were kept at the target temperature for 4 h and then returned to 5°C at a rate of 6°C h^–1. Finally, the dormant buds were removed from the twigs and placed in test tubes with 25 mL deionized water. The test tubes were placed in the refrigerator at 5°C for 24 h. The conductivity (S₁, i.e., electrolyte leakage after treatment) of the water in the test tubes was measured with a conductivity meter, and then the electrical conductivity (S) was measured again after autoclaving the samples at 120°C for 30 min. The relative conductivity (E) was calculated by combining the conductivity of deionized water (S₀) with the following equation:

The relationship between E and the target temperature (T) was analyzed by logistic regression:

where E(T) is the relative conductivity at temperature T, and i, j, and K are model parameters. The temperature at a relative conductivity of 50%, i.e., the semi-lethal temperature (LT₅₀), represents cold hardiness (Vitasse et al., 2014).

To describe the seasonal variation in cold hardiness, the cold hardiness (LT₅₀) of each species (averaged from five individuals) with sampling dates was plotted. Because the cold hardiness initially increased and thereafter decreased during the study period, a double logistic function was used to fit the curves:

where x is the sampling date, y is the cold hardiness at date x, and k₀, k₁, and k₂ are parameters estimated using a least-squares method.

To simulate the seasonal variation in cold hardiness, the cold hardiness model proposed by Leinonen (1996) was modified and a normalized cold hardiness index (CHI) was constructed. CHI = 0 represents the lowest cold hardiness and CHI = 1 represents the strongest cold hardiness. We use [R(t)] to represent the CHI on day t. [R(t)] on day t depends on the CHI on the preceding day [R(t - 1)], the potential cold hardiness on day t [R_c(t)] and on day t - 1 [R_c(t - 1)]:

where w is the weight of potential cold hardiness R_c(t) on actual cold hardiness. w was set to 5, because the R(t) corresponded better with the measured cold hardiness when w = 5 compared with other values (2, 3, 4, and 6, 7, 8). The calculation of CHI started on 1 October. Because the outdoor temperature was still approximately 18°C on 1 October, it was assumed that cold acclimation had not started yet. The R(t) on 1 October was set to its minimum value (zero), and then R(t) on the next day was calculated with a daily scale.

The potential cold hardiness [R_c(t)] was determined by the effect of temperature [R_T(t)], the effect of photoperiod [R_p(t)], and hardening capacity (C_R). The process of cold acclimation before the winter solstice is driven by temperature and photoperiod (Kalberer et al., 2007), whereas deacclimation is mainly driven by temperature (Leinonen et al., 1997). Therefore, R_c(t) was calculated as:

where R_T(t) is the effect of temperature on cold hardiness, and R_P(t) is the effect of photoperiod on cold hardiness. The equations for calculation of R_T(t) and R_P(t) are:

where T_up is the maximum temperature limit, which was set to 10°C, T_low is the minimum temperature limit, which was set to −16°C (Leinonen, 1996), and T(t) is the daily mean temperature on date t.

where P_up is the maximum daylength limit, which was set to 14 h, P_low is the minimum daylength limit, which was set to 8 h, and P(t) is the photoperiod on date t, which is calculated based on the latitude and the day of the year (Forsythe et al., 1995).

C_R is assumed to increase linearly from 1 October to the winter solstice and thereafter decrease linearly from 1 to 0 from the winter solstice to the onset of spring events. Thus, the equation for calculation of C_R is:

Regarding the initial value of C_R on 1 October, a set of values (0, 0.1, 0.2, 0.3, and 0.4) were tested, and the correlation coefficient between CHI and measured cold hardiness was largest at C_R = 0.2. Therefore, the initial value of C_R was set to 0.2.

After calculation of the daily CHI for each species, a linear regression model was calibrated to simulate cold hardiness (CH):

where m is the slope and n is the intercept, estimated with the least squares method. CH is the cold hardiness measured in the experiment. CHI is the daily cold hardiness index [R(t)].

Using the daily mean temperature data during 1952–2021, CA from 1 November to each day until the end of April was calculated for each species based on Equation 3. The course of dormancy depth (GDD_R) over the last 70 years was estimated through Equation 4, and the interannual spring phenology was simulated as the date when GDD_R was attained. We also obtained the long-term FLD and FFD data of five species in Beijing from the China Phenological Observation Network (CPON) during 1963–2018. The historical phenological data were collected at the Summer Palace, about 12 km from the sampling site of the experiment. We compared the simulated spring phenology with the long-term observed data to verify the robustness of phenological simulation.

The daily CHI from 1 October to 30 April was calculated for each year and Equation 13 was used to transform CHI into cold hardiness of each species. Finally, the interannual changes in dormancy depth (on 1 March, because after this date the daily mean temperature is always >5°C) and maximum cold hardiness were analyzed.

Based on the simulated cold hardiness, we identified historical frost events (1952–2021) from the previous winter (1 November) to spring as the dates when daily minimum temperatures equal or below the estimated cold hardiness.

Furthermore, we adopted a widely used false spring index (FSI) to quantify the likelihood of damage to occur after budburst (Marino et al., 2011; Chamberlain et al., 2019).

where LFD is the last frost date in spring (the last day when daily minimum temperature <-2.2°C). BBD represents the budburst date, i.e., FLD or FFD of each species. Negative values indicate no-risk situations, whereas a positive FSI indicates cold damage occurred to young flowers or leaves.

For most species, budburst did not occur in the earliest-collected samples with sufficient GDD (871°C⋅days, equivalent to the outdoor GDD from January to May in Beijing), indicating that the twigs were in the maximum dormancy. For S. babylonica, the samples collected before 21 October underwent budburst, indicating that dormancy induction was not complete. However, budburst failed to occur in the samples collected between 21 October and 8 November, indicating that the buds entered deep dormancy during this period. Overall, maximum dormancy depth generally occurred between 1 October and 3 December. Thereafter, all species showed a clear exponential decrease in dormancy depth (Figure 2). An exponential function was used to fit the change in GDD_R with sampling dates. This function simulated the variation in dormancy depth with R² ranging from 0.59 (F. chinensis) to 0.94 (F. suspensa). Therefore, the exponential function accurately reflected the change in dormancy depth over time.

The CA was correlated negatively and exponentially with GDD_R (Figure 3). For all species, the GDD_R decreased rapidly with increase in CA. However, when CA attained a certain threshold, the effect of CA on GDD_R decreased. The R² of the exponential function between GDD_R and CA ranged from 0.61 to 0.95, indicating that CA was a good indicator of dormancy depth.

The cold hardiness of all species initially increased from autumn to mid-winter and thereafter decreased from mid-winter to spring (Figure 4). The maximum cold hardiness was observed between 29 November (S. babylonica) and 22 January (J. nudiflorum and F. chinensis). F. chinensis had the strongest cold hardiness of −11.45°C on 1 October, which increased to −27.68°C on 23 December. Jasminum nudiflorum had the weakest cold hardiness of −5.16°C on 1 October, which only increased to −18.34°C on 30 December. The temporal variation in cold hardiness for the five species was expressed by a double logistic function with R² ranging from 0.83 (J. nudiflorum and F. chinensis) to 0.91 (S. oblata). The maximum cold hardiness (Figure 4) was observed later than the maximum dormancy depth (Figure 2) for all species.

The cold hardiness measured by the electrical conductivity method was correlated significantly and negatively with CHI for all species (Figure 5). The Pearson correlation coefficients (r) ranged from 0.82 (S. babylonica) to 0.93 (F. suspensa). In the processes of both cold acclimation and deacclimation, CHI corresponded well with cold hardiness, indicating that CHI may be used as an indicator of actual cold hardiness.

Over the past 70 years, the CA from the preceding 1 November to the mean date of spring events decreased significantly with a linear trend of 3.0 days decade^–1 (Figure 6A). Accordingly, the GDD_R of spring events (on 1 March) significantly increased for all species, and ranged from 0.06°C⋅days decade^–1 (J. nudiflorum) to 2.17°C⋅days decade^–1 (S. babylonica) (Figures 6B–D). Therefore, dormancy depth in early spring increased in strength for all species in the context of climate warming.

The maximum CHI showed a trend of significant decrease over time for all species (Figure 7A). The changes in cold hardiness were simulated by a linear regression function between CHI and the measured cold hardiness for each species (Figure 7B). All species showed a significant decrease in cold hardiness with trends ranging from 0.22°C decade^–1 (J. nudiflorum) to 0.39°C decade^–1 (S. babylonica). Therefore, the reduced duration and magnitude of low winter temperatures have weakened the cold hardiness of all species to varying degrees.

To identify the historical frost events before or after budburst, we first simulated the annual FLD or FFD using the phenological model based on CA-GDD_R functions. The results showed that the simulated budburst dates of five species advanced by 2.12 days decade^–1 to 2.44 days decade^–1 over the past 70 years (Supplementary Figure 1). Then, we validated the simulated spring phenology with the long-term observation data at a nearby site in Beijing. The results showed that the simulated spring phenology correlated positively and significantly with the observations (Figure 8), suggesting that the phenological models in this study well reflect the interannual changes in the spring phenology of all experimental species. However, the linear regression lines between the simulated and observed dates deviated from the 1:1 line for three species (J. nudiflorum, F. suspensa and F. chinensis), possibly because of the phenological differentiation between populations and the difference in microclimate between two locations.

When comparing daily minimum temperatures and cold hardiness simulated by CHI, we found that the frost events before budburst rarely occurred in the past 70 years. Frost events only occurred several times in 4 years (1966, 1969, 1971, and 1979/1981) for J. nudiflorum and S. babylonica, but occurred only in one year (1966) for the other 3 species (Supplementary Table 1). Thus, the species investigated here have adapted to the seasonal temperature variation of the study area to a great extent.

We also assessed the post-budburst frost events in the past seven decades. The results showed that frost damage (positive FSI) only occurred in three species. Frost damage on young leaves or flowers occurred in 2, 1, and 1 years (all before 1970) on J. nudiflorum, S. oblata and S. babylonica, respectively (Figure 9). Furthermore, there were no significant trends in FSI due to the synchronous advance of the last frost date and budburst date in response to climate warming (Supplementary Figure 1).

After endodormancy induced by the shortened daylength and lower temperatures in autumn, budburst in temperate trees and shrubs may not occur even in a growth-promoting environment (Campoy et al., 2011; Basler and Körner, 2012; Singh et al., 2017). In the present study, the maximum dormancy generally occurred between 1 October and 3 December for all five study species. This result is in accordance with previous findings. For example, the maximum dormancy depth of peach in Georgia, USA, was observed between late October and early December (Weinberger, 1950). Similarly, Campoy et al. (2011) reported that the maximum dormancy depth for apple usually occurred 1–2 months after the autumn equinox, which was consistent with the present results.

The release of endodormancy requires a certain amount of CA. In the current study, the dormancy level gradually decreased with increase in CA after plants attained the maximum dormancy. This result confirmed the effects of CA observed in previous studies. For example, the GDD requirements for budburst of seven temperate species gradually decrease with increase in chilling duration in Ontario, Canada (Man et al., 2017). The GDD requirement for budburst of European beech (Fagus sylvatica) is negatively correlated with duration of chilling (Vitasse and Basler, 2013). The effects of CA on dormancy depth are also evident in observational data for 30 perennial species in Europe (Wang et al., 2020b). Although chilling plays a similar role in the dormancy release of temperate species, the different sensitivity of GDD requirements to CA among species may lead to varying phenological responses to climate warming (Hu et al., 2022).

The cold hardiness of the five temperate species increased from the beginning of October, with the maximum cold hardiness observed at the end of December. This pattern has been reported in previous studies. For example, experiments on two Japanese cultivars of tea plant (‘Yabukita’ and ‘Yutakamidori’) in Kagoshima showed that cold hardiness gradually increases from the beginning of December to a maximum value in mid-February of the following year and thereafter steadily decreases (Kimura et al., 2021). The cold hardiness of twigs of 11 woody plant species at different altitudes in central France increases from early November to late winter (January–February) and then gradually decreases until budburst (Charrier et al., 2013). The cold hardiness of four woody plant species in Switzerland gradually decreases from the end of January to the beginning of April (Vitra et al., 2017). During the transition from autumn to winter, plant tissues start to harden owing to the decrease in temperature and the shortening of daylength (Beck et al., 2004), whereas cold deacclimation is mainly driven by progressively warmer temperatures in spring. Maximum cold hardiness differed among the five species, ranging between −18°C and −28°C, in the current investigation. The maximum cold hardiness of fruit-oriented walnut can attain −28°C, whereas that of wood-oriented walnut is −35°C (Charrier et al., 2011). Nordic Betula species and other cold-adapted species can tolerate temperatures lower than −70°C at full dormancy (Körner, 2012).

The CA and CHI provided sufficient indications of seasonal variations in dormancy and cold hardiness. The present simulations indicated that dormancy depth is increasing owing to decreased chilling. Similar results have been observed in previous studies. For instance, Wang et al. (2022) reported that insufficient chilling under climate warming has led to an increase in dormancy depth, which has approximately halved the rate of advance of spring phenology caused by increased forcing.

Temperate species with long life cycles have evolved the ability to cope with low seasonal temperatures within their distribution range. This explains why we detected very few frost events before budburst in the past 70 years (Supplementary Table 1). Since lower temperatures are more effective for cold acclimation, the present study simulated a decrease in cold hardiness from the previous autumn to spring under climate warming. The decreased cold hardiness may cause a higher risk of frost damage. However, the occurrence of extremely low-temperature events decreased significantly under warming (Griffiths et al., 2005). Therefore, combined with decreased cold hardiness and frequency of low-temperature events, frost risk prior to budburst would not change significantly. The results of this study confirmed this assumption, as no frost events were simulated to occur in recent years.

The budburst date of temperate plants is now earlier owing to the influence of climate warming (Menzel et al., 2006; Ge et al., 2015). As a result, plants may be at increased risk of frost damage after budburst (Chamberlain et al., 2019). However, the earlier last frost date in a warmer climate tends to reduce the spring frost risk (Park et al., 2021). Thus, considering both the changes in spring phenology and the frequency of frost events, the post-budburst frost risk may vary among regions. This study found that the frost risk to young leaves and flowers did not change significantly. In Switzerland, the risk of spring frost exposure and subsequent potential damage has increased for four tree species at the vast majority of stations located at elevations higher than 800 m while remaining unchanged at lower elevations (Vitasse et al., 2018). A recent study showed that the frequency of frosts in late spring has increased to a greater extent in Europe and East Asia than in North America from 1959 to 2017 (Zohner et al., 2020).