Dominant forest tree species are potentially vulnerable to climate change over large portions of their range even at high latitudes

Climatic data

Normalized (1961–1990) monthly surfaces of total precipitation and average, maximum, and minimum temperatures were downloaded from the USDA Forest Service Rocky Mountain station website (http://forest.moscowfsl.wsu.edu), as were other derived climatic variables (see Rehfeldt et al. (2006) for more details). Data were obtained at a spatial resolution of 0.0083 decimal degrees (≈1 km) and averaged for each 20-km × 20-km grid cell of the modelling area. Climate variables tend to be highly correlated, so we used the VARCLUS procedure in SAS 9.2 (SAS Institute Inc, 2008) to find groups of variables that were as correlated as possible among themselves and as uncorrelated as possible with variables in other clusters. This analysis led to three clusters; in each cluster, we then selected one climate variable influencing plant survival and growth. They were mean annual temperature (TEM), mean annual precipitation (PRE) and useful precipitation (i.e., the ratio of the sum of June, July, August and September monthly precipitation to total annual precipitation; PRATIO).

Ouranos (http://www.ouranos.ca/en/), a consortium on regional climatology and adaptation to climate change, provided different climate simulations using output from 12 general and one regional coupled atmosphere–ocean general circulation models. Each of these was coupled with 1, 2 or 3 projected greenhouse gas emissions scenarios (scenarios B1, A1B and/or A2, based on the Special Report on Emissions Scenarios, or SRES; http://www.ipcc.ch/ipccreports/sres/emission/index.php?idp=0). This generated a total of 70 climate simulations, which are a subset of the 86 climate simulations (Logan et al., 2011) made available from phase 3 of the Coupled Model Intercomparison Project (Meehl et al., 2007).

For each climate simulation, future (2071–2100) TEM, PRE and PRATIO values were obtained using the “change field” method (IPCC, 2001). Monthly mean differences between the baseline period model run (1961–1990) and the future climate model run (2071–2100) were calculated and then combined with baseline values of the observed monthly climate data set. However, due to the relatively coarse spatial resolution of the climate simulations (45 km per cell side for the regional coupled atmosphere–ocean model, and ∼250-km per cell side for the general coupled atmosphere–ocean circulation models), monthly delta values for the centroids of each 20-km × 20-km grid cell (6,418 in all) were interpolated using a linear triangle-based interpolation method (De Berg et al., 2008) between climate model grid cell centroids. Climate simulations for each month were then created by applying interpolated delta values to each observed grid cell value.

To maintain a range of variability in climate projections while reducing computation time, we selected 7 of the 70 available climate simulations as drivers (Table S1), using an objective approach that employs a k-means clustering approach to obtain a good coverage of overall future uncertainty (Casajus et al., 2016). We considered all selected scenarios as equiprobable in this analysis.

Note that emission scenarios are now represented by four Representative Concentration Pathways (RCPs), which became available with the IPCC fifth assessment report and after this study was initiated. The RCPs span a larger range of stabilization, mitigation and non-mitigation pathways than the range covered by the SRES we used (Table S2). As a result, the RCPs now estimate a larger range of temperature increase than the SRES (Rogelj, Meinshausen & Knutti, 2012). Moreover, climate models have also been developed since phase 3 of the Coupled Model Intercomparison Project—CMIP3, particularly by including the representation of biogeochemcal cycles (Flato et al., 2013). Some models do perform better than others for certain climate variables, but no individual model is clearly the best overall. The comparison of median model capability in reproducing historical climate shows relatively modest improvement between CMIP3 and the current CMIP5 generation (Flato et al., 2013). Each generation exhibits a range in performance, with CMIP5 showing fewer ‘bad models’ than CMIP3, but the species distribution modelling uses a consensus approach similar to looking at an ensemble mean of General Circulation Models. As well, for climate projections, we used the delta method (see above) where only relative changes from the GCMs are calculated and then applied to observed data. In this case, a similar range in changes of temperature and precipitation between CMIP5 and CMIP3 is what matters most in terms of species distribution model results. A detailed comparison of CMIP3 and CMIP5 is beyond the scope of this paper, but interested readers can consult Flato et al. (2013), particularly their Figs. 9.44 and Fig. 1 of FAQ 9.1.

Topographic and soil data

Elevation data were provided by the Climate Change Tree Atlas database for the eastern United States portion of the modelling area (Landscape Change Research Group, 2014), whereas for Quebec it was obtained from the Canadian Surface Model Mosaic (http://geogratis.gc.ca/api/en/nrcan-rncan/ess-sst/3A537B2D-7058-FCED-8D0B-76452EC9D01F.html) at a resolution of ca. 20 m and averaged to match our grid. For the eastern United States part of the modelling area, we obtained soil characteristics data (surface deposit and drainage class; Table 1) from the NRCS Soil Survey data (version 2.1, scale 1:24,000; http://websoilsurvey.nrcs.usda.gov/). For the Quebec part, we used soil data from the 3rd decennial permanent and temporary surveys (1:20 000 scale) of the ministère des Forêts, de la Faune et des Parcs, Quebec (available on request at http://www.mffp.gouv.qc.ca/forets/inventaire/donnees-inventaire.jsp). For each grid cell, we computed the percentage of the 20-km × 20-km cell occupied by each level of each edaphic variable.

Species distribution modelling

We computed the geographical distribution of suitable climatic and edaphic conditions—or habitat, as defined by these particular dimensions of the niche—for each of the target tree species, following an ensemble procedure (Araújo & New, 2007) with the BIOMOD 1.1 modelling package (Thuiller et al., 2009) implemented in R (R Development Core Team, 2010). We considered both a baseline period (1961–1990) and a future period (2071–2100, hereafter referred to as 2080).

We used species occurrence data and environmental predictors to build species distribution models using eight modelling techniques: three regression methods (generalized additive models, GAM; generalized linear models, GLM; multivariate adaptive regression splines, MARS), two classification methods (mixture discriminant analysis, MDA; classification tree analysis, CTA) and three machine learning methods (artificial neural networks, ANN; generalized boosted models, GBM; random forest, RF). All models were produced using default BIOMOD parameters where appropriate (Thuiller et al., 2009). Further parameters were as follows: GLMs were generated using quadratic terms and a stepwise procedure with the AIC criteria; GAMs were generated with a spline function with three degrees of smoothing; GBMs were built with a maximum of 2,000 trees; ANNs were produced with five cross-validations (see Marmion et al., 2009a for further details on these modelling techniques). For each species, we built the eight species distribution models using a random subset of data containing 70% of the 20 × 20-km cells (i.e., 4,493 cells). We used the remaining 30% (i.e., 1,925 cells) to evaluate the predictive performance of the models. We repeated this split-sample procedure ten times, thus calibrating 80 different statistical models for each species. We simulated suitability under climate change (future suitability) by projecting each of the 80 projections under each of the seven climate simulations for 2080. This generated a total of 560 probabilities (ten repetitions × eight modelling techniques × seven climate simulations) of habitat suitability for each species for the 2080 period. We combined the different probabilities of habitat suitability (P) based on the area under the receiver-operating characteristic (ROC) curve (AUC) values; we assigned the AUC values from each modelling technique as the weights of the weighted average in order to enhance the contributions of models with higher performance values: (1)${\displaystyle {\mathrm{WAP}}_{i_{\mathrm{baseline}}}=\frac{\sum _{j=1}^8\sum _{k=1}^{10}\left({\mathrm{AUC}}_{jk}\times P_{ijk}\right)}{\sum _{j=1}^8\sum _{k=1}^{10}\left({\mathrm{AUC}}_{jk}\right)}}$ (2)${\displaystyle {\mathrm{WAP}}_{i_{2080}}=\frac{\sum _{j=1}^8\sum _{k=1}^{10}\sum _{l=1}^7\left({\mathrm{AUC}}_{jk}\times P_{ijkl}\right)}{7\times \sum _{j=1}^8\sum _{k=1}^{10}\left({\mathrm{AUC}}_{jk}\right)}}$ where WAP is the weighted average probability of habitat suitability, i is the index of the grid cell (1, …, 6418), j is the modelling technique (GAM, GLM, MARS, CTA, MDA, ANN, GBM, RF), k is the repetition (1, …, 10) and l is the climate simulation (1, …, 7). Averaged projections resulted in a single projection at each grid cell for each species (hereafter referred as the “average model”) for the baseline period (WAP_ibaseline; Eq. 1) and the 2080 period (WAP_i2080; Eq. 2). This method is considered to be more robust than other model fusion methods or single model projections (Marmion et al., 2009b).

Transforming probabilities of suitability to binary values

To transform continuous probabilities of suitability into binary (0/1) values, we calculated a common threshold (cut-off) value for both the baseline period and the 2080 period using a binary vector of observed occurrence and a vector of probability of occurrence from the average model (WAP_i2080). We searched for the threshold which jointly maximized sensitivity and specificity (Liu et al., 2005). This approach is considered among the most reliable for choosing a threshold (Freeman & Moisen, 2008).

Model evaluation

The predictive model performance was evaluated using area under the receiver operating characteristic curve (AUC; Fielding & Bell, 1997) as an accuracy measure. The area under the ROC function (AUC) is usually taken to be an important index because it provides a single measure of overall accuracy that is not dependent upon a particular threshold. Suggested AUC values for classifying the accuracy of models using AUC are: 0.90–1.00 = excellent; 0.80–0.9 = good; 0.70–0.80 = fair; 0.60–0.70 = poor; 0.50–0.60 = fail (e.g., Virkkala et al., 2010 adapted from Swets, 1988). Sensitivity (true positive fraction) and specificity (false positive fraction) values were also reported for each species (Lobo, Jiménez-Valverde & Real, 2008).

Agreement between the average future projection in each cell and the single projections

To measure the level of confidence in our average future projection for a given cell, we also calculated the percentage of the 560 single projections for that cell that agreed with the average projection (hereafter referred as “agreement value”).

Identifying vulnerable habitats under future climates

We focused on Quebec’s productive forest territory to evaluate whether predicted future conditions remained suitable for a species within its baseline range. For this purpose, the baseline range of a species was defined as the set of grid cells within Quebec productive forests where the baseline average model predicted a suitable habitat (WAP_ibaseline ≥ threshold value), as defined by climatic, edaphic and topographic variables. Note that a ‘suitable habitat’ does not necessarily mean an ‘optimal habitat’, since a species can be found on sites with suboptimal conditions. Cells modelled as suitable habitat under baseline climatic conditions, but which became unsuitable under future climate conditions, were classified as unsuitable habitat (UH). Cells modelled as ‘suitable’ under both baseline and future climate further subdivided as:

Less Suitable Habitat (LSH): (3)${\displaystyle \left[{\mathrm{WAP}}_{i_{2080}}-{\mathrm{WAP}}_{i_{\mathrm{baseline}}}<0\&,\left|\right.{\mathrm{WAP}}_{i_{2080}}-{\mathrm{WAP}}_{i_{\mathrm{baseline}}}\left|\right.\ge 0.15\right]}$

Persistent Habitats (PH): (4)${\displaystyle \left[{\mathrm{WAP}}_{i_{2080}}-{\mathrm{WAP}}_{i_{\mathrm{baseline}}}<0\&\left|{\mathrm{WAP}}_{i_{2080}}-{\mathrm{WAP}}_{i_{\mathrm{baseline}}}\right|<0.15\right]}$

OR (5)${\displaystyle \left[{\mathrm{WAP}}_{i_{2080}}-{\mathrm{WAP}}_{i_{\mathrm{baseline}}}\ge 0\right].}$ LSH reflects predicted probabilities of habitat suitability that decrease over time, but not to the point of unsuitability like UH. We used the threshold of a 15% change of probabilities of habitat suitability (WAP_i2080 − WAP_ibaseline) to select the proper subcategory for each cell. This threshold was chosen after examining spatial predictions for 2050 and comparing them with predictions for 2080. The majority of cells classified as less suitable at the 15% threshold in 2050 became unsuitable in 2080. A sensitivity analysis, where the threshold value varied from 5% to 25%, showed how forecasts are affected when this value changes (Fig. S7).

For each species, we reported trends in relation to the entire productive forest territory, the baseline range of the species in Quebec, and each of five vegetation domains.

Black spruce (Table 2; Fig. 2A)

The baseline range for black spruce in the study area essentially covers all five bioclimatic domains. Overall, 78% of the baseline range of black spruce in Quebec’s productive forest is projected to shift towards unsuitable (18%) or less suitable (60%) conditions compared to baseline conditions (agreement value = 68%). Shifts in suitability are projected largely within the sugar maple domain (89% of baseline spruce habitat in that domain shifting to unsuitable), the balsam fir–yellow birch domain (13% shifting to unsuitable), and the balsam fir–white birch domain (2% shifting to unsuitable). Moreover, all the remaining baseline habitats in these domains are projected to become less suitable for black spruce compared to baseline climatic conditions. In the spruce–moss domain, 52% of suitable habitats are projected to become less suitable for the species.

Balsam fir (Table 2; Fig. 2B)

The baseline range for balsam fir covers more than 97% of Quebec’s productive forests. Overall, 59% of the baseline range of balsam fir is projected to shift towards unsuitable (21%) or less suitable (38%) climatic conditions (agreement value = 69%) with climate change. Essentially, all baseline sites over the entire sugar maple domains and the balsam fir–yellow birch domain are projected to become unsuitable or less suitable. Further north, in the balsam fir–white birch domain, shifts towards unsuitability are projected on 1% of the range, while less suitable conditions are projected on another 87%.

White birch (Table 2; Fig. 2C)

White birch is widely distributed in the study area with a baseline range covering 94% of Quebec’s productive forests. Overall, 62% of the baseline range of balsam fir is projected to shift towards unsuitable (14%) or less suitable (48%) climatic conditions (agreement value = 71%) with climate change. In the sugar maple domains, unsuitability is projected on 63% of the baseline range, with the remainder projected as less suitable compared to baseline conditions. Only 2% of habitats shifts towards unsuitability in the balsam fir–yellow birch domain, but less suitable habitats are projected in 67% of the balsam fir–yellow birch domain, 79% of the balsam fir–white birch domain, and 15% of the spruce–moss domain.

Yellow birch (Table 2; Fig. 2D)

The baseline range for yellow birch covers 44% of Quebec’s productive forests. Shifts towards unsuitability (5%) or less suitability (19%) are projected over 24% of the baseline range (agreement value = 78%). All unsuitable areas are in the sugar maple domains (13%), as are most habitats projected as less suitable (48%).

Sugar maple (Table 2; Fig. 2E)

The baseline range of sugar maple covers 31% of Quebec’s productive forests, essentially in the south. Shifts towards unsuitability (8%) or less suitability (1.3%) are projected for 9.3% of the sugar maple baseline range (agreement value = 60%). All sites shifting to unsuitable conditions are in the sugar maple domains. The more northern domains are predicted to maintain their current habitats for sugar maple.