Strong latitudinal gradient in temperature-growth coupling near the treeline of the Canadian subarctic forest

Boreal forests are experiencing severe climatic changes that vary widely across the broad geographic distribution of the biome. The changes are greatest near the subarctic treeline where trees often exhibit high climatic sensitivity because climatic conditions approach the limits of their physiological tolerance. Despite the importance of subarctic boreal forests, the lack of field-acquired growth data remains a critical issue that limits the generalization of forest productivity models across the entire boreal biome. Using tree-ring chronologies from remote stands distributed along three latitudinal gradients ranging from 65 to 102°W, we investigated recent trends in black spruce growth and their relationships with recent climate warming near the subarctic treeline in eastern Canada. Our results show a generally positive effect of temperature and a negative effect of precipitation, both indicating that black spruce growth is temperature-limited near its northern range limit. However, we observed a strong gradient in temperature-growth coupling within a small latitudinal gradient (about one degree of latitude), where strong temperature constraints appear limited to the northernmost, coldest stands. Moreover, the positive growth response to temperature decreased from wetter to dryer sites and climate-growth coupling declined over the study period in the driest sites. These results suggest that the growth increase associated with warmer temperature may be limited by reduced precipitation and potential moisture limitation. Lastly, our results suggest that acute climatic events have the potential to induce abrupt shifts in tree climate-growth relationships. Such results indicate that the expected beneficial effect of warming on high latitude tree growth may be less generalized and more complex than previously thought in northeastern Canada, perhaps due to factors other than temperature, which might confound the climate-growth coupling southwards. Thus, our results highlight the need for a better understanding of additional growth drivers in these poorly studied regions and for physiologically informed definitions of acute climatic events, in order to refine broad-scale forest productivity modeling.

Boreal forests are experiencing severe climatic changes that vary widely across the broad geographic distribution of the biome. The changes are greatest near the subarctic treeline where trees often exhibit high climatic sensitivity because climatic conditions approach the limits of their physiological tolerance. Despite the importance of subarctic boreal forests, the lack of field-acquired growth data remains a critical issue that limits the generalization of forest productivity models across the entire boreal biome. Using tree-ring chronologies from remote stands distributed along three latitudinal gradients ranging from 65 to 102 • W, we investigated recent trends in black spruce growth and their relationships with recent climate warming near the subarctic treeline in eastern Canada. Our results show a generally positive effect of temperature and a negative effect of precipitation, both indicating that black spruce growth is temperaturelimited near its northern range limit. However, we observed a strong gradient in temperature-growth coupling within a small latitudinal gradient (about one degree of latitude), where strong temperature constraints appear limited to the northernmost, coldest stands. Moreover, the positive growth response to temperature decreased from wetter to dryer sites and climate-growth coupling declined over the study period in the driest sites. These results suggest that the growth increase associated with warmer temperature may be limited by reduced precipitation and potential moisture limitation. Lastly, our results suggest that acute climatic events have the potential to induce abrupt shifts in tree climategrowth relationships. Such results indicate that the expected beneficial effect of warming on high latitude tree growth may be less generalized and more complex than previously thought in northeastern Canada, perhaps due to factors other than temperature, which might confound the climate-growth coupling southwards. Thus, our results highlight the need for a better understanding of additional growth drivers in these poorly studied regions and for physiologically informed definitions of acute climatic events, in order to refine broad-scale forest productivity modeling.

