Significance Limited knowledge about the mechanistic drivers of forest growth and responses to environmental changes creates uncertainties about the future role of circumpolar boreal forests in the global carbon cycle. Here, we use newly acquired tree-ring data from Canada’s National Forest Inventory to determine the growth response of the boreal forest to environmental changes. We find no consistent boreal-wide growth response over the past 60 y across Canada. However, some southwestern and southeastern forests experienced a growth enhancement, and some regions such as the northwestern and maritime areas experienced a growth depression. Growth–climate relationships bring evidence of an intensification of the impacts of hydroclimatic variability on growth late in the 20th century, in parallel with the rapid rise of summer temperature.

Abstract Considerable evidence exists that current global temperatures are higher than at any time during the past millennium. However, the long-term impacts of rising temperatures and associated shifts in the hydrological cycle on the productivity of ecosystems remain poorly understood for mid to high northern latitudes. Here, we quantify species-specific spatiotemporal variability in terrestrial aboveground biomass stem growth across Canada’s boreal forests from 1950 to the present. We use 873 newly developed tree-ring chronologies from Canada’s National Forest Inventory, representing an unprecedented degree of sampling standardization for a large-scale dendrochronological study. We find significant regional- and species-related trends in growth, but the positive and negative trends compensate each other to yield no strong overall trend in forest growth when averaged across the Canadian boreal forest. The spatial patterns of growth trends identified in our analysis were to some extent coherent with trends estimated by remote sensing, but there are wide areas where remote-sensing information did not match the forest growth trends. Quantifications of tree growth variability as a function of climate factors and atmospheric CO 2 concentration reveal strong negative temperature and positive moisture controls on spatial patterns of tree growth rates, emphasizing the ecological sensitivity to regime shifts in the hydrological cycle. An enhanced dependence of forest growth on soil moisture during the late-20th century coincides with a rapid rise in summer temperatures and occurs despite potential compensating effects from increased atmospheric CO 2 concentration.

Circumpolar boreal forests are estimated to store ∼53.9 Pg of carbon or ∼14% of terrestrial vegetation biomass (1). These regions are currently experiencing accelerated changes, including warmer and longer growing seasons, tree line expansion, species migration, increased frequency and severity of drought, and increases in the frequency and severity of disturbances (2⇓⇓⇓⇓⇓⇓⇓–10). These changes create uncertainty about the boreal forests’ future role in the global carbon cycle (11). Adding to this uncertainty is the discrepancy over recent changes in the productivity of boreal and other northern latitude forests. Some empirical evidence suggests increases in the forest productivity (12⇓–14), whereas other studies suggest decreasing productivity over the last decades (7, 8, 15⇓–17). Furthermore, inversion and process-based ecosystem models indicate large carbon sinks (7, 8), whereas field-based bottom-up approaches suggest smaller carbon sinks or small carbon sources (3, 18), or large sinks (19). Quantifying the response of boreal forests to environmental changes and the subsequent effects on the global carbon cycle thus remains a pending, interdisciplinary scientific challenge (11, 18, 20).

Observations of vegetation productivity at northern latitudes may be obtained using either ground- or satellite-based observations, which may represent different components of forest productivity, that is, stem vs. leaf level. Satellite-derived indices of photosynthetic activity, such as the normalized difference vegetation index (NDVI), are now widely used for quantifying responses of forested areas to environmental changes (21, 22). However, a recent comparison of commonly used NDVI datasets showed marked differences in their representation of mean seasonal vegetation productivity and decadal-long trends (21), which stresses the necessity for large-scale, ground-based observations. Ground-based observations of vegetation productivity may originate from permanent sample plots (13, 17) or from tree-ring analyses (12, 14, 16, 23), which differ in the temporal resolution (decadal vs. annual) and spatial scale (stand vs. tree). Both are prone to bias because of unbalanced spatial sampling. Permanent sample plots may be preferentially established in the more productive forests (24). On the other hand, tree-ring data have been typically focused on the dry or cold marginal limits of forest distribution (23, 25). Furthermore, uncertainty in determining the representative vegetation productivity trends for northern latitudes is exacerbated by short observational time series inherent to many of these data sources. The resulting limited quantification of growth, and subsequent uncertainties in mechanistic drivers and predictions of the ecosystem’s fate, emphasize the need for continuous, highly resolved, and spatially extensive datasets that are statistically representative and constitute an unbiased estimate of the forested regions.

