Joint effects of climate, tree size, and year on annual tree growth derived from tree-ring records of ten globally distributed forests.

Tree rings provide an invaluable long-term record for understanding how climate and other drivers shape tree growth and forest productivity. However, conventional tree-ring analysis methods were not designed to simultaneously test effects of climate, tree size, and other drivers on individual growth. This has limited the potential to test ecologically relevant hypotheses on tree growth sensitivity to environmental drivers and their interactions with tree size. Here, we develop and apply a new method to simultaneously model nonlinear effects of primary climate drivers, reconstructed tree diameter at breast height (DBH), and calendar year in generalized least squares models that account for the temporal autocorrelation inherent to each individual tree's growth. We analyze data from 3811 trees representing 40 species at 10 globally distributed sites, showing that precipitation, temperature, DBH, and calendar year have additively, and often interactively, influenced annual growth over the past 120 years. Growth responses were predominantly positive to precipitation (usually over ≥3-month seasonal windows) and negative to temperature (usually maximum temperature, over ≤3-month seasonal windows), with concave-down responses in 63% of relationships. Climate sensitivity commonly varied with DBH (45% of cases tested), with larger trees usually more sensitive. Trends in ring width at small DBH were linked to the light environment under which trees established, but basal area or biomass increments consistently reached maxima at intermediate DBH. Accounting for climate and DBH, growth rate declined over time for 92% of species in secondary or disturbed stands, whereas growth trends were mixed in older forests. These trends were largely attributable to stand dynamics as cohorts and stands age, which remain challenging to disentangle from global change drivers. By providing a parsimonious approach for characterizing multiple interacting drivers of tree growth, our method reveals a more complete picture of the factors influencing growth than has previously been possible.

INTRODUCTION

Tree rings provide a long‐term record of annual growth increments that is invaluable for understanding forests in an era of global change (Amoroso et al., 2017; Fritts & Swetnam, 1989; Zuidema et al., 2013). Spanning time scales of decades to centuries or even millennia, they provide by far the most robust method for characterizing the interannual climate sensitivity of tree growth (Bräker, 2002; Fritts, 1976). Combined with forest censuses, they can be used to estimate forest woody productivity (Davis et al., 2009; Dye et al., 2016; Graumlich et al., 1989) and its climate sensitivity (Helcoski et al., 2019; Klesse et al., 2018; Teets et al., 2018). Tree rings also provide the long‐term perspective necessary for understanding how slowly changing environmental drivers including rising atmospheric carbon dioxide (CO2) concentrations, changing climate, and other anthropogenic and natural changes are influencing tree growth and forest productivity (e.g., Levesque et al., 2017; Mathias & Thomas, 2018; Walker et al., 2021). This information is critical to predicting forest responses to anthropogenic changes—particularly climate change—and thereby reducing the large uncertainty surrounding future contributions of Earth's forests to the global carbon cycle (Arora et al., 2020). Yet, collection and analysis of dendrochronological records has traditionally been optimized to detect climate signals rather than to understand variation among trees—including size‐related variation in climate sensitivity (e.g., Bennett et al., 2015; McGregor et al., 2021; Rollinson et al., 2021)—or to predict forest productivity, its climate sensitivity, and how it may be changing (Babst et al., 2018; Cherubini et al., 1998; Klesse et al., 2018; Nehrbass‐Ahles et al., 2014; Wilmking et al., 2020). As a result, prevailing approaches hold a number of limitations for using tree rings to address pressing questions concerning forest carbon sequestration in the current era of rapid environmental change.

To realistically estimate forest woody productivity, it is necessary to measure or model the growth rate of individual trees within a stand based on the primary biotic and abiotic drivers. Specifically, what is needed is an analysis framework that can capture the additive and interactive effects of climate, tree size (most commonly diameter at breast height, DBH), and other environmental drivers, all of which may be best described by nonlinear functions (Muller‐Landau et al., 2006; Rollinson et al., 2021). Although multifactorial and sometimes nonlinear individual‐based analysis frameworks have been applied in tree‐ring analysis (e.g., Evans et al., 2017; Klesse et al., 2020; Rollinson et al., 2021; Zuidema et al., 2020), their use has been relatively limited, and none simultaneously account for climate, DBH, and calendar year (a proxy for slowly changing environmental drivers). Below, we outline major hypotheses regarding the influence of these factors on tree growth that can best be addressed using a multifactorial, nonlinear approach to tree‐ring analysis (Table 1). We then develop such a framework and apply it to test these hypotheses.

TABLE 1

Summary of hypotheses and specific predictions tested using the method developed here, along with the frequency at which they were supported in our analyses of tree‐ring data from ten globally distributed forests

Hypotheses and specific predictions	Frequency observed
Interannual climate variation a
Drought limits growth, but water responses are nonlinear.
Growth responds positively to water,	93% (42/45 SSC)
…but positive responses decelerate or decline at high precipitation.	76% (32/42 SSC)
High temperatures (T) limit growth, often nonlinearly.
Growth responses to T are predominantly either negative…	29% (13/45 SSC)
…or non‐linear concave down.	40% (18/45 SSC)
However, there are cases where growth increases with T.	15% (7/45 SSC)
Climate sensitivity often varies with tree diameter at breast height (DBH).
Water and DBH have an interactive effect on growth.	44% (16/36 SSC) b
Temperature and DBH have an interactive effect on growth.	38% (12/32 SSC) b
Diameter (DBH) c
DBH—ring width (RW) relationships depend upon the light environment.
RW declines with DBH for light‐demanding species,	46% (6/13 SSC)
…but initially increases with DBH for shade‐tolerant species.	73% (8/11 SSC) d
Basal area and biomass increments reach maxima at intermediate DBH.
Basal area increment (BAI) peaks at intermediate DBH.	95% (41/43 SSC)
Biomass increment (∆AGB) peaks at intermediate DBH.	98% (42/43 SSC)
Calendar year e
Size‐corrected growth rates decline with time since severe disturbance.
In secondary or disturbed forests, growth rates of most species have declined.	92% (23/25 sp. at 7 sites)
In forests dominated by >100‐yr‐old trees, growth trends are mixed.
In older forests, growth rates of some species have declined,	50% (6/12 sp. at 3 sites)
….whereas others have increased.	25% (3/12 sp. at 3 sites)

Abbreviation: SSC, species–site combination.

Results summarized here are for climate‐only models with RW as response variable.

Refers to SSC with significant (p<0.05) main effect of climate on RW.

Results summarized here are for models without year.

100% (9/9 SSC) for models including year.

Results summarized here are for models with BAI as response variable.

Key hypotheses on tree growth

