Studying predator–prey dynamics requires an understanding of both predator diet and prey availability (Fryxell and Lundberg 1994, Křivan and Sikder 1999). Predators can affect other species directly through predation events or indirectly through competition or trophic facilitation with other predators, trophic cascades and apparent competition (Lima and Dill 1990, Schmitz et al. 2000, Wilmers et al. 2003, Kortello et al. 2007, Wittmer et al. 2013). Insights into predator–prey interactions are important for understanding the factors behind prey decline, particularly for small populations at risk of extinction (Sinclair et al. 1998, Wittmer et al. 2005, Mech and Fieberg 2014).

Gray wolves Canis lupus have been studied to examine both the direct and indirect effects of predation on other species, as well as conflicts with humans through predation on livestock (Wilmers et al. 2003, Kortello et al. 2007, Wittmer et al. 2013, Nelson et al. 2016, Santiago-Avila et al. 2018). Wolf predation has been studied on wild ungulates such as deer (Odocoileus spp.), elk Cervus elaphus, and moose Alces alces, and may influence population dynamics of prey (Huggard 1993a, Hebblewhite et al. 2002, Smith et al. 2004). Less-commonly studied are interactions between wolves and plains bison Bison bison bison, as there are fewer areas where they co-exist (Carbyn and Trottier 1988, Smith et al. 2000, Jung 2011, MacNulty et al. 2014, Tallian et al. 2017). Historically, bison numbered in the tens of millions but were nearly driven to extinction in the late nineteenth century due to overhunting and habitat loss (Samson and Knopf 1994, Isenberg 2001). Bison have since increased in abundance, but most herds consist of farmed or captive animals with cattle gene introgression (Halbert and Derr 2006, Freese et al. 2007). Therefore, studies investigating wolf predation on wild, genetically-pure bison are limited to only a few areas in North America (Samson and Knopf 1994, Smith et al. 2000, Freese et al. 2007, Harvey and Fortin 2013).

The diet of wolves can be estimated using direct observation of kills (Boyd et al. 1994, Sand et al. 2005), scat and stomach content analysis (Floyd et al. 1978, Ciucci et al. 1996, Merkle et al. 2009), and stable isotope analysis (SIA; Szepanski et al. 1999, Derbridge et al. 2012, Stanek et al. 2017, O'Donovan et al. 2018). SIA quantifies the change in isotopic ratios as nutrients are consumed, metabolized and reorganized at each trophic level, which allows for the relative proportion of each prey item to be determined for a consumer (DeNiro and Epstein 1978, 1981, Peterson and Fry 1987). SIA can examine diet at the individual, group or population level (Urton and Hobson 2005, Derbridge et al. 2012, O'Donovan et al. 2018). As well, tissues used in SIA can provide diet history spanning weeks (e.g. blood) to months (e.g. hair) to lifetimes (e.g. bone collagen), allowing researchers to explore temporal differences in diet (Chisholm et al. 1982, Tieszen et al. 1983, Hilderbrand et al. 1996, Darimont and Reimchen 2002, Hall-Aspland et al. 2005, Gómez et al. 2018).

We used SIA to examine seasonal and inter-annual variation in the diet of wolves for two packs that overlap with plains bison range in and around Prince Albert National Park (PANP), Saskatchewan. The Sturgeon River plains bison (SRPB) are one of only a few wild populations of plains bison in their historical range in Canada, and have declined by greater than 50% since 2005, when the population was estimated to be around 500 individuals (Merkle et al. 2015, Cherry et al. 2019). Current population estimates are lower than the target management threshold required to prevent losses to genetic diversity (Cherry et al. 2019). Disease and annual harvests have contributed to past bison mortality (Shury et al. 2009, Merkle et al. 2015), but the role of wolf predation in the SRPB decline is unclear. Population simulations for the SRPB indicate unsustainable harvest is likely the main factor limiting population recovery; however, these models assume predation rates similar to other bison populations (Cherry et al. 2019, Simon and Fortin 2019). No studies have directly investigated the role of predation by gray wolves on the SRPB population. As free-ranging populations of plains bison are reintroduced into their historical range, it is important to understand the effects of predation on bison population dynamics (Steenweg et al. 2016).

