Cold-season freeze frequency is a pervasive driver of subcontinental forest growth.

As northern latitudes experience rapid winter warming, there is an urgent need to assess the effect of varying winter conditions on tree growth and forest carbon sequestration potential. We examined tree growth responses to variability in cold-season (November–April) frequency of freeze days (FFD) over 1951 to 2018 using tree-ring data from 35,217 trees and 57 species at 4,375 sites distributed across Canada. We found that annual radial growth responses to FFD varied by species, with some commonalities across genera and clades. The growth of gymnosperms with late spring leaf-out strategies was negatively related to FFD; years with high FFD were most detrimental to the annual growth of Pinus banksiana, Pinus contorta, Larix lyalli, Abies amabilis, and Abies lasiocarpa. In contrast, the growth of angiosperms with early leaf-out strategies, namely, Populus tremuloides and Betula papyrifera, was better in the coldest years, and gymnosperms with intermediate leaf-out timing, such as widespread Picea mariana and Picea glauca, had no consistent relationship to FFD. Tree growth responses to FFD were further modulated by tree size, tree age, regional climate (i.e., mean cold-season temperature), and local site conditions. Overall, our results suggest that moderately warming winters may temporarily improve the growth of widespread pines and some high-elevation conifers in western Canada, whereas warming winters may be detrimental to the growth of widespread boreal angiosperms. Our findings also highlight the value of using species-specific climate-growth relationships to refine predictions of forest carbon dynamics.

The boreal and temperate conifer-dominated forests that cover most of Canada play a significant role in global carbon capture and storage (1). These northern biomes, particularly boreal forests, are experiencing some of the world’s most rapid increases in temperature, with disproportionate warming occurring in winter (2). In parallel with mean temperature increases of 1.7 °C, the freeze-free season has lengthened by more than 20 d, on average, between 1948 and 2016, with the season starting earlier by about 9 d and ending later by about 12 d (3). Overall, the number of freeze days has decreased by an average of 15 d, with the greatest reduction found in western Canada south of 60°N (3). Since trees are often highly attuned to climate, these changes could impact tree growth and survival across vast areas. While much recent research suggests tree growth and carbon sequestration may be increasingly limited by summer drought (4), it is yet unclear how winter warming, and changes to seasonal freeze trends, might affect these ecological processes.

Tree species of the northern hemisphere have evolved strategies to cope with freezing temperatures. In late summer and early autumn, decreasing photoperiod (day length) and cooling temperatures trigger the cessation of both shoot and stem growth and the setting of buds and other physiological processes that enable resistance to freeze damage (5, 6). This cold acclimation period is followed by dormancy in winter, when trees reach maximum cold-hardiness (resistance to freeze). Trees require a period of chilling—daily mean temperatures between −2 °C and 10 °C—to break winter dormancy and resume growth in response to increasing day length and warming with the onset of spring (7), at which point tree resistance to freezing temperatures declines (5, 6). The timing of this growth and dormancy cycle is species specific, “hard-wired” by genetics, and modulated by tree vigor and local environmental conditions (8). As the global climate warms, and becomes more variable, temperature cues may become increasingly mismatched with these hard-wired phenological cycles, making trees more vulnerable to freezing conditions (8, 9).

Trees exposed to severe cold that exceeds their freeze resistance can be damaged in multiple ways, with consequences for tree growth. Unseasonably warm spells, followed by below-freezing temperatures, can cause injuries to tissues during midwinter and early spring (9–12). These injuries may include damage to the apical and cambial meristems, photosynthetic apparatus, and roots. For example, photoinhibition caused by the stress of freezing temperatures on the photosynthetic apparatus suppresses photosynthetic capacity and quantum-use efficiency in spruce, with recovery taking up to 13 d after the freeze (13). Unusual early-freeze events in autumn may cause leaf abscission before leaves are fully senescent; this may reduce opportunities for remobilization of leaf nutrients and be detrimental to carbon and nutrient storage (8). Xylem cavitation, which compromises water transport to foliage in the following growing season, may also increase during harsh winters through freeze-thaw embolisms (8, 14). Furthermore, although tree roots may tolerate short-term low temperatures in their freeze-hardy state (15), mild but long-lasting soil freezing can damage fine roots in saplings or mature trees (16–20). Because snow cover insulates and prevents soil freezing in winter and buffers against cold snaps in autumn and spring, root injuries are most severe when snow cover is low (21–23). In the same way the legacy effects of summer water stress can cause episodes of mass tree mortality, freezing can lead to growth loss, permanently damaging trees and potentially resulting in death several years later (13, 24–28). While the most extreme freezes can damage many species, and trees of all sizes, over vast areas, freeze damage often varies both among and within species (6, 11, 29, 30).

