Combining US and Canadian forest inventories to assess habitat suitability and migration potential of 25 tree species under climate change

To evaluate current and future dynamics of 25 tree species spanning United States and Canada.


| INTRODUC TI ON
Climate change has already impacted many natural systems (Scheffers et al., 2016) and is projected to have increasingly serious impacts as the century progresses. Forest ecosystems are of particular concern due to the high levels of biodiversity that they maintain (Aerts & Honnay, 2011) and the slow migration rates of trees (Corlett & Westcott, 2013;Feurdean et al., 2013). Consequently, understanding the response of forests to climate change has been an ongoing focus of scientific research for several decades (see Iverson & Prasad, 1998;Keenan, 2015). A key aspect of this response is the extent to which long-distance tree migration can track the climatic shifts projected for the coming decades (Aitken, Yeaman, Holliday, Wang, & Curtis-McLane, 2008;Iverson, Schwartz, & Prasad, 2004;Prasad, Gardiner, Iverson, Matthews, & Peters, 2013).
Modelling tree migration under climate change is most effective at broad spatial extents which incorporate entire species ranges.
One reason for this is that the climate is not changing evenly across the planet. For example, the boreal forest in North America has already experienced more than a 1.5°C increase in mean annual temperature since 1900 (Price et al., 2013), compared to an increase of 0.7°C across the conterminous United States (USGCRP, 2017). These changes, in combination with projections for modest, regionalized changes in precipitation (Gauthier, Bernier, Kuuluvainen, Shvidenko, & Schepaschenko, 2015) indicate that many species could experience different climatic drivers across their range (Bouchard, Aquilué, Périé, & Lambert, 2019;Housset, Girardin, Baconnet, Carcaillet, & Bergeron, 2015;Matthews, Sadler, Kubota, Woodall, & Pugh, 2019).
Efforts to model range shifts of tree species under climate change have evolved through time Iverson & Prasad, 1998;Iverson, Prasad, Matthews, & Peters, 2008). Early studies tended to focus solely on climate envelope shifts, identifying suitable climate habitat under current and future conditions (McKenney, Pedlar, Lawrence, Campbell, & Hutchinson, 2007a, 2007bMcKenney, Pedlar, Rood, & Price, 2011;Shafer, Bartlein, & Thompson, 2001;Sykes, Prentice, & Cramer, 1996). Other studies have attempted to generate better predictions by incorporating key drivers of tree suitability, including edaphic variables (Chambers, Périé, Casajus, & de Blois, 2013;Iverson & Prasad, 1998;Iverson et al., 2008;, migration constraints Prasad et al., 2013) and competitive interactions (Araujo & Luoto, 2007). Regardless of the exact approach taken, most efforts incorporate habitat models that relate species abundance or occurrence to a suite of environmental variables. While many studies have employed occurrence data for this step (e.g. Bell & Schlaepfer, 2016;Case & Lawler, 2017), the use of abundance data allows for a more detailed and useful assessment of habitat suitability (Howard, Stephens, Pearce-Higgins, Gregory, & Willis, 2014;Iverson et al., 2008). However, in North America, seamless coverage of tree abundance data has been challenging to obtain due to national-level differences in data collection methodology and availability in the United States and Canada (see Beaudoin et al., 2014;Forest Inventory Analysis, 2017). Consequently, attempts to model tree migration using abundance-based habitat models have typically been regional in scope (Chambers et al., 2013;Coops & Waring, 2011;Gray & Hamann, 2013), with minimal attempts at cross-border analyses using abundance data owing to several challenges related to the merging of different datasets.
In order to overcome this limitation, we have combined United States' Forest Inventory Analysis (FIA) data (Forest Inventory Analysis, 2017) with spatial datasets derived from Canada's National Forest Inventory (NFI) data (Beaudoin et al., 2014;Beaudoin, Bernier, Villemaire, Guindon, & Guo, 2017) to produce a spatial database of relative abundance for 25 ecologically and commercially important tree species that span the two countries (Table 1). We use this database to model climatic habitat quality (HQ, also referred to as habitat suitability or climatic habitat suitability) under current, and future climate (RCP 4.5 and 8.5 emission scenarios- Moss et al., 2008) using a multimodel ensemble approach (Prasad, 2018). We further employ a migration simulation model (Prasad et al., 2013) to compute colonization likelihoods (CL, also referred to as migration potential) associated with future suitable habitats. Finally, we combine HQ and CL to visualize and assess potential tree shifts under climate change. Specifically, our goals are to a) examine the potential spatial employ range-wide data, evaluate colonization potentials and enhance cross-border collaborations.

