Quantifying the trade-off between spring phenology and lethal frost risk: a meta-analysis

Domain: chemistry
Published: 2026-03-05
Editorial outcome: major_revision

Abstract
The timing of spring phenology represents a critical trade-off between growing season extension and frost risk avoidance, the balance of which governs plant fitness and shapes plant species’ distributions. However, few quantitative assessments exist, indicating that this trade-off is largely untested empirically. Here, we present a global meta-analysis encompassing 193 plant species across 126 study sites, finding consistently high freezing resistance (to -12°C) with remarkably low lethal frost risk (safety margin of 17°C) during spring emergence. Under projected climate warming scenarios, predictions indicate that advances in spring phenology do not affect the lethal frost risk under light and moderate warming scenarios. This is likely due to the reduced temperature sensitivity, consistent with trade-off predictions. By contrast, in the high-warming scenario, the freezing safety margin is predicted to expand (by ~1.7°C), resulting in a lower lethal frost risk. This expansion may compensate for the diminished freezing resistance predicted under extreme warming. Our findings challenge the conventional view that simple environmental factors moderate spring phenology, and highlight the necessity of incorporating biotic factors (e.g., freezing resistance) and biotic processes (e.g., early growth versus frost risk trade-offs) into next-generation phenological models to enhance prediction accuracy under climate change.

Introduction
Plant phenology, or the timing of recurrent biological events during the year, plays a pivotal role in shaping species distributions, regulating water and nutrient cycles1,2,3. In recent decades, advanced spring phenology with warmer springs is widely reported4,5. However, it remains unclear whether this phenological advance, accompanied by warmer temperatures, alters the risk of lethal frost (LFR).

Theoretically, spring phenology represents an adaptive compromise between maximizing the growing season and avoiding lethal frost damage, which reflecting an evolutionary adaptation shaped by long-term selection pressures6,7,8,9,10. The timing of spring phenology under climate warming arises from adaptive trade-offs between two competing phenological strategies. Emerging earlier provides a competitive advantage by extending the photosynthetic period, maximizing light-use efficiency and photosynthetic activity11,12, reducing herbivore/pathogen pressure13,14, and providing priority access to ecological niches15, all of which may enhance plant fitness and reproductive success16. On the other hand, potentially lethal spring frost events favor delayed emergence, which protect vulnerable tissues17 and enabling greater resource allocation to growth and reproduction instead of frost tolerance18,19 (Fig. 1a, Supplementary Table 1). Ultimately, natural selection favors a balance between these two opposing pressures, the benefits of early emergence versus the risk of lethal frost events, which allows individual species to optimize their spring phenology based on local environmental conditions. This evolutionary balance, known as the evolutionarily stable strategy (ESS), provides a theoretical framework for predicting plant phenological responses to a changing climate.

Specifically, this balance manifests as a buffering mechanism: safeguarding against lethal frost damage during spring phenological stages while simultaneously maximizing individual fitness. Such trade-offs broadly shape plant evolutionary strategies and explain spring phenological variation across species and regions18,20. Therefore, the optimal phenological strategy for plants under climate change is to maintain this evolutionary balance. Lenz et al.9 propose a metric termed the ‘safety margin’, defined as the difference between the daily minimum temperature and freezing resistance during spring emergence, which operationalizes and quantifies this evolutionary balance (Fig. 1a). The key advantage of this dynamic metric lies in its incorporation of biotic factors (e.g., plant freezing resistance), rather than relying solely on abiotic factors21,22,23,24,25. In fact, plants must cope with stochastic freeze events throughout the year6. The capacity to withstand such events, termed freezing resistance, peaks during deep physiological dormancy in winter and declines sharply with spring leaf emergence7. Consequently, when both biotic and abiotic factors are considered simultaneously, a larger safety margin indicates greater tolerance of freezing (and lower frost risk), while narrower margins correspond to greater vulnerability (Fig. 1a).

Climate-change-induced shifts in spring phenology may fundamentally alter this risk-reward trade-offs in natural plant populations22,26, leading to three potential distinct risk scenarios (Fig. 1b). In the high-risk scenario, earlier springs confer transient competitive advantages, but a single extreme climatic event beyond a critical thermal threshold could result in local extinction7 (Fig. 1b). The low-risk scenario allows for near-complete frost avoidance, but comes at the cost of a competitive disadvantage that may lead to replacement by other aggressive species (Fig. 1b). The stable-risk scenario, which best aligns with the EES, illustrates the fundamental trade-off between advancing growth and avoiding lethal frost damage (Fig. 1b). However, empirical data supporting the existence of a spring phenology risk-reward trade-offs under warming conditions remain scarce, rendering this theoretical framework highly speculative.

In this study, we investigate (i) plant spring freezing resistance and associated safety margin, as well as their environmental drivers; (ii) how warming-induced spring phenological shifts may modify risk-reward trade-offs. To this end, we perform a global meta-analysis encompassing 193 plant species belonging to 336 orders, 54 families, and 122 genera across 126 sites (Fig. 2a, c). Data are collected from 88 peer-reviewed articles (Supplementary Data 1) published before September 1st, 2024 (Supplementary Fig. 1). The meta-analysis quantifies lethal frost risk during spring emergence and identifies key environmental drivers predicting freezing resistance and the safety margin. We also implement eight process-based phenological models to evaluate the trade-off between early emergence and frost risk under climate warming, as projected by three Shared Socioeconomic Pathways (SSPs). We show that biotic factors and the trade-off mechanisms play a critical role in regulating spring phenology.

Results
Low lethal frost risk during the spring phenology
Overall, study species exhibited strong spring freezing resistance (-12.42°C ± 5.19°C, mean value ± SD) relative to the daily minimum temperature during spring phenology (Fig. 2d), resulted in wide safety margin (17.14°C ± 5.58°C) which ensured an low lethal frost risk (Fig. 2e). The observed patterns of both freezing resistance (K = 0.059, p = 0.144) and safety margin (K = 0.048, p = 0.409) did not differ significantly from the null expectation under a Brownian motion model, indicating that the influence of phylogenetic relatedness did not differ from what would be expected by chance (Fig. 2d, e). Furthermore, convergent cross mapping (CCM) analysis identified a strong mutual causality between freezing resistance and safety margin, providing evidence for a fundamental adaptive trade-off within plants, balancing cold hardiness against growth opportunity (Fig. 2f, Supplementary Fig. 2).

Across studies, after accounting for potential sample effects (e.g., smaller n in shrubland, n = 17), plant species within the same ecosystem showed similar levels of freezing resistance and safety margin (Fig. 3). Neither freezing resistance nor the safety margin varied significantly among different life forms (i.e., trees, herbs, and shrubs) or across life cycles (i.e., perennial and annual species). However, early-spring species showed higher freezing resistance than medium-spring species, while medium-spring species exhibited a lower safety margin than both early- and late-spring species (Fig. 3). Surprisingly, freezing resistance estimates were sensitive to methodological variation. Values from differential thermal analysis (DTA) were systematically overestimated, leading to a consequent overestimation of the safety margin (Fig. 3). However, after excluding these observations, we found our results remain unchanged (Supplementary Fig. 3, 4).

Environmental and ecological drivers
For freezing resistance, random forest modeling identified the three most important predictors as annual maximum temperature (Tmax) (24.34%, p < 0.01), diurnal temperature range (DTR) (21.62%, p < 0.01), and mean annual temperature (Tmean) (17.62%, p < 0.01) (Fig. 4a). It is worth noting, that the relationship between mean freezing resistance and DTR (slope = 0.16, p = 0.33) was not significant in the linear mixed model analysis (Supplementary Fig. 5e). Absolute latitude (slope = -0.06, p < 0.01) and frost day frequency (FRSE) (slope = -0.02, p = 0.01) exhibit strong and moderate negative correlations, respectively. Daily minimum temperature during phenology (Phe-Tmin) (slope = 0.37, p < 0.01), Tmean (slope = 0.34, p < 0.01), Tmax (slope = 0.35, p < 0.01), evapotranspiration (Eva) (slope = 0.01, p < 0.01), solar radiation (slope = 1.18, p < 0.01) and sunshine duration (Ssh) (slope = 0.01, p < 0.01) show strong positive correlations, while Tmin (slope = 0.14, p = 0.01) showed a moderate positive correlation (Fig. 4b, c, d, e, f, j). For forest-inhabiting species, the most significant predictors of freezing resistance were Tmax, DTR, and Tmean, whereas for grassland species, precipitation and elevation were most influential (Supplementary Table 4).Freezing resistance in early spring species were more dependent on solar radiation and water availability, while mid- and late-spring species relied more on Tmax (Supplementary Table 4).

For safety margin, freezing resistance (FR) (60.42%, p < 0.01) was the strongest predictor, followed by Phe-Tmin (27.65%, p < 0.01) and FRS (17.90%, p < 0.01) (Fig. 5a). Linear mixed models revealed strong positive correlations with Phe-Tmin (slope = 0.63, p < 0.01), Tmean (slope = 0.24, p < 0.01) and Tmin (slope = 0.21, p < 0.01), and strong negative correlations with minimum temperature variance (Tmin-var) (slope = -0.68, p < 0.01) and FRS (slope = -0.03, p < 0.01) showed the strong negative correlations (Fig. 5b, c, d, e, f, g). Our analysis revealed consistent safety margin across ecosystem types, life cycles, life forms, and phenological strategies (Supplementary Table 5).

Spring phenology and safety margin under future warming scenarios
Based on projected future climatic data from the Coupled Model Intercomparison Project Phase 6 (CMIP6) (Supplementary Fig. 6), we compared spring phenology predictions (2020-2090) from eight models (Supplementary Fig. 7). Predictions were generally consistent across the models (Supplementary Fig. 8, 9). Given this agreement, we generated an ensemble by averaging the outputs of all eight models. The results showed climate warming significantly advanced spring phenology, at a rate of -7.14, -6.82, and -5.71 days per degree Celsius under the low-, medium-, and high-warming scenarios, respectively (Fig. 6a). To account for these shifts when calculating the freezing risk, we extracted projected daily minimum temperature data from CMIP6 corresponding to the updated onset timings of spring phenology. Regarding future freezing resistance, numerous studies have demonstrated that warming leads to reduced chilling accumulation, which may impair freezing resistance8,19,28,29,30. However, the actual impact depends on the future climate scenarios31, data sources, and the chilling model used32. Therefore, we projected chilling accumulation for 2021-2099 using four climate scenarios and four process-based chilling models. Most models indicate that chilling accumulation was not always diminished under warming, however, half of the results show a strong decline only under the SSP5-8.5 scenario (Supplementary Fig. 10). Additionally, our in situ observations show no significant change in the freezing resistance of the 15 perennial herbs under moderate warming (Supplementary Fig. 11, 12, see methods). Despite the inherent methodological limitations of such manipulative studies, our results suggest that plant freezing resistance may remain stable at least two warming scenarios.

Subsequently, we calculated the future safety margin using the same methods described previously (Supplementary Fig. 14, see methods). Interestingly, under both the low- and medium-warming scenarios, the safety margin predictions remained stable (p = 0.51 and 0.11, respectively; Fig. 6b), which aligns with the stable-risk scenario outlined in Fig. 1b. In contrast, under the high-warming scenario, the safety margin increased significantly, even though spring phenology exhibited the lowest temperature sensitivity (Fig. 6a). To address potential methodological biases from differential thermal analysis (DTA), we repeated all analyses after excluding DTA-derived observations. The similarity of the results reinforced confidence in our findings (Supplementary Fig. 14).

