Predicting Climate Change Impacts to the Canadian Boreal Forest

Climate change is expected to alter temperature, precipitation, and seasonality with potentially acute impacts on Canada's boreal. In this research we predicted future spatial distributions of biodiversity in Canada's boreal for 2020, 2050, and 2080 using indirect indicators derived from remote sensing and based on vegetation productivity. Vegetation productivity indices, representing annual amounts and variability of greenness, have been shown to relate to tree and wildlife richness in Canada's boreal. Relationships between historical satellite-derived productivity and climate data were applied to modelled scenarios of future climate to predict and map potential future vegetation productivity for 592 regions across Canada. Results indicated that the pattern of vegetation productivity will become more homogenous, particularly west of Hudson Bay. We expect climate change to impact biodiversity along north/south gradients and by 2080 vegetation distributions will be dominated by processes of seasonality in the north and a combination of cumulative greenness and minimum cover in the south. The Hudson Plains, which host the world's largest and most contiguous wetland, are predicted to experience less seasonality 134 and more greenness. The spatial distribution of predicted trends in vegetation productivity was emphasized over absolute values, in order to support regional biodiversity assessments and conservation planning.


Introduction
Latitudes north of 40° are expected to experience the greatest temperature increases due to climate change [1].In northern locations, impacts of climate change are expected to be acute and the ecological response to climate change will vary spatially [2].Alterations to biodiversity are one of the expected responses to climate change [3].Most commonly, climate change is expected to reduce climatic constraints to plant growth [4].Warmer, wetter conditions will result in increased vegetation productivity, which has been demonstrated as an indirect indicator of biodiversity that correlates with geographic variation in species richness [5].
Canada's boreal is a unique northern environment where climate change is expected to have some of its largest impacts (e.g., [2,3,6,7].To date, changes to fire regimes [8] and periodicity and impacts of insect infestations [9] have been documented and are unevenly distributed due to their dependence on local factors including terrain and moisture conditions.A large and ecologically continuous landscape [10], Canada's boreal has been documented as essential for maintaining biodiversity [11][12][13].Using measures of forest contiguity and functional traits, 80% of Canada's boreal occurs in contiguous blocks [14,15] and provides a range of ecosystem services, such as carbon storage [16] within a stable ecosystem [17]. Over expansive areas such as Canada's boreal, biodiversity management requires large area mapping and monitoring of vegetation conditions.In general, highly productive areas are hypothesized to support more species by providing ample food and habitat resources to partition among species [18,19].Research also indicates that migratory fauna respond to seasonal productive pulses, whereas resident fauna preferentially select areas of consistent annual productivity, demonstrating the importance of modelling integrated annual (and variance) measures of vegetation production [20].
Remote sensing-derived productivity indicators have been shown to provide indirect representations of biodiversity in several areas [5,[20][21][22][23][24][25], including Canada's boreal [26] and have been related to vegetation and wildlife species richness [26][27][28].In addition to providing spatially extensive and continuous data for indirect mapping of biodiversity indicators, remote sensing archives are temporally rich.Historical data from the Advanced Very High Resolution Radiometer (AVHRR) satellite sensor have been used to map Canada's vegetation productivity annually, from 1981-2007, at a 1 km spatial resolution using imagery captured monthly [29].
In our research, productivity mapping is based on the fraction of photosynthetically active radiation (fPAR) absorbed by a vegetated canopy [23,29].Radiation spanning the visible electromagnetic spectrum from 400-700 nm is absorbed by plant pigments for photosynthesis.Integrating maximum monthly absorption can be used as a proxy for production, as growth is indicative of the rate at which vegetation absorbs energy [19].For an area where fPAR remains high over a growing season, high vegetation productivity can be inferred [23].
Several studies have documented a general link between vegetation productivity and biodiversity [5,28,30,31].Though the associations are nuanced and complex, at the most simplistic level increased vegetation productivity generates more energy available to support species and therefore higher biodiversity [32][33][34][35][36][37].More specifically, derivatives of fPAR have been related to species distribution and diversity [28,38].
Given the links between climate, vegetation, and ultimately biodiversity, quantifying possible responses of biodiversity to climate change provides important information to decision makers.Through integration with historical climate data, temporal archives of vegetation productivity provide an opportunity to quantify past trends between climate and productivity, and indirectly biodiversity [39].Modelled relationships between historical climate and productivity data can be integrated with future climate scenarios to map possible future spatial distributions of vegetation productivity.For instance, climate modelling has been used to show that vegetation primary productivity may increase due to climate change [40].
The goal of this research was to understand how climate change might alter vegetation productivity in Canada's boreal forest, and in turn how biodiversity will be impacted.To meet this goal we: (1).Modelled possible spatial distributions of vegetation productivity across Canada's boreal in 2020, 2050, and 2080, using historical climate data, linked to historical remote sensing derived indirect indicators of vegetation productivity from 1987-2007, which are then projected using a range of future climate scenarios; (2).Assessed the performance of predictive models by comparing observed and predicted maps of vegetation productivity using a fuzzy map comparison approach, and (3).Mapped temporal variability in productivity predictions (from 2020, 2050, and 2080) to determine the locations where vegetation appears most and least sensitive to climate change.

Study Area and Data
The boreal region of Canada covers approximately 537 million km 2 and represents one quarter of all naturally functioning forests globally [41] (Figure 1).In contrast to the southern locations of the boreal zone where commercial forestry and human land uses are common, the northern boreal is shaped by active natural disturbances, principally fire [42][43][44].Insect outbreaks are another important disturbance agent in Canada's boreal [45].Typically, insect outbreaks cause spatially localized disturbances.However, climate change is leading to unprecedented insect outbreaks that are having regional effects across North America [46,47].The geographic distribution of productivity within the boreal forest is closely related to climatic regimes with the forests characterized by long, cold, dry winters and short, cool, moist summers.

Ecodistricts and Ecozones
Due to the large spatial extent of the study area, ecozones and ecodistricts were used to provide geographic context and describe landscape characteristics.Ecozones and ecodistricts are two different landscape scales with homogenous and unique interactions between geologic, landform, soil, vegetation, climatic, water and human factors [48].There are 15 ecoregions in Canada's boreal (Table 1).The southern portion of the Boreal Plains ecozone is characterized by trembling aspen (Populous tremuloides) and balsam poplar (Populus balsamifera), as well as open parklands.The central portion of the boreal is covered largely by closed-canopy forest and is comprised of three distinct ecozones: the Boreal Plains, Boreal Shield, and Boreal Cordillera.In this research, as in others (e.g., [16,49]) the Boreal Shield and Taiga Shield are split into east and west units due to differences in precipitation and disturbance regimes.The northern boreal consists of sparse forest and open barrens extending to the northern limit of trees and is the least developed part of the Canadian boreal [15,50].Within Canada's boreal there are 592 ecodistricts with a minimum size of ~100 ha and these formed the basis of all mapping and predictions.

Indirect Biodiversity Indicators
The dynamic habitat index (DHI) is a representation of vegetation productivity and an indirect indicator of biodiversity.The DHI is acquired from three indices, calculated from monthly fPAR composite images, and derived from MODIS, MERIS or AVHRR reflectance values [5].fPAR is a measure of the amount of radiation absorbed by the plant canopy during photosynthesis (400-700 nm wavelengths) [51].A proxy for landscape vegetation primary productivity, fPAR is utilized as a surrogate measure for biodiversity [52].For this research the DHI indices were computed from AVHRR data from 1987-2007 with a spatial resolution of 1 km [29], which is the only boreal-wide satellite dataset available over this time period.The three fPAR indices representing vegetation productivity that were combined into the DHI were: (i) annual cumulative greenness; (ii) seasonal variation in greenness; and (iii) minimum annual cover [5,38].Annual cumulative greenness is an indication of overall vegetation productivity and the potential landscape productivity at any given location [19].Annual cumulative greenness was estimated by summing monthly fPAR observations for each year [23].Seasonal variation in greenness is directly related to local climate and geography and characterizes annual variations in vegetation productivity at a given location.Calculated as the coefficient of variation, seasonal variation in greenness is sensitive to extreme changes in vegetation cover such that the index is low in highly productive evergreen forests and high in areas exhibiting seasonal snow covers in the winter and vegetative growth in the spring.Minimum annual cover indicates the lowest amount of vegetated cover at a specific location within a year.Minimum annual cover provides information about the landscape's capacity to support adequate levels of green vegetation cover and permanent resident species throughout the year [23].Baseline data on annual cumulative greenness, seasonal variation in greenness, and minimum annual cover are provided in Figure 2. We set the baseline to 1990, as future predictions are provided at 30 year intervals for 2020, 2050, and 2080.The DHI in 1990 represents conditions 30 years prior to the first prediction and provides data at consistent time intervals.

Historical Climate Data
Historical climate data were provided by the Pacific Climate Impacts Consortium (PCIC) and were derived from models produced by the National Centers for Environmental Prediction (NCEP) North American Regional Reanalysis (NARR) project.The NARR project provides a long-term, dynamically consistent, high-resolution (32 km), high-frequency, atmospheric dataset for North America [53] that is available as point data.In order to spatially integrate the NARR data with the DHI, monthly climate datasets from 1987-2007 were interpolated using natural neighbours and stored in 1 km by 1 km grid cells.Daily climate variables were summarized annually and include precipitation (prec), maximum temperature (tmax), minimum temperature (tmin), and the mean annual growing season index (GSI) [54,55], which is an estimation of the mean vapour pressure deficit (iVDP) and mean suboptimal temperatures (iTmin).The vapour pressure deficit of the atmosphere serves as an index of vegetation evaporative demand where water stress constrains vegetation growth due to stomatal closure and reduced photosynthesis [6,54,56].Suboptimal temperatures affect phenological development of vegetation during the first two to three months of the growing season, whereby cooler conditions retard and warmer conditions accelerate development [57].The growing season index is computed by multiplying the suboptimal temperature, vapour pressure deficit, and photoperiod [54] and serves as indicator of the relative constraints to phenological development of vegetation due to climatic limits.

Scenarios of Future Climate
Three Intergovernmental Panel on Climate Change (IPCC) future climate scenarios were used: B1, A1B, and A2.The B1 scenario represents the least extreme climate change, the A1B (scenario represents business as usual, and the A2 represents the most extreme climate change.Based on the IPPC scenarios, maps of future climate for Canada's boreal were obtained from the Canadian Centre for Climate Modelling and Analysis (CCCma) and represent the thirty-year climate averages centered on 2020, 2050, and 2080.By using three models we aimed to represent climate change conditions over a full range of possible scenarios [58].Mote et al. [59] compared the ability of 10 global climate models to track recently recorded trends in temperature and precipitation in the Pacific northwest and reported that the Canadian model predicted trends consistent with other models, including a relatively rapid rate in warming and increases in precipitation when compared to actual observations.

Predicting Future Spatial Distributions of Forest Productivity
Described in more detail below, the general approach to prediction was to first quantify relationships between historical climate and the DHI.Using a Random Forest data mining technique that was used for each of the DHI components (annual cumulative greenness, seasonal variation in greenness, and minimum annual cover) we generated three 1987-2006 climate-productivity models.Using model fits, built from historical data, and future climate scenarios (A1B, A2, B1) each component of the DHI index was forecasted for the years 2020, 2050, and 2080, resulting in 27 future vegetation productivity maps at the spatial resolution of ecodistricts.In addition, the DHI indices were visualized in a composite image by simultaneously displaying seasonal variation in greenness, annual cumulative greenness, and minimum annual cover, which enabled interactions between the three DHI components to be explored.
For each of the three DHI components (annual cumulative greenness, seasonal variation in greenness, and minimum annual cover) a Random Forest model was generated to quantify relationships with historical climate data.Multiple fully grown regression trees ("the forest") were used to create a model fit.These singular trees are formed by hierarchical binary partitions of singular predictor variables (ecoregion and climate indices) to explain the variation in the response variable (cumulative, seasonal and minimum fPAR levels) [60].Predictor partitions were selected using iterative ANOVA analysis to ensure the greatest decrease in the sum of the squares within the vegetation productivity groups (see [61]).The variable selected at each split was chosen from a random subset of the predictors.Predictions of vegetation productivity at the terminal splits were averaged over all trees in the forest to predict the final annual productivity values [60].In our case, 999 singular regression trees, created from bootstrapped samples, produced stable predictions.
The Random Forest algorithm was utilized for prediction as it does not make assumptions about the shapes of relationship and has been demonstrated to be a flexible model for ecological prediction [60,61].Using an ensemble of trees is known to reduce error rates by smoothing the severity of the prediction boundaries created by singular regression trees [61].In addition, the implementation of a second randomization step, of selecting a subset of predictors to be chosen from for the binary splits, is known to reduce the effects of redundant variables and allowed indices that may have been masked by dominant predictors to be selected [61,62].The second randomization step also introduced a greater display of interaction effects between climate variables that may otherwise have been hidden [63].
Random Forest model fits for annual cumulative greenness, seasonal variation in greenness, and minimum annual cover, were assessed by the mean square error and percentage of variance explained, which is a generalized error rate calculated from out of the bag data [61], or data not used in the bootstrapped sample.To interpret which climate variable had the greatest influence on the productivity indices' predictions, node purity was also used to indicate the average total reduced sum of the squares attained by all partitions on the environmental variable of interest [61].

Model Validation
For model validation the Random Forest models were applied to data from 1987-2006 to predict the three DHI indices in 2007.The 2007 observed and predicted maps of annual cumulative greenness, seasonal variation in greenness, and minimum annual cover were then compared using a technique known as fuzzy numerical map comparison [64].The fuzzy numerical map comparison considers the similarity in modelled and observed values by quantifying how similar a DHI value is for a pixel, within the pixel's spatial neighbourhood, and the mean of the map or study area.By comparison with the mean, the fuzzy numerical statistic identifies where the local pattern of similarity is unique [65]  Where s i (A,B) is the one-way similarity between map A and B at cell i.Index j iterates through all N cells and the function w(d) defines the neighborhood, which defined by distance or, in this case, adjacent pixels.S i (A,B) combines the two one-way similarities into an overall similarity.S(A,B) is the overall map similarity, it is the mean of overall n locations.The function f(a,b) determines the similarity of two values.The fuzzy numerical statistic values range from zero to one, where zero indicates no similarity and one indicates complete similarity.

Predicting Future Spatial Distributions of Forest Productivity
For the annual cumulative greenness and seasonal variation in greenness, the models explained 79.3% and 79.0% of the variance respectively.Large increases in node purity and increases in mean square errors (after variable value permutations) indicated climate indices with the greatest importance to the overall prediction of annual vegetation productivity characteristics (Table 2).The variables with the strongest impact on node purity, a measure of the variable importance for cumulative vegetation production, were maximum temperature (tmax) and mean suboptimal temperatures (iTmin).Precipitation (prec) and mean vapour pressure deficit (iVDP) were the variables that when included reduced mean square error the most.For the minimum cover, the model explained 72.0% of variance and, as with the other models, maximum temperature (tmax) and mean suboptimal temperatures (iTmin) increased node purity.For minimum cover precipitation (prec) and mean annual growing season index (GSI) had the greatest impact on mean square error.
Table 2. Random Forest model fits.Percent increases in MSE (mean squared error) are indicative of the degree of the model fit based on the data used to train the Random Forest.Node purity is a summarized measure of how much the prediction was improved when the variable was selected (i.e., sum of the decrease in root mean square error across all trees in the forest).As models of climate change became more extreme, or as climate change continued through time, increased annual cumulative greenness shifts northward through Canada's boreal and intensifies.The most extreme climate model scenario (A2) predicted the highest expansion of cumulative greenness into areas of the Hudson Plains wetlands and Taiga Plains tundra by 2080 (Figure 3).More conservative climate change scenarios, the A1B and B1 models, also predicted increased greenness, but the southern boreal forest was predicted to be more strongly impacted.

DHI Indice Annual Cumulative Greenness
Figure 4 illustrates that as the climate changes, seasonal variation in greenness decreases.Future scenarios of seasonality indicated that the region around the Great Lakes is predicted, by the moderate (A1B) and high climate change (A2) scenarios, to be less seasonally coupled with an increase in cumulative productivity by 2080 (Figure 4).In contrast, the conservative climate change model (B1) indicated a reduction in the lowest seasonality class in the south boreal regions.
In general, minimum annual cover will decrease under increasing climate change.Forecasts showed a strong decrease in the southern portion of the Boreal Plains (Figure 5).However, throughout the Boreal Shield an increase in minimum cover was predicted, with the most extreme model (A2) predicted to have high minimum cover in the Hudson Plains' region by 2080.
Simultaneous display of the three DHI indices for 1990 and the A1B scenario for 2080 allowed for more complex interpretation of habitat shifts due to the interaction of vegetation productivity processes (Figure 6).Results predicted reduced heterogeneity in vegetation productivity, particularly west of Hudson Bay as the complex colors shown in 1990 become more homogenous aqua and red in 2080.Reduced vegetation heterogeneity indicated that in the south, vegetation productivity is the expression of a combination of annual cumulative greenness and minimum annual cover, while in the north seasonal variation in greenness is the dominant process.

Model Validation
Comparison of the 2007 observed and predicted DHI generated fuzzy numerical values for the cumulative greenness ranging from 0.7-1.0 with a mean of 0.9 suggesting strong agreement.Fuzzy numerical values for the coefficient of variation or seasonality were similar, ranging from 0.6-1.0 with a mean of 0.9.The minimum cover had less agreement between observed and predicted with fuzzy numerical values ranging from 0.1-1.0 and a mean of 0.6.
The level of model agreement varied spatially (Figure 7).For annual cumulative greenness, the lowest agreement between modelled and observed 2007 data occurred primarily in the Southern Arctic and Taiga Plains and agreement was highest in the Hudson Plains.The Atlantic Maritime Ecozone, which is small in area, also showed poor agreement for annual cumulative greenness.For seasonal variation in greenness the lowest agreement occurred in the Hudson Plains, while the Taiga ecozones had the highest agreement.Minimum annual cover, overall, had much lower agreement, and the least agreement occurred in the Hudson Plains and Taiga Shield.All DHI indices had poor agreement in the Boreal Cordillera, likely due to the small area and thus the relationships that govern productivity may not be well represented in the model.

Discussion
Developing climate/vegetation relationships from 1987 to 2006 and then predicting 2007 vegetation spatial distributions (as DHI indices) allowed assessment of model performance beyond variance explained.Both model variance explained and model performance were quite high, 72.0%-74.0%and a mean of 0.6-0.9respectively.Though, minimum annual cover model performance was highly variable and as low as 0.1.High variability in predictive success of minimum annual cover may be related to low variance, as minimum annual cover could be low for several different environments [39].For example, minimum annual cover will be low in deciduous forests as well as grasslands.Given the large spatial extent of our study area it was essential that the interaction between climate processes and vegetation vary geographically.By building a spatially explicit model that included ecozones, we did not assume spatial non-stationarity.Rather, we allowed modelled relationships between climate and vegetation to be developed based on regional dynamics [61].
Using geographically explicit map comparison, we were also able to quantify the geographic nature of uncertainty and better understand ecological factors influencing model confidence [66].Annual cumulative greenness was typically low in areas with long winters and cool summers (Southern Arctic and Taiga Plains), possibly reflecting the importance of temperature in the model (Figure 3).The Atlantic Maritime Ecozone, which also showed poor agreement, may be impacted by oceanic influence that mediates climate [3] and reduced associations between temperature and vegetation productivity [67,68].The lower agreement between modelled and predicted results for seasonal variation in greenness in the Hudson Plains may reflect the presence of the largest contiguous wetland in the world [26,69].Wetlands are a complex mosaic of vegetation (trees, shrubs, and herbs) interspersed with water bodies of various sizes, posing a difficult environment for consistent capture with coarse spatial resolution remote sensing [26,70].Much of the vegetation in the Hudson Plains is low and easily covered with temporally dynamic snow and ice [5].In contrast to forest dominated environments, wetlands express cover and seasonality characteristics that are not as well captured using the remote sensing methods, with further investigation warranted.Minimum annual cover produced much lower levels of agreement between observed and modelled maps and was more difficult to model due to numerous near zero values, which can occur in environments ranging from rock and ice to locations exhibiting vegetation senescence cycles.The lowest agreement occurred in the Hudson Plains, which again could be the result of the unique response of wetlands [5,26,70,71], and in the Taiga Shield, where a patchwork of wetlands, forests, meadows, and shrub occurs in a cold climate [69].The strength of mapping minimum annual cover is best realized when it is utilized in conjunction with annual cumulative greenness.
Maps of modelled change in individual DHI indices indicated that overall cumulative greenness will increase, seasonality will decrease, and minimum cover will increase.Similar trends have been predicted for other regions of Canada [39].The amount of change will vary spatially and increase through time with more extreme climate scenarios (Figures 3-5).A north/south gradient of change was apparent and matches the observed and experienced temperature increases north of 40° latitude [1].These mid to high latitudes have had diurnal decreases in temperature ranges and decadal increases in autumn and winter precipitation contributing already to a 10% decrease in snow cover and ice extent exacerbated in the most northern extents [72].Using the composite image we can identify, using the legend, how the combination of DHI components are expressed spatially.The composite images of DHI from 1990 and 2080 (Figure 6) revealed that, when compared to 1990, vegetation productivity modelled to 2080 (using the business-as-usual (A1B) climate change scenario) will be most impacted by the increased seasonality predicted in the north and a combination of cumulative greenness and minimum cover in the south.More inter-and intra-annual variability in precipitation and snowfall is expected [73].Areas in the west and east were modelled as changing the most.For example, Newfoundland is currently influenced by multiple productivity components and predicted to be strongly impacted by a combination of seasonality and minimum cover with a reduction in cumulative greenness.In the northwest, the Boreal Cordillera, Taiga Cordillera, and western portion of the Taiga Plains will also undergo substantive change influenced by seasonality more than any other DHI indices.
The relatively large predicted change in vegetation in the Hudson Plains and associated wetlands (Figures 3-6) is in line with other research that has indicated climate change is a major threat to the integrity of wetlands [74].Relative to other wetlands worldwide, the high latitude of Canada's wetlands make them particularly vulnerable to climate change [75].Previous research has suggested that drier and warmer conditions may lead to a 70%-90% loss of wetlands in the Hudson Plains [76].Reduced duration of snow cover, associated with predicted increases in minimum greenness and decreased seasonality, will negatively impact some duck species in boreal wetlands [77].For example, reduced flexibility in timing of breeding has been predicted to impact population levels of species such as scaups (Aythya spp.) and scoters (Melanitta spp.).
The Taiga has similar climate characteristics but different geology, flora, and fauna [71] to the Hudson Plains and was predicted to undergo substantial change in vegetation (west for annual cumulative greenness; east for seasonal variation in greenness and minimum annual cover).The Taiga is characterized by discontinuous permafrost and annual temperatures that are below freezing.Observing rapid change the Taiga may indicate that climate change will surpass a threshold that will cause a substantive response in vegetation [78].For instance, in the western Taiga, Olthof and Pouliot [78] predicted conifer recruitment and growth would occur in the north.
A striking result of the composite images was the loss of heterogeneity in the DHI.Increased homogeneity may lead to reductions in species richness with fewer niches for species to select [79].Examples of species of particular concern is the yellow rail (Coturnicops noveboracensis), an at-risk wetland bird whose narrow breeding habitat and important boreal breeding ranges south of Hudson and James Bay lowlands makes it especially vulnerable.Yellow rail breeding habitat is tightly constrained to shallow water levels (i.e., <15 cm for initial nesting) and the availability of overlying dry dead mats of vegetation, environments akin to wet marshy areas dominated by grass-like vegetation [80].As such, climate-induced habitat loss or degradation through wetland reduction or shifts in wetland configuration and altered hydrological conditions, such as flooding or drought, represent an important threat to this species [80].
Research has previously documented that wide-ranging and sensitive species, such as birds and butterflies, follow vegetation productivity trends modelled by remotely sensed imagery [20][21][22][23]27,28,81]. Notably, seasonal trends captured by the fPAR index have been able to distinguish spatial variation in bird species richness, including species' responses to seasonal and minimum annual productivity levels [23] and associations between composite DHI heterogeneity [81].Demonstrated associations between vegetation productivity, represented by fPAR and the DHI index, suggest that predicting future possible spatial distributions of vegetation productivity provides guidance on how species spatial distributions may change.The strength of our approach is that it is general, covers a large area, and is repeatable.We provide a landscape perspective on predicted trends in the response of vegetation to climate change.The response of vegetation to climate is anticipated to correlate with change in species richness and distributions.However, responses will be complex and impacted by the interaction between multiple species [2] and nonlinear responses of species to climate change [82].Our results are most appropriately interpreted as broad scale and relatively trends on spatial variation in the response of vegetation productivity and species distribution to climate change [9].
While our projections showed that climate change may increase the overall potential for greenness, other factors were not accounted for.For instance, soils will not change at the same rate as climate and may become a stronger limiting factor for vegetation growth.For example, our models showed potential increased greenness in higher elevations but thin soils may limit vegetation development.As well our models do not capture episodic events and climate extremes, such as droughts and flooding [72] or forest fire [43], which are all expected to increase in either magnitude or frequency due to climate change [83].

Conclusions
We have integrated historical climate and vegetation productivity data with future IPCC scenarios of climate change to map future spatial distributions of indirect indicators of biodiversity for Canada's boreal.Our results indicate that vegetation productivity, an indirect indicator of biodiversity, will become more homogenous, particularly west of Hudson Bay.We can also expect climate to impact biodiversity along north/south gradients and by 2080 vegetation distributions will be dominated by processes of seasonal variation in greenness in the north and a combination of annual cumulative greenness and minimum annual cover in the south.The Hudson Plains are modelled as having the potential to experience substantive change becoming more productive or green with less seasonality.
In this research, we highlighted opportunities for space-time modelling of indirect indicators of biodiversity.Remote sensing provides datasets with continuous spatial coverage and long temporal extents and these are beneficial for mapping and modelling how vegetation has and will change.Though all models suffer uncertainty, managers tasked with planning and conserving biodiversity, particularly in large areas like Canada's boreal, benefit from information on possible trends and changes in vegetation.The maps generated thorough this research are not the final word.One of the strengths of our methodology is that models and maps can be updated as new and better data come online.However, representations of future trends in vegetation productivity provide important context for long-term decision making and possible adaptation opportunities in Canada's boreal.

Figure 2 .
Figure 2. The three fraction of photosynthetically active radiation (fPAR) indices that represent vegetation productivity and are combined to produce the dynamic habitat index (DHI): annual cumulative greenness; seasonal variation in greenness; and minimum annual cover.Annual cumulative greenness is an indication of overall vegetation productivity and the potential productivity at any given location.Seasonal variation in greenness is calculated as the coefficient of variation and is sensitive to extreme changes in vegetation cover.Minimum annual cover indicates the lowest amount of vegetated cover at a specific location within a year.Maps shown represent conditions in 1990.

Figure 6 .
Figure 6.Composite DHI model results for observed 1990 and 2080 (A1B).Maps are a simultaneous display of the three DHI components where the red, green, and blue color bands represent seasonality, cumulative greenness, and minimum cover, respectively.Dark colored regions have low values in all DHI components while light colored areas have high values.Areas of definitive color indicate a dominant habitat indicator; for example a region that is purely red denotes high seasonal variation in greenness.As indicators shift in 2080, such as when seasonal variation in greenness decreases, other components such as annual cumulative greenness or minimum annual cover fill into newly favorable habitats.Common DHI blending occurred with light blue, yellow, and darker colorations like purple or dark orange.The light blue regions represent the most productive habitats with high minimum annual cover, high annual cumulative greenness, and low seasonal variation in greenness.Darker purple and orange colorations indicate low annual cumulative greenness, low seasonal variation in greenness, and low minimum annual cover.

Figure 7 .
Figure 7. Map comparison results comparing modelled and observed DHI components in 2007.Values are an average of map comparison values for each ecozone, where high values indicate more similarity between observed and modelled maps.