We used the locations of moose killed by wolves and a spatial index of moose density to calculate a relative predation rate in a spatially heterogeneous and highly disturbed landscape. We tested whether wolves kill a spatially variable proportion of moose by estimating the relationship between relative moose density and the frequency of kills across gradients of natural and anthropogenic landscape features. By including both natural and anthropogenic features, we were able to compare the strength of each source of spatial heterogeneity on predation rate dynamics. Human‐induced rapid ecological change can disrupt predator–prey spatial interactions by adjusting either the type or behavior of predator (Latham et al. 2011 , Sih et al. 2011 ). We therefore predicted that because human disturbance is a novel source of heterogeneity, disruption to the wolf aggregative response would lead to stronger changes in predation rate near oil sands mining features than near or in natural habitats.

Alberta's Athabasca oil sands region (AOSR) is a region of boreal forest in the Canadian western sedimentary basin with extensive deposits of bitumen. The AOSR is characterized by extensive human disturbance including open pit mines, tailings ponds, and industrial facilities (Fig. 1 ; Schindler 2010 ), with a total footprint of approximately 650 km 2 . Mining activity in the AOSR has removed large swaths of habitat with considerable variation in human density and intensity of use across the footprint. The AOSR is home to spatially overlapping moose ( Alces alces ) and wolf ( Canis lupus ) populations (Fuller and Keith 1980 , Wasser et al. 2011 ). Across the boreal forest and mountainous regions of North America, wolves respond to human disturbance depending on the types and intensity/frequency of use of the disturbance. When the intensity of use by humans is low, wolves use human disturbance such as linear features but reduce their use with increasing human use (Whittington et al. 2005 , Hebblewhite and Merrill 2008 , Houle et al. 2010 , Rogala et al. 2011 ). Several studies have shown that moose use areas avoided by their predators near human disturbance to reduce the probability of predation (Edwards 1983 , Stephens and Peterson 1984 , Dussault et al. 2005 , Berger 2007 , Latombe et al. 2014 ).

Both predators and prey can respond to spatial heterogeneity caused by novel human disturbance such that the proportion of prey killed by predators is a function of proximity to human disturbance. However, the direction and magnitude of such effects are not always predictable. Increased prey density at human‐made forested edges can elicit a predator functional response leading to increased predation rate (Gates et al. 1978 ). Predators can aggregate in areas disturbed by humans such as along linear features to increase movement rates (Dickie et al. 2017 ), potentially increasing predation rates in the absence of prey density changes (James and Stuart‐Smith 2000 , Decesare 2012 ). When a predator avoids human disturbance (Frid and Dill 2002 ) more than their prey, predation rates will be decreased near disturbance due to a prey refuge effect (Hebblewhite et al. 2005 , Berger 2007 ). When predation influences prey population dynamics, it is important to assess if and how predation rates vary with novel human disturbance.

Predation rate, or the proportion of prey killed by a predator, is a fundamental component of understanding the effect of predation on prey population dynamics (Holling 1959 , Messier 1994 ). Predation rate varies when predator abundance (numerical response), per predator kill rate (functional response) or both, varies with prey density (Solomon 1949 , Holling 1959 , Vucetich et al. 2011 ). Predator abundance changes temporally in response to prey density via predator reproduction and spatially due to a predator's aggregative response, by which individual predators hunt more in some areas than others in response to prey clumping (Holling 1959 , Hassell 1978 ). Therefore, landscape heterogeneity causes predation rate to vary spatially when prey aggregate in habitat with increased access to food or mates and predators allocate disproportionately more hunting effort there than other areas.