Methods
Data compilation
In this study, we present a global meta-analysis of plant freezing resistance and its relationship to lethal frost risk in the context of spring phenology. We searched for journal articles on the ISI Web of Science and the China National Knowledge Infrastructure using the following key word combinations: 'freezing resistance*frost resistance*cold resistance' OR 'freezing tolerance*frost tolerance*cold tolerance' OR 'freezing hardiness*frost hardiness*cold hardiness' OR 'hardening*de-hardening' OR 'cold acclimation*cold de-acclimation*cold re-acclimation' AND 'leaf*foliage*bud*'. We selected studies that satisfied the following criteria: (i) data were collected at the species level under natural conditions, excluding garden or glasshouse studies to avoid human management; (ii) freezing tests were conducted during spring phenology on new leaves/buds using standardized methods (Supplementary Fig. 1), excluding studies based solely on visual scores; and (iii) the study site location was clearly documented. For experimental studies, our meta-analysis included only data from control groups. Means, standard errors and/or standard deviations, and sample sizes were extracted from the study text, tables, and/or figures using the GetData (V. 2.26) image tool. To ensure the independence of data from individual studies, we averaged freezing resistance measurements for each species within a study site. In the final dataset (Supplementary Data 1), the selected studies spanned 110.53° of latitude (from 65.45°N to 54.06°S), 1°-170.23° of absolute longitude (Supplementary Table 2), and almost all terrestrial biomes (Fig. 2b).

We quantified the lethal frost risk using the safety margin (SFM)9, which is defined as the difference between plant freezing resistance and the absolute daily minimum temperature during spring phenology (Fig. 1a):
SFM = Phe-Tmin - LT50 (1)
where Phe-Tmin is the daily minimum temperature during spring emergency, LT50 means the freezing resistance (i.e., the lethal temperature of which 50% of samples died) (Supplementary Table 3). In context of climatic instability and unpredictable extreme cold events, a larger safety margin translates into a lower lethal frost risk. Based on this negative relationship, we quantified the species-specific lethal frost risk in spring (Fig. 1a).

From each study, we compiled ancillary site information spanning three categories: site spatial distribution, temperature regime, and water availability. To characterize the sites spatially, we recorded their geographic coordinates (latitude, longitude), elevation, and key solar parameters (solar radiation and sunshine duration). These data are critical because photoperiod cues can protect plants from a potentially fatal consequences, initiating growth at the wrong time by relying solely on ambient temperature49. We compiled key temperature parameters, including the annual means, extremes, variances, diurnal range, frost frequency, and daily minimum temperature during spring phenology. These variables are established drivers of cold acclimation/de-acclimation processes and the associated biochemical activities that determine freezing resistance37,38. We compiled site-level moisture data, including drought severity, evapotranspiration, and precipitation, due to their established role in governing freezing resistance mechanisms53 (Supplementary Table 3). Where possible, we extracted site-level data from the individual studies, missing data were obtained from the WorldClim Database and/or the Climatic Research Unit (CRU TS Version 4.06, provides daily climatic data with a spatial resolution of 1° from 1981 to 2021). Additionally, we recorded the ecosystem type (forest, grassland, or shrubland)54, life cycle (perennial or annual plant species), and life form (tree, herb, or shrub), these data were gleaned from the papers themselves, with missing data collected from online databases such as the USDA Plants Database or the China Flora/Fauna Handbook (Supplementary Table 2). Study species were classified into early-, mid-, and late-spring groups based on their emergence dates, reflecting distinct phenological strategies. For Southern Hemisphere observations, dates were standardized by subtracting 180 days, following the typical green-up period in the latter half of the year.

A central objective of this study is to project the future risk of lethal frost. This is achieved by integrating projected shifts in spring phenology and the potential changes in freezing resistance under warming scenarios. Such an integrated approach is critical, as the direct impact of climate warming on freezing resistance itself remains poorly understood. To evaluate the freezing resistance during spring phenology under warming conditions, we selected 15 perennial herbaceous species at the Qinghai Haibei National Field Research Station of Alpine Grassland Ecosystems (37°29′-37°45′N, 101°12′-101°23′E, 3200m a.s.l., see our previous studies for site details55,56,57). The warming treatment increased air temperature by an average of approximately 1.8℃ over a ten-year period. Plant cutting samples (leaves) were cooled to each target temperature (using commercial freezers) and subsequently thawed at a rate of approximately 6℃ h-1 until room temperature was reached (Supplementary Fig. 11a). Initial assessment of freezing damage was conducted through visual inspection, such damage typically manifested as tissue discoloration due to polyphenol oxidation, or as a characteristic odor resulting from protoplast decompartmentalization and autolysis9. Following visual inspection, cutting samples were transferred to small beakers containing distilled water for electrical conductivity (Ec) measurement. These Ec assays were conducted hourly until values stabilized, a process typically requiring over ten hours (the slow leakage of electrolyte). Concurrently, fresh samples were boiled at 100℃ for 30 minutes and then cooled to room temperature to determine the maximum electrolyte conductance (Ecmax). In this study, freezing resistance was defined as the lethal temperature (LT50, ℃) at which electrolyte leakage reached 50% of ECmax 58,59,60, and this LT50 value serves as a common proxy for the temperature at which 50% of the sample cells died9 (Supplementary Fig. 11b).

Phenological models and future climatic data
We employed eight phenological models, four one-stage and four two-stage, to project future spring phenology. One-stage model assume that spring phenological events are triggered when plants have experienced sufficient warming (forcing temperature) over a fixed period (Eq. 2). Note the dormancy phase is not accounted for, as the requirements for dormancy release (i.e., sufficient chilling) are assumed to have been met before the starting date61. The state of forcing (Sf) is modeled as:
Sf(D) = ∑R(T) (2)
where R determines the rate of change in Sf, t is the first day of the forcing period, and D is the day of the year. Spring phenological events are triggered when Sf reaches a forcing threshold F*:
Sf(D) ≥ F* (3)
In this study, we employed four forcing models to estimate the growing degree days (GDDs) needed for the initiation of spring growth32.

The first model F1, assumes that the rate of forcing is linearly correlated with the temperature above a particular threshold. This widely accepted model adopts a temperature threshold of 0℃ in calculating the accumulated GDDs (Eq. 4)62:
T0 = max(T, 0) (4)
where T0 is the rate of forcing for model F1 and T is the daily mean temperature (℃).

The second model F2, assumes that the rate of forcing represents a logistic function of temperature. This model allows a larger forcing rate in relatively warm spells. Here, we used a logistic function to describe this linear function as proposed by Hänninen63:
Tlf = 28.4 / (1+e-0.185(T-18.5)) (5)
where Tlf is the rate of forcing for model F2 and T is the daily mean temperature (℃).

According to Piao et al.62, the onset of spring phenological events is more strongly influenced by daytime than nighttime temperatures. Thus, we estimated a third model F3, based on the maximum temperature with a threshold of 5℃:
Tmax_R = max(Tmax - 5, 0) (6)
where Tmax_R is the rate of forcing for model F3, and Tmax is the daily maximum temperature (℃).

The fourth model F4, considers both daytime and nighttime temperatures, as Fu et al.64 demonstrated that the impact of daytime temperatures on spring phenology was approximately threefold higher than that of nighttime temperatures. Thus, we utilized a temperature-mixed model with two parameters (0.75 and 0.25) to weight the impact of daytime versus nighttime temperatures (℃) on GDD accumulation:
Tmix = 0.75 × max(Tmax - 5, 0) + 0.25 × max(Tmin - 5, 0) (7)
where Tmix is the rate of forcing for model F4, and Tmax and Tmin are the daily maximum and minimum temperatures (℃), respectively.

The two-stage phenological models are based on the premise that spring phenology is initiated only after sufficient chilling accumulation. We applied the commonly used Parallel model, a chilling-based framework, to simulate future spring phenology32. The state of chilling (Sc) is modeled as:
Sc(D) = ∑Rc(T) (8)
where Rc determines the rate of change in Sc, t is the first day of the chilling period, and D is the day of the year. The triggering of spring phenological events occurs when both Sf (Eq. 2) and Sc (Eq. 8) reach their specific thresholds. We adopted the commonly used model for calculating Sf (Eq. 4), and employed four chilling models to quantify the chilling accumulation required for spring phenology32:
if Sf(D) ≥ F* then leaf induction is completed
if Sc(D) ≥ C* then dormancy is completed (9)

The first model C1, also known as the Utah Model, incorporates differential temperature weighting (Eq. 10). It was originally proposed for quantifying the chilling requirements of peach trees65:
CU1 = 0 (T≤1.4), 0.5 (1.4<T≤2.4), 1 (2.4<T≤9.1), 0.5 (9.1<T≤12.4), 0 (12.4<T≤15.9), -0.5 (15.9<T≤18), -1 (T>18) (10)

After removing the negative contributions of warm temperatures to the accumulated chilling, we obtained the modified Utah Model66, designated as model CU2:
CU2 = 0 (T≤1.4), 0.5 (1.4<T≤2.4), 1 (2.4<T≤9.1), 0.5 (9.1<T≤12.4), 0 (T>12.4) (11)
where CU1 and CU2 is the rate of chilling for model C1 and C2, respectively. T is the daily mean temperature (℃).

Note that models C1 and C2 operate on the premise that no relevant plant physiological activity occurs below freezing. Therefore, sub-zero temperatures are omitted from their chilling calculations. However, Zhang et al.67 fitted the triangular model to observational data from 24 plant species and found that a mean optimal chilling temperature of 0.2℃ and an upper chilling threshold of 6.9℃ were most effective. Consequently, we adopted a triangular function with these parameters (0.2℃ and 6.9℃) as our model C3:
CU3 = (T+6.5)/(6.9-0.2) (-6.5<T≤0.2), 0 (T≤-6.5 or T≥6.9), (6.9-T)/(6.9-0.2) (0.2<T≤6.9) (12)

Model C4 also adopts a triangular form, as defined by Hänninen63 through fitting to previous experimental results for Finnish birch seedlings:
CU4 = 3.13 * ((T+4.66)/10.93)^2.10 * e^(-((T+4.66)/10.93)^3.10) (T≤-4.7 or T≥16), 1 (2.5<T≤7.4), 0 (else) (13)
where CU3 and CU4 is the rate of chilling for model C3 and C4, respectively. T is the daily mean temperature (℃). Note that both model C3 and model C4 permit chilling accumulation to occur at sub-zero temperatures.