Understanding the climate sensitivity of tree growth is critical to predicting forest dynamics and productivity as the climate changes. The classic dendrochronological approach to characterizing the climate sensitivity of tree growth describes linear relationships between the primary growth‐limiting climate factor (moisture or temperature) and population‐level growth responses captured in ring‐width index chronologies (Cook, 1985; Fritts, 1976; Speer, 2010). While invaluable for applications such as reconstructing past climates (e.g., Buntgen et al., 2011), accurately representing the climate sensitivity of forest productivity requires an analysis of sensitivity to multiple climatic variables and how this varies across trees in a forest stand (Babst et al., 2018). Precipitation and temperature can have additive or interactive effects on growth (Foster et al., 2016; Meko et al., 2011; Sánchez‐Salguero et al., 2015; Vlam et al., 2014; Zuidema et al., 2020), and we hypothesize that both influence the growth of most species, often over different seasonal windows (Table 1). In addition, we hypothesize that nonlinear climate responses, already known to be widespread within forest settings (Rollinson et al., 2021; Wilmking et al., 2020; Woodhouse, 1999), are in fact the predominant form of response in forests around the world (Table 1). Finally, the influence of diameter at breast height (DBH) is typically removed through detrending (Cook & Peters, 1997), which eliminates the potential to directly model its influence on climate sensitivity, but we hypothesize that interactive effects of DBH and climate are, in fact, quite common in forest settings and fundamental to understanding climate change responses of forests (Table 1; Bennett et al., 2015; McGregor et al., 2021; Rollinson et al., 2021; Trouillier et al., 2019).

Tree DBH scales with numerous traits affecting tree growth (e.g., height, crown size and position, root mass, hydraulic architecture) and, therefore, is itself linked to growth (e.g., Anderson‐Teixeira, McGarvey, et al., 2015). Yet, there remain inconsistencies across disciplines (dendrochronology, forest ecology) as to what is considered a typical growth pattern. Dendrochronological records from trees that established in high‐light environments commonly show a pattern in which radial stem growth increments (ring width, RW) are initially large and decline with DBH (Fritts, 1976), whereas censuses of globally distributed forests indicate that diameter increments most commonly increase with DBH (Anderson‐Teixeira, McGarvey, et al., 2015; Muller‐Landau et al., 2006). We hypothesize that this discrepancy is primarily a distinction between trees that establish in the open versus in the understory (Table 1, Lorimer et al., 1988). Building upon observed ontogenetic patterns in RW, dendrochronology studies often adopt a null hypothesis that the annual basal area increment (BAI) fluctuates around a constant mean after a juvenile growth phase (Biondi & Qeadan, 2008; Fritts, 1976)—a pattern that we would not expect to hold in understory‐established trees (Table 1). Finally, there is debate as to whether the aboveground biomass increment (ΔAGB) increases continuously with DBH (Foster et al., 2016; Meakem et al., 2018; Sillett et al., 2010; Stephenson et al., 2014) or peaks at intermediate DBH and then plateaus or declines as trees divert carbon to other functions such as reproduction and respiration (Sheil et al., 2017; Thomas, 2011; West, 2020). Following the finding that the latter pattern is common for individual trees whereas the former emerges in “cross‐sectional” analyses of forest stands (Forrester, 2021), we hypothesize that ΔAGB—and also BAI—peaks and declines as DBH increases (Table 1). Discerning these ontogenetic growth trends is essential not only for predicting the growth rate of any given tree but also for standardizing for DBH to deduce the influence of slowly changing environmental drivers (see next paragraph, Peters et al., 2015), with the reliability of such analyses contingent upon accurate assumptions of ontogenetic growth patterns.

Beyond the direct effects of climate, other factors, such as rising atmospheric CO2 concentrations, changes in atmospheric deposition of sulfur dioxide (SO2) and nitrogen oxides (NOx), and the indirect effects of climate change all potentially influence tree growth (Belmecheri et al., 2021; Levesque et al., 2017; Mathias & Thomas, 2018; Maxwell et al., 2019; Takahashi et al., 2020; Walker et al., 2021). Understanding these effects is central to predicting the future of the terrestrial carbon sink (Walker et al., 2021). Yet, characterizing how tree growth and forest productivity are responding to slowly changing environmental drivers is challenging and uncertain. Ontogenetic patterns in tree growth must be accounted for, yet two of the most commonly used methods of standardizing for tree size, conservative detrending and basal area correction (Peters et al., 2015), assume certain growth patterns unlikely to be universal in forest settings, as discussed above. Approaches that combine cross‐sectional with temporal analyses to correct for growth ontogeny, such as regional curve standardization, perform better at growth trend detection (Peters et al., 2015). Yet, even after correcting for ontogeny, growth trend detection remains subject to various potential sampling and analysis biases, which result in a limited potential of a contemporary set of tree cores to represent the growth history of a population (Bowman et al., 2013; Brienen et al., 2012, 2017; Cherubini et al., 1998; Duchesne et al., 2019; Hember et al., 2019; Nehrbass‐Ahles et al., 2014; Sullivan et al., 2016). Tree growth rates are sensitive to stand dynamics, with competition—the intensity of which tends to increase as forests mature—reducing woody growth rates (e.g., Zhang et al., 2015). Similarly, ecosystem‐level carbon allocation to woody growth—as opposed to leaf or fine root production, reproduction, defenses, etc.—has been shown to decline as forest stands age (Collalti et al., 2020; DeLucia et al., 2007; Goulden et al., 2011; Pregitzer & Euskirchen, 2004; West, 2020). Thus, we hypothesize that size‐corrected growth rates of tree populations sampled from within secondary or severely disturbed stands (i.e., those with large recruitment pulses within the past century) will generally decline, whereas populations sampled from older, relatively undisturbed stands will display mixed growth trends that are more dependent on external environmental drivers (Table 1).

We address the above hypotheses (Table 1) across 10 forested sites spanning 52° latitude, using a new method that allows simultaneous consideration of the effects of primary climate drivers (i.e., the most influential climate variables and the seasonal window over which they operate), DBH, and calendar year on annual tree growth.

MATERIALS AND METHODS

Data sources and preparation

We analyzed tree‐ring data, most of which were collected for earlier studies (see references in Table 2), from 10 sites ranging from 9.15° to 61.30°N latitude and representing a wide range of forest and tree types: tropical broadleaf deciduous and evergreen, temperate broadleaf deciduous and conifer, and boreal conifer (Table 2, Tables S1 and S2). Nine of these sites (exception: LT) were co‐located with large forest dynamics plots of the Forest Global Earth Observatory (ForestGEO, Anderson‐Teixeira, Davies, et al., 2015; Davies et al., 2021). Trees were cored within the ForestGEO plots (n = 5 sites) and/or nearby within similar forest types (n = 5 sites), following a variety of sampling protocols designed to meet the objectives of the original studies (Tables S1, S3). In using this diversity of data sources, we ensured that our approach could handle challenges presented by varying methodologies and forest types.

TABLE 2

Sites included in this analysis. Here and throughout, sites are ordered by descending mean annual temperature. Additional site information is provided in Appendix S1 and Table S1, and tree species and sampling details are detailed in Tables S2 and S3

