Impacts of Climatic Variation on the Growth of Black Spruce Across the Forest-Tundra Ecotone: Positive Effects of Warm Growing Seasons and Heat Waves Are Offset by Late Spring Frosts

Climate strongly limits the physiological processes of trees near their range limits, leading to increased growth sensitivity. Northeastern North America is experiencing considerable warming, so the growth of trees near the northern treeline represents a key indicator of forest responses to climate change. However, tree-ring series and corresponding climatic data are scarce across the forest-tundra ecotone when compared to southern boreal regions, resulting in fewer studies on growth-climate relationships focused on this ecotone. Using daily climatic data, we identified trends in growing season heat accumulation and the intensity of acute climatic events over the last several decades in the southern and the northern parts of the forest-tundra ecotone in northeastern North America, and investigated their influence on black spruce radial growth. We found that black spruce trees responded positively to the increase in growing season temperatures and heat wave intensity, suggesting that growth is currently limited by suboptimal temperatures. While tree growth in the southern region generally benefited from warm spring temperatures, vulnerability to late spring frosts reduced tree growth in the northern region and increased probability of abrupt growth decline. In this region, late spring frosts offset approximately half of the additional growth that would otherwise occur over the course of a warm growing season. This vulnerability of northern trees may result from local adaptations to short growing seasons, which initiate biological activities at colder temperatures in the spring. Overall, our results highlight the need to explicitly incorporate acute climatic events into modeling efforts in order to refine our understanding of the impact of climate change on forest dynamics.

Climate strongly limits the physiological processes of trees near their range limits, leading to increased growth sensitivity. Northeastern North America is experiencing considerable warming, so the growth of trees near the northern treeline represents a key indicator of forest responses to climate change. However, tree-ring series and corresponding climatic data are scarce across the forest-tundra ecotone when compared to southern boreal regions, resulting in fewer studies on growth-climate relationships focused on this ecotone. Using daily climatic data, we identified trends in growing season heat accumulation and the intensity of acute climatic events over the last several decades in the southern and the northern parts of the forest-tundra ecotone in northeastern North America, and investigated their influence on black spruce radial growth. We found that black spruce trees responded positively to the increase in growing season temperatures and heat wave intensity, suggesting that growth is currently limited by suboptimal temperatures. While tree growth in the southern region generally benefited from warm spring temperatures, vulnerability to late spring frosts reduced tree growth in the northern region and increased probability of abrupt growth decline. In this region, late spring frosts offset approximately half of the additional growth that would otherwise occur over the course of a warm growing season. This vulnerability of northern trees may result from local adaptations to short growing seasons, which initiate biological activities at colder temperatures in the spring. Overall, our results highlight the need to explicitly incorporate acute climatic events into modeling efforts in order to refine our understanding of the impact of climate change on forest dynamics.