There is a strong need for real-world, in situ, multispecies comparisons to quantify the differential impacts of climate change on northern forests. Although significant research has been conducted on the impacts of freeze damage on growth and carbon sequestration potential (13, 23, 31–34), scaling research results across northern forests can be challenging. Practical constraints to research sometimes lead to limited site replication and short time spans. To quantify freeze impacts, often controlled conditions in the field or laboratory are required, but a large number of observations are needed to complement these experimental findings. These constraints therefore leave considerable gaps in our understanding of cold acclimation (6, 35, 36). Consequently, freeze impacts are not included in many vegetation models or are limited to acting primarily on photoinhibition processes (13) rather than incorporating tree damage and impacts on growth (37, 38). Quantifying the response of tree growth to freeze could improve our understanding of forest growth responses under climate change, which will benefit model predictions that form the basis of climate change mitigation and adaptation strategies and sustainable forest management.

Here, we provide a comprehensive assessment of the response of forests to the frequency of freeze days (hereafter, FFD), based on a well-replicated growth increment network covering Canada’s forests over the period from 1951 to 2018. Our dataset represents extensive sampling from boreal, temperate conifer, and temperate deciduous forests. We define FFD as the percentage of days with a minimum temperature below freezing (0 °C) from November to April. Using mixed-effects models, we quantified the relationship between tree radial growth increment and FFD, as well the relative predictive power of other freeze variables. We then explored mechanistic reasons for variability in the effect of FFD on growth, including factors such as species, leaf-out strategy, regional climate and local site conditions, and tree age and size. Finally, we determined how FFD and growth are changing over time for the most responsive and widespread species.

Results

Seasonal and Spatial Patterns in the FFD.

We summarized monthly and spatial patterns in the cold-season FFD (Fig. 1). Across Canada, monthly FFD begins to increase in September and plateaus from December to March, as 90 to 100% of days fall below the freezing point during that period. FFD then decreases in May. The period from November to April, when mean daytime temperatures (Tday) are generally near or just below freezing, is hereafter referred to as the cold season. Between 1951 and 2018, cold season FFD varied from <80 to 100% across Canada. Forests in coastal regions, and north of the Great Lakes, had the lowest percentage of FFD. FFD typically averaged ∼100% in the forests of mountainous regions in western Canada and over most of the boreal forest (north of about 50 to 55°N), with the exception of central Canada’s boreal area where FFD was lower (∼90 to 94%). Cold-season FFD decreased exponentially with rising Tday (), dropping below 50% when Tday exceeds 4 °C, and then decreasing by ∼8% for each 1 °C temperature rise (Fig. 1).

Fig. 1.

Climate characteristics of the study area. (A) Long-term average (1951 to 2018) cold-season FFD (% of all days below zero °C from November to April). Circles indicate sites from which we obtained tree-ring data (n = 4,375). (B) Monthly mean Tday (brown) and monthly FFD (blue) from July of the year previous to growth (lowercase text) to June of the year contemporaneous to tree growth for the period 1951 to 2018. The boundaries of the box denote the 25th and 75th percentiles, and the line within the box marks the median of the distributions. Error bars above and below the box indicate the 90th and 10th percentiles. Gray shading delineates the cold season. (C) Cold-season FFD versus mean cold-season Tday.

Relationships between FFD and Basal Area Increment (BAI).

Using tree-ring data from 35,217 trees (57 species) across Canada, we constructed 5,911 species-by-site Generalized Additive Mixed Models (GAMMs) to assess the sensitivity of annual BAI to FFD (Fig. 1 and ). The sites from which we obtained tree-ring data spanned the range of variability in FFD across Canadian forests. Convergence was successful for 78% of the species-by-site GAMM models. Among these species-by-site models, 35% had at least one tree where the t-value for the relationship between annual BAI and cold-season FFD reached a threshold of P < 0.05. The t-value indicates the strength of the statistical associations, where negative t-values mean that high FFD is detrimental to tree growth and positive t-values indicate high FFD favors tree growth. Overall, a significant relationship to FFD was reached in 11,218 tree samples. About two-thirds of these statistically significant t-values were negative, indicating that years with high cold-season FFD were often detrimental to tree growth (). Negative t-values were concentrated in western Canada and north of the Great Lakes (Fig. 2 and ). The remaining one-third of the statistically significant t-values (positive) suggest that high FFD is sometimes also associated with increases in annual growth.