Here, we examine the aboveground growth responses of the Canadian boreal forest to climate change over the past ∼60 y based on a new ground plot network that provides a much more representative sampling of this major North American biome. Information was obtained for a total of 19 boreal tree species from annually resolved and absolutely dated ring width measurements from 2,807 trees at 598 ground plots of Canada’s National Forest Inventory (NFI) (Fig. 1 and SI Appendix, part A). We examine growth trends before and during the satellite era and compare them with trends observed from NDVI datasets. We relate seasonal temperature, soil moisture, and atmospheric carbon dioxide concentration ([CO 2 ]) variation to long-term changes in annual tree growth rates. We also evaluate whether the observed growth trends are consistent with the expected growth rate responses to climate change. Across the vast territory under study, mean annual temperatures increased between 0.5 and 3.0 °C during the 20th century along with global increases in atmospheric CO 2 concentrations (9). Given the cold climatic conditions, it could be expected that these changes should have led to increases in forest growth in this subcontinental region.

Fig. 1. Distribution of Canada’s National Forest Inventory (NFI) sample plots for tree-ring analysis and of the main terrestrial ecozones under study. The nonboreal Pacific Maritime, Mixedwood Plains, and Atlantic Maritime ecozones are included in the study owing to coexisting species with the boreal ecozones.

Conclusions Our analyses of a new methodologically standardized tree-ring dataset covering Canada’s boreal forest provide insights into the growth responses of this ecosystem to climate change. Although revealing no overarching “growth enhancement” or “growth decline” in recent years, results do point to significant regional- and species-related trends in growth. The observed link between climate variation and growth variability revealed unique evidence of an intensification of the impacts of hydroclimatic variability on growth late in the 20th century, in parallel with the rapid rise of summer temperature. Such response can be attributed to annual growth variability in these forests being mainly driven by negative sensitivity to summer temperature (warmer summers leading to less growth) and positive sensitivity to summer soil moisture (more moisture leading to more growth), albeit species-specific variations in these responses can be found. We caution that neither mortality nor disturbance impacts are addressed here and that therefore extrapolation of future carbon storage in these forests is neither straightforward nor simple to achieve (3, 6, 8, 42). Furthermore, the complex interplay of different biotic and abiotic drivers of boreal forest productivity, including unpredictable wildfires, insect population outbreaks, and other disturbances, may yield counterbalancing effects on the overall net carbon sequestration (6, 42, 43). Notably, insect-induced collapses in growth in western and eastern Canadian forests played a major role in defining growth trajectories at the turn of the 21st century (3, 44, 45), but their relative effects on growth are only partially, if at all, captured by this current NFI tree-ring analysis. Although remote-sensing products can provide insights into the impacts of these phenomena and of climate change, our results suggested extensive areas of disagreement between forest growth trajectories and remotely sensed NDVI trends: the accelerated growth over large regions was not necessarily correlated with greening and, inversely, with browning where trees experienced a slower growth. Complex effects of climate on growth rates, unpredictable future disturbances, and uncertainties in trends assessed from remotely sensed forest productivity data provide justification for ground-based monitoring of species-specific trends, disturbance impacts, and responses to climate change at local to regional scales.

