Acquiring knowledge about habitat selection for threatened species is crucial to their recovery, allowing for essential habitats to be identified and protected. This may be particularly acute for species reintroduced to landscapes where they have long been absent, because the response of founding individuals to potentially novel environments may be uncertain (Frair et al. 2010, Osborne and Seddon 2012). Moreover, knowledge of habitat preferences is critical for anticipating the expansion of recovering populations (Jung 2017), and can guide the selection of suitable sites for future reintroductions (Le Gouar et al. 2012).

After approaching the brink of extinction, American bison Bison bison have been reintroduced to parts of their historic range; yet, they currently occupy a small fraction of the landscape they once did (COSEWIC 2013). Bison restoration efforts in the North American boreal forest are focused on the survival and growth of 11 small, geographically isolated populations that occur across landscapes that vary in habitat composition. As bulk feeders (Reynolds et al. 1978, Jung 2015), the preservation of foraging areas, such as graminoid-dominated habitats (e.g. wet sedge meadows and xeric grasslands) may be critical to bison recovery. Similar to European bison B. bonasus (Kuemmerle et al. 2018), American bison exhibit behavioral plasticity that has allowed them to persist in diverse landscapes, likely resulting in disparate habitat selection strategies among distinct populations. They readily adapt their diets based on forage availability (Hecker et al. 2021), enabling reintroduced populations of bison to persist in landscapes that may differ in vegetation composition from historical habitats (Rivals et al. 2007, Kowalczyk et al. 2019). In the boreal forest, reintroduced bison are established in a variety of dissimilar ecozones that vary with respect to elevation, terrain ruggedness, dense forest cover and graminoid-dominated habitats. The adaptability of bison to local conditions, combined with their relatively recent reintroduction to a heterogeneous boreal landscape, has resulted in a limited understanding of their habitat selection behaviour.

Research on several bison populations has revealed both similarities and differences in diet and habitat use. Bison select seasonal habitats based on food availability (Larter and Gates 1991, Fortin et al. 2003, Jung et al. 2018a), but their choices may be constrained by the energetic costs of locomotion (Fortin et al. 2003, Belanger et al. 2020), abundance of biting insects in summer (Melton et al. 1989, Belanger et al. 2020) and predation risk (Fortin et al. 2009, Harvey and Fortin 2013). Both forage quantity and quality are strong determinants of bison habitat use (Hudson and Frank 1987, Larter and Gates 1991), as with other grazing or browsing ungulates (Hebblewhite et al. 2008, VanBeest et al. 2010). During winter, bison may prefer graminoid-dominated wetlands (Larter and Gates 1991, Jung et al. 2018a, Belanger et al. 2020) where they can obtain the sedges and rushes that dominate their winter diet (Fortin et al. 2003, Larter and Allaire 2007, Strong and Gates 2009, Jung 2015); however, some populations also use woody browse (Campbell and Hinkes 1983, Waggoner and Hinkes 1986, Jung et al. 2015). In summer, bison tend to avoid wetlands due to soft footing and an abundance of biting flies (Belanger et al. 2020), instead adopting a broader diet found in alternative habitats such as shrub meadows and eskers (Larter and Gates 1991, Larter and Allaire 2007, Belanger et al. 2020), or alpine tundra in mountainous regions (Jung et al. 2015, 2018b). Bison may spend little time in coniferous and deciduous forests (Jung et al. 2018a) where there is minimal forage (Strong and Gates 2009, Belanger et al. 2020), but forests may be used periodically for resting, ruminating and avoiding biting flies and predators (Jensen et al. 2004, Fortin et al. 2009), or when lichen is consumed during times of food scarcity (Larter and Gates 1991). Predation by wolves Canis lupus – the primary predators of bison – may be an important limitation for long-established (> 50 years) bison populations (Larter et al. 1994, Carbyn et al. 1998), but is rarely observed in recently reintroduced populations (Fortin et al. 2003, Jung 2011) for which the role of predation in habitat selection is poorly understood.

Figure 1.

Reintroduced bison Bison bison from the Nahanni population foraging near the end of an airstrip in northwestern Canada. Photo by Cameron D. Eckert.