Fig. 2.

Tree species growth responses to cold-season FFD. (A) Strength of relationships between tree BAI (mm2) and the FFD during the cold-season (November to April) for the period (1951 to 2018), expressed as t-values obtained from GAMMs. Negative t-values (red) indicate areas where high FFD is detrimental to tree growth, t-values close to zero denote an absence of relation between BAI and FFD, and positive t-values (blue) indicate areas where high FFD favors tree growth. The map was obtained by applying an inverse distance weighted interpolation to the vector of georeferenced t-values (n = 28,115 trees). Sampling for interpolation was set at n = 100 trees with a search radius of 10 degrees and resolution of 0.1 degree (ArcGIS v10.5.1). Circles denote sites where a significant statistical association was detected for at least one tree. (B and C) Same as above but for the combined responses of two evergreen gymnosperm species (P. banksiana and P. contorta) and two deciduous angiosperm species (P. tremuloides and B. papyrifera) associated with the most strongly negative and positive growth responses to FFD, respectively. (D) Distribution of t-values for all 57 species (gray) and for species with large sample sizes (i.e., n > 250 trees, colored). (E) Number of tree-ring series per species. (F) The percentage of significant t-values (P < 0.05, positive and negative) by species ().

GAMMs constructed using cold-season FFD outperformed models using 16 other freeze metrics to predict BAI ( and Table S3), indicating the importance of using FFD to predict cold exposure risk. In addition, winter FFD (December to February) yielded a high percentage (30%) of significant relationships compared to FFD in the previous fall (September to November; 17%) and spring (March to May; 22%) ( and Table S3). Also, day-of-year when accumulated growing degrees reached 100 °C (typically during mid-May) reproduced the FFD distribution of t-values well, albeit with a much lessened percentage of significant t-values (). In general, indicators meant to quantify freezing after growth initiation (e.g., the mean temperature or sum of growing degree days (GDD) prior to the last freeze day in spring) did not perform as well as FFD ( and Figs. S3 and S4).

Growth Responses to FFD by Clade and Species.

The effect of FFD (expressed as t-value) on tree BAI varied by clade (gymnosperm vs. angiosperm) and by species. Among gymnosperms, 30% of trees had significantly negative relationships between BAI and FFD (); significantly positive relationships were less common (∼11%). Angiosperms exhibited the opposite trend, with BAI for 23% of trees having a significant positive relation to FFD (this relationship with BAI was negative in 9% of trees). Among species, FFD had the strongest effect (|t-values| <2) on the BAI of well-replicated and widely distributed gymnosperm species, namely, Pinus banksiana and Pinus contorta, as well as Larix lyalli, Abies amabilis, and Abies lasiocarpa (Fig. 2 and ). For all of these species, more than 50% of trees had a high percentage of significant negative relationships with FFD (Fig. 2 and and Table S2).

A relatively high percentage of trees (30 to 50%) having negative relationships between BAI and cold-season FFD was also observed for gymnosperm Pseudotsuga menziesii, Picea engelmannii, and Thuja occidentalis (Fig. 2 and and Table S2). Among several gymnosperm species with large geographic ranges—Picea glauca, Picea mariana, Abies balsamea, and Larix laricina—the effect of FFD on BAI was considerably less and the percentage of significant relationships was more evenly partitioned among positive and negative t-values (Fig. 2). Negative BAI responses to FFD in P. mariana and P. glauca were of higher magnitude north of the Great Lakes ().

Growth of well-sampled angiosperm species with the largest geographic ranges, namely, Populus tremuloides and Betula papyrifera, responded mostly positively to high cold-season FFD (Fig. 2 and and Table S2). However, the strength of these positive relationships (t-value and percentage of trees with significant t-values) was inferior to the negative relationships shown by the boreal pine species (Fig. 2). Of note was that for 20 to 23% of P. tremuloides and B. papyrifera samples, growth was improved when the day of year when cumulative GDD of 100 °C occurred later, supporting the idea that they have lower growth following shorter winters and when the last freeze date occurred earlier (), implying possible injury from spring frost. For other angiosperms, our ability to provide species-specific summaries for the effect of FFD on BAI is tenuous because of limited site replications and sample sizes (). Grouping of samples by genera provided indications that samples from the Acer and Quercus genera had small percentages (<25%) of significant BAI and FFD relationships. In contrast, the Fraxinus genus had a high percentage (86%) of significantly negative BAI and FFD relationships ( and Table S2).

