To better understand the mechanisms driving the distribution of a threatened carnivore, wolverine Gulo gulo, in the southern Columbia Mountains, we contrasted four hypotheses; climate, food, human disturbance and trapping harvest. We used non-invasive hair snagging methods to examine wolverine occupancy with respect to these factors, collectively and by sex, at a 5 km radius scale around sampling sites and at a larger 10 km radius to evaluate consistency over changes in scale. For all analyses our top models included food and human disturbance. Of the four food items examined; caribou Rangifer tarandus, mountain goat Oreamnos americanus, moose Alces alces and hoary marmot Marmota caligata, wolverine occurrence was most closely related to hoary marmot habitat. With regards to human disturbance, we documented a negative association with forestry road density and a positive association with protected areas. The importance of climate was low compared to the food and disturbance hypotheses. Somewhat surprisingly, recent harvest was positively associated with subsequent wolverine presence. Our top 10 km scale models were similar to the 5 km scale, but with a stronger positive link to caribou distribution. We show that marmot habitat is important to wolverine in winter and suggest that management actions for conservation prioritize factors related to female occurrence, as these were more clearly defined than male. We demonstrate that human disturbance is a major driver of wolverine distribution. Protected areas appear to be providing secure habitat, and reducing road density or mechanized use of roads in winter should be considered for species conservation. A positive spatial relationship with recent harvest, in the context of a heavily harvested population, suggests that harvest data alone may not be useful in detecting population declines.
Effective conservation of threatened species demands an understanding of habitat factors that link to demographic processes or population regulation such as food resources, reproductive requirements and mortality risk (Nielsen et al. 2010, Festa-Bianchet and Apollonio 2003). Habitat quality may be further modified by human disturbance or mortality (Nielsen et al. 2010). Understanding the processes linking life history to demography allows land managers to focus conservation effort on key areas such as food resources or habitat security which should lead to better evaluation of risks and more certain conservation outcomes.
Wolverines Gulo gulo have been assessed as a species of conservation concern in both British Columbia (BC) and in Canada (BCCDC 2017, SARA 2018). These listings were due to low population densities and reproductive rates, while identified threats include climate change, caribou decline, human disturbance, habitat loss and mortality from trapping (COSEWIC 2014). Since divergent factors may limit species ecological relationships across landscapes with different disturbance regimes (Shirk et al. 2014), we investigated four mechanistic factors that have been deemed important, but never explored together to explain wolverine winter distribution in the Southern Columbia Mountains. These factors were climate, food, anthropogenic disturbance and trapping harvest. The effect of trapping as a mechanism influencing distribution in North America (Krebs et al. 2004, Squires et al. 2007) has not been examined previously. The study area, proximal to the southern edge of the current species range, provides a link to US populations, which are at considerable risk (USFWS 2016).
Wolverines are a circumboreal species associated with deep and persistent snowpack (Copeland et al. 2010). The relationship with persistent spring snow has been attributed to the need for snow tunnels in which to den (Magoun and Copeland 1998, May et al. 2012, Magoun et al. 2017), to facilitate movement (Lofroth et al. 2007, Schwartz et al. 2009), and food caching (Inman et al. 2012). Snow cover may be less important in colder, flatter, more northerly latitudes (Copeland et al. 2010, Webb et al. 2016), but we expected that spring snow cover would limit wolverine occurrence in our study area, as it has in other alpine environments (Aubry et al. 2007, Schwartz et al. 2009, Copeland et al. 2010, Heim et al. 2017).
Wolverines are facultative predators and scavengers, relying mostly on small- and medium-sized mammals and ungulates for sustenance in winter in BC (Lofroth et al. 2007). Collectively, caribou Rangifer tarandus, mountain goat Oreamnos americanus, moose Alces alces, hoary marmot Marmota caligata and porcupine Erithizon dorsatum made up 96% of prey items in scat samples in the Columbia Mountains, immediately north of our study area (Lofroth et al. 2007). Winter food resources influence habitat selection (Krebs et al. 2007) and the level of scavenging (Mattisson et al. 2016), reproductive rates (Persson 2005, Rauset et al. 2015) and survival (Krebs et al. 2004). We predicted important winter wolverine food items would influence distribution.
Wolverines are thought to be sensitive to anthropogenic disturbance. Road density (Rowland et al. 2003, Krebs et al. 2007, Scrafford et al. 2018), winter recreation (Krebs et al. 2007), and infrastructure (May et al. 2006, May et al. 2012) have been used as surrogates for disturbance in previous studies. Other studies have shown that wolverine distribution is negatively associated with linear industrial features (Fisher et al. 2013, Heim et al. 2017). While Scrafford et al. (2017) suggested that individual wolverines were attracted to such features in circumstances when there was low human use, he also noted roadside habitats may be population sinks due to high mortality. Forestry roads in our study area are used in winter for snow machine recreation and forest harvest, and rarely, access to industrial sites or homes. We predicted that wolverine occurrence would be negatively related to road density, to reduce encounters with humans or with apex predators such as wolves that use the compacted snow on roads to facilitate movement (Dickie et al. 2017).
Recent analyses of wolverine distribution in Canada's western mountain ranges have documented distributional margins for wolverine that coincide with increased human activity (Fisher et al. 2013, Heim et al. 2017). However, these associations also occur along a longitudinal gradient of decreasing elevation and topography, concurrent with changes in climate and legislative protection. Our study provides a contrasting geographic background to address similar questions, as factors relating to human disturbance occur as a mosaic throughout our study area and across several mountain ranges. For example, while human settlements and road networks tend to occur in low elevation valley bottoms, forestry roads with mechanized recreational and industrial activities are widespread at higher elevations. Wilderness areas such as parks provided a patchwork of low human activity dispersed throughout our sampled area.
Large protected areas (260–2000 km2) in the South Columbia Mountains share similar elevational gradients and topography as surrounding areas with more human use and may also contain traplines and roads. However, unlike their surroundings, these protected areas have difficult access and are largely devoid of human activity in winter; snow machine and helicopter use are prohibited or strictly controlled. Thus, we hypothesized that wolverine occurrence would be associated with protected areas, and that this association would give us the strongest evidence for anthropogenic avoidance.
Consideration of top down mortality sources may be necessary to describe distribution when the level of mortality is substantial (Mowat 2006, Nielsen et al. 2010). Since wolverines are territorial (Bischof et al. 2016) with low population growth rates, mortality can be a strong limiting factor. Some of the areas we sampled were heavily harvested (Lofroth and Ott 2007, Mowat et al. unpubl.) and trapping mortality may result in vacant wolverine territories (Vangen et al. 2001). Consequently, we expected that recent trapping mortality would be negatively associated with wolverine presence (Krebs et al. 2004).
Male and female behavioural differences with respect to parental investment may be reflected in sex-specific factors related to distribution. For example, winter habitat selection differed by sex in the North Columbia Mountains; with male habitat associations best explained by food items, while females were positively associated with food and negatively associated with predation risk and human disturbance (Krebs et al. 2007). With greater restrictions to home range use due to denning, we predicted female occurrence would be linked to spring snow and a more constrained suite of prey items that can be scavenged, captured or cached and accessed during this period of restricted movement. Since females must leave their vulnerable kits in order to forage, they are sensitive to human disturbance near maternal dens (May et al. 2012), and even the presence of skiers has led to den abandonment (Pulliainen 1968). We predicted that female wolverine presence would be related to areas of minimal human disturbance (Krebs et al. 2007, May et al. 2012). Conversely, with lower parental investment, we expected males to be more tolerant of disturbance and more at risk to trapping-induced breaks in distribution.
We considered that wolverine would be distributed according to climate-adapted behaviour and reproductive requirements, in addition to food. Additionally, we hypothesized that anthropogenic impacts and trapper harvest would have negative influences. Our objectives were to simultaneously examine the factors that limit the distribution of wolverine in southeast BC to try to understand the relative importance of environment, disturbance and mortality.
Study area
Our 26 000 km2 study area in southeast BC included parts of the Selkirk, Purcell and Monashee Mountains (Fig. 1). The study area is bisected east–west by Highway 3, which crosses two high elevation passes. Smaller highways follow low elevation valleys. Extensive forest harvest has occurred throughout the area and mining was widespread historically; both industries built many roads. Registered traplines occur throughout the area, and wolverine trapping is allowed from 1 November to 31 January with no harvest quota (Province of BC 2018). Twenty-nine wolverines were trapped in 20 traplines within the study area boundaries during our sampling (2012–2016). Winter recreation (snow machine use, ski resorts, helicopter or snowcat-access skiing, ski lodges and backcountry skiing) is common. Provincial parks and protected areas occupy ∼20% of the area.
The Monashee, Selkirk and Purcell Mountains have elevations ranging from 540 m to 3000 m and receive 150–450 mm of winter precipitation (MacKillop and Ehman 2016). Low elevation forests are composed of western redcedar Thuja plicata, western hemlock Tsuga heterophylla, Douglas fir Pseudotsuga menziesii, Ponderosa pine Pinus ponderosa, lodgepole pine Pinus contorta, trembling aspen Populus tremuloides and western larch Larix occidentalis. At higher elevations, Englemann spruce Picea engelmannii and subalpine fir Abies lasiocarpa forests transition to treeless alpine meadows, rock and ice (MacKillop and Ehman 2016). Potential wolverine foods available through predation or scavenging included large ungulates; mountain goats, mountain caribou, moose, elk Cervus elaphus, mule deer Odocoileus hemionus, white-tailed deer O. virginianus and small mammals; hoary marmots, Columbian ground squirrels Spermophilus columbianus, snowshoe hares Lepus americanus, American pika Ochotona princeps and porcupine (Lofroth et al. 2007).
Methods
Sampling
We partitioned our study area into 10 by 10 km cells, excluding primarily residential and agricultural areas, and sampled each cell using non-invasive genetic techniques during one winter between 2012 and 2016. All 227 cells contained one or sometimes two bait sites, sampled twice in ∼25-day sampling intervals, between January and April. Site locations were selected based on accessibility, topography and local knowledge of wolverine movements. Sites averaged 7 km apart. They consisted of a tree wrapped in barb wire for hair collection with bait (deer head or beaver) nailed 2 m from the ground or snow surface to entice the animal to climb (Fisher et al. 2013). A scent lure was applied to a cloth and hung from a nearby tree. During each check, the barb-wire was examined for hairs and the bait and lure replenished. Hair was collected and stored in paper envelopes in a dry environment.
Hair samples were submitted to Wildlife Genetics International (WGI) for genetic identification analysis. DNA was extracted using QIAGEN DNeasy Tissue kits, following the manufacturer's instructions (Qiagen Inc.). Species identification was based on a sequence-based analysis of a segment of the mitochondrial 16S rRNA gene (Johnson and O'Brien 1997). For samples that yielded wolverine DNA, a ZFX/ZFY microsatellite marker identified sex and a group of 8 microsatellite markers (Gg-4, Gg-7, Ggu-216, Ma-2, Tt-4, Gg-14, MP0182, MP0197) was used to assign individual identity as described in Mulders et al. (2007).
Table 1.
Description of variables, their data sources and their hypothesized influence on wolverine occupancy in winter in the southern Columbia Mountains 2012–2016.
Landscape attributes
Factors affecting animal distributions are scale-dependent and a multi-scale approach should generate more robust conclusions about species–environment relationships (Lindenmayer 2000, Ciarniello et al. 2007, Mayor et al. 2009). Fisher et al. (2011) determined a characteristic scale for wolverine at a 5 km radius around detection points, though larger scales (10 km) may be better (Heim et al. 2017). Using ArcGIS (ESRI ArcGIS 10.4.1), we generated both 5 and 10 km radius buffers (area = 78.5 km2 and 314 km2) around each sampling site. These correspond to the minimum and average home range sizes for females in similar habitats to our study area (Hornocker and Hash 1981, Krebs et al. 2007). We present our 10 km scale data with the caveat that, although modelling terms and parameters should be unaffected, measures of variance and p-values may be underestimated because of the autocorrelation caused by the overlap among adjacent buffers (x̄ = 22%).
Within each buffer we measured the proportion, density or count of each landscape feature. Large lakes were excluded as non-habitat. We produced the spring snow layer using moderate resolution imaging spectroradiometer (MODIS) satellite data (Hall and Riggs 2007, Copeland et al. 2010; Table 1, Fig. 2). We used 17 years (2000–2016) of daily ‘MOD10A1 C6’ snow product images, from 24 April to 15 May. For each image, all values (0–100) representing snow albedo were reclassified as 1 and the value (125) for land as 0. These were combined for each year, on the condition that a pixel retained a value of 1 if it was snow covered for all days, otherwise 0. We summed the number of years (0–17) that a pixel had persistent spring snow coverage, and averaged this value over the extent of the buffer to measure SPRINGSNOW.
We mapped the distribution and habitat quality of preferred wolverine prey to index that food source. We measured caribou winter range (CARIBOU) based on relocations of radio-collared animals. We used the 95% fixed kernel estimates of annual sub-population ranges of resident herds from 1990 to 2000 (Wittmer et al. 2005: Table 1, Fig. 2). Forest-dwelling caribou are known to have high fidelity to home ranges (Faille et al. 2010), however, this measure may overestimate caribou distribution, given the decline in caribou over the interim between that study and ours.
Figure 2.
Thumbnail images of habitat variables in the South Columbia Mountains: (a) SPRINGSNOW, (b) FSR (grey)/OTHERROAD (red, same image), (c) PROTECTED, (d) MOOSE, (e) MARMOT, (f) GOAT, (g) CARIBOU and (h) HARVEST5.
We considered habitat surrogates for the remaining prey species (Table 1, Fig. 2) using locally derived habitat models. For mountain goat (GOAT), we mapped winter habitat using a buffer of 500 m around escape terrain (100% slope) and included south facing slopes (135–185°) between 1330 and 2320 m (Poole and Heard 2003, Poole et al. 2009). We modeled moose winter habitat (MOOSE) as the area below 1300 m within 500 m of riparian areas with canopy cover <50% and shrub cover >40% (Poole and Stuart-Smith 2006). We considered marmot habitat (MARMOT) to include talus or moraine landforms adjacent to alpine meadows on south aspects (135–185°) with slope less than 60% (Turnock et al. 2017). Porcupines are habitat generalists (Snyder and Linhart 1997) and appear to be declining in western North America (Appel et al. 2017). As there were no regionally relevant habitat studies available; we were unable to adequately model porcupine habitat.
We modeled anthropogenic disturbance using both the density of roads (km km–2) and the percent overlap with protected areas (Table 1, Fig. 2). We were interested in the disturbance effect of winter recreation because snow machine, snowcat and helicopter skiing use occurs in high elevation subalpine and alpine areas associated with wolverines. Consequently, we separated our road layer into two components: forestry road (FSR) density, including decommissioned and overgrown roads, and the density of all other roads (OTHERROAD), including paved and unpaved highways and secondary roads. These latter roads are typically at low elevations and associated with residential, commercial, agricultural and transportation infrastructure. Forestry roads are gravel roads built to facilitate forest harvest and extend to high elevations up the drainage systems of almost all named creeks in the area. These roads are only cleared of snow during active logging operations and are the primary conduits for winter recreation in the form of snow machine use. Although the traffic levels are generally low and can vary greatly, we had no way to account for this variation.
Snow machine and other mechanized winter recreation also occur in cut blocks and alpine areas but quantifying these spatially extensive use patterns is difficult. Winter recreation tenures are up to 1500 km2 in size and cover ∼20% of the study area. For heli-skiing especially, the actual area used is highly dependent on the temporally and spatially variable snowpack conditions. In lieu of monitoring these directly, we used protected areas as a proxy for low human use. Protected areas included provincial parks, ecological reserves and privately owned conservation properties (Table 1, Fig. 2). We measured the percent overlap of the 5 and 10 km buffer areas with protected areas as an alternative index of human disturbance which was not correlated to road density or persistent spring snow.
Harvest was estimated for the study area using a combination of compulsory reporting, fur auction reports and carcasses submitted by trappers, most of which did not include a specific kill location (Table 1, Fig. 2). All pixels within a trapline were weighted by the number of animals trapped and for each buffer we calculated the average value of the pixels. HARVEST5 included all harvest within the November–January trapping season immediately previous to our sampling, five months is the longest possible time span between a wolverine removal and the final check of our hair traps. We considered that this short time frame had the highest potential to detect distribution gaps caused by a trapping kill, however, we also evaluated time frames that included >1 harvest season; of 17 (HARVEST17) and 65 months (HARVEST65) previous to sampling efforts, to assess consistency. In the context of the analysis, this time period spanned 2008–2016 and comprised a harvest of 7, 16 and 45 wolverine for 5, 17 and 65 months of harvest prior to sampling, respectively.
Occupancy modelling
To facilitate comparison of effect sizes among variables, we standardized each of our covariates by subtracting the mean value and dividing by the SD (Steel and Torrie 1980). We assessed collinearity and multicollinearity using Spearman's correlation coefficients (Glasser and Winter 1961) and variance inflation factors (VIF; Zuur et al. 2010) using cutoff values of ρ > 0.6 or VIF > 3, respectively. We retained most correlated covariates in the model sets because these variables were of particular interest but we did not include variables that were highly correlated in the same model (Dormann et al. 2013).
We used occupancy as a state variable to model wolverine distribution as a function of habitat characteristics while accounting for imperfect detection (MacKenzie et al. 2002) in a stacked single season framework. We modeled the response variable, wolverine presence, using four sets of data; a 5 km scale for combined sex, female and male models and a 10 km scale for combined sex. We considered only two detection covariates a priori in all models; the number of nights a hair trap was in place (TRAPNIGHT; average 26 nights in each session; range 11–56) and the session in which a trap was checked (SESSION 1 or 2). We expected that the longer the trap was in place, the greater our chance of detecting wolverine, as for other species (Mowat 2006, Lamb et al. 2016). We also expected wolverine detectability to increase with each sampling session (Fisher et al. 2013).
We fit occupancy models using the ‘unmarked’ package in R (Fiske and Chandler 2011, < www.r-project.org>). Our global model included the covariates SPRINGSNOW, CARIBOU, GOAT, MARMOT, MOOSE, FSR, PROTECTED, OTHERROAD and HARVEST (Table 1).
We checked for potential non-linear relationships between occupancy probability and its covariates by plotting univariate response curves. Where a non-linear relationship was suggested, we transformed the original data and selected either the log10 or original variable form, retaining the form of the covariate having the lowest Akaike information criteria for small sample sizes (AICc) in all subsequent analyses (Burnham and Anderson 2002). Next, we constructed candidate models by evaluating all subsets of the global model using the R package MuMIn (Doherty et al. 2012, Barton 2016). We removed candidate models that included variables having parameter estimates opposite to the predicted response, but these variables are acknowledged and examined. We regarded all models within 2 AICc points of the top ranked model as having empirical support (Burnham and Anderson 2002), while possibly containing uninformative parameters (Arnold 2010). We used AICc weights and log-likelihoods to compare the relative support for each model and we calculated parameter estimates for the top-ranked model. We assessed goodness of fit for the top ranked models following MacKenzie and Bailey (2004), using a parametric bootstrap χ2 test with 5000 simulations (R package AICcmodavg, Mazerolle 2017).
Results
We detected wolverine at 71 of 235 sites (Fig. 1) and identified 37 individuals. 16 sites had detections in both sessions. Overall estimates of occupancy and detectability in the study area were 0.50 (SE = 0.09) and 0.36 (SE = 0.07), respectively. Both detection covariates: SESSION and TRAPNIGHT, were positively related to detectability. In general, SESSION was strongly significant and TRAPNIGHT marginally so, supporting our a priori decision to retain these in all models.
Some occupancy variables had collinear relationships. In particular SPRINGSNOW was correlated with MOOSE, GOAT and FSR at both scales as well as OTHERROAD at 10 km. FSR was correlated with GOAT, MOOSE and OTHERROAD at 10 km (Supplementary material Appendix 1 Table A1). We excluded models containing combinations of correlated variables. None of the occupancy variables were correlated with the detection covariate TRAPNIGHT. Additionally, the variables MARMOT, FSR, PROTECTED and GOAT were log10 transformed to linearize their relationship with occupancy.
For all four analysis, top ranked global models fit much better than the constant (NULL) occupancy models (Table 2, see Supplementary material Appendix 1 Table A2, A3 for a complete list). For the climate hypothesis, SPRINGSNOW did not appear as a variable in the top ranked models, although in the univariate SPRINGSNOW model the beta term was positive and strongly significant at both scales and sexes (Table 3). For the food hypotheses, MARMOT was the most important food variable and the coefficient was positive and significant in all top models (Table 2, 4, Fig. 3). CARIBOU appeared in the top models for all data sets except female, the relationship was strongest at the 10 km scale (Table 2, 4). GOAT was present in top models for males and at the combined sex 5 km scale. Human disturbance variables were present in all top models (Table 2). The response to FSR was negative, and significant for every analysis except male. PROTECTED was positive, and significant for all analyses except female, in which the covariate may be uninformative (Table 2, 4, Fig. 3). OTHERROAD also featured in the top ranked models for 5 km male and combined sex scale, although the relationship was weak. The variables HARVEST and MOOSE consistently presented a sign on the covariate that was not congruent with our hypothesis, in both univariate and multivariate models and at multiple time scales for HARVEST (Table 3). We did not include them in our final global model sets.
Male only models had many more terms and the choice of top model is less clear than the female models (Table 2). In particular, the negative association with FSR density was less pronounced, although the positive association with PROTECTED areas was much greater and showed a significant relationship (Table 4). While not in our top models, females had a stronger univariate association with SPRINGSNOW than males (Table 3).
In summary, food and disturbance terms feature prominently for all analyses. MARMOT was present in all top models as was either FSR or PROTECTED or both. SPRINGSNOW was a strong univariate term but not a component of any top models while MOOSE and HARVEST did not match our predicted direction of influence. Our top models all fit the data well: 5 km combined sexes (χ2 = 0.07, p = 0.68, ĉ = 0.17), male (χ2 = 0.12, p = 0.50, ĉ = 0.34), female (χ2 = 0.16, p = 0.25, ĉ = 1.09) and 10 km (χ2 = 0.07, p = 0.67, ĉ = 0.17).
Table 2.
Top occupancy models (ΔAICc ≤ 2) for wolverine in the South Columbia Mountains 2012–2016 at 5 km (overall, sex-specific) and 10 km scales (overall). Degrees of freedom (df), log-likelihood value (LL), Akaike information criteria values for small sample sizes (AICc) and AICc weights (wi) are reported. Models are ranked based on relative AICc values (ΔAICc). All models include p (SESSION +TRAPNIGHT), only occupancy variables are shown.
Table 3.
β parameter estimates with standard errors (SE) and p values for wolverine occupancy variables in univariate climate, food, disturbance and harvest models in the South Columbia Mountains 2012–2016.
Discussion and conclusions
Knowledge of factors shaping wolverine distribution are important for conservation and species management. Sampling from a large mountainous landscape at the periphery of the species contiguous range, we provided multi-scale and sex-specific decomposition of factors influencing wolverine distribution in winter, a time when female have specialized needs related to reproduction. Nonetheless, we found clear consistencies in wolverine response to climate, food, human disturbance and harvest across both spatial scales and sexes.
Climate
Climate influenced wolverine distribution but did not explain occupancy as well as the combination of food and disturbance variables. SPRINGSNOW was correlated with FSR at both scales of analysis, consequently the diminished importance of SPRINGSNOW is likely because the strong FSR variable incorporated information on both climate and disturbance. Notably, if the FSR variable is removed from the analysis, SPRINGSNOW appears in both the female and 10 km top models. Hence, our results do not reject the hypothesis that wolverine occurrence is constrained by an obligate association with persistent spring snow (Aubry et al. 2007, Copeland et al. 2010), but do suggest the alternative explanation that the relationship between spring snow and wolverine distribution could be functionally related to the distribution of food, disturbance or mortality risk.
Table 4.
β parameter estimates with standard errors (SE) and p values for the top ranked wolverine detection and occupancy models in the south Columbia Mountains 2012–2016. Dash (–) indicates that the variable was not contained in the top model.
Food
We found marmot habitat to be a predictor of wolverine occurrence in every top model. Rodent abundance has been associated with reproductive success of wolverine (Landa et al. 1997, Persson 2005) and marmots made up 16–67% of prey items in the scat of females at den sites in BC (Lofroth et al. 2007). As marmots are quiescent during winter, perhaps wolverine were accessing previously cached carcasses or excavating hibernating individuals (Hornocker and Hash 1981). Our observations may also reflect a spurious correlation to another functional process such as den selection, although the similarity between male and female models make this unlikely.
Scavenged ungulates comprise the majority of wolverine winter diet in North America (Lofroth et al. 2007, Inman and Packila 2015, Scrafford and Boyce 2018) and wolverine presence in our combined and male only models were weakly associated with large ungulates. In the North Columbia Mountains, caribou were an important element of the diet of reproductive females (Lofroth et al. 2007), but they did not feature prominently in our female analysis. Caribou numbers have declined precipitously since this earlier work (currently <50 total), which probably limits their availability to many resident wolverine (Fig. 2). However, in spite of low numbers, and restricted spatial extent in the South Columbia Mountains, caribou remained in the top models at the 5 and especially 10 km scales for combined sex and for males. While some other feature of caribou habitat could be driving this relationship, wolverine are important predators on caribou neonates (Gustine et al. 2006), and our data suggest some wolverine still rely on caribou as a winter food source. As the southern mountain population of caribou is assessed as endangered provincially and federally (BCCDC 2017) and the wolverine is also classified as threatened (COSEWIC 2014, SARA 2018), a greater understanding of the amount of caribou in the diet of wolverine is warranted.
Figure 3.
Predicted wolverine occupancy (both sexes) as a function of the occupancy covariates: (a) MARMOT, (b) FSR, (c) PROTECTED, in the top ranked 5 km model. Grey lines indicate 95% confidence intervals. Values along the x-axis have been back-transformed to reflect actual landscape measurements.
Given the paucity of caribou in our study area, mountain goat is the major large ungulate overlapping with high elevation territories in winter. Mountain goats made up 15% of the diet in winter in the North Columbia Mountains during the 1990s (Lofroth et al. 2007). Mountain goat habitat was only weakly associated with wolverine occurrence, but considering the marginal correlation of this variable with forestry road density, we can make no clear conclusion about the relative importance of mountain goat as a winter food source for wolverine.
In our models, moose habitat was negatively associated with wolverine occurrence. Moose comprised the bulk of wolverine winter diet in the North Columbia Mountains (Lofroth et al. 2007) but occur at lower densities in our project area. We expected wolverine occurrence to be related to moose habitat in winter because they are known to eat moose. Our results did not support this prediction so we removed this covariate from the model selection exercise because it obviously was not indexing a functional relationship that we understand at this time. Stronger road and climate variables that are also correlated with moose may be driving this association. Other studies have suggested that a negative relationship with moose habitat may be due to wolf avoidance (Krebs et al. 2007), but we had no means to corroborate this.
Disturbance
Most roads are concentrated around human development at low elevations with less snow, consequently a paucity of wolverine proximal to these features in itself is not clear evidence of adverse anthropogenic effects (Copeland et al. 2007). We did find highways and secondary roads were negatively related to wolverine presence, but these roads were largely peripheral to our sampled areas. However, forestry roads extend into high elevation wolverine habitat (Fig. 2). Density of forestry roads was a strong negative predictor of wolverine distribution, and in particular of females. Wolverine may travel on forestry roads in regions of low human use (Copeland et al. 2007, Scrafford et al. 2017), but because of the high level of snow machine operation in our study area, we believe the negative relationship with these roads is indicative of anthropogenic disturbance. Most sampling sites that were adjacent to a forestry road showed evidence of recent mechanized activity. We support Krebs et al. (2007) in cautioning that forestry roads and the associated winter recreation could be impacting habitat quality by disturbing wolverine. Although use of packed roads, such as snow machine trails, by other large carnivores for travel may also play a role (Dickie et al. 2017), wolves and cougars are seldom detected on roads above 1200 m during winter. For this reason we do not believe wolf or cougar avoidance is a major factor driving this relationship. Further, trapping kills were not negatively related to distribution, hence our results suggest the explanation for this pattern of distribution was habitat avoidance due to disturbance, rather than trapping.
Wolverine avoidance of mechanized human disturbance was also corroborated by our protected areas variable. Unlike road density, protected areas were not correlated with any other variable. Protected areas were strongly and positively related to wolverine presence in the top-ranked models for both scales. Since we also found a positive and weaker response to harvest, it seems unlikely that this effect was principally due to reduced levels of trapping in protected areas. Although factors such as lower harvest rates and minimal habitat alteration may play a role, in this study area the primary difference between protected areas and the surrounding landscapes in winter is lower human use. We consider our results to add to an accumulating body of evidence documenting human disturbance impacts on wolverine (May et al. 2006, Krebs et al. 2007, Fisher et al. 2013, Heim et al. 2017, Scrafford et al. 2018).
Harvest
Hunting and trapping are major sources of mortality for wolverine populations in North America (Krebs et al. 2004, Squires et al. 2007), and mortality rates were likely not sustainable during our study (at least 10%/year; Mowat et al. unpubl.). Harvested wolverine populations depend on immigration from un-trapped areas for persistence (Krebs et al. 2004, Sæther et al. 2005, Squires et al. 2007, Dalerum et al. 2008) but our index of mortality did not support a negative influence of trapping on the short-term distribution of wolverine. It is possible that the small number of trapped individuals and relatively large trapline areas, 115–550 km2, precluded sufficient resolution to detect distribution gaps, although our results were consistently positive at multiple scales. Obviously, wolverine must be initially present for harvest to occur. A positive relationship with occupancy may indicate that wolverine presence was driving the harvest variable rather than the converse. Since traplines are seldom larger than wolverine home ranges, and wolverine home ranges are exclusive to same sex conspecifics (Persson et al. 2010), persistent wolverine presence in recently trapped areas implies either: 1) where male and female territories overlap, one of the pair has been trapped, with the other remaining; 2) the presence of juvenile offspring, or; 3) the rapid recolonization of the empty territory. Some combination of all three is to be expected, however, since trapper harvest is strongly biased toward young males (Kukka et al. 2017) but the relationship of recent harvest with occupancy was not female biased, it doesn't appear that we are typically detecting the female of a resident pair, where the male has been harvested, as in 1). Consequently, we suspect a combination of 2) and 3) is more probable. Since harvest likely exceeds productivity, and females seldom breed in consecutive years (Persson et al. 2006), while only a few traplines were responsible for exerting heavy mortality in a limited area over successive years, we further suspect that 3) is weighted more heavily than 2). We speculate that rapid recolonization is a consequence of high quality vacant habitat in these regularly trapped areas, thereby creating population sink dynamics (Delibes et al. 2001, Battin 2004). This is congruent with observations by Aronsson and Persson (2018) that suggest strong intraspecific competition for high quality territories. If the posited mechanism is correct, the observed trend towards a weaker positive relationship between harvest and occurrence over longer time scales suggests a shrinking distribution, particularly for females, i.e. recent harvest is a much better spatial predictor of wolverine occurrence than harvest from as much as five years ago. Clearly, the relationship between wolverine distribution, harvest, habitat quality and population dynamics merits further investigation.
Sex
Our data were mostly consistent with our hypotheses regarding differences between male and female distribution. Females did show a stronger association with spring snow extent and had a more constrained association with prey items, as predicted from denning requirements. Female presence was also strongly related to areas with minimal human disturbance, as measured by forestry road density. We expected males to be more tolerant of human disturbance; although the relationship with forestry roads was only weakly negative, surprisingly, they were strongly associated with our other human disturbance variable, protected areas. As roads occur at finer scales than protected areas, we surmise that these differences could be related to sex-based differences in home range size (male: 800 km2 versus female: 300 km2; Krebs et al. 2007). Regardless, areas of low human disturbance were strong predictors of wolverine presence for both sexes. We expected male distribution to have a more negative association with recent harvest than females, as males are harvested at approximately 3× the rate of females (Kukka et al. 2017). However, this relationship was positive and similar between the sexes. Overall, female distribution was characterized by a few well-defined factors, whereas males had numerous candidate variables in the top models with none clearly better that the other.
Our top models indicate robust relationships to multiple processes that were hypothesized to influence the distribution of wolverine. At both scales and sexes, an association with persistent spring snow was largely eclipsed by a strong link to marmot distribution and a negative relationship with anthropogenic disturbance. A consistent response across two different measures of disturbance, protected areas and roads, lends convincing support for the anthropogenic disturbance hypothesis. In particular, the forestry road variable was compelling as a measure of disturbance within preferred subalpine habitat, as typically human disturbance occurs in low elevation habitat peripheral to core home ranges (Copeland et al. 2007). In a heavily harvested population, a positive association with harvest suggests rapid recolonization rates and source–sink dynamics.
Management implications
Although distribution was explained by multiple factors, the obvious management strategies would involve human disturbance and harvest. First, factors influencing female wolverine distribution were more defined than those of males, and due to requirements for denning and reproduction, these needs should be prioritized by resource managers. Mapping female habitat associations in late winter could be a useful tool in identifying candidate areas for protection. Since parks and conservations areas in the region positively impacted wolverine occurrence, protection of vulnerable wolverine habitat is likely to reap observable benefits. Second, this paper substantiates mounting evidence of road density impacts on wolverine winter distribution, suggesting that access management of forestry roads is likely to have positive implications for species conservation (Lamb et al. 2018). Precluding or limiting mechanized vehicle use on forestry roads may be effective in restoring wolverine habitat use during the later-winter denning period. This could entail decommissioning and restoring forestry roads after resource extraction, or closing roads seasonally or permanently. Our work suggests that priority should be put on restricting access in areas with marmot habitat and that access restrictions to reduce the disturbance of caribou will likely benefit wolverine. Our results also point to the need for consideration of the cumulative impacts of the many uses on the landscape. Future research should focus on wolverine movements with respect to the intensity of snow machine and other recreational use.
Finally, wolverine are currently monitored on the basis of annual fur harvest and we describe a positive relationship between recent harvest and subsequent wolverine occurrence in a heavily harvested population, suggesting that harvest data alone may not be useful in detecting population declines. Additional methods should be employed to monitor wolverine populations. If trapped areas are indeed functioning as mortality sinks in attractive habitat, careful monitoring is critical for population sustainability.
Acknowledgements
The authors would like to thank the regional trapping community for wolverine carcasses, bait and assisting in field operations. Thanks to B. Dorsey and K. McGuinness for GIS work. Thanks to field technicians B. Phillips, T. Abraham, J. Robbins, C. Lehman, S. Forrest, D. Fear, S. Himmer, A. Page, A. Bourelle, D. Lynch, C. Hiebert and J. McCullough and many volunteers. We thank Jeff Parker of Kootenay Valley Helicopters Ltd. for transportation and field assistance. Thanks to staff at WGI for genetic analysis and L. Smit for preparing the harvest data and Fig. 1. We thank J. Krebs and A. Clevenger for reviewing the manuscript.
Funding – Funding was provided by the Columbia Basin Trust and the Fish and Wildlife Compensation Program, as well as BC Ministry of Forests, Lands and Natural Resource Operations (MFLNRO) and the Wolverine Foundation. MFLNRO, BC Provincial Parks, Idaho Fish and Game and Idaho Panhandle National Forest also provided logistical support and assistance.
Permits – The authors extend their thanks to MFLNRO, BC Parks and the Nature Conservancy of Canada for research approvals in their respective areas.
References
Appendices
Supplementary material (available online as Appendix wlb-00480 at < www.wildlifebiology.org/appendix/wlb00480>). Appendix 1.