Some bison populations occupy wilderness regions, whereas others persist in industrialized landscapes, and many have highways bisecting their range. For some ungulates, anthropogenic disturbance causes detrimental habitat loss or increased predation rates (Neilson and Boutin 2017, Plante et al. 2018), but bison may benefit from some degree of disturbance. Linear features such as seismic lines and roads may reduce the energetic costs of travel (Bruggeman et al. 2006, Jung 2017, DeMars et al. 2020) and support the growth of key forage plants such as graminoids and shrubs (Strong et al. 2013, Finnegan et al. 2018). Resource and residential developments also create grazing habitats that may attract bison, such as clearcuts (Redburn et al. 2008), airstrips (Fig. 1), fields and lawns (Larter and Allaire 2007). However, a limiting factor to bison recovery is the negative attitudes of local people towards bison as a result of their use of roads and communities, where they can become safety hazards or damage property (Doney et al. 2018, Bath et al. 2021). For example, lack of social acceptance has resulted in resistance to the establishment and growth of some reintroduced bison populations (Clark et al. 2016, Klich et al. 2018, Jung 2020). In addition, roads may become ecological traps for recovering bison due to mortality from vehicle collisions (COSEWIC 2013). Consequently, there is a need to better understand the selection of roads and other anthropogenic features so that wildlife managers can work to find socially acceptable solutions that will improve attitudes towards bison.

Our objective was to evaluate seasonal habitat selection by a transboundary population of reintroduced bison in northwestern Canada, known as the Nahanni population, in order to inform conservation planning efforts. This population has a strong association with riverine habitat along a major river (i.e. Liard River; Larter and Allaire 2007, Thomas et al. 2021). The population's range has a low to moderate density of linear features associated with resource exploration and extraction, and bison are known to frequent local communities and a highway that bisects their range (Fig. 2), causing conflict with local residents and hazards to people and bison. Thus, we also assessed bison selection of anthropogenic features to better understand the potential for human–bison conflict. We predicted that bison would use different habitat in winter and summer, selecting forage-rich graminoid-dominated wetlands in winter but occupying more diverse habitats (e.g. riverine habitats such as gravel bars and islands) in summer, when they broaden their diet and avoid wetlands with poor footing and high densities of biting insects. We also predicted that bison would select roads year-round, ostensibly for forage, insect relief (summer) and ease of travel (winter), and would be particularly attracted to the abundant forage found near communities (e.g. lawns and airstrips) during summer. We predicted that bison would select anthropogenic clearings associated with forestry and natural gas development, as these may provide more forage than the surrounding forest. Because habitat selection may differ between individuals in the same population, and these individual differences may be concealed by analyses that pool all individuals (Gillingham and Parker 2008, Nielsen et al. 2009, Anderson and Johnson 2014), we accounted for variation in behaviour among individuals in our analyses. We predicted that individual-based resource selection would differ from population-level resource selection derived by pooling animals.

Figure 2.

Land cover and linear disturbances within the study area in Northwest Territories, Canada, which comprises the home ranges (100% MCPs) of all GPS-collared bison (Bison bison; n = 10) in this study.

Methods

Study area

The Nahanni bison population is located in northwestern Canada, adjacent to Nahanni National Park and spanning the borders of Northwest Territories, Yukon and British Columbia (Fig. 2). Bison were reintroduced to the region in 1980, 1989 and 1998 as part of a national recovery effort (Larter and Allaire 2007, COSEWIC 2013). From approximately 160 animals in 1998, the population grew to an estimated 962 animals in 2017 (Larter 2021). The population is ‘wild by nature’ (Jung 2020); that is, free-ranging (unfenced) and subject to ecological and evolutionary processes such as competition and predation. The Nahanni population is geographically isolated from other bison populations.

The range of the Nahanni population is located in the taiga plains ecozone, characterized by low-lying forested plains with the Mackenzie Mountains and foothills along the western boundary. Low elevation areas are dominated by dense boreal forest (Fig. 2), consisting of spruce (Picea glauca, P. mariana), pine Pinus banksiana, tamarack Larix laricina, aspen Populus tremuloides, poplar P. balsamifera and birch Betula papyrifera. The forest is punctuated by patchily-distributed wetlands, including bogs and graminoid-dominated fens, the latter of which are often associated with lakes. The population is strongly associated with the Liard River, which contains numerous islands and gravel bars that are used by foraging bison (Larter and Allaire 2007, Thomas et al. 2021). The range is subject to resource development and contains a network of seismic lines and secondary service roads (Fig. 2). Otherwise, human presence on the landscape is minimal; there are two small villages and one unpaved highway. Climate is continental, with short, warm summers and long, cold winters. Relative to other reintroduced populations in northern Canada (Larter and Gates 1991, Redburn et al. 2008, Jung et al. 2018b, Belanger et al. 2020), the range of the Nahanni population is densely forested, lacking large patches of wet sedge meadow, regenerating forest or relict boreal grasslands.