INTRODUCTION
As the rate of atmospheric warming increases at high latitudes of the northern hemisphere (ACIA, 2004;IPCC, 2014), arctic and subarctic ecosystems are undergoing major changes (ACIA, 2004). Over recent decades, a general increase in primary productivity has been observed at the circumpolar scale, where a significant expansion of shrub species has been linked to warmer temperatures, longer growing seasons, and increased precipitation (Myers-Smith et al., 2011;Elmendorf et al., 2012;Tremblay et al., 2012;Fraser et al., 2014;Moffat et al., 2016). However, the effect of climate warming on trees growing at high latitudes remains unclear, as revealed by the contrasting responses of treelines (Payette, 2007;Dufour-Tremblay et al., 2012) and tree growth (Vallée and Payette, 2004;D'Arrigo et al., 2008;Andreu-Hayles et al., 2011;Dufour-Tremblay et al., 2012). Considering that around 30% of global terrestrial carbon is contained in boreal forests (Mcguire et al., 2009), a refined understanding of their response to climate change is essential, especially as it may have significant implications for the regulation of future climate (Bonan, 2008).
Tree physiological processes are strongly limited by climatic conditions near their range limits, resulting in an increased climatic sensitivity near the treeline (Holtmeier and Broll, 2005). Accordingly, the growth response of trees to a changing climate is considered a key indicator for anticipating how forest structure and function may be altered by climate change (Payette et al., 1989;Callaghan et al., 2002). However, because of their remote location, tree-ring series and corresponding climatic data are scarce in northern regions, so there are fewer studies on growthclimate relationships in the forest-tundra ecotone compared with the southern boreal forest (i.e., Goldblum and Rigg, 2005;Huang et al., 2010;Girard et al., 2011;D'Orangeville et al., 2016D'Orangeville et al., , 2018Ols et al., 2018;Marchand et al., 2019;Marquis et al., 2020). This is particularly true for northeastern North America, where very few studies have investigated the relationships between tree growth and climatic conditions at such high latitudes (but see Dufour-Tremblay et al., 2012;Nicault et al., 2015). Near the northern treeline, Dufour-Tremblay et al. (2012) reported nonsignificant relationships between growth of black spruce (Picea mariana Mill. B.S.P.), the dominant forest-tundra tree species, and both monthly temperature and precipitation. This apparent lack of climate-growth relationships is surprising considering the strong warming trends observed over the last decades in northeastern North America , and could be due to the low temporal resolution (i.e., monthly) of the climate data used in the analyses, which may obscure the effects of more discrete climatic events (Graumlich, 1993;Ols et al., 2018;Moreau et al., 2020). Important shifts in forest composition, productivity and mortality regimes have been observed worldwide as a result of the increasing occurrence and intensity of acute climatic events (Allen et al., 2010;Carnicer et al., 2011;Anderegg et al., 2015;Dorothea et al., 2015;Zhang et al., 2018). In temperate forests, late spring frosts, defined as freezing that occurs after a substantial accumulation of warmth during spring, cause injuries that influence the range limits of tree species (Chamberlain et al., 2019;Zohner et al., 2020). Severe heat waves can also induce important physiological stress through an increase of the vapor pressure deficit, which in turn is strongly related to droughtinduced tree growth decline and mortality (Grossiord et al., 2020). Finally, severe cold waves have been identified as a factor potentially responsible for the progressive loss of foliage and mortality of tree buds, thereby leading to growth decline and tree mortality in northeastern North America (Payette, 2007). A weakening of the tree growth response to growing season temperatures has also been observed in subarctic regions, a phenomenon that could be attributed to temperature-induced drought stress counteracting the favorable growth conditions of warmer summers (Barber et al., 2000;D'Arrigo et al., 2008;Andreu-Hayles et al., 2011). As the frequency and intensity of acute climatic events are expected to be exacerbated by further climate change (Bell et al., 2004;Dai, 2013;Zohner et al., 2020), the integration of acute climatic events in climategrowth modeling is necessary to better understand and anticipate forest responses to climate change. However, such events are highly variable in duration and intensity, which brings additional challenges related to their detection and to the quantification of their potential effects on forest ecosystems.
In this study, we evaluated the influence of growing season heat accumulation and the occurrence of acute climatic events on the radial growth of black spruce near its northern distribution limit. Our goal was to gain novel insights into tree growth dynamics in subarctic environments by pairing high-resolution daily climatic data with tree-ring series from remote forest locations in Nunavik (Québec, Canada). More precisely, we aimed to (1) identify climatic trends over the last decades, including growing season heat accumulation and the occurrence of acute climatic events, i.e., late spring frost, heat waves and cold waves, and (2) quantify their respective effects on annual radial tree growth in the southern and the northern parts of the forest-tundra ecotone in Nunavik.

Study Area
Sampling was conducted in two remote regions of the Taiga Shield Ecoregion (Stralberg et al., 2020) in eastern Nunavik (Québec, Canada). The first region is located near the southern limit of the forest-tundra ecotone (Payette et al., 2001;Scheffer et al., 2012), about 125 km northeast of Schefferville (55 • 17-55 • 51N; 064 • 48'−065 • 16 W). The climate in this area is subarctic with mean annual temperatures ranging from −5.8 to −3.4 • C and mean annual precipitation of 620-898 mm during our study period (Environment Canada, 2019;CEN, 2020). The second region is located ca. 400 km further north, near the northern limit of the forest-tundra ecotone (57 • 50-58 • 36N; 064 • 46'−065 • 49 W). The climate in this region during the study period was characterized by slightly lower temperature and precipitation, with mean annual temperatures ranging from −6.5 to −3.9 • C and mean annual precipitation of 430-620 mm (Environment Canada, 2019;CEN, 2020). The whole study area is underlain by granitic and gneissic rock of Precambrian age constituting a vast plateau that is eroded by major rivers, creating several valleys (Payette, 2007). These lowlands allow continuous lichen woodlands to stretch up to the Arctic ocean in the Ungava Bay, near the northern tree line in Nunavik. Lichen woodlands are characterized by an open structure (above 40% of tree cover) that includes only a small number of dominant trees (Payette and Delwaide, 2018). The main surface deposits in these valleys are shallow tills and the dominant tree species are black spruce, white spruce (Picea glauca (Moench) Voss) and eastern larch (Larix laricina (DuRoi) K. Koch). As reported by Payette (2007), this study area is characterized by an absence of recent wildfires and insect outbreaks.

Data Collection and Tree-Ring Chronologies
Like most of Nunavik, neither region is accessible by road, so the data were collected as part of a canoeing expedition. Before field work, potential study sites of lichen woodland dominated by mature black spruce stands were identified along the Depas, George and Koroc rivers, based on the 2018 Québec Northern Vegetation Map produced by the Ministère des Forêts, de la Faune et des Parcs (MFFP) of the Government of Québec. Potential sites were visited in August 2018 and only sites meeting the following three conditions were retained for sampling: (i) established on well-drained sand terraces, (ii) dominated by large black spruce trees, and iii) characterized by the presence of latesuccessional lichen species (Cladonia stellaris Opiz). Accordingly, two sampling sites were retained in the southern region, and four in the northern region (Figure 1). For each of the sampling sites, wood disks were collected 45 cm aboveground from 10 dominant black spruce trees that were randomly selected. Wood disks were gradually air-dried and finely sanded to allow complete delineation of the latewood boundaries. For each tree, annual ring widths were measured on two opposite radii on the crosssection of the disk using a Velmex micrometer (±0.002 mm; see Supplementary Figure 1 for site-specific mean tree-ring series). One tree from the sixth site had to be removed from the sample because of advanced deterioration induced by decay that prevented accurate ring-width measurement. Each individual tree-ring series was standardized using a spline function with a 66% frequency response to remove the age-related trend and highlight the climatic signal (Cook and Peters, 1981). Individual standardized master chronologies were constructed for each study site and each region, using all respective tree-ring series. Master chronologies were then used for detection of missing or false rings, which was conducted using the COFECHA software (Holmes, 1983), and to assess the quality of the climatic signal. The high site and regional expressed population signal (EPS ; Table 1) values indicated adequate statistical robustness among the tree-ring series, and a strong climatic signal (Wigley et al., 1984).

Tree-Ring Analysis
In addition to the standardized growth chronologies described above, we assessed the occurrence of abrupt growth reductions for each individual tree-ring series based on a method proposed by Das et al. (2007). An abrupt growth decline was defined as a year-to-year negative percentage growth change (PGC) greater than one standard deviation below the mean of all the negative year-to-year PGCs in an individual tree-ring series. The mean (± SD) threshold used to identify abrupt growth declines across all individual series was a growth reduction of 24.8 ± 4.1% between subsequent years.

Climate Data
Daily mean, minimum and maximum temperatures were obtained from Environment Canada and Nordicana D (CEN, 2020). We used two main long-term climatic series that were the closest to our study sites and with comparable elevation: ′′ W) and Schefferville were used to fill the Kuujjuaq dataset. In both regions, temperature data from the paired weather stations were strongly correlated to the main series, with coefficients of determination (R 2 ) >0.9. The period covered by the reconstructed daily temperature series was from 1953 to 2017 for the southern region and from 1947 to 2017 for the northern region. Daily precipitation data are very scarce in northern Québec and only showed weak correlations between associated weather stations (R 2 < 0.3), indicating important regional variations. Because of the limited accuracy of the available data, we chose not to include precipitation values in our analyses.
For each region and each year, we first identified the start and end dates of the growing season. We defined the start date of the growing season as the fifth consecutive day during which the daily mean temperature was > 4 • C, while the end date was set as the fifth consecutive day after August 1st for which the daily minimum temperature was < 0 • C (adapted from Huang et al., 2010). Then, for each growing season, we calculated the cumulative growing degree-days (GDD) when the daily mean temperature reached >4 • C, which was used as a potential predictor variable for the interannual variation in growth (see below).