Modulating Influences of Tree Traits, Site Characteristics, and Site Climate.

Our Generalized Additive Models (GAMs) indicated that variability in BAI response to cold-season FFD (e.g., t-values from above) was significantly influenced by species-specific spring leaf-out timing, tree development stage (tree age and size), and climatological means of Tday and total snowfall at the tree-ring sites (Fig. 3 and ). The effect of FFD on BAI was positive among species having early leaf out and negative among species leafing out later in the spring (after day of year 125) (Fig. 3 and ). As indicated by higher F-values, tree size had a greater effect on BAI sensitivity to FFD than tree age (Fig. 3). Predicted t-values increased with increasing tree basal area until 4.0e+4 cm2 and leveled-off afterward. Predicted t-values increased with age among trees older than 200 y. FFD effects on BAI were strongly positive among eastern sites (east of –80°W) and strongly negative among northwestern sites (west of –110°W and north of 50°N; Fig. 3). We also predicted a negative relationship between BAI response to cold-season FFD and cold-season mean Tday (Fig. 3), indicating that FFD had the most detrimental influence on annual BAI at warmer sites (i.e., mean Tday greater than 10 °C over the period 1951–2018). At colder sites, higher-than-average cold-season FFD was beneficial to annual growth. Moreover, at sites where mean cold season snowfall was >200 mm, the response of BAI to cold-season FFD was negative and became increasingly stronger in regions with increasing total snowfall (Fig. 3). In addition to variation in FFD effects along climate gradients, and due to tree phenology and ontogeny, we found that geographic location and random site effects were also very important predictors of t-values (Fig. 3 and and Table S4).

Fig. 3.

Variation in predicted growth response to cold-season (November to April) FFD (i.e., t-value) relative to the following explanatory variables: mean climate conditions (A and B), joint effects of latitude and longitude (C), and tree traits (age, size, leaf-out timing [D–F, respectively]). A GAM of relationships between t-values (see Fig. 2) and the aforementioned explanatory variables generated the predictions (n = 28,115 trees). Predicted freeze response t-values from the GAM are shown as black lines with two-sided 95% confidence intervals shown as blue shading for all panels except C, which illustrates a multivariate t-value surface without confidence intervals. Tree age is cambial age. Leaf-out day was inferred from an experiment conducted in 2012 at eight common garden studies in the northern hemisphere (9); for graphing, the mean leaf-out day is delineated for those genera having more than 10 sampling sites (arbitrarily defined). The F-test and associated P values are shown for each factor in their respective panels. For reference, a term for random site effect was included in the GAM and produced statistically significant results (F = 16.8 and P < 0.001, ).

Changes in FFD and Growth of Widespread Pines.

We examined changes in the annual growth of two widespread pine species that responded strongly and negatively to annual variability in FFD, namely, P. banksiana and P. contorta (see Fig. 2). Both species occur in areas that have undergone important declines in FFD since 1951 (Fig. 4). We found that percent growth changes (GC) in terms of BAIs, averaged by species, were linearly and inversely related to FFD (Fig. 4). From the slope of this relationship, we estimated that for about every percent increase in cold-season FFD, there was an ∼3% decrease in mean BAI throughout the geographic ranges of these pines (Fig. 4). This response is substantial, considering the average FFD across species’ range varies annually by about 8% (i.e., 14 d). Increases in tree growth from 1951 to present were apparent in both species, paralleling with decreases of cold-season FFD over the same period (Fig. 4).

Fig. 4.

Changes in the annual growth of P. contorta and P. banksiana relative to cold-season (November to April) FFD (%). (A) Distribution of P. contorta and P. banksiana sample sites on a map of temporal trends in FFD over the period 1951 to 2018. Contour lines delineate the 95% significance level. (B) FFD (blue line) and percent GC (i.e., change in BAI, red line) averaged across all sites; the shaded band represents the 95% confidence interval. Annual GC are expressed as the percent deviation from predicted values generated by the GAMMs of relationships between BAI and the effects of tree size and age, as well as growth legacies effects been accounted for (Material and Methods, Eq. ). Only years in which the number of sampled sites exceeded 100 were used. (C) Correlative relationships (r) between GC and FFD (years where sites >100); the 95% confidence intervals around r were computed from a bootstrapping technique that accounted for autocorrelation and trend in data (78), and regression slope (β1) from ordinary-least-square regressions.