Table 1.

Summary of location data collected and predictive accuracy of pooled RSF models and individual RSF models for GPS-collared bison Bison bison (n=10) from the Nahanni population in northwestern Canada, during summer and winter. Fix success rates were calculated as the number of successful fixes by the total number of possible fixes. Predictive accuracy of RSF models is based on k-fold cross-validation procedures. r_s=Spearman's rank correlation coefficient of habitat bin rank by area-adjusted frequency of used locations per habitat bin.

Habitat data

We selected habitat variables anticipated to have an influence on bison habitat selection (Table 2). These variables may relate to the abundance of forage, ease of travel and the presence of biting insects (summer only). Habitat data were primarily derived from the NASA Arctic Boreal Vulnerability Experiment (ABoVE) land cover product documenting vegetation cover in 2014 at 30-m resolution. We aggregated the existing dataset into eight habitat classes with relevance to bison: coniferous forest, deciduous forest, fen, bog, shrub, herbaceous, fluvial and anthropogenic (Fig. 2; defined in Table 2). Our fluvial habitat class, including gravel bars and islands, was modified from the existing ‘barren’ land cover class within and along major rivers. Anthropogenic disturbances (e.g. logging, mining, airstrips) were not adequately delineated by existing maps, so we manually digitized disturbance polygons using Google Earth imagery. We reasoned that bison select habitat at larger scales than our pixel size (30 m), and may occupy habitat adjacent to a feature of interest (i.e. foraging patch) while resting and ruminating. Thus, we calculated the percent cover of each habitat class within 100 m of a given 30-m pixel using a moving window analysis in ArcGIS (ver. 10.8, ESRI, Redlands, CA). We chose 100 m buffers because in a previous study bison responded more to land cover at 100 m than at 1000 m scales (Jung et al. 2018b).

Table 2.

Descriptions of predictor variables used to evaluate habitat selection in GPS-collared bison Bison bison (n = 10), Northwest Territories, Canada.

Linear feature data, including roads and seismic lines, were provided by the Government of Northwest Territories and included disturbances up to 2013. Given our study area was largely wilderness (0.1% anthropogenic disturbance), we assumed disturbances (linear and non-linear) were relatively static throughout our period of study, and historic satellite images showed little change. We used the Line Density tool in ArcGIS to calculate linear feature densities (km km^–2). All distance variables (distance to primary roads, secondary roads, communities, waterbodies and the Liard River) were transformed using an exponential decay function because we expected landscape features to have a stronger influence on bison at closer distances. The function took the form e^–αd, where d was the distance to a feature and α was a constant (Nielsen et al. 2009). We set the value of α to 0.02 so the effect of landscape features dropped to < 15% beyond 1 km, and declined to 0% at a distance of 5 km. We derived elevation data from a digital elevation model at 30-m resolution, and calculated a terrain ruggedness index (TRI) based on the methods proposed in Riley et al. (1999) using the spatialEco package in R.

Habitat selection models and validation

We developed resource selection function models to assess habitat selection by bison in summer and winter based on a used/available design (Manly et al. 2002). Used habitats were defined by GPS collar locations. To derive available locations, we established annual home ranges for each bison using 100% minimum convex polygons (MCPs; Johnson and Gillingham 2008). We considered this an appropriate method for defining available habitat because MCPs did not contain large areas that were unused by bison or inaccessible to them (Northrup et al. 2013). We generated random points within each home range at a ratio of five available locations for each used location (Milakovic et al. 2012, Morrison et al. 2014), resulting in densities of 2–11 available points km^–2. Thus, we assessed resource selection within individual home ranges (third-order selection, DeCesare et al. 2012) and used observations of individuals to make inferences about the population (design III from Manly et al. 2002). We considered beta coefficients (i.e. selection coefficients) to be representative of relative habitat selection.