Acute Climatic Events
Warm spring temperatures are environmental cues that initiate cold dehardening of tree tissues during budburst (Chuine, 2010). The amount of warming that precedes a late spring frost therefore increases the vulnerability of newly emerging shoots to freezing temperatures (Bourque et al., 2005;Chamberlain et al., 2019;Zohner et al., 2020). Accordingly, we characterized potential late spring frost damage by defining an index (hereafter referred to as spring frost index), which was calculated as the accumulated degree-days starting when the daily maximum temperature reached 4 • C and ending when the daily minimum temperature subsequently decreased below −4 • C (Bourque et al., 2005;Zohner et al., 2020). The upper threshold approximates the temperature at which biological activity is activated (Braathe, 1995(Braathe, , 1996, while the lower threshold approximates the FIGURE 1 | Location of the six sites where black spruce wood disks were collected in 2018 (black dots, 1-6) in Nunavik (Québec, Canada), and the two weather stations from which daily temperatures were compiled for analyses (black triangles). Schefferville data were used for the southern region (sites 1 and 2), and Kuujjuaq data were used for the northern region (sites 3 to 6). Quebec's bioclimatic zones follow the map produced by the Direction des inventaires forestiers (MFFP, unpublished). RBAR is the correlation between series, SNR is the signal-to-noise ratio and EPS is the expressed population signal.
Frontiers in Forests and Global Change | www.frontiersin.org temperature that exceeds the freezing tolerance of most tree species during budburst (Chamberlain et al., 2019). For each individual year of the study period in both the southern and northern regions, the most intense event (maximum accumulated GDD) was retained as a potential predictor of the yearly growth variation. We defined heat waves and cold waves as periods during which the daily temperatures were extreme compared to the average conditions in each of the southern and northern regions. To identify heat waves, we used a threshold corresponding to one standard deviation above the mean of the overall daily maximum temperature of all growing seasons for which we had tree-ring series (adapted from Oswald et al., 2018). We measured the intensity of heat waves as the accumulated degree-days for all periods during which the daily maximum temperature remained above that threshold. Similarly, we used one standard deviation below the mean daily minimum temperature for the coldest months of the year (December to March, inclusive) as the threshold for cold waves. The intensity of cold waves was measured as the accumulated absolute degree-days for all periods during which the daily minimum temperature remained below that threshold. For each individual year of the study period in both the southern and northern regions, the most severe event (maximum accumulated degree-days) of both heat and cold waves were retained as potential predictors of yearly growth variation. These two metrics were defined as the cold wave and heat wave indices.

Statistical Analyses
Separate growth models were developed for the northern and southern regions using the standardized tree-ring series. All candidate explanatory variables were standardized prior to fitting models to facilitate the comparison of their respective effects (Gelman, 2008). We first fit a reference model for each region using accumulated GDD during the growing season (Huang et al., 2010). The reference models were linear mixed-effects models with a tree-level random effect (Moreau et al., 2020), as well as a continuous autoregressive (AR1) correlation structure to account for the interannual correlation between tree-ring measurements (Hartmann et al., 2008). We used mixed-effects models to allow for heterogenous responses to acute climate events, which can be captured by a tree-level random effect. We then performed model selection procedures to assess the influence of acute climatic events on annual growth in each region over and above that of accumulated GDD during a given growing season. We built a total of eight additional models by adding and combining intensity indices for spring frosts, heat waves and cold waves, including a full model with all three indices. The resulting models were systematically compared to an intercept-only (null) model. The selection of the best model was based upon the second-order Akaike's information criterion (AICc). Because the model selection procedure did not identify a best single model for explaining the response variable (AICc weights < 0.90), we used model averaging to compute modelaveraged coefficients with a 95% unconditional confidence interval (CI) to evaluate the effects of the different variables (Mazerolle, 2006). We considered variables with confidence intervals excluding zero to be good predictors of annual tree growth (Mazerolle, 2006). Models were developed using the "lme" function of the "nlme" R package (Pinheiro et al., 2018) and model selection was performed using the "AICcmodavg" package (Mazerolle, 2015).
We also conducted a separate analysis to quantify the probability of an abrupt growth decline (defined above) following the occurrence of an acute climatic event, using mixed-effects logistic regressions applied to individual tree-ring series. We used a binary variable indicating the occurrence of an abrupt growth decline as the response and tested for the influence of the acute climatic events that had previously been identified as negative predictors of annual indexed tree growth using their corresponding intensity indices. Again, models were developed independently for the southern and the northern regions. Each logistic model included a tree-level random effect. The "glmmTMB" function of the "glmmTMB" package (Brooks et al., 2017) was used to develop the models. For both the linear and the logistic models, visual inspection of the model residuals was used to confirm that the model assumptions were met, i.e., homogeneity of variance, normality of residuals, presence of outliers and over-dispersion. The variance inflation factor (VIF) was also calculated between candidate variables with VIF < 5 used as a threshold to avoid multicollinearity issues (Zuur et al., 2010). All analyses were performed in the R statistical programming environment, version 3.5.2 (R Core Team, 2019).

Climatic Trends
Growing season length and accumulated GDD increased over time in both the northern and southern sites. The growing season length showed a significant increasing trend at a rate of 3 days per decade in the north (p = 0.004), while a marginally significant increasing rate of 2 days per decade was observed in the south (p = 0.057). Accumulated GDD during the growing season showed similar trends, with a 3.1 • C.days/year increase in the north (p < 0.001) and a 3.2 • C.days/year increase in the south (p < 0.001, Figure 2A).
The intensity of heat waves showed increasing trends over time, with greater differences in the north than in the south. Heat wave intensity increased by 0.12 • C.days/year in the north (p = 0.040; Figure 2C) compared to a 0.061 • C.days/year increase in the south (p = 0.470; Figure 2C). Changes in spring frost intensity also showed increasing trends over time, but these were not statistically significant in either region (p = 0.213 in the north and p = 0.679; Figure 2B). Changes in the intensity of cold waves showed a negative but non-significant trend during the observation period (p = 0.168 in the north and p = 0.359; Figure 2D).

Variables Influencing Annual Tree Growth
In the northern region, the model selection process indicated that indexed growth was best explained by the full model, with an AICc weight of 0.58 ( Table 2). The second most plausible model was the model that included accumulated GDD during the growing season, and spring frost and heat wave indices FIGURE 2 | (A) Accumulated growing degree-days (GDD) during the growing season (GS) for the northern region (left panel; linear regression: y = 516.69 + 3.05; p < 0.0001) and the southern region (right panel; y = 626.40 + 3.16x; p < 0.0001); (B) accumulated degree-days during the most severe spring frost of each year for the northern region y = 3.11 + 0.04x; p = 0.21) and the southern region (y = 8.19 + 0.03x; p = 0.68); (C) accumulated degree-days during the most severe heat wave of each year for the northern region (y = 15.43 + 0.12x; p = 0.04) and the southern region (y = 24.42 + 0.06x; p = 0.47): and (D) absolute accumulated degree-days during the most severe cold wave of each year for the northern region (y = 7.79 -0.02x; p = 0.17) and the southern region (y = 9.11-0.017x; p = 0.36). Confidence interval predictions (shaded area) were calculated with α = 0.05. Solid lines indicate a significant relationship while dashed lines indicate a non-significant relationship. as predictors, with an AICc weight of 0.42. Model averaging indicated that growing season GDD had a positive effect on indexed growth (model-averaged β ± unconditional SE = 4.54 ± 0.67; 95% CI: 3.22, 5.85), spring frost events had a negative effect (-4.96 ± 0.07; 95% CI: −6.33, −3.58), and heat waves had a positive effect (3.72 ± 0.70; 95% CI: 2.30, 5.15). The effect of cold waves was not significant (0.64 ± 0.74; 95% CI: −0.82, 2.10). In the north, the effect size of spring frost events was of sufficient magnitude to offset half of the additional growth that would otherwise have occurred over the course of a warm growing season (Figure 3A). For example, in years with little to no thawing, the indexed growth is predicted to increase from 0.9 Tested variables are accumulated growing degree-days (GDD) during the growing season (GS), accumulated GDD during a spring frost (SF), accumulated degree-days during a heat wave (HW) and accumulated degree-days during a cold wave (CW). K represents the number of estimated parameters, LL the log-likelihood, AICc the delta AIC, and AICcWt the AICc weight. All models included an intercept among the fixed effects as well as a tree-level random intercept. Frontiers in Forests and Global Change | www.frontiersin.org to 1.1 across the full observed range of accumulated GDD, but a single late spring frost reduces growth down to 1.0 in the warmest years. However, since heat waves had a positive effect of roughly similar magnitude (Figure 3B), these acute events effectively offset each other if they occur in the same year ( Figure 3C).
In the southern region, model selection indicated that the full model was the most plausible model for explaining indexed growth, with an AICc weight of 0.71 ( Table 3). The secondbest model had an AICc weight of 0.20 and included growing season GDD, as well as spring frost and heat wave indices. Model averaging indicated that growing season GDD had a positive effect on indexed growth (3.12 ± 0.80; 95% CI: 1.56, 4.68), as did late spring frosts (2.18 ± 0.95; 95% CI: 0.32, 4.04), and heat waves (2.16 ± 0.92; 95% CI: 0.36, 3.97), while the effect of cold waves was again not significant (−1.51 ± 1.11; 95% CI: −3.68, 0.66). Overall, the effect size of the acute climatic events was lower in the southern region compared with the northern region (Figure 3). However, since spring frosts also had a positive effect in the south, their combined effect was sufficient to double the additional growth that occurs over the course of a warm growing season (from about 0.9 to 1.0 without acute events, vs. 1.0 to 1.1 with acute events).

Abrupt Growth Decline
Several abrupt growth declines were detected during the study period (Figure 4). Among these, high proportions of trees from the northern region experienced abrupt growth declines in 1953 (49% of all trees), 1999 (23%), 2013 (38%), and 2017 (67%), all of which coincided with the occurrence of severe spring frosts (>1 SD above the mean of all spring frost indices). In the southern sites, high proportions of trees experienced abrupt growth declines in 1958 (35%), 1986 (35%), and 1999 (45%), which also coincided with years of severe spring frosts (Figure 4). In this region, additional growth decline events occurred during the 1960-1970 period, but could not be related to climatic stress events (Figure 4). According to the logistic regression analyses, an increasing spring frost intensity was associated with a significant increase in the probability of an abrupt growth decline for both the northern (p < 0.001) and the southern (p = 0.001) sites ( Figure 5). Conversely, increasing GDD during the growing season reduced the probability of a tree experiencing an abrupt growth decline for both the northern (p = 0.004) and the southern (p = 0.004) sites. Overall, the probability of abrupt growth decline was much greater in the northern region, where it reached about twice that found in the south with only half of the frost index intensity.

Regional Growth Dynamics
Across the northern and the southern limits of the foresttundra ecotone, our results showed that black spruce growth responded positively to warm growing seasons. These positive growth responses support previous findings that black spruce trees growing near their northern distribution limit are subject to suboptimal temperatures for growth (Thomson et al., 2009;Pedlar and McKenney, 2017;Fréjaville et al., 2020). Accordingly, our results also revealed a positive impact of heat waves on tree growth, indicating that growth is also promoted by periodic heat waves regardless of the overall temperatures of the growing season. This suggests that periods of unusual warmth did not enhance stand respiration and evapotranspiration sufficiently to induce adverse moisture deficits, thereby confirming the prevalence of temperature as a limiting factor for growth at the forest-tundra ecotone. This is in line with previous studies that found little effect of summer precipitation on tree growth in Nunavik (Gamache and Payette, 2004;Huang et al., 2010;Dufour-Tremblay et al., 2012;Nicault et al., 2015), but contrasts with the increasing precipitation limitation observed in northwestern North America (D'Arrigo et al., 2008 in Alaska;Beck et al., 2011;Girardin et al., 2016 in Canada). While northwestern North America experienced an extensive warming from the mid-1970's to the early 2000's (Alaska; Bieniek et al., 2014), extensive warming in Nunavik only started in the mid-1990's (Allard et al., 2007). Hence, the vegetation response to the increased temperatures, such as the accelerated shrub expansion Naito and Cairns, 2015), is asynchronous between those two regions. These contrasting responses to climate change between the northwestern and northeastern parts of the continent may thus result from a transitory positive effect of warming on tree productivity, which can eventually be counteracted by limited water availability D'Orangeville et al., 2018). Tested variables are accumulated growing degree-days (GDD) during the growing season (GS), accumulated GDD during a spring frost (SF), accumulated degree-days during a heat wave (HW) and accumulated degree-days during a cold wave (CW). K represents the number of estimated parameters, LL the log-likelihood, AICc the delta AIC, and AICcWt the AICc weight. All models included an intercept among the fixed effects as well as a plot-level random intercept.  While trees showed positive responses to warm climatic conditions, the occurrence of cold waves was unrelated to tree growth dynamics in our study region. Our results confirm the greater influence of summer rather than winter temperatures on black spruce growth in northeastern Nunavik, as previously reported by Nicault et al. (2015). Although extreme cold conditions were associated with growth reduction and mortality events in our study area in the Nineteenth century (Payette, 2007), the warming observed in Nunavik over the last decades appears to be associated with a decline in cold wave intensity. Even though our data did not provide a statistical confirmation of the significance of this trend, it is possible that warmer conditions have attenuated the impacts of cold waves on tree growth.

Growth Responses to Late Spring Frosts
Despite the positive effects of warm climatic conditions, we found a strong negative impact of late spring frosts on northern tree growth, as spring frost intensity decreased annual growth and substantially increased the probability of abrupt growth decline. Moreover, some of the years with significantly reduced radial growth corresponded to the years with the most severe spring frosts, suggesting that they are also important determinants of tree growth in the forest-tundra ecotone. While late spring frosts have previously been related to severe episodes of forest dieback and growth reductions in temperate forests of North America (Auclair et al., 1996;Bourque et al., 2005;Moreau et al., 2020) and Europe (Menzel et al., 2015;Vanoni et al., 2016;Príncipe et al., 2017;Gazol et al., 2019;Vitasse et al., 2019), this study provides, to the best of our knowledge, the first report of their impact on tree growth dynamics in subarctic environments. This finding also corroborates the potential negative impact of spring frosts on boreal forest growth dynamics (Kullman, 1989;Suvanto et al., 2016Suvanto et al., , 2017Marchand et al., 2019), and marks a new milestone toward a better understanding of the limiting climatic conditions encountered in the forest-tundra ecotone.
Unlike trees from the northern region, trees from the southern region appear to have generally benefited from late spring frosts, as revealed by their overall positive impact on annual growth. Several processes may contribute to the contrasting effects of spring frosts between our two study regions. First, the timing of spring phenological events may partly explain the greater vulnerability of northern trees to spring frosts (Suvanto et al., 2016;Marchand et al., 2019). Because of their adaptation to short growing seasons, northern trees initiate biological activities at colder temperatures in the spring (Westin et al., 2000;Aitken and Hannerz, 2001;Kalliokoski et al., 2012). While this adaptation optimizes the growing season length in cold environments, it also leads to greater chances of frost injuries following a warming that had previously dehardening of tree tissues (Kreyling et al., 2014;Suvanto et al., 2016;Chamberlain et al., 2019). By maintaining winter dormancy until higher temperature thresholds are reached, southern trees may have a significantly reduced vulnerability to early spring frost events, allowing them to benefit from warm spring temperatures thereafter (Kreyling et al., 2014;Suvanto et al., 2016). However, it is worth mentioning that southern trees were not completely invulnerable to late spring frosts, as suggested by the increasing growth decline probability in years of high intensity spring frosts recorded in the southern region.
The increasing trend in the intensity of spring frost events over the observation period, although not confirmed statistically, appeared to be more prevalent in the north than in the south, which may have amplified the vulnerability of northern trees to such events. Indeed, the increase in spring frost intensity occurred over a relatively short period, allowing little time for tree-level acclimation processes, such as adaptive cold stress memory (Walter et al., 2013). The resulting maladaptation to changing climatic conditions may also be amplified in locallyadapted range margin populations, where reduced gene flow prevents rapid adaptations (Nunes et al., 2009;Kreyling et al., 2014). Lastly, in cases where the snow cover has rapidly melted during a warm spring, the return to freezing conditions may injure shallow roots through the freezing of tissues and the development of hinder root pressure that disrupts root processes (Cox and Zhu, 2003;Gaul et al., 2008;Kreyling, 2010;Comerford et al., 2013). The lower amount of precipitation recorded in the northern region during our study period is likely to have increased the vulnerability of northern trees to such injuries.

Perspectives in the Context of Global Change
Over the last few decades, rapid changes in climatic conditions have occurred over our study area at the forest-tundra ecotone. These include increasing growing season length and temperatures, as well as an increase in the intensity of heat waves, which are in line with the climatic trends observed across Canada (Zhang et al., 2011;Seneviratne et al., 2012;Donat et al., 2013). Moreover, in accordance with the greater amplitude of climate change at higher latitudes recorded at the global scale (ACIA, 2004;IPCC, 2014), we observed greater changes in climatic conditions in the northern part of our study area. With projected increases in mean annual temperatures ranging from 5 to 10 • C by 2100, Nunavik is expected to experience the greatest warming in the province of Québec (Ouranos, 2015), which will likely induce major changes in northern ecosystem dynamics.
Since black spruce inhabits colder than optimal environments across the northern and the southern limits of the forest-tundra ecotone, further warming should continue to be beneficial in the near future. However, this beneficial effect on boreal forest tree growth may be transitory, as climate warming may continue beyond the optimal temperature range and induce substantial growth declines through limited water availability (D'Orangeville et al., 2018). Eventually, a transition from temperature-limited to precipitation-limited ecosystems, such as that observed in northwestern North America may occur, although when this transition might occur remains unknown. The ongoing increase in heat wave intensity may also exacerbate this process and affect northern ecosystem dynamics by promoting periodic drought-induced growth declines. Ultimately, such a shift from a positive to a negative effect of heat waves may interact with late spring frost and generate cumulative climatic stresses that are known to promote long-term shifts in forest ecosystem dynamics (Anderegg et al., 2015;Nolet and Kneeshaw, 2018;Moreau et al., 2020). However, an increasing precipitation trend has been observed over the last decades in northern Québec and is expected to continue until 2100 (Ouranos, 2015), which may mitigate the potential negative effects of extensive warming.

CONCLUSION
By using daily-level climate data, we identified relevant climatic variables influencing black spruce growth and demonstrated the influence of acute climatic events on trees growing near the northern treeline in Nunavik, Canada. We also report an increasing trend in the intensity of climatic stress events, which underscores the importance of explicitly incorporating acute climatic events into modeling efforts to better evaluate the impact of ongoing climate change on forest growth dynamics. Toward this end, further studies should pursue the identification of additional biologically meaningful climatic metrics related to tree growth, even though the scarcity of high temporal resolution climatic data, particularly in subarctic regions, represents an additional challenge. The relationships between environmental cues triggered by climatic stressors and tree physiology also need to be carefully investigated at a fine scale, across ecosystems and between species, in order to better identify key thresholds likely to affect tree growth dynamics. Such information is required to provide novel insights into the ongoing tree growth responses to rapid climate change, and to ultimately help support future adaptation, wherever possible.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author, Guillaume Moreau, upon reasonable request.

AUTHOR CONTRIBUTIONS
GM and AA contributed to conception and design of the study. GM collected the data. CC and DA organized the database and performed the statistical analysis. GM and CC wrote the first draft of the manuscript with contributions from JC. All authors contributed to the interpretation of the results and to the manuscript revision. All authors read and approved the submitted version.