K E Y W O R D S
colonization likelihood, Forest Inventory Analysis, habitat suitability, long-distance migration, multimodel ensemble, National Forest Inventory, random forest, range-wide analysis, stochastic gradient boosting distribution of future suitable habitats across the entire range of each species, b) assess the likelihood that newly gained suitable habitats will be colonized in the future, c) evaluate habitats colonized, gained, lost and retained across the two countries and d) provide a template to evaluate current and future management options for important tree species. We use a multistage modelling approach (Prasad, Iverson, Matthews, & Peters, 2016) to minimize the bias and confoundment of the model by modelling HQ as a purely climate response Price, Cooke, Metsaranta, & Kurz, 2015) and allowing for other dependencies, such as land cover and elevation, to be evaluated as optional post-model filters, where appropriate.

| ME THODS
Our overall methodology and the data used to model habitat suitability and colonization potentials are illustrated in Figure 1.

| Combining FIA and NFI-kNN
We use FIA data in the eastern United States to derive relative importance values (RIV, also known as relative abundance) of tree species (Iverson et al., 2008;Peters et al., 2019). This index is calculated as the sum of the (1) abundance (i.e. number of stems) and (2) dominance (i.e. basal area) of each species relative to all trees within the plot (~83K plots in the eastern United States-for details on the FIA plot design please refer to the FIADB user guide available online, and O' Connell et al., 2017). Because the FIA plot locations are unequally sampled across the United States, we aggregated relative importance values over grid cells (20 x 20 km) TA B L E 1 Combined R-square values from the random forest and gradient boosting models for the 25 tree species in this study  (Gillis, Boudewyn, Power, & Russo, 2010). The remote sensing component consists of 2 x 2 km sample units located on a systematic national sampling grid of 20km, with a total of 13,158 sample units in the first current 10-year (2008-2017) re-measurement survey.
Recently, Beaudoin et al. (2014), Beaudoin et al. (2017) have produced gridded estimates of forest attributes from the NFI, using a k-nearest neighbour (kNN) approach to impute cell values from a combination of NFI photo plot data (Gillis, Omule, & Brierley, 2005), and 26 geospatial data layers including MODIS spectral imagery, climate and topographic variables. One of the products derived from this project was a set of species composition maps that gave estimates of species' relative abundance.
We merged grids of relative abundance from the two data sources (FIA and NFI-kNN) to obtain continuous tree abundance estimates for each tree species. Both indices were scaled from 0 to 100 to indicate the relative per cent composition of each species. One consequence of this merger was that the coarser estimates from the kNN-based approach tended to produce smoother patterns of species relative abundance patterns, likely influenced by the secondary processing. Although admittedly not seamless, the merger provided continuous estimates of relative abundance up to the northern edge of each species range, which was vital for our regression-based ensemble modelling approach. We aggregated the abundance data for the 25 species to a 20 x 20 km resolution that smoothed differences in relative abundance methodologies across the USA-Canada border and allowed a computationally feasible framework for looking at macroscale patterns across the entire range.
Of the 25 species, there were only two (quaking aspen and paper birch) whose ranges extended significantly into the western United States and thus were not fully represented by our abundance data, which was comprised of FIA plots east of the 100th meridian.
Additionally, there were some data gaps in the NFI-kNN data for quaking aspen in Ontario, Canada.

| Climate data
Predictors for the HQ model (described below) were climate normals from the AdaptWest database (AdaptWest Project, 2015;Wang, Hamann, Spittlehouse, & Murdock, 2012) for current  and future (2071-2100) periods. We obtained data for future climate scenarios (Moss et al., 2008), RCP 4.5 (medium emission pathway) and RCP 8.5 (high emission pathway), in order to examine the degree to which habitat projections differed under these alternative pathways (Wang, Hamann, Spittlehouse, & Carroll, 2016). Future climate projections were based on an ensemble of 15 general circulation models (AdaptWest Project, 2015). The original data, at 1 km resolution, was aggregated to 20 km to match the abundance variables for the HQ model.

| Multimodel ensemble (MME)
Habitat quality (HQ) for current and future climate scenarios was predicted using a multimodel ensemble (MME) technique (Prasad, 2018) at a cell resolution of 20 x 20 km. The HQ model was developed using current climate values at tree abundance cells and then projected into the future under the RCP 4.5 and 8.5 emissions scenarios. We screened 27 bioclimatic climate variables available in the AdaptWest database using random forest and stochastic gradient boosting (see MME below) to find the best fit for all the 25 species. We selected a subset of 21 climate variables (  (Cutler et al., 2007;Prasad, Iverson, & Liaw, 2006). The version of RF with random number of splits, known as extremely randomized forest (ERF), takes this one step further by creating p random predictor splits in each node (in contrast with RF, which chooses the 'best' split), independent of the response variable and then the split with the best gain (mean-squared error for regression) is chosen. The rationale is that by randomizing the split the variance is further reduced, but with a slight increase in model bias. To compensate, ERF typically uses the entire learning sample instead of a bootstrapped one to keep the bias low (Geurts, Ernst, & Wehenkel, 2006).
Boosting is a method of iteratively converting weak learners (regressors that are weakly correlated with the true model) to strong ones using decision trees as the base learner by reweighting observations that have higher errors. Stochastic gradient boosting (gbm module in R) reduces variance by shrinkage (regularization) and stochasticity in which each newer iteration of the model learns from the previous one by minimizing residuals based on the shrinkage parameter (Friedman, 2002). Xgboost is a slightly different approach to boosting which implements an enhanced regularization technique to limit overfitting (Chen & Guestrin, 2016).
The MME approach is best used where prediction uncertainty needs to be stabilized to yield more robust predictions that reflects the consensus (Jones & Cheung, 2015;Martre et al., 2015) by averaging individual model errors (Zhang, Huang, et al., 2015;Zhang, Liu, et al., 2015). It further minimizes bias and variance (Hastie, Tibshirani, & Friedman, 2009) as well as prediction error while extrapolating current prediction to future climates. Only predictions where there was 'consensus' among all models (i.e. where all models predicted nonzero values) were averaged to ensure that the main trends were captured. For the multimodel approach to work well, the models should be based on a similar framework (in our case, decision trees) but should adopt structurally different approaches (the different approaches to random forest and boosting) so that the final ensemble averages these heterogeneous approaches (Tebaldi & Knutti, 2007).
The MME also helps mitigate the effects of spurious model artefacts at the low end of the abundance spectrum. Here, the 'consensus' approach averages prediction signals that are present in all four individual models of the ensemble, thus highlighting the dominant patterns of prediction in the ensemble (Prasad, 2018).

| MME validation
The four models used in the ensemble have built-in validation methods that ensure that the models do not overfit and provide estimates of their predictive accuracy. For example, in Random Forest, the training is performed using a bootstrap sample (with replacement), which included 2/3 of the data. At each iteration of the bootstrap sample, the remaining data (called out-of-bag, OOB) are used to measure the predictive accuracy of regression trees grown from the bootstrap sample; these errors are aggregated to obtain the OOB error rate (Svetnik et al., 2003). With large datasets like ours, where the number of observations far exceeds the number of predictors, and with large numbers of decision trees (1,000 in our case), the OOB error rate matches very closely with the cross-validation error rate from separate training and test datasets (Hastie et al., 2009).
For boosting, cross-validation was performed (we used 10-fold) to calculate an estimate of the generalization error using 1,000 trees.
There was close correspondence between the error rates of Random Forest and stochastic gradient boosting, and we therefore calculated the average R-square to estimate overall model performance (Table 1).
To ensure that we chose the optimal parameters for the random forest and boosting variants in our ensemble, we used a package called caret in R which streamlines the model training process for complex regression and classification problems by doing a stratified random split of the data into training and test sets. These two sets are used to evaluate the effect of tuning parameters on model

TA B L E 2
The 21 climate variables used as predictors in the HQ (MME) model performance and choose the optimal model across these parameters (Kuhn, 2008).
Finally, our approach to average only the consensus of the four models in the ensemble ensured we chose the most dominant trend across these models, further boosting our confidence in the final result.

| Estimation of colonization potentials
Colonization likelihoods (CL) were computed for each species using a long-distance migration model that incorporates current abundance, historical migration rates and current habitat fragmentation (Schwartz, 1993). Simulation of long-distance migration was implemented via a fat-tailed inverse power function at a cell resolution of 1 km. The likelihood of an unoccupied cell becoming colonized by a particular species during each generation is a function of that species' abundance in the surrounding cells, the habitat quality of the unoccupied cell (must have at least 10% per cent forest cover, to ensure that it is a habitable cell) and a search window distance function. The colonization likelihood for each unoccupied cell is calculated by summing over all occupied cells at each generation. The stochastic nature of the tail of the long-distance migration is simulated by drawing a random number from an even distribution and comparing it with the calculated likelihood to determine if the cell gets colonized. Simulation using multiple historical migration rates and application to multiple species was made feasible by using convolution and Fast Fourier Transforms to keep the computation time to a minimum (Prasad et al., 2013). For this study, we use an optimistic, but historically defensible migration rate of 50km/century for all species and a generous search window of 500 km (to accommodate stochastic long-distance dispersals) as the limit for migration within each generation. These estimates are generally at the higher end of reported tree migration rates (Davis, 2001;McLachlan & Clark, 2004;McLachlan, Clark, & Manos, 2005; but see Snell & Cowling, 2015;Ordonez & Williams, 2013). We deliberately chose not to parameterize individual species migration rates owing to large uncertainties in life histories and dispersal syndromes among the 25 tree species (Iverson et al., 2004).
A recent improvement to this model (cf. Prasad et al., 2013)

| Combining suitable habitats and colonization likelihoods
The habitat quality (HQ) based on relative abundance is estimated by the MME model for the RCP 4.5 and RCP 8.5 scenarios. Outputs are scaled 0 to 100, where 0 represents absence and 100 represents the maximum where the cell contains monotypic stands of that species.
These outputs were reclassified into three classes: low (1-5), medium (6-15) and high (16-100). The current distribution of relative abundance based on FIA and NFI (scaled 0-100) was also classified into three classes-low, medium and high-based on the HQ categories described above.

| Temporal matching of HQ and CL
As can be expected, combining suitable habitats and colonization likelihoods through time poses a problem because of the approximate temporal periods of the various data and model outputs. The current distribution of the species was based on data that spanned from approximately 2008-2017 from the FIA and NFI-kNN data. The species were allowed to migrate for ~100 years and CL was calculated relative to this current distribution. Because future climatic habitats (HQ) for RCP 4.5 and 8.5 were based on the 2071-2100 equilibrium scenarios, and the current climate spanned 1981-2010, future temporal matching of HQ and CL is approximate. We have therefore assumed that allowing the species to migrate ~100 years with a slight temporal overshoot is justified for combining HQ and CL rather than a premature termination of colonization ( Figure 1). Further, we found that this approximation does not alter the outcome of this work to any appreciable degree. Additionally, we illustrate the potential HQ lost, gained and retained and report the per cent colonization of the newly suitable habitats. We also assess the overall trends of HQ and CL for the 25 species. The maps for all 25 species are available in the Data Availability section.

| HQ Model performance
The combined R-square for the HQ models, based on an average of the random forest and stochastic gradient boosting models, shows that model reliability varied widely among the 25 species (Table 1)

| Habitats lost, gained and retained
Sugar maple (Acer saccharum) is currently distributed primarily in the eastern United States (Figure 2), but also extends into southern Ontario and Quebec as well as the coastal regions of eastern Canada. Little's range extent (Little, 1971, used extensively for tree species ranges in North America; red line in Figure

| Quantifying Habitat Quality and Colonization Likelihood
The combination of HQ & CL (Figure 4a

| Habitat dynamics, all 25 species
In addition to the species-specific maps and graphs, assessing all 25 species together provides insights into general patterns in habitat and colonization dynamics. This analysis also allows comparison of the HQ gained and lost for the two climate scenarios (RCP 4.5 and scenario compared to the RCP 4.5 scenario. This is consistent with the greater warming, and consequently larger northward shifts, associated with RCP 8.5 scenario. However, the degree to which species exhibited differential responses to the two emissions scenarios appear tied to current range location, with boreal species (e.g. black spruce) that are already close to the northern boundaries of the  5 and 6) is often an order of magnitude larger compared to habitat lost (Figures 7 and 8).

F I G U R E 9
The percentage of newly gained habitat that was projected to be colonized by the end of the current century for 25 tree species under RCP 4.5 and RCP 8.5 emissions scenarios were projected to have some probability of colonization within ~100 years under future climates (e.g. coloured cells in Figure 4), as compared to cells that had essentially no chance of colonization (grey cells in Figure 4). For each species and emissions scenario, we calculated the percentage of new habitats (depicted as shades of green-GainLow, GainMedium and GainHigh, in Figure 3) that have the potential to be colonized (Figure 9). In general, the proportion of habitats colonized is much higher under RCP 4.5 than for RCP 8.5, reflecting the greater amount-and more northerly

| Changes in suitable habitats and management implications
Changes in climatically suitable habitats provide useful metrics to understand how habitat expansion under future climates can impact the United States differently compared to Canada. The future habitat projections also provide a visual map of the altered climatic landscapes in which these species will exist in the future.
Many species that are currently in the eastern United States are projected to gain suitable habitat in Canada. The maps of habitats gained, lost and retained ( Figure 3, and Data Availability section), and the multispecies summaries for each country separately ( Figures 6 to 9), show how these dynamics play out for individual species, which is useful for forest managers in both countries.
Species whose current range is primarily in the United States (red maple-Acer rubrum, silver maple-Acer saccharinum, sugar maple-Acer saccharum, white ash-Fraxinus americana, bigtooth aspen-Populus grandidentata, northern red oak-Quercus rubra) and especially those that just reach the southern edge of Canada  (Figures 6 to 9).
All 25 species gained more habitat than they lost under future warming scenarios, often by an order of magnitude ( Figures   6, 7 versus Figures 8, 9). The species that are projected to have the highest gains in the United States include silver maple (Acer saccharinum), black walnut (Juglans nigra), bur oak (Quercus macrocarpa) and American elm (Ulmus americana)-these species have relatively low R-square values (Table 1) In the United States, almost all species lost more habitats under the harsher RCP 8.5 scenario compared to RCP 4.5 (Figure 7). In Canada, this trend is reversed, with the notable exception of black spruce (Picea mariana; Figure 8). This pattern is consistent with the notion that habitats will shift further northward under more extreme warming scenarios. For many species in our study, this translates into habitat gains in Canada, though northern species that already extend to nearly the northern land limits of the continent (e.g. black spruce -Picea mariana) would be expected to gain little habitat from such shifts.
The loss of the high HQ class is greater relative to low and medium HQ (Figure 7) for the United States, which is cause for concern. The extent of this trend varies among species, but should be factored into management protocols. More focused management in areas where the models are depicting this trend may be warranted.
In Canada, although the amount of habitat lost is lower compared to the United States, black spruce, balsam poplar and quaking aspen all lose sizeable amounts of habitat. In fact, quaking aspen (Populus tremuloides) loses a significant amount of high-quality habitat in both the United States and Canada. While these losses in habitat are cause for concern, species like black spruce (Picea mariana) are also quite opportunistic and can colonize relatively rapidly under changing conditions (Joyce & Rehfeldt, 2017;Truchon-Savard et al., 2019).
The magnitude of the habitat losses reported here is relatively small compared to those reported in some earlier studies. For example, McKenney, Pedlar, Lawrence, Campbell, and Hutchinson (2007b) reported that over half of the tree species in their study were projected to have smaller climatic niches in the future (compared to none in the current study)-with the majority of losses happening at species' southern range limits. One reason for these different findings may be the use of abundance (as opposed to occurrence) data in the current study, which allows for a more detailed (and nuanced) delineation of suitable habitat. In fact, for many of the species examined here, southern range limits are projected to remain relatively intact, albeit with reduced habitat quality. Thus, while leading-edge habitats are projected to expand rapidly due to greater warming at these latitudes, trailing-edge populations may remain viable, but subject to greater climatic stress. Findings such as this illustrate the value of the range-wide abundance data employed here. This finding reflects the rapid rate of warming projected for the north-western region of North America (Bintanja, 2018;Coumou, Di Capua, Vavrus, Wang, & Wang, 2018;Goosse et al., 2018). This gain in HQ via Alaska signals future potential possibilities for the United States, but, as indicated by our colonization models, many species would have essentially no chance of reaching this region by means of natural migration by the end of the century, as discussed in the next section.

| Colonization of suitable habitats-HQCL combinations
While our HQ models provide valuable insights into future habitat shifts, it is well understood that only a portion of future suitable habitats will be naturally colonized by the end of the century.
Therefore, the combination of HQ with CL allows for a more realistic assessment of the future habitat occupancy by incorporating colonization likelihoods of future suitable habitats based on an optimistic historical migration rate of ~50 km/century (Figure 4, online maps).
As might be expected, the area of HQ colonized is much smaller than that uncolonized for most species (grey shades dominate over coloured ones in Figure 4, especially under RCP 8.5). However, for species with scattered, spotty distributions (e.g. bitternut hickory- These findings (Figure 9, Figure 4 and rest of the maps in Data Availability section) can help to identify suitable conservation strategies for the various tree species examined here. For example, for the species listed above with considerable potential for infilling migration, efforts could focus on facilitating within-range movements, while for those species relying primarily on migration at range edges, efforts aimed at assisted range expansion may be appropriate. Our results suggest that differing approaches could be featured in the two countries-for example, a focus on within-range conservation in the United States and range expansion in Canada. These contrasting situations require different approaches, but a common paradigm of international cooperation to combat the consequences of rapid climate change.

| Role of assisted migration
As outlined above, this study identifies a number of situations in which various forms of assisted migration could be considered as potential tools for promoting forest conservation and adaptation under climate change (Pedlar et al., 2012 Regardless of the exact approach taken, efforts to conserve forest habitats across United States (including Alaska) and Canada will call for close cooperation between the two countries, not only for mitigating ill effects, but also for realizing the opportunities these changes offer.

| Post-model filter
Even though a sizeable portion of the gained climatic habitats are colonized according to our models (Figure 9), these are estimates under a full climatic response and optimistic historical migration rate According to our HQ model, under the RCP 8.5 scenario (Figure 3,   right), approximately 4.31 million km 2 of HQ-low, 2.5 million km 2 of HQ-medium and 0.38 million km 2 of HQ-high become available by the end of the current century. These amounts, with the land cover mask overlaid, shrink by 35% for HQ-low; 20% for HQ-medium; and 9% for HQ-high. Additionally, in line with our multistage scheme, regions where topography, soil depth or type (Carteron et al., 2020) are major modifiers can be further filtered or modelled by more accurate regional databases, depending on the interest/expertise of the manager, to produce locally relevant maps.

| Limitations and strengths of the FIA + NFI approach
One of the novel aspects of this study is the combination of relative importance data (derived from the USA's FIA dataset) with per cent composition data (derived from Canada's NFI-based data product), which allowed us to model the entire range of almost all the species considered here. Though there are differences between the two data sources, and the merger of FIA and NFI data while not optimal, can be justified in terms of 'information gain' derived from our MME approach for predicting HQ. Additionally, because FIA data were not available for interior Alaska that area could not contribute to the MME; however, these data would have included only a few species, such as black spruce and quaking aspen, and not changed our overall conclusions.
One of the consequences of post-processing of the NFI data (via MODIS and k-nearest neighbour approach) is that the relative abundance is smoothed (because it tends to eliminate the tails of the distribution) compared to the FIA plot-based data. This may have resulted in more conservative migration potentials, but our optimistic assumption of the migration rate of 50 km/century should compensate for this deficiency (Davis, 2001;McLachlan & Clark, 2004;McLachlan et al., 2005; but see Snell & Cowling, 2015;Ordonez & Williams, 2013).

| Future directions
An important consideration with respect to tree responses to climate change is that populations may be adapted to climate in only a portion of the range (i.e. local adaptation). This means that, for widespread species, populations may respond differently to climate across the species' range (Leites, Rehfeldt, & Steiner, 2019;Peterson, Doak, & Morris, 2019). In addition, natural populations tend to occur in climates colder than where they grow the best-that is, populations occur in climates suboptimal for their growth and productivity-this discrepancy increases with the degree of winter cold of the provenance (Rehfeldt, Leites, Joyce, & Weiskittel, 2018). However, capturing this variability across the entire species' range via habitat suitability models is not possible without data from provenance studies and/or genetic analysis, which tend to be sparse or non-existent for many species. One way to approximately assess this variability in the absence of these data is to split the range into three regions based on plant hardiness zones (PHZ): leading cold-adapted, a trailing warm-adapted and the 'optimal' middle (Prasad, 2015). We are in the process of using this approach for the 25 species studied here to evaluate intraspecific variability in conjunction with the range-wide analysis to understand how the leading-edge and trailing-edge responses can be different. Our current multistage approach, when combined with the analysis of intraspecific variability, will have direct relevance to future climatic changes because the leading and trailing regions can be under greater climatic stress, although they may possess genotypes better adapted to extreme climates (Hampe & Petit, 2005;Isaac-Renton et al., 2018;Rehm, Olivas, Stroud, & Feeley, 2015).
In conclusion, our study underscores the importance of evaluating range-wide habitat suitability and colonization potential of tree species when examining climate change impacts on future tree distributions. For the species examined here, habitat losses were primarily experienced along southern range limits, while habitat gains were associated with northern range limits-though there was also a significant amount of infilling migration within the ranges of species with discontinuous distributions. However, even with optimistic tree migration rates, our models predict that, for most of the species examined here, only a small portion of the climatic habitat generated by climate change will be colonized naturally by the end of the current century. Significant management efforts, including assisted migration of populations and species, may need to be considered in order to conserve and promote tree species that are deemed important for economic or ecological reasons. Our work highlights the need for increased cross-border cooperation as management demands can differ between the United States and Canada owing to the differing species-specific objectives and responses under future climates.

ACK N OWLED G EM ENTS
We are grateful to the field crews of the Forest Inventory Analysis program in the United States and the National Forest Inventory in Canada-this paper would not have been possible without their tireless efforts. Support in obtaining and preparing the kNN-NFI data was provided by Kevin Lawrence at Great Lakes Forestry Centre.
We would like to thank Stephen Handler and Paul Berrang for their incisive reviews and to the anonymous reviewers whose comments greatly improved the paper.

CO N FLI C T O F I NTE R E S T S
We have no competing interests.

DATA AVA I L A B I L I T Y S TAT E M E N T
Can be accessed via: https://doi.org/10.5061/dryad.qz612 jmbn (36 MB). 1) Unzip the downloaded file to a folder, 2) go to the 'Data' folder and 3) open file index.html. This will allow rapid visualization of three sets of maps for the 25 species: 1) abundance distribution; 2) habitats lost, gained and retained (RCP 4.5 and RCP 8.5); 3) combined maps of HQ and CL (RCP 4.5 and RCP 8.5).