We constructed six candidate models based on a priori predictions about factors that would influence bison habitat selection in summer and winter, including a global model and a null model (Table 3). We also evaluated a ‘combined’ model with a selection of predictors representing land cover, water features, anthropogenic disturbance and topography, because multiple factors may influence bison habitat selection. To facilitate comparison of selection coefficients, all continuous predictor variables were standardized to have a mean of zero and a standard deviation of one. To prevent multicollinearity, we calculated variance inflation factor (VIF) scores and dropped any variables with VIF > 3 (Zuur et al. 2010). We compared models using Akaike's information criterion (AIC) scores. We considered any model with an AIC weight ≥ 0.9 to be the single best model, with very little model-selection uncertainty, and based any subsequent inference on that model alone (Burnham and Anderson 2002, Arnold 2010).

In resource selection studies, data from multiple animals are often pooled together without adequate consideration for inter-animal variation and how this affects population-level estimates (Gillies et al. 2006, Gillingham and Parker 2008). Thus, we took two approaches to RSF modeling. First, we conducted an RSF on pooled data from all collared animals within each season, using mixed-effects logistic regression with animal ID as a random intercept to account for correlated location data from the same animal (Gillies et al. 2006). We used the lme4 package in R to run mixed-effects models. Because female bison are gregarious, we tested for independence among collared bison by calculating the minimum and mean distance between individuals at regular time intervals during which collars were simultaneously active. We considered GPS-collared bison to have formed a group when their locations were < 100 m apart (Fortin et al. 2009). Second, we calculated separate RSFs for each individual animal during each season via logistic regression. This approach was important because we anticipated substantial inter-animal variation as a result of unequal sample sizes and inconsistent land cover across home ranges. Such variation can result in unreliable estimates of selection coefficients when data are pooled, and the precision of such estimates may be overestimated (Gillies et al. 2006, Gillingham and Parker 2008). We then estimated population-level habitat selection from individual RSFs via a two-stage approach (Nielsen et al. 2009, Fieberg et al. 2010). After estimating selection coefficients for each animal, we averaged individual coefficients (taken from the top model for each animal) using an inverse variance weighting (IVW) method (Gillies et al. 2006, Sawyer et al. 2006, Nielsen et al. 2009, Takahata et al. 2014). The IVW method accounts for variation in the precision of individual coefficient estimates by affording higher weight to those animals with larger sample sizes and lower standard errors. This approach is analogous to a pooled model with random slopes for each predictor (Gillies et al. 2006, Nielsen et al. 2009). We calculated the variance of IVW selection coefficients using a conservative method that accounted for variation among animals as well as uncertainty surrounding individual estimates (Marzluff et al. 2004, Sawyer et al. 2006, Anderson and Johnson 2014).

Table 3.

Model structure and selection results for mixed-effects resource selection function (RSF) models with pooled data from GPS-collared bison Bison bison (n = 10) from the Nahanni population in northwestern Canada, during summer and winter. Predictor variables are described in Table 2.

We evaluated the same set of models for both the pooled and individual RSFs, in both seasons. However, some land cover variables (e.g. bogs) were dropped from individual RSF models because those habitats were completely or nearly completely avoided by the animal (e.g. < 5 used points) during a given season, resulting in perfect separation between used and available habitat such that selection coefficients could not be estimated (Gillingham and Parker 2008). For distance variables, perfect separation occurred when there were no used locations within 5 km of a feature (i.e. the maximum decay distance). Consequently, complete habitat avoidance was not captured in individual RSFs, and we noted all cases where this occurred. Lastly, we compared selection coefficients from the pooled RSF to averaged IVW coefficients from individual RSFs. We calculated 95% confidence intervals for selection coefficients, and considered those with confidence intervals that did not overlap zero to have a significant effect on habitat selection.

We used the averaged selection coefficients to predict and map habitat selection probabilities across the study area. We defined our study area as a polygon encompassing all of the individual minimum convex polygon home ranges of GPS-collared bison. We applied the following equation to raster habitat layers with standardized values:

where β_k is the estimated selection coefficient for predictor variable x_k. Predictions were classified into 10 quantile bins that represented relative habitat selection (Boyce et al. 2002, Morris et al. 2016). We evaluated predictive accuracy via k-fold cross-validation (Boyce et al. 2002, k=5) of the top model for each individual bison and the pooled data. GPS locations were randomly assigned to one of five folds for cross validation. The validation procedure was performed five times, each time withholding a different test set and validation results were averaged. Our metric of predictive accuracy was the Spearman's rank correlation (r_s) between the habitat bin rank and the area-adjusted frequency of used locations per habitat bin.

Results