Discussion

Our cross-continental assessment shows that tree growth in Canada’s forests is significantly influenced by exposure to freezing. Using annual growth rings from over 35,000 trees and 52 species across Canada, the FFD during the cold season (November to April) was the most powerful predictor of annual BAI among the multiple variables tested. This suggests that freezing conditions have been a pervasive driver of interannual variability in tree growth across Canada’s forests over the last seven decades. While some studies have examined forest growth patterns across large areas, like the boreal zone (4, 39), we have comprehensively linked yearly fluctuations in freeze exposure over time to variations in species’ annual growth. The patterns of variation in the direction and magnitude of the relationships between growth and freeze exposure are important findings of our study. A large part of this variation can be attributed to clade, but we also observed large variations in growth sensitivity to FFD that were species specific and related to spring leaf-out timing, tree size, and tree age, as well as local climate and site factors.

Tree growth over the last 70 y was strongly related to FFD over vast areas of Canada. However, the direction and strength of the relationship between annual BAI and FFD varied substantially among clades and species. The annual growth of widely distributed gymnosperms P. banksiana and P. contorta and the growth of widely distributed angiosperms P. tremuloides and B. papyrifera were most sensitive to yearly fluctuations in FFD, with gymnosperm growth responding negatively to higher than average FFD and angiosperm growth responding more positively. The contrasting growth responses of these angiosperms and gymnosperms and high sensitivity of pines to FFD are consistent with other continental- to global-scale assessments of the effect of temperature changes (and drought) on tree growth (40, 41). Also, our findings that years with above-average FFD had especially adverse effects on P. contorta and P. banksiana growth are similar to other studies of freeze effects on pine in Canada and Europe (25, 26, 42). Furthermore, P. menziesii, an important timber species in western North America, is notoriously sensitive to freezing conditions (29). This was reflected in our observations of its moderately negative growth responses to high FFD as well. Boreal species P. glauca and P. mariana, as well as high-latitude species L. laricina and, to some extent, A. balsamea, were less sensitive to yearly fluctuations in FFD. Having evolved strategies to tolerate freezing temperatures, tree species of boreal forests and in high-latitude environments are highly adapted to freezing temperatures (6). Therefore, it is not surprising that the annual growth increment of these species was less sensitive to yearly fluctuations in FFD than the growth of some lower latitude pines, other more southern gymnosperms, or the angiosperms we studied. Overall, our analyses using tree-rings indicated diverse growth responses to FFD that broadly align with long-recognized trends in the relative cold tolerance of boreal and temperate tree species (43).

The positive responses of angiosperm growth to higher-than-average yearly FFD but negative growth responses among evergreen gymnosperms are likely due to a continuum of physiological, phenological, and anatomical trait differences. The growth of two widespread pine species, namely, P. contorta and P. banksiana, was more vulnerable to higher-than-average FFD than most other gymnosperms. Extracellular freezing of water, and associated cell desiccation during winter, commonly causes pine needle damage that is expressed by reduced growth among trees (12, 44) and may partly explain our finding. Also, in midseral pine stands where stem density is very high, colder-than-average winters, combined with strong winds, force tree collisions that can cause substantial mechanical damage whereby twig/branch breakage and needle loss adversely impact growth potential in the months to follow (45).