We calculated a study area as a minimum convex polygon calculated for the 95% isopleth of summed kdes for each wolf pack. This generated an 8122 km 2 study area polygon used for modeling. We compared the locations of kills to random locations (10 km −2 ) in the study area using mixed‐effects logistic regression with a random intercept for wolf pack to account for varying GPS data collection durations. We fit multiple models containing all combinations of five variables (Table 2 ) interacted with moose density such that interaction coefficients estimated how predation rate changes with each natural and anthropogenic variable. For our inference, we selected the top model using AIC. No covariates were collinear >0.7, and we scaled all continuous variables so they were centered on zero. For ease of interpretation, we re‐calculated the distance to disturbance variables as their difference from the max distance such that the largest values were closest to the disturbance. We compared the odds ratios of interaction terms (predation rate) in the top model. All modeling was conducted in R 3.3.0 (R Core Team 2016 ). Code used to summarize the clusters in the GPS data is included in the supplementary materials (Appendix S2 ).

The distribution of linear features was delineated by the Alberta Biodiversity Monitoring Institute at a 1:15,000 scale using 2012 SPOT imagery. We used proximity to the city of Ft. McMurray, located in the south of our study area, as a surrogate for intensity of human recreational use on linear features. Using network analyst in ArcMap 10.1 (ESRI 2011 ), we calculated the density of seismic, transmission, and pipelines as km/km 2 weighted with network distance to Ft. McMurray. We re‐calculated the network distances as their difference from the max distance such that the largest values were closest to the Ft. McMurray.

We defined the mine footprint as mining excavation including pit mines and tailings ponds, and the facilities footprint as buildings, oil sands upgraders, processing plants, work camps, and parking lots. To delineate the borders of mines and facilities in AOSR, we modified 2009 landuse shapefiles supplied by the Regional Aquatics Monitoring Program ( 2011 ) using 2012 SPOT satellite imagery. We used only facilities that were larger than one km 2 for analysis and clipped facility polygons from mines. We removed locations falling inside mines and facilities and calculated the distances to mines and facilities for all locations outside mines.

We calculated the kernel density of moose sightings using kernel density in ArcGIS 10.1 (ESRI, Redlands, California, USA). The kernel density estimator was weighted by the number of moose sighted at each location, using a pixel size of 1 km 2 and a search radius of 8 km, corresponding to the average moose home range diameter (7.15 ± 0.55 SE km) in our study area. We rounded the search radio to 8 km to be more inclusive of the moose in our study area and to create a smoother density surface. Moose home range areas were calculated using moose GPS telemetry from collared moose. We collared moose throughout the study area in February 2010 (University of Alberta, ACUC Study ID AUP00000102). Twenty‐five moose were captured and outfitted with a GPS (Lotek 7000MU) collar (University of Alberta, ACUC Study Id. AUP00000102). Moose GPS collars were programmed to fix the moose's location every three hours. We calculated each moose home range area from the 95% isopleth of a kernel density estimate (kde) using the least‐squares cross‐validation (Seaman et al. 1999 ) smoothing factor in Geospatial Modelling Environment (Beyer 2012 ).

We calculated a spatial index of moose density using Alberta Environment and Parks moose surveys of our study area. Alberta Environment and Parks conducts ungulate surveys using a random stratified block survey in which they fly initial surveys along lines of latitude (~1.8 km spacing) to stratify the area by moose density. We used the sightings from these initial stratification flights in order to capture uniform survey effort over our entire study area. Stratification flights are flown with a Cessna 206 fixed‐wing aircraft at 100 km/h allowing sightings of moose 250–300 m on either side of the aircraft (Morgan and Powell 2010 ). The four wildlife management units (WMU) in our study area were flown in the winters of 2013, (WMU 518), 2010 (WMU 530), 2009 (WMU 531), and 2008 (WMU 519).