We captured and collared 12 bison from 2007 to 2017. Two collars either malfunctioned or fell off within 30 days of deployment; data from the remaining 10 animals (9♀, 1♂) were used for habitat selection analysis. Collars were active intermittently from 2007 to 2019 with variable lifespans (x̄ = 481 ± 311 [SD] days, range = 159–964 days; Table 1). Fix success rates were also variable and all were < 100% (x̄ = 80.8% ± 14.5%, range = 55.1–95.5%; Table 1), likely as a result of bison damaging collars (Jung and Kuba 2015, Jung et al. 2018b). We collected 17 260 used locations (after data screening) and 86 300 available locations across both seasons, with 9169 used locations in summer (x̄ = 917 ± 814 per bison, range = 77–2575) and 8089 in winter (x̄ = 809 ± 710, range = 214–2274). Collared bison typically travelled independently of one another. The mean distance between pairs of bison with simultaneously active collars, at any time, ranged from 5.2 to 79.2 km (overall x̄ = 30.4 ± 29.2 km). Five bison pairs (involving six individuals) travelled together intermittently for short periods of time, demonstrating the unstable fission–fusion dynamics that are typical of bison populations (Fortin et al. 2009, Merkle et al. 2015, Jung 2020). For four of these pairs, < 5.8% of locations were within 100 m of one another during the same time period. For one pair, 31.6% of locations were within 100 m. Because dyads were not common in our data, we retained all location data for further analyses.

The top RSF model for pooled data was the global model for both summer and winter (Table 3). This model included all predictor variables excluding coniferous forest, which was dropped as a result of a high VIF score. Bog was also omitted from the pooled model in winter due to near-perfect avoidance (i.e. almost all used locations were in areas where bog cover was 0%) and consequent convergence issues. The global model received an AIC weight of 1.00, indicating this was our best model, and the only one with adequate support. The top model for each individual bison in both summer and winter was also the global model (excluding coniferous forest), which received an AIC weight of 1.00 in all cases; thus, we only report selection coefficients from global models. Although all habitat variables except coniferous forest and bog were included in the pooled models, additional variables were dropped from several individual RSF models because the animal completely (or nearly completely) avoided that habitat type. In summer, 5 of 10 collared bison completely avoided bogs, two avoided communities, one avoided the Liard River and one avoided anthropogenic habitat. In winter, all bison showed complete or near-complete avoidance of bogs, and seven bison strongly avoided communities. Both primary roads and anthropogenic habitat were avoided by five bison in winter, fluvial habitat and the Liard River were avoided by three, and secondary roads were avoided by two. We did not calculate averaged IVW coefficients unless a predictor variable was included in a minimum of five individual RSF models. Consequently, we did not calculate IVW coefficients for bog and distance to communities during winter, but we calculated coefficients for all predictor variables in summer.

During summer, averaged IVW selection coefficients indicated that bison selected shrubby and herbaceous land cover, fens, fluvial habitat and anthropogenic disturbances (Fig. 3, Table 4). Bison also tended to avoid bogs in summer, but there were substantial inconsistencies among individuals. Bison had an affinity for areas near primary and secondary roads and used areas close to the Liard River, and selected lower elevation areas. In winter, bison again selected shrubby habitats and showed stronger selection of fens than in summer. Bison used anthropogenic disturbances, but to a lesser extent than in summer. Bison showed relatively strong selection for areas near waterbodies, and preferred secondary roads.

The direction of selection (i.e. the sign of selection coefficients) estimated by the pooled RSF and the IVW method generally agreed, with the exception of distance to primary roads, distance to Liard and elevation for winter models (Table 4, Fig. 3). At a population level during winter, bison appeared to avoid primary roads, avoid the Liard River and select lower elevations, while the IVW method suggested the opposite behaviours. Because the IVW method adequately accounted for variation in individual selection, IVW coefficients were typically smaller than pooled RSF coefficients when inter-animal variation was high (e.g. bogs in summer, elevation in winter; Table 4). Further, the estimated precision of selection coefficients was much lower with the IVW method, resulting in fewer predictors that were statistically significant (Table 4, Fig. 3). This is likely a more accurate reflection of the uncertainty associated with population-level coefficients.