We hypothesized that wolf diet would vary seasonally, between summer and winter, due to variation in prey availability and vulnerability. We predicted that wolves would primarily consume wild ungulates in winter, as ungulate movements are inhibited by deep snow and body condition is generally poorer, making prey more susceptible to predation (Sweeney and Sweeney 1984, Telfer and Kelsall 1984). Pack cohesiveness is also higher during winter (Benson and Patterson 2014), potentially allowing wolves to consume larger prey such as bison (Zimen 1976, Metz et al. 2011). During summer, we expected beaver Castor canadensis to contribute more to wolf diet, as they are more accessible outside of their lodges. In summer, there is also overlap between pack territory and summer pasture of cattle, which may increase the occurrence of predation on livestock. We also hypothesized that diet would differ between wolf packs. Since wolves are territorial, and buffer zones that separate adjacent territories can be as wide as 1–2 km, overlap between wolves from different packs is uncommon and often results in inter-pack conflict (Mech 1977, White et al. 1996). We predicted bison would constitute a higher proportion of diet for wolves that show greater overlap with the SRPB range (Bergeson 1993).

Material and methods

Stable isotope analysis

SIA of blood and hair samples were performed at the Chemical Tracers Lab (University of Windsor, Windsor, ON, Canada). Blood samples were freeze-dried, homogenized and lipids were removed, since variation in lipid concentrations can influence measurements of carbon isotope ratios (Rau et al. 1992, Post et al. 2007, Logan et al. 2008). Two milliliters of a 2:1 chloroform methanol mixture was added to dried blood samples, which were then vortexed for two seconds and placed in a water bath for 24 h at 30°C. Samples were centrifuged, and the chloroform methanol solution was drained. Another two milliliters of the chloroform methanol solution was then added, and samples were centrifuged and drained again. We removed underfur from hair samples and rinsed samples under a ventilation hood in a 2:1 chloroform/methanol solution to remove fine debris and oils (Darimont et al. 2007). All samples were air dried before weighing.

Figure 1.

The approximate range and territory size of plains bison and gray wolves, respectively, in Prince Albert National Park. The summer and winter range of bison was estimated using 95% minimum convex polygons. Location data for bison were collected from GPS-collars deployed from 2011 to 2017 over both summer/autumn (June–October) and winter (November–March). Location data for wolves were collected from GPS-collars deployed from 2013 to 2017 in winter (November–March).

Samples were weighed (∼1 mg) and placed into tin capsules for continuous-flow mass spectrometry analysis using a Delta V Advantage mass spectrometer coupled to a Costech 4010 elemental combustion system and a ConFlo IV gas interface. During continuous-flow mass spectrometry analysis, samples were combusted, resulting in the separation of CO₂ and N₂, which we used to quantify isotopic ratios (Fry 2006). Isotope values are expressed in delta notation as:

where X is ¹³C or ¹⁵N, and R is ¹³C/¹²C or ¹⁵N/¹⁴N. The standards used in SIA are Pee Dee Belemnite limestone for carbon, and atmospheric N² for nitrogen (DeNiro and Epstein 1978, DeNiro and Epstein 1981).

Statistical analysis and stable isotope mixing models

We used guard hairs and red blood cells (RBCs) for SIA of wolves. Isotopes in hair do not turnover, and thus reveal the diet of an individual when it is growing (Roth and Hobson 2000). Wolves undergo an annual moult in the spring, and guard hair is grown from summer to autumn (Young and Goldman 1964, Darimont and Reimchen 2002, Darimont et al. 2009). Therefore, guard hair reflects the summer and autumn diet of wolves. Cellular components of blood differ in their isotopic half-lives, with plasma reflecting diet integrated approximately one to two weeks before collection, and RBCs representing the diet of the previous several months (Hilderbrand et al. 1996, Thomas and Crowther 2015, Rode et al. 2016). Since we used RBC samples collected from wolves in November to April, we expected δ¹³C and δ¹⁵N values to be more representative of wolf diet in late autumn and winter.