We assessed the ability of the selected kill and species model to distinguish among kills and beds, and species, using k ‐fold cross‐validation (Boyce et al. 2002 ), for which we re‐fit the model to a 60% subset of the data, then calculated predicted values from the remaining 40%. From the predictions, we calculated the area under the receiver‐operating characteristic (ROC) curve (AUC) and optimal threshold to maximize both the model classifier specificity and sensitivity. We bootstrapped the above evaluation 1000 times. Using the mean optimal cutoff for both the kill and species models, we predicted first whether each GPS cluster from the full dataset of clusters summarized between October and March, 2012–2014 was a kill and then among the kills, the prey species. Once GPS clusters were predicted for each wolf, we combined GPS cluster centroids from all wolves within the 300 m and within a max of 30 d by retaining only the first GPS cluster of the first wolf at the GPS cluster.

We plotted the location of these ground‐truthed (GT) GPS clusters and selected from the full set of summarized GPS clusters the closest centroid in space and time that began before our GT visit. When two wolves were present at the same GPS cluster, we selected the summarized GPS cluster from the first individual present to maintain a single wolf‐kill sampling unit for modeling. Using this subset of 49 summarized GPS clusters, we modeled the effect of GPS cluster characteristics on whether the GPS clusters were beds or kills (kill model), and among kill clusters whether the prey was deer or moose (species model) using a two‐step model selection (Webb et al. 2008 , Knopff et al. 2009 ). We coded GPS clusters known to be beds as 0 and kills as 1, then subset the kills from the GT GPS clusters, and coded deer as 0 and moose as 1. We developed a candidate set of models with all combinations of GPS cluster characteristics (Table 1 ) using logistic regression modeling. We excluded any models containing covariates with a correlation coefficient >0.7. Latitude was included in the species model selection because deer densities are expected to decrease with latitude (Dawe and Boutin 2016 ). For both the kill and species models, we selected the most parsimonious model within 2 delta Akaike's information criterion (AIC c ) of the top ranking model.

We visited 49 GPS clusters in the field in the winters of 2013 and 2014. At each, we searched for a prey carcass, blood, or other remains and then scored whether the location was a kill or a bed distinguished by areas in the snow where wolves had rested. At kills, we scored the species, sex, age class, and condition of the prey. As we were interested only in the components of GPS clusters when wolves were at a prey carcass, we did not distinguish between kills and scavenge locations. However, only two carcasses near roads appeared to have been caused by means other than wolves.

We measured a suite of variables at each GPS cluster (Table 1 ). The duration was calculated as the sum of all times between locations less than six hours at a GPS cluster. We tallied the number of returns to a GPS cluster as the number of consecutive locations separated by more than 6 h. The mean distance and maximum (radius) distance between all points and the GPS cluster centroid were calculated. We tallied the number of times the wolf returned and created a new cluster within the space window (300 m) of another cluster within 30 d. We attributed the pack size and the latitude to each GPS cluster.

We estimated the locations of wolf‐caused moose mortalities using GPS clusters in wolf telemetry (Webb et al. 2008 ). We followed the method developed by Knopff et al. ( 2009 ), which established GPS cluster centroids as the geometric mean of points within a space‐time window. Seed centroids were established when a minimum of three points and three hours (given variable fix rates) had past. Once a single centroid was established, points were added chronologically and the centroid mean was adjusted. As the centroid moved, we corrected which points were included in the GPS cluster using the space‐time window. When two centroids from the same wolf fell inside the sampling window, they were combined into a new GPS cluster centroid using all points from each. We used a 300 m radius and four‐day space‐time window (Webb et al. 2008 , Lake et al. 2013 ). We allowed GPS cluster durations to extend beyond the four‐day temporal sampling window to a maximum of thirty days (Webb et al. 2008 ) as long as points added chronologically were not more than four days past the end of the GPS cluster (Knopff et al. 2009 ).

In the winter of 2012–2013, collars were replaced where they had failed and new wolves were captured where individuals had died. In the second winter, GPS locations of the previous year were used to establish pack boundaries to assess whether or not all packs in the area had a collared wolf. We discovered two distinct gaps: one situated between Fort McMurray and the oil sands mines west of the Athabasca River and the second south of oil sands mines and north of the Steepbank River on the east side of the Athabasca River. Both areas were searched using two helicopters for a full day. We did find limited wolf tracks indicating that wolves had moved through the area but could not find any sign of wolves using this habitat as territory. Based on our extensive efforts, we are confident that we placed collars on wolves from all packs in our study area.