Results from individual RSFs suggest there was variation in selection behaviour that was not always captured by the pooled RSF, which likely overestimated precision. However, there was consistent selection or avoidance of some habitat types (Table 4). For example, in summer all 10 bison selected shrub, and nine bison selected herbaceous habitat, anthropogenic disturbances and primary roads. During winter, 8 of 10 animals strongly selected fens and nine selected areas near waterbodies. In addition, there were two cases where complete avoidance of a habitat type resulted in discrepancy between individual RSFs and the pooled and IVW results (Table 4). Notably, seven bison completely avoided communities in winter, but the pooled RSF, which obtained estimates from the remaining three individuals, indicated significant selection for communities (i.e. a negative distance coefficient), suggesting that those animals that selected communities likely stayed there for longer periods. Five bison avoided anthropogenic habitat in winter, but the remaining five bison selected it, resulting in significant positive coefficients for both the pooled RSF and IVW results.

IVW averaged coefficients derived from individual RSFs were used to predict and map the relative probability of bison habitat selection in summer and winter (Fig. 4). Individual RSF models generally had good predictive capacity, although summer models ( = 0.87, range = 0.63–0.97) were typically more accurate than winter models ( = 0.79, range = 0.38–0.97; Table 1). Of the individual models, only the winter model for the lone male in our study had poor predictive ability (r_s = 0.38); however, this bison's collar was active for only ∼2 months during the winter (Table 1). The pooled RSF models had moderate predictive success, but values were lower than the mean accuracy of individual RSFs (summer: r_s = 0.66, winter: r_s = 0.78).

Discussion

A key finding from our study was the variation in habitat selection behaviour among individual bison, particularly in winter. Some bison apparently occupied low elevation areas along the highway corridor and visited communities, while others dispersed away from the highway to higher elevations, where they sought fen complexes and secondary industrial roads. Broader dispersal in winter may have alleviated inter-specific competition for limited forage (Plumb et al. 2009). We would have overlooked this variation in habitat selection strategies, had we relied solely on a pooled RSF and not investigated individual behaviour. The variation that we uncovered in individual behaviour may indicate subpopulation structure, similar to that found in caribou Rangifer tarandus (Nagy et al. 2011), Yellowstone bison (Olexa and Gogan 2007, Halbert et al. 2012) and elk Cervus canadensis (Walter et al. 2010). However, we did not have data adequate to detect subpopulation structure, which may be defined by spatial distribution, social interaction, diet and habitat selection. In the future, detection of subpopulation structure could be aided by analysis of a larger number of collared bison with simultaneously active collars (Nagy et al. 2011), as well as dietary and genetic data (Walter et al. 2010, Halbert et al. 2012). Regardless, similar to other studies (Gillingham and Parker 2008, Nielsen et al. 2009, Anderson and Johnson 2014), our analyses broadly demonstrate how habitat models estimated by pooling data from multiple animals may conceal individual behaviour, which may have implications for management decisions. We advocate for further use of individual-based approaches when investigating habitat selection behaviour, particularly when few animals provide location data or subpopulation structure is apparent.

With regard to bison habitat selection, our main finding was that Nahanni bison selected habitat in a manner that was generally consistent with what is known about bison ecology in the boreal forest, yet we also discovered new insights. As predicted, Nahanni bison showed divergent habitat selection behaviour in winter and summer, reflecting seasonal changes in forage availability and other constraints. During winter they selected wetlands (fens) and waterbodies associated with them, despite their relative scarcity on the landscape (4.5% of the study area) and their patchy distribution within a matrix of forest (88% of the study area; Fig. 2). Fens provide an abundance of high-quality graminoid forage (e.g. sedges and rushes; Larter and Gates 1991, Fortin et al. 2003, Strong and Gates 2009) and are easily accessible to bison in winter when frozen. Nahanni bison also selected fens in summer, but to a lesser extent than in winter, consistent with the expectation of tradeoffs between forage, footing and biting insects that may cause bison to seek drier foraging habitats (Belanger et al. 2020) and a broader diet in summer (Jung et al. 2015, Hecker et al. 2021). One unexpected behaviour by Nahanni bison was their strong selection for shrub-dominated habitats, particularly in winter when > 80% of the diet of other boreal bison populations is graminoids (Larter and Gates 1991, Fortin et al. 2003, Jung 2015). This may reflect the low availability of graminoids in our study area, which is supported by consumption of horsetail (Equisetum spp.), despite the damage it causes to their teeth (Larter and Allaire 2007). In Alaska, the Farewell Lake population provides a similar example of dietary plasticity to persist in a graminoid-poor environment, where shrubs comprise up to 94% of their summer diet (Waggoner and Hinkes 1986).

