Wolverines in winter: indirect habitat loss and functional responses to backcountry recreation

. Outdoor recreation is increasingly recognized to impact nature and wildlife, yet few studies have examined recreation within large natural landscapes that are critical habitat to some of our most rare and potentially disturbance-sensitive species. Over six winters (2010 – 2015) and four study areas ( > 1.1 million ha) in Idaho, Wyoming, and Montana, we studied the responses of wolverines ( Gulo gulo ) to backcoun-try winter recreation. We ﬁ t Global Positioning System (GPS) collars to 24 individual wolverines and acquired > 54,000 GPS locations over 39 animal-years during winter (January – April). Simultaneously, we monitored winter recreation, collecting ~ 6000 GPS tracks ( ~ 200,000 km) from backcountry recreationists. We combined the GPS tracks with trail use counts and aerial recreation surveys to map the extent and relative intensity of motorized and non-motorized recreation. We integrated our wolverine and backcountry recreation data to (1) assess patterns of wolverine habitat selection and (2) evaluate the effect of backcoun-try recreation on wolverine habitat relationships. We used resource selection functions to model habitat selection of male and female wolverines within their home ranges. We ﬁ rst modeled habitat selection for environmental covariates to understand male and female habitat use then incorporated winter recreation covariates. We assessed the potential for indirect habitat loss from winter recreation and tested for functional responses of wolverines to differing levels and types of recreation. Motorized recreation occurred at higher intensity across a larger footprint than non-motorized recreation in most wolverine home ranges. Wolverines avoided areas of both motorized and non-motorized winter recreation with off-road recreation eliciting a stronger response than road-based recreation. Female wolverines exhibited stronger avoidance of off-road motorized recreation and experienced higher indirect habitat loss than male wolverines. Wolverines showed negative functional responses to the level of recreation exposure within the home range, with female wolverines showing the strongest functional response to motorized winter recreation. We suggest indirect habitat loss, particularly to females, could be of concern in areas with higher recreation levels. We speculate that the potential for backcountry winter recreation to affect wolverines may increase under climate change if reduced snow pack concentrates winter recreationists and wolverines in the remaining areas of persistent snow cover.