We attempted to collar at least two wolves in every pack in the area covered by moose GPS telemetry. Wolves were captured and equipped with Iridium GPS collars (Lotek, Newmarket, Ontario, Canada) in winter of 2011/2012 using aerial net‐gunning following the Wildlife Animal Care Committee Class Protocol #761 (University of Alberta, ACUC Study ID AUP00000040). A total of 41 wolves were captured and collared. Wolf collars were programmed to fix the wolf's location at variable intervals depending on the time of year. In winter, the fix rates varied from 10 min to three hours. We estimated the number of wolves in the pack using counts of unique individuals during the collaring process.

The AOSR is comprised of the broader Athabasca watershed surrounding the Athabasca oil sands mines. The mines lie along the Athabasca River approximately 20 km north of the town of Ft. McMurray between 56.9° and 57.4° N and −111.0° and −112.0° E. The mining footprint consists of large extraction pits (325 km 2 ), tailings ponds (248 km 2 ), and associated large facilities (66 km 2 ). The area surrounding the mines has densities of linear features including seismic, transmission, and pipelines as high as 18 km/km 2 . Forest surrounding the mines is made up of peat ( Sphagnum sp.)‐forming wet areas with variable black spruce ( Picea mariana ) and tamarack ( Larix laricina ; 33%), uplands of aspen ( Populus tremuloides ), white spruce ( Picea glauca ), and jackpine ( Pinus banksiana ; together 30%). Forestry is minimal in AOSR. Moose densities, estimated with random block surveys (Gasaway et al. 1986 ), are low in the AOSR, but vary spatially between 0.04 and 0.15 km −2 (90% CI ranged from 21% to 41% of the estimate) in the three wildlife management units overlapping our study area as measured between 2008 and 2013 (Morgan and Powell 2008 , 2009 , 2010 , Alberta Environment and Sustainable Resource Development 2013 ).

Model predicted effect of relative moose density on the frequency of kills (relative predation rate) as a function of distance to rivers. Effects were estimated using logistic regression that compared the locations of kills to random available locations across the study area against an interaction of distance to mines and a spatial index of moose density. Distances were transformed such that close distances had the largest values. Grey areas are 95% confidence intervals.

Model predicted effect of relative moose density on the frequency of kills (relative predation rate) as a function of distance to mines. Effects were estimated using logistic regression that compared the locations of kills to random available locations across the study area against an interaction of distance to mines and a spatial index of moose density. Distances were transformed such that close distances had the largest values. Grey areas are 95% confidence intervals.

Our top‐ranked predation rate model contained interactions between moose density and rivers, distance to mines, weighted linear feature density, and forest cover type (Appendix S1 : Table S3). Wolves killed moose more frequently with decreasing distance to mines and rivers (Table 3 ) such that the relative predation rate increased significantly near mines (Fig. 3 ) and rivers (Fig. 4 ). Wolves killed moose more frequently in upland forest (Table 3 ). Fifty‐nine percent of kills occurred in upland forest, which only made up 45% of the available area, whereas only 19% of kills occurred in wetland forest, which made up 43% of the study area. Frequency of kills also decreased with increasing density of weighted linear feature density (Table 3 ). However, when compared to moose density, the number of kills in upland forest and along the linear feature gradient translated into very weak (insignificant effect size) changes to the relative predation rate.

Map of the Athabasca oil sands clipped by our study area. Locations of moose killed by wolves were estimated from wolf GPS data and ground‐truthed kill locations. The relative density of moose was calculated as a kernel density estimate (kde) of sightings from Alberta Environment and Parks moose density stratification flights. The kde was calculated with an 8 km search radius to match the mean moose home range diameter in our study area.

