Winter hunting behavior and habitat selection of wolves in a low-density prey system

The functional response is the relationship between food intake rates and prey density, and is shaped by factors including handling time, predator speed, habitat or prey movement. For many predator—prey systems, the density-dependent functional response is represented by a type II or type III functional response. Determination of the relationship type is important, as managers can often predict the response of predators to changing prey densities. In wolf—moose (Canis lupus—Alces alces) systems with relatively high prey density, the functional response often follows a predicted type II functional response. However, in a very low prey-density system, wolves have previously been shown to escape the density dependent phase of the functional response and demonstrated kill rates mimicking high prey-density systems. We conducted a study to evaluate winter wolf movements between moose kills in the Yukon Flats, Alaska where moose exist at densities <0.2 km-2. Our research objectives were to understand whether habitat selection when moving and specific behaviors could be mechanisms used by wolves to maintain kill rates that mimic those in high prey density systems and if those behaviors may allow wolves in our study system to escape a density dependent functional response. We used GPS collars to characterize wolf travel paths between kills to estimate wolf travel speed, movement distance, time between kills, and handling time of each kill. Our results demonstrated selection for frozen river corridors by wolves and provided new information on long-distance movements in a low prey-density system. These adaptations may influence the functional response by moderating the effect of low prey densities.

Within predator-prey systems the relationships between the consumption rate of prey by predators and the prey density is known as the functional response and are generically represented as type I, II or III functional response curves (Holling 1959, Solomon 1949. For many predator-prey systems the density-dependent functional response is represented by a type II or type III functional response, but it is possible for predator-prey relationships to fall between type II or type III at different life stages (Streams 1994). Traditionally type III curves are associated with large-bodied predators and prey (Hassel et al. 1977). Functional response is primarily moderated by prey kill rate and handling time (Hassel 1978), but a variety of other factors may be influential such as prey size, attack rate, encounter rate, multiple predators, learning, adaptation or prey-switching (Abrams 1990, McCoy et al. 2012, Streams 1994, Van Leeuwen et al. 2013. Kill rate, and hence the functional response, is also influenced by both predator and prey movement rates (Streams 1994). Predators select for specific habitat characteristics, high prey density, prey age or prey vulnerability to increase kill rates or efficiency (Hebblewhite et al. 2005, Sand et al. 2012, Montgomery et al. 2014, Kittle et al. 2015. Using a predator-prey system that consists of wolves Canis lupus relying on a low-density population of a single prey species, moose Alces alces, we explored how predator movement may shape the functional response. More specifically, our research provided new insight on how the adaptive capacity of a predator's movement behavior and habitat selection may uncouple the relationship between kill rate and prey density. Wolves are coursing predators that hunt continuously while on the move and simultaneously maintain territory boundaries (Mech 1970, Mech andBoitani 2010). While hunting, they may modify their speed, travel distances, amount of area searched, prey selection, terrain or habitat selection to maintain kill rates or increase prey encounters (DeCesare 2012, McPhee et al. 2012a, Mech and Cluff 2011, McKenzie et al. 2012, Moffatt 2012, Vander Vennen 2016 Mech et al. 2015). Wolf selection for linear corridors results in low density prey being at higher risk of predation in environments with densities of linear corridors (McKenzie et al. 2012).
In wolf-moose systems, a type II functional response has been predicted and observed (Hayes and Harestad 2000, Messier 1994, Zimmermann et al. 2015. However, Lake et al. (2013) found in their system of very low moose densities that a type II functional response was not observed because there was no density-dependent response. They report that wolves maintained a kill rate in a low prey-density system (0.2 moose km -2 ) comparable to high prey-density systems. In that system, wolves may have adjusted pack size to accommodate for low prey density. However, the functional response was likely influenced by other factors, and Lake et al. (2013) speculated that wolves were selecting corridors to facilitate travel or changing their movement characteristics. We extend on that study by analyzing movement behavior of predators in a low-density ungulate system. We analyzed travel paths when hunting, speed, distances traveled and underlying habitat characteristics of wolves from six packs in the Yukon Flats of Interior Alaska (Fig. 1), and compared those characteristics to movements of wolves in systems of higher prey density reported in the literature. Our research objective was to investigate underlying drivers of the functional response by analyzing the movement behavior and modeling habitat selection of wolves while traveling (e.g. excluded resting and kill site behavior) in a low-density prey system. We hypothesized that wolves in our system were traveling farther than wolves in high prey-density systems to make kills and were maintaining high kill rates by utilizing landscape characteristics that aid efficient travel, such as non-forested areas or river corridors. Such use of river corridors may also affect the functional response as greater numbers of prey may be encountered because moose may preferentially forage in these areas in winter (Baigas et al. 2010, MacCracken et al. 1997, McKenzie et al. 2012.

Study area
We conducted our study in the western Yukon Flats of Interior Alaska (Fig. 1). The Yukon Flats is bounded by the Brooks Range to the north and the White Mountains to the south. Elevations within our study area range from 91 to 912 m, but most of the area is low and flat. The Yukon River bisects the region and at its center is the confluence of the Yukon, Porcupine and Chandalar rivers. The Yukon Flats National Wildlife Refuge (Yukon Flats NWR) covers approximately 34 000 km 2 (8.6 million acres) and a majority of our study area. It stretches approximately 350 km from east to west and 190 km from north to south. Based on the 2001 National Land Cover Dataset, the Yukon Flats is 67% boreal forest and 33% riparian areas. Boreal forests and riparian species include white spruce Picea glauca and black spruce P. mariana, white birch Betula papyifera, aspen Populus tremuloides and poplar P. balsamifera, alder Alnus spp. and willow Salix spp. (Homer et al. 2007).
The climate of the Yukon Flats is classified as sub-arctic and characterized by long cold winters (November-March) and short dry summers (May-August). Temperatures are seasonably variable, and the mean temperature is -28.5°C in January and 16.7°C in July. The dry climate generates snow depths much less than 90 cm, which is considered a threshold that results in changes in moose movement and survival (Coady 1974, Gasaway et al. 1983, 1992. During our study period, snow depths at two snow stations averaged 69 and 48 cm. The 10-year average at those stations was 52 and 64 cm, respectively (Natural Resources Conservation Service 2015).
Aerial estimates during the study period indicated  0.2 moose km -2 in the western and eastern Yukon Flats (Lake et al. 2013). Wolf densities in the Yukon Flats were estimated at 3.4-3.6 1000 km -2 (Lake et al. 2015). Moose densities are thought to remain at a low-density equilibrium due to high calf mortality from bears, Ursus americanus and U. arctos, and adult mortality from wolves, combined with illegal harvest of adult females (Bertram andVivion 2002, Gasaway et al. 1992). Within the Yukon Flats, moose are the primary food source for wolves, with occasional takes of snowshoe hare Lepus americanus or beaver Castor canadensis (Lake et al. 2013). Caribou Rangifer tarandus are not common in the area.

Data collection
Wolves were chemically immobilized by darting from a helicopter (US Fish and Wildlife Service Region 7 Animal Care Protocol no. 2008022), beginning in November 2009 in the region of Beaver, Alaska (Fig. 1). Further details of wolf immobilization are described in Lake et al. (2013). Nine wolves from six packs were marked with Telonics model TGW-3580 GPS radio collars. The GPS collars recorded locations at three-hour intervals and had a life expectancy until May 2010. All data were accessed from the collar following recapture in April 2010.
Moose kill site locations were determined by Lake et al. (2013), using aerial surveys coupled with an analysis of location clusters. Webb et al. (2008) reported that a four-hour GPS interval was sufficient to identify 100% of kill sites by wolves on large-bodied prey, such as moose. Lake et al. (2013) used three-hour intervals and reported no errors related to incorrectly classifying a kill as a non-kill. Hence, it was unlikely that they omitted any kills. At the conclusion of their study and based on their cluster modeling, they identified thirteen location clusters (19% of confirmed kills) of seven locations or more where flights did not confirm if a kill existed; if the location clusters were classified as a kill, but were instead a rest site, a commission error (i.e. classifying a rest cluster as a kill cluster) may occur. The six packs were monitored for different amounts of time.  (Table 2). All GPS collars demonstrated a high fix success (mean = 98%, range = 96-99%) rate, which was attributed to flat terrain and lack of canopy (Lake et al. 2013).

Dataset preparation
In packs with two collared individuals, we observed that the collared individuals traveled within 25, 50, 75, 100 and 200 m of each other 66, 80, 83, 86 and 90% of the time respectively. Since they traveled together (50 m) a majority of the time, we chose one wolf from each pack to represent all movements of the pack. We further justify that decision based on high likelihood of a carcass being attended by both wolves in packs with two collared individuals during the winter (Metz et al. 2011). Second, no pack maintained two operating collars for the entire winter due to mortalities, collar slippage or collar failure. We used locations from the breeding male or female from each pack except Hodzana Mouth, where locations from a juvenile were used because the breeding female collar failed prematurely. In the final dataset, we standardized GPS data for each pack by removing points from capture up to their first kill and after the date of their last kill to collar retrieval. In the analysis, we included kills (n = 68) that were confirmed through aerial observation and location clusters that lasted longer than one day (n = 10 fixes) where the model of Lake et al. (2013) predicted the cluster was a kill. Errors of omission were zero for the model of Lake et al. (2013), but we acknowledge a commission error could have occurred. Such a commission error would have resulted in rest sites mistakenly being classified as kills. This would have decreased true search distance or decreased time to kill. For each individual, we characterized the resulting GPS location data into four distinct behavioral classes that have been used in previous studies to characterize wolf movement (DeCesare 2012, McPhee et al. 2012a, hereafter, referred to as 'path characterization'. They included presence at kill site, resting, kill-site revisits, and traveling. Lake et al. (2013) located kill sites through aerial surveys, clustering and tracking. We characterized the first kill-site point as the first time that a wolf arrived at a kill location and all locations were considered to be associated with the kill site until the wolf left for more than 24 h (eight locations). Once the wolf left for more than 24 h, we used the location closest to the kill site as the last kill-site location. We characterized rest locations as any time two-or-more consecutive locations within travel paths did not change more than 26 m from the last location (i.e. in 3 h). This distance was the approximate maximum accuracy of the GPS location (Adams et al. 2013), hence any locations that did not move more than that could be considered the same location. Revisits included all locations where a wolf returned to a kill and remained there 6 h (i.e. two fixes) or more. Traveling included paths between kill sites, but excluded rest locations and all but the first revisit location at the kill. We maintained the first revisit location to keep the travel path intact.

Data analysis
We derived several descriptive statistics from our path characteristics that could provide insight to the functional response. We chose these parameters because they are quantifiable, comparable to previous literature, and hypothesized to be behaviors that wolves can modify to adapt to a low preydensity system. These statistics included mean and standard deviation of handling time (days), median and maximum days spent traveling (time between kills), median and maximum distance (km) traveled, and median and maximum travel speed (km h -1 ). Days spent traveling, travel distance, and travel speed were strongly skewed right and reporting their mean would be inappropriate. Handling time is the amount of time a pack of wolves requires to consume an animal after it is killed. Within our dataset, we determined handling time to be the interval between when the kill began and the first time the wolf left the kill for more than 24 h. We calculated the travel speed by dividing the segment length (i.e. distance between two GPS locations) by the total time elapsed during that segment. A log 10 transformation was used to normalize the distribution of handling time, time between kills and travel distances. To examine pack differences, we used an analysis of covariance (ANCOVA) to test for differences while controlling for pack size. If a difference was detected by the ANCOVA, we used a t-test with a Bonferroni adjustment to determine which packs differed. All models were checked to ensure that they met basic assumptions of normality and homogeneity of variance.
We examined underlying habitat selection during travel to aid with inference of our movement statistics output. We hypothesized that wolves were utilizing corridors such as rivers or habitats with minimal travel barriers to enable efficient travel (i.e. travel with the least amount of energy expended), and our covariates were chosen to test corridor usage. We defined a corridor as a landscape or habitat feature that enhances the movement of an animal (Bennett 1999). Previous studies have associated wolves with linear corridors that may increase speed up to 2.8 times over forested habitats (James 2000). In the winter, rivers become frozen and hard packed with snow, reducing energetic expenditure during travel. For instance, river corridors used intensively by wolves may be avoided by prey as an anti-predator strategy (Bergerud and Page 1987).
We assessed habitat selection while traveling using a step selection function (SSF). SSFs are an effective method that use movements of animals during discrete time steps to quantify fine-scale selection patterns (Thurfjell et al. 2014). Matched sets of used and available steps are compared using conditional logistic regression, taking the same generalized exponential form as a resource selection function with a log-link function (Fortin et al. 2005). Five available steps were generated for each used location by randomly drawing step length and turn angles from two distributions established from observations of monitored individuals. Steps can be characterized by the line segments between locations, the average continuous habitat variables along the step, the proportion of habitat along each step, or by the environmental characteristics at the endpoint of each step. We used the endpoints of each step for our analysis because selection of linear features (e.g. travel corridors) can be underestimated when landscape variables are measured along the lines between steps (Thurfjell et al. 2014). Five available steps were generated for each used location by randomly drawing step length and turn angles from two distributions established from observations of monitored individuals.
In order to maintain statistical power, we only chose biologically plausible covariates that could be related to wolf travel paths (Table 1). We measured distance to rivers and waterbodies as the distance from a location to the nearest river or waterbody of the high-resolution National Hydrography Dataset at a scale of 1:24 000 (United States Geological Survey 2015). We measured distance to ridges as the distance of locations to ridges derived from a 17-m resolution digital elevation model. We used the National Landcover Dataset (NLCD) from 2001 to generate underlying categorical habitat variables. We grouped NLCD into four broader categories based on habitat height. NLCD 1 was the water NLCD class, which includes water bodies or rivers greater than 30  30 m in width or area. NLCD 2 included shrub land cover of medium height. NLCD 3 included tall tree classes, and NLCD 4 included riparian or wetland classes with short or grassy vegetation. We assumed that some NLCD categories created efficient travel corridors in open habitats (e.g. NLCD 1 , NLCD 4 ) and some created barriers to travel in tall or medium vegetation-height habitats (e.g. NLCD 3 , NLCD 4 ) (James 2000).
We used a two-stage modeling approach that fits models separately for each individual animal and then averages regression parameters across individuals to quantify Table 1. Summary of covariates utilized in step selection function (SSF)analysis of wolf Canis lupus habitat selection during winter in the Yukon Flats, Alaska. Categorical variables were grouped together as: NLCD 1 (11 -water), NLCD 2 (31 -barren, 52 -shrub scrub, 51 -dwarf scrub), NLCD 3 (41 -deciduous forest, 42 -evergreen forest, 43 -mixed forest), NLCD 4 (72 -sedge/herbaceous wetlands, 90 -wood wetlands, 95 -emergent wetlands). The group description describes the continuous and categorical variables, and GIS layer derived from describes the data surface or derived data surface. population-level patterns for wolf packs (Fieberg et al. 2010). We fit conditional logistic regression models for each individual wolf with matched sets of used and available locations using the coxph package in R ver. 3.2.0 ( www.rproject.org ). We then averaged logistic coefficients and standard errors across individual wolf packs as an estimate of the population-level effect of predictor variables on the relative probability of use (Sawyer et al. 2009). To normalize the distance to water and ridge variables, we used a log 10 transformation. Prior to modeling, we tested for multicollinearity bases on Pearson's pairwise correlation analyses, and did not find any highly correlated variables (|r| 0.70). We did not use an information theoretic approach for model selection because these methods lack standardized approaches to keep the animal as the experimental unit and build a population-level model from a common set of predictor variables (Sawyer et al. 2009). Instead, we used a t-statistic to test whether coefficients averaged across individuals were significantly different from zero (α  0.05), and included only those significant variables in the population-level model (Hosmer andLemeshow 2000, Squire et al. 2013).

Results
Our analysis dataset contained 5561 locations from six packs, and the number of locations from each pack ranged from 499-1123 (Table 2). All kills were of moose. We analyzed 68 unique paths to kills during our study period. The number of kill paths was less than the number of kills as some paths led to multiple kills at one site. The number of kills varied by pack, and the largest packs made the most kills ( Table 2). The maximum handling time was 16 days and mean handling time was 4.0 days (SD = 2.5) for all packs (Table 2). While controlling for pack size, the logtransformed handling time of Beaver Creek was significantly different from Hodzana Mouth (ANCOVA, p = 0.03, df = 5, F = 2.73, t-test, p = 0.04). Average days between kills (ANCOVA, p = 0.17) and travel distance between kills (ANCOVA, p = 0.39) were not significantly different among packs. The maximum number of days between kills was 17.8, the median was 5.6, and the mean days between kills was 5.9. The maximum distance between kills was 263 km and median travel distance was 53 km. Log-transformed travel distances and times between kills were highly correlated (p  0.01, r 2 = 0.85, df = 56, F = 340.9). Median travel speed was 0.4 km h -1 and mean travel speed was 0.6 km h -1 . From 2009 to 2010, we modeled how traveling wolves selected landscape characteristics during winter. SSF coefficients averaged across individuals that were significantly different from zero (α  0.05) included effects of distance to water, and NLCD 2 and NLCD 3 (Table 3). While traveling, wolf habitat selection during winter was characterized by shrub (NLCD 2 ) and forest habitat (NLCD 3 ) types, and a decreased distance to water. Wolves did not exhibit selection for open water (NLCD 1 , p = 0.18) and wetland (NLCD 4 , p = 0.13) habitats. Additionally, wolves did not show strong selection for distance to ridgelines (p = 0.47).

Discussion
Our results supported our hypothesis that wolves were traveling further to make kills and were selecting for river corridors. This may demonstrate how habitat selection and movement distance of wolves in a low prey-density system can provide insight into how plastic behavioral response of a predator may theoretically moderate the functional response.
The results of the SSF suggest that wolves were selecting for frozen river corridors. Our results were consistent with other studies that report selection for linear corridors, including rivers or seismic lines (James and Stuart-Smith 2000, McKenzie et al. 2012, McPhee et al. 2012a. McKenzie et al. (2012) found wolves used seismic lines for travel, and surmised that the functional response in low prey-density systems was strongly influenced by prey clustering around seismic lines, as this would increase the encounter rate. The Yukon Flats is laden with an intricate network of streams that create a high density of linear corridors. Selection for shrub and forest habitats did not align with our hypothesis that wolves would avoid barriers to travel. Wolves select for regions of higher prey density, and may select shrub or forest habitat along riparian corridors because moose utilize the corridors and adjacent habitat in the winter (Baigas et al. 2010, MacCracken et al. 1997, McPhee et al. 2012b). However, selection for shrub habitat should be interpreted cautiously as we ran a parallel analysis using a Resource selection function (RSF) and found that nearly all variables significant to the SSF were also significant to the RSF, except that shrubs were selected against in the RSF (Supplementary  material Appendix 1 Table A1-A2). We speculate that the RSF may have selected against shrubs because that approach did not account for travel paths between point locations. Along braided river corridors, buffered travel paths may include more sandbars and recently disturbed areas (e.g. seasonal flooding) dominated by shrubs. Selection for linear corridors also may facilitate longer travel distances between kills in our low prey-density system. Moffatt (2012) reported distance and time to kill in a multiprey system, characterized by moose and caribou (Table 4). In his study system, moose density ranged from 0.12-0.25 moose km -2 . While an estimate of woodland caribou densities could not be found for his study area, total prey densities were probably higher than the density of moose (0.2 moose km -2 ) in the Yukon Flats. Our result of 53 km between kills was 1.9 times greater than the search distance reported by Moffatt (2012) of 27 km. Average time to kill reported in our analysis (5.9 days) was 1.4 times greater than in Moffatt (2012). Although the GPS fix interval in Moffatt (2012) was 5 h and the fix interval in this study was 3 h, we can gain confidence that wolves in the Yukon Flats were traveling farther to kill prey than reported by Moffatt (2012) because of the longer travel time. McPhee et al. (2012b) reported that average time to event in their highdensity system was 5.3 days, which is similar to time to event in our low-density system (5.9 days). Sand et al. (2005) reported in a very high prey-dense system of moose that time to kill was approximately 4 days. Our results also corroborate with Lake et al. (2015) who documented large wolf territories in the Yukon Flats, and we speculate that long travel distances are necessary to encounter sufficient vulnerable prey (Mech et al. 2015). Compared to low preydensity systems, we speculate that wolves travel less distance in high prey-density systems because vulnerable prey are encountered at greater rates. During long searches, wolves may be either not finding prey at all, or they may not be encountering vulnerable prey.
Long time-to-kill intervals may require that wolves optimize kills by consuming everything possible, thus adding to long handling times. Handling time was longer in our study system than previously reported by several studies with largebodied prey (Eriksson 2003, Sand et al. 2005, although shorter than the results of Messier and Crête (1985) (Table 4). Pack size did not explain handling time in our study. We were unable to differentiate between adult and calf moose kill sites, but prey size did not seem to be a tenable explanation, as Hayes et al. (2000) found that handling time for adult or calf moose (2.9 and 2.6 days, respectively) was similar. Longer handling time may be a function of prey density and we speculate long handling times reflect the effort required to secure a kill in a low prey-density system. Long handling times may indicate that wolves were food-limited and they completely consumed the kill before initiating another hunt. Investigations on the ground at a sample of kills by Lake et al. (2013) found that kills were completely consumed before wolves revisited. This is in contrast to Eriksson (2003) who documented partially consumed moose in their high prey-density system. On Isle Royale, Michigan, partial moose consumption may be an optimal foraging strategy for wolves and associated with severe winters (Vucetich et al. 2012). Partial moose consumption in Scandinavia may be linked to high disturbance rates of wolves at kills by humans (Sand et al. 2005). In contrast, disturbance of wolves by humans in the Yukon Flats is low and unlikely to be a factor.
The inverse relationship of prey abundance and wolf territory size is well documented in the literature (Messier 1985, Fuller et al. 2003, Kittle et al. 2015. However, to our knowledge an underlying mechanism has not been well described. Our finding (i.e. long travel distance to make a kill) presents a mechanism to explain why territories are so large; as wolves travel longer distances in response to low availability of vulnerable prey. Wolves in the Yukon Flats exhibit some of the largest territories and lowest densities reported in North America (Lake et al. 2015). If prey density or vulnerability were to further decline, we hypothesize that wolves may further adjust travel distances to locate vulnerable prey. If the distance and time traveled increases beyond what wolves can survive, then wolves may starve, reduce pack size (i.e. small packs) or disperse from the territory to find suitable prey availability. Thus, wolf persistence in our study area appears to be limited by the energetic constraints This study # Kill rate = a study designed to specifically determine kill rate, Movement = a study designed to look at movement parameters, Time to kill = a study designed to look at event timing.
associated with locating a sufficient number of vulnerable moose in this low prey-density system.