Site code	Site name	Location	July T (°C) a	Jan T (°C) a	MAP (mm) a	Vegetation type(s)	n species	n cores	Original publication(s)
BCNM	Barro Colorado Nature Monument	Panama	26.6	25.5	2,627	BD, BE	3	84	Alfaro‐Sánchez et al., (2017)
HKK	Huai Kha Khaeng	Thailand	25.7	22.4	1,428	BD, BE	4	470	Vlam et al., (2014)
SCBI	Smithsonian Conservation Biology Institute	Virginia, USA	24.3	0.9	1,018	BD, C	14	704	Bourg et al., (2013); Helcoski et al., (2019)
LDW	Lilly Dickey Woods	Indiana, USA	24.0	−2.2	1,099	BD	6	170	Maxwell et al., (2016)
HF	Harvard Forest	Massachusetts, USA	21.6	−5.1	1,104	BD, C	4	366	Dye et al., (2016); Alexander et al., (2019); Finzi et al., (2020)
ZOF	Žofín Forest Dynamics Plot b	Czech Republic	18.1	−2.0	731	C, BD	4	2,059	Šamonil et al., (2013)
NIO	Niobrara	Nebraska, USA	23.4	−6.5	520	BD	1	138	Bumann et al., (2019)
LT	Little Tesuque b	New Mexico, USA	16.2	−3.1	608	C	2	34
CB	Cedar Breaks b	Utah, USA	13.8	−6.2	842	C, BD	7	187	Birch et al., (2020a), Birch et al., (2020b), Birch et al., (2020c), Birch et al., (2020d)
SC	Scotty Creek	Northwest Territories, Canada	16.5	−24.7	373	C	1	443	Sniderhan and Baltzer, (2016)

Abbreviations: BD, Broadleaf Deciduous; BE, Broadleaf Evergreen; C, Conifer; MAP, mean annual precipitation; T, mean monthly temperature.

Refers to 1950–2019 mean climate.

Older forest, with majority of sampled trees established before 1900.

All tree cores were cross‐dated and measured by the original researchers using standard dendrochronological practices (Stokes & Smiley, 1968). From among the full set of original RW measurements, we excluded cores for which we detected technical errors (e.g., labeling inconsistencies, obvious dating errors) that could not be resolved before finalizing the analysis. We also excluded records with small sample sizes or highly anomalous growth patterns, including (a) species with <7 cores, (b) cores with <30 years of record, (c) contiguous portions of cores containing large outliers (RW >mean plus 5 x SD of RW for the entire core), and (d) the final 20 years prior to death for trees cored when dead. The final criterion was implemented to avoid periods of growth decline and potentially altered climate sensitivity prior to death (Cailleret et al., 2017; DeSoto et al., 2020). From analyses including DBH (see below), we further excluded (a) trees for which we lacked data required to reconstruct DBH, (b) trees for which there was a significant inconsistency between measured DBH and the sum of RW’s across the core (Appendix S2), and (c) poorly represented tails of the DBH distribution, starting where reconstructed DBH (see below) included <3 conspecific trees. In total, this resulted in inclusion of 4655 cores from 3811 trees, 4513 of which (from 3705 trees) could be included in analyses with DBH (Table S3).

For each year in the tree‐ring records, we reconstructed (i.e., back‐calculated) DBH, as detailed in Appendix S2. We applied allometric equations for bark thickness to account for changes in bark thickness as the tree grew (Appendix S2; Tables S2, S4). Once DBH had been reconstructed, we estimated basal area (BA =  r 2, where r is radius) and aboveground biomass (AGB). Biomass allometries for temperate and tropical species were calculated using the R packages allodb (Gonzalez‐Akre et al., 2021) and BIOMASS (Réjou‐Méchain et al., 2017), respectively. We then calculated basal area increment (BAI = Bay+1‐BAy, where y is year) and aboveground biomass growth increments (ΔAGB = AGBy+1‐AGBy).