Using the top kill model, we distinguished between bed and kill GPS clusters across our entire GPS dataset, yielding 988 kills. From this subset, we applied our species model, yielding 199 moose kills. After combining GPS clusters from wolves at the same kill, we were left with 153 unique GPS clusters describing moose kills across the 10 wolf packs. After subsetting to our study area polygon, we were left with 129 locations of moose killed by wolves (Fig. 2 ).

Of the 49 GPS clusters visited in the field, 29 were beds. Twelve of the 20 kill GPS clusters were of moose and seven of deer, with one large‐bodied carcass un‐identified. In the selected kill model (Appendix S1 : Table S1), GPS clusters defining locations of ungulates killed by wolves were distinguished from locations of wolf beds by longer durations (Appendix S1 : Fig. S1) and larger mean distance of points to the GPS cluster centroid (Appendix S1 : Fig. S2). The bootstrapped mean AUC ROC score was 0.87 yielding an optimal classifying cutoff of 0.42. The kill model successfully classified an average 84% of kills from beds and 89% of beds from kills in bootstrapped k ‐fold evaluation. The selected species model (Appendix S1 : Table S2) distinguished clusters of moose kills from those of deer kills by longer durations (Appendix S1 : Fig. S3) and smaller mean distance of points to the GPS cluster centroid (Appendix S1 : Fig. S4). The model with equal parsimony to our selected species model, which contained duration and cluster radius, was also within two AIC c of the top model (Appendix S1 : Table S2). To maintain a more general model for predicting kill type, we selected the model with the higher likelihood as opposed to model averaging. With a bootstrapped mean AUC ROC score of 0.85 and an optimal classifying cutoff of 0.49, the species model successfully classified 85% deer from moose kills and 90% of moose from deer kills.

Discussion

Comparing locations of moose killed by wolves to an index of moose density revealed dynamics in the relative predation rate of moose driven by natural and anthropogenic landscape features. Our approach allowed inference concerning both the distribution of kills, which revealed areas where wolves aggregated while hunting, and how that distribution contributes to the relative predation rate of moose across space. Our results indicate that the influence of predation on the moose population is not uniform across the Athabasca oil sands, with areas of relatively high and low rates of predation, and that these areas have changed with building of mining features. The wolf aggregative response demonstrated avoidance of high density of linear features near Ft. McMurray but not to the extent that the predation rate was altered in these areas. On the other hand, areas in proximity to oil sands mines and tailings ponds exhibited increased wolf kills above that predicted by moose density, indicating that predation rates have increased near human‐disturbed areas.

As predicted, variation in the predation rate was better explained by gradients in anthropogenic than natural landscape features. Predators and prey compete in an evolutionary arms race to maximize or minimize spatial overlap (Sih 2005). In a system with relatively stable population dynamics, such as wolves and moose in AOSR (Fuller and Keith 1980, Hauge and Keith 1981, Wasser et al. 2011), it is unlikely that either species wins the race over long time spans. Long‐standing features such as forest cover and rivers are likely therefore to exhibit relatively uniform predation rates across space. For instance, despite the larger frequency of moose kills in upland forest cover over wetlands, the predation rate was not different between these two forest types. The two landscape gradients over which predation rate did vary spatially, distance to rivers and mines, had comparable effect sizes but the effect on distance to mines on predation rate dynamics mines was stronger. The imposition of oil sands mines has provided wolves with an advantage in the arms race with potential population consequences for moose.

