Functional shifts in natural forests under environmental change over the last 65 years are faster in colder regions

Global environmental changes have signicantly impacted plant diversity and composition over many decades. Maintaining biodiversity and composition is critical for sustainability of ecosystem functioning and related services. While global environmental changes have modied plant species and functional compositions in forest ecosystems, it remains unclear how temporal shifts in functional composition differ across regions and biomes. Utilizing extensive spatial and long-term forest inventory data (17,107 plots monitored 1951–2016) across Canada, we found that functional composition shifted toward fast-growing deciduous broadleaved trees and higher drought tolerance over time; notably, this functional shift was more rapid in colder regions. Further analysis revealed that the functional composition of colder plots shifted toward drought tolerance more rapidly with rising CO2 than warmer plots, which suggests the vulnerability of the functional composition of colder plots against global environmental changes. Future ecosystem management practices should consider spatial differences in functional responses to global environmental change, with particular attention to colder plots that experience higher rates of warming and compositional changes.

to affect the temporal trends of biomass changes in boreal forests. For example, the growth of boreal forests in colder regions that experienced lower, or no changes in water availability, were less negatively affected by long-term climate change compared to boreal forests in the warmer regions of Eastern Canada 13 . Moreover, boreal forests in more humid regions suffered a lower extent of biomass loss under long-term changes in climate in Western Canada 14 . If long-term global environmental change favours tree species with traits that are better adapted to new climate realities (while causing higher mortality for species with unfavoured traits), these differences in demographic changes may induce spatially divergent shifts in functional composition. However, exactly how temporal shifts in functional traits are associated with spatial gradients of the baseline climate have rarely been tested, with insights being limited only to a water availability gradient in tropical forests 7 . Yet, we are not cognizant of how the rate and directionality of compositional shifts are dependent on larger scale environmental contexts across regions and biomes, particularly as relates to baseline temperature. Moreover, no study as yet has determined the relative contributions of these regionally dependent environmental change drivers to shifts in functional composition of natural forests at a macroecological, cross-biome scale.
For this paper, we explored how the directionality and rate of temporal functional changes of natural forests in response to persistent long-term global environmental change have been dependent on baseline climatic conditions over 65 years across multiple biomes in Canada. We hypothesized that, overall, the compositions of forests across Canada would shift toward resource acquisitive deciduous broadleaved trees in response to rising CO 2 and temperatures 3,15 . Further, the functional composition of forests in drier plots would respond more quickly to long-term global environmental change than humid plots, by increasing community-level drought tolerance due to their greater susceptibility to decreasing water availability 14 . Moreover, the functional shifts of colder plots would be more prominent than their warmer plots due to the higher propensity for increased temperatures at higher latitudes 16 (Fig. 1). To test these hypotheses, we surveyed the data for 17,107 permanent sampling plots (with 1,471,165 trees naturally regenerated following wild re and unmanaged) of temperate and boreal forests monitored between 1951 and 2016 across Canada (Fig. 1). We represented functional composition by the community-weighted mean (CWM) 17 of eight traits associated with competitive and tolerative abilities. These included leaf nitrogen and phosphorus content per leaf dry mass, speci c leaf area, wood density, shade tolerance, drought tolerance 18 , leaf habit, and leaf structure (Fig. 2). To compress the dimensionality, we used the rst and second axes of the principal component analysis (PCA). The rst axis was correlated with deciduous broadleaved trees vs. conifers (CWM PC1 , positively associated with leaf nitrogen and phosphorus content, speci c leaf area, and wood density), whereas the second axis was negatively correlated with drought tolerance (CWM PC2 × -1; converted to make it positively associated with drought tolerance; Fig. 2). Subsequently, we analyzed temporal trends in these functional composition metrics over 65 years (representing changes in atmospheric CO 2 concentration, temperature, and water availability), while simultaneously accounting for the in uences of stand development 19 and spatial variations in the baseline climate (i.e., long-term averages of mean annual temperature (MAT ave ) and climate moisture index (CMI ave ) 13,14 . Further, we examined the relationships between functional shifts and temporal trends in atmospheric CO 2 concentrations, anomalies of mean annual temperature (ATA), and the climate moisture index (ACMIA).

Results And Discussions
Across all study plots and temporally repeated measurements, stand age (0.2 to 379 years; succession following stand-replacing disturbances) accounted for more variations in both functional composition metrics, and it had greater effect sizes on them than background temporal changes (i.e., calendar year) ( Fig. 3). Overall, the functional composition shifted toward conifers and lower drought tolerance (or higher shade tolerance; see Fig. 2) with stand age (Fig. 4a). Increased baseline water availability (CMI ave ) was the most strongly associated with conifers and lower drought tolerance (Fig. 4b), while baseline temperature (MAT ave ) showed the strongest association with deciduous broadleaved trees and lower drought tolerance (Fig. 4c).
Even after factoring out the potent in uences of stand age and baseline climate covariates, we found that functional composition shifted toward deciduous broadleaved trees and higher drought tolerance (or lower shade tolerance) over the calendar year (see the main black lines in Fig. 5). This was consistent with a previous study of the boreal and temperate forests of Western Canada 15 , Eastern USA 9 , and Europe 29 . Furthermore, our analysis revealed that the functional shift toward drought tolerance over the calendar year was strongly modulated by baseline temperature (Fig. 5d), while the shift speed was nearly consistent regardless of water availability (Figs. 3, 5c). Speci cally, drought tolerance increased (or shade tolerance decreased) with the calendar year more rapidly in colder plots than warmer plots (Fig. 5d), and even the slope direction was opposite in the case of warmer plots. The temporal shift toward deciduous broadleaved trees was consistent across the baseline climate (Figs. 3, 5a). To understand the functional shift process, we examined temporal changes in the relative abundance of major tree genera (Fig. S1).
We found that the observed trends in functional composition, as related to the baseline temperature, were due to greater temporal increases in the relative abundance of deciduous broadleaved trees and earlysuccessional conifers (particularly Betula and Pinus) with a reduction in drought intolerant conifers (Picea; see Table S2 for genus-level trait values) in colder plots (Fig. S1b). This suggested the higher sensitivity of functional composition in boreal forests at higher latitudes ( Fig. 1b) under long-term global environmental change.
Over the 65 years, the atmospheric CO 2 concentrations (Fig. 1e) and temperature (ATA) increased across the study area (Figs. 1d, f, g). However, the temperature rose more quickly in wetter and colder sites, in contrast to drier and warmer sites (Figs. 1f, g). Across the study area, water availability (ACMIA) showed a convex curve, which increased and then decreased over the calendar year (Figs. 1f, g). However, drier sites experienced more substantial temporal variations in water availability, although the most humid sites exhibited a gradual concave curve, which decreased and then increased over the calendar year (Fig. 1f). The temporal trend in water availability was consistently concaved with the baseline temperature ( Fig.   S1g).
To investigate whether the expedited compositional shifts at colder sites were the result of more rapidly warming temperatures, we replaced the baseline climate variables in eqn. 1 with temporal change rates in temperature (ATA ChangeRate ) and water availability (ACMIA ChangeRate ) (Figs. 1c, d, S2). We found that the functional shift toward drought tolerance was more prominent in plots with higher change rates in temperature, ATA ChangeRate, than those in water availability (Fig. S3). Thus, the faster shifts in functional composition towards drought tolerance/shade intolerance in colder plots were likely due to more rapid warming than the change in water availability in those plots (Figs. 1b, g).
We then tested whether the association of functional composition with these global environmental change drivers would be dependent on the baseline climate. We employed two alternative approaches (one driver at a time, and all three drivers simultaneously, using a linear mixed-effect model, respectively) to model the main and interactive effects of individual drivers and baseline climate on functional composition. These approaches yielded similar coe cient estimates (Fig. S4). As interpreted from the results of the linear mixed models with all three drivers modelled simultaneously, CO 2 had the greatest effect size on both types of functional composition (Fig. S4).
The association of functional composition with rising CO 2 levels largely mirrored that of the calendar year, due to their high correlation (r 2 = 0.99) (Figs. 6a, d). Similar to previous studies 3, 15 , increased CO 2 concentraiton was associated with deciduous broadleaved trees and higher drought tolerance across the study area (see black average lines in Figs. 6a, d). However, our new nding was that although the response slope of broadleaves vs conifers to rising CO 2 was consistent across the spatial gradient of baseline climate (Fig. 6a), the positive relationship between rising CO 2 levels and the functional shift toward drought tolerance (or shade intolerance) was more prominent in colder plots (the lower panel of Fig. 6d). This was attributable to the trends that rising CO 2 levels were also related to the reduced relative abundance of drought-intolerant Picea spp. (Table S2) and increased abundance of droughttolerant/shade-intolerant Pinus spp. and early successional Betula spp. (Table S2) in colder plots (Fig.  S5b). Although there was also a negative relationship between rising CO 2 levels and the relative abundance of Picea spp. in humid plots (Fig. S5a), their association was weaker. Thus it was not translated to ecologically meaningful trends in community-level functional shifts. Speci cally, our study enhanced previous ndings by showing that drought-tolerant and early-successional (resource acquisitive or fast-growing species 3,21 ) might have bene tted more from rising CO 2 levels in colder plots, which was likely due to the extended growing season.
This might have augmented the growth of shade-intolerant (i.e., early successional) species in colder plots; however, in turn, it might have also facilitated mortality 42,43 , leading to faster changes in composition toward drought tolerance/shade intolerance. In addition to rising CO 2 , other environmental drivers, such as nitrogen deposition, might have also in uenced functional shifts. However, in most parts of Canada, nitrogen deposition occurs at low levels 44 , thus it is not very likely to drive signi cant changes in composition. Therefore, it is likely that rising CO 2 was at least partially responsible for the baseline climate-dependent shifts in functional composition over the study period. Nevertheless, other factors may have also contributed to the observed functional shifts.
Across the study areas, warming and temporal changes in water availability had negligible effects on functional composition ( Fig. S4; average effects are also shown as black lines in Figs. 6b, e). Temporal increases in temperature were associated with the higher relative abundance of drought intolerant Abies spp. and with decreased Betula spp. in humid plots (Fig. S5a), and with increased Betula spp. in warmer plots. However, warming had no association with both types of community-level functional composition regardless of the baseline climate. Our results contradicted a previous study, which showed that warming had positive effects on the growth of Picea mariana, particularly at higher latitudes (i.e., colder areas) in Eastern Canada 13 . This suggested that such trends may not be pervasive at larger scales across Canada, but rather regionally speci c.
Temporal changes in water availability consistently had no relationship with both types of functional composition throughout the baseline climate gradient (Figs. 6c, f). Although a decrease in the relative abundance of Picea spp. (with the temporal reduction in water availability) was more prominent in historically drier plots (Fig. S5a), and an increase in the relative abundance of Pinus spp. was greater in colder plots (Fig. S5b), these changes were not translated to shifts in functional composition. Our result was consistent with a previous study of the boreal forests of Western Canada, as it also showed no signi cant in uence of temporal variations in water availability to life history-based composition 15 .
Previous studies in temperate and boreal forests at the regional scale revealed that the functional composition of forests shifted toward fast-growing and drought-tolerant identity (or early-successional and deciduous traits) with rising CO 2 levels, increasing temperatures, or decreased water availability 9,15,30 . However, we found clear patterns in functional composition as relates to the baseline climate: colder plots experienced more rapid functional shifts toward drought tolerance under rising CO 2 levels than in warmer plots. Moreover, we also found that the functional shifts toward deciduous broadleaved trees, as well as the regional trends observed in previous studies 29,34 , were persistent and pervasive across biomes in North America. Thereby, our study generalized the ndings from regional observations 3,4,7,9,15 to larger spatial networks across multiple biomes (i.e., boreal forest, temperate broadleaved and mixed forests, and temperate coniferous forests 45 ) in North America.
While global environmental change is anticipated to intensify, our study suggests the vulnerability of the composition (higher rates in compositional changes from their original states 46 ) of colder plots after experiencing, and as a result of, these environmental changes). As these compositional shifts are likely to impact the functioning of forest ecosystems (e.g., net changes in biomass through growth and mortality) 9 by altering their functional identity 47 , our ndings of baseline climate-dependent functional shifts may partially elucidate the spatial variations in the impacts of global environmental on forest ecosystem functions 13,14,16 . Since such macro-scale interactions remain largely elusive and need exploratory analysis, the applications of arti cial intelligence can be helpful for pattern discovery 48,49 .
These environmental change-induced functional shifts could have signi cant impacts on temperate and boreal forests. Speci cally, greater increases in the capacity for resource acquisition (or earlysuccessional functional identity) might consequently be translated to increased productivity and mortality 42,43,50 , while increases in drought-tolerant abilities could result in reduced productivity and mortality in the face of changes brought about by global warming 9,50,51 . To ensure the sustainable functioning of forest ecosystems, future ecosystem management strategies should consider spatial differences in the response of forest composition to global environmental change, with a particular emphasis on colder forests experiencing higher rates of warming and compositional changes.

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 t the following criteria. The plots must: (i) be unmanaged, with a  (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).
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 rst (CWM PC1 , explained 60% of the variation) and second axes (CWM PC2 , explained 22% of the variation) as proxies for functional composition. The CWM PC1 collectively represented traits associated with deciduous broadleaved trees and higher resource acquisition 18,21,23,24 , being positively related with CWM Nmass , CWM Pmass , CWM SLA , CWM Habit (i.e., deciduous), CWM Struct (i.e., broadleaves), and CWM WD . On the contrary, the CWM PC2 represented traits associated with environmental tolerance, being negatively associated with CWM DT and positively related with CWM ST (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 CO 2 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 de ned 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 (CMI ave ) and MAT (MAT ave ) for each plot, as proxies for site-speci c 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: (CWM PC1 ) ijkl or (CWM PC2 ) ijkl = β 0 + β 1 × (PS) j + β 2 × (Prov) j + β 3 × f(SA) ij + β 4 × (Year) i + β 5 × (CMI ave ) j + β 6 × (MAT ave ) j + β 7 × f(SA) ij × (Year) i + β 8 × f(SA) ij × (CMI ave ) j + β 9 × f(SA) ij × (MAT ave ) j + β 10 ×(Year) ij × (CMI ave ) j + β 11 × (Year) ij × (MAT ave ) j + β 12 × (CMI ave ) j × (MAT ave ) j +π j + ε (1) where i and j are ith census and jth plot; CWM PC1 and CWM PC2 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 in uences 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 (p j ) to take locally unique conditions (site-specific disturbance histories; e.g., short-term climate events, insect outbreaks, non-catastrophic small re/wind/ ooding 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 in ation factor (VIF) was 2.77 for the CWM PC1 model and 2.79 for the CWM PC2 model, indicating that multicollinearity was not an issue.
As the measurement interval (years) varied between censuses, we performed variation partitioning involving those xed variables above and distance-based Moran's eigenvector maps (dbMEM) 36 to explicitly factor out the in uences 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 CWM PC1 model and 11 dbMEMs for the CWM PC2 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 signi cant temporal autocorrelation (examined by autocorrelation function estimation using the acf function in the stats package). Note that the coe cient 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 con rmed 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 coe cient comparison.
We also examined the temporal trends in atmospheric CO 2 concentrations, ATA, and ACMIA, and how they were associated with the CMI ave and MAT ave via linear xed effects models. To investigate the associations between the CWMs and rates of global environmental change, we replaced baseline climate variables (CMI ave and MAT ave ) in eqn. 1 with temporal change rates in ACMIA (ACMIA ChangeRate , cm/yr;  where CO 2 , ATA, and ACMIA are the atmospheric CO 2 concentration, anomalies of mean annual temperature, and climate moisture index, with which the CO 2 and ATA in our data was positively correlated (r 2 = 0.19). As the maximum VIF in this model was = 4.33 for the CWM PC1 model and 4.34 for the CWM PC2 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 coe cient 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 R 2 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 (m 2 /ha) by major tree genus. Similar to a previous study 19 , major tree genus was de ned 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 xed-effects parameters used in eqn. 1, as well as the three global environmental change drivers with the same xed effects used in eqn. 2.
For the interpretation of all analyses, we focused on not the statistical signi cance (i.e., p-value) but the 'ecological signi cance'; 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 signi cant 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 signi cance level, we advocate that statistical signi cance does not necessarily equal to ecological     Temporal and spatial trends in functional composition. The main effects of stand age (a), the long-term average of climate moisture index (CMIave, b), and the long-term average of mean annual temperature (MATave, c) on community-weighted mean of trait values (CWMPC1 and -CWMPC2). CWMPC1 is a functional composition associated with deciduous broadleaved trees (higher value) vs conifers (lower value), while CWMPC2 is related to environmental tolerance (higher value = higher drought tolerance (DT); see Methods). Dots and error bars re ect the means and their 95% Bayesian con dence intervals. Blue lines are tted main effects with their 95% Bayesian con dence intervals shown as shaded areas. Based on AIC, stand age was transformed by the squared root.

Figure 5
Temporal trends in community weighted-mean of traits associated with deciduous broadleaved trees vs conifers (CWMPC1) and drought tolerance (-CWMPC2). Trends dependent on the long-term averages of the climate moisture index (CMIave, a) and mean annual temperature (MATave, b). Values are means and their 95% Bayesian con dence intervals. CMIave and MATave were binned from 1.0 to 104.8 (cm) and from -1.6 to 7.1 (°C) (their 5th and 95th percentiles) for four levels.