Materials and Methods Tree-Ring Data. The tree-ring width dataset is based on increment cores sampled during the establishment of Canada’s NFI plot survey (24), representing an unprecedented degree of sampling standardization for a large-scale dendrochronological study. The sampled plots (n = 598) were designed to be representative of the distribution of species and growing conditions in the managed forests of Canada (Fig. 1). At each NFI plot, samples were collected from up to 10 living dominant and codominant trees (diameter at breast height, >5 cm; median n per plot was four trees and varied depending on the year of survey, with more recent surveys having larger n; see SI Appendix, part A, for further details). Samples were processed, crossdated, and measured using standard dendrochronological methods (SI Appendix, part A). We only used core samples covering a minimum of 20 y. Over 80% of the samples either contained the innermost ring (pith) or were estimated to lie within a centimeter of it. A high proportion of the samples originated from plots located on moraine (44% of plots) and lacustrine (19% of plots) deposits, followed by bogs (8% of plots; one-half of them north of 60°N), and unspecified organic materials (e.g., sedge, sphagnum; 10% of plots). The remaining 19% of plots represented other types of parent materials (e.g., fluvioglacial, bedrock). Ring width measurement series were scaled to provide annual estimates of tree BAIs (in square centimeters per year) using the following: B A I = π R t 2 − π R t − 1 2 , [1]where R t and R t − 1 are the cumulative stem radial increments (in centimeters) at the end and beginning of a given annual ring increment, respectively. BAI and derived metrics thereof (long-term trends in BAI) are more directly related to tree productivity and stem biomass than to diameter increment (46). Vegetation Index Data. Processed NDVI satellite data based on advanced very high-resolution radiometer (AVHRR) raw data were obtained from the latest version (3g) of the Global Inventory for Mapping and Modeling Studies (GIMMS) (47), which covers July 1981 to December 2012 at a temporal resolution of ∼2 wk (24 scenes per y per pixel) and a spatial resolution of 0.083° (∼9 km). For photosynthetically active vegetation, the NDVI ranges from 0 to 1. In this study, we extracted NDVI data for pixels in Canada where the percentage of tree cover was >75% (based on ref. 48), and averaged them by pixel and year to create annual NDVI time series. Climatic and CO 2 Data. For each plot, daily weather data (maximum and minimum temperature; in degrees Celsius), precipitation (sum; in millimeters), relative humidity (in percentage), and vapor pressure deficits (in kilopascals) were obtained for the period of 1950−2010 using BioSIM, which interpolates site-specific estimates from historical weather observations (49) as described in ref. 50. The quantity of available soil moisture was estimated for each month using the quadratic-plus-linear (QL) formulation procedure described in ref. 51, which accounts for water loss through evapotranspiration (simplified Penman–Monteith potential evapotranspiration) and gain from precipitation. Parameter values for maximum and critical available soil water were set at 300 and 400 mm, respectively; the number of weather stations for interpolation was set to 8. We used annual mean atmospheric concentrations of CO 2 recorded at Mauna Loa observatory since 1958 (52); the data were extended to 1950 using estimates from ice cores (53). Quantification of Past Annual Tree Growth. Eliminating the intrinsic age- and size-related growth trend in BAI makes it possible to address annual growth variability independently of biological age or tree size, and furthermore permits comparison of long-term growth trends among species and regions. Departures from the detrended growth rates are interpreted as being caused by climate variability, insect disturbances or other external drivers. Here, we used GAMMs (54, 55) to remove internal biological effects for each species-by-plot combination (873 analyses in total). This approach is based on the modeling of BAI as a function of cambial age and tree basal area. Compared with a continental-scale model approach, the species-by-plot-level models were better able to handle the variation in species responses across climate gradients, varying plot disturbance histories, unequal distributions of species-by-plot sample sizes, and regional climate trends. A logarithmic transformation was first applied to deal with the skewed distribution of BAI (LBAI). The fitted GAMM took on the following form: L B A I i j k t = β i j ⁡ log ( B A i j k ) + s ( cambial age i j k t ) + Z i j k B i j k + ν i j k + ε i j k t , [2]where i represents the species, j represents the plots, k represents the tree, and t represents the year. Cambial age refers to an estimation of tree age based on the ring counts (age n); we did not adjust age n for missing innermost rings as such error, which concerns a minority of samples, does not cause significant bias in tree-ring chronologies (56). Tree identity ( Z i j k B i j k ) was considered as a random effect. We also included an error term ( ν i j k ) with an AR1 (p = 1, q = 0) correlation structure, and the residuals of the model (ε). The smooth terms s of the GAMM were represented using cubic regression splines which degree of smoothness was determined through an iterative fitting process (see ref. 54). The growth model was fitted using the mgcv package, version 1.8–4 (57) in R (58). For quantification of past annual growth variability, a back-transformation was applied to estimated LBAI (59): B A I i j k t ∧ = c f i × exp ( L B A I i j k t ∧ ) , [3]where c f i is a correction factor of species i computed as follows: c f i = exp ( M S E i 2 ) , [4]and M S E i is the mean-squared error computed as: M S E i = ∑ j , k , t ( L B A I i j k t − L B A I i j k t ∧ ) 2 N i , [5]where N i is the number of observations. Annual growth variability (expressed in percent-difference from the long-term mean) was computed as follows: G C i j k t = B A I i j k t − B A I i j k t ∧ B A I i j k t ∧ × 100 . [6]The results of these calculations are hereafter referred to as the GAMM NFI growth data. The ability to detect effects of low-frequency climate variability or rising [CO 2 ] on growth may be limited by the use of detrending techniques that remove long-term growth trends, and by biased sampling of trees that produces spurious trends in growth rates (39, 40). To assess the extent to which the choice of approach for removing age, size, and competition effects might influence our conclusions, we additionally applied a tree level-based statistical procedure, and two uniform ecozone-level procedures to generate three different forms of tree growth indices. Ecozones are areas representative of large and generalized ecological units characterized by interactive and adjusting abiotic and biotic factors (60). The detection of similar trends in all analyses would enhance our confidence that variability is an inherent feature of the tree-ring data. First, a tree-level GNE statistical detrending procedure was used. Here, the ring-width measurement series were rescaled using a power transformation method (61) and detrended using linear (L), negative exponential (NE), or generalized negative exponential (GNE) methods (SI Appendix, part B). The residuals calculated as observed minus fitted estimates are hereafter referred to as the GNE growth data. Note that this GNE approach may remove a greater percentage of low-frequency variance in the LBAI data compared with the GAMM approach. This can limit our ability to track subtle long-term growth trends, particularly when dealing with young trees. Second, for the uniform procedures, we used GAMM models and the RCS technique (56) applied at the scale of ecozones (respectively, GAMM eco and RCS, detailed in SI Appendix, part B). A requirement for the uniform procedures was that at least 20 trees of the same species were available for analysis within an ecozone. Analysis of Forest Growth Trends. Species-by-plot GAMM NFI chronologies were computed using robust averaging of the detrended individual tree BAI data, and then spatially averaged to 1° × 1° grids to obtain continuous yearly raster maps of growth covering 1950–2002 (i.e., the period of maximum sample replication). We used ordinary Kriging in ESRI ArcGIS 10.3 with a 5° radius and an exponential model to fit the semivariograms. The Kriging procedure is essentially a weighted average of annual tree growth around the grid to be estimated. Further testing was conducted to confirm that the results of this approach were robust with respect to the inclusion of rare tree species and chronologies with low tree sample replication (SI Appendix, part G). Ordinary least-square regression analyses (using a linear term alone) were conducted for each grid location, and on NDVI data for the periods 1982−2002 (to correspond to the available NDVI time series) and 1950−2002 (for the growth-based maps). Slope estimates from these regressions were then mapped. Furthermore, NDVI data were used for pairwise correlation with gridded growth time series. These procedures were repeated with the GNE, GAMM eco , and RCS data to ensure that our main findings from the GAMM NFI procedure remain valid under alternative treatments of data. In addition, we examined tree growth variability since the mid-20th century after averaging species-by-plot GAMM NFI chronologies at the level of Canada’s terrestrial ecozones; estimates and their 90% confidence intervals were obtained using bootstrap resampling (61). The procedure was repeated for GNE, GAMM eco , and RCS data. Linear trends for 1950–2002 at the species-by-ecozone level were examined, and direction was interpreted for the level of certainty: “very likely,” if the sign of the trend in GAMM NFI was matched by the three other methods; “likely,” if the sign of the trend in GAMM NFI was matched by two other methods; and “uncertain,” if inconsistency was found between methods; the term “inconclusive” was used when a trend in GAMM NFI was between 0 and ±0.10%⋅y−1. Analysis of Climate Sensitivity. Linear mixed models were used to explore the climate effects on species-by-plot averaged residuals of Eq. 2 ( G A M M L R e s i j k t ) G A M M L R e s i j t = ∑ n = 1 8 β i j n * C l i m j t n + ν i j . [7] An error term with an AR1 (P = 1, q = 0) correlation structure was included in the models. Eight explanatory variables (Clim) were tested: mean summer (June to August) soil moisture (in millimeters) and temperature (in degrees Celsius) of the previous year and current year to ring formation, previous fall (September to November) mean temperature, current spring (March to May) mean temperature, cool season (December to May) total snowfall (in millimeters), and mean annual atmospheric CO 2 concentrations ([CO 2 ] in parts per million by volume). The period of analysis ranged from 1950 to 2010 depending on the length of the individual G A M M L R e s i j t records. Multicollinearity among these variables is moderate to low (SI Appendix, part F). Multimodel selection was performed using the MuMIn package (61). In this approach, we ranked all of the potential models that could be generated with the different explanatory variables according to the second-order Akaike information criterion corrected for small sample sizes (AICc). For each model, we considered its Akaike weights (i.e., normalized model likelihoods) for ranking. Among all ranked models, those with cumulative Akaike weights within 95% were considered the best candidate models for multimodel inference. Parameter and error estimates were then computed through averaging of these multiple models (62). We used four different metrics to report the climate sensitivity: the t statistics (i.e., the coefficient divided by its corresponding unconditional SE), the variable relative importance (sum of “Akaike weights” from the selected models), the 95% adjusted confidence intervals for the parameter coefficients of the average model (species-by-plot), and the Pearson correlation between the estimation from the average model and the response variable for each species-by-plot.

Acknowledgments We acknowledge the important contribution of Graham Stinson and the National Forest Inventory program, from which most of the data were obtained. We thank Christine Simard, Julie Fradette, David Gervais, Catherine McNalty, and Thierry Varem-Sanders for valuable contributions to the laboratory work. This work was made possible thanks to the financial and in-kind support provided by the Canadian Forest Service of Natural Resources Canada.

Footnotes Author contributions: M.P.G., O.B., W.K., X.J.G., and J.B. designed research; M.P.G., O.B., and X.J.G. performed research; M.P.G., O.B., E.H.H., N.E.Z., J.M.M., R.d.J., D.C.F., J.E., and U.B. contributed new reagents/analytic tools; M.P.G., O.B., J.M.M., R.d.J., and X.J.G. analyzed data; and M.P.G., O.B., E.H.H., W.K., N.E.Z., J.M.M., D.C.F., U.B., and X.J.G. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1610156113/-/DCSupplemental.