Given the varying observation years in our dataset (up to 2020), we harmonized all projections to a 2020 basecline using the aforementioned phenological models. For future climatic data (2020-2099), we used daily near-surface mean, minimum, and maximum temperatures from a multi-model ensemble based on CMIP6. This ensemble integrated four climate models (CAMS-CSM1-0, CAS-ESM2-0, CESM2-WACCM, MRI-ESM2-0) across three Shared Socioeconomic Pathways (SSP1-2.6, SSP2-4.5, and SSP5-8.5). These three scenarios belong to low, medium and high emission scenarios, respectively, and represent forcing levels with a peak of 2.6, 4.5, and 8.5 W m-2 by the end of this century. The projected mean, maximum, and minimum temperature anomalies under these three emission scenarios range from 1.2℃-2.2℃, 2.2℃-3.2℃, and 4.6℃-5.7℃, respectively (Supplementary Fig. 6). We posit that plant frost resistance is likely to remain stable under moderate climate warming scenarios (< 4℃)31. This stability is attributed to the relatively stable chilling accumulation under such conditions (Supplementary Fig. 10), a conclusion also supported by our in situ manipulation experiment (Supplementary Fig. 12). Therefore, future safety margins under different warming scenarios can be derived as the difference between current freezing resistance and the projected daily minimum temperature during spring phenology:
Future-SFM = Future-Phe-Tmin - LT50 (14)
where Future-Phe-Tmin is the daily minimum temperature during spring phenology and LT50 is the freezing resistance (Eq. 1).

Statistical analysis
We obtained a species-level phylogeny using the V.PHYLOMAKER package16, which utilizes the largest datced phylogeny for vascular plants (GBOTB.extended.tre). Next, we assessed whether evolutionary history was a key predictor of freezing resistance and the safety margin using Blomberg's K (R package 'phytools', version 4.1.3)68. Blomberg's K is a statistical metric that quantifies the phylogenetic signal in species traits by comparing their observed variance against a Brownian motion model of evolution:
Blomberg's K = (concordance between trait and phylogenetic relationships) / (concordance expected under a Brownian motion model) (15)
Under this model, a value of K > 1 indicates that closely related species are more similar than expected under Brownian motion, reflecting a strong phylogenetic signal. Conversely, a value of K < 1 signifies that trait variation is more random than expected, making the pattern indistinguishable from randomness and demonstrating a weak or absent phylogenetic signal. However, the K value alone is insufficient and must be accompanied by a significance test to determine if it deviates from random expectation. We therefore employed a permutation test (n = 5000) to generate a null distribution of K values (Fig. 2d, e). The observed K was compared against this distribution, with significance assessed at 0.1 level, a more lenient threshold applied consistently across all subsequent analyses in this study. We applied convergent cross mapping (CCM)69, an algorithm designed to detect causal relationships, to explore the potential causality between freezing resistance and the safety margin. The method is grounded in the principle that if X causally influences Y, then the state of X can be reconstructed from the historical record of Y, because the evolutionary trajectory of Y embeds information about X70,71. We performed the analysis using the ‘rEDM’ package in R. We tested library sizes from 31 to 304 (τ = 1) to evaluate the convergence of prediction correlation (Rho). The optimal embedding dimension (E = 5) was determined via simplex projection over a range of dimensions (E = 1-20) (Supplementary Fig. 2). We assessed convergence using the twin-surrogate method with 10 bootstrapped iterations, and considered a causal relationship significant at 0.1 level.

