Forest types outpaced tree species in centroid-based range shifts under global change

Discussion: Our ﬁ ndings provide an important scienti ﬁ c basis for adaptive forest management and conservation, which primarily depend on individual species assessment, in mitigating the impacts of rapid forest transformation under climate change.


Introduction
Treesexcept in the seed life stageare individually immobile organisms, but tree species worldwide are found to undergo substantial changes in geographic distributions through altered dispersal, growth, and survival patterns under global change.Drivers of species range shifts include climate change, associated threats such as insect and disease outbreaks (Anderegg et al., 2015), and other anthropogenic activities (Vanderwel and Purves, 2014).In North America, while consistent upward shifts have been widely observed in recent years along an altitudinal gradient (Kelly and Goulden, 2008;Lenoir et al., 2008;Chen et al., 2011), more diverse patterns have been observed in horizontal shifts in the latitudinal and longitudinal dimensions; Some tree species have moved to higher latitudes to track increasing temperature (Woodall et al., 2009(Woodall et al., , 2010;;Boisvert-Marsh et al., 2014;Sittaro et al., 2017), a smaller proportion of species have moved southward (Woodall et al., 2010;Zhu et al., 2012), while many species have moved longitudinally in response to changes in precipitation patterns (Fei et al., 2017).Little is known, however, about whether and how this substantial shifting and reshuffling of tree species ranges may cause an overall type of foresta distinctive assemblage of tree species distributed across a wide geographical extentto shift its range.
Since certain tree species are often found together in forest communities, various assemblages of tree species imply different types of forest communities which support different types of plants, wildlife (Perry et al., 2022), and microbiomes (Steidinger et al., 2019).A forest type represents a distinctive assemblage of tree species distributed across a wide geographical range (Perry et al., 2022), which in some studies is also referred to as a forest region (Rowe, 1972;Dyer, 2006), a tree species assemblage (Costanza et al., 2017), or a forest community (Knott et al., 2020).The classification of forest types and quantification of their dynamics thus provide important references for forest management, conservation, climatechange mitigation, and restoration (Perry et al., 2022).That said, a lack of consistent classification of forest types at a continental scale has been a major obstacle to the understanding of the patterns of forest type range shifts.For over a century, forests have been classified based on tree species composition and structural characteristics (Küchler, 1964;Rowe, 1972;Eyre, 1980;Dyer, 2006;Ruefenacht et al., 2008;Costanza et al., 2017;Knott et al., 2020).With the recent advancement of forest data availability (Liang and Gamarra, 2020) and computational capacity, new data-driven forest type classification schemes minimize subjective biases and exhibit great potential to replace conventional approaches (Dyer, 2006;Costanza et al., 2017;Knott et al., 2020).With statistical methods like cluster analysis (Dyer, 2006;Costanza et al., 2017) and allocation models (Knott et al., 2020), data-driven classification is capable of effectively identifying forests with similar characteristics of species relative abundance and dominance.
Although past research has significantly advanced our understanding of tree species range shifts (Kelly and Goulden, 2008;Lenoir et al., 2008;Woodall et al., 2009Woodall et al., , 2010;;Zhu et al., 2012;Boisvert-Marsh et al., 2014;Fei et al., 2017;Sittaro et al., 2017) and patterns of tree fecundity and recruitment (Sharma et al., 2022), the range shifts of forest types and how they differ from those of tree species remain largely unknown.In community ecology, Gleason's individualistic viewpoint has gained more acceptance than Clements' organismic concept of plant communities (Clements, 1916;Gleason, 1926).Gleason (1926) states that species distributions are determined by the unique ecological niches and requirements of each species, and the assemblages they form are merely a result of chance.Therefore, communities do not behave as "superorganisms," as Clements (1916) contends, and could be highly dynamic (Guisan and Zimmermann, 2000;Cavender-Bares et al., 2009).As further exemplified in Markowitz's portfolio theory of investment (Markowitz, 1952), the change of an ensemble can differ from the changes of its constituents.An ecological hypothesis (Chesson et al., 2005;Schindler et al., 2015;Hui et al., 2017), derived from the portfolio theory (hereafter, portfolio hypothesis), postulates that the dynamics of an ecological community are not the same as the dynamics of its constituent species.According to the portfolio hypothesis, the speed and direction of forest type range shifts can differ substantially from the average speed and direction of tree species range shifts.However, this hypothesis has never been formulated and tested for forest type dynamics.
Past studies widely vary in their methods to quantify range shifts (Iverson and McKenzie, 2013).With similar and sufficient sample distributions between the two time periods across the study area (e.g., Boisvert-Marsh et al., 2014;Fei et al., 2017), range shifts can be quantified without interpolation or extrapolation.However, interpolation or extrapolation is necessary for studies with temporally differing sample distributions (e.g., Iverson et al., 2019a).In this case, correlational species distribution models are useful for comparing species ranges over time (Iverson and McKenzie., 2013).While the former approach provides observational, and thus more realistic, patterns of range shifts, it is often limited in terms of spatial resolution (e.g., Fei et al., 2017) and area covered (e.g., Boisvert-Marsh et al., 2014).In contrast, the latter approach provides larger spatial and temporal scales, but it is subject to uncertainties associated with interpolation or extrapolation (e.g., sampling bias, uncertainty in predictor variables, model assumptions, incorporation of dispersal ability, etc.; see Araujo et al. (2019) and Beale and Lennon (2012) for details).Furthermore, range shifts can be represented by shifts in the centroids (e.g., Fei et al., 2017;Iverson et al., 2019a) or boundaries (e.g., Boisvert-Marsh et al., 2014).In the context of forest type and its constituents (i.e., species), forest type boundary shifts simply depict the minimum of each constituent species' boundary shift.As such, they will often reflect the slowest species' range boundary shift.Instead, the range centroid of a forest type reflects the distribution of its most representative and abundant/dominant species, capturing the primary function and service of this forest type.
Here, we formulated and tested the portfolio hypothesis on the divergent range shift rates between forest types and constituent tree species between 1970-1999 and 2000-2019.To do this, we collated a continental-scale forest inventory database containing in situ measurements of more than 20 million trees in 596,282 sample plots located across North America.First, we used autoencoder neural networks and K-means cluster analysis to group the data based on similarities in tree species composition for each period.By modifying established classification algorithms to accommodate larger datasets, the resulting forest types were comparable to those described by Costanza et al. (2017) and Dyer (2006).In parallel, we developed correlative species distribution models to estimate the range shift of each species based on geographic centroids.Finally, based on the portfolio hypothesis, we compared the speed of forest type range shifts and the average of their constituent species.

Materials and methods
We first formulated the portfolio hypothesis, describing the theoretical relationship between forest type range shift and the average range shift of its constituent species.To test this hypothesis, we prepared tree data with an importance value index (IVI) as a proxy for species relative abundance for periods 1970-1999 and 2000-2019 (Figure 1).We created interpolated maps of species distribution ranges using a species distribution model based on predictor variables encompassing climate, topography, soil, and anthropogenic characteristics.In parallel, we classified forest types using the same tree data, which were further mapped across the study area.Finally, for each pair of forest types identified in both periods, we calculated the forest type range shift and the average range shift of its constituent species based on species composition (Figure 1).Species range shift was represented by the geographic centroid shift of its interpolated distribution range.

The portfolio hypothesis of range shifts
We derived the following hypothesis from Markowitz's portfolio theory of investment (Markowitz, 1952) to quantify the difference in range shifts between forest types and constituent tree species.Tree species range shift was represented by the geographic centroid shift of its modeled distribution range.
Let a i (x) be the relative abundance of species i at spatial location x (demarcated by coordinates of latitude and longitude).The cumulative relative abundance of species i over space x is ox a i (x) = A i , and the geographic centroid of species i's range can be calculated as the weighted mean of x, c i = ox x • a i (x)=A i .The relative abundance of a forest type at location x is the sum of its constituent species' relative abundance at the location, a(x) = oi∈G a i (x), where G = (1,...,n) represents the set of all species of this forest type.The cumulative relative abundance of this forest type can be computed as ox a(x) = A; let p i = A i =A, and the Overview of the methods.Tree data was prepared for each period at the 0.025°grid level (approximately 3km) with species importance value index (IVI).Using the tree data, we built species distribution models to create interpolated maps of species distribution ranges and classified forest types.For each pair of forest types that were identified in both periods, we estimated the range shift of the forest type and the average of constituent species based on the portfolio hypothesis.Species range shift was represented by the geographic centroid shift of the interpolated distributions range.Abbasi et al. 10.3389/fevo.2024.1366568Frontiers in Ecology and Evolution frontiersin.orggeographic centroid of the forest type can be calculated as c = ox x • a(x)=A = oi∈G c i • p i .By the Leibniz product rule of calculus (see details in Supplementary Text 1.1), we have Because oi∈G p i = 1 and thus oi∈G Dp i = 0, the latter term of Equation 1 equals the covariance between species' centroids and the change in their cumulative relative abundance: cov(c i , Notably, the shift in a forest type's centroid (Dc) is driven not only by the weighted mean of centroid shifts at the species level ( oi∈G Dc i • p i ), but also, counterintuitively, by the covariance (cov(c i , Dp i )) that inflates or deflates the centroid shift of the forest type (Figure 2).This covariance term, which we call the portfolio effect, can be attributed to a difference between forest type range shift and the range shift of its constituent tree species.We therefore made the following hypothesis: when the range shift of a forest type and its constituent tree species are measured along a given direction, a positive covariance (cov(c i , Dp i )) will inflate the magnitude and velocity of the forest type range shift along this direction, while a negative covariance will reduce the magnitude and velocity along this direction (Figure 2).In most cases, portfolio effects are expected to be positive because species at the front of a forest type range (i.e., in the direction of forest type range shift) often have increasing abundances due to preferential allocation to dispersal and reproduction, while species at the rear of forest type range are decreasing in abundance, reflecting the compounded effect of constituent species' mortality and recruitment on the centroid shift (Burton et al., 2010;Rumpf et al., 2018).

Data integration
For this study, we compiled and integrated in situ forest-tree data from independent and standard forest inventories.Data for the United States came from the Forest Inventory and Analysis (FIA) (Burrill et al., 2021) and the Cooperative Alaska Forest Inventory (CAFI) (Malone et al., 2009).Data for Canada came from two independent sources: permanent sample plot networks (Paquette and Messier, 2011;Chen et al., 2016) and Canada's National Forest Inventory ground plot network (National Forest Inventory, 2011;Zhang et al., 2017).See Supplementary Text 1.2 for a detailed description of each data.
We derived the following data integration protocol to harmonize the different forest inventory datasets described above into consistent continental data frames.From each dataset, we obtained tree-level information for all the trees with a minimum diameter at breast height (DBH) of 1 cm.We grouped these treelevel records by the year of inventory and compiled one data frame for 2000-2019 and another data frame for 1970-1999.For each period, we summarized tree-level information into a plot-level species abundance matrix.We calculated, for each sample plot, the importance value index (IVI hereafter) for a species, which is the sum of the percent number of stems and the percent basal area for the species.Frequently used in forestry research as a typical measure of species abundance (Dyer, 2006;Costanza et al., 2017;Iverson et al., 2019a), IVI equally weighs the number of stems and basal area of a particular species, and ranges from 0 to 200.
The final continental data frames consisted of plot identification and coordinates, as well as the IVIs of all tree species present on each plot.The plots were widely distributed across the forested areas of the continent (Supplementary Figure 1).For the 1970-1999 data frame, because some trees in the genera of Aesculus, Amelanchier, Carya, Crataegus, Halesia, Malus, and Salix were recorded only to the genus level, we also calculated the IVIs of these genera (Supplementary Table 1).To mitigate the impact of sampling bias (Araujo et al., 2019), we derived the average species IVI of all plots located within a 0.025°by 0.025°(approximately 3 by 3 km) grid cell in each past and present plot-level dataset, which is a reasonable aggregation regardless of the distribution of species IVI (see Supplementary Text 1.3).
We also compiled 38 predictor variables for the supervised learning of species and forest type distributions, consisting of 17 climate variables (Fick and Hijmans, 2017;Karger et al., 2017 Range shift dynamics of forest type and constituent tree species based on the portfolio effect.Forest type shift is controlled by constituent species centroid shifts (Δc i ) and changes in cumulative relative abundance (Dp i ).The covariance between a species' centroid and the change in its cumulative relative abundance (cov(c i , Δp i )) determines whether forest type shift outpaces or lags tree species shifts.(A) In this simple scenario where tree species and forest type shift along the same direction, if covariance > 0, forest type shift outpaces tree species shifts, but (B) if covariance< 0, forest shift is outpaced by tree species shifts.et al. 10.3389/fevo.2024.1366568Frontiers in Ecology and Evolution frontiersin.orgTrabucco and Zomer, 2019;Karger et al., 2022), 13 topographic variables (Amatulli et al., 2018), seven soil variables (Batjes, 2016), and human footprint (Venter et al., 2016).These predictor variables were derived from open-access satellite-based remote sensing and ground-based survey data layers, all of which have a nominal resolution of 1 km.Detailed information on the predictor variables is available in Supplementary Table 2. Different climate and human footprint data were used for past and present periods to match the represented years (see Supplementary Table 2).We extracted predictor variables to the centroid of each grid cell using the "sf" or "raster" packages (Pebesma, 2018;Hijmans, 2023;Hijmans and Van Etten, 2022) in R (R Core Team, 2021) and used the "Hmisc" package in R to impute missing data in those predictor variables (Harrell, 2023).Our final training data encompassed more than 114,000 grid cells for the present period and 97,000 grid cells for the past period containing species IVI and 38 predictor variables.

Abbasi
For mapping purposes, we prepared another 0.025°by 0.025°g rid, covering forest extent across North America with all predictor variables.Each grid cell had a minimum 10% canopy cover based on the global forest range map (Hansen et al., 2013), in accordance with FAO's definition of "forest" (FAO, 2020).Species or forest type distributions were predicted for each grid cell of this mapping grid.Our study region encompassed 1,004,358 grid cells of forested area across North America, with a total of ~5 million km 2 .The tropical regions of North America, i.e., Mexico, Central America, and the Caribbean, were not included in this analysis due to a lack of remeasured in situ data.
Our study region covered 92 terrestrial ecoregions (Olson et al., 2001) across the United States and Canada.These ecoregions were grouped into three distinct arch-biomes based on the definition of the eastern United States in previous studies (e.g., Knott et al., 2020) and biome map (Olson et al., 2001): West (39 ecoregions), East (33 ecoregions), and Boreal (20 ecoregions, Supplementary Figure 1).For each arch-biome and time frame (2000-2019 and 1970-1999), we classified forest type (see 2.4 Forest type classification) and quantified distribution ranges of forest types and constituent tree species separately.

Range shifts of tree species
Following the formulation of the portfolio hypothesis, the first step to quantify forest type range shifts is to determine constituent species' range shift patterns, which consisted of two steps: creating a spatially continuous map of species i's relative abundance (a i (x)) at location x (i.e., grid latitude and longitude) across the continent and quantifying its geographic centroid shift (Dc i ).To estimate a species' relative abundance, we first created interpolated maps of species IVI across the 4.9 million-km 2 study region using random forests model and 38 predictor variables, which took similar steps as correlative species distribution models (Araujo et al., 2019;Iverson et al., 2019a).Since we do not have sufficient samples in some areas, we used interpolated relative abundance map instead of observations.For each arch-biome (West, East, and Boreal) and time frame (2000-2019 and 1970-1999), only species with sufficient sample size (≥60 grid cells where species IVI > 0) in both time periods were included (van Proosdij et al., 2016; Supplementary Table 1).
Random forests are a non-parametric ensemble learning approach (Breiman, 2001), which combines a variant of decision trees and an additional level of randomness by bootstrapping subdata and different sets of predictor variables.Random forests are commonly used in species distribution models (Araujo et al., 2019) as it mitigates the multicollinearity issues that most statistical models face (James et al., 2013).We used the "randomForest" package in R (version 4.0.4)(Liaw and Wiener, 2002) to train a separate model for each species, period, and arch-biome.Following Iverson et al. (2019a), we reported the mean predicted IVI of all decision trees for each species or zero for species with zero median and a coefficient of variation no less than 2.75 among all predicted values of decision trees.For each species, we built 20 random forests models to calculate an average IVI in each grid cell in the mapping grid.
Based on the estimated species IVI, we calculated species relative abundance (a i (x)) along each latitudinal and longitudinal gradient by calculating percent IVI for each species in each grid cell.We then calculated the cumulative relative abundance (A i = ox a i (x)) and geographic centroid (c i = ox x • a i (x)=A i ) for each species.The direction and velocity of species range shift (Dc i ) were calculated based on the displacement between the past and present geographic centroids using the "sp" and "sfsmisc" packages in R (Pebesma and Bivand, 2005;Maechler et al., 2023).In this study, the direction was measured in azimuth, the angle between past and present geographic centroids around the same horizon (i.e., altitude was not considered), ranging from 0 to 360°measured from the North direction.We also determined the total area of each species' range as the sum of the grid cell area weighted by the species' relative abundance.Grid cell area was estimated using the "raster" package in R (Hijmans and Van Etten, 2022).

Forest type classification
We modified the existing forest type classification algorithm (Dyer, 2006;Costanza et al., 2017) to determine forest type definition in each arch-biome and period.Specifically, instead of hierarchical cluster analysis used in previous research (Dyer, 2006;Costanza et al., 2017), we used K-means cluster analysis due to the size of the data.We first used an autoencoder neural network to calculate a latent space representation of the original grid tree IVI data, which is often more suitable for K-means cluster analysis than the original data (Song et al., 2014).Instead of the interpolated maps of species IVI, we used the grid-level observational data to maximize the performance of cluster analysis.See Supplementary Text 1.4 for a detailed description of autoencoder neural networks and the architecture used in this study.To avoid potential bias caused by insufficient sample sizes, we excluded the species that are present in less than 60 grid cells (where IVI > 0) (Supplementary Table 1).
Based on the reduced dimensional representation, we then conducted a K-means cluster analysis using the built-in function "kmeans" in R (R Core Team, 2021).We set the number of starts to 50 and the maximum iterations to 100.To fine-tune hyperparameters for the autoencoder neural networks (the number of dimensions) and K-means cluster analysis (the number of clusters; Supplementary Figure 2), we calculated the silhouette width using the "cluster" package in R (Maechler, 2018).Silhouette width is an indicator of between-cluster heterogeneity; With a range between −1 and 1, positive silhouette width values indicate that a given member of a cluster is closer to its own cluster's centroid than to the nearest cluster's centroid.Negative values indicate that a given member is closer to the nearest cluster's centroid than to the centroid of its own cluster.Generally, higher silhouette width values indicate greater between-cluster heterogeneity.
With the optimal hyperparameter values, we conducted K-means cluster analysis 20 times to derive the final classification results due to the random nature of the analysis.Since forest types (i.e., clusters) were defined independently in each of the 20 repetitions, we manually matched the same forest type based on similarities in mean species IVI (i.e., species composition for each forest type) between all the combinations of forest types generated from the 20 repetitions.When 10 or more repetitions identified the given forest type, we recognized the forest type as a final forest type.Moreover, since we classified forest types for three arch-biomes separately (West, East, and Boreal), there were potential overlaps of forest types between them.To identify and merge potential overlaps, we calculated the Euclidean distance of all combinations of the final forest types in terms of mean species IVI.If an Euclidean distance was less than 60, forest types across arch-biomes were merged.One exception was that western aspen-mixed conifer (W-J) and boreal quaking aspen-spruce (B-A) remained separated due to the large expanse of quaking aspen (Populus tremuloides).
Finally, with the final set of forest types for each period, we identified forest types that were present in both periods.To do this, we calculated the Euclidean distance of all combinations between past and present forest types in terms of mean species IVI.Pairs were considered matching when the forest type of minimum distance was the same between the past-and-present pair.For example, if and only if present forest type X's closest past forest type is Y, and past forest type Y's closest present forest type is also X, they were considered matching.Finally, we manually grouped final forest types into larger forest biome categories based on Eyre (1980).

Forest type mapping
We predicted the distribution ranges of forest types by first considering two candidate imputation models: random forests and support-vector machines.Support-vector machines are supervised learning models which construct a hyperplane or set of hyperplanes in a high-or infinite-dimensional space to help analyze data for classification and regression analysis (Vanpik and Cortes, 1995).We used the "e1071" package in R with the default hyperparameter setting (Meyer et al., 2019).For random forests, we used the default hyperparameter setting of the "randomForest" package in R (Liaw and Wiener, 2002) and the same set of predictor variables as described in 2.2.Data integration.
To assess the performance of the imputation model in mapping forest types across the continent, we conducted a rigorous 80/20 cross-validation using bootstrapping.In each iteration, we used stratified sampling to split the entire dataset into the training (80%) and testing (20%) sets and conducted the combination of undersampling and oversampling of the training set for both random forests and support-vector machines to balance the sample size for all classes (i.e., forest types).Stratified sampling was conducted using the "caret" package in R (Kuhn et al., 2020), and undersampling and oversampling were conducted using the "UBL" package (Branco et al., 2016).Based on five random iterations with sample replacement in each of the 20 repetitions, we calculated the 95% confidence interval of classification accuracy, the Kappa statistic, and elements of the confusion matrix.For each candidate imputation model, the output was a matrix of class probability from five iterations.We chose the forest type of majority vote from the five iterations, and thus, our final output was a matrix of class probability from the 20 repetitions.Based on how many repetitions, out of 20, returned the given forest type, we calculated percent forest type in each grid cell.Every grid cell has a sum of 100 for all forest types, with the forest type with the highest percentage being the final assigned forest type.(e.g., The final forest type assigned to a grid cell consisting of 20% forest type A, 30% forest type B, and 50% forest type C is forest type C).The percent forest type was used to visualize the geographic patterns of forest type range shifts, specieslevel range shifts, and the difference between them.

Range shifts of forest types
Following the formulation of the portfolio hypothesis, we calculated the weighted sum of species-level shifts ( oi∈G Dc i • p i ), the covariance (the portfolio effect, cov(c i , Dp i )), and the sum of the two (i.e., forest type-level shift (Dc)) for each forest type that was present in both periods along the latitudinal and longitudinal gradient.For each forest type, G represents the set of all species present in the forest type (mean species IVI > 0).Along each respective gradient, Dc i represents the centroid shift of species i, D p i represents the change in p i from past to present periods, and p i and c i take the mean values of past and present periods.Thus, we have the species-level shifts, covariance, and forest type-level shift for each latitudinal and longitudinal gradient for each forest type, where the weighted sum of species range shifts ( oi∈G Dc i • p i ) and covariance term (cov(c i , Dp i )) precisely matches a forest type range shift (Dc).
In a two-dimensional space with latitude and longitude combined, we also derived the velocity of a forest type range shift in kilometers per decade and azimuth angle by calculating the distance between past and present forest type centroids using "sp" and "sfsmisc" packages in R (Pebesma and Bivand, 2005;Maechler et al., 2023).This forest type range shift can be visualized as a vector in a two-dimensional space of latitude and longitude, as a resultant of two composing vectors: species-level vector and covariance vector.Quantifying the length of these two vectors (i.e., velocity) is not practical due to the earth's curvature, yet we aimed to approximate it by calculating the distance between past forest type centroid and the point to where the vector heads, both in degrees and kilometers (see Supplementary Figure 3 for a visualized example).Finally, to expand the velocity of forest type and species range shifts to a grid level for visualization purposes, we weighted the velocity of each forest type by the percent forest type in each grid cell.Average percentage of past and present forest types was used.

Results
Our forest type classification identified eight forest biomes and 43 forest types across North America (Figure 3, Supplementary Figures 4-11, Supplementary Table 3).The mean silhouette widths from our Kmeans cluster analyses were significantly greater than zero for all forest types (p< 0.001) in the West for the present dataset (Supplementary Table 4).Eighteen out of 19 forest types in the West, 22 out of 26 in the East, and all six forest types in the Boreal region were significantly greater than zero in the mean silhouette width.In summary, 90% of the forest types classified here were significantly distinct from one another in terms of species composition (Supplementary Table 4).
For mapping the ranges of classified forest types, the random forests model was 10-17% and 11-20% more accurate in terms of overall accuracy and the Kappa statistic, respectively, compared with the support vector machine model (Supplementary Figure 12).Therefore, we selected random forests as the final imputation model.The confusion matrices based on random forests models were based on the number of cases in class prediction, standardized in percentage (Supplementary Figures 13, 14).For the present dataset, the coastal redwood-tanoak forest (W-P) had the highest classification accuracy (88%; Supplementary Figure 13), and the red maple-hardwood forest (E-F) had the lowest one (18%; Supplementary Figure 13) among all forest types.
At the individual species level, Sitka spruce (Picea sitchensis) had the greatest velocity in range shifts (480.4 km•decade -1 ), followed by balsam fir (Abies balsamea; 438.6 km•decade -1 ) and gray alder (Alnus incana; 354.6 km•decade -1 ) (Figure 4A, Supplementary Table 1).In contrast, pond cypress (Taxodium ascendens) had the lowest velocity We classified forested areas across North America into 43 forest types in eight forest biomes and three arch-biomes (East, West, and Boreal) following the existing algorithms.See Supplementary Table 3 for the definition of forest types and constituent tree species.Colors in the circular dendrogram corresponds to those in the map of forest type distribution range.Here, we show a map of forest types for the present period, based on the distribution of forest types estimated by random forests (see 2.5 Forest type mapping).

B A FIGURE 4
The velocity of (A) tree species range shifts and (B) forest type range shifts.Velocity of shifts (here in logarithmic scale) is defined as the distance between past and present centroids of range in kilometers per decade (km•decade -1 ).See Supplementary Table 3 for the definition of forest types and constituent tree species. of all (1.5 km•decade -1 ), followed by water-elm (Planera aquatica; 4.5 km•decade -1 ) and sweetgum (Liquidambar styraciflua; 4.5 km•decade -1 ).The velocities of temperate tree species range shifts found in our study are consistent with those reported in the previous studies (Supplementary Table 5).Few boreal tree species have ever been assessed in terms of shifting velocity, and here, we found that boreal tree species' ranges are shifting much faster than temperate ones.In terms of direction, 36 out of 150 species shifted northwards, 34 eastwards, 27 southwards, and 53 westwards during the studied period (Supplementary Table 1).
At the forest type level, we found that the range of Sitka sprucewestern hemlock forest (W-A) shifted with the highest velocity at 327.8 km•decade -1 (Figure 4B).Among the top six fast-shifting forest types, three were in the boreal forest biome (B-A, B-B, and B-D, see Supplementary Table 3 for forest type names), two were in the eastern mixed forest biome (E-A and E-K), and one was in the Pacific-coastal forest biome (W-B) (Supplementary Table 4).The remaining forest types shifted at a velocity lower than 100 km•decade -1 .In terms of the direction of shift, nine out of 43 forest types shifted westwards, 16 eastwards, 11 southwards, and seven northwards during the studied period (Supplementary Table 4).
Overall, forest type range shifts differed substantially from tree species range shifts (i.e., the average range shift of constituent tree species) in terms of velocity (Figures 5, 6).On average, forest type ranges shifted at 86.5 km•decade -1 , more than three times as fast as the weighted average of their constituent tree species (28.8 km•decade -1 ) across the continent at the grid level (Figures 5B, 6C).For more than 75% of forest types, the range shifts at the forest type level substantially outpaced the average range shifts of constituent tree species, and only 10 out of 43 forest types moved more slowly than their constituent tree species in terms of range shifts (Figure 6C, Supplementary Table 4).In the boreal and Great Lakes regions, the velocity of forest type range shifts exceeded that of tree species range shifts by 200 km•decade -1 or more (Figure 5B).Along the Rocky and Appalachian Mountains, forest type ranges shifted with a lower velocity than the ranges of their constituent tree species (Figure 5B).

Discussion
While existing forest type classifications for North America are limited to either Canada or the United States (Rowe, 1972;Dyer, 2006;Ruefenacht et al., 2008;Costanza et al., 2017;Knott et al., 2020), our study supports ongoing and future international collaborations in forest management and climate actions (FAO and UNEP, 2020), with a consistent continental-wide data-driven forest type classification scheme.Nevertheless, our classified forest types are generally compatible with existing forest type classifications for the North American continent (Rowe, 1972;Dyer, 2006;Ruefenacht et al., 2008;Costanza et al., 2017;Knott et al., 2020).
Our findings provide strong support for the portfolio hypothesis.The marked difference in the range shift velocity between forest type and constituent tree species is mainly attributable to the covariance term (cov(c i , Dp i )), which we call the portfolio effect.Positive portfolio effects provide a mechanistic underpinning of our finding that, across the North American continent, forest type range shifts generally outpaced tree species range shifts (Figures 5B, 6C).In some cases, negative covariance terms pervade, leading to reduced forest type range shift, such as in some forests across the Rocky and Appalachian Mountains (Figure 5B).In these montane forest types, constituent species increase their relative abundances in the rear edge, where investments in competitive traits are likely to dominate, even if the species centroids drift towards the front edge (Figure 2).
Various factors in complex interactions are attributable to species abundance dynamics, leading to positive or negative portfolio effects.While climate change (e.g., increasing temperature and precipitation changes) is commonly considered the major driver of species abundance dynamics and associated range shifts (Boisvert-Marsh et al., 2014;Fei et al., 2017;Sittaro et al., 2017), climate change-induced disturbances (e.g., fires and disease outbreaks), human disturbances (e.g., logging and land use changes), and biotic interactions influencing the dispersal and The velocity of forest type range shift and the difference in velocity between forest types and the weighted mean of their constituent species across the continent, mapped at a 0.025°resolution.(A) The velocity of forest type at the grid level.(B) The difference in shift velocity between forest types and the weighted mean of their constituent species at the grid level, with warm colors representing areas where forest type range shifts outpaced species range shifts, and cold colors representing areas where species range shifts outpaced forest type range shifts.Forest type range shift and the weighted mean of its constituent species were first quantified for each forest type.Here, we used interpolated maps of forest types to visualize them based on percent forest type in each grid cell.et al. 10.3389/fevo.2024.1366568Frontiers in Ecology and Evolution frontiersin.orgsurvival rates can further exacerbate the dynamics of species abundance (Vanderwel and Purves, 2014;Anderegg et al., 2015;Gauthier et al., 2015;Brice et al., 2019Brice et al., , 2020;;Seidl et al., 2020).

Abbasi
Although it is out of the scope of this study to assess the determinants behind species abundance dynamics, we utilized non-climate predictor variables, including human footprint (Venter et al., 2016) (Supplementary Table 2), to account for the impacts of potential human disturbances.
There are profound impacts of forest type range shifts on forest biodiversity and associated ecosystem functioning, food, water, energy security (Gitay et al., 2002;Neilson et al., 2005;Wellstead and Howlett, 2017;Kremen and Merenlender, 2018), human wellbeing (IPBES, 2018), and socioeconomic value (Hanewinkel et al., 2013).Forest type range shifts can jeopardize the sustainability of local forest industries, making them more vulnerable to timber price fluctuations (Zhou, 2021).It can also inflate timber procurement ranges and increase transportation costs, causing significant downstream financial implications with serious welfare and economic consequences comparable to the impact of COVID-19 on transportation and logistics (Polzin and Choi, 2021).Furthermore, the collective human experiences of rural communities embedded within these forested landscapes have strong ties to surrounding forest types.From the Sitka spruce-western hemlock forests in the Pacific Northwest to the oak-pine forests along the Appalachians (Supplementary Table 4), the change in native forest ecosystems may be threatening the customs, identities, and culture of indigenous (Chamberlain et al., 2018) and other local communities while jeopardizing non-timber forest products supply and overall environmental justice (Fleetwood, 2020).The rapid shift of forest types may place urgency on human communities, especially rural populations, to adapt their cultural norms and relationships with surrounding forests.
Our finding that a majority of forest types shifted faster than their constituent tree species due to positive portfolio effects suggests that the impacts of tree species' range shifts under global change on forest ecosystem functioning and services may be underestimated.Forest ecosystem functioning (Loreau, 2000 The difference between forest type range shifts and tree species range shifts.Scatter plots show the velocity of forest type vs. tree species range shift velocity, in terms of latitudinal (A) and longitudinal (B) shift of geographic centroids, with positive values representing northward (A) and eastward (B) shifts and negative values southward (A) and westward (B) shifts.Vertical line segments represent the difference between the velocity of forest type range shifts and tree species range shifts, along latitude (A) and longitude (B).The length of these segments is identical to the portfolio effect.(C) A comparison between the velocity of forest type and tree species range shift velocity, regardless of the direction.The horizontal dotted lines represent the mean velocity of forest type and tree species shists in pink and blue, respectively.The mean velocity represents a grid-level velocity (i.e., Figure 5A for forest type).To ease identification of various forest types, axis text colors of the Panel (C) are consistent with Figures 3, 4B.Cardinale et al., 2011), productivity (Liang et al., 2016), phenology, and population turnover (Zhu et al., 2014;Pecl et al., 2017) are directly related to tree species composition and other forest characteristics (Loreau, 2000;Cardinale et al., 2011;Paquette and Messier, 2011;Liang et al., 2016).Therefore, existing adaptive forest management regimes, which are based primarily on individual species range projections and associated environmental and social aspects (Millar et al., 2007;Iverson et al., 2019b), may have underestimated the impacts of global change (including climate change, land use change, invasive species regimes, habitat fragmentation, and forest degradation).Compositional changes of species into novel forest types imply that forest managers may need to address effects on forest ecosystems that are unique in space and time (Iverson and McKenzie, 2013), pursuant to a need for more dynamic silvicultural practices.For instance, in the central United States, a diminishing supply of various white oak species, such as white oak (Quercus alba) and bur oak (Q.macrocarpa), is threatening the bourbon industry (Conrad et al., 2019), a staple of American culture and tradition.These oak species are supported by a variety of co-existing tree species; For example, deep roots of hickory species (Carya spp.) improve soil structure (Smith, 2024), eastern redbud (Cercis canadensis) enhance soil nutrient availability through nitrogen fixation, pine species (Pinus spp.) co-exist without intense competition (Nowacki and Abrams, 2008), and species like black cherry (Prunus serotina) support a wide range of wildlife species, which play a critical role in acorn dispersal.Shifts and contracting ranges of these co-existing species communities can negatively affect the quality and volume of oak species for local industry.As forest type-level change (25% reduction in Appalachian oak-pine forest range, Supplementary Table 4) was far more apparent than species-level change (e.g., 12.5% reduction in Q. alba range and 6% increase in Q. macrocarpa range, Supplementary Table 1), missing the forest for the trees for their range shift patterns may reduce the capability and preparedness of local forest industry and communities to face global change through disruptions in timber supply chains and ecosystem benefits.
Positive portfolio effects enable us to scale global change-induced range shifts up to the levels where decision-making processes often take place.Because of the portfolio effect, climate change and other factors of global change, such as anthropogenic effects, land use change, as well as associated insects/disease outbreaks and invasive species regimes, will result not only in the migration of forest types, but also in the emergence of new forest types (Williams and Jackson, 2007).Due to unequal changes in relative abundance among species, new assemblages of species (i.e., new forest types that are unknown today) may appear.To this end, additional research could outline the relative importance of various global change factors on forest migration and tree species range shift patterns.Further spatial analyses could also disentangle directional and non-directional factors influencing the absolute difference in range shift velocity between forest types and constituent tree species.
In this study, we quantified the velocity of range shifts using interpolated maps over a large geographic area to best assess the difference between forest type range shifts and constituent species range shifts.As such, some species in our results exhibited higher velocities than previously reported.This is especially prominent in the boreal region, where sample density is limited (Supplementary Figure 1) and no previous studies have ever assessed at this spatial scale.The contrast between our reported velocities and previous studies mainly results from three sources: uncertainty associated with interpolation, limitations in incorporating biotic factors into species distribution, and the spatial scale of this study.While estimating the distributions of species or similar entities through interpolation or extrapolation has been a popular method in ecology (Araujo et al., 2019), predictions are subject to various uncertainties (Beale and Lennon, 2012), which could affect the reported velocities.Nevertheless, we minimized the uncertainty associated with sampling bias through grid aggregation, model selection by considering several candidate models, and extrapolation by analyzing three arch-biomes separately.However, the incorporation of biotic interactions (competition, dispersal ability, etc.) remains an enormous challenge in correlative distribution modeling (Araujo et al., 2019), and thus, the interpolated maps do not reflect realized niche but rather suitable areas for species or forest types in terms of climate, topography, and soil characteristics.Our study also covers a larger geographic area than previous studies, and we were able to capture the entire range of species or forest types without being limited by political boundaries.Nevertheless, our main objective was the comparison between forest type range shifts and constituent species range shifts, and the use of interpolated maps enabled an effective comparison at a large geographic scale.
The divergent range shift patterns between forest types and tree species observed here represent only a snapshot of a more prominent trend seen in geological time scales.Forests, because of their sensitivity to changes in tree species compositions, have over the millennia exhibited shorter life spans than individual species (Williams et al., 2004).Positive portfolio effects may provide a critical scientific basis for adaptive forest management and conservation in mitigating the impacts of rapid forest transformation under climate change.