Introduction
Boreal forests cover nearly 30% of the Earth's forested area and play an essential role in the global carbon cycle by storing about 14% of its terrestrial biomass (Pan et al., 2011;FAO, 2020). This biome is experiencing some of the most severe changes in climate on the planet, including rapid increases in temperature and heterogeneous changes in precipitation (Price et al., 2013;Wang et al., 2014). The length of the growing season has also increased in many boreal regions due to warmer spring temperatures and earlier snowmelt, as well as warmer fall temperatures (Vincent et al., 2018). These dramatic changes are likely to influence forest productivity by enhancing tree growth and survival where temperature is limiting (D'Orangeville et al., 2016;Mirabel et al., 2022;Wang et al., 2023). Conversely, in regions where water availability limits tree growth, forest productivity is projected to decline as increasing temperature increases evaporative demand (McKenney et al., 2007;D'Orangeville et al., 2016;Mirabel et al., 2023;Wang et al., 2023). Indeed, due to the broad geographic distribution of the boreal biome, climate change and associated tree growth responses are expected to exhibit considerable spatial heterogeneity, limiting our ability to anticipate how forests will respond (Chagnon et al., 2022).
Boreal forests are dominated by few species with broad geographic distributions. For example, black spruce (Picea mariana (L.) Mill B.S.P.) is the most common boreal tree species in North America, and its range encompasses most of the boreal biome on the continent. Thus, assessing climatic constraints on the growth of black spruce will go a long way toward anticipating the response of boreal forests to climate change. Two recent studies have used tree-ring series spanning broad latitudinal gradients to show that the growth of black spruce is limited by temperature at high latitudes (>49 • N) and precipitation at low (<49 • N) latitudes (D'Orangeville et al., 2016;Chagnon et al., 2022). Black spruce forests at high latitudes are therefore expected to respond positively to projected increases in temperature (Chagnon et al., 2022). Moreover, since northern Canada is also characterized by a longitudinal precipitation gradient, greater precipitation in the east is likely to compensate for increasing evapotranspiration caused by warmer temperatures (Gauthier et al., 2015).
While broad-scale studies such as these have significantly improved our understanding of climatic constraints on forest productivity, they have generally excluded the northernmost part of the boreal biome and overlooked the impact of acute climatic events on tree growth. At higher latitudes, closed-canopy boreal forests transition into the forest-tundra ecotone, a mosaic of lichen woodlands consisting of open forest stands (10-40% tree cover) characterized by a small number of dominant trees and a lichen-dominated understory (Payette and Delwaide, 2018). In North America, these subarctic forests cover about two million km 2 (Johnson and Miyanishi, 1999), thereby forming the second largest ecosystem in the boreal biome (Payette and Delwaide, 2018). Yet, few tree-ring chronologies have been assembled in these remote forests, limiting our understanding of climatic constraints on growth near the northern treeline (D'Arrigo et al., 2009;Nicault et al., 2015;Moreau et al., 2020b). Nevertheless, northernmost trees are generally expected to exhibit higher climatic sensitivity because climate conditions approach the limits of their physiological tolerance (Holtmeier and Broll, 2005). Thus, quantifying growth responses to a rapidly changing climate at high latitudes may provide critical insights that will help anticipate how further climate change will alter boreal forest productivity (Payette et al., 1989;Callaghan et al., 2002).
In this study, we made use of recent sampling efforts conducted in remote forests near the northern treeline in Canada to gain novel insights on the influence of climate on the growth of subarctic black spruce. We paired high-resolution daily climatic data with 118 individual tree-ring series sampled in stands distributed among three latitudinal gradients in eastern Nunavut and Nunavik, which together span a much broader longitudinal gradient from ∼65 to 102 • W. The objectives of the study were to (1) characterize variation in biologically relevant climatic variables over the last decades, including growing-season heat accumulation and the occurrence of acute climatic events, (2) quantify the effect of the climate variables on annual growth at sites spanning a latitudinal temperature gradient as well as a longitudinal precipitation gradient, and (3) assess the temporal stability of these growthclimate relationships.

Study area
The study area is located in the Taiga Shield Ecoregion in northern Canada (Figure 1) and is dominated by mature lichen woodlands. The climate is subarctic and characterized by a latitudinal gradient in temperature (increasing by 0.5-1.6 • C from north to south; Table 1) and a longitudinal gradient in precipitation (increasing by ca. 250 mm from west to east). Between 2018 and 2021, three to four black spruce-dominated mature lichen woodlands (Figure 1) were sampled along latitudinal gradients within each of three study regions, for a total of 10 sampling sites. Each study region is located in the vicinity of the northern treeline, i.e., near the northern boundary of the forest-tundra ecotone (Payette et al., 2001). From west to east, the three study regions include the Kazan River watershed (three sites; thereafter Kazan; visited in 2021), the Hudson Bay watershed (three sites; Hudson, 2018, 2019), and the Ungava Bay watershed (four sites; Ungava; 2018; Figure 1). Study sites in Ungava and Kazan were located along the Kazan, George and Koroc rivers, and were accessed during canoe expeditions. Sites in Hudson were located within a five-kilometer radius around the research stations of Boniface (57 • 45 N, 76 • 10 W) and Lac-à-l'Eau-Claire (56 • 20 N, 74 • 27 W) affiliated with the Center for Northern Studies.
For the Hudson and Ungava regions, potential sites were identified a priori using the Quebec Northern Vegetation Map (Ministère des Forêts, de la Faune et des Parcs (MFFP), 2018), while the Kazan sites were identified in the field. Sites had to meet the following criteria to be included in this study: (i) established on well-drained sand terraces, (ii) dominated by black spruce, and (iii) characterized by the presence of late-successional lichen species (i.e., Cladonia stellaris Opiz). Kazan study sites were located on surficial deposits consisting of thick and continuous tills whereas the study sites in the Hudson and Ungava regions were on shallow tills.  MAT is the mean annual temperature and MAP is the mean total annual precipitation. RBAR is the correlation between the tree ring series within a site, SNR is the signal-to-noise ratio, EPS is the expressed population signal, and start-end is the annual time span of the tree ring chronologies. Wood samples were collected at 130, 30, and 45 cm aboveground in the Kazan, Hudson and Ungava regions, respectively.

Sampling and tree-ring chronologies
Due to their difficulty of access, each study region was visited during separate expeditions, which explains some of the variations in data collection protocols. A total of 125 black spruce trees were sampled in the 10 study sites; 15 trees were sampled in each of the Kazan and Hudson sites, whereas 10 trees were sampled in the Ungava sites. To be selected for sampling, black spruce trees had to be dominant or co-dominant and free of immediate competition. In the Kazan sites, two cores (opposite radii) were sampled, while only one core was collected in the Hudson sites. Cores were sampled 130 and 30 cm aboveground in the Kazan and Hudson sites, respectively. In the Ungava sites, wood disks were sampled 45 cm aboveground. All samples were air-dried and sanded to allow a precise delimitation of latewood boundaries. For both cores and wood disks, annual ring widths were measured using a Velmex micrometer (±0.002 mm). Two opposite radii were measured on the cross section of each disk. Cross-dating of the treering series was done using COFECHA (Holmes, 1983) and master chronologies were created for each site to detect missing or false rings. Each individual tree-ring series was detrended using a spline function with a 0.50 frequency response at a wavelength of 30 years to highlight interannual variability (Klesse, 2021; Supplementary  Figure 1). Detrending was done in the R environment (R Core Team, 2019) using the "dplR" package (Bunn, 2008). Regional expressed population signal (EPS; Wigley et al., 1984), mean correlation between the series (RBAR), and signal-to-noise ratio (SNR) were calculated using the "dplR" package to characterize the tree-ring chronologies. Over the whole dataset, seven trees were removed from the chronologies due to decay in the wood sample that prevented accurate measurements or weak correlations with the master chronology. A total of 118 individual tree-ring series was thus used for the following analyses.

Climate data
For each study site, daily climate was extracted with the R implementation of BioSIM (Fortin et al., 2022), using the "Climatic_Daily" model and the six nearest weather stations (Supplementary Figure 2) for the 1950-2020 period. BioSIM is a software developed by Natural Resources Canada which computes daily climate variables based on daily records of Environment Canada historical weather observations. These observations are then corrected for the site-specific latitude, longitude, elevation, slope, aspect, and distance from the ocean (Régnière et al., 2017). Using the resulting daily climate data, we computed several variables specific to the growing season, defined as starting on the fifth consecutive day during which maximum temperature reaches 4 • C, and ending on the fifth consecutive day during which minimum temperature falls below 0 • C (but only after August 1st, to avoid considering spring frosts as the end of the growing season; Huang et al., 2010;Moreau et al., 2020b). These include growingseason length (GSlength), growing-season precipitation (GSprec), mean growing-season temperature (mean GST), and growingseason degree-days (GSdd), which was calculated by summing the degree-days ( • C) above 4 • C ( Table 2). Finally, we calculated the yearly snowfall (cm) by summing the amount of precipitation that fell while the daily maximum temperature was below 0 • C, after the end of the previous growing season and before the start of the current growing season. These climatic variables were detrended using the same method as the tree-ring series (Klesse, 2021;Ols et al., 2023), i.e., with a 30-year spline and a 50% cut-off.
We also quantified the impact of two types of acute climatic events: late spring frosts and winter thaws. Late spring frosts were defined as temperature excursions below 0 • C that occurred after the start of the growing season, when frost may induce the freezing of parenchyma cells in newly emerging shoots that have initiated cold dehardening following spring warming (Chuine, 2010). We quantified the impact of frost using two indices ( Table 2): the number of days during which the temperature remained below the 0 • C threshold between the beginning of the growing season and the last frost event (LSFlength), and the cumulative freezing (degree-days below the 0 • C threshold) during the frost events (LSFfreezing). The amount of warming that precedes a late spring frost may also increase the vulnerability of tree tissues to freezing temperatures (Chamberlain et al., 2019;Zohner et al., 2020a), so we also quantified the impact of frost using a third index that accounts for the cumulative warming (degree-days above 4 • C; Zohner et al., 2020b) prior to the last frost event (LSFwarming). Finally, in addition to using 0 • C as the threshold for defining the beginning and end of a frost event, we also used two other thresholds (−2 and −4 • C; Vitasse et al., 2019;Zohner et al., 2020b) to calculate each of the three indices ( Table 2).
We defined winter thaws as positive temperature excursions occurring between the end of one growing season and the start of the next, when the daily maximum temperature reached 4 • C. Winter thaws may affect tree growth through several physiological processes. First, winter thaws followed by freezing temperatures have been connected to decreased xylem conductivity induced by winter cavitation (Zhu et al., 2001;Cox and Zhu, 2003). Second, shallow root systems, such as that of black spruce (Elie and Ruel, 2005), are susceptible to thaw-freezeinduced injuries, particularly during extended thaws that melt the snow cover (Comerford et al., 2013). We quantified the impact of winter thaws using two indices ( Table 2): the number of days the temperature remained above 4 • C between the end of one growing season and the start of the next (WTlength) and the cumulative warming (degree-days above 4 • C; WTwarming) during these thaws. These two acute climatic events were not detrended because of their stochastic nature.

Statistical analysis
In a preliminary model selection analysis, we used the individual tree-ring series (n = 118) from the three regions to identify the candidate variables that best quantify the effect of growing-season climate, late spring frost, and winter thaw on growth for the 1950-2020 period ( Table 2). We built individual models for each candidate variable using the tree-level annual indexed ring width (RWI) as the response variable. The best variable was identified using the second-order Akaike Information Criterion (AICc), computed with the "AICcmodavg" R package (Mazerolle, 2006), and only the best variables were selected for further analysis. There were no alternative candidate variables for growing-season precipitation and snowfall, so these were advanced for further analysis. Temporal trends for each of the selected variables were then evaluated at the regional scale using mixed linear models to quantify the relationship between climate variables and year, including a random intercept to allow the climate trends to vary among the study sites (Supplementary Figure 3).
We then used the selected variables to develop separate models for relating annual growth to climate using the individual indexed tree-ring widths series (RWI). Mixed linear models were built independently for each study site since we expected the growth response to vary among sites (see Chagnon et al., 2022). A treelevel random intercept and a first-order autoregression correlation structure taking the form ∼ year | treeID were added to each sitespecific model in order to address temporal autocorrelation, which was detected through autocorrelation and partial autocorrelation functions on the preliminary models' residuals. Models were built using the "nlme" package in the R environment using the following formula: RWI ∼ GSdd + GSprec + GSdd_lag + GSprec_lag + Snowfall + LSFwarming + WTwarming All climatic variables were centered on the sites-specific mean and divided by the overall standard deviation of the variable across all sites prior to modeling to allow comparison between the effect sizes. Using the regression coefficients estimated for each climatic variable and for each site, we then investigated if there was a trend in the growth response related with the site-specific climate. To do so, we computed linear models with "base" R functions using the regression coefficients as the response variable and the site-specific climatic normal (1981-2020) for the growing season temperature and precipitation and snowfall as predictors.
Finally, to assess the stability of the growth-climate relationships over time, we partitioned our dataset using a 30-year moving window and repeated the site-specific investigation of the effect of each predictor on growth for each 30-year window. We then plotted the regression coefficients for each predictor variable against time to detect temporal trends or instability in the magnitude of the growth response. Linear models were then used to confirm temporal trends over the study period for each predictor and for each site. We considered a temporal trend to be present when the slope of the relationship between the regression coefficients and the year (first year or the 30-year moving window) was significant (α = 0.05). Finally, we also computed the partial correlation for each combination of predictor, site and 30-year window (Supplementary Figures 7-9) using the package "rockchalk" (Johnson, 2022).

Spatial variation in growth-climate relationships
Across the study sites, black spruce growth was positively related to current year temperature (GSdd, Figure 2) in all but the warmest site (Kazan 3), but negatively related to current year precipitation (GSprec, Figure 2). The effect sizes of growth response to temperature followed clear spatial patterns with the positive effect of temperature decreasing from north to south. The relationship between the growth response to temperatures and sitespecific climate supported this spatial pattern, with the regression coefficients significantly decreasing from colder to warmer sites (p = 0.008; Figure 3A). Growth response to temperature was also affected by precipitation and decreased from drier to wetter sites (p = 0.02; Figure 3B), in accordance with the north-tosouth gradient. In contrast to temperature, there were no spatial patterns for the growth response to precipitation, nor for the lagged variables, which were not significantly influenced by local climate normals. The effect of the lagged variable was the opposite of that of the current year variables (Figure 2): the previous year's temperature (GSdd_lag) had a negative effect on growth, and the previous year's precipitation (GSprec_lag) had a positive effect on growth.
Snowfall had a negative effect on growth in most sites (Snowfall; Figure 2), and the negative effect seems to increase from west to east along the precipitation gradient. Again, this spatial pattern was supported by the relationship between correlation coefficients and site-specific climate, which indicated a marginal relationship with snowfall (p = 0.08; not shown). Accordingly, the negative effect was greatest where snowfall is greatest (Ungava) and negligible where snowfall is least (Kazan). The effect of late spring frosts and winter thaws (Figure 2) was inconsistent across sites, with both positive and negative effects on growth.

Temporal variability
Growth-climate relationships estimated using 30-year moving windows exhibited various temporal trends across the study period (Figure 4 and Supplementary Figures 4-6). In the Hudson region, the positive effect of temperature on current year's growth began Growth responses (standardized regression coefficients ± SE) to current and previous growing season degree-days (GSdd, GSdd_lag), growing season precipitation (GSprec, GSprec_lag), snowfall (cm), late spring frost (LSF) and winter thaws. For each climate variable, the regions (Kaz, Kazan; Hud, Hudson; Ung, Ungava) are plotted from west to east, and the sites within each region are plotted from north to south, with the northernmost site (1) on the left and the southernmost site (3 or 4) on the right. Significant (α = 0.05) regression coefficients are in blue.
increasing at the end of the 1960s (Figure 4 and Supplementary  Figure 5). The lagged effect of temperature on growth exhibited the opposite temporal trend, becoming more negative over time.
Overall, the effect of the current year growing season temperature on black spruce growth remained greater than that of the previous year. The effect of current year precipitation oscillated over the study period without obvious patterns (with positive, negative, and non-significant slope observed), and remained mostly The effect of the climatic normal (1981-2020) for growing season degree-days (A; R 2 = 0.43) and precipitation (B; R 2 = 0.25) on black spruce growth response to current-year growing season degree-days. The growth response is represented by the standardized regression coefficient of the relationship between indexed ring width and current-year growing season DD for the 1950-2020 period. The lines and shaded areas show the model estimation and associated 95% CI. The symbols represent the study sites, with Kazan sites represented as triangles, Hudson sites as squares, and Ungava sites as circles.
non-significant throughout the observation period. The negative effect of snowfall increased over the study period.
In the wetter Ungava region, similar trends were detected for both current year and previous year temperature, although the magnitude of their effects was less pronounced than in Hudson (Figure 4 and Supplementary Figures 5, 6). The effect of current year precipitation was mostly non-significant over time but switched to significant and positive values in the most recent years. The effect of previous year precipitation, which was positive at the beginning of the study period in all four sites, decreased through time and became non-significant or weakly negative over time (Supplementary Figure 6). The effect of snowfall showed weak trends but remained negative over the whole study period.
Across all Hudson and Ungava sites, the effect of late spring frosts and winter thaws was very unstable, with both positive and negative values recorded for each site and heterogenous temporal trends (Figure 4 and Supplementary Figures 4, 5). This instability was even greater in the Kazan region, where the effect of acute climate events on growth shifted abruptly around 1980 (Supplementary Figure 4), coinciding with particularly extreme winter thaw and late spring frost events in 1980 (Supplementary  Figure 3). The temporal variability of other climatic effects on growth displayed similar shifts (Supplementary Figure 4), although less pronounced than for the acute events. These shifts led to less obvious temporal trends that were often non-linear, and thus less appropriately described with our linear models. Nevertheless, it is worth noting the effect of multiple climatic variables, and especially precipitation (i.e., current and previous year precipitation, snowfall), weakened over time in the Kazan region, with small and often non-significant coefficients in the most recent years (Supplementary Figure 4).

Climatic drivers of black spruce growth
In line with recent broad-scale assessments of growth-climate relationships (D'Orangeville et al., 2018;Marchand et al., 2019;Chagnon et al., 2022;Mirabel et al., 2022), our results indicate that black spruce growth in eastern Canada is constrained by cold temperatures and associated short growing seasons near its northern distribution limit. However, as the positive effect of temperature decreases abruptly from north to south, our results suggest that the beneficial effect of increasing temperatures on tree growth expected at latitudes >49 • N in northeastern Canada (D'Orangeville et al., 2018;Chagnon et al., 2022) may be less generalized and more complex than previously thought. Indeed, the 0.5-1.6 • C increase in MAT at the southern edge of the latitudinal gradients studied here reduced the positive effect of the temperature by more than half, suggesting that strong temperature-growth coupling is spatially limited to the very coldest sites. Vapor pressure deficit, i.e., atmospheric dryness, may be involved in limiting growth in the southernmost sites as it increases from north to south in our study area (e.g., Mirabel et al., 2023). Moreover, non-climatic factors whose effect is masked by temperature limitation in the far north may become more limiting to the south as temperature limitation is alleviated. However, we can only speculate about such factors, which may include soil depth, nutrient availability or other factors acting more locally, such as microtopography (Camarero et al., 2021). Although similar gradients have been observed in other parts of the boreal biome (see Hellmann et al., 2016;Hartl et al., 2021 in Eurasia), to our knowledge, this study is the first Temporal variability of the climatic effects on growth in the northernmost site of each region. The standardized regression coefficients (and 95% CIs) quantify the effect of climate on growth in a given 30-year window, beginning the year plotted on the x-axis. Significant regression coefficients (α = 0.05) are in blue. Note that the y-axis scales differ between plots of different climatic variables.
to report such a steep gradient in black spruce climate-growth relationships near the Canadian treeline.
The growth response to precipitation was negative across the study area, in accordance with recent findings in the northeastern part of the Canadian boreal forest, where increased vapor pressure deficit benefitted tree growth (Mirabel et al., 2022(Mirabel et al., , 2023. If such findings do not seem to indicate a current moisture stress, the lagged positive response to precipitation may indicate that greater water availability has an overall positive effect of tree growth (Mirabel et al., 2023), as the previous-year effect of precipitation is greater than the current-year effect across our study area. Moreover, the positive growth response to temperature decreases from wetter to dryer sites, suggesting that reduced precipitation may limit the growth increase associated with warmer temperature (Mirabel et al., 2023). Similarly, temporal trends in growth-climate relationships indicate a decreasing climate-growth coupling in the last decades at Kazan (Supplementary Figure 4). In northwestern North America, a similar weakening of the temperature effect on spruce growth has been previously observed (i.e., the "divergence problem"), and attributed in part to increasing water deficit (D'Arrigo et al., 2008(D'Arrigo et al., , 2009Andreu-Hayles et al., 2011). The limited precipitation in Kazan may thus partly restrain the growth response to the most recent warming. When added to the decreasing negative effect of precipitation over time, this effect points toward an upcoming moisture limitation for black spruce growth in this part of the study area. Nevertheless, opposing temporal trends were observed across most of our study area, with increasing temperatures being associated with increasing growth in the last decades in Ungava and Hudson. These contrasting temporal trends suggest that future tree growth will depend on local climate and challenge the generalization of the "divergence problem" to the entire subarctic forest, which appears limited to the driest sites (Juday et al., 2015).
Our results also indicate that the growth response to climatic conditions in the previous growing season was the opposite of that observed for the current growing season, as warm and dry seasons negatively impacted growth in the following year. The lagged effect of temperature and precipitation on black spruce growth are mostly consistent across the boreal forest (e.g., Huang et al., 2010;Girardin et al., 2016b;Ols et al., 2018;Marchand et al., 2019;Mirabel et al., 2023). Several explanations have been suggested, including reduced carbohydrate accumulation caused by heat stress (Girardin et al., 2016a;Mirabel et al., 2023), increased resources allocation to roots to improve water assimilation capacity (Lapenis et al., 2013), or production of reproductive structures following a good year of growth (Hacket-Pain et al., 2018;Tumajer and Lehejèek, 2019), all of which may limit resource allocation to radial growth.
Finally, we found a negative effect of snowfall on black spruce growth, which likely results from a delayed start of the growing season induced by late soil thawing and associated root activity under the insulated snow cover (Vaganov et al., 1999;Girard et al., 2011). Along with other recent work , our results underline the necessity to include winter variables into tree growth modeling. This is especially true as changes in winter conditions are projected to surpass those of growing season climate, with a decreasing number of frost-days and a higher proportion of precipitation falling as snow in northern latitudes (Vincent et al., 2018).

Potential for acute events to disrupt growth-climate relationships
Our results suggest that acute climatic events have the potential to disrupt the growth response to climate, as indicated by the abrupt shifts in the Kazan region coinciding with the extreme winter thaws and late spring frosts that occurred around 1980. While the acute climatic events showed a much higher frequency and severity in the Kazan region than elsewhere, the lesser snowfall in the region may lead to an insufficient snow cover to effectively insulate soils and protect tree roots from freezing during late spring frosts or winter thaws (Comerford et al., 2013). The higher severity of late spring frost in the Kazan region, i.e., the greater amount of degree-days cumulated before the last frost event, is also likely to be associated with a more advanced dehardening process of tree tissues than in other regions, leading to more severe damage (Chamberlain et al., 2019). Yet, other factors such as climate oscillations could act on both the occurrence of acute climatic events and tree growth, leading to a confounding effect. For example, the Arctic Oscillation and El Niño-Southern Oscillation were related to both climatic anomalies and altered tree-growth (Buermann et al., 2003), whereas the Pacific Decadal Oscillation was related to decreasing growthclimate coupling in Alaska (Ohse et al., 2012).
Previous studies have also suggested that acute climatic stressors may severely affect growth-climate relationships during the current and following years (Chamberlain et al., 2019;Vitasse et al., 2019;Moreau et al., 2020a,b). However, the evidence remains marginal, probably because of the ecological and computational uncertainties associated with such events. Indeed, climate indices describing acute climatic events that are biologically meaningful for tree growth are difficult to compute for several reasons. First, while evidence of frost injuries associated with late spring frosts have been reported, the timing of cold dehardening and the temperature thresholds leading to tree damage remains unclear (e.g., 0 • C, Zohner et al., 2020a; −2 • C, Vitasse et al., 2019; −4 • C, Moreau et al., 2020a,b, herein), and varies intra-and inter-specifically (Zohner et al., 2017(Zohner et al., , 2020a, which complicates the identification of such events. Moreover, few studies investigated the tree physiological processes associated with frost events in situ (Mayr et al., 2020), as most were conducted in laboratory (e.g., Zhu et al., 2001;Mayr et al., 2003). Second, although gridded climatic data modeling has considerably improved in the recent years, weather stations density and temporal coverage vary spatially and influence the quality of the derived gridded data, which is of particular concern at high latitudes. When computing the climatic variables used herein, we observed that small variations in climatic data led to important changes in the detection of climatic events that use thresholds. Accordingly, we detected far fewer and less severe late spring frosts in Ungava than did Moreau et al. (2020b). This highlights the sensitivity of acute climatic events indices to small variation in climate data, and the necessity to inform and refine the climatic thresholds used to define such stressors through in situ, speciesspecific physiological studies.

Conclusion
Black spruce growth from remote open stands near the northern Canadian treeline follows the broad-scale patterns previously identified, with black spruce growth being temperaturelimited at high latitudes. However, the large gradient in temperature-growth coupling within a rather small area that we observed was unexpected and highlights the fact that factors other than temperature are influencing tree growth near the treeline, leading to less homogenous growth patterns than that suggested by previous broad-scale studies for northeastern Canada. Thus, this variability should be accounted for in broad-scale forest productivity modeling, and a generalized positive effect of increasing temperature on northern tree growth should not be assumed. Overall, by highlighting the spatial and temporal variability of climate-growth relationships and the potential of acute climatic events to induce abrupt shifts of these relationships, our results demonstrate the need for a better understanding of additional growth drivers and tree physiological processes in these poorly studied regions.

Data availability statement
The original contributions presented in this study are included in the article/Supplementary material, further inquiries can be directed to the corresponding author.

Author contributions
CC, GM, AA, and LD'O contributed to the conception and design of the study. CC, GM, and J-PL-F collected the data and organized the dataset. CC performed the statistical analysis and data visualization. CC and GM wrote the first draft of the manuscript with contributions from JC. LD'O, AA, and J-PL-F substantially contributed to the manuscript revision. All authors read and approved the submitted version.