Size-, species-, and site-speci ﬁ c tree growth responses to climate variability in old-growth subalpine forests

. Tree-ring data have become widely used to model tree growth responses to climate variability and gain insight about the potential effects of global warming on forests. We capitalized on a rare opportu-nity to develop growth – climate models using tree-ring data collected from all trees ( > 4 cm in diameter at breast height) within 50 × 50 m plots established in subalpine old-growth forests of western Canada. Our objective was to determine how tree growth responses to climate vary among tree size classes, species, and sites. We modeled relationships between times series of annual basal area increment ( Δ BA) and yearly climate variables for individual trees; this approach obviated key statistical criticisms of “ traditional ” tree-ring analysis methods. Time series of annual basal area increment were detrended a priori for size, age, legacy, and competition effects. We found that the overall climate signal in our time series of Δ BA was weak; < 6% of the interannual variance was explained by climate variables. Nevertheless, there were clear patterns in climate – growth relationships related to tree size and species. Relationships between Δ BA and ﬁ ve climate variables increased in strength with tree size class; large trees were most sensitive to annual climate ﬂ uctuations and accounted for ~ 71% of the overall climate effect on growth across all trees and sites. In all stands, Δ BA variance explained by climate variables was stable over the 20th century for large trees but decreased in the 1940s for small trees, indicating a temporal reduction in sensitivity to annual ﬂ uctuations in ﬁ ve climate variables. In coastal forests, relationships between Δ BA and climate for Callitropsis nootkatensis were signi ﬁ cantly different in direction and magnitude than those of co-occurring Pinaceae species. The effect of climate on tree growth was idiosyncratic among stands and could not be discriminated by forest type (coastal vs. interior). Our individual-tree modeling approach adds to a growing body of research providing novel insights about the complexities of tree growth responses to climate variability and the challenges associated with predicting future tree growth and forest productivity using tree-ring data.