During summer, bison consistently selected fluvial habitats such as riverine gravel bars and islands, where sedges and willows grow (Larter et al. 2003). Nahanni bison frequent the Liard River corridor and regularly cross the river to reach riparian foraging areas, despite the risk of drowning (Larter et al. 2003, Thomas et al. 2021). The strong association of Nahanni bison with riverine habitat is a unique aspect of their ecology, and is likely an adaptation to low forage availability. By comparison, other reintroduced bison populations that border major rivers (e.g. Mackenzie and Ronald Lake populations) are not known to use riparian areas or regularly cross these rivers (Larter and Gates 1991, DeMars et al. 2016). Foraging habitat is relatively abundant in the range of the Ronald Lake population, where graminoid, shrub and herbaceous plants are prevalent (∼15% cover, DeMars et al. 2016, compared to 6% cover in the Nahanni region) and where burned forest provides a rich source of forage (grasses and forbs) across most of the population range. In contrast, the Nahanni region is dominated by dense, mature boreal forest, which offers little forage for bison (Redburn et al. 2008, Strong and Gates 2009). Thus, fluvial habitats may provide much-needed forage for Nahanni bison, given the lack of other suitable foraging habitat within their range. Bison may act as ecological indicators of suitable foraging habitat for other grazing mammals in the boreal forest, and may have a role in maintaining these open foraging habitats, similar to that in more southern locales (Ranglack and du Toit 2015). However, avoidance of biting insects and predators are alternative explanations for their selection of riverine habitat (Larter and Allaire 2007, Thomas et al. 2021), and are not mutually exclusive with forage availability. Biting insects are known to influence habitat selection by forest-dwelling bison (Belanger et al. 2020), but predator avoidance is unlikely to be a strong factor for the Nahanni population given limited observations of predation events.

Figure 3.

Averaged selection coefficients and 95% confidence intervals from resource selection function models of GPS-collared bison Bison bison (n = 10) in Northwest Territories, Canada, during summer and winter (2007–2019). (A) is models for individual bison, while that for (B) is all individuals pooled together. Coefficient averages were derived by an inverse variance weighting method, and confidence intervals were calculated conservatively using a method developed by Marzluff et al. (2004). Abbreviations for model predictors are explained in Table 1.

Table 4.

Comparison of selection coefficients (β) and standard errors (SE) from individual RSF models, inverse variance weighted (IVW) averages (derived from individual models) and mixed-effects models with pooled data for GPS-collared bison Bison bison (n = 10) in the Northwest Territories, Canada, during summer and winter. Individ + and – show the number of individual models with significant positive and negative coefficients, respectively (models with non-significant coefficients are not included in counts). Asterisks (*) indicate IVW or pooled coefficients with 95% confidence intervals that do not overlap zero. Predictor variable abbreviations are explained in Table 1.

Figure 4.

Maps of predicted relative habitat selection by GPS-collared bison Bison bison (n = 10) in Northwest Territories, Canada, during summer and winter. Predictions are calculated from the averaged inverse variance weighted (IVW) selection coefficients of individual RSF models. Predicted habitat selection values for bison are ranked from low (1 – green) to high (10 – red).

Bison generally selected roads during both seasons, as predicted. Roadsides represent patches of high forage biomass (i.e. grasses and forbs), which are rare in the boreal forest, including introduced agronomic species that are attractive to grazing ungulates (Rea 2003, Leverkus 2011, Strong et al. 2013). Additionally, bison may use roads to enhance their movement efficiency (Jung 2017, DeMars et al. 2020), particularly during winter when plowed roads provide a reprieve from deep snow (Bruggeman et al. 2006). In support of forage as a driver of habitat selection, Nahanni bison consistently selected a primary road corridor (Highway 7) in summer when wetlands were less used as foraging habitat, likely because the latter were partially submerged and offered poor footing and high densities of biting insects. In contrast, half of the GPS-collared individuals moved away from the road in winter, instead occupying higher elevation areas with fen complexes. Regardless of the mechanism for road attraction, major roads can be a significant source of mortality for small populations of reintroduced bison. In some years, vehicle collisions have killed up to 15% of the Nordquist population, which occupies a highway corridor in British Columbia and Yukon (Leverkus 2011, Jung 2017), and Nahanni bison are killed each year along Highway 7 (Government of Northwest Territories 2019). Thus, the strong seasonal attraction to primary roads may be an ecological trap.