Although gymnosperms have evolved many physiological mechanisms to mitigate the impacts of freezing, foliar damage, and subsequently reduced growth, can still occur in extremely cold years. The deciduous habit of the angiosperms we studied may confer an advantage. Our finding that the deciduous conifer L. laricina also had better growth during years of high FFD supports this hypothesis. In general, our findings of contrasting growth responses to FFD between angiosperms and gymnosperms support the view that angiosperm traits permit more plastic responses to extreme cold or drought (40, 46). Differences in spring leaf-out timing may play a role in the observed FFD response differentiation among gymnosperms and angiosperms (41). Earlier leaf-out, which marks the beginning of significant photosynthesis in angiosperms species (47), as well as higher rates of carbon transfer (41), may confer advantages to Populus and Betula. However, early spring leaf-out also comes with high risks of foliar damage from spring freezes (48); our analyses showed negative impacts of a late last freeze day in ∼22% of samples for both these species (cf, and Table S3). In contrast to Populus and Betula, the growth of trees with late-spring budburst (particularly Pinus spp.) responded negatively during years of higher-than-average cold-season FFD. Differences in growth responses to FFD among Populus, Betula, and Pinus likely reflect the relative capacity of angiosperms and gymnosperms to recover and allocate carbon to growth increment following long, harsh winters. However, the relatively high percentage of angiosperm trees that responded positively to FFD could also indicate adverse effects of lower-than-average cold-season FFD on bud development; without adequate chilling, budburst can be delayed and shoot growth reduced (49, 50). This could be particularly important in central Canada, where Populus and Betula have expressed the greatest positive sensitivity to FFD (Fig. 2). It is also possible that positive responses seen in BAI of angiosperms to FFD reflect an energy deficit (increased respiratory losses) being accrued by these species due to prolonged warm springs (50, 51). Leaf out might not be advancing at the same rate as warming (52), implying that a depletion of carbon reserves may be caused by a metabolism being activated earlier than photosynthesis is renewed.

Tree age and size were also important predictors of variations in the direction and magnitude of growth responses to FFD. The predicted annual growth increment of large (and old) trees responded more positively to high FFD than the growth of small trees. Global studies of critical temperatures for xylogenesis show that the onset of cambial activity can vary among trees of different sizes and ages (53). The insulating effect of thicker bark on large, old trees protects the cambium and shifts critical temperatures for xylogenesis, permitting large old trees to recover growth increment earlier than small trees (54). These differences may be particularly apparent during very cold years when the stems of young trees may also be more susceptible to tissue damage (55).