As part of an importance ranking analysis, we used an inverse-variance weighting method to account for heteroscedasticity across study observations. This method assigned higher weights to cases where estimates had greater precision and lower variance72. The weights were computed as:
Weighted (X) = (∑k i=1 wixi) / (∑k i=1 wi) (16)
where X is an observation, k the number of cases, and w the weighting factor.
The weighting factor took both between-case (τ2) and within-case variance in X (Vx) into account and was calculated as73:
wi = 1 / (τ2+Vx(i)) (17)
Given the absence of heterogeneity in our dataset, we restricted our weighting scheme to within-case variance components and set τ2 = 0 for all subsequent analyses. To assess the relative importance of our study predictors, we used random forest models as implemented in R (R package 'randomForest', version 4.7-1.2)74. Variable importance was quantified using the percent increase in the mean squared error (%IncMSE), with higher values indicating greater contributions. Negative %IncMSE values indicate improved model predictions when excluding the predictor in question, suggesting that the predictor should be removed from the model. Many ecological processes covary, and the sheer number of possible combinations (of predictors) in our analyses made predictor ranking and the interpretation of significant interactions complicated, while also increasing the risk of model over-fitting34,75. Thus, we retained all predictor variables (without implementing stepwise backward elimination) in the following analyses, as most explanatory variables were statistically significant (p < 0.1) (Supplementary Tables 4, 5). Next, we evaluated associations between the two response variables (i.e., freezing resistance and the freezing safety margin) and seventeen predictors (i.e., five measures of site spatial distribution, eight temperature-related traits, three moisture-related traits, and one biotic trait (freezing resistance) using linear mixed models, with random intercepts among species and sites. To interpret the mixed model results, instead of significance testing with an arbitrary p-value cutoff, we implemented the gradual language of evidence recommended by Muff et al.76 (Supplementary Tables 4, 5, 6, 7), this approach allowed us to communicate model findings in a more nuanced way16.

To determine the phenological temperature sensitivity (that is, days per ℃), we used ordinary least squares linear regression41 for the entire future period (2020-2099). We then compared future temperature sensitivity to predicted trends in spring phenology and freezing safety margins under future warming scenarios77. To assess differences in spring phenology and freezing safety margins among diverse warming scenarios, we applied linear mixed model with random intercepts among species and sites. We conducted all our statistical analyses using MATLAB R2022b and R 4.4.278.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "The authors provide interesting and important new information on how frost risk affects spring phenology, especially under climate change scenarios. Such information is important to increase the understanding of plant tolerance to extreme weather events and how this affects biodiversity. However, despite the relevance of the study, I believe that the authors still need to make significant improvements to the text to increase the clarity of several aspects discussed. I also recommend an English revision. Below I provide comments that I hope will help the authors improve the text.\n# Introduction\nL1-L24: I recommend restructuring the paragraph, prioritizing a more orderly sequence of the ideas presented. The structure with which different magnitudes of processes are presented (e.g., species distribution, spring phenology, climate change, trade-off) is a bit confusing. The authors should follow a more orderly direction in describing the different processes and phenomenons and the interconnection between them. In other words, lead the paragraph from smaller-scale processes to larger-scale processes, ending by highlighting the practical implications of these processes.\nNote: correct Supplementary Table 1: 3. Less damage by heibivory…. It is herbivory.\nL33-34: ambient minimum temperatures. I recommend adopting another term since this one can lead to some misinterpretations.\n# Results\nL72-75: Exactly. It is what I mentioned about the caption of your supplementary table 2.\nFigure 1: Define ‘DOY’ of the X-axis in the caption.\nFigure 2: Since the phylogeny of the evaluated groups is very important information for the study, I recommend that you separate the cladogram from this figure, giving it its own figure (it can be supplementary).\nFigure 3: P < 0.01? It is unclear.\n# DiscussionL\n159-161: You have not found strong evidence for evolutionary conservation of frost tolerance. Be careful with that statement.\nL162-193: This paragraph is just too big and confusing. You mix different ideas here rather than present a logical order of reasoning. For instance, you start explaining the methodological limitations of DTA, but then start to mix other concepts without a clear connection to it. I recommend you split the paragraph for clarity.\nL194: Firstly, as I mentioned before, your discussion in the previous paragraph is confusing. Second, after improving the mentioned paragraph, stress how and what your discussion exactly outlines the predictors of freezing resistance. Give the reader a clear sense of the novelty with your results.\nL218: What is this biochemical foundation required for freezing resistance? Clarify it.\nL219-221: How can it explain the mentioned negative relationship? Clarify it for the reader.\nL237-238: What is the main implication of that for the performance and reliability of your results? Make it clear to the reader.\nL263: Make it clearer how and why it remains equivocal.\nL286-288: Exemplify shortly these approaches.\nL288-290: I recommend the authors to create a section “Caveats, limitations, and implications” for addressing all blindspots or limitations of this study, as well as indicating future prospects that can assist in solving these limitations.\n# Methods\nIn the section 'Phenological models and future climatic data' you mention: \"We assumed that freezing resistance will exhibit high evolutionary conservatism ...\". But it was not the case according to the results presented in your supplementary table 2. Therefore, rephrase the sentences about that.\nIn the section \"Statistical Analysis\", at the end of the first paragraph you mentioned: \"When Blomberg's K was statistically significant (p < 0.05), this indicated that species’ evolutionary history contributed to variance in freezing resistance and the freezing safety margin\". Be careful with this statement. Not necessarily K would indicate that evolutionary history explains variance in FR. Rather, it only suggests that phylogenetic relatedness might influence variance in a trait or niche above the expected under a random fashion. However, your results don’t point to that, but that FR variance related to the phylogeny of the species is not different from the expected by chance. Further, you don’t mention here which evolutionary model you used to compute K., i.e., Brownian Motion? Ornstein-Uhlenbeck (OU)? The type of evolutionary model assumed to compute K also can be very influential in the output of its p-values. You must account for that in your study while informing and justifying clearly your choice to the reader. A last thing I mention, so you reflect about it and also account for that in your study: were average values of FR computed from sufficient numbers of individual values? That is, a reliable number of repetitions to compute average per species. Highly variable traits or niches in plants’ life history must be calculated from a robust sampling effort, so it can reflect consistent average conditions of them for each species. Therefore, you must also account for this and include such information in the manuscript as well.\nFurther in the \"Statistical Analysis\" section you mention: \"To interpret the mixed model results, instead of significance testing with an arbitrary p-value cutoff, we implemented the gradual language of evidence recommended by Muff et al., this approach allowed us to communicate model findings in a more nuanced way\". Explain more about this approach since it is a vital information for your study. Include a simple sentence explaining in slightly more details.# Data availability and Code\navailabilityThe URL links provided don't work. Update the links.\n# Supporting Information\nTable 2: Your p-values in the table don’t confirm that. Rather, the non-significant p-values you have here, related to low K values, indicate that phylogenetic relatedness influence on FR and SFM doesn’t differ from the expected by chance.\n(Remarks on code availability)\nThe codes are functional. However, they only provide commands to generate figures shown in the manuscript, but there is not enough code to portray the statistical background used.\nI recommend that authors provide all the code routines used in data processing and treatment (along with the data), as well as the statistical modeling, thus ensuring the replicability of their work.", "reviewer_2": "This manuscript addresses the important trade-off between extending the growing season and avoiding frost risk, mediated by spring phenology. Furthermore, this manuscript investigated how warming-induced phenological shifts may alter risk-reward trade-offs under future climate scenarios. This work appears to provide novel findings about the frost risk of plants, but the conclusions rely on several highly questionable assumptions, such as the LT50 remaining constant across seasons, plant tissues, and future climate scenarios. These assumptions significantly undermine the credibility of the paper's conclusions. Specific comments are detailed below:\n(1) The assumption that LT50 is temporally constant across seasons contradicts empirical observations.\nIn this study, LT50 is assumed to be constant over time (depicted as the fixed blue dashed line in Figure 1). However, experimental evidence demonstrates that freezing resistance undergoes continuous seasonal variation. Induced by short photoperiods and low temperatures in autumn, freezing resistance increases progressively—a process termed cold acclimation. Freezing resistance typically reaches its peak (i.e., lowest LT50) in December or January and subsequently declines, reaching its minimum level (i.e., highest LT50) during the active growing season. This seasonal fluctuation in LT50 often exceeds 10 to 20°C. Applying the same LT50 value to assess freezing resistance at different spring phenology timing is therefore unsound and directly renders the paper's conclusions unreliable.\n(2) The plant tissues used for LT50 measurements are ambiguous.\nThe meta-analysis compiled LT50 values from diverse locations and plant species. A critical question arises: Which specific plant tissues were used for LT50 measurements? Was it buds, flowers, leaves, or twigs? Crucially, freezing resistance differs markedly between tissues. Twigs generally exhibit higher resistance compared to leaves and flowers. The significance of freezing risk during spring phenology lies in the vulnerability of newly developing tissues, such as leaves or flowers. Consequently, the defined safety margin (SFM) lacks clarity: for which specific tissue is this margin measured? Furthermore, this study defined spring phenology in terms of the emergence of new tissues. Do new tissues refer to leaf, flowering, or bud? The definition of spring phenology (e.g., leaf unfolding or flowering) should be explicitly stated, and the LT50 measurements used to calculate the SFM must correspond to the specific vulnerable tissues associated with the phenological event.\n(3) The models used for spring phenology prediction may overestimate phenological advancement.\nThe authors employed one-phase models to predict future shifts in spring phenology. A major limitation of these models is that they overlook potential future changes in chilling conditions, which are crucial for releasing dormancy. Under future climate warming, reduced chilling accumulation in winter is expected to increase forcing requirements in spring, which could substantially mitigate the advancing trends of spring warming. Relying solely on one-phase models is, therefore, methodologically unsound and likely overestimates the magnitude of future spring phenological advancement. This overestimation impacts the validity of the study's predictions and conclusions regarding future spring frost risk.\n(4) LT50 is likely to change under future climate warming.\nThe assumption that LT50 remains constant under future warming is merely an assumption and is potentially flawed. Warming conditions may lead to insufficient cold acclimation in plants, increasing the LT50, which indicates a reduction in freezing resistance. This reduction would directly act to narrow the safety margin, a factor that is not accounted for in this study.\n(5) The frequency distribution of freezing resistance values shown in Fig. 2c warrants careful examination. The reported mean is -12.42°C. However, the bulk of the distribution appears to lie significantly to the right of this mean value.\n(Remarks on code availability)"}, "author_response": "Point-to-point responses to reviewer comments\nBelow comments from reviewers are copied in italics, followed by our responses\nin roman font. Additional or altered text in the revised manuscript, when cited in this\nresponse document, is marked in blue.\nTo Reviewer #1\n[General Comment] The authors provide interesting and important new information\non how frost risk affects spring phenology, especially under climate change scenarios.\nSuch information is important to increase the understanding of plant tolerance to\nextreme weather events and how this affects biodiversity. However, despite the\nrelevance of the study, I believe that the authors still need to make significant\nimprovements to the text to increase the clarity of several aspects discussed. I also\nrecommend an English revision. Below I provide comments that I hope will help the\nauthors improve the text.\n[Response] We are grateful to the reviewer for his/her valuable feedback, which has\nsignificantly improved our manuscript. Our point-by-point responses are provided\nbelow, and we hope the revisions fully address all concerns.\nComments in Introduction\n[Comment 1] L1-L24: I recommend restructuring the paragraph, prioritizing a more\norderly sequence of the ideas presented. The structure with which different\nmagnitudes of processes are presented (e.g., species distribution, spring phenology,\nclimate change, trade-off) is a bit confusing. The authors should follow a more\norderly direction in describing the different processes and phenomenons and the\ninterconnection between them. In other words, lead the paragraph from smaller-scale\nprocesses to larger-scale processes, ending by highlighting the practical implications\nof these processes.\n[Response] In response to the reviewers' comments,, we have restructured the\nintroduction. First, we now clearly articulate the current knowledge gap:“In recent\ndecades, advanced spring phenology with warmer springs has been widely reported4,5.\nHowever, it remains unclear whether this phenological advance, accompanied by\nwarmer temperatures, alters the risk of lethal frost (LFR).”(Line 25-27 on Page 3).\nSecond, we introduce the evolutionary trade-off that plants face, “The timing of\nspring phenology under climate warming arises from adaptive trade-offs between two\ncompeting phenological strategies. Emerging earlier may provide a competitive\nadvantage by extending the photosynthetic period and maximizing light-use\nefficiency and photosynthetic activity11,12, reducing herbivore/pathogen pressure13,14,\nand providing priority access to ecological niches15, all of which may enhance plant\nfitness and reproductive success16. On the other hand, potentially lethal spring frost\nevents favor delayed emergence, which protect vulnerable tissues17 and enabling\ngreater resource allocation to growth and reproduction instead of frost tolerance18,19\n(Fig. 1a, Supplementary Table 1).”(Line 30-39 on Page 3), and quantified this\nevolutionary balance, “Specifically, this balance manifests as a buffering mechanism:\nsafeguarding against lethal frost damage during spring phenological stages while\nsimultaneously maximizing individual fitness. Such trade-offs broadly shape plant\nevolutionary strategies and explain spring phenological variation across species and\nregions18,20. Therefore, the optimal phenological strategy for plants under climate\nchange is to maintain this evolutionary balance. Lenz et al.9 proposed a metric termed\nthe ‘safety margin’, defined as the difference between the daily minimum temperature\nand freezing resistance during spring emergence, which operationalizes and quantifies\nthis evolutionary balance (Fig. 1a). The key advantage of this dynamic metric lies in\nits incorporation of biotic factors (e.g., plant freezing resistance), rather than relying\nsolely on abiotic factors21,22,23,24,25.” (Line 45-55 on Page 4). In the end, we introduce\nthe potential risk scenarios arising from climate-change-induced phenological shifts\nand warmer springs, which directly defines the critical knowledge gap our research\ntargets. “Climate-change-induced shifts in spring phenology may fundamentally alter\nthis risk-reward trade-offs in natural plant populations22,26, leading to three potential\ndistinct risk scenarios (Fig. 1b). In the high-risk scenario, earlier springs confer\ntransient competitive advantages, but a single extreme climatic event beyond a critical\nthermal threshold7 could result in local extinction (Fig. 1b). The low-risk scenario\nallows for near-complete frost avoidance, but comes at the cost of a competitive\ndisadvantage that may lead to replacement by other species (Fig. 1b). The stable-risk\nscenario, which best aligns with the EES, illustrates the fundamental trade-off\nbetween advancing growth and avoiding lethal frost damage (Fig. 1b). However,\nempirical data supporting the existence of a spring phenology risk-reward trade-off\nunder warming conditions remain scarce, rendering this theoretical framework highly\nspeculative.” (Line 62-73 on Page 4-5).\n[Comment 2] Note: correct Supplementary Table 1: 3. Less damage by heibivory….\nIt is herbivory.\n[Response] Noted and corrected.\n[Comment 3] L33-34: ambient minimum temperatures. I recommend adopting\nanother term since this one can lead to some misinterpretations.\n[Response] We have replaced ‘ambient minimum temperatures’ with ‘daily minimum\ntemperature’ in the relevant passages.“Lenz et al.9 proposed a metric termed the\n‘safety margin’, defined as the difference between the daily minimum temperature\nand freezing resistance during spring emergence, which operationalizes and quantifies\nthis evolutionary balance (Fig. 1a).”(Line 50-53 on Page 4).\nComments in Results\n[Comment 1] L72-75: Exactly. It is what I mentioned about the caption of your\nsupplementary table 2. 159-161: You have not found strong evidence for evolutionary\nconservation of frost tolerance. Be careful with that statement.\nIn the section 'Phenological models and future climatic data' you mention: \"We\nassumed that freezing resistance will exhibit high evolutionary conservatism ...\". But\nit was not the case according to the results presented in your supplementary table 2.\nTherefore, rephrase the sentences about that.\nIn the section \"Statistical Analysis\", at the end of the first paragraph you mentioned:\n\"When Blomberg's K was statistically significant (p < 0.05), this indicated that\nspecies’ evolutionary history contributed to variance in freezing resistance and the\nfreezing safety margin\". Be careful with this statement. Not necessarily K would\nindicate that evolutionary history explains variance in FR. Rather, it only suggests\nthat phylogenetic relatedness might influence variance in a trait or niche above the\nexpected under a random fashion. However, your results don’t point to that, but that\nFR variance related to the phylogeny of the species is not different from the expected\nby chance. Further, you don’t mention here which evolutionary model you used to\ncompute K., i.e., Brownian Motion? Ornstein-Uhlenbeck (OU)? The type of\nevolutionary model assumed to compute K also can be very influential in the output of\nits p-values. You must account for that in your study while informing and justifying\nclearly your choice to the reader. A last thing I mention, so you reflect about it and\nalso account for that in your study: were average values of FR computed from\nsufficient numbers of individual values? That is, a reliable number of repetitions to\ncompute average per species. Highly variable traits or niches in plants’ life history\nmust be calculated from a robust sampling effort, so it can reflect consistent average\nconditions of them for each species. Therefore, you must also account for this and\ninclude such information in the manuscript as well. Table 2: Your p-values in the\ntable don’t confirm that. Rather, the non-significant p-values you have here, related\nto low K values, indicate that phylogenetic relatedness influence on FR and SFM\ndoesn’t differ from the expected by chance.\n[Response] We are grateful for these constructive comments. In response, we have (1)\nrevised the incorrect statement as suggested, and (2) provided a more detailed\nexplanation of Blomberg's K in the Methods.\n(1) “ The observed patterns of both freezing resistance (K = 0.059, p = 0.144) and\nsafety margin (K = 0.048, p = 0.409) did not differ significantly from the null\nexpectation under a Brownian motion model, indicating that the influence of\nphylogenetic relatedness did not differ from what would be expected by chance (Fig.\n2d, e).” (Line 93-97 on Page 6).\n(2) “Next, we assessed whether evolutionary history was a key predictor of freezing\nresistance and the safety margin using Blomberg's K (R package 'phytools',\nversion\n4.1.3)63. Blomberg's K is a statistical metric that quantifies the phylogenetic signal in\nspecies traits by comparing their observed variance against a Brownian motion model\nof evolution (Eq. 15):\nconcordance between trait and phylogenet ic relationsh ips\nBlomberg' s K  (15)\nconcordance expected under a Brownian motion model\nUnder this model, a value of K > 1 indicates that closely related species are more\nsimilar than expected under Brownian motion, reflecting a strong phylogenetic signal.\nConversely, a value of K < 1 signifies that trait variation is more random than\nexpected, making the pattern indistinguishable from randomness and demonstrating a\nweak or absent phylogenetic signal. However, the K value alone is insufficient and\nmust be accompanied by a significance test to determine if it deviates from random\nexpectation. We therefore employed a permutation test (n = 5000) to generate a null\ndistribution of K values (Fig. 2d, e). The observed K was compared against this\ndistribution, with significance assessed at 0.1 level, a more lenient threshold applied\nconsistently across all subsequent analyses in this study.” (Line 555-568 on Page\n31-32).\n[Comment 2] Figure 1: Define ‘DOY’ of the X-axis in the caption.\n[Response] Done as suggest.\n[Comment 3] Figure 2: Since the phylogeny of the evaluated groups is very\nimportant information for the study, I recommend that you separate the cladogram\nfrom this figure, giving it its own figure (it can be supplementary).\n[Response] Following your suggestion, we have rebuilt Fig. 2 and employed\nConvergent Cross Mapping (CCM) analysis to evaluate the potential causal\nrelationships between freezing resistance and the safety margin.“We applied\nConvergent Cross Mapping (CCM)64, an algorithm designed to detect causal\nrelationships, to explore the potential causality between freezing resistance and the\nsafety margin. The method is grounded in the principle that if X causally influences Y,\nthen the state of X can be reconstructed from the historical record of Y, because the\nevolutionary trajectory of Y embeds information about X65,66. We performed the\nanalysis using the ‘rEDM’ package in R. We tested library sizes from 31 to 304 (τ = 1)\nto evaluate the convergence of prediction correlation (Rho). The optimal embedding\ndimension (E = 5) was determined via simplex projection over a range of dimensions\n(E = 1-20) (Supplementary Fig. 2). We assessed convergence using the twin-surrogate\nmethod with 10 bootstrapped iterations, and considered a causal relationship\nsignificant at 0.1 level.”(Line 568-579 on Page 32).\nFig. 2 | Global dataset of plant freezing resistance and safety margins during\nspring emergence. a-b Geographic distribution and environmental characteristics of\nstudy sites. Circle size indicates species number per site. c Phylogenetic relationships\namong study species. Branch colors represent taxonomic orders; blue and orange bars\nshow freezing resistance and safety margin for each lineage. d-e Frequency\ndistributions of freezing resistance and safety margin for all samples; arrows indicate\nmean values. f Convergent Cross Mapping (CCM) analysis showing cross-mapping\ncorrelation coefficients (Rho) versus library size. Lines and shading represent mean ±\nSD from 10 bootstrap iterations. Convergence of p-values (twin-surrogate method)\nindicates strength and direction of causal relationships; p < 0.1 indicates significant\npotential causality. See Methods for details on Blomberg's K and CCM analysis.\n[Comment 4] Figure 3: P < 0.01? It is unclear.\n[Response] We have updated Fig. 3 to display the actual p-values for all groups.\nFig. 3 | Freezing resistance and safety margin across data subsets. Boxplots show\nmean (solid line), quartiles (25th-75th percentiles, colored boxes), and range\n(whiskers). Jittered points indicate individual samples. Sample sizes shown in\nparentheses. P-values from linear mixed models; bold indicates p < 0.1 (see Methods).\nAbbreviations: REL: relative electrolytic leakage; DTA: differential thermal analysis;\nVOB: visual observation (see Supplementary Table 2).\nAll analyses in this study used a significance threshold at 0.1 level, as stated in the\ntext. “...with significance assessed at 0.1 level, a more lenient threshold applied\nconsistently across all subsequent analyses in this study.” (Line 567-568 on Page 32).\nComments in Discussion\n[Comment 1] L162-193: This paragraph is just too big and confusing. You mix\ndifferent ideas here rather than present a logical order of reasoning. For instance,\nyou start explaining the methodological limitations of DTA, but then start to mix other\nconcepts without a clear connection to it. I recommend you split the paragraph for\nclarity.\nL194: Firstly, as I mentioned before, your discussion in the previous paragraph is\nconfusing. Second, after improving the mentioned paragraph, stress how and what\nyour discussion exactly outlines the predictors of freezing resistance. Give the reader\na clear sense of the novelty with your results.\n[Response] Following your suggestion, we have revised the text to address the\nconcerns, “We found that T , DTR, and T were key predictors of freezing\nmax min\nresistance; however, the linear mixed model showed no significant relationship with\nDTR. This discrepancy is likely because distinct freezing-resistance strategies (e.g.,\nfreezing avoidance versus tolerance) entail different biochemical requirements6.\nVariations in methods used to assess freezing resistance can lead to inconsistent\nresults. For instance, in Differential Thermal Analysis (DTA), the rapid biophysical\nchanges during metastable equilibrium may compromise measurement reliability.\nThis underscores the need for further research to validate and reconcile these\nmethodological discrepancies28.” (Line 181-189 on Page 10), and we have revised the\ncorresponding section as you recommended. “The correlation between freezing\nresistance and various climatic factors, including temperature, evapotranspiration, and\nsolar radiation, fits within an eco-physiological framework. This framework posits\nthat the development of freezing resistance is shaped by the interplay of intrinsic\nphysiological constraints and extrinsic ecological constraints. In terms of intrinsic\nfactors, plant biochemistry, particularly the concentrations of non-structural\ncarbohydrates (NSCs), such as antifreeze proteins (AFPs) and soluble sugars, is\nessential for freezing resistance. Notably, NSC concentrations in plant tissues are\nhigher in colder regions29 and also in areas with greater diurnal temperature\nvariation30. This physiological pattern offers a mechanistic explanation for why both\nT and DTR emerged as strong predictors of freezing resistance in our analysis19,31.\nmean\nIn addition, plant species inhabiting colder regions and those with earlier leaf-out\ndates typically experience more severe environmental conditions, necessitating\ngreater freezing resistance7,22. This may explain the observed positive correlations\nbetween freezing resistance and both T and Phe-T . Extrinsic factors influence\nmin min\nfreezing resistance primarily through cold acclimation and de-acclimation processes.\nFor example, non-lethal low temperatures promote cold acclimation and enhance\nfreezing resistance, whereas warmer spells trigger de-acclimation, reducing freezing\nresistance6,32. This explains why T was identified as a key predictor of freezing\nmax\nresistance in our study. We also observed positive correlations between freezing\nresistance and both evapotranspiration and solar radiation. The evapotranspiration\nrelationship likely reflects a physiological trade-off, higher transpiration rates can\nincrease tissue water content and metabolic activity, thereby reducing freezing\nresistance19, 22. Similarly, under high light conditions, plants preferentially allocate\nmore resources to growth versus defense mechanisms33,34, potentially explaining the\nnegative effect of solar radiation on freezing resistance.” (Line 190-215 on Page\n10-11).\n[Comment 2] L218: What is this biochemical foundation required for freezing\nresistance? Clarify it.\nL219-221: How can it explain the mentioned negative relationship? Clarify it for the\nreader.\n[Response] We have clarified these points as suggested. “Frequent frost exposure can\nreduce photosynthetic productivity24, potentially undermining the biochemical\nfoundation (photosynthates, such as sugars, proteins, enzymes, and other\nnon-structural carbohydrates) necessary for freezing resistance. This relationship\nbetween frost exposure and photosynthetic rate may explain the strong-to-moderate\nnegative relationships we observed between the safety margin and both freezing\nresistance and FRS.” (Line 238-244 on Page 12-13).\n[Comment 3] L237-238: What is the main implication of that for the performance\nand reliability of your results? Make it clear to the reader.\nL263: Make it clearer how and why it remains equivocal.\nL286-288: Exemplify shortly these approaches.\nL288-290: I recommend the authors to create a section “Caveats, limitations, and\nimplications” for addressing all blindspots or limitations of this study, as well as\nindicating future prospects that can assist in solving these limitations.\n[Response] We have incorporated additional analysis using two-phase phenological\nmodels to prevent potential bias and enhance the robustness of our findings (Line\n494-529 on Page 28-30). Furthermore, following your constructive suggestion, we\nhave added a new section titled \"Caveats, Limitations, and Implications\" to\ncomprehensively address all raised concerns. “Our study represents a significant step\ninto the trade-off between spring phenology and lethal frost risk, providing a fresh\nperspective on predicting future spring phenology. Nevertheless, several limitations\nmust be acknowledged. First, plant phenology and freezing resistance data remain\nscarce for regions such as Africa, the Eurasian Steppe, and the Southern Hemisphere\n(Fig. 1a), limiting our ability to generate robust global predictions of terrestrial plant\nphenology under future climate warming scenarios. Second, fundamental challenges\nin modeling spring phenology persist due to two key limitations: (1) discontinuous\nphenological series can undermine the accuracy of future phenological predictions; (2)\noverly empirical parameterizations in current phenological models47. Third, and most\nimportantly, there is a notable lack of data on how warming affects plant freezing\nresistance, especially in trees, as implementing simulated warming on trees is\nparticularly challenging. Addressing such constraints will require innovative\napproaches. Promising directions such as long-term monitoring of freezing resistance\ndynamics, which integrate controlled manipulations with advanced physiological\nmeasurements.” (Line 294-309 on Page 15-16).\nComments in Methods\n[Comment 1] Further in the \"Statistical Analysis\" section you mention: \"To interpret\nthe mixed model results, instead of significance testing with an arbitrary p-value cutoff,\nwe implemented the gradual language of evidence recommended by Muff et al., this\napproach allowed us to communicate model findings in a more nuanced way\". Explain\nmore about this approach since it is a vital information for your study. Include a simple\nsentence explaining in slightly more details.# Data availability and Code\navailabilityThe URL links provided don't work. Update the links.\n[Response] We are grateful for your meticulous review. In response, we have added\ncitations after the method description to indicate which figures and/or tables\nsubsequently employ each approach. “To interpret the mixed model results, instead of\nsignificance testing with an arbitrary p-value cutoff, we implemented the gradual\nlanguage of evidence recommended by Muff et al.72 (Supplementary Tables 4, 5, 6, 7),\nthis approach allowed us to communicate model findings in a more nuanced way16.”\n(Line 607-610 on Page 34).\n[Remarks on code availability] The codes are functional. However, they only\nprovide commands to generate figures shown in the manuscript, but there is not\nenough code to portray the statistical background used.\nI recommend that authors provide all the code routines used in data processing and\ntreatment (along with the data), as well as the statistical modeling, thus ensuring the\nreplicability of their work.\n[Response] We have updated the repository link, which now hosts all code used in\nthis study. This includes the code to generate all main and appendix figures,\nimplement all models, reproduce the statistical analyses, and the raw data. We believe\nthis will fully address your concerns and enhance the replicability of our work:\nhttps://doi.org/10.6084/m9.figshare.28877072.v3.\nTo Reviewer #2\n[General Comment] This manuscript addresses the important trade-off between\nextending the growing season and avoiding frost risk, mediated by spring phenology.\nFurthermore, this manuscript investigated how warming-induced phenological shifts\nmay alter risk-reward trade-offs under future climate scenarios. This work appears to\nprovide novel findings about the frost risk of plants, but the conclusions rely on\nseveral highly questionable assumptions, such as the LT50 remaining constant across\nseasons, plant tissues, and future climate scenarios. These assumptions significantly\nundermine the credibility of the paper's conclusions. Specific comments are detailed\nbelow:\n[Response] We sincerely thank the reviewer for the meticulous review and valuable\ncomments. To strengthen the manuscript, we have (1) incorporated results from our in\nsitu warming experiment to support our central hypothesis of stable freezing\nresistance; (2) incorporated four chilling-based phenological models to strengthen and\nsolidify our major findings, and (3) added a dedicated \"Caveats, Limitations, and\nImplications\" section to address the study's scope and limitations. Our point-by-point\nresponses are provided below.\n[Comment 1] The assumption that LT50 is temporally constant across seasons\ncontradicts empirical observations.\nIn this study, LT50 is assumed to be constant over time (depicted as the fixed\nblue dashed line in Figure 1). However, experimental evidence demonstrates that\nfreezing resistance undergoes continuous seasonal variation. Induced by short\nphotoperiods and low temperatures in autumn, freezing resistance increases\nprogressively—a process termed cold acclimation. Freezing resistance typically\nreaches its peak (i.e., lowest LT50) in December or January and subsequently\ndeclines, reaching its minimum level (i.e., highest LT50) during the active growing\nseason. This seasonal fluctuation in LT50 often exceeds 10 to 20°C. Applying the\nsame LT50 value to assess freezing resistance at different spring phenology timing is\ntherefore unsound and directly renders the paper's conclusions unreliable.\n[Response] We thank the reviewer for this valuable suggestion. We acknowledge that\nour original text may have insufficiently emphasized the study's primary goal. To\nclarify, we did not assume LT to be constant; our objective was specifically to\n50\nquantify the safety margin during spring phenology. Accordingly, our meta-analysis\nexclusively used freezing resistance data from this period, as detailed in the Methods.\n“(ii) freezing tests were conducted during spring phenology on new leaves/buds using\nstandardized methods (Supplementary Fig. 1)”. (Line 383-385 on Page 23). We have\nremoved the blue dashed line from Fig. 1a to enhance clarity.\n[Comment 2] The plant tissues used for LT50 measurements are ambiguous.\nThe meta-analysis compiled LT50 values from diverse locations and plant\nspecies. A critical question arises: Which specific plant tissues were used for LT50\nmeasurements? Was it buds, flowers, leaves, or twigs? Crucially, freezing resistance\ndiffers markedly between tissues. Twigs generally exhibit higher resistance compared\nto leaves and flowers. The significance of freezing risk during spring phenology lies in\nthe vulnerability of newly developing tissues, such as leaves or flowers. Consequently,\nthe defined safety margin (SFM) lacks clarity: for which specific tissue is this margin\nmeasured? Furthermore, this study defined spring phenology in terms of the\nemergence of new tissues. Do new tissues refer to leaf, flowering, or bud? The\ndefinition of spring phenology (e.g., leaf unfolding or flowering) should be explicitly\nstated, and the LT50 measurements used to calculate the SFM must correspond to the\nspecific vulnerable tissues associated with the phenological event.\n[Response]We thank the reviewer for this comment. Our assessment specifically\ntargeted newly emerging tissues (leaves and buds), as we sought to capture the freezing\nresistance during spring phenology. Overwintering organs like branches were\ndeliberately excluded due to their typically higher freezing resistance, as stated in our\nmethods: “(ii) freezing tests were conducted during spring phenology on new\nleaves/buds using standardized methods (Supplementary Fig. 1).” (Line 383-385 on\nPage 23).\n[Comment 3] The models used for spring phenology prediction may overestimate\nphenological advancement.\nThe authors employed one-phase models to predict future shifts in spring\nphenology. A major limitation of these models is that they overlook potential future\nchanges in chilling conditions, which are crucial for releasing dormancy. Under\nfuture climate warming, reduced chilling accumulation in winter is expected to\nincrease forcing requirements in spring, which could substantially mitigate the\nadvancing trends of spring warming. Relying solely on one-phase models is, therefore,\nmethodologically unsound and likely overestimates the magnitude of future spring\nphenological advancement. This overestimation impacts the validity of the study's\npredictions and conclusions regarding future spring frost risk.\n[Response] We are grateful for this valuable suggestion. In response to the point\nregarding chilling's key role in spring phenology, we now simulate spring phenology\nusing the chilling-based Parallel model, incorporating four different models to\ncalculate chilling accumulation. “The two-stage phenological models are based on the\npremise that spring phenology is initiated only after sufficient chilling accumulation.\nWe applied the commonly used Parallel model, a chilling-based framework, to\nsimulate future spring phenology39. The state of chilling (S ) is modeled as:\nc\nD\nS (D)  R (T) (8)\nc c\ntt1c\nwhere R determines the rate of change in S , t is the first day of the chilling\nc c 1c\nperiod, and D is the day of the year. Spring phenological events are triggered when S\nf\n(Eq. 2)and S (Eq. 8) reaches threshold (Eq. 9):\nc\nif S(D) ≥ F* then leaf induction is completed\nf\nif S (D) ≥ C* then dormancy is completed (9)\nc\nIn this study, we employed four chilling models to quantify the chilling\naccumulation required for spring phenology47.\nThe first model C , also known as the Utah Model, incorporates differential\n1\ntemperature weighting (Eq. 10). It was originally proposed for quantifying the chilling\nrequirements of peach trees60.\n0 T 1.4\n\n0.5 1.4 < T  2.4\n\n1 2.4 < T 9.1\n\nCU\n1\n 0.5 9.1< T 12.4 (10)\n\n0 12.4 < T 15.9\n\n0.5 15.9 < T 18\n\n1 T >18\n\nAfter removing the negative contributions of warm temperatures to the\naccumulated chilling, we obtained the modified Utah Model, designated as model C\n2\n(Eq. 11)61:\n0 T 1.4\n\n0.5 1.4 < T  2.4\n\n\nCU\n2\n 1 2.4 < T 9.1 (11)\n\n0.5 9.1< T 12.4\n\n0 T >12.4\n\nwhere CU and CU is the rate of chilling for model C and C , respectively. T is\n1 2 1 2\nthe daily mean temperature (℃).\nNote that models C1 and C2 operate on the premise that no relevant plant\nphysiological activity occurs below freezing. Therefore, sub-zero temperatures are\nomitted from their chilling calculations. However, Zhang et al.62 fitted the triangular\nmodel to observational data from 24 plant species and found that a mean optimal\nchilling temperature of 0.2℃ and an upper chilling threshold of 6.9℃ were most\neffective. Consequently, we adopted a triangular function with these parameters (0.2℃\nand 6.9℃) as our model C (Eq. 12):\n3\nT 6.5\n-6.5< T 0.2\n\n6.9-0.2\n\nCU\n3\n 0 T -6.5 or T 6.9 (12)\n\n6.9-T\n 0.2< T 6.9\n6.9-0.2\nModel C also adopts a triangular form, as defined by Hänninen58 through fitting\n4\nto previous experimental results for Finnish birch seedlings (Eq. 13):\n1 2.5< T 7.4\n\n 2.10 T4.66 3.10\nCU\n4\n\n\n3.13(\nT4.66\n)\ne-(\n10.93\n)\nT -4.7 or T 16 (13)\n10.93\n\n\n0 else\n\n\nwhere CU and CU is the rate of chilling for model C and C , respectively. T is\n3 4 3 4\nthe daily mean temperature (℃). Note that both model C and model C permit\n3 4\nchilling accumulation to occur at sub-zero temperatures.” (Line 494-529 on Page\n28-30).\nThe results are consistent with our previous findings (Supplementary Fig. 7, 12),\nthereby reinforcing our major conclusion (Fig. 6).\nSupplementary Figure 7 | Projections of future spring phenology and models\ncomparison. Projections use CMIP6 general circulation models for three Shared\nSocioeconomic Pathways: SSP1-2.6 (low warming), SSP2-4.5 (medium warming),\nand SSP5-8.5 (high warming). Significant findings (at 0.1 level) are highlighted in\nbold and marked with an asterisk, p-values were determined using linear mixed model\ncompared to zero, see methods for more details.\nSupplementary Figure 12 | Projections of future safety margin and models\ncomparison. Projections use CMIP6 general circulation models for three Shared\nSocioeconomic Pathways: SSP1-2.6 (low warming), SSP2-4.5 (medium warming),\nand SSP5-8.5 (high warming). Significant findings (at 0.1 level) are highlighted in\nbold and marked with an asterisk, p-values were determined using linear mixed model\ncompared to zero, see methods for more details.\n[Comment 4] LT50 is likely to change under future climate warming.\nThe assumption that LT50 remains constant under future warming is merely an\nassumption and is potentially flawed. Warming conditions may lead to insufficient\ncold acclimation in plants, increasing the LT50, which indicates a reduction in\nfreezing resistance. This reduction would directly act to narrow the safety margin, a\nfactor that is not accounted for in this study.\n[Response] We thank the reviewer for this constructive comment. While future\nwarming may alter plant freezing resistance, direct empirical evidence is scarce. To\naddress this and test our hypothesis of stable freezing resistance, we conducted an in\nsitu experiment in 2023. The results from 15 perennial herbs show no significant\nchange in freezing resistance (Supplementary Fig. 11), providing direct support for\nour hypothesis.\nSupplementary Figure 11 | Freezing resistance under ambient and warming\nconditions for fifteen perennial herbs in Qinghai Haibei National Field Research\nStation. Jittered points represent independent samples for each species, ns indicate\nnon-significante at 0.1 level using permutation test. Pb, Potentilla bifurca (p = 0.45);\nBv, Bistorta vivipara (p = 0.69); Tm, Taraxacum mongolicum (p = 0.3); Lt, Lancea\ntibetica (p = 0.89); Vb, Viola bulbosa (p = 0.78); It, Iris tectorum (p = 0.33); Ll,\nLeontopodium leontopodioides (p = 0.56); Th, Tibetia himalaica (p = 0.34); Tl,\nThermopsis lanceolata (p = 0.6); Ac, Anemone cathayensis (p = 0.75); Rn,\nRanunculus nephelogenes (p = 0.93); Oa, Oxytropis alpina (p = 0.62); Gs, Gentiana\nstraminea (p = 0.75); Aa, Aster alpinus (p = 0.74); Sn, Saussurea nigrescens (p =\n0.58).\nWe have expanded the Methods section to include more detailed information\nabout our in-situ experiment. “A central objective of this study is to project the future\nrisk of lethal frost. This is achieved by integrating projected shifts in spring phenology\nand the potential changes in freezing resistance under warming scenarios. Such an\nintegrated approach is critical, as the direct impact of climate warming on freezing\nresistance itself remains poorly understood. To evaluate the freezing resistance during\nspring phenology under warming conditions, we selected 15 perennial herbaceous\nspecies at the Qinghai Haibei National Field Research Station of Alpine Grassland\nEcosystems (37°29′-37°45′N, 101°12′-101°23′E, 3200m a.s.l., see our previous\nstudies for site details50,51,52). The warming treatment increased air temperature by an\naverage of approximately 1.8℃ over a ten-year period. Plant cutting samples (leaves)\nwere cooled to each target temperature (using commercial freezers) and subsequently\nthawed at a rate of approximately 6℃ h-1 until room temperature was reached\n(Supplementary Fig. 10a). Initial assessment of freezing damage was conducted\nthrough visual inspection, such damage typically manifested as tissue discoloration\ndue to polyphenol oxidation, or as a characteristic odor resulting from protoplast\ndecompartmentalization and autolysis9. Following visual inspection, cutting samples\nwere transferred to small beakers containing distilled water for electrical conductivity\n(Ec) measurement. These Ec assays were conducted hourly until values stabilized, a\nprocess typically requiring over ten hours (the slow leakage of electrolyte).\nConcurrently, fresh samples were boiled at 100℃ for 30 minutes and then cooled to\nroom temperature to determine the maximum electrolyte conductance (Ec ). In this\nmax\nstudy, freezing resistance was defined as the lethal temperature (LT , ℃) at which\n50\nelectrolyte leakage reached 50% of EC 53,54,55, and this LT value serves as a\nmax 50\ncommon proxy for the temperature at which 50% of the sample cells died9\n(Supplementary Fig. 10b).” (Line 428-452 on Page 25-26).\nGiven that the response of tree freezing resistance to warming remains poorly\nunderstood, we have acknowledged this gap in a new, dedicated section.\n“Caveats, limitations, and implications\n...Third, and most importantly, there is a notable lack of data on how warming\naffects plant freezing resistance, especially in trees, as implementing simulated\nwarming on trees is particularly challenging. Addressing such constraints will require\ninnovative approaches. Promising directions such as long-term monitoring of freezing\nresistance dynamics, which integrate controlled manipulations with advanced\nphysiological measurements” (Line 304-309 on Page 15-16).\n[Comment 5] The frequency distribution of freezing resistance values shown in Fig.\n2c warrants careful examination. The reported mean is -12.42°C. However, the bulk\nof the distribution appears to lie significantly to the right of this mean value.\n[Response] Thank you for pointing out this error. The erroneous marking in Fig. 2\nhas been corrected."}, {"round": 2, "reviewer_comments": {"reviewer_2": "I appreciate the authors’ effort to revise the manuscript. The attempt to (1) explicitly define the timing and plant tissues associated with freezing resistance measurements, (2) incorporate more rigorous two-stage phenological models that account for chilling accumulation, and (3) include an experimental evaluation of freezing resistance stability during spring phenology under warming conditions has improved the manuscript’s clarity and structure.\nHowever, the revised manuscript still assumes spring freezing resistance remains stable under warming scenarios. The experimental evidence provided to support this assumption includes only herbaceous species. This is insufficient to generalize across all life forms. In woody species, abundant literature demonstrates that freezing resistance starts to decline in late winter or early spring (Vitra et al., 2017; Hu et al., 2022; Grossman, 2023). Moreover, both chilling accumulation and forcing temperature strongly influence the rate of deacclimation in trees (Kovaleski et al., 2018; North et al., 2022). Considering the combined effects of phenological advancement, insufficient chilling, and intensified thermal forcing under warming scenarios, the assumption of invariant freezing resistance lacks empirical support and may lead to inaccurate estimations of safety margins.\nIn addition, the forcing model component employed within the two-stage phenological framework was not clearly specified. This information should be explicitly provided to ensure methodological reproducibility.\nIn summary, although the authors have made notable improvements, the central issue remains unresolved. The assumption of constant freezing resistance under future climate warming still lacks sufficient empirical evidence, thereby undermining the robustness and reliability of the study’s primary conclusions.\nReferences\nGrossman JJ. 2023. Phenological physiology: seasonal patterns of plant stress tolerance in a changing climate. New Phytologist 237: 1508-1524.\nHu Z, Lin S, Wang H, Dai J. 2022. Seasonal variations of cold hardiness and dormancy depth in five temperate woody plants in China. Frontiers in Forests and Global Change 5: 1061191.\nKovaleski AP, Reisch BI, Londo JP. 2018. Deacclimation kinetics as a quantitative phenotype for delineating the dormancy transition and thermal efficiency for budbreak in Vitis species. AoB PLANTS 10: ply066.\nNorth M, Workmaster BA, Atucha A. 2022. Effects of chill unit accumulation and temperature on woody plant deacclimation kinetics. Physiologia Plantarum 174: e13717.\nVitra A, Lenz A, Vitasse Y. 2017. Frost hardening and dehardening potential in temperate trees from winter to budburst. New Phytologist 216: 113-123.\n(Remarks on code availability)"}, "author_response": "Point-to-point responses to reviewer comments\nBelow comments from reviewers are copied in italics, followed by our responses in\nroman font. Additional or altered text in the revised manuscript, when cited in this\nresponse document, is marked in blue.\nTo Reviewer #1\n[General Comment] I appreciate the authors' effort to revise the manuscript. The\nattempt to (1) explicitly define the timing and plant tissues associated with freezing\nresistance measurements, (2) incorporate more rigorous two-stage phenological\nmodels that account for chilling accumulation, and (3) include an experimental\nevaluation of freezing resistance stability during spring phenology under warming\nconditions has improved the manuscript's clarity and structure.\nIn summary, although the authors have made notable improvements, the central issue\nremains unresolved. The assumption of constant freezing resistance under future\nclimate warming still lacks sufficient empirical evidence, thereby undermining the\nrobustness and reliability of the study’s primary conclusions.\n[Response] Thank you for acknowledging our efforts in revising the manuscript, and\nfor your positive feedback on our work. We fully understand your concerns regarding\nthe assumption of constant freezing resistance under future climate warming, which is\nbased on the potential decline in chilling accumulation under future warming\nscenarios.\nTo address this, we conducted additional analyses to quantify chilling\naccumulation using four chilling models and four earth system models. Our results\nalign with those of previous studies, confirming that chilling accumulation is highly\ndependent on the choice of both chilling and climate models. We will discuss this\npoint in detail below and have made corresponding revisions to fully address your\nconcern.\n[Comment 1] However, the revised manuscript still assumes spring freezing\nresistance remains stable under warming scenarios. The experimental evidence\nprovided to support this assumption includes only herbaceous species. This is\ninsufficient to generalize across all life forms. In woody species, abundant literature\ndemonstrates that freezing resistance starts to decline in late winter or early spring\n(Vitra et al., 2017; Hu et al., 2022; Grossman, 2023). Moreover, both chilling\naccumulation and forcing temperature strongly influence the rate of deacclimation in\ntrees (Kovaleski et al., 2018; North et al., 2022). Considering the combined effects of\nphenological advancement, insufficient chilling, and intensified thermal forcing under\nwarming scenarios, the assumption of invariant freezing resistance lacks empirical\nsupport and may lead to inaccurate estimations of safety margins.\n[Response] Thank you for your professional feedback and for sharing the relevant\narticles, which are closely related to our study and have now been cited in our\nmanuscript. We will discuss your concerns in detail below and present the\ncorresponding revisions we have made to address them.\nConclusion first, the assertion that climate warming with earlier spring\nphenology inevitably reduces chilling accumulation is highly empirical and lacks\nrigorous testing. For example, Wang et al. (2020) have indicated that more than half\nof current chilling models are invalid and have demonstrated opposing results\nregarding the interaction between chilling accumulation and heat requirement\n(Response Fig.1).\nResponse Fig.1\nWang, H. J. et al. Overestimation of the effect of climatic warming on spring\nphenology due to misrepresentation of chilling. Nat Commun. 11, 4945 (2020).\nDOI: https://doi.org/10.1038/s41467-020-18743-8.\nFollowing this line of reasoning, we have calculated chilling accumulation using\nthe four chilling models and four Earth System Models (ESMs) employed in our\nstudy (Supplementary Fig. 10). The predictions reveal inconsistent trends in chilling\naccumulation under the three future warming scenarios. Specifically, accumulation\nremains stable under SSP1-2.6 and SSP2-4.5 in most models. Under SSP5-8.5, 7 out\nof 16 projections indicate stable chilling, while the other models predict a reduction\n(Supplementary Fig. 10).\nSupplementary Figure 10 | Projected chilling accumulation under future climate\nscenarios (2021-2099). Projections use CMIP6 general circulation climate models\n(M1-M4) and four process-based chilling models (CU1-CU4) under three Shared\nSocioeconomic Pathways: SSP1-2.6 (low warming), SSP2-4.5 (medium warming),\nand SSP5-8.5 (high warming).\nOur results now provide evidence that freezing resistance remains stable under at\nleast two warming scenarios, as assumed. Furthermore, Ettinger et al. (2020)\nindicated that the stalling effects of chilling generally appear above 4 ℃ (Response\nFig.2). Since the SSP1-2.6 and SSP2-4.5 scenarios do not exceed this threshold\n(Supplementary Figure 6 in appendix), their findings support our statement.\nResponse Fig.2\nEttinger, A. K. et al. Winter temperatures predominate in spring phenological\nresponses to warming. Nat. Clim. Chang. 10, 1137-U119 (2020).\nDOI: https://doi.org/10.1038/s41558-020-00917-3.\nRegarding the SSP5-8.5 scenario, we maintain our hypothesis and cautiously\nconclude that freezing resistance may decline under this condition. This decline may\ncompensate for the extended safety margin, align with our major point, although the\nmagnitude of this effect requires further investigation. In line with this reasoning, we\nhave revised the relevant sections of the manuscript.\nWe have highlighted that the potential reduction in chilling accumulation under\nfuture warming may compromise freezing resistance. Subsequently, we quantified\nchilling accumulation using the four chilling and climate models employed in our\nstudy, thereby providing evidence to support our assumption. “Regarding future\nfreezing resistance, numerous studies have demonstrated that warming leads to\nreduced chilling accumulation, which may impair freezing resistance8,19,28,29,30.\nHowever, the actual impact depends on the future climate scenarios31, data sources,\nand the chilling model used32. Therefore, we projected chilling accumulation for\n2021-2099 using four climate scenarios and four process-based chilling models. Most\nmodels indicate that chilling accumulation was not always diminished under\nwarming, however, half of the results show a strong decline only under the SSP5-8.5\nscenario (Supplementary Fig. 10). Additionally, our in situ observations show no\nsignificant change in the freezing resistance of the 15 perennial herbs under moderate\nwarming (Supplementary Fig. 11, 12, see methods). Despite the inherent\nmethodological limitations of such manipulative studies, our results suggest that plant\nfreezing resistance may remain stable at least two warming scenarios.” (Line 152-164\non Page 8-9).\n“To our knowledge, only two studies to date have directly demonstrated that warming\nreduces freezing resistance, with the effect being contingent on the magnitude of\nwarming51,52. For instance, Ettinger et al. (2019) showed that the stalling effects of\nchilling appear only above 4℃ of warming31, a pattern that aligns with our findings.”\n(Line 285-289 on Page 14-15).\nWe have revised the relevant sentences in the “Caveats, Limitations, and\nImplications” and “Materials and methods” section to better describe the\nlimitations of our study.\n“Third, and most importantly, there is a notable lack of data on how warming affects\nplant freezing resistance, especially in trees under extreme warming (e.g. SSP5-8.5\nscenario), as implementing simulated warming on trees is particularly challenging.”\n(Line 312-315 on Page 16).\n“ We posit that plant frost resistance is likely to remain stable under moderate climate\nwarming scenarios (< 4℃)31. This stability is attributed to the relatively stable chilling\naccumulation under such conditions (Supplementary Fig. 10), a conclusion also\nsupported by our in situ manipulation experiment (Supplementary Fig. 12). ” (Line\n550-554 on Page 31).\n[Comment 2] In addition, the forcing model component employed within the two-\nstage phenological framework was not clearly specified. This information should be\nexplicitly provided to ensure methodological reproducibility.\n[Response] Thank you very much for your meticulous review and valuable\nsuggestion. We have now added a more detailed description of the forcing model used\nin the two-stage phenological models. “The triggering of spring phenological events\noccurs when both S (Eq. 2) and S (Eq. 8) reach their specific thresholds. We adopted\nf c\nthe commonly used model for calculating S (Eq. 4), and employed four chilling\nf\nmodels to quantify the chilling accumulation required for spring phenology32” (Line\n510-513 on Page 28-29).\n[Additional revisions]\n[No. 1] We have now replaced the term 'species' with 'aggressive species' to\nemphasize the fitness loss incurred by adopting a more conservative strategy. “The\nlow-risk scenario allows for near-complete frost avoidance, but comes at the cost of a\ncompetitive disadvantage that may lead to replacement by other aggressive species\n(Fig. 1b).” (Line 66-68 on Page 5).\n[No. 2] We have now added the definition of LT in Supplementary Table 3. “LT ,\n50 50\nLethal temperature at which 50% of the samples died15” (Supplementary Table 3 in\nappendix).\n[No. 3] We have stated that the lower risk of lethal frost (high freezing resistance) is\ndriven by extreme events rather than by regular temperature fluctuations. “First, we\nfound that plants sustain enhanced freezing resistance during the spring emergence,\nresulting in an exceptionally low lethal frost risk, which may have been driven by the\nstress of extreme events.” (Line 182-185 on Page 10)."}, {"round": 3, "reviewer_comments": {"reviewer_2": "The authors have adequately addressed my previous concern regarding the assumption of constant freezing resistance during spring phenology under future warming. Two issues remain to be addressed.\n(1) The chilling models adopted by the authors impose a fixed lower effective temperature threshold (ranging from −6.5 to 1.4 °C). However, one study found that freezing temperatures down to −10 °C can still effectively promote dormancy release in temperate woody species (https://doi.org/10.1093/aob/mcae112). Therefore, chilling models adopted in this paper neglect the contribution of freezing temperatures (below −6.5 °C) to chilling. Under warming conditions, temperatures that were previously below the model’s lower bound may shift into the effective range, potentially leading to an underestimation of the projected decline in chilling. To compare the influence of a lower temperature threshold on chilling estimates, I recommend including an additional chilling model without a lower threshold to test whether chilling accumulation remains stable under future climate scenarios. For example, a model assuming temperatures below an upper threshold (e.g., 10 oC) have equal effects on dormancy release.\n(2) A minor issue remains in the equation of Model C4 (L537). “T≤-4.7 or T ≥16” should be placed after “0”.\n(Remarks on code availability)"}, "author_response": "Point-to-point responses to reviewer comments\nBelow comments from reviewers are copied in italics, followed by our responses in\nroman font. Additional or altered text in the revised manuscript, when cited in this\nresponse document, is marked in blue.\nTo Reviewer #2\n[General Comment] The authors have adequately addressed my previous concern\nregarding the assumption of constant freezing resistance during spring phenology\nunder future warming. Two issues remain to be addressed.\n[Response] Thank you for your professional, meticulous comments and for your\npatience. Your suggestions have been invaluable, significantly improving the\nmethodological rigor of our manuscript during the revision process.\nWe will provide a detailed justification for the effective temperature threshold\nrange (−6.5 to 1.4°C) used in our study. Concerning the study you cited, we believe its\nconclusions are not directly generalizable to our specific research context and dataset,\nfor reasons we will elaborate.\nOnce again, we are deeply appreciative of the significant time and effort you\ninvested in reviewing our manuscript. We feel truly fortunate to have had you as our\nreviewer.\n[Comment 1] The chilling models adopted by the authors impose a fixed lower\neffective temperature threshold (ranging from −6.5 to 1.4°C). However, one study\nfound that freezing temperatures down to −10 °C can still effectively promote\ndormancy release in temperate woody species (https://doi.org/10.1093/aob/mcae112).\nTherefore, chilling models adopted in this paper neglect the contribution of freezing\ntemperatures (below −6.5 °C) to chilling. Under warming conditions, temperatures\nthat were previously below the model’s lower bound may shift into the effective range,\npotentially leading to an underestimation of the projected decline in chilling. To\ncompare the influence of a lower temperature threshold on chilling estimates, I\nrecommend including an additional chilling model without a lower threshold to test\nwhether chilling accumulation remains stable under future climate scenarios. For\nexample, a model assuming temperatures below an upper threshold (e.g., 10°C ) have\nequal effects on dormancy release.\n[Response] Thank you for this excellent and insightful suggestion. The topic of\nchilling accumulation at subfreezing temperatures is indeed both fascinating and\nimportant, as it has not been extensively studied in the existing literature.\nWe would like to clarify two points in response. First, it is important to note that\nthis specific issue currently lacks rigorous physiological testing in the literature.\nSecond, while the study you cited is valuable, its conclusions may be influenced by\nmethodological constraints related to sampling and design. A core strength of our\nmeta-analysis is its ability to mitigate these biases by synthesizing a more\ncomprehensive and representative dataset.\n How do plants maintain physiological activity at subfreezing temperatures?\nIt must be emphasized that, whether for forcing or chilling phases, should be\ninterpreted as a proxy for plant physiological processes, and not simply as a purely\nmeteorological abstraction. It is well known that pure water freezes at 0°C. This phase\ntransition results in a loss of molecular mobility, preventing the diffusion and\nencounter of ions and thus halting all metabolic activity. However, plant cell sap is not\npure water; it contains various solutes that depress the freezing point. For instance,\nKörner et al. (2016) reported that a 1-molar solution depresses the freezing point by\n1.8 K. Notably, the cells of most plants function with a cell sap concentration at or\nslightly below 1 molar.\nWhen a cold snap occurs, this basal freezing tolerance becomes inadequate. In\nthe short term, as the extracellular fluid (which typically has a lower solute\nconcentration than the cytosol) begins to freeze, the resulting ice creates a osmotic\ngradient, subjecting the cell to dehydration risk. In response, cells accumulate\nnon-structural carbohydrates, further raising the cytoplasmic solute concentration.\nThis acclimation mechanism may explain why fruits from colder regions (or those\nexperiencing greater diurnal temperature variation) and scallions harvested in late\nautumn are often perceived as sweeter (Tan et al., 2023; Blumstein et al., 2022). In the\nlong term, sustaining freeze resistance through constant solute production is\nunsustainable due to its high metabolic cost. Consequently, plants trigger a\ndormancy-prep dehydration response. This dehydration concentrates the cytoplasm,\nboosting overall freezing tolerance. Thus, in nature, the physiological\nlow-temperature limit during brief cold events is generally maintained within the -5°\nC to -10°C range.\nMetabolism is critically dependent on enzyme activity. Although cold\nacclimation confers tolerance via supercooling, maintaining liquid cytosol even at\ntemperatures as low as -40°C (Körner et al., 2016). Low temperatures severely limits\nmetabolic and metabolism ceases to contribute meaningfully to chilling accumulation.\nThis fundamental physiological limitation is precisely what our CU4 model\nquantifies.\nThe aforementioned mechanisms account for plant survival under extreme cold\nbut also impose a severe constraint on chilling accumulation. Furthermore, the\ninactivation of antioxidant enzymes at low temperatures further exacerbates stress by\nincreasing the risk of reactive oxygen species damaging the cell wall (Grossman,\n2022). This damage, in turn, contributes to a reduction in photosynthetic capacity of\nup to ~70% compared to optimal conditions (Körner et al., 2016). Cumulatively, this\nevidence demonstrates that chilling accumulation is profoundly limited at low\ntemperatures (e.g., below -10°C), justifying its exclusion from the model. Therefore,\nwe contend that the temperature range used in our chilling model (-6.5 to 1.4°C) is\nboth conservative and physiologically well-founded.\nIn the following section, we will discuss the specific context of the study you\ncited and why its findings may offer valuable insights for that particular system but\nmay have limited generalizability to the broader ecological and methodological\nframework of our analysis.\n What to consider when studying freezing resistance？\nAs established, extra cold acclimation can substantially enhance freezing\ntolerance. To isolate the natural trade-off between spring phenology and lethal frost\nrisk, studies involving human management, such as common garden experiments or\nartificially bred individuals, must be excluded. This criterion has been rigorously\napplied in our meta-analysis. “(i) data were collected at the species level under natural\nconditions, excluding garden or glasshouse studies to avoid human management”\n(line 390-392 in page 23).\nRegarding the study you cited, its sampling was conducted in the Olympic Forest\nPark in Beijing. It is important to clarify that this park is a large-scale, artificially\nconstructed landscape developed between 2003 and 2008 for the Olympic Games. As\na result, the tree, shrub, and floral communities within it consist predominantly of\ncultivated varieties that were artificially selected and transplanted. Their physiological\nresponses thus may not represent the natural adaptive strategies of wild populations\nthat are central to our research question. Furthermore, as a managed urban landscape,\nthe park's vegetation was specifically selected and maintained for cold hardiness.\nProtective measures such as painting trunks white, applying insulation wraps, and\ninstalling windbreaks are routinely implemented for vulnerable species. This intensive\nmanagement artificially buffers the plants from the full selective pressures of natural\nfrost events. Consequently, the observed phenotypes and survival rates in this setting\ndo not reflect the natural trade-offs and evolutionary constraints that are the focus of\nour meta-analysis on wild populations. For such intensively managed species,\nestimates of chilling accumulation at subfreezing temperatures are likely artificially\ninflated.\nMoreover, the cutting-sample methodology employed in that study introduces a\nsecond layer of methodological bias (Vitasse and David, 2014; Ambroise et al., 2019).\nBy placing excised twigs in a controlled greenhouse with artificially optimal\nconditions, such as extended photoperiods, unlimited water, and no wind stress, the\nexperimental design likely overstimulates their physiological response to chilling\ntreatment. This further decouples the observed results from the resource-constrained\nand multi-stress reality of plants in natural field environments, which is the context\nour meta-analysis seeks to understand.\n Conclusion\nIn summary, after careful consideration, we have decided not to incorporate a\nnew chilling model without an upper temperature threshold in the current manuscript.\nWe believe that such a model lacks robustness given the current physiological\nunderstanding and available data. Moreover, introducing an additional,\nless-constrained model would not significantly enhance the core contributions of our\npaper. Instead, it might confuse readers who are not specialists in this specific area,\nthereby hindering the clarity and broad communication of our findings.\nWe thank the reviewer for raising this fascinating and important question, which\nwe agree points to a key future research direction. We appreciate this insightful\nsuggestion and will explore it in subsequent work. Once again, we are grateful for the\nreviewer's dedicated effort and constructive comments.\nReference\nKörner, C. 2016. Plant adaptation to cold climates. F1000Research.\nTan, X. Q. et al. 2023. haracterization of the Difference between Day and Night\nTemperatures on the Growth, Photosynthesis, and Metabolite Accumulation of\nTea Seedlings. Int. J. Mol. Sci.\nBlumstein, M. et al. 2022. Global variation in nonstructural carbohydrate stores in\nresponse to climate. Glob Chang Biol.\nGrossman, J. J. 2022. Phenological physiology: seasonal patterns of plant stress\ntolerance in a changing climate. New Phytol.\nVitasse, Y. & Basler, D. 2014. Is the use of cuttings a good proxy to explore\nphenological responses of temperate forests in warming and photoperiod\nexperiments? Tree Physiol.\nAmbroise, V. et al. 2019. The Roots of Plant Frost Hardiness and Tolerance. Plant\nand Cell Physiology.\n[Comment 2] A minor issue remains in the equation of Model C4 (L537). “T≤-4.7\nor T ≥16” should be placed after “0”.\n[Response] Done as suggested.\n0 else\n\n2.10 T4.66 3.10\nCU\n4\n\n\n3.13(\nT4.66\n)\ne-(\n10.93\n)\nT-4.7 or T16 (13)\n10.93\n\n1 2.5<T7.4\n"}]

Decision letter