INTRODUCTION
Fostering societal appreciation for the conservation of nature partly relies upon individuals connecting to nature during leisure activities, which includes participating in outdoor recreation activities (Teisl andO'Brien 2003, Gifford andNilsson 2014). Snow-based recreation during the winter is an important component of that outdoor recreation. In recent years, technological advancements in over-snow equipment including more powerful snowmobiles and lightweight backcountry ski gear provide increasing opportunity for winter recreation enthusiasts to access previously remote backcountry landscapes and have resulted in an increase in human presence in a landscape that has previously been de facto winter wilderness. Indeed, backcountry winter recreation has become valuable both economically and culturally for many small communities (Scott et al. 2008).
Unfortunately, recreation activities can negatively impact landscapes and the wildlife that reside in them (Steven et al. 2011, Sato et al. 2013, Larson et al. 2016. The most commonly reported wildlife responses to recreation are behavioral and physiological, including elevated stress hormones and avoidance or displacement from areas of disturbance (Harris et al. 2014, Arlettaz et al. 2015, Larson et al. 2016. Avoidance of disturbed areas may lead to indirect habitat loss (Patthey et al. 2008, Polfus et al. 2011, Coppes et al. 2017b, the impacts of which could be compounded during winter seasons if animals face increased energetic demands for thermoregulation and travel over snow with limited food supplies (Telfer and Kelsall 1979, Parker et al. 1984, Neumann et al. 2009). Habitat displacement and indirect habitat loss from winter recreation activities have been documented in a few montane and alpine species. In Europe, for example, high elevation forest grouse (Tetrao sp.) are negatively impacted by backcountry winter recreation including habitat displacement as well as energetic and physiological effects (Patthey et al. 2008, Braunisch et al. 2011, Arlettaz et al. 2015, Coppes et al. 2017b). Many species of large herbivore (e.g., red deer, Cervus elaphus; mountain caribou, Rangifer tarandus caribou; bighorn sheep, Ovis canadensis; mountain goat, Oreamnos americanus; moose, Alces alces) have exhibited negative physiological or behavioral responses including indirect habitat loss through avoidance of motorized and non-motorized winter recreation (Seip et al. 2007, Neumann et al. 2009, Courtemanch 2014, Richard and Cote 2016, Coppes et al. 2017a, Lesmerises et al. 2018. Although useful, many of the previous studies assessing the effects of winter recreation on wildlife have been limited spatially and temporally, and most were focused within a single study area and on a single form of winter recreation (Larson et al. 2016). As backcountry winter recreation grows in intensity and spatial extent, coupled with the potential concentration of activities due to climate change-driven reductions in snow area and season (Dawson et al. 2013, Rutty et al. 2015, there is a growing need to understand the impacts of recreation on wildlife species, and particularly on those that are sensitive, snowassociated, and occupy alpine habitats.
Large carnivores are globally threatened and have experienced negative effects of humancaused habitat loss and fragmentation throughout their range (Ripple et al. 2014). In North America, the Rocky Mountains represent a large carnivore hotspot (Noss et al. 1996, Laliberte andRipple 2004), where some species are restricted to high elevation habitat. The wolverine (Gulo gulo) is limited to northern latitudes across its circumpolar distribution and is closely associated with snow and boreal forests, subalpine or alpine habitats (Magoun and Copeland 1998, Aubry et al. 2007, Copeland et al. 2010. Consequently, there is high potential for overlap and interactions between wolverines and backcountry winter recreationists because they both frequent similar areas, that is, areas with deep and persistent snow. Wolverines are also a species of conservation concern throughout much of their expansive range, further highlighting the importance of assessing interactions between wolverine and winter recreation. Wolverines may be vulnerable to direct and indirect impacts of recreation during winter, as they naturally occur at low densities, have low reproductive rates, and remain active through the winter (Hash 1987, Persson 2005, Persson et al. 2006, Copeland et al. 2017). There has been no effort focused on understanding wolverine responses to winter recreation, though some research suggests they are sensitive to human activities and infrastructure (May et al. 2006, Krebs et al. 2007, Stewart et al. 2016, Heim et al. 2017. Females enter reproductive dens within deep snowpack during the winter recreation season with kits born in mid-February to early March, and they occupy these dens through late April or mid-May (Hash 1987, Magoun and Copeland 1998, Persson et al. 2006, Copeland et al. 2010, Inman et al. 2012b). The potential impact of backcountry winter recreation to denning females is of primary concern (Carroll et al. 2001, May et al. 2006, Krebs et al. 2007). In Canada, wolverine status was changed to Special Concern in 2014 with increased winter recreation use combined with sensitivity of denning females among the considerations (COSEWIC 2014). In the United States, wolverines are being considered for listing under the Endangered Species Act, with the most recent status review (U.S. Fish and Wildlife Service 2013) indicating a lack of evidence to assess potential effects of winter recreation.
Understanding the responses of elusive, lowdensity wildlife species to relatively novel human uses such as backcountry winter recreation require innovative approaches that capture the spatio-temporal variability inherent in human activity and the responses of animals to this disturbance (Tablado andLukas 2017, Squires et al. 2018). Over six years, we monitored the movements and habitat use of wolverines in four different study areas in the Rocky Mountains of Idaho, Wyoming, and Montana. Simultaneously, we tracked and monitored winter recreation to characterize the spatial extent and relative intensity of recreation across the landscape. We predicted that wolverine responses to winter recreation would be influenced by the type, spatial extent, and intensity of the recreation. We developed resources selection analyses to both understand wolverine habitat selection within home ranges and to test wolverine responses to winter recreation. These analyses allowed us to evaluate the potential for indirect habitat loss due to winter recreation (Johnson et al. 2005, Polfus et al. 2011, Hebblewhite et al. 2014. While resource selection analyses provide an estimate of average responses, they tell us little about how wolverine responses may change based on the level of exposure to winter recreation (Mysterud andIms 1998, Hebblewhite andMerrill 2008). Functional responses in habitat selection can provide important insight concerning behavioral changes in animals as they experience differing levels of a resource or disturbance (Hebblewhite and Merrill 2008, Moreau et al. 2012, Holbrook et al. 2017. We tested for functional responses in habitat use of wolverines by evaluating how wolverines changed their use of increasingly recreated areas. The goals of our research were threefold (1) characterize fine-scale (i.e., third-order home range scale, Johnson 1980) habitat use and selection of male and female wolverines; (2) assess the importance of motorized and non-motorized winter recreation in influencing wolverine habitat selection and predicted use; and (3) test whether the responses of wolverines to winter recreation were dependent upon the relative intensity of the recreation within individual home ranges.

Overview
We fit GPS collars on wolverines to monitor responses to winter recreation and other resources in mid-and late winter (January-March) and concurrently sampled the spatial patterns of winter recreationists. We developed wolverine resource selection functions (RSF) with a use: availability design to estimate the relative probability of selection (Manly et al. 2002, Johnson et al. 2006, McDonald 2013 including models with and without winter recreation covariates. Based on the selected models, we assessed the effect of winter recreation on wolverine habitat selection and evaluated indirect habitat loss from winter recreation. Finally, we tested whether wolverines showed functional responses to winter recreation based on the relative intensity of winter recreation to which they were exposed. We used ArcGIS (ArcGIS Desktop: Release 10.1-10.5; ESRI, Redlands, California, USA) and R (R Core Team 2016) for data management and analyses.
McCall study area (Payette NF, northern Boise NF); Sawtooth study area (including portions of the Sawtooth NF, southern Boise NF); West Yellowstone study area (including portions of the Caribou-Targhee NF, Custer-Gallatin NF, and Beaverhead-Deerlodge NF), and the Teton study area (including portions of the Caribou-Targhee NF, Bridger-Teton NF, and the Grand Teton National Park). Each study area was a popular backcountry winter recreation destination with backcountry snowmobiling, skiing, or both occurring in the range of wolverines. Each study area also contained large areas without intense winter human activity. Study areas were primarily U.S. Forest Service lands, but also contained a mix of other state and federal land designations. Topography was mountainous with alpine dominated by rock, ice, and low-growing herbaceous vegetation, transitioning into more open conifers with open rocky or subalpine shrub, grass, and herbaceous vegetation. Mid-elevation vegetation was dominated by coniferous forests, with interspersed deciduous tree and shrub communities. The lower boundaries of the study areas were defined by the lower limits of wolverine use, typically near the lower limit of forested habitats, with rare agricultural and sagebrush steppe near these margins.
Infrastructure supporting backcountry recreation varied across the study areas. All study areas had maintained parking areas for backcountry access at trailheads or along plowed roads, and some study areas had a network of groomed snowmobile trails. Within wolverine home ranges, roads were almost exclusively secondary roads that were not plowed for vehicle travel though some were groomed for snowmobile use. The few plowed roads occurred near home range boundaries. All roads were snow-covered during our study, and motorized and non-motorized recreation use was allowed on most roads regardless of whether they were groomed for recreation use. Winter recreation activities varied in the number of recreationists and types of recreation, and each study area had a unique combination of backcountry recreation including snowmobile, ski (including snowboards), snowmobile-accessed ski/board (hybrid), cat-ski, heliski, and yurt-supported ski. The McCall, Sawtooth, and Teton study areas also had developed ski resorts which allowed for backcountry or outof-bounds skiing.

Wolverine capture and monitoring
Between 2010 and 2015, we captured wolverines from early January through April using modified box traps built from logs (Lofroth et al. 2008) baited with road-kill deer or trappercaught beaver and a skunk-based lure. Each trap was equipped with a satellite device that notified us when the trap was triggered (Vectronics trap transmitters TT2, TT3; Vectronic Aerospace GmbH, Berlin, Germany), as well as a VHFbased trap trigger (Telonics trapsite transmitters, TBT series; Telonics, Mesa, Arizona, USA); traps were visited immediately if triggered and maintained every 3-5 d. Traps were closed late February to late March to avoid capturing a reproducing female and re-opened in late March through April for collar removal. Wolverines were anesthetized using a 10 mg/kg ketamine hydrochloride and 0.075 mg/kg medetomidine mixture (Fahlman et al. 2008) delivered by a jab stick. A GPS collar (either WildCellSL collar from Lotek Wireless, Newmarket, Ontario, Canada, or Quantum 4000 collar from Telemetry Solutions, Concord, California, USA) was attached and programmed to collect a location every 20 min on weekends (Saturday, Sunday) and mid-week (Tuesday, Wednesday), which we expected to differ in intensity of human use. Collars were modified with a cotton strip designed to rot away within 4-6 months if we were unable to recapture the animal. Trapping and handling procedures were approved through the University of Montana Institutional Animal Care and Use Committee (IACUC; Permit #055-10MHECS-113010) and the National Park Service IACUC under a research permit (GRTE-2015-SCI-0003). We also obtained research permits through Idaho Department of Fish and Game (IDFG Scientific Research Permit #091210) and Wyoming Game and Fish (WGF Collection Permit #33-928). We monitored the status of wolverines through aerial telemetry flights, including location, denning status, survival, and confirming collar function.

Resource selection function analyses
Resource selection functions compare covariates at used GPS locations with random locations (putatively available) to identify covariates that are used disproportionately more (i.e., selected) or less (i.e., avoided) than available or proportional to availability (lack of selection: Manly et al. 2002). We used general linear mixed-effects models with a logit link function (GLMM) and animal-year as a random effect to control for repeated sampling of individuals (Gillies et al. 2006). The mixed-effects RSF model therefore takes the form: where x n are covariate values for location i of animal-year j with the fixed regression coefficient b n ; c 0j is the random intercept for animal-year j and is e ij is the residual variance within each animalyear. Logistic regression (Hosmer et al. 2013) was used to fit the exponential approximation to an inhomogeneous spatial-point process model, but without the intercept because in used-available designs the true amount of non-use is unknown (McDonald 2013). Thus, the resultant probability is best considered a relative probability of selection or use (Boyce et al. 2002, Lele et al. 2013). Animal and random (available) locations were attributed with the environmental and winter recreation covariates (see Environmental variables, Winter recreation sampling and mapping), which were then standardized ((value À mean)/SD) to support model fitting and allow for comparisons between model coefficients (Hosmer et al. 2013).

Location data and winter season home range analyses
We defined available habitat by estimating winter season home range or seasonal use area boundaries using a local convex hull non-parametric kernel method (Getz et al. 2007) with a fixed "k" number of nearest neighbors. We buffered boundaries by the sex-specific median step length (331 m for females, 441 m for males) to account for habitat immediately available to the animal. We included an individual animal-year for each wolverine with ≥5 weeks of GPS monitoring for any given winter. Data for individuals that exhibited localized or home range-type movements but were monitored for <5 weeks were withheld for model validation; subadults exhibiting exploratory or dispersal type ❖ www.esajournals.org movements were removed from all analyses. Within each home range, we estimated available habitat with random locations generated at a ratio of 2:1 random:use with random locations forced to be ≥30 m apart.
The time wolverine spent under snow and structures resulted in low GPS fix-rates and potential behavioral or habitat-induced bias (Frair et al. 2004, Mattisson et al. 2010. To account for behavior-based missed locations, we developed a modification of Knopff et al. (2009) to identify clusters of wolverine locations based on their spatial (within 25 m of each other) and temporal (within 24 h of each other) proximity. Missed locations were associated with a known cluster site if the location before or after the failed GPS attempt was within a cluster, and the cluster centroid was imputed for their location (Frair et al. 2004). Locations <100 m of an active trap site were censored given the effect of baited traps.

Environmental variables
We evaluated land cover, topographic, snow, climate, and anthropogenic covariates (Appendix S1: Table S1) that may be important predictors of wolverine resource selection at the third order. First, we identified the spatial scale at which each potential covariate was most strongly selected by wolverines (DeCesare et al. 2012; Appendix S1: Table S1). Second, we screened covariates for collinearity (|r| ≥ 0.6), and the covariate with the lowest univariate Akaike's information criterion (AIC) was retained (Hosmer et al. 2013). We also evaluated potential interactions. Finally, we evaluated covariates for potential non-linear relationships using general additive models (Hilbe 2015) and by testing potential non-linear models, keeping the form of the covariate with the lowest AIC (Hosmer et al. 2013). This resulted in slope being included in a quadratic form.

Winter recreation sampling and mapping
We developed spatially explicit maps of winter recreation by sampling backcountry recreation using three methods: GPS tracking of volunteer recreationists , infra-red trail use counters, and aerial surveys. We combined spatial information from GPS tracks with counts of recreational use from trail counters to develop maps of winter recreation intensity. We used the aerial surveys to validate recreation maps (Appendix S2).
To collect GPS tracks of recreation, we sampled recreationists at known recreation access points during mid-week (Tuesday, Wednesday) and weekend (Saturday-Sunday) days from mid-January through mid-April. We sampled recreation systematically, not in proportion to recreation use at access points or across study areas. We asked recreation groups to carry one GPS unit (Qstarz International, Taiwan, ROC, model BT-Q1300, 1 location/5 s, position accuracy <10 m) per ≤4 people in the group, and we recorded the type of winter recreation and the group size per GPS unit. We also distributed GPS units to backcountry guide, heli-ski, and cat-ski operations, with guides carrying the GPS units and recording their group size. To estimate the number of recreationists accessing each study area, we installed infra-red trail counters (Trafx Research Ltd, Canmore, Alberta, Canada) at constriction points on backcountry snowmobile and ski/ snowboard access routes. If the access route was used by both outgoing and incoming recreationists, the counts were divided by two to estimate the one-way traffic.
We developed maps of different types of backcountry winter recreation, including linear travel (primarily access routes along forest roads) and dispersed (off-road) use. We calculated the relative density or intensity of dispersed use (meters of track/100 m 2 ) based on the GPS tracks of recreationists. To account for differences in overall use within and between study areas, we weighted each GPS track based on the proportion of the estimated total recreation use it represented from each trailhead or access point, with total use estimated from the trail use counters associated with the access point (Appendix S2). The GPS tracks of recreationists that use motorized access (e.g., snowmobile, cat-ski, heli-ski) to undertake non-motorized activities were split into their motorized and non-motorized components. For heli-ski GPS tracks, we used only the non-motorized portions of GPS tracks and discarded the track associated with the helicopter transport; any helicopter-based disturbance was not accounted for in our analyses. To test for wolverine responses to spatial pattern and intensity of winter recreation, we developed maps of recreation that became candidates for inclusion ❖ www.esajournals.org as covariates in the wolverine habitat models: (1) the recreation footprint as a binomial characterization of recreation extent that includes roadbased and dispersed recreation; (2) linear recreation along roads and groomed trails; (3) the relative intensity of all winter recreation; (4) the relative intensity of off-road or dispersed recreation (tracks >30 m from a road or groomed route) recreation; (5) the relative intensity of dispersed motorized; and (6) the relative intensity of dispersed non-motorized recreation (Appendix S2).

Model selection
To assess wolverine responses to winter recreation, we first developed RSFs (habitat models) based on environmental covariates not including recreation, which predicts potential habitat quality in the absence of recreation based on relative probability of use (Polfus et al. 2011, Trainor andSchmitz 2014). Then, we added winter recreation covariates to the potential habitat model(s) to test for responses of wolverines to different characteristics of winter recreation (e.g., recreation footprint, relative intensity, recreation type) and to identify the best model to predict "realized" habitat quality accounting for effects of winter recreation on wolverine habitat selection. We followed a two-step process to identify the environmental predictors of wolverine habitat use for all animals combined (global model), for females (female model), and for males (male model). To identify the most predictive of the potential covariates and covariate interactions, we used fixed-effect least absolute shrinkage and operator selection (LASSO) logistic regression (Tibshirani 1996, Reineking andSchr€ oder 2006) implemented using the glmnet package in R (Friedman et al. 2010) for male, female, and global (male and female combined) models. We removed covariates that were not within the selected covariate set (penalty strength set within one standard error [SE] of the minimum cross-validated error; Friedman et al. 2010). In the second step, we used the covariates identified in step one to developed RSF global, female, and male models using GLMM with animal-year as a random effect using the lme4 package in R (Bates et al. 2015). To determine whether a single global model or separate sex-based models were supported, we compared the summed AIC scores of the male and female RSF models to the global RSF AIC; this is possible because the combined male and female data are exactly the full global data (Burnham and Anderson 1998). To include winter recreation effects, we then developed five additional RSF models that included the potential habitat RSF covariates and different combinations of the six winter recreation covariates from our winter recreation maps. We selected the model with the lowest AIC to best represent realized wolverine habitat use in areas that also have winter recreation. For the selected models of potential habitat and realized habitat, we used 10-fold cross-validation to assess the goodness of model fit (Boyce et al. 2002). We also validated the models using out-of-sample GPS location data from wolverine animal-years not used in the development of habitat models to determine how our models predicted the frequency of wolverine use (DeCesare et al. 2012, Holbrook et al. 2017).

Comparing potential and realized habitat quality
We estimated habitat degradation due to winter recreation by calculating the reduction in habitat quality between the potential habitat and realized habitat models (Johnson et al. 2005, Polfus et al. 2011, Hebblewhite et al. 2014). This may underestimate the influence of winter recreation on wolverines because we assume the influence of winter recreation is independent of environmental variables and did not confound modeled relationships. To assess this assumption, we calculated the relative percent change between the potential and realized environment coefficients and identified those covariates with greater than a 20% change in value (Hosmer et al. 2013: equation 3.9). If model coefficients were stable in the potential and realized models (<20% change), this suggests that recreation and the environmental covariates were not confounded.
Each model was mapped at a 30 m resolution, and mapped values were binned into 10 quantiles from low to high quality (i.e., relative probability of use). We classified habitat quality into three groups: (1) the top 30% of the area (bins 8-10) as high-quality habitat, (2) the next 30% (bins 5-7) as moderate quality habitat, and (3) the lowest 40% of habitat values (bins 1-4) as low-quality habitat. We did not include areas where gaps in winter recreation monitoring information did ❖ www.esajournals.org not allow us to predict the probability of use. Indirect habitat loss was calculated as the spatially explicit reduction in habitat quality when comparing the realized habitat maps to the potential habitat maps (Johnson et al. 2005, Polfus et al. 2011. We calculated the degree of habitat degradation by the number of classes reduced, with the most severe degradation indicated by high-quality habitat that is degraded by two classes to low-quality habitat.

Functional responses to winter recreation
We tested whether wolverines exhibited a functional response to the relative intensity of motorized and non-motorized dispersed winter recreation by evaluating how habitat use of recreated areas changes with availability of these areas. If there is no functional response, habitat use of recreation changes in proportion to availability, while deviations from proportional use indicate a functional response (Holbrook et al. 2017(Holbrook et al. , 2019. We calculated the mean recreation intensity at used (animal) and available (random) locations for each animal-year home range and used these data in the following model: where R indicates the recreation type (motorized or non-motorized); U R i = the average recreation intensity at used locations of each animal-year i; b 0 = y-intercept, b R = slope of the functional response; and A R i = the average recreation intensity at available locations within the home range of animal-year i. The null expectation is b R = 1 (proportional use), while b R < 1 indicates decreasing use and b R > 1 indicates increasing use as availability increases. We limited the scope of our functional response analyses to wolverine responses to recreation type and intensity as the primary focus of this work; functional responses to other covariates may also exist.

Wolverine trapping and location data
We captured and GPS-collared 24 individual wolverines (11 females, 13 males) over five years of live-trapping (Fig. 1). We did not capture any female wolverines in the Tetons study area. Each wolverine was monitored for 1-4 yr for a total of 39 animal-years. We obtained >5 weeks of data from 18 (10 females, 8 males) animals over 25 animal-years, averaging 2101 locations/animalyear between mid-January and end of March (Table 1). An additional nine animal-years (five female animal-years and four male animal-years) were used for model validation. The cluster analysis identified groups of animal locations with an average (AESD) distance between an animal location and the cluster center of 18 AE 24 m, and we estimated missed locations associated with a cluster as the cluster center. Raw fix-rates were 75.8%, yet 78% of failed GPS attempts were associated with clustered behavior and were thus imputed. Our corrected fix-rate was 94.7%, providing 53,301 locations used in the spatial modeling and 6603 for model validation. The average size of female winter home ranges was smaller than male winter home ranges (Table 1).

Recreation monitoring
Study areas had two years of GPS-based recreation tracking, infra-red trail use counts, and aerial surveys, though areas without successful wolverine identification may have had less effort. We recorded 5899 GPS tracks (i.e., trips by recreationists) of combined length of 198,019 km (Table 2). While we recorded a diversity of backcountry recreation types (Appendix S3: Table S1), snowmobiling was the most popular motorized backcountry recreation while skiing was the most popular non-motorized recreation. Over 90% of non-motorized recreation tracks were collected in the Teton study area, with localized areas of non-motorized recreation in other study areas (Table 2). Snowmobiling was a common recreation activity across all study areas, and snowmobile tracks were longer (average of 60 km) than ski tracks (average of 10 km); snowmobile tracks constituted 82% of our total track length. Heli-ski recreation only occurred within our Sawtooth study area, and cat-ski recreation was only present in the McCall study area. We established trail use counters at 25 sites. The total estimated recreation visits varied across our study areas from <7000 to >23,000 ( Table 2). The proportion of recreationists sampled using GPS tracking also varied, based partly on the total recreation use and on localized access patterns, from 15% to 42% (Table 2).
Winter recreation occurred in 12.5% of our combined study areas (as shown in Fig. 1), and ❖ www.esajournals.org the spatial extent and relative intensity of both motorized and non-motorized winter recreation varied notably within and across individual study areas (Fig. 2). In all study areas except the Tetons, motorized recreation represented the majority of the recreated footprint. The lowest overall levels of winter recreation occurred across much of our Sawtooth study area with <5% disturbance from each of motorized and non-motorized recreation activities though recreation did have areas of high localized intensity (Fig. 2). The highest overall winter recreation levels were in the southern Tetons where we recorded >50% of this area with winter recreation, primarily as non-motorized winter recreation (Fig. 2).
The spatial extent and relative intensity of backcountry winter recreation also varied within and across wolverine home ranges (Fig. 2). Motorized recreation, on average (AEstandard deviation [SD]), occurred in 22% AE 19% and 14% AE 15% of female and male home ranges, respectively, but varied greatly from <1% to 50%. Non-motorized winter recreation covered <5% of home ranges on average, and two females were not exposed to non-motorized recreation. The male monitored in the Teton study area had more non-motorized recreation than all other wolverines. Within home ranges, average recreation intensity of motorized recreation ranged from 0.0 to 42.2 m tracks/100 m 2 and average non-motorized recreation intensity value ranged from 0.1 to 9.3 m tracks/100 m 2 .

Potential habitat models: environment-only resource selection functions
The summed AIC score of the male and female potential habitat models (i.e., environment-only models) was notably lower than the AIC of the global model with ΔAIC of 1669, thereby justifying sex-specific models (Appendix S3: Table S1). The male model uniquely included covariates for distance to roads and the proportion of lower elevation grass and shrub land cover types. Alternatively, the female model included talus, persistent spring snow cover and forest edge:area covariates, which were not identified as important predictors of male habitat use. All covariates were Notes: SD, standard deviation. Home range areas were estimated using a local convex hull non-parametric kernel method (Getz et al. 2007).
† Animal-years indicates the total number of winter seasons cumulatively monitored accounting for multiple seasons of monitoring of some individual animals. Table 2. The number (%) of motorized and non-motorized recreation GPS tracks collected in our study areas, the annual average number of recreationists sampled (carrying or in a group with a GPS), the average annual trail use counts from infra-red trail use counters, and the estimated proportion of total use that we sampled (total people represented by GPS tracks/total use).  statistically significant. The models shared several covariates including topographic position index (TPI), quadratic form of slope, distance to forest edge, solar insulation and the percent cover of forest, riparian, and montane open cover types (Appendix S3: Table S1).
Cross-validation of female and male potential habitat models had similar Spearman rank correlations (r S ) of 0.92 and 0.91, respectively. Out-ofsample data validation similarly showed strong validation (female r S = 0.86, male r S = 0.95).

Realized habitat models: environment and winter recreation resource selection functions
Of the six models developed for male wolverines, Model 4 (combined recreation intensity) had the lowest ΔAIC (Table 3) and defined our realized habitat model for male wolverines (Appendix S4: Figs. S1-S4). There was a significant avoidance of areas with higher recreation Notes: AIC, Akaike's information criterion. Model 1 for male and female are the environment-only models. Models 2-6 use the environment covariates identified in Model 1 and winter recreation covariates to assess the responses of wolverines to different characteristics of winter recreation. Models 2-6 were developed separately for males and females. NA indicates not applicable.
† The realized models (Models 2-6) for males included recreated roads in the recreation covariates so the distance to road covariate in the Male Potential Model was redefined as distance to unrecreated roads in these models. intensity (b male = À0.06, SE = 0.01) though the overall importance of this was relatively low (ranked 9 out of 12 covariates) compared to other coefficients in Model 4 (Table 4). Ten-fold crossvalidation of this model showed high support for the model (r S = 0.91), and the out-of-sample male locations also validated very well (r S = 0.90).
The best-supported habitat model for female wolverines was Model 6 (Table 3; Appendix S4: Figs. S5-S8), with three significant (P-value < 0.01) winter recreation covariates: distance to linear recreation, intensity of dispersed motorized recreation, and intensity of dispersed nonmotorized recreation. Beta coefficients of Model 6 show females strongly avoided dispersed motorized winter recreation (b female = À0.31, SE = 0.02), and this covariate is the second ranked covariate (Table 4). Females also strongly avoided dispersed non-motorized winter recreation (b female = À0.19, SE = 0.01; ranked fifth in importance). Females avoided areas near recreated roads and groomed routes as indicated by the positive coefficient (b female = 0.08, SE = 0.01), and this covariate ranked 10 out 14. Similar to the male model, both the cross-validation and out-of-sample model validation showed strong support (r S = 0.91, r S = 0.83, respectively).
Model 6 did not provide the best overall predictor of male resource selection, but it allowed us to evaluate male wolverine responses to different forms of winter recreation (Table 4). All covariates in Model 6 were significant (or nearly so) in predicting male wolverine habitat use (Table 4). Similar to females, males avoided areas of dispersed motorized recreation (b male = À0.07, SE = 0.01), dispersed non-motorized recreation (b male = À0.15, SE = 0.02), and areas close to recreated roads and groomed routes (b male = 0.02, SE = 0.01) but the relative importance of winter recreation to males was less than for females. The importance of dispersed motorized recreation to male wolverine resource selection ranked 10 out of 13, while avoidance of dispersed non-motorized recreation was similar to females at a rank of 6. Avoidance of linear recreation by male wolverines was marginally insignificant (P = 0.056) and of lowest importance (Table 4).

Indirect habitat loss
Comparing the potential (Appendix S3: Table S1) and realized (Table 4) habitat models coefficients, there is very little evidence of confounding between the environmental covariates and the winter recreation covariates. Nine of the 12 environmental covariate coefficients in the female wolverine models were stable when comparing potential and realized models, including the top 7 ranked covariates (Table 4). Similarly, 9 of the 10 male model environmental coefficients were stable between models (Table 4).
Winter recreation resulted in indirect habitat loss of moderate and high-quality wolverine habitats as measured by areas transitioning to a lower class when comparing the realized habitat map to the potential habitat map (Fig. 3). On average (AESD), 14.1% AE 9.4% of female habitat and 10.9% AE 4.1% of male habitat was degraded to lower habitat classes, ranging from <10% to >70% within individual home ranges. This represented an average (AESD) area degraded by winter recreation within home ranges of 42 AE 36 km 2 for female wolverines (average home range 289 AE 92 km 2 ) and 118.2 AE 55.6 km 2 for males (average home range 1273 AE 471 km 2 ). Both the amount and severity of indirect habitat loss were related to the relative intensity of winter recreation within home ranges. The incremental effect of higher levels of winter recreation was large across home ranges with relatively low winter recreation levels (i.e., substantial habitat loss for each unit of recreation intensity), but the amount of indirect habitat loss tended to plateau across home ranges with the highest levels of recreation use (Fig. 4A). Female wolverines experienced more degradation to high-quality habitat, represented by a reduction in high-quality habitat to low-quality habitat (change of two classes; Fig. 4B). An average of 9.6% of available female high-quality habitat was degraded to low quality across the study area, while only 0.2% of available high-quality habitat for males was reduced to low quality.
These responses translated into more pronounced indirect habitat loss for females compared to males within the same landscapes. For example, a male and female that resided in the same landscape had similar average recreation intensity within their respective home ranges (Fig. 3). The female experienced indirect habitat losses of 36% and 38% of her high and moderate quality habitats, respectively, and 21% of the high-quality habitat was predicted to be degraded to low-quality habitat. In contrast, the male experienced predicted habitat degradation to 20% of high and moderate quality habitats, with only 0.9% of high-quality habitats predicted to be degraded to low-quality habitat.

Functional responses to winter recreation
Wolverines displayed negative functional responses in habitat use as the intensity of both motorized and non-motorized winter recreation increased. Use of areas with motorized recreation decreased as the average intensity of motorized recreation increased (Fig. 5A) within male and female home ranges, with slopes of 0.22 (R 2 = 0.40) and 0.38 (R 2 = 0.72), respectively. Similarly, both males and females showed negative functional responses to non-motorized winter recreation, even at the relatively low intensities of this recreation type. Habitat use of areas with non-motorized recreation declined as the availability of these areas increased within wolverine home ranges (Fig. 5B), with slopes significantly <1: 0.32 (R 2 = 0.89) and 0.10 (R 2 = 0.13) for males and females, respectively. The male functional response was driven by the high average intensity of non-motorized recreation that one male (2 animal-years) experienced in the Tetons. If the Teton animal was removed, male wolverines did not show a significant functional response to non-motorized winter recreation (Table 5). Additionally, the low R 2 of the female functional response to non-motorized recreation indicates high variation and a comparatively weak relationship.

DISCUSSION
We found that male and female wolverines showed some notable differences in the select-❖ www.esajournals.org 13 February 2019 ❖ Volume 10(2) ❖ Article e02611 ion for environmental covariates and that their selection for these covariates appeared to be independent of the potential effects of winter recreation. The realized habitat models that included winter recreation further showed that male and female wolverines responded negatively to increasing intensity of winter recreation within home ranges. Dispersed or off-road recreation activities elicited a stronger response than recreation along roads and groomed routes, with females showing more sensitivity to disturbance than males. The functional responses to dispersed recreation, particularly to motorized dispersed recreation, suggests that avoidance results in potentially important indirect habitat loss when a significant portion of an animal's home range receives recreation use, as it is exactly those animals exposed to higher levels of recreation that are most strongly displaced from these areas. Wolverines exposed to lower levels Fig. 4. The proportion of habitat degraded (A) and the proportion of habitat severely degraded (B) across home ranges of male (N = 12) and female (N = 13) wolverines (Gulo gulo) with varying levels of winter recreation intensity. Degradation is defined by the proportion of high and moderate quality habitat that degrades by at least one class (A; female R 2 = 0.93, male R 2 = 0.64), while severe degradation is measured by the proportion of the degradation that is high-quality habitat dropping two classes to low-quality habitat (B; female R 2 = 0.93, male R 2 = 0.44). Fig. 5. Functional responses of male and female wolverines (Gulo gulo) habitat use to the available relative intensity of (A) motorized (male N = 12, female N = 13) and (B) non-motorized (male N = 12, female N = 11) winter recreation in individual home ranges. The y-axis shows the average relative intensity of recreation at wolverine locations for each monitored wolverine, and x-axis shows the average recreation intensity within the animal home range. The dotted 1:1 slope line indicates the null expectation of random use. Responses below the 1:1 line indicate that use is lower than expected based on availability.
of winter recreation exhibit weaker avoidance based on the functional responses, which should result in relatively less indirect habitat loss. Also, the weak avoidance of areas near linear access used by winter recreationists suggests wolverines may be less sensitive to these linear disturbances.

Wolverine habitat selection
Previous habitat analyses in the Rocky Mountains for wolverines have been mainly at the firstor second-order of selection (Aubry et al. 2007, 2010, Fisher et al. 2013, Inman et al. 2013, identifying characteristics that predict the distribution or presence of wolverines. These efforts have indicated that wolverine are found at higher elevations , 2010, Krebs et al. 2007, Inman et al. 2013, in areas associated with late spring snowpack (Aubry et al. 2007, Copeland et al. 2010, Inman et al. 2013, and in alpine and subalpine habitats (Aubry et al. 2007) with higher topographic ruggedness (Krebs et al. 2007, Fisher et al. 2013, Inman et al. 2013 compared to the broader landscape. In contrast to the broader association to more rugged terrain, our analyses at the third order showed wolverines select less extreme topography characterized by concave or drainage bottoms (negative coefficient of TPI and slope covariates) within their home ranges. Additionally, our analyses showed selection for riparian habitats and forested edge habitats, which may represent good travel paths or more productive habitats (Scrafford et al. 2017) within a generally low productivity, high elevation landscape.
We expect that the habitat selection of our females was influenced by reproductive denning as 7 of 13 female animal-years represented denning females. Wolverine reproductive dens, particularly in the southern portion of their distribution, have been linked to deep and persistent snowpack and high structure such as talus boulders (Magoun and Copeland 1998). We found areas that support persistent spring snow as well as talus habitat were selected by female wolverines. In addition, females also selected for cold areas (negative solar radiation covariate), which also would support the selection for areas with persistent snow. Female habitat selection is complex, including characteristics that may be linked to some of the coldest and snowiest habitats as well as characteristics that may represent some of the more productive areas. Indeed, Krebs et al. (2007) proposed female selection was driven by a multitude of factors including food, predator, and human avoidance.

Influence of winter recreation on wolverine habitat
Wolverines maintained multi-year home ranges within landscapes that support winter recreation and some resident animals had >40% of their home range within the footprint of winter recreation. This suggests that at some scale wolverines tolerate winter recreation. However, within home ranges, wolverine avoided all forms of winter recreation and showed increasing avoidance of areas as the amount of off-road winter recreation increased, resulting in indirect habitat loss or degradation of moderate-or highquality habitats. Krebs et al. (2007) also found that wolverines, particularly females, avoided areas with winter recreation. Habitat displacement from winter recreation activities has been Table 5. Functional responses of male and female wolverines (Gulo gulo) to dispersed motorized (male N = 12, female N = 13) and non-motorized (male N = 12, female N = 11) winter recreation measured as the proportional use of recreation intensity compared to the average recreation intensity across home ranges of individual animals.
Model documented in other montane and alpine species. Endangered mountain caribou (R. tarandus caribou) in southern British Columbia have been displaced from high-quality winter habitat due to high levels of snowmobile recreation (Seip et al. 2007). In the Teton Mountains of Wyoming, backcountry ski recreation resulted in a 30% loss of high-quality winter habitat to bighorn sheep (Courtemanch 2014), and mountain goats avoided otherwise high-quality habitat associated with a developed ski area near Banff, Alberta (Richard and Cote 2016). Additional responses to winter recreation include changes in movement rates and temporal patterns, as was found in Canada lynx (Lynx canadensis) in response to winter recreation (Olson et al. 2018).
It can be challenging to identify animal responses to existing anthropogenic infrastructure and disturbance given the limited ability to control for these factors. One approach is to develop models that capture theoretical situations of no disturbance and compare these models to realized models that include the disturbance effect, which is the technique previous studies have used. For instance, Polfus et al. (2011) compared habitat models with and without human infrastructure covariates to assess indirect habitat loss to northern woodland caribou in northern British Columbia. Using a similar approach, Hebblewhite et al. (2014) modeled Amur tiger (Panthera tigris altaica) habitat with and without human-related covariates to evaluate anthropogenic habitat degradation. Patthey et al. (2008) used a regression approach to predict the potential abundance of capercaillie (Tetraago urogallus) if alpine ski recreation developments were not present, which they compared to the actual population estimate to assess the effects of winter recreation on the endangered Eurasian grouse.
Applying this approach to wolverines, we demonstrated that winter recreation had a stronger influence on female wolverine habitat selection than the habitat selection of males, as was also found by Krebs et al. (2007). Scrafford et al. (2018) also found that females are more sensitive than males to disturbances from industrial activities. Avoidance of areas with winter recreation degraded an average of 14% of moderate and high-quality female wolverine habitat, with 10% of high-quality habitat degraded two habitat classes to low quality. An average of 10% of male wolverine moderate-and high-quality habitat was degraded, and <1% of high-quality habitat degraded to low-quality habitat. While wolverine home ranges may be notably large, we expect female home ranges, in particular, represent the minimum spatial requirement necessary to provide the resources for the individual as well as offspring and kin as expressed by the resource dispersion hypothesis (Macdonald andJohnson 2015, Copeland et al. 2017). Rauset et al. (2015) found that wolverine reproductive success is related to habitat quality within their home ranges, suggesting factors that cause habitat degradation for reproductive females could translate into reduced fitness. A series of studies on mule deer responses to oil and gas development in Wyoming found avoidance of habitat surrounding oil and gas wells translated directly to declines in population size, empirically linking avoidance of habitat and fitness consequences (Sawyer et al. , 2017. We did not have the information required to assess demographic or fitness effects of winter recreation on wolverine. Our approach to estimate the indirect effects of recreation on habitat quality assumes independence between recreation and other environmental covariates. Our evaluation suggests minimal bias based on (1) our efforts to screen collinear, and hence, confounded variables in the development of RSF models, (2) the stability of the majority and most influential coefficients when comparing potential and realized models, and (3) 77% of our wolverine locations were outside the winter recreation footprint where confounding would not have affected the coefficient estimates for the potential model. Nevertheless, despite these precautions and caveats, our approach explicitly underestimates the potential effect of recreation on wolverines if recreation activities negatively influenced how wolverines used other environmental covariates.

Responses to recreation type
Male and female wolverine avoided both motorized and non-motorized winter recreation and avoided recreation occurring on and off roads. Females showed the strongest avoidance of off-road motorized winter recreation, which was the second most important predictor of female habitat use in areas where this recreation ❖ www.esajournals.org occurred, and they show a functional response of increasingly strong avoidance as exposure to dispersed motorized recreation increases within their home range. Dispersed or off-road motorized winter recreation also represented the largest proportion of the recreation footprint across our study areas, as well as occurring at much higher intensities than non-motorized recreation. These characteristics of dispersed motorized recreation and female response to it likely result in higher levels of indirect habitat loss experienced by females with higher levels of motorized recreation within their home range than our averaged population model indicates.
Both male and females also showed a strong avoidance of areas with dispersed non-motorized recreation, though these areas were limited within home ranges (<5% of home ranges affected on average). We recorded the highest and most extensive backcountry non-motorized recreation in the Teton study area, but we only captured one male wolverine in this study area. He exhibited strong avoidance of non-motorized recreation and was influential in our functional response analysis (Table 5). This suggests that the strength of avoidance exhibited by male wolverines to non-motorized recreation might depend on the intensity of recreation within their home ranges, similar to the functional response of female wolverines to dispersed non-motorized recreation. Given our limited sampling of male and female wolverines exposed to higher levels of backcountry non-motorized winter recreation, it would be useful to perform additional monitoring in areas with abundant backcountry, nonmotorized recreation.
Research examining wolverine responses to human infrastructure has suggested wolverines avoid roads, roaded areas, and development (May et al. 2006, Fisher et al. 2013, Inman et al. 2013, Stewart et al. 2016, Heim et al. 2017, Scrafford et al. 2018. Within home ranges and during winter when roads are covered in snow, we found human use of roads may be more important than the existence of the road itself in determining wolverine responses. Male wolverines were found closer than expected to unused roads but both male and female wolverines avoided areas near roads and groomed routes with winter recreation. Recent research in northern Canada also found that both males and female wolverines avoided active winter roads and that their movement rates increased with increased traffic volume (Scrafford et al. 2017(Scrafford et al. , 2018. In our research, the avoidance of recreated roads was significant but relatively weak compared to avoidance of off-road recreation areas, suggesting that spatially predictable or confined recreation travel patterns may be perceived by wolverines as less risky. Harris et al. (2014) also reported less disturbance to northern ungulates from road-based recreation as compared to recreation that is unpredictable in space or time.

Cumulative impacts of climate change and winter recreation
Both wolverines and backcountry winter recreation are expected to be affected by climate change, potentially resulting in a funnel effect where the overlap between winter recreation and wolverine distribution increases as they both respond to declining snow extent, depth, and the snow season. In the southern portion of their North American range, wolverines appear to be tightly linked to the area defined by the presence of persistent spring snow (Aubry et al. 2007, Copeland et al. 2010, Inman et al. 2013. The underlying ecological requirements that drive this close relationship may include denning requirements Copeland 1998, Copeland et al. 2010), a dependence on scavenging large ungulate carcasses effectively preserved within and under the snowpack (Mattisson et al. 2016), caching food (Inman et al. 2012a), and competitor or predator avoidance (Mattisson et al. 2016). Heim et al. (2017) suggested that the association of wolverines to persistent spring snow makes them vulnerable to climate changes, and McKelvey et al. (2011) predicted a 67% loss of wolverine habitat in the western United States by 2059 due to loss of snowpack.
The demonstrated loss of snow pack and reduced length of winter (Mote et al. 2005) may also have profound impacts for winter recreation in the future (Bowker et al. 2012, White et al. 2016, Wobus et al. 2017. While the reductions in winter length are predicted to cause a decline in per capita participation in winter recreation, human population growth may counter these declines and most projections of winter recreation are stable or increasing (Bowker et al. 2012, White et al. 2016, Wobus et al. 2017. Winter recreationists will likely need to adapt when and where they recreate to adjust to shortened snow season and reduction of winter recreation areas due to snow loss (Dawson et al. 2013, Rutty et al. 2015. Winter recreation may become more concentrated and intense in both space and time (Dawson et al. 2013, Rutty et al. 2015, especially during the mid-to late winter period when snowpack is predicted to be the most consistent (Mote et al. 2005). This is also the time period when female wolverines are entering reproductive dens. Predictions of winter recreation distribution and intensity would likely suggest even more severe indirect habitat loss than our current assessment indicates. Our results underscore the importance of managers to consider growth of the recreation industry concurrent with declining habitat for winter recreation, which could exacerbate conflicts between recreation and wildlife.

CONCLUSION
Balancing the many positive benefits of outdoor recreation with the impacts it may have on natural systems is a growing field of study. Our research into the effects of backcountry winter recreation on wolverines represents information at spatial and temporal scales rarely achieved in other disturbance research. Habitat quality has been linked to reproductive success in wolverines (Rauset et al. 2015), and sufficiently high levels of indirect loss of high-quality habitats through disturbance would affect the reproduction and survival of animals. However, thus far we do not have the information to assess the population level effects of winter recreation on wolverines. Here, we have shown significant avoidance by wolverines of areas used by backcountry winter recreationists and that this results in habitat degradation, particularly for female wolverines. Given the low density and fragmented nature of wolverines in the contiguous United States, impacts to the relatively few reproductive females should be of concern.
Our results suggest that winter recreation should be considered when assessing wolverine habitat suitability, cumulative effects, and conservation. We found that the effects of winter recreation on wolverine habitat are dependent upon the intensity of recreation and that winter recreation patterns are highly variable at the scale of wolverine home ranges such that some animals may experience higher levels of indirect habitat loss while adjacent animals may experience little. Our research provides land managers with a more detailed understanding of important habitat characteristics used by wolverines and should inform management of wolverine habitats across the extensive landscapes they use. These backcountry landscapes represent critical habitats for wolverines, important and highly valued areas for people to connect with nature, and are economic drivers for the small communities that surround them. Solutions to finding a balanced approach to sustaining the diverse values of these wild landscapes require creative approaches and collaboration between land managers, stakeholders, and wildlife professionals.