A high level of uncertainty remained in the attribution of the freeze responses in growth to environmental factors. This was indicated by the strong contribution of geographic location and random site effects to explaining variability in BAI responses to FFD. Strong negative effects of FFD on BAI in southcentral and southwestern Canada and strong positive effects in northwestern Canada—as well as relatively insensitive growth in other areas of the country—imply regional ecotypic differences in freeze tolerance and acclimation to fluctuating conditions (56). However, broad regional climate differences (total snowfall and mean daily temperature) accounted for a relatively small amount of variability in growth sensitivity to FFD compared to geographic location. Several other factors are likely influential, including microtopography (cold air drainages) (57), microclimate (e.g., exposure to wind and cold (45), soil nutrient and moisture conditions, and stand composition and structure. Within tree species distributions, populations can also express different responses due to genetic adaptation to local environments (29, 30, 58). Such population structure and genetic adaptation across landscapes could affect our results.

Canada’s forests are important contributors to the global carbon balance. As such, there is a growing interest in describing long-term changes in tree growth and forest productivity in the context of climate change. Across Canada, FFD decreased over the past 70 y (3), and likewise, the cold-season freeze trends affected tree growth and carbon sequestration potential across a variety of species in diverse ways over that period. Many studies have focused on quantifying forest growth and productivity trends in Canada using multiple data sources, including remote sensing and tree rings, but report very different results (59). Given our findings that tree growth responded diversely to variability in FFD, inconsistencies in reports of temporal changes in forest growth/productivity among these studies is not surprising. The high spatial variability of growing conditions found in Canada’s forests, as well as the multiple interactive effects and feedbacks of the environmental variables present, contribute to a significant diversity in the response of trees to climate. This being said, long-term impacts of freezing dynamics on growth can be exemplified from the results of this study for the two evergreen gymnosperms P. banksiana and P. contorta. These two species stand out particularly strong in their sensitivity to FFD, and they both have significant ecological implications owing their wide distribution and contrasting climatic environment. The effect of FFD is of broad significance in these two species, and increases in growth in parallel with decreases of cold-season FFD during the last 70 y are observable. Our findings suggest positive outcomes for tree health and growth of these pine species as the climate warms and damage associated with winter freeze decreases. This is true, however, to the extent that other growth-detrimental factors, such as drought and insect disturbances, do not counteract the benefits.

Conclusion

Our subcontinental analyses of tree-ring data indicate highly variable tree growth responses to yearly fluctuations in cold-season freeze frequency across Canada and among species. Leaf-out timing and tree ontogeny predicted tree annual increment in responses to yearly fluctuations of FFD better than regional climate. Projected increases in winter temperatures imply overall decreases of cold-season FFD in the 21st century in many high latitude and altitude regions (36, 60). Mechanistic representations of the interaction between climate and physiological processes in forests are desirable to improve the precision of simulations describing the implications of these changes on carbon dynamics (27). In that regard, our results demonstrate that the use of a metric based on the frequency of extreme conditions performs better than average or cumulative climate metrics, such as mean temperature and degree days, to predict tree growth. Although our results may suggest that diminishing freeze exposure may temporarily increase the BAI of gymnosperms with late-spring leaf-out strategies, other seasonal conditions like summer aridity and damage by insects and pathogens may become increasingly growth limiting in the future. Additionally, warming winters could negatively impact angiosperms species in need of sufficient chilling or through the negative impacts of more frequent winter thaw events (61). These considerations have merit in a context where stress factors compound each other, affecting the health of trees and forests (62). Our empirical assessment of tree growth response to cold-season freeze in boreal forests can support efforts to represent tree growth processes in global vegetation models, which is urgently needed for improved nature-based solutions to mitigate, and adapt to, climate change (4, 63).

Materials and Methods

Tree-Ring Data.

Tree-ring width data were extracted from CFS-TRenD (64), a data repository containing increment core samples obtained during forest inventory surveys and targeted collections from old trees on more ecologically marginal sites. The data in the repository have already undergone quality-control procedures (64). A significant proportion of the data in the repository come from tree cores sampled during the establishment of Canada’s National Forest Inventory (NFI) program (), while other important data contributions are more regional in scope (). These tree-ring data also contain several contributions from smaller-scale research projects, including tree-ring datasets from the International Tree-Ring Data Bank (65). The dataset we compiled from all these sources contained annual ring-width measurements from 35,217 trees (greater than 30 y old to ensure sufficient degrees of freedom in our analyses), collected at 4,594 sites across Canada. The data represent 57 tree species and all of Canada’s forested biomes, with the largest being the pan-Canadian boreal forest.

Climate Data.

As per the Glossary of Meteorology of the American Meteorology Society (66), we characterized freeze as the condition that exists when, over a widespread area, the surface temperature of the air remains below freezing (0 °C) for long periods. Accordingly, we computed the annual FFD at each site as the percentage of days with a minimum temperature of <0 °C from November of the year previous to April of the year contemporaneous with tree growth (i.e., the cold season). BioSIM (67) provided daily minimum temperatures (Tmin, °C) for the period 1950 to 2018. This software tool generates site-specific estimates of Tmin (and other air temperature and precipitation variables) from historical weather observations by interpolating daily data from the four closest weather stations, adjusting them for elevation and location differentials with regional gradients, and averaging values using a 1/d2 weight, where d is distance. In addition to cold-season FFD, we also calculated monthly and seasonal (spring, winter, autumn) FFD for each year between 1951 and 2018.

Through graphical analysis, we depicted spatial and temporal patterns in FFD. We also explored relationships between FFD and other characteristic climate features for the cold season (November to April), as follows: Tday (0.75*Tmax+0.25*Tmin), cold-season snowfall (mm of water), cold degree days (i.e., cold-season sum of daily mean temperature below 0 °C), GDD (i.e., cold-season sums of daily mean temperature above 0 °C), and the accumulating warming before the last spring freezing (i.e., the sum of GDD from January first until the last freeze day [9]). We considered Tday, rather than mean daily temperature, because it better reflects the temperature experienced during photosynthesis (i.e., weighted by Tmax) (68).

Statistical Analyses.

We used GAMMs (69, 70) to explore the effect of FFD on tree growth. To do this, we modeled BAI (cm2/yr−1) as a function of tree ontogeny (cambial age and tree size; i.e., basal area [BA], cm2), and the explanatory climate variable FFD. Prior to analyses, annual ring-width measurement series were scaled into annual estimates of BAI. A logarithmic-power transformation (log) was applied to BAI and BA because of the skewed distribution of their values. We then fitted species by site GAMM that took on the form:where h represents the site, k represents the tree, t represents the year, and the ε residual error. Cambial age refers to an estimation of tree age based on the ring counts (age i). We explored the tree-specific FFD effect by setting tree identity as categorical predictor in the model (). We considered tree identity (Z) as a random effect. We also included an error term () with a first-order autoregressive (AR1 [P = 1, q = 0] correlation structure, which is usually sufficient to account for legacy effects on current year growth increment (71). The smoothing term, s, represents a cubic regression spline for which the degree of smoothness was determined through an iterative fitting process (69). We quantified tree growth sensitivity to FFD using t-value statistics (i.e., model parameter estimate divided by its corresponding unconditional SE) and their associated P values. We quantified the importance of FFD effect by the percentage of significant t-values (P < 0.05, positives and negatives) by clades, genera, and species. GAMMs were also fitted to 16 other freeze-related thermal metrics to predict BAI. These included among others calendar-season FFD, cold-season snowfall, cold degree-days, growing degree-days, fifth percentile of coldest days, day-of-year of last freeze event, and the accumulating warming before the last spring freezing. The performances of these metrics were also measured accordingly with their percentage of significant t-values. The GAMMs were fitted using the mgcv package version 1.8-4 (70) in R (72).

To gain further insight into the potential mechanisms driving growth responses to freezing, we modeled FFD effects on the growth of each tree (quantified by t-value statistics above) in relation to mean cold-season climate at each site, tree ontogeny, and a leaf-out trait. Mean cold-season climate was characterized by two variables—mean cold-season Tday (°C) and total amount of snowfall (mm of water)—which were calculated annually from 1951 to 2018 and then averaged across all years. Variables characterizing tree ontogeny included tree age and BA. The leaf-out trait, i.e., the average day of the year budburst is initiated, was determined for each species from observational data of leaf-out timing in spring 2012 at eight common garden studies in the northern hemisphere (9). Despite this single year of vegetation phenology data, this experimental setup minimizes the influence of phenotypic differences resulting from environmental forcing in native habitats and genotypic variation in leaf-out dates. Leaf-out data were not available for all the species we studied; for those without data, we used mean leaf-out values calculated at the genus level. We employed a forward-stepwise procedure to select the most influential predictor variables to include in the GAM. A variable selected at each step satisfied three criteria, as follows: 1) P values of all the variables in the model are less than 0.05, 2) they had the smallest Akaike information criterion (AIC) in each step, and 3) AIC decreased more than it had in the previous step. Collinearity among the candidate variables was assessed using the variance inflation factor (VIF), with the maximum VIF limit set at 10 at each step. We fit GAMs (73) that took the following form:where the effect of FFD on BAI (expressed as a t-statistic) at site h, tree k, and species p was cast as a function of clim, all the climate variables (i = 1 to N) at site h; ontogeny, all the ontogeny variables (i = 1 to N) at site h, for tree k: and, leaf outp, a genus-specific (p) day to budburst. The term spatial explicitly models the joint effects of latitude and longitude, and site identity, Z, was accounted for the site random effects. We used penalized cubic regression splines as smooth functions for the climate, ontogeny, and leaf-out trait variables and splines on the spheres for the spatial effects. Because the dataset was very large, we used the BAM function for memory-efficient estimation and the associated fast REML (fREML) method to fit the model (73, 74). Residuals were tested for the presence of spatial dependencies using Moran’s I (spdep R package v1.1-5 (75)). Goodness-of-fit of smoothed terms were assessed using an F-test and associated P values. The GAM was fitted using the mgcv package v1.8-33 (70) in R v4.0.4 (72). We reported model results graphically using the SigmaPlot v14 and ArcGIS v10.5.1 for Windows.

