Study area and forest inventory data
To examine the temporal compositional shifts, we used a large network of permanent sampling plots (PSP), which were established by the provincial governments of British Columbia, Alberta, Saskatchewan, Manitoba, Ontario, Quebec, Nova Scotia, Newfoundland, and Labrador between the 1950’s and 1980’s (Fig. 1). We selected the PSPs that fit the following criteria. The plots must: (i) be unmanaged, with a known stand age (year); (ii) have all trees tagged and repeatedly measured; (iii) have all trees marked with their diameter at breast height (DBH). A total of 17,107 plots (914.21 ha; 43°47’N–60°00’ N, 52°81’W–133°71’ W) were selected for our analyses, with 1,438,577 trees measured during the monitoring period of from 1951 to 2016. The average measurement interval was 9.47 years with 4.70 census times, where the initial and final census years varied as 1951-2011 and 1956-2016, respectively. The plot sizes ranged from between 20 m2 and 2,023 m2 (Table S1). The mean annual temperature and precipitation in the area varied between -3.91 °C and 12.26 °C, and between 291 mm and 3,884 mm (1951-2016), respectively. The elevation ranged from 0.1 m to 2,355 m above sea level (Table S1).
Functional composition
To quantify functional composition, we employed eight key functional traits related to growth and competitive abilities, as well as environmental tolerance capacities, based on previous studies 18,20,21: ‘leaf nitrogen content per leaf dry mass’ (Nmass, mg g-1), ‘leaf phosphorus content per leaf dry mass’ (Pmass, mg g-1), ‘specific leaf area’ (SLA, mm2 mg-1; i.e., leaf area per leaf dry mass), ‘wood density’ (WD, g cm-3), ‘shade tolerance’ (ST, categorical class 1–5 18), ‘drought tolerance’ (DT, categorical class 1–5 18, ‘leaf habit’ (‘deciduous’ = ‘1’ vs ‘evergreen’ = ‘0’), and ‘leaf structure’ (‘broadleaves’ = ‘1’ vs ‘coniferous’ = ‘0’). These trait values were extracted from the TRY database 22 and other published sources 18,23-26. We obtained trait data for > 92% of all species. For those minor species of which trait data was absent, we used genus-level trait values (averaged for the genus) 27,28.
The functional composition was represented by the CWM 4,17,29 that weighs trait values according to the relative abundance of each species based on DBH 15. Similar to previous studies 30, we performed a PCA with the CWMs of the eight traits to obtain a comprehensive functional identity to represent them, as these values were highly correlated with each other (Fig. 2). We employed the first (CWMPC1, explained 60% of the variation) and second axes (CWMPC2, explained 22% of the variation) as proxies for functional composition. The CWMPC1 collectively represented traits associated with deciduous broadleaved trees and higher resource acquisition 18,21,23,24, being positively related with CWMNmass, CWMPmass, CWMSLA, CWMHabit (i.e., deciduous), CWMStruct (i.e., broadleaves), and CWMWD. On the contrary, the CWMPC2 represented traits associated with environmental tolerance, being negatively associated with CWMDT and positively related with CWMST (Fig. 2).
Stand age
Stand age (SA) represents changes in stem density and composition associated with forest succession 15,19) The SA of each plot was determined through dendrochronological aging based on the average age of the oldest species in the stand. We employed SA to account for the effects of forest development processes on functional composition.
Global environmental change drivers
Similar to previous studies 4,15, we used the calendar year (Year), which represented the effects of global environmental change overall on functional composition. For global environmental change drivers, we derived CO2 measurements from the Mauna Loa Earth System Research Laboratory in Hawaii (http://www.esrl.noaa.gov/gmd/ccgg/trends/co2_data_mlo.html), and annual mean temperature, as well as annual mean precipitation and potential evapotranspiration, using BioSIM 11 software 31. BioSIM generates plot-level climates, based on the simulation using daily observations and monthly historical statistics from the sampled points (latitude/longitude), being adjusted by differences in elevation. Therefore, the generated climate data was unique to each plot. Subsequently, we calculated the annual climate moisture index (CMI; mean annual precipitation minus potential evapotranspiration 32). Following a previous study 15,19, we calculated the anomalies of annual mean temperature (ATA) and climate moisture index (ACMIA), which were defined as a deviation from their long-term means between 1951 and 2016 33. CMI is extensively employed as an indicator of drought conditions in Canada 15,19,32. Negative values indicate drier conditions, while positive values denote wetter conditions.
Baseline climate
Following previous studies 13,14, we calculated the long-term averages between 1951 and 2016 of annual CMI (CMIave) and MAT (MATave) for each plot, as proxies for site-specific baseline climates (i.e., local historical climate; Table S1; Fig. 1).
Statistical analysis
To examine the temporal trends of functional shifts associated with spatial variations in baseline climates, we employed the following linear mixed models:
(CWMPC1)ijkl or (CWMPC2)ijkl = β0 + β1 × (PS)j + β2 × (Prov)j + β3 × f(SA)ij + β4 × (Year)i +
β5 × (CMIave)j + β6 × (MATave)j + β7 × f(SA)ij × (Year)i +
β8 × f(SA)ij × (CMIave)j + β9 × f(SA)ij × (MATave)j +
β10 ×(Year)ij × (CMIave)j + β11 × (Year)ij × (MATave)j +
β12 × (CMIave)j × (MATave)j +πj + ε (1)
where i and j are ith census and jth plot; CWMPC1 and CWMPC2 are community-weighted mean of ‘broadleaves vs conifers traits’ 30 and ‘stress-tolerance traits’ 18, respectively; β are the coefficients to be estimated; SA is the stand age being transformed by a square root function f based on Akaike Information Criterion (AIC); Year is the calendar year representing long-term global environmental change effects 15,19. To control for potential influences of plot size on composition 34, we included plot size (PS) in the model as a covariate. We also inculded province (Prov) to account for differences in sampling methods (e.g., DBH threshold for tree measurement 35) among provinces. The identities of each plot was included as a random effect (pj) to take locally unique conditions (site-specific disturbance histories; e.g., short-term climate events, insect outbreaks, non-catastrophic small fire/wind/flooding disturbances) and spatial autocorrelation structures into account. ε was a random error. All of the two-way interaction terms were included, as the model that included these showed a consistently lower AIC according to our preliminary analysis. The maximum variance inflation factor (VIF) was 2.77 for the CWMPC1 model and 2.79 for the CWMPC2 model, indicating that multicollinearity was not an issue.
As the measurement interval (years) varied between censuses, we performed variation partitioning involving those fixed variables above and distance-based Moran’s eigenvector maps (dbMEM) 36 to explicitly factor out the influences of temporal autocorrelation on functional composition. A dbMEM matrix was generated based on the calendar year as explanatory variables of temporaly correlated structures, using the adespatial package 37 in R. We initally selected 24 MEMs that well represented temporaly correlated structures, using Moran’s I statistic with 1,000 random permutation 37. We then selected nine dbMEMs for the CWMPC1 model and 11 dbMEMs for the CWMPC2 model by stepwise selection and added them to eqn. 1, as well as subsequent analyses with global environmental drivers (eqn. 2), as covariates to account for temporal autocorrelations 36. After modelling with the dbMEMs, there was no significant temporal autocorrelation (examined by autocorrelation function estimation using the acf function in the stats package). Note that the coefficient estimates of the dbMEMs are shown separately for brevity since it is not our intension to understand the importance of autocorrelative structure (Fig. S6).
To account for uncertainties in sampling, models, and parameters we employed Bayesian Markov chain Monte Carlo methods for linear mixed models, using the MCMCglmm package 38. To obtain a reliable posterior distribution with a satisfactorily effective sample size (i.e., the size of an uncorrelated sample), we used a thinning interval with a lag of 50 (examined by autocorr.diag function in the MCMCglmm package). Thus, we ran the models for 53,000 iterations with a burn-in period of 3,000 and thinning interval of 50 to achieve the recommended >1,000 effective sample size 39 (checked by effectiveSize function in the MCMCglmm package). We estimated the posterior distribution with a sampling of 1,000 in accordance with the default and confirmed that the performance stabilized with no autocorrelation 40. All explanatory variables were centred and scaled (mean = 0, SD = 1) prior to analysis to allow a coefficient comparison.
We also examined the temporal trends in atmospheric CO2 concentrations, ATA, and ACMIA, and how they were associated with the CMIave and MATave via linear fixed effects models. To investigate the associations between the CWMs and rates of global environmental change, we replaced baseline climate variables (CMIave and MATave) in eqn. 1 with temporal change rates in ACMIA (ACMIAChangeRate, cm/yr; Fig. 1c) and ATA (ATAChangeRate, °C/yr; Fig. 1d). Furthermore, we explored the roles of global environmental drivers on the CWMs using the following model (simultaneous modelling with three environmental change drivers rather than the Year term in eqn. 1):
(CWMPC1)ijkl or (CWMPC2)ijkl = β0 + β1 × (PS)j + β2 × (Prov)j + β3 × f(SA)ij + β4 × (CO2)i +
β5 × (ATA)i + β6 × (ACMIA)i + b7× (CMIave)j + β8 × (MATave)j +
β9 × f(SA)ij × (CO2)i + β10 × f(SA)ij × (ATA)ij +
β11 × f(SA)ij × (ACMIA)ij + β12 × f(SA)ij × (CMIave)j +
β13 × f(SA)ij × (MATave)j + β14 × (CO2)i × f(CMIave)j +
β15 × f(ATA)ij × (CMIave)j + β16 × f(ACMIA)ij × (CMIave)j +
β17 × (CO2)i × f(MATave)j + β18 × f(ATA)ij × (MATave)j +
β19 × f(ACMIA)ij × (MATave)j + β20 × f(CMIave)j × (MATave)j +
πj + ε (2)
where CO2, ATA, and ACMIA are the atmospheric CO2 concentration, anomalies of mean annual temperature, and climate moisture index, with which the CO2 and ATA in our data was positively correlated (r2 = 0.19). As the maximum VIF in this model was = 4.33 for the CWMPC1 model and 4.34 for the CWMPC2 model, concerning the multicollinearity, we also modelled CWMs with each of the climate drivers separately (one driver at a time). This preliminary attempt showed that the coefficient estimates did not qualitatively differed among these models and eqn. 2 (Fig. S4). Therefore, similar to a previous study 19, we focused on the outcomes from the simultaneous model with eqn.2. Conditional and marginal R2 for eqns. 1 and 2 are shown in Table S3.
To understand the functional response processes to the calendar year, and global environmental change drivers, we calculated genus-level relative abundance (%) by sub-setting the basal area (m2/ha) by major tree genus. Similar to a previous study 19, major tree genus was defined as those that accounted for >5% of the total basal area across all of the plots during the entire census, and occurred in all the provinces: Picea spp. (26.7%); Abies spp. (11.8%); Populus spp. (8.5%); Acer spp. (7.9%); and Pinus spp. (15.4%) (Table S2). The basal areas of individual stems were summed to obtain the overall basal area. The relative abundance of each major genus was calculated as the proportion of its basal area to the total basal area of the stand at each census for each plot, which was then multiplied by 100 to obtain an abundance percentage 15. Similar to a previous study 15, we then examined the responses of the relative abundance of each genus to the calendar year, with the same fixed-effects parameters used in eqn. 1, as well as the three global environmental change drivers with the same fixed effects used in eqn. 2.
For the interpretation of all analyses, we focused on not the statistical significance (i.e., p-value) but the ‘ecological significance’; that is, the effect sizes and the directionality and steepness of slopes (positive/negative/neutral directionalities). If these elements were substantially different, we interpreted as an ecologically meaningful trend in functional shifts, while we considered statistically significant but small effect sizes that result in qualitatively similar slopes as negligible difference. Although such evaluation scheme cannot offer an exact threshold for conclusion in comparison to statistical significance level, we advocate that statistical significance does not necessarily equal to ecological significance, and statistical assessment based on effect size should be the standard, according to the recent statement by the American Statistical Association 41.