Aerial wildlife surveys, scat analysis and wolf kill site visits were used to select prey to include in our mixing models (Szepanski et al. 1999, Darimont and Reimchen 2002, Derbridge et al. 2012, O'Donovan et al. 2018). We used nonparametric Mann–Whitney U tests to examine if species were isotopically distinct, and accounted for error associated with multiple comparisons using a Bonferroni correction (Bland and Altman 1995, Cabin and Mitchell 2000). We used diet-hair and diet-blood discrimination actors (hair: δ¹³C = 2.6%, δ¹⁵N = 3.2%; blood: δ¹³C = 0.7%, δ¹⁵N = 2.6%) from a captive feeding study on red foxes Vulpes vulpes (Roth and Hobson 2000) to account for isotopic change through trophic levels (Del Rio and Anderson-Sprecher 2008, Parnell et al. 2013, McLaren et al. 2015). To augment our study, we used data analyzed from wolf scat as priors for our mixing models (Merkle et al. 2014). Wolf scat was collected opportunistically during summer and winter of 2012–2013 and analyzed for the percent frequency of occurrence of prey. Mammalian hairs in each scat were identified by microscopic examination of the cuticular pattern, the medulla and the cross section. Prey biomass estimates were calculated using a regression equation that accounted for the ratio of indigestible to digestible remains (Floyd et al. 1978, Weaver 1993; Table 1).

We ran Bayesian stable isotope mixing models to estimate wolf diet during summer and winter from 2011 to 2017 (Moore and Semmens 2008, Semmens et al. 2009). We included a process^*residual error term in our models (Stock and Semmens 2016a). For each model we ran three chains that were 100 000 iterations long, with a burn-in of 50 000 and thinned every 50th iteration. We used Gelman–Rubin diagnostic tests to assess convergence for each model, not allowing more than one value above 1.1 for all variables in the model, as well as trace plots produced by MixSIAR (Stock and Semmens 2016b). We also used a pairs plot to examine tradeoffs between any strongly correlated species before interpreting diet results (Stock and Semmens 2016b). All models were run using JAGS and R software (Plummer 2003, < www.r-project.org>) and the R package MixSIAR (Stock and Semmens 2016b). We presented mixing model results as median ranges in the text of the results section and the 95% credible intervals (CI) for the median ranges in Table 3, 4.

Table 1.

Percent of total biomass (%) of prey items from wolf scats. Wolf scats (n = 465) were collected opportunistically during summer and winter of 2012–2013 in the southwest corner of Prince Albert National Park, Saskatchewan. Mammalian hairs in each scat were identified by microscopic examination of the cuticular pattern, the medulla and the cross section. Biomass estimates were calculated using a regression equation, which accounted for the ratio of indigestible to digestible remains (Floyd et al. 1978, Weaver 1993).

Results

SIA of summer and winter wolf diet

For combined pack data, the percentage of bison, deer/elk and moose consumed in summer was consistently high for all years (Table 3; bison: median range = 26–39%; deer and elk: 21–24%; moose: 16–33%). Wide 95% CIs for the posterior distribution of models resulted in considerable overlap between prey. Cattle (3–7%) and beaver (1–9%) comprised a lower percentage of wolf diet in summer for all years. Both wolf packs had similar prey contributions in summer. The contribution of bison to the diet of the Amyot pack was highest in 2011, 2013 and 2014 (Table 3).

Table 2.

Means (x̄) and standard errors (SE) of δ¹³C and δ¹⁵N values estimated from hair and blood tissue of wolves, and hair tissues from diet sources. Hair and blood samples were collected from 2011 to 2017 in Prince Albert National Park, Saskatchewan.

In winter, we found tendencies towards higher percentages of deer/elk (40–49%) consumed by wolves for all years, compared to summer (Table 4). We found a similar percentage of bison consumed in winter (25–45%) compared to summer, and almost no contribution from cattle or beaver in winter. Bison comprised a higher percentage of diet for the Amyot pack (26–40%) compared to the Nesslin pack (20–23%), while the median percentage of moose in wolf diet (19–35%) was higher for the Nesslin pack in 2013, 2014 and 2016.

Discussion

We found wolf diet differed between summer and winter, with a more even contribution of wild ungulates to wolf diet in summer. Deer/elk comprised the highest percentage of wolf diet in winter, supplemented by bison and moose. While we found evidence that wolves consumed beaver and cattle in summer, diet contributions were comparatively small. Relative proportions of prey contributions were similar among years. Bison comprised a higher percentage of the diet for the Amyot pack in both summer and winter, while the contribution of moose in winter was higher in some years for the Nesslin pack.

Wolves were centered on ungulates in the isotope mixing space, which is supported by other studies on wolf diet in the region Of the wolf kill sites visited (n = 270) in PANP from 2013 to 2017, 215 had prey remains: 57.2% (n = 123) were white-tailed deer, 18.2% (n = 39) were moose, 8.4% (n = 18) were deer for which genetic testing was inconclusive, 7.9% (n = 17) were mule deer, 5.1% were bison (n = 11), 1.9% were black bears (n = 4) and 1.4% (n = 3) were elk (Irvine 2019). Urton and Hobson (2005) performed SIA on wolves in central Saskatchewan, which overlapped with our study area in PANP, and found elk and deer dominated wolf diet, although they did not include bison in their mixing model.

The main prey of wolves shows considerable variation across regions, and includes moose (Ballard et al. 1987, Mech and Fieberg 2014), deer (Ballard et al. 1987, Huggard 1993a), elk (Huggard 1993a, Smith et al. 2000, 2004, Hebblewhite et al. 2002, Kortello et al. 2007, Atwood et al. 2007), and caribou Rangifer tarandus (Ballard et al. 1987, Merkle et al. 2017). Wolves may also supplement their diet with smaller prey that show a minor contribution to wolf diet when included in mixing models (Szepanski et al. 1999, Milakovic and Parker 2011, Derbridge et al. 2012, O'Donovan et al. 2018). The primary prey of wolves can be dependent on habitat overlap between predator and prey, prey density, and thus the encounter rate between wolves and prey (Huggard 1993a).

Wolves are opportunistic and prey switching can occur as a result of a number of factors (that are not mutually exclusive), including in response to seasonal and environmental changes (Nelson and Mech 1986, Szepanski et al. 1999), the occurrence of other predators (Kortello et al. 2007, Merkle et al. 2017), prey vulnerability (Bergman et al. 2006, Garrott et al. 2007, Metz et al. 2012), or changes in the abundance of prey (Garrott et al. 2007, Fortin et al. 2015, Sand et al. 2016). We found that winter diet for wolves centered on wild ungulates, particularly deer, which is what we had predicted. Deep snow (>70 cm) can inhibit ungulate movement, restricting access to forage and increasing energetic output, making prey more susceptible to predation by wolves (Telfer and Kelsall 1984, Huggard 1993b, Mech et al. 2001, Smith et al. 2004). We expected bison contributions to be higher in winter relative to bison contributions in summer, as wolf pack cohesiveness is higher (Benson and Patterson 2014), and deeper snow may make bison more susceptible to predation by wolves. However, bison harvests most often occur in the late summer and autumn, and it is likely that wolves scavenge gut piles and/or consume bison that are wounded or stressed during unsuccessful hunting attempts (Tallian et al. 2017), which may have increased diet contributions of bison during summer. Beaver and cattle are less available to wolves in winter, as beavers reside in lodges and forage under the ice (Benson and Patterson 2014), and cattle are on their winter pasture. In summer when there was greater overlap between wolves and cattle pasture, we found little evidence of cattle in wolf diet. Other studies have shown wolf predation on livestock occurs infrequently when wild ungulates are readily available (Meriggi and Lovari 1996, Oakleaf et al. 2003, Treves et al. 2004).

Variation in diet may also be observed between wolf packs due to differences in prey availability. Prey abundance can differ between packs as a result of habitat and landscape heterogeneity between regions, which can influence foraging and migration patterns of prey (Ballard et al. 1987, Gustine and Parker 2008, Fortin et al. 2015, Muposhi et al. 2016). Variation in habitat type/quality not only influences prey availability and density, but it can also affect prey vulnerability by modifying prey detection, access, and/or the success of an attack (McPhee et al. 2012a, b, Torretta et al. 2017). We found bison comprised a higher percentage of diet for the Amyot pack compared to the Nesslin pack. There is greater overlap between the SRPB population and the Amyot wolf pack (Bergeson 1993), which would have increased the encounter rate between wolves and bison. In addition, landscape characteristics within the SRPB range, including flat, open meadows, may have increased vulnerability of bison to wolf predation relative to forested areas (Bergeson 1993, McPhee et al. 2012b, Torretta et al. 2017). Finally, the Amyot pack was larger than the Nesslin Pack, which may increase the acquisition and consumption of larger prey, such as bison, for wolves in the Amyot pack, due to benefits of cooperative hunting (MacNulty et al. 2014, Tallian et al. 2017).

Table 3.

Posterior median estimates and 95% credible intervals (CI) of summer diet proportions for Amyot, Nesslin and combined wolf packs, from Bayesian stable isotope mixing models. Wolf diet was estimated from 2011 to 2017 in Prince Albert National Park, Saskatchewan. Years that data were unavailable are indicated by N/A.

Overall, SIA, kill site investigation and scat analysis reveal that wolves are consuming bison in PANP. There was considerable variation in the posterior distribution of our stable isotope mixing models, which confounds our interpretation of prey contributions. In our mixing models, we incorporated uncertainty in both our isotope and red fox discrimination factors. While controlled feeding studies have been performed on wolves (Derbridge et al. 2015, L'Hérault et al. 2018), red fox discrimination factors have been commonly-used for wolf diet studies (Ballard et al. 1987, Derbridge et al. 2012, Rode et al. 2016). Derbridge et al. (2015) found that specifying wolf fractionation rates made little practical difference when estimating wolf diet. The red fox discrimination factors placed wolves within the prey hypervolume in iso-space, suggesting that it is accurately representing wolf diet-tissue discrimination. We also incorporated priors into our models using results from wolf scat analysis, which often improves mixing model outputs by guiding parameter estimates and reducing the variance of prey contributions (Moore and Semmens 2008). However, there was considerable overlap among the posterior distributions of our mixing models, which was most likely a result of the similarity of our sources in the isotope mixing space and is common when using SIA to reconstruct carnivore diet (Derbridge et al. 2012).

Table 4.

Posterior median estimates and 95% credible intervals (CI) of winter diet proportions for Amyot, Nesslin and combined wolf packs, from Bayesian stable isotope mixing models. Wolf diet was estimated from 2011 to 2016 in Prince Albert National Park, Saskatchewan.

Deer and elk could not be isotopically separated, as they display high dietary overlap of browse and grass species, particularly during summer (Hansen and Reid 1975, Krysl and Bryant 2016). While moose and beaver spend a large portion of time foraging in aquatic habitats, beaver may consume greater quantities of browse, leading to higher δ¹³C values (Drucker et al. 2003). Bison had higher δ¹³C, indicating they are foraging primarily on C₃ plants, which is expected in an aspen parkland ecotone (Chisholm et al. 1986). The observed overlap between δ¹³C and δ¹⁵N signatures is normal when prey occupy a similar ecological niche, particularly when compared with δ¹³C and δ¹⁵N values from other ecosystems (Chisholm et al. 1982, Szepanski et al. 1999, Webb et al. 2017). Therefore, caution should be used when reconstructing the diet of a terrestrial carnivore that shows a considerable overlap in prey isotopic signatures, as this may make conclusions difficult to discern. Using other methods of diet analysis (e.g. direct observations of kills and/or scat or stomach content analysis) to support or refute conclusions drawn from SIA can be helpful in this regard.

In conclusion, bison constitute a smaller proportion of wolf diet when compared to other ungulates, which is consistent with results from other wolf–bison predator–prey systems (Smith et al. 2000, 2004, Jędrzejewski et al. 2002). However, predation in combination with other mortality factors, such as disease and harvest, may lead to population decline and increase extinction risk of the SRPB population (Joly and Messier 2004, Sigaud et al. 2017, Simon and Fortin 2019). Extinction risk for the SRPB was highest in population simulations when the effects of harvest and disease were combined (Cherry et al. 2019). But, small populations can be more vulnerable to extinction caused by predation (Festa-Bianchet et al. 2006). Therefore, management of the SRPB population should focus on reducing harvests, rather than wolf predation, to assist with the population recovery of the SRPB.