Monthly climate data for 1901–2019 were obtained from CRU v.4.04 (Harris et al., 2014, 2020), and in a few cases corrected based on higher‐resolution or local records (Appendix S3). Variables considered here included average daily minimum, maximum, and mean temperatures (T min, T max, T mean, respectively); precipitation (PPT); and, when deemed reliable (Appendix S3), potential evapotranspiration (PET) and precipitation day frequency (PDF). For the one riparian site, NIO, we tested for a relationship with stream flow, for which we obtained data for the Sparks, Nebraska station (station code: 06461500; 42°54’14“N, 100°26’13”W) from the U.S. Geological Survey (USGS) National Water Information System (https://waterdata.usgs.gov/nwis/uv/?site_no=06461500&agency_cd=USGS&referred_module=sw). All ForestGEO climate records used here are archived in the ForestGEO Climate Data Portal, v1.0 (https://doi.org/10.5281/ZENODO.3958215).

Data analysis

Data analysis consisted of two main steps: (a) identifying the primary climate drivers (i.e., variables and seasonal windows over which they are most influential on tree growth) and (b) combining these climate drivers, DBH, and year into a multivariate model (Figure 1). The analysis was run separately for each site (step 1), site–species combination (step 2), and each response variable estimating different measures of tree growth (RW, BAI, or Δ AGB). We note that the decision to identify primary climate drivers at the level of site, as opposed to species, was motivated by the expectation that differences in the most influential climate drivers across species in one site would be small compared with cross‐site differences (Figure 2); however, analyses focused on interspecific differences could optimize species‐specific climate sensitivity estimates by fitting individually by species.

FIGURE 1

Schematic illustration of the analysis process. In step 1, the R package climwin (van de Pol et al., 2016) is used to identify the primary climate drivers in water and temperature variable groups for each site, defined as the variable‐seasonal window combination that are most strongly correlated to the residual variation around splines fit to trends in growth (here, ring width, RW) for all cores sampled at the site. In step 2, a generalized last squares model is used to produce a combined model with the previously identified drivers, reconstructed diameter at breast height (DBH), and year. Example plots show raw data and partial effects of each variable from the best multivariate model for Pinus ponderosa P. Lawson & C. Lawson at Little Tesuque, New Mexico, USA

FIGURE 2

Example comparison of climate sensitivity derived via traditional methods (a) and our approach (b‐f). Example is for the sensitivity of 14 species at SCBI to potential evapotranspiration (PET). Panel (a) shows a matrix of Pearson correlations between ring‐width index and monthly climate variables (produced using the bootRes package in R, Zang & Biondi, 2013). Black rectangle represents the period selected by climwin as the most influential window. Panels (b‐d) give statistics for seasonal windows tested in climwin, where window open and close refer to the start‐ and end‐months of the window tested, expressed as months prior to current August, and cells across the lower diagonal indicate single‐month tests (akin to panel a). Panels (b) and (c) give values of linear and quadratic terms for each seasonal window, and (d) gives the difference in Akaike information criterion for small sample sizes ΔAICc for each. The seasonal window with the minimum ΔAICc (1–3 months prior to August, or May‐July; black circles), was identified as the most influential window. Panel (e) shows the correlation of individual‐level residuals to PET, with the function fit in climwin. Finally, panel (f) shows the generalized least squares model output, where PET was a candidate driver variable (along with PPT; DBH not included in this model). Plotted are responses of species for which PET was included in the top model, with best‐fit polynomials plotted with solid lines when both first‐ and second‐order terms are significant, dash‐dotted lines when only one term is significant, and dotted lines when neither is significant. Transparent ribbons indicate 95% confidence intervals. Species names corresponding to the codes are given in Table S2

Step 1: Identifying primary climate drivers

We used the climwin package in R (van de Pol et al., 2016) to identify the most influential climate variable (i.e., that most strongly correlated with annual growth) and the seasonal window over which its effect was strongest for each of two categories of variables: a temperature group (T min, T max, T mean, and PET) and a precipitation group (PPT, PDF). To remove low‐frequency variation that most likely represents responses to non‐climatic drivers (e.g., growth and aging of the tree, change in competitive dynamics, atmospheric pollution), we detrended the response variables by fitting penalized thin‐plate regression splines in generalized additive models (GAM, functions gam and s in the R package mgcv, Wood, 2011) to individual growth records (RW, BAI, or ΔAGB) from each core, and extracted the residual variation for each observation. The smoothing parameters were automatically selected by the gam function with generalized cross‐validation. Our application of the thin‐plate regression splines acts similar to more traditional a priori detrending methods using a two‐third spline commonly used in dendrochronology studies and results in similar predictor variable selection (Appendix S4; Cook & Peters, 1997; Rollinson et al., 2021).

We then used climwin to identify the climate drivers that most strongly correlated with the individual tree‐level residuals of the growth variables, RW, BAI, or ΔAGB, specifying linear and quadratic terms to allow for potential nonlinearities in the climate response. Within climwin, we specified a mixed‐effects model in which the fixed effects were the climate variables and the random intercepts were species (when n ≥3) and core identity (noting that these effects should be minimal given that residuals are centered around zero). For each climate variable, we ran permutations for all possible combinations of consecutive months within a 15‐month period ending near the time of growth cessation of each annual ring (Table S1). Climwin runs all potential models to select the best fit (lowest Akaike information criterion corrected for small sample size, AICc), and does k‐fold cross‐validation in its computation of AICc to guard against over‐fitting (van de Pol et al., 2016). For each group of candidate climate variables (water and temperature; Figure 1), we selected the variable ‐ seasonal window combination with the lowest AICc as a candidate climate variable for the multivariate models.

We tested whether this process identified similar seasonal windows and direction of response as would be identified using traditional methods for four species (detailed in Appendix S4). Furthermore, we explored alternate methods of climate driver selection for the two sites that have undergone the most rapid changes in climate and tree growth: LT, where increasingly warm drought has dramatically reduced growth (Touchan et al., 2011; Williams et al., 2013), and SC, where rapidly rising temperatures are causing permafrost thaw, which limits access to soil moisture during summer months and drives growth declines (Sniderhan & Baltzer, 2016). We ultimately determined that the method described above captured these sources of variation (Appendix S5).

Step 2: Combining drivers in generalized least square model

We next combined the primary climate drivers in temperature and precipitation variable groups (included in all models) and DBH (included in models with DBH and DBH–climate interactions) as candidate variables in linear mixed‐effects models (function lme in the R package nlme, Pinheiro et al., 2021). In all models, we included core identity as a random intercept and year as a continuous time covariate for the within‐group correlation structure (function corCAR1) to account for temporal autocorrelation (similar to how detrending would). We will refer to this model as a generalized least squares (GLS) model (Figure 1). Within the GLS models, our response variables were raw, log‐transformed growth estimates (as opposed to residuals): log[RW], log[BAI], or log[ΔAGB]. Prior to running the models, we checked for collinearity among the candidate variables using the vifstep function (Naimi et al., 2014). Our analysis code was programmed to remove any variable with a variance inflation factor >3, but none required removal.

For each species independently, we ran GLS models including every possible combination of independent fixed‐effect variables (i.e., candidate climate drivers, DBH, and year), including both first‐ and second‐order terms for each. For climate response, we allowed concave‐down fits but ignored any concave‐up fits on the basis that exponential functions would be captured by a linear fit to the log‐transformed growth variables, while u‐shaped fits are not expected biologically.

As an example, a full model for log[RW] responses to PPT, T max, and DBH would look like this in R: where x is a complete data set for one species at one site (all records after excluding cores as described above, and with no missing values). The method is set to maximum likelihood (ML) during the fixed effect model selection phase, but to restricted maximum likelihood (REML) for parameter estimation with the best model.

For models including interactive effects of climate and DBH, we tested for interactions between first‐order linear terms for climate drivers and DBH.

To test for year effects, we limited the analysis to species with reasonable coverage of the DBH x year matrix. Specifically, we required that the species be represented by cores from ≥3 trees and that the core record spanned ≥40% the total DBH range for ≥2/3 of the total time range analyzed. To avoid severe big‐tree selection bias (Brienen et al., 2012), we also required that the minimum DBH sampled be ≤25 cm (exception: Abies alba Mill. at ZOF, where mature trees <50 cm DBH are extremely rare). Species that failed to meet these criteria (n = 8; Table S3) were excluded from the analysis of temporal trends, but were included in analyses of climate and DBH and their interactions. We then ran models as described above, including a first‐order linear effect of year. We note that the random effect of tree may reduce analytical biases arising from persistent growth differences among individuals that are not accounted for by DBH or year (Brienen et al., 2012, 2017). To verify that GLS model trends for year were not an artifact of inherent covariation between DBH and year within each core, we compared GLS results with an analysis of DBH–growth relationships by decade.

Within each of three categories of models run (climate only, climate + DBH, climate × DBH, climate + DBH + year), we selected as the top model that with the lowest AICc. We did not run models of precipitation × temperature, climate × DBH + year, or climate × year, which would be possible in the GLS model framework but would require additional statistical and conceptual validation beyond the scope of the current analysis.

RESULTS

Validation of the method

Our process identified similar primary climate drivers to those identified via established dendrochronological analysis methods for identifying climate signals (Figure 2; Figures S1‐S4; Table S5; Appendix S4). Although one‐to‐one correspondence of estimated slope coefficients describing the response of tree growth to interannual climatic variation was neither expected nor observed, estimates obtained using the two methods were correlated and rarely differed significantly from one another (Appendix S4; Figures. S1‐S4).

Trends with year, when assessed, were generally consistent with those observed in a separate analysis of DBH–growth relationships by year (Figure S5).

Full model results overview

When a precipitation variable and a temperature variable (each selected using climwin; e.g., Figure 2; Figures S6‐S9), DBH, and calendar year were all included as candidate variables in the GLS models, the most common outcome was that all four were included in the top model and had statistically significant effects (p ≤ 0.05), regardless of the growth metric used (Figure 3). In general, DBH and calendar year explained more variation in growth rates than did climate, but their relative importance varied across growth metrics and sites (e.g., Figures S10–S13). Climate responses were generally similar regardless of the other variables included, although some of the weaker climate responses were not consistently included in top models (e.g., T min responses at BCNM; Figure 4; Figure S10). In contrast, effects of DBH and year often interacted such that the shape of the DBH response curve or its inclusion in the top model were frequently modified by the inclusion of calendar year (Figure 3).

FIGURE 3

Summary of top models for ring width as a function of climate (water and temperature variables), diameter at breast height (DBH), and calendar year. Arrow shapes approximate responses detailed in Figures 4 and 6, and S16. Each symbol indicates one species, and species are ordered alphabetically within each site. Overlapping arrows for the same species indicate that the response shape changed when Year was included in the model. For species on which the effect of Year was tested, the presence of only an arrow representing models without Year indicates that the effect was not included in the top models with Year. Interactive effects of climate and DBH are not shown. BCNM through SC are codes for ten forested sites spanning 52° latitude (Table 2); sites are ordered by descending mean annual temperature

FIGURE 4

Species‐level responses of ring width to climwin‐selected variables in precipitation and temperature variable groups for 10 sites. Models presented here include only climate variables as fixed effects. Primary climate drivers are coded on the x‐axes as the climate variable abbreviation followed by the range of months (p, previous year; c, current year) over which it is most influential. PPT, precipitation; PDF, precipitation day frequency; PET, potential evapotranspiration. For each species (color‐coded as in Figure 6), relationships are plotted if included in the top model. For each relationship shown, other terms in the model are held constant at their medians. Best‐fit polynomials are plotted with solid lines when both first‐ and second‐order terms are significant (t‐test's p‐value <.05), dash‐dotted lines when only one term is significant, and dotted lines when neither is significant. Transparent ribbons indicate 95% confidence intervals. Vertical grey lines indicate the long‐term mean for the climate driver over the analysis period; shading indicates 1 SD

FIGURE 6

Growth sensitivity to diameter at breast height (DBH): (a) ring width (RW), (b) basal area increment (BAI), (c) Δ aboveground biomass (ΔAGB). Models presented here include climate variables and DBH as fixed effects. Relationships for tree species in 10 sites (Table 2) are plotted when included in the top model. Other terms in the model are held constant at their medians. Best‐fit polynomials are plotted with solid lines when both first‐ and second‐order terms are significant (t‐test's p‐value <.05), dash‐dotted lines when only one term is significant, and dotted lines when neither is significant. Transparent ribbons indicate 95% confidence intervals. Species authorities are given in Table S2

Climate sensitivity

Most influential climate drivers

At each site, the three metrics of growth (RW, BAI, and ΔAGB) exhibited similar patterns in the direction of response, and relative strength of correlation, to climate variables across the range of potential seasonal windows. However, the seasonal window exhibiting the strongest climatic effect on growth, and even the most influential climate variable, sometimes differed among the growth metrics. For eight of 20 site–variable group combinations (i.e., water and temperature, each at 10 sites), both the most influential climate variable and seasonal window were identical across growth metrics (e.g., PPT at SCBI; Figure S7). For nine site–variable group combinations, climwin identified the same climate variable and overlapping seasonal windows (e.g., PET at SCBI Figure S8), and in one case (at HKK) different variables (Tmax and Tmean) were selected with overlapping seasonal windows (Figure S6). For just two site–variable group combinations (both variable groups at HF, where climate had only weak effects and mixed responses among species in the final models), climwin identified completely different seasonal windows and, for precipitation, different variables (PPT and PDF; Figure S9). Henceforth, unless otherwise noted, we focus on the climate sensitivities identified using RW as the growth metric and for the full set of cores (i.e., including those for which DBH could not be reconstructed).

Precipitation amount (PPT) was selected over precipitation frequency (PDF) as the top variable in five of the eight sites for which both variables were available (but had no significant main effect at one site, NIO), and was the only option at two sites (LT and CB). The most influential seasonal windows were most commonly long (≥3 months at 7 of the 9 sites with significant main effects) and coincided at least partially with months of active growth in the current year (Figure 4; Table S1): year‐round in the tropics (BCNM and HKK) or late spring/ summer outside of the tropics (n = 5 of 7 sites with significant main effects). In the tropics, the long time‐windows over which precipitation was influential (12 mo at BCNM, 9 mo at HKK) also included the majority (BCNM) or all (HKK) of the dry season months (<100 mm rainfall/month). Outside of the tropics, the most influential windows at three sites included the current growing season and extended back to the previous fall (LT, CB) or summer (SCBI), whereas they were limited to the current spring and early summer at LDW. At three sites (HF, ZOF, and SC), precipitation of the previous growing season was the most influential variable.

Within the temperature group (Figure 1), the most commonly selected variables were Tmax and PET, which were identified by climwin as the top temperature‐related driver at six and three of the 10 sites, respectively, noting that PET was not available for two sites. Tmin was identified as the top driver at BCNM, where its effects were only marginally significant for one species (Figure 4). Tmean was never selected as the top driver. The most influential seasonal windows for temperature tended to be shorter than those of precipitation (≤3 months at 9 of 10 sites). They most commonly occurred during the current growing season (n = 5 of 10 sites), but there were cases where the most influential windows occurred during the preceding dry season (BCNM), late winter/early spring (HF, ZOF), or the previous growing season (NIO, CB). Temperature and precipitation variables were rarely influential over the same seasonal window (exception: LDW).

Climate responses

Analyses of species‐specific responses at each site used the GLS model to test for first‐ and negative second‐ order linear effects of both a precipitation and a temperature variable. Both a precipitation and a temperature variable were included in the top model for 78% (n = 36 of 46) of site–species combinations (Figure 4). There were seven site–species combinations for which only a precipitation term was significant (two at BCNM, three at SCBI, and two at LDW), two for which only a temperature term was significant (Chukrasia tabularis A. Juss. at HKK and Betula papyrifera Marshall at NIO), and none with no significant climatic effects on RW. Below, we summarize the precipitation and temperature variables included in these models and their direction of response.

Responses to precipitation amount (PPT) and frequency (PDF) were included in the best model for all but two species–site combinations, and were predominantly positive (Table 1, Figure 4). Specifically, there were positive first‐order linear terms for precipitation for all but one species–site combination (Tsuga canadensis (L.) Carrière at HF; Figure 4). Negative second‐order terms were commonly included in the best model (32 of 42 with positive first‐order terms), generally resulting in a deceleration or decline at the highest levels of precipitation, but occasionally producing a unimodal (e.g., several species at SCBI) or predominantly negative response (e.g., Betula alleghaniensis Britton at HF; Figure 4).

A temperature variable was included in the best model for all but eight site–species combinations, with predominantly negative responses, particularly at the higher end of the temperature range (81%; 34% with negative first‐order term, 47% with positive first‐order term but negative second‐order term; Figure 4). Within the tropics, there was minimal effect of temperature at BCNM and a negative effect of wet season Tmax for three of four species at HKK. For temperate sites with the most influential seasonal windows covering the current and/or past growing season, responses were universally negative (i.e., negative first‐order linear or unimodal, peaking at temperatures lower than the long‐term mean). In contrast, there were positive effects of Jan‐March Tmax for all three species at ZOF and of March PET for T. canadensis at HF, the latter contrasting with a negative response of the three deciduous species analyzed at HF (Figure 4). At the highest‐latitude site (SC), which has undergone rapid warming and permafrost melt, Picea mariana (Mill.) Britton, Sterns & Poggenb. responded positively (but with wide 95% CI on the slope) to temperature over the full analysis period (1903–2013); however, responses were predominantly positive prior to 1970 and predominantly negative afterwards (Figure S14).

Variation in climate sensitivity with diameter at breast height

Interactive effects of climate and DBH were found for 91 of the 203 (45%) species‐site‐variable (RW, BAI, or ΔDBH) combinations for which they were tested. For precipitation variables, interactions were significant for 16 of the 36 (44%) interactions with RW as the growth metric (Figure S15) and for 17 of the 36 (47%) with BAI as the growth metric. The majority of these interactions were positive (75% for RW; 65% for BAI), indicating that larger trees generally respond more positively to precipitation or its frequency (Figure 5). Among the exceptions to this pattern (n = 4 for RW, 6 for BAI), only a minority (n = 1 for RW, 4 for BAI) occurred in species responding positively to precipitation in the current growing season.

FIGURE 5

Examples of climate–diameter at breast height (DBH) interactions for three tree species at three sites. Shown are modeled response functions at the minimum and maximum and maximum tails of the DBH distribution. Other terms in the model are held constant at their medians. Transparent ribbons indicate 95% confidence intervals. Vertical gray lines indicate the long‐term mean for the climate driver over the analysis period; shading indicates 1 SD. PPT, precipitation; PDF, precipitation day frequency; PET, potential evapotranspiration. Species authorities are given in Table S2

Temperature variable ×DBH interactions were significant for 38% of cases with RW as the growth metric (Figure S15) and for 50% with BAI as the growth metric. Directions of these interactions were mixed, with 5 of 12 significant interactions negative with RW as the growth metric and 10 of 16 significant interactions negative when BAI was the growth metric. For both RW and BAI, the majority of significant negative interactions (i.e., more negative/ less positive response of larger trees to higher temperatures) occurred in cases where the main effect temperature was negative (e.g., HKK, LT, CB; Figure 5), whereas positive interactions were more common when the main effect of temperature was positive (e.g., HF, ZOF).

Variation with diameter at breast height

Growth rate—whether measured as RW, BAI, or ΔAGB—varied with DBH for most species at all sites (Figure 6). Because the effects of calendar year could not be evaluated for all species (Figure 3), the DBH trends described here refer to models without year. Relationships between population‐mean growth rate and DBH were best described by models with second‐order terms for the majority of site–species combinations (81–98% depending on growth metric; Figure 6).

For RW, DBH was included in the best model for 81% of species–site combinations (n = 35 of 43), and the majority of best models also included a significant second‐order term (n = 26, 21 of which were negative). There was substantial variation in these trends, with patterns mixed across both forests and species within a single stand (Figure 6). On one end of the spectrum, some species exhibited maximum RW at low DBH, followed by fairly rapid declines in RW with increasing DBH. This pattern was common among light‐demanding species (6 of 13 site–species combinations; Table 1, S2; e.g., Melia azedarach L. at HKK, Populous tremuloides Michx. at CB) and in relatively open stands (e.g., both species at LT, P. mariana at SC; Figure 6). At the other end of the spectrum, some species had low RW at small DBH, increased to peak RW at intermediate DBH, and subsequently declined. This pattern was common among shade‐tolerant species (8 of 11 site–species combinations; Table 1; e.g., Trichilia tuberculata and Tetragastris panamensis at BCNM; Fagus spp. at SCBI and ZOF, Picea spp. at ZOF and CB; Table S2), and universal for shade‐tolerant species in models accounting for calendar year (Figure 3).

Trends in both BAI and ΔAGB were far more consistent across sites and species, typically increasing to a peak at intermediate DBH and then declining (Table 1, Figure 6). Best models for BAI included DBH and DBH2 for 42 of 43 species (exception: Acer rubrum L. at HF), with a positive coefficient for DBH in 40 (exceptions: non‐significant negative coefficients for Pinus ponderosa P. Lawson & C. Lawson at LT and Pinus longaeva D.K. Bailey at CB, whose reconstructed DBHs did not extend down to 0 cm within the time frame analyzed) and near‐universally negative coefficients for DBH2 (exception: P. longaeva at CB). For Δ AGB, models were even more consistent, with the best models for 98% of species containing a positive coefficient for DBH and a negative coefficient for DBH2 (exception: P. longaeva at CB).

Effects of year

There was a significant effect of year in the GLS models for 31–32 (depending on growth metric) of the 37 species–site combinations tested (Figure 7; Figure S16). In 90–91% of cases (depending on growth metric), the growth trend over time was negative. Declines were particularly prevalent in secondary or disturbed forests, occurring in 92% of species–site combinations (100% of all species with significant year effects) at the seven disturbed sites (Table 1, Figure 7). In older forests (ZOF, LT, CB), growth trends were mixed (Table 1, Figure 7, Figure S16). Significant positive growth trends were observed for only three species (consistently across all three growth metrics), Fagus sylvatica L. at ZOF, Picea pungens Engelm. and Pinus flexilis E. James at CB, and all were modest compared with the steep negative trends observed for some species. Growth rate was consistently independent of year for only four species: C. tabularis A. Juss. at HKK, Pinus strobiformis Engelm. at LT, and Picea engelmannii Engelm. and P. longaeva at CB.

FIGURE 7

Effect of year, when included in the best model, on basal area increment. Models presented here include climate variables, reconstructed diameter at breast height, and year as fixed effects for 10 sites (Table 2). For each tree species (all listed), relationships are plotted if the year effect could be analyzed and was included in the top model. Other terms in the model are held constant at their medians. Best‐fit polynomials are plotted with solid lines when both first‐ and second‐order terms are significant (t‐test's p‐value <.05), dash‐dotted lines when only one term is significant, and dotted lines when neither is significant. Transparent ribbons indicate 95% confidence intervals. Species authorities are given in Table S2

Effects of year and DBH interacted such that inclusion of year in models altered the shape of DBH responses, typically resulting in less pronounced growth declines with increasing DBH (Figure 3; Figure S11 and S12).

DISCUSSION

The long‐term growth records contained in tree rings provide an exceptional tool for understanding past drivers of growth and anticipating future forest changes, yet traditional dendrochronological analysis methods were not designed to disentangle multiple, simultaneously acting drivers of tree growth, nor their implications for whole‐forest productivity. Our novel method provides a powerful approach to elucidate how tree growth is simultaneously shaped by climate, tree size, slowly changing environmental conditions, and their interactions. Analyzed with respect to each of these drivers individually, our method yields results that are consistent with what would be obtained using established methods. Beyond this, because our approach considers these factors simultaneously, it allows analyses of their joint and interactive effects. Applied across a wide range of forest types and species distributed globally across 10 sites, we have shown that tree species vary in the shapes of their functional responses with respect to size‐related sensitivity to different climate variables. Dissecting these species‐specific long‐term responses is essential to understanding the drivers of variability and directional changes in tree growth over the past century and to predicting changes in forest composition and function in the future.

Across diverse climates and forest types (Table 2), growth rates of 40 tree species usually responded positively to water availability (PPT or PDF)—at least up until the long‐term mean—and negatively to temperature (usually T max or PET), with the exception of several positive responses at times and in places where temperature was limiting (Table 1, Figures 3 and 4). These findings are generally consistent with current understanding of global‐scale patterns in climate sensitivity (Babst et al., 2019; Rozendaal & Zuidema, 2011): outside of the wet tropics (where there are few tree‐ring records), the majority of forests are moisture limited and respond negatively to temperature, with a shrinking area of temperature‐limited forests in cold, humid regions (with SC falling near the transition zone). Within warmer regions, warm winter or early spring temperatures in humid climates may advance the growing season (Keenan et al., 2014) and increase annual growth (Babst et al., 2019; Tumajer et al., 2017), as we show for all three species at ZOF and one species at HF (Figure 4). However, the predominantly negative temperature responses (Figure 4) imply that warmer temperatures are likely to reduce growth across the wide range of forest types and climates represented here. The primary mechanism underlying growth decreases at high temperatures is presumably increased evaporative demand (PET or vapor pressure deficit, VPD) and ensuant exacerbation of observed water limitations (Humphrey et al., 2021; López et al., 2021; Novick et al., 2016). This effect occurs in addition to the effects of precipitation (Figure 4), highlighting the fact that temperature and associated VPD increases limit growth even under conditions of high soil moisture (Novick et al., 2016), and occurs over shorter time‐frames (usually ≤3 mo) than the effects of precipitation (usually ≥3 mo.; Table 1, Figure 4). This suggests that relatively short periods of anomalously high temperatures and evaporative demand, themselves caused in large part by soil dryness (Humphrey et al., 2021), add to effects of prolonged periods of reduced precipitation to shape forest drought responses.

Our analysis differed fundamentally from most conventional approaches in testing for nonlinear responses of growth to climate, finding that the majority of climate responses were nonlinear (Table 1, Figure 4). This result, which is consistent with physiological expectations (e.g., Kumarathunge et al., 2019; Wilmking et al., 2020), indicates that the majority of tree‐ring records examined here cover climate variation beyond the range over which the response is approximately linear. The nonlinear form of most climate growth responses implies that as the climate changes such that high temperatures and strong precipitation anomalies become more common (IPCC, 2021), nonstationary climate responses, already common (Wilmking et al., 2020), could become more prevalent (Germain & Lutz, 2020).

We found that interactions between climate drivers and DBH were common (45% of total cases analyzed; Table 1, Figure 5; Figure S15). The most coherent pattern observed in this analysis was a tendency for larger trees to be more sensitive to precipitation and high temperatures (Figure 5), consistent with widespread observations that larger trees are more sensitive to drought (e.g., Bennett et al., 2015; Gillerot et al., 2020; Hacket‐Pain et al., 2016; McGregor et al., 2021; Pretzsch et al., 2018). An analytical structure such as ours that can account for this pattern and other DBH–climate interactions (e.g. Rollinson et al., 2021; Rossi et al., 2007) will be critical to using tree‐ring records to understand and forecast the effects of climate on tree growth and forest productivity.

There was substantial variation across species–site combinations in the population mean relationship between DBH and growth rate (Figure 6). Variation was most pronounced when RW was considered as the growth metric, as would be expected based on basic geometric principles given that RW patterns are most variable at small DBH. This variation was driven by two primary, interrelated factors: species ecology and stand history. As hypothesized, on average RW declined with DBH for nearly half of light‐demanding species, but most commonly RW initially increased across the lower end of the DBH range for shade‐tolerant species (Table 1), particularly when the effects of calendar year were accounted for in the model (Figure 3). However, species shade tolerance alone did not explain variation in RW–DBH relationships; rather, we observed instances where, on average, RW declined with DBH for a shade‐tolerant species growing in a relatively open stand (P. mariana at SC) or initially increased with DBH for shade‐intolerant species growing at sites where competition for light was likely more intense (e.g., Afzelia xylocarpa (Kurz) Craib at HKK, Liriodendron tulipifera L. at LDW). These results imply that while species that typically grow in high‐light conditions commonly display dendrochronology's “textbook” pattern of declining RW with DBH—in part attributable to the geometric constraint that new growth is spread around an ever‐growing circumference (Biondi & Qeadan, 2008; Fritts, 1976)—the majority of trees within forest settings exhibit hump‐shaped patterns of RW in relation to DBH. This latter pattern is consistent with the observation that when contemporary growth rates are compared across individuals within a closed‐canopy stand (e.g., a “cross‐sectional” analysis of census data), average RW increases continuously across most of the DBH range (e.g., Anderson‐Teixeira, McGarvey, et al., 2015; Muller‐Landau et al., 2006), or increases and subsequently decreases (Schelhaas et al., 2018). Our finding that on average across populations, BAI and ΔAGB generally saturate or decline with increasing DBH (Table 1, Figure 6) contrasts with findings of cross‐sectional analyses of forest census data showing that mean ΔAGB increases continuously with DBH (Meakem et al., 2018; Stephenson et al., 2014). In large part, this discrepancy can be explained by differences between cross‐sectional analyses and “longitudinal” patterns of individual trees through time (Forrester, 2021; Sheil et al., 2017), consistent with the principle that individual‐scale growth patterns and environmental responses do not necessarily match population‐ or stand‐level average responses (Clark et al., 2003, 2011). Declines in BAI and ΔAGB at larger DBH may be in part attributable to increasing carbon allocation to functions other than woody growth, such as reproduction (Thomas, 2011), and possibly to physiological declines as trees age (Forrester, 2021; Qiu et al., 2021). Apparent declines in ΔAGB at large DBH (or old age) may also be driven by shifts towards proportionally greater wood production within the crown (e.g., branch production, Sillett et al., 2010, 2021) that are not adequately captured by biomass allometries based on DBH and sometimes height (Disney et al., 2020; Goodman et al., 2014). Growth declines may also be linked to slowly changing environmental conditions (e.g., successional changes in stand structure, climate change). Notably, inclusion of year in the GLS models reduced the frequency of RW declines with DBH (Figure 3) and tended to reduce the magnitude of BAI and ΔAGB declines at larger DBH (Figures S11 and S12), suggesting that some of the growth declines at large DBH (Figure 6) are more properly attributed to recent environmental changes than to large DBH.

Changing growth rates

Our analytical framework reconstructs growth changes in a sampled tree population over time while accounting for climate, DBH, and persistent growth differences among individuals (Figure 1), thereby addressing some important challenges to obtaining unbiased estimates of growth trends attributable to non‐climatic environmental drivers. First, we account for changes in climate that may drive directional growth trends. For example, dramatic growth declines at LT, driven by a strong regional warming and drying trend (Touchan et al., 2011; Williams et al., 2013), are in part factored out by accounting for the primary climate drivers, such that a significant decline over time was detected for only one of the two species (Figure 7). Second, we show that growth rate—by any metric—varies nonlinearly with DBH and with patterns dependent upon the species and environmental context (Figure 6), reinforcing the concept that growth trend analyses should incorporate cross‐sectional analyses to correct for growth ontogeny (e.g., through regional curve standardization, Peters et al., 2015). Our method does this, differing from the conceptually parallel method of regional curve standardization in that we standardize relative to DBH rather than age, correct for any trends in the most influential climate drivers, and include random effects of tree to account for persistent growth differences among individuals. The latter addresses a third important challenge, as those growth differences among individuals can bias estimated growth trends in positive or negative directions (Brienen et al., 2012, 2017; Groenendijk et al., 2015; Nehrbass‐Ahles et al., 2014; van der Sleen et al., 2017). For instance, older trees, which provide the only records available for the earliest decades, may be competitive winners that had above‐average growth rates within their cohorts (Aubry‐Kientz et al., 2015), which would upwardly bias average growth rate estimates for early decades (“juvenile selection effect,” Groenendijk et al., 2015). In contrast, the oldest age classes being dominated by trees with below‐average growth rates (e.g., Piovesan et al., 2019) could downwardly bias average growth rate estimates for early decades (“a slow‐grower survivorship bias,” Brienen et al., 2012; Duchesne et al., 2019). By including a random effect of tree, our approach likely reduces the most severe potential biases associated with persistent growth differences across individuals (Bowman et al., 2013; Brienen et al., 2012, 2017), yet observed trends nevertheless represent only the sampled population of trees, as opposed to tree populations at all points throughout the time frame analyzed. Within this context, signals of changing growth rate over time are attributable to some combination of stand dynamics (e.g., recruitment and succession, changing stand structure) and environmental drivers (e.g., indirect effects of climate change, rising atmospheric CO2, deposition of SO2 and NOx).

In all seven sites dominated by trees less than 100 years old, population‐mean growth trends were universally negative, when significant (Table 1, Figures 3 and 7). Although attribution of these trends is difficult, it is likely that some trends are caused by stand dynamics as cohorts and stands develop over time. Such negative trends are fairly typical of mixed‐species stands that experience vertical stratification (Oliver & Larson, 1996). For species exhibiting a pulse of recruitment in the past followed by little subsequent recruitment (e.g., A. rubrum and B. alleghaniensis at HF), persistent differences in growth rates among individuals could produce a trend of declining growth, as faster‐growing individuals reach various size thresholds earlier (Brienen et al., 2017; see also van der Sleen et al., 2017). Particularly in secondary stands where many of the sampled species recruited in pulses that were followed by low recruitment (e.g., SCBI, HF; Appendix S1), growth declines are consistent with the tendency for faster tree growth during early stand development (Lorimer & Frelich, 1989; Lorimer et al., 1988; Oliver & Larson, 1996), and with increasing competition and declining woody productivity as young stands mature (e.g., Goulden et al., 2011; Pregitzer & Euskirchen, 2004; West, 2020). Gradual shifts in abiotic drivers (e.g., indirect effects of warming) likely also play a role at some of these sites. At Scotty Creek, in northern Canada, rapid warming is thawing permafrost and altering hydrologic conditions (Baltzer et al., 2014), resulting in high mortality, growth declines, and low recruitment of P. mariana (Dearborn et al., 2020; Sniderhan & Baltzer, 2016). At this site, we attribute pronounced negative growth trends to a combination of successional declines and indirect climatic stress.

Within the three older forests (ZOF, LT, CB), population‐mean growth trends were mixed (Table 1, Figures 3 and 7), probably reflecting some combination of successional changes, changing mortality rates, and shifting competitive advantages, perhaps in part driven by changing environmental conditions (Furniss et al., 2017; Vrška et al., 2009) or the lack of intermediate disturbances giving rise to increasing crowding (e.g., Lutz et al., 2009). In particular, light‐demanding species that establish in gaps (e.g., P. tremuloides Michx. at CB; Table S2) would tend to experience an increasingly competitive environment through time. At Zofin, size‐corrected growth rates were lowest in the 1970s and 1980s, consistent with other studies from central Europe showing dramatic growth reductions due to acid deposition during this period (Elling et al., 2009; Šamonil & Vrška, 2008). Nonlinear trends such as this would be more accurately described by a nonlinear response function to year, or incorporation of data on pollution, but that is beyond the scope of the current analysis. Notably, there were only three species—F. sylvatica at ZOF and P. pungens and P. flexilis at CB—whose population‐mean growth rate increased significantly over the analysis time frame (Figure 7).

The rarity of positive growth trends observed here indicates that any growth benefit from elevated CO2 was outweighed by some combination of demographic changes and chronic environmental shifts. This aligns with the preponderance of studies using tree rings to infer growth responses to rising CO2 (e.g., Girardin et al., 2016; Groenendijk et al., 2015; Hararuk et al., 2019; Walker et al., 2021), including previous analyses from HKK (Groenendijk et al., 2015; Nock et al., 2011; van der Sleen et al., 2015, 2017), although some studies have detected growth increases (e.g., Hember et al., 2019; Voelker et al., 2006). A growth benefit of increasing atmospheric CO2 concentration is expected, based on physiological mechanisms, under water‐limited conditions and has been observed in young forests in experimental settings (Walker et al., 2021). However, significant woody growth stimulation by elevated CO2 has not been observed in experimentally manipulated mature forests (Walker et al., 2021), and increasing CO2 does not appear to be a dominant growth driver for the trees in natural forest settings analyzed here.

CONCLUSIONS

As global change pressures intensify and the need to understand changing forest dynamics becomes increasingly urgent (McDowell et al., 2020; Thom et al., 2017), we expect that the approach presented here will prove valuable to understanding drivers of tree growth and forest change. Multiple elements of global change—including changing atmospheric composition, warming, drought, changing disturbance regimes, and thawing permafrost—are simultaneously influencing forests worldwide (e.g., Anderson‐Teixeira, Davies, et al., 2015) and are expected to interact with each other, differentially affecting trees of different size and species (Table 1, Figures 3, 4, 5). Simultaneous evaluation of multiple drivers can improve analysis, understanding, and predictions of global change effects on tree growth and forest productivity. The method we present and apply here allows doing so.

AUTHORS’ CONTRIBUTIONS

Kristina J. Anderson‐Teixeira, Valentine Herrmann, Christine R. Rollinson, M. Ross Alexander, and Camille Piponiot conceived the ideas and designed methodology; Neil Pederson, Craig D. Allen, Raquel Alfaro‐Sánchez, Tala Awada, Jennifer L. Baltzer, Joseph D. Birch, Sarayudh Bunyavejchewin, Paolo Cherubini, Ryan Helcoski, Jakub Kašpar, James A. Lutz, Ellis Q. Margolis, Justin T. Maxwell, Pavel Šamonil, Anastasia E. Sniderhan, Alan J. Tepley, Ivana Vašíčková, Mart Vlam, and Pieter A. Zuidema collected the data; Valentine Herrmann, Bianca Gonzalez, Erika B. Gonzalez‐Akre, Cameron Dow, and Neil Pederson organized and analyzed the data; Kristina J. Anderson‐Teixeira led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.

Supplementary Material

Click here for additional data file.