Previous work has demonstrated that human disturbance increases wolf hunting efficiency. Wolves use linear features (James and Stuart‐Smith 2000) to move faster and further in a day (Dickie et al. 2017). The edges of mines are similar to linear features in that they are open and can be long and straight, characteristics previously shown to facilitate wolf movement (Dickie et al. 2017). For a coursing predator that uses a large area to encounter prey, such features facilitate movement and should increase encounters. For moose living near a mine, the edge of a pit, tailings pond, or fence presents a consistent feature past which there are reduced escape opportunities. We speculate that the large barriers such as mines provide a hunting advantage to wolves because unlike the edge of a territory, which is a boundary to wolves (Mech 1977) but not their prey, mines constrain both species, thereby increasing overlap near mines. This effect has been observed in hunting by African wild dogs (van Dyk and Slotow 2003, Davies‐Mostert et al. 2013) and demonstrated experimentally with predatory shrimp in which co‐aggregation with their prey along microcosm walls increased predator attack rates (Bergstrom and Englund 2004). Muhly et al. (2015) speculated that restriction to the movement of woodland caribou (Rangifer tarandus) due to impermeability of human disturbance increased the probability of encounter with predators. Additional work examining predator and prey movement along boundaries is needed to understand this effect in vertebrate systems.

Predicting the relative predation rate as a function distance to mines (Fig. 3) demonstrated increased frequency of kills with increasing moose density within several km of mines. Given that, on average, a moose home range in our area is approximately 8 km in diameter, individual moose near mines will experience increased exposure to predation across their entire home range. Increased predation rates near mines will reduce moose densities there, driving dispersal of juvenile moose to the area (Pulliam 1988). Such a novel source–sink would be unsustainable if mortality near mines is higher than can be replaced by dispersal, putting downward pressure on the moose population in AOSR. Consequences for the moose population may be particularly concerning considering the large footprint of the mine edge (975 km). However, the top predation model demonstrated that distances far from mines and rivers actually exhibit a negative predation rate such that wolves are not killing moose in areas with high relative moose density away from mines (Fig. 3). It is possible that the imposition of the mines has further aggregated wolves in some areas of the AOSR, decreasing the relative predation rate in others. Our results do not predict specific changes to the absolute predation rate of moose in AOSR over time, and additional work is needed to assess moose population dynamics over longer time periods.

Wolves may be trading off food and perceived risk in areas near high‐intensity human disturbances. Overall, there were fewer kills in areas of high linear feature density near the city of Ft. McMurray, but wolves killed moose there in proportion to other areas. We suggest that wolf use of areas near human disturbance in AOSR is context dependent (Haswell et al. 2016). Wolves will use areas in proximity to human presence where moose density is highest. Avoidance of human disturbance by predators can be thought of as anti‐predator behavior (Frid and Dill 2002). However, the benefits of such behavior will be weighed against costs (Fryxell 1991). If moose near human disturbance are unavailable to wolves, wolves using those areas potentially incur reduced access to food. If that loss limits wolf survival for packs near Ft. McMurray, they would likely begin using areas near human disturbance because the advantage of avoiding the disturbance does not outweigh the disadvantage of lost hunting opportunities (Hebblewhite and Merrill 2009). Our results demonstrate that for wolves, access to prey is more important than avoiding human disturbance.

The aggregative response of wolves is an important component of predation rates, and measuring how it varies in space provides a deeper understanding of how the total proportion of wolves’ prey is killed. Several recent studies have demonstrated that wolves adjust aggregation behavior in response to human disturbance by either avoiding (Hebblewhite et al. 2005, Muhly et al. 2011, Rogala et al. 2011) or using (James and Stuart‐Smith 2000, Latham et al. 2011, Dickie et al. 2017) human disturbance, particularly linear features. By estimating the relative predation rate of wolves on moose in the AOSR, we have illustrated that removal of large areas of habitat, which creates boundaries and alters the amount of space available to wolves and their prey, facilitates wolf predation with potential prey population and trophic consequences. Considering that oil sands mines are large and in places, used intensely by humans, the expected outcome may be wolf avoidance of mining features that creates refugia for moose. Our results highlight the un‐predictability of the effects of human disturbance and the importance of investigating novel types and magnitudes of disturbances on predator–prey interactions.