INTRODUCTION
Global climate change is expected to alter tree growth around the world (Babst et al. 2019, Cailleret et al. 2019. Changes in growth affect ecological processes at multiple scales, including species interactions, rates and trajectories of forest succession, forest productivity, and the biogeochemical cycles that regulate landatmospheric feedbacks and the global climate.
Dynamic global vegetation models often predict substantial increases in forest productivity and have generated considerable debate about the capacity of rising atmospheric CO 2 concentrations to enhance forest growth and mitigate greenhouse gas emissions (Büntgen et al. 2019, Hararuk et al. 2019. These models frequently predict increases in forest productivity that is inconsistent with field observations of tree growth (Girardin et al. 2016, D'Orangeville et al. 2018, Giguère-Croteau et al. 2019, Hararuk et al. 2019. One reason for this is that large-scale model outcomes are often constrained by a lack of information about tree growth responses to climate variability at local scales (Babst et al. 2012, Huber et al. 2018, Klesse et al. 2018, Marchand et al. 2018. Rising global temperatures are expected to increase growth rates in cold-humid regions as cell development becomes more rapid and growing seasons lengthen (Rossi et al. 2014). This benefit of warming could be transitory though if thresholds for atmospheric water demand are crossed (D'Orangeville et al. 2018, Klesse et al. 2018. In regions that are already warm and dry, continued warming is expected to decrease growth-and ultimately kill trees-as water stress compromises wood production and hydraulic conductivity (Adams et al. 2017). However, general conclusions about the expected effect of climate change on tree growth are difficult to make because species and size (age) classes can respond differently, and these responses can be modified by factors that differ among sites, such as competition, nutrient availability, and herbivory (Babst et al. 2012, Marchand et al. 2018, Cailleret et al. 2019. Nonlinear, or threshold, growth responses to climate complicate the matter further (Lloyd et al. 2013, Buechling et al. 2017, Klesse et al. 2018). This complexity makes predictions about future tree growth and forest productivity challenging.
Tree-ring width chronologies offer an expedient way to examine how trees might respond to future changes in climate (Mann et al. 2002, Smith andLewis 2007) and numerous tree-ring studies report important modifications to growth related to climate variation (Vicente 2001, Martinelli 2004, Luckman and Scott 2007, Lapointe-Garant et al. 2010. However, several factors should be considered when using tree-ring studies to make inferences about climate influences on growth. Small sample sizes and sampling biases, such as a focus on the largest living trees in a stand or targeted collections from trees on growing on harsh sites, limit the scope of study conclusions (Nehrbass-Ahles et al. 2014, Büntgen et al. 2019, Duchesne et al. 2019, Trouillier et al. 2019. The outcomes of growth-climate analyses may be dependent upon the detrending method used to first remove non-climate influences on growth (i.e., size/age trends, or competition effects) (Esper et al. 2003, McShane and Wyner 2010, Peters et al. 2015, Schofield et al. 2016, Wood et al. 2016, Lee et al. 2017, Dietrich and Anand 2019. Traditional analytical approaches that estimate climate-related effects on tree growth from a "site-level" time series (obtained by averaging individual-tree times series of standardized values of annual growth) ignore statistical problems associated with averaging time series of individual trees that exhibit variance heteroscedasticity and different covariance structures (Esper et al. 2003, McShane and Wyner 2010, Schofield et al. 2016, Nothdurft and Vospernik 2018. These traditional approaches to modeling growth-climate relationships using tree rings have also rested on the assumption that tree growth remains sensitive to the same aspects of climate over time (Babst et al. 2019, Wilmking et al. 2020. New approaches to analyzing tree-ring data, as well as tree-ring data sets that better represent whole stands and a range of site conditions, could help overcome these challenges and generate more robust inferences about the potential effects of climate change on tree growth and forest productivity. In this study, we examined the influence of climate variability over the twentieth century on the annual basal area growth increment of trees in subalpine old-growth forests where short growing seasons are a constraint on forest productivity (Churkina andRunning 1998, Parish et al. 1999). We utilized a uniquely large set of tree-ring data (basal disks) obtained from most living and standing dead trees (>4 cm diameter at breast height, DBH) within plots established in old-growth forests that were intensively surveyed and then logged. Our primary objective was to understand subalpine tree growth sensitivity to climate variability. Given the findings of recent tree-ring studies in other forest environments (Babst et al. 2012, De Frenne et al. 2013, Harwood et al. 2014, Kelsey et al. 2018, Marchand et al. 2019, we hypothesized that annual tree growth responses to a set of key climate variables are dependent on tree size class, species, and site. To test this hypothesis, we modeled relationships between annual basal area increment and climate for individual trees-as opposed to modeling a mean time series of standardized growth indices-and then summarized individual model results by tree size class, species, and site. With this individual-tree approach, we bypass statistical problems associated with averaging time series with diverse trends, variance, and covariance structures (Schofield et al. 2016). Also, in contrast to most tree-ring studies, we evaluated the outcomes of different detrending methods before choosing a particular function to separate intrinsic growth effects from external climate effects, and we tested all individual-tree models for stationarity in growth-climate relationships. With large sample sizes and having addressed key challenges associated with traditional approaches to the analysis of tree-ring data, we expected an improved basis for understanding and predicting climate effects on tree growth in mesic subalpine old-growth forests of western North America.

Study sites
We used tree-ring data from four study sites in British Columbia (BC), Canada (Fig. 1); two in the Coastal Mountain Ranges-Mt. Cain (1000 m asl.) and Mt. Elphinstone (1100 m asl)-and two in the Columbia Mountain Ranges-Adams Lake (1900 m asl.) and Damfino Creek (1700 m asl). All sites represent forests characteristic of mesic edaphic conditions. Sites in the Coastal Mountains (coastal sites) have a maritime climate, whereas sites in the Columbia Mountains (interior sites) have a more continental climate with larger differences between mean winter and mean summer temperature (i.e., 16°vs. 23°C) (Wang et al. 2016) (Fig. 2). Mean annual temperature (MAT) at coastal sites is 4°C, about 3°C higher than at interior sites. Mean annual precipitation (MAP) exceeds 2400 mm/yr at coastal sites, which is more than double that at interior sites. Deep snowpacks accumulate at these sites and typically persist until late spring or early summer; average snow depth is~20% greater at coastal sites than interior sites and lowest at Damfino Creek. Growing-season climate moisture deficits, which are most intense in July, occur more frequently at interior sites than coast sites.
All four study sites represent old-growth, closed-canopy forests that have remained largely undisturbed for centuries. Typical of such forests, they were dominated by long-lived, shade-tolerant tree species and contained trees within most size and age classes (Antos and Parish 2002a, b, Parish and Antos 2004, 2006, Parish et al. 2008). The two interior sites were in high-elevation forests composed primarily of Abies lasiocarpa (Hook.) Nutt. and lesser amounts of Picea engelmannii Parry ex Engelm. The forest at Adams Lake was at least 462 yr old and contained no evidence of disturbance whereas the somewhat younger forest at Damfino Creek (~330 yr old) contained evidence of partial canopy disturbances caused by episodic spruce budworm outbreaks (Parish and Antos 2002). The two very old forests sampled in the coastal mountains contained four tree species: Abies amabilis (Dougl. ex Loud) Dougl. ex J. Forbes, Tsuga mertensiana (Bong.) Carrière, Tsuga heterophylla (Raf.) Sarg., and Callitropsis (Chamaecyparis) nootkatensis (D. Don) Spach, with trees exceeding 1300 yr in age Antos 2004, 2006).

Field sampling and tree-ring measurements
The old-growth forests sampled were typical of the region and scheduled for commercial clear-cut logging. We selected stands scheduled for logging so that basal disks could be collected after trees were felled. Prior to logging, four 50 × 50 m plots were established at each site and all living and standing dead trees ≥4.0 cm in diameter at breast height (DBH) were mapped and tagged. The DBH (cm) of all trees was measured, and then, basal disks were cut from trees too small to be of commercial value (generally trees <15 cm DBH). After logging, basal disks were cut as close to the ground as possible from the stumps of the larger trees. Logging occurred in 1996 at Adams Lake, in 1998 at Damfino Creek and Mt. Elphinstone, and in 2000 at Mt. Cain.
Ring widths on the basal disks were measured to the nearest 0.01 mm using a Measuchron digital positiometer (Micro-Measurement Technology, Bangor, Maine, USA). A radius that appeared representative of average growth, and avoided rot pockets, was used for the measurement of annual rings. A year was assigned to each tree ring and dating, or measurement errors were detected using COFECHA (Holmes 1983). Details of the ring measurement and sampling procedures are reported in Antos and Parish (2002a, b) and Antos (2004, 2006). Across all sites, dated chronologies of tree-ring widths were obtained from 3685 trees (living and dead) of varying sizes and ages (Table 1). We used chronologies from all four plots at Mt. Elphinstone and Adams Lake but only from two of the four plots at Mt. Cain and Damfino Creek because we found errors in the field records for the spatial locations of trees in these plots. This pool of tree-ring chronologies represented nearly 80% of all tagged trees in the plots (some tagged trees were lost during harvest operations). Almost all tree-ring width chronologies spanned the entire period used for building growth-climate models (i.e., from 1901 to the year before the stand was logged); trees that established or died during this period had shorter chronologies. Time series of basal area increment (cm 2 ) were calculated from these chronologies of tree-ring width.

A model for annual basal area increment
Given tree-to-tree variation in temporal growth trends can be substantial, particularly in old-growth forests Parish 2002b, Parish andAntos 2002), we modeled the basal area v www.esajournals.org increment of each tree separately. For any tree k in plot j, on site i, basal area increment between year t − 1 and t was cast as a function of last year's tree age (AGE), basal area (BA), and basal area increment (ΔBA) an indicator of competitive pressure (COMP), plus a function of (centered) climate variables (CLIM 1 to CLIM p ), and a residual error term (ϵ) as follows: The function f ijk is a "detrending function" to capture tree size-/age-and competition-related effects on annual basal area increment. Because tree size and age are not as highly correlated in old-growth forests as they are in even-aged stands-our stands contain many small shadetolerant understory trees that are very old (Antos et al. 2005)-we include both BA t and AGE t as terms in f ijk to account for both aspects of ontogenetic development on annual basal area increment, ΔBA ijkt (cf. Peters et al. 2015, Dietrich andAnand 2019). Growth in previous years may strongly influence current-year growth increment, and these legacy effects, often regulated by the abundance of stored carbon reserves (Hoch and Körner 2012), are represented by ΔBA tÀ1 in Eq 1. Competition for above-and below-ground resources may affect annual basal area increment in any year throughout the lifetime of a tree (Fritts 1976, Marchand et al. 2018, regardless of tree size or age. Thus, we also included a competition index, COMP t , in f ijk (details below). There are several ways the function f ijk could be formulated to extract size, age, and competition trends from times series of annual basal area increment. After comparing the nearly identical outcomes of four options (see Appendix S1), we selected a parametric function because of its transparency and to avoid poorly supported tree-specific decisions about frequency cutoff points for the extraction of non-parametric time series trends (Fritts 1976, Cook and Kairiukstis 1990, Magnussen and Alfaro 2012, Schofield et al. 2016, Nothdurft and Vospernik 2018. The selected parametric detrending function is where the parameters α 0 to α 7 were estimated by nonlinear least squares (Gallant 1987) with the observed basal area increments as the dependent variable. Initial search starting values for the parameters came from ordinary least squares estimates with a linear model on a log-log scale and a very small positive value added to each value of ΔBA to avoid taking the logarithm to 0. For each tree, the competition index (COMP) was computed for each year of observation as the inverse distance-weighted average basal area of its 6 (or fewer) nearest competitors within a circle of radius R12, which is the radius of a circle holding an average of 12 living trees (sensu Radtke et al. 2003, Rozas 2015. To calculate the v www.esajournals.org competition index, dead trees were given a BA of zero in the years when their status was dead. Given the detrending function, f ijk , the marginal effects of climate on ΔBA were assumed to follow a single linear model for all sites, plots, and trees such that where g ijk is a linear function of climate variables in year t and the climate variables lagged by one year. For each tree, parameter estimates β (current year) and γ (previous year) are the linear effect of selected climate variables l to p. We assume ϵ ijkt in (1) is a Gaussian random variable with an expectation of 0 and a weakly stationary (Daley 1971) parametric covariance function Γ ɛ ijks , ɛ ijkt jΘ À Á capturing both variance heteroscedasticity and temporal correlation among tree-specific residuals at time s and t (Θ denotes a vector of time-varying parameters, in the covariance function Γ).
We undertook a variable selection process to determine the climate variables to include in g ijk (see Appendix S2 for details). Yearly climate data were obtained from ClimateBC (Wang et al. 2016, v.5.4) because no continuous records from 1901 to 2000 were available from weather stations near our study sites. Based on published records about the influences of climate on tree growth (Fritts 1976, Briffa et al. 2004, Martinelli 2004, Speer 2010, we initially selected data for 22 of the 60 climate variables calculated by Cli-mateBC. Nine were related to monthly or seasonal temperature (°C), three to growing degreedays (degree-days > 5°C), four to precipitation (mm), four to measures of evapotranspiration (mm) (Hargreaves 1974), and two to solar radiation (MJÁm −2 Ád −1 ). Using a randomly selected subset of time series (N = 653 trees), we constructed regression models to test the marginal association between detrended basal area increment and climate variable. Each of the 22 candidate climate variables was tested separately for each of the 653 trees. A climate variable was determined to have a statistically significant effect if the probability of a likelihood ratio test statistic under the null hypothesis was ≤0.05 (adjusted for family-wise error using Holm's sequential rejection test, Holm 1979). Climate variables that had the highest frequency of statistically significant effects across all 653 × 22 individual-tree models were considered the most influential climate variables. An additional criterion for retaining a variable was that no pairwise correlation among the retained variables exceeded 0.7. The five climate variables retained for inclusion in our model were as follows: (1) average July temperature (JULT,°C); (2) minimum winter (December-February) temperature (MINWT,°C); (3) mean summer (May-September) precipitation (MSP, mm); (4) Hargreaves reference evaporation during the seasons of spring, summer, and autumn (EREF, mm); and (5) precipitation as snow in winter prior to a growing season (PAS, mm).
Estimation of climate effects (g ijk ), in the detrended empirical residuals (ΔBA ijkt Àf ijk ), was done with generalized least squares (GLS) for several tree-level candidate time series models, including commonly used autoregressive processes (AR1, AR2, . . . etc.) (Cook and Kairiukstis 1990), to capture temporal autocorrelation structures. The maximum likelihood of candidate time series models was computed, and we selected the one with the lowest corrected AIC (Fan and Yao 2005); a diverse array of autocorrelation functions was selected, but low-order AR processes were dominant. Following an estimation of the regression coefficients in the climate model, we conducted tests on the residuals for weak stationarity and temporal autocorrelation (Ljung-Box test) to determine whether GLS residuals were independently and identically distributed. Almost all individual-tree models (>95%) passed the test of weak stationarity, and most (87%) had no statistically significant temporal autocorrelation.
The models developed for nearly 3700 trees provided information about growth-climate response variability among sites, tree size classes, and species. Our estimated generalized least squares regression coefficients for each tree model correspond to mean-centered values of JULT, MINWT, MSP, EREF, and PAS, and their one-year lagged effects (JULT -1 , MINWT -1 , MSP -1 , EREF -1 , and PAS -1 ). To interpret variability in the annual detrended basal area increment due to yearly climate variable fluctuations, we calculated mean regression coefficients (and standard errors) for each climate variable for trees grouped by site, by size class, and by species. Trees were assigned to one of six relative size classes by dividing the basal area of a tree by the basal area of the largest tree in a plot, multiplying it by six, and then rounding this value upwards to the nearest integer. The null hypothesis of no significant differences in mean regression coefficients among size classes and sites was tested with a least significant test criterion at the 5% level of significance (Snedecor and Cochran 1971). Hotelling T 2 permutation tests (Good 1993), with 2000 permutations of species labels, were used to test null hypotheses of no differences in pairwise species comparisons of multivariate mean regression coefficients for the five climate variables and their one-year lags considered jointly. Species comparisons indicating a significant difference in mean regression coefficients were followed by a Student's t test of differences in each of the set of two mean regression coefficients linked to a climate variable. To explore temporal variability in growth-climate relationships, we plotted time series of estimated annual basal area increment attributed to the five climate variables (and their one-year lags) over all trees, by tree size class and by species. The variance attributed to the climate effects in detrended growth residuals was taken as an indicator of climate signal strength in a time series.

Detrending time series of basal area increment
We used a parametric function, f ijk , to extract the size, age, and competition trends from each time series of annual basal area increment (see Appendix S1 for comparisons of detrending methods). The detrending function explained a large percentage of the annual variability in basal area increment at all sites (Fig. 3). At coastal sites, it accounted for 66% of the annual variability in basal area increment at Mt. Cain and 69% at Mt. Elphinstone. At interior sites, the detrending function accounted for more of the annual variability in annual growth increment: 77%.

Climate variable effects on detrended basal area increment
We constructed 3685 models to explore the effects of yearly climate variable fluctuations on detrended annual basal area increment. Across all these models, the five climate variables and their one-year lags explained, on average, <6% of the interannual variance in detrended annual basal area increment (hereafter, referred to simply as "annual basal area increment"), suggesting a relatively weak climate signal in the tree rings at each site (Table 2). They accounted for 3% of the variance in basal area increment at all sites except Mt. Cain where climate explained about 6% of variance in basal area increment. Near maximum (90th percentile) variance explained by climate ranged from 10% to 66% across all sites. The climate variables had a significant effect on basal area increment in 8 to 10% of trees, except at Damfino Creek where this percentage was lower.
Annual basal area increment was most strongly influenced by yearly fluctuations in JULT and MINWT with the mean model regression coefficients for these climate variables being orders of magnitude larger than those for MSP, EREF, or PAS (Appendix S3: Table S1). While JULT and MINWT had the greatest influences on basal area increment overall, their effects varied across sites (Fig. 4). JULT (in the current and previous year) had the greatest influence on annual basal area increment at Mt. Cain and Damfino Creek, with above-average JULT being associated with increased growth increment. JULT had smaller effects on tree growth at Mt. Elphinstone and Adams Lake, with above-average JULT being associated with below-average basal area increment. JULT −1 had no significant effect on growth at these sites. MINWT had the strongest effect on tree growth at Mt. Cain. At all sites except Damfino Creek, higher than average MINWT was associated with below-average basal area increment. MINWT −1 had significant, but opposite effects on annual basal area increment at Mt. Elphinstone and Adams Lake.
The annual basal area increment of large trees was most sensitive to yearly fluctuations in influential climate variables (Fig. 5). At all sites, mean regression coefficients for JULT, MINWT, and their one-year lagged effects were much larger for trees in BA rank classes 5 and 6 than for trees in smaller BA rank classes. Large trees accounted for most (mean = 71%, range = 60-86%) of the total climate effect on annual basal area increment (i.e., the sum of regression coefficients for BA rank classes 5 and 6 divided by the sum of regression coefficients over all classes; Fig. 5; v www.esajournals.org Appendix S3: Table S2). Significant mean responses (i.e., mean regression coefficient/standard error >2) to yearly climate variable fluctuations were observed in several size classes at Mt. Elphinstone, Adams Lake, and Damfino Creek, but only in the largest size classes at Mt. Cain.
Using model estimates of annual basal area increment from 1900 to 2000, we found pronounced size-dependent temporal trends in growth responses to climate variables. Variance in the proportion of annual basal area increment explained by climate variables was greatest early in the century and decreased substantially by the middle of the century for trees in basal area rank classes 1 and 2 at all sites (Figs. 6, 7) and basal area rank classes 3 and 4 at Mt. Elphinstone (Fig. 6). The mean also decreased at Mt. Cain and Mt. Elphinstone (Fig. 6). The mean proportion of annual basal area increment explained by the climate variables (as well as the variance) was comparatively constant in the larger basal area rank classes (Figs. 6, 7). Among the 35 pairwise species comparisons undertaken with Hotelling T 2 permutation tests for equality in multivariate regression coefficients, 13 indicated statistically significant differences in species growth responses to yearly climate variable fluctuations; they all occurred at coastal sites (Mt. Cain and Mt. Elphinstone), and 11 of them involved C. nootkatensis (Table 3; Appendix S3:  Table S3). These differences occurred over all basal area classes at Mt. Elphinstone but were limited to the largest trees at Mt. Cain (Table 3). At Mt. Cain, JULT (and its lagged effect, JULT −1 ) generated the greatest species differences in annual basal area increment of large trees (BA rank class 6); above-average JULT was associated with basal area increment decreases in C. nootkatensis but increases in A. amabilis and Tsuga spp. (Table 3). At Mt. Elphinstone, significant differences between C. nootkatensis and A. amabilis were predominantly due to their growth responses to MINWT -1 ; both species had a negative relationship with MINWT -1 but the response of C. nootkatensis was much stronger. Response differences between C. nootkatensis and T. mertensiana occurred among BA rank classes 3, 4, and 6 and were strongly influenced by JULT −1 and MINWT −1 . The basal area increments of A. amabilis and T. mertensiana had significant negative relationships with MINWT −1 but the effect of aboveaverage winter temperature was more detrimental to the growth of A. amabilis.
Temporally varying differences in coastal species sensitivity to climate variability were also evident. For example, at Mt. Elphinstone, where differences were most pronounced (Appendix S3: Table S3), the influence of climate variables on annual basal area increment varied over the 20th century and the temporal patterns in this variability were different for C. nootkatensis and A. amabilis (Fig. 8). During years when climate variables had above-average influences on annual basal area increment of A. amabilis (i.e.,ΔbBA CLIM > 0), the same variables often contributed less to the annual increment of C. nootkatensis and vice versa. Between 1920 and 1940, the contribution of climate variables to basal area increment was below average for C. nootkatensis and above average for . ‡ Variance at the 90th percentile of a distribution of percentage variance explained by CLIM l-p for each model. § No. of signif. = the number of models with significant CLIM l-p effects. ¶ ϵ ijkt is the residual variance after accounting for CLIM l-p effect son detrended basal area increments. Fig. 4. Mean regression coefficients, β (multiplied by 10 4 ), for influential climate variables-July temperature (JULT,°C), and its one-year lag effect (JULT −1 ), as well as, minimum winter temperature (MINWT,°C) and its one-year lag effect (MINWT −1 )-averaged across all trees at each site. Asterisks indicate significant mean regression coefficients, where mean β/mean standard error >2.
A. amabilis. The opposite was true for the periods 1945-1950 and 1955-1965, during which climate variables had more influence on the basal area increment of C. nootkatensis than on the growth increment of A. amabilis.

DISCUSSION
Given different detrending approaches could generate disparate findings about climate influences on basal area increment (Schofield et al. Fig. 5. Mean regression coefficients, β (multiplied by 10 4 ), for climate variables-July temperature (JULT,°C ) with its one-year lag effect (JULT −1 ), and minimum winter temperature (MINWT,°C) with its one-year lag effect (MINWT −1 ). Results are by relative tree size class (i.e., relative basal area rank class, 1 = smallest trees, and 6 = largest trees in a stand). Coefficients for each climate variable from each individual-tree model were averaged across all models within each of six size classes. Asterisks indicate significant mean regression coefficients, where mean β/mean standard error >2. v www.esajournals.org Estimates were averaged each year, over all individual-tree models grouped by tree size class. To simplify the graph, we have pooled results from classes 1 and 2, 3 and 4, and 5 and 6. Rank classes 1 and 2 represent the smallest trees in a stand, while classes 5 and 6 represent the largest trees. The mean age of trees in a given year (in brackets along the x-axis) indicates a gradient in tree age. The gray areas represent the upper and lower limits of 95% CI. Estimates were averaged each year, over all individual-tree models grouped by tree size class. To simplify the graph, we have pooled results from classes 1 and 2, 3 and 4, and 5 and 6. Rank classes 1 and 2 represent the smallest trees in a stand, while classes 5 and 6 represent the largest trees. The mean age of trees in a given year (in brackets along the x-axis) indicates a gradient in tree age. The gray areas represent the upper and lower limits of 95% CI. 2016), we compared the outcomes of a parametric detrending method to three non-parametric detrending alternatives: a Gaussian smoother (Wood et al. 2016), a semi-parametric truncated power series (Ruppert et al. 2003), and a smoothing function via Wiener filtering (Wiener 1949). We found no significant differences in the resulting residual time series (Appendix S1, Appendix S2: Table S1) and used the parametric method to detrend time series for subsequent analyses of climate effects on basal area increment. This detrending method had the advantage of transparency; it avoids heuristics when choosing frequency response cutoffs (i.e., bandwidth) for commonly used non-parametric smoothing splines (Cook and Holmes 1986, Cook and Kairiukstis 1990, Bunn et al. 2019) and enabled intuitive summaries of climate effects (by tree size, species, and site) using parametric model statistics. While we sought the benefits of a parametric detrending function in this study, we acknowledge that non-parametric methods are more appropriate for applications with data from even-aged stands.
We constructed 3685 parametric models quantifying the effect climate variables on detrended annual basal area increment. Across all models, climate variables explained, on Notes: ABAM, Abies amabilis; CHNO, Callotropsis nootkatensis; TSHE, Tsuga heterophylla; TSME, Tsuga mertensiana. † A Hotelling's T 2 permutation test statistic was used to identify significant species differences in multivariate mean regression means (see Appendix S3: Table S3 for all 35 species comparison tests). Pairwise comparisons of mean regression coefficients for the single most influential climate variable are presented here (see Appendix S3: Table S4 for selection of these climate variables).
‡ Number of individual-tree models for each species used in comparisons, regression coefficients for each climate variable are averaged over all models, by species. Sp1 = the first species in the pairwise comparison, Sp2 = the second species in the pairwise comparison.
§ Regression coefficients are multiplied by 10 4 . ¶ Significant differences in mean regression coefficients (for each climate variable) were tested using Student's t test. average, only~4% of the annual variation in detrended basal area increment of the subalpine trees we sampled in western Canada. Moreover, 10%, or less, of the trees at each site exhibited a statistically significant growth response to climate variables. This small climate signal in the tree-ring chronologies may seem surprising given the purported climate sensitivity of trees in high-altitude forests (Pepin et al. 2015, Seddon et al. 2016, and the obvious importance of climate in controlling the growth, abundance, and distribution of species (Churkina and Running 1998). Dendro-climate researchers have long recognized that the strongest relationships between radial growth and climate occur in open forests at the cold and dry margins of tree species' geographic distributions (Fritts 1976). None of the species we examined in closedcanopy forests were near their geographic or edaphic distributional limits, where we might expect them to be most sensitive to interannual fluctuations in climate parameters (Buechling et al. 2017, Klesse et al. 2018. However, other tree-ring studies in closed-canopy subalpine and coastal montane forests of western North America report somewhat larger climate influences on interannual growth variability (Buechling et al. 2017, Kelsey et al. 2018, Legendre-Fixx et al. 2018). Such differences in climate influences among studies could be due to differences in site conditions, species studied, and various methodological factors. Our tree-ring samples were collected near the base of harvested trees rather than at breast height (i.e., 1.3 m), or higher, which might contribute to lower sensitivity to climate (Chhin andWang 2005, Millar et al. 2018). Detrending to extract the influences of allometry (size), age, and competition could have removed some climate effects on growth increment (Fritts 1976). Also, growth-climate relationships may be confounded by interannual variability in the storage and remobilization of carbon, including the reallocation of carbon to reproductive tissues instead of radial growth (i.e., seed masting years) , or the use of stored reserves during environmentally adverse years (Litton et al. 2007, Hoch and Körner 2012, Sala et al. 2012, Locossellii and Buckeridge 2017. However, these latter two challenges are common to all tree-ring studies.
We found pronounced and predictable relationships between annual basal area increment and climate, even though the overall climate signal in our tree-ring series was weak. Among the five climate variables examined, JULT and MINWT had the strongest influences on basal area increment at all sites; many tree-ring studies in temperate forests (using various methods) find these two aspects of climate to be key drivers of tree growth (e.g., Biondi et al. 1999, Briffa et al. 2002, McCoullough et al. 2017. Given the findings of other North American studies in subalpine forests (Buechling et al. 2017, Legendre-Fixx et al. 2018, we anticipated that warmer-than-average current-year JULT would increase basal area increment, but at two of our four sites (one coastal and one interior site), we found a negative relationship between basal area increment and JULT. Many subalpine tree-ring studies report reductions in annual growth on warm, dry sites due to summer moisture stress (e.g., Peterson et al. 2002, Büntgen et al. 2006, Nakawatase and Peterson 2006, Kelsey et al. 2018 but moisture stress can also reduce growth on some mesic sites (Splechtna et al. 2000). Summer moisture stress might account for the negative relationship we found at the interior site, but this seems a less likely explanation in the wetter coastal climate of Mt. Elphinstone. Above-average MINWT had a positive effect on annual basal area increment at Damfino Creek but a negative effect on growth increment at the other three sites. The benefits of mild winters on subalpine tree growth are commonly documented and typically linked to longer growing seasons (Vitas 2006, Millar et al. 2018. The negative effects of mild winters reported in tree-ring studies are often attributed to the effects that reduced winter snow accumulation and early spring snowmelt have on growing-season soil moisture availability (Sanmiguel-Vallelado et al. 2019, Harley et al. 2020. The relatively weak relationships we found between growth increment and PAS imply other factors play a more important role. Some studies reporting negative relationships between winter temperature and annual growth suggest the rapid consumption of stored carbohydrates needed for early reactivation of cambial growth may be a limiting factor (Sala et al. 2012, Locosselli and Buckeridge 2017, Nothdurft and Vospernik 2018. The largest trees in a stand were most sensitive to year-to-year climate variable fluctuations and the strength of climate effects on annual basal area increment progressively decreased with tree size (i.e., basal area rank class). We observed this trend at all sites and even when we grouped all species together, which suggests a degree of generality regarding size-related climate sensitivity. This generality is supported by another study in a different subalpine forest type (Carrer and Urbinati 2004) but in contrast to others reporting opposite trends, or no sizerelated trends at all (Esper et al. 2008, Vieira et al. 2009). Size-related trends in tree growth sensitivity to climate variability can be attributed to physiological processes mediated by tree size (Mencuccini et al. 2005). Limitations on the upward transport of water and nutrients-that are imposed by tree height (Bond 2000)-may make large trees more sensitive than small trees to climatic changes. This could account for sizerelated differences in sensitivity on sites where we found negative relationships between temperature and growth, but it is a less likely explanation for differences in sensitivity where we observed positive temperature-growth relationships. On these sites, the benefits of warmer temperatures (longer growing seasons and increased photosynthetic rates) may disproportionately benefit large trees that are in more direct contact with the macroclimatic environment whereas smaller trees are buffered by cooler microclimates created by a forest canopy (De Frenne et al. 2013, Harwood et al. 2014. Regardless of the eco-physiological causes, patterns in size-dependent climate sensitivity documented in the literature are highly variable, and this topic remains an active area of tree-ring research (Mérian and Lebourgeois 2011, Wunder et al. 2013, Trouillier et al. 2019).
For small trees, the proportion of annual basal area increment explained by climate was highly variable at the beginning of the 20th century and this variability in growth response decreased substantially at all sites in the latter half of the century. At some sites, the mean trend in annual increment explained by climate also decreased. This implies that small trees became progressively less sensitive to climate over time.
Developmental-or age-related decreases in climate sensitivity are reported in other tree-ring studies (Dorado-Liñán et al. 2011, Konter et al. 2016. Many functional and physiological processes that change over a tree's lifespan can influence growth sensitivity to climate (Bond 2000, Thomas 2011, Rossi et al. 2014, and young trees may exhibit greater physiological plasticity to climate fluctuations than older trees (Delagrange et al. 2004, Fischer et al. 2014. Age-and size-dependent influences on tree sensitivity to climatic variability are often difficult to disentangle in tree-ring research (Konter et al. 2016). The results of our study in uneven-aged old-growth forests suggest both age and size play a role in determining tree sensitivity to climate variability. Whereas tree size effects on climate sensitivity are related to hydraulic limitations imposed by tree height, age-related differences in climate sensitivity are possibly the outcome of functional and physiological processes under genetic control (Greenwood 1995, Mencuccni et al. 2005. Other plausible explanations for the decreasing temporal trend in climate sensitivity of small trees include an early-century landscape-scale disturbance that opened forest canopies and temporarily exposed small trees to macroclimates normally experienced by canopy trees, or the presence of compression wood (e.g., from mechanical stress due to snow loads) that can generate considerable noise near the beginning of tree-ring series. Artifacts introduced by detrending may also have played a role. We have no evidence of synchronous landscape disturbances at our study sites (Antos and Parish 2002a, Parish and Antos, 2002, 2006, and a recent study indicates compression wood may have little effect on climate signals in tree rings (Janecka et al. 2020). However, we cannot rule out the possibility that the decreasing trend in small-tree sensitivity is an artifact of our detrending method, particularly if our competition index did not adequately represent competition intensity early in the century (i.e., we were unable to document early-century competitors if they had fallen, died, and decayed before the stands were sampled).
Differential species responses to climate variables were pronounced at the coastal sites, primarily due to the distinct growth-climate relationships of C. nootkatensis. Trait differences v www.esajournals.org among C. nootkatensis and co-occurring species could explain differential species growth responses. For example, the indeterminate leaf growth (i.e., lack of overwintering buds) may give C. nootkatensis greater flexibility to respond to favorable climate conditions than co-occurring species. While this flexibility may not translate directly into basal area increment, it may affect carbon allocation dynamics (Locosselli and Buckeridge 2017) and explain why we found temperatures in previous years were the primary drivers of species differences. The shallow root systems of C. nootkatensis can also make it more sensitive than co-occurring species to cold soils. Root exposure to lethal freezing events, due to recent decreases in snow cover, is the primary explanation for the recent decline of this species in forests further north (Buma et al. 2017). At Mt. Elphinstone, where species differences were most pronounced, climate variables often contributed most to the annual increment of C. nootkatensis during years when they contributed least to the annual growth of A. amabilis. It is difficult to determine the reasons for these temporal differences in sensitivity. Between 1920 and 1940, when JULT and MINWT temperatures were warmer than average, climate had more of an influence on the annual basal area increment of A. amabilis than it had on C. nootkatensis but during a slightly warmer period later in the century (around 1960), the opposite was true. Changes in the aspects of climate affecting growth and nonlinear threshold responses to the climate variables might explain this finding. Unexpected growth responses to climate found in other studies were attributed to interactions between interannual fluctuations in temperature and precipitation and decadal climate oscillations (Linderholm 2002, Millar et al. 2018. In contrast to the coastal sites, tree species at the interior sites, P. engelmannii and A. lasiocarpa, responded similarly to interannual climate variable fluctuations. Other studies report divergent growth responses to climate variability at the southern limits of P. engelmannii-A. lasiocarpa forests and attribute this to differential species sensitivity to drought (Buechling et al. 2017, Kelsey et al. 2018.
We anticipated that subalpine tree growth increment would respond similarly to interannual climate variation within each of the two climatic regions (i.e., coastal and interior forests). However, both the magnitude and direction of growth-climate relationships were site-specific within each region. Several other tree-ring studies report variable growth responses among sites and attribute this to spatial differences in regional climate and local site factors, such as topoedaphic conditions, underlying geology, stand structure, species composition, competition intensity, and genotype (Babst et al. 2012, Buechling et al. 2017, Kelsey et al. 2018, Legendre-Fixx et al. 2018, Millar et al. 2018, Nothdurft and Vospernik 2018, Marchand et al. 2019, Sáenz-Romero et al. 2019. A greater abundance of C. nootkatensis at Mt. Elphinstone likely contributes to differences in stand-level responses at the coastal sites. However, at the interior sites, stand species composition and structure (i.e., tree size-class distributions) are virtually identical (Antos and Parish 2002a, b), suggesting other factors, such as regional climate differences, are more important in shaping site-specific differences in growth responses to climate variability. Regardless of underlying causes, a highly differentiated sitespecific response suggests that local populations are likely to respond idiosyncratically to future climate change (Chen et al. 2010, Buechling et al. 2017. Most tree-ring studies assume stationary growth-climate relationships despite a long history demonstrating they may vary over time (Briffa et al. 1998, Wilmking et al. 2005, Carrer and Urbinati 2006, D'Arrigo et al. 2008, Visser et al. 2010, Leonelli et al. 2011, Lloyd et al. 2013, Babst et al. 2019, Hofgaard et al. 2019. Wilmking et al. (2020) suggest that non-stationarity likely characterizes the fundamental nature of long-term relationships between tree growth and climate. In our study, we did not see strong evidence for temporally varying regression coefficients. This, in addition weak stationarity tests on model residuals, suggests growth-climate relationships were relatively stable over the period of analysis. Nevertheless, we agree with Wilmking et al. (2020) that failing to account for non-stationarity impedes our capacity to accurately estimate the effects of climate variability on tree growth. Regression analyses that account for temporally varying regression coefficients (Tibshirani et al. 2005, Stock and Watson 2019) could improve the v www.esajournals.org predictive capacity of models when non-stationarity in growth-climate relationships occur.
Tree-ring studies like ours, characterized by large number of samples from a wide range of tree sizes and ages, are rare. Both the fieldwork to obtain these samples and laboratory work preparing and measuring tree rings are time consuming, costly, and impractical to implement on a wide scale. Existing collections of tree-ring data that are spatially extensive-such as those from regional forest monitoring plots or the International Tree-Ring Data Bank-are often comprised of a small number of samples from the largest trees in a sampled population. Using these data to make inferences about forest-level sensitivity to climate variability could be misleading (Nehrbass-Ahles et al. 2014, Klesse et al. 2018, Duchesne et al. 2019. We found large among-tree variability in basal area increment responses to interannual climate variable fluctuations-including many trees that exhibited no significant response-emphasizing a need for relatively large sample sizes that represent the range of tree sizes, and species in a study population. In our uneven-aged old-growth forests, we determined that a sample size of about 60 trees would suffice to protect us with a probability of 0.9 (power) to accept the null hypothesis of no climate effects when it is false (see Appendix S4); similar sample sizes have been advocated for historical treering reconstructions of stand growth . While tree growth responses to climate were highly variable, we also found that the growth increment of the largest trees in a population was most sensitive to yearly climate variable fluctuations, accounting for about 70% of the total effect of climate on growth across all trees in a population. Because large trees were most sensitive to climate variability and represent a disproportional amount of stand biomass, it is possible that tree-ring data from the largest trees in our old-growth subalpine forests are sufficient to estimate the effects of climate variability on stand biomass and suitable for use in large-scale models of forest productivity. But this may not be the case for all forests (Esper et al. 2008, Vieira et al. 2009). Because the potential impacts of a large-tree bias on estimates of stand biomass sensitivity to climate variability may be unknown until tree size-dependent differences in climate responses have been quantified, the risk of limiting a sample to larger trees in a population should not be taken lightly.
Our study demonstrated highly diverse subalpine tree growth responses to interannual climate variable fluctuations. This finding adds to a growing body of tree-ring research providing novel insights about the complexity of tree growth responses to climate variability (Girardin et al. 2016, Buechling et al. 2017, D'Orangeville et al. 2018, Housset et al. 2018, Klesse et al. 2018, Marchand et al. 2018, Millar et al. 2018, Nothdurft and Vospernik 2018, Babst et al. 2019, Büntgen et al. 2019, Cailleret et al. 2019, Correa-Díaz et al. 2019, Hararuk et al. 2019, Sáenz-Romero et al. 2019. Interactions among ecological factors -tree size and life stage, species traits, and site conditions-and interannual fluctuations in temperature and precipitation, as well as decadal/multidecadal climate variability, can influence tree growth responses, sometimes in unexpected ways (Linderholm 2002, Millar et al. 2018). Further complexity is imposed by directional changes in the global climate (Babst et al. 2019, Hofgaard et al. 2019. Both JULT and MINWT have increased at our study sites since the tree-ring data were collected, and as climate regimes move outside the envelope of past conditions, the capacity of retrospective models to accurately predict future tree growth is likely to decrease. New tree-ring data collections that sufficiently capture variability in recent growth, and iterative model improvements that consider the implications of different detrending options and changing growth-climate relationships, will be necessary when we strive to predict the effects of climate on future tree growth.