Linear trends of FFD were examined for 1951 to 2018 using least-squares linear regressions. The significance of the slope was tested against the null hypothesis that the trend is different from zero, by using a variant of the t test with an estimate of the effective sample size that takes into account the presence of serial persistence in data (76). Results were reported graphically using SigmaPlot v14 and ArcGIS v10.5.1 for Windows.

To visually illustrate the relationship between freeze frequency and annual fluctuations in growth over time, we removed tree age and size trends from BAI using the same GAMM approach described by Eq. but eliminated climate variables from the model. The fitted species by site GAMM hence took the form:

Annual GC, expressed as the percent deviation from predicted values generated by the GAMM, were then computed following Girardin et al. (39). Next, for each calendar year, we calculated the mean values of all samples, by species, to show temporal variability in species-specific GC since 1951. For uncertainty in GC estimations, calculated as species-specific 95% confidence intervals for each calendar year, we first corrected the effective degrees of freedom (n′) based on first-order (i.e., lag = 1) autocorrelation estimates of Moran’s I (77) (moran.test function). Then we limited our analysis to those years during which n′ >100, considering the temporal coverage of tree-ring data and allowing us to focus on the period during which the uncertainty band was relatively stationary. Correlative relationships (r) between GC and FFD were computed from a bootstrapping technique that accounted for autocorrelation and trend in data (78) and regression slope (β1) from ordinary least-squares regressions.