Many reintroduced bison populations are located in protected areas, including those most extensively studied (e.g. Wood Buffalo, Prince Albert and Yellowstone national parks). Consequently, bison response to landscapes affected by resource extraction is not well understood. The Ronald Lake population occurs in an industrialized landscape and has been shown to avoid anthropogenic disturbances (DeMars et al. 2020), possibly reflecting a greater prevalence of natural foraging habitat within their range. In contrast, Nahanni bison may rely on anthropogenic disturbances to supplement their nutritional needs. Nahanni bison consistently selected artificial clearings associated with forestry and oil and gas development. Although these disturbances were the rarest type of land cover in the study area (0.1%), bison may have been attracted to them because they were likely rich in forage species compared to the surrounding forest (Redburn et al. 2008).

Based on a history of human–bison conflict in local communities, where bison forage on lawns and airstrips (Larter and Allaire 2007, Government of Northwest Territories 2019; Fig. 1), we expected that bison would be strongly attracted to communities. However, we found no evidence of an association with communities at a population level. Bison behaviour was highly variable; half of the collared bison selected communities in summer, whereas the other half either avoided communities or had neutral selection patterns. During winter, most bison showed complete avoidance, never coming within several kilometers of either community. These results imply that only some individuals and groups occupy communities, rather than the population at large, although more collared individuals would be required to affirm this.

Our study had several limitations. Most importantly, inferences were based on a small number of collared bison (n = 10) with relatively short collar lifespans (x̄ = 481 days). A larger sample size would provide greater confidence in our conclusions about population-level selection behaviour, given the substantial inter-animal variation we now know to exist. Unfortunately, GPS collar data is perhaps more difficult to acquire for bison than for other ungulates. Most collars deployed on bison malfunction in < 2 years due to bison behaviour (Jung and Kuba 2015). We discarded data from two individuals as a result of immediate collar failure, and only 3 of 10 individuals had collars that were active for ≥ 2 years (Table 1), which is their battery life expectancy. As such, studies relying on GPS collar data to examine resource selection by bison are challenged to obtain multi-season or multi-annual data from individuals, which is a constraint for modeling an individual's choices over longer time spans. However, the small number of GPS-collared bison in our study likely represented selection by a greater number of animals than represented solely by the collared individuals because forest-dwelling bison are gregarious. In an adjacent bison population mean group size was 18, with significantly larger groups occurring during the summer, often > 50 animals (Jung 2020). In our study area, the mean group size of bison during July was 15.1 ± 13.2 (SD, range = 1–65, n = 174 groups), based on annual surveys from 2002 to 2018 (Larter and Jung, unpubl.).

Two additional limitations of our study are noteworthy: First, all but one of our collared individuals were female, leaving a substantial gap in our knowledge about sex-specific habitat selection for this population. This may be important, as males and females in forested environments typically do not travel together outside of the rut (Jung 2017, 2020). Obtaining additional data about habitat selection by male bison may be exceedingly difficult given males often destroy collars within a few weeks of deployment (Jung and Kuba 2015, Jung et al. 2018b). Lastly, our study was limited by considering only two broad seasons: summer and winter. Previous studies have found important distinctions between habitat selection and diet by bison in early versus late winter (Jung 2015, Jung et al. 2018a), and bison diet composition may vary on a monthly basis (Larter and Gates 1991). Our small sample size prevented us from further data division, but future work should evaluate habitat selection at a finer temporal resolution.

In conclusion, our study provides new information on habitat selection behaviour by threatened bison in a boreal landscape that is seemingly poor in forage and subject to natural resource development. We report individual and seasonal variation in habitat selection, and demonstrate that some land cover types, such as graminoid-dominated wetlands, fluvial habitats, roads and other anthropogenically disturbed areas may be disproportionately used relative to their availability. Our predictive habitat suitability maps serve as an important first approximation of areas of conservation interest for Nahanni bison, and should be considered in land use planning and environmental assessment processes within their range. Habitat selection by Nahanni bison highlight the importance of graminoid-dominated habitats, particularly in winter, which may be limiting in some boreal landscapes. Areas considered for reintroductions of new populations or expansion of existing populations should prioritize identifying and protecting these key habitats to support bison population growth and recovery. More broadly, our results indicate that individual behaviour is an important factor when considering bison habitat selection, as well as the potential for human–bison conflicts. Finally, our work provides another example of determining habitat selection behaviour from a small number of GPS-collared animals that may be used to initially identify areas of highest conservation interest, where resources do not allow for a larger sample of instrumented animals.