Minimizing wolf-livestock conflicts requires identifying conditions placing livestock at risk and focusing adaptive management at a local scale. Gray wolves Canis lupus began recolonizing Washington in 2008. We used generalized linear mixed models to investigate characteristics of wolf pack territories in Idaho and Montana from 1991–2008 (n = 137) and predict cattle Bos taurus depredation risk for current and probable wolf-occupied areas in Washington. Cattle depredation risk increased with cattle abundance and if the pack depredated cattle the previous year. When models were applied to wolf pack territories in Washington from 2008–2016 (n = 43), 3 of 7 (43%) depredating pack territory/years were predicted at ≥61% depredation risk. During the summer grazing season (1 May – 31 October) when most cattle depredations (97%; n = 34) occurred in Washington, cattle east of the Cascade Mountains were often on grazing allotments whereas cattle west of the Cascade Mountains were located on small, private farms. Thus, relative cattle abundance per grazing allotment and county likely represented cattle depredation risk east and west of the Cascade Mountains, respectively. County-wide and allotment cattle abundance forecasted 10.3% and 1.4% of probable wolf-occupied areas at ≥ 61% cattle depredation risk, respectively. These risk models and maps provide locations for federal and state wildlife managers to focus depredation prevention measures and a template for future analyses as wolves continue to recolonize Washington.
Large carnivore attacks on livestock and subsequent carnivore removal have led to conflict between farming communities and conservation practitioners worldwide (Woodroffe et al. 2005, Suryawanshi et al. 2013). Over the past few decades, human–carnivore conflict has received increasing attention as most large carnivore populations continue to decline (Ripple et al. 2014). Identifying areas of potential human– carnivore conflict where mitigation efforts would be most beneficial is an emerging practice for informing livestock management and carnivore conservation (Miller 2015).
The gray wolf Canis lupus (hereafter wolf) has one of the largest global ranges of any mammal (Ripple et al. 2014) and a considerable history of human–carnivore conflict. Livestock depredation was a primary factor in motivating the extirpation of wolves from most of the conterminous United States by the early–mid 1900s (Young and Goldman 1944). Following protection through the Endangered Species Act, wolf populations have increased in the north-western, southwestern, and Great Lakes regions of the United States (Wydeven et al. 2009, US Fish and Wildlife Service 2010, US Fish and Wildlife Service et al. 2016) and are federally delisted in some of these areas. In the northern Rocky Mountains (NRM; including Idaho, Montana and Wyoming) of the United States, wolves recolonized portions of their former range through natural dispersal from Alberta and British Columbia, Canada in the early 1980s (US Fish and Wildlife Service 1987) and reintroduction into Yellowstone National Park and central Idaho in the mid-1990s (US Fish and Wildlife Service et al. 1996). By the late 2000s dispersing wolves from the NRM and British Columbia, Canada successfully established packs in Washington and Oregon (US Fish and Wildlife Service et al. 2010).
Livestock depredations increased as wolf range expanded in the NRM (US Fish and Wildlife Service et al. 2013), however, confirmed depredations remained small (< 0.1%; Bangs and Shivik 2001, Wiles et al. 2011) in comparison to total livestock losses. From 1987–2015 (29 years) the NRM wolf population grew to an estimated 1704 individuals, during which time 2274 cattle and 4514 sheep depredations by wolves were documented in the region (US Fish and Wildlife Service et al. 2016). Depredation events involving cattle usually included 1–2 animals whereas sheep depredations often involved multiple sheep per incident (Muhly and Musiani 2009, US Fish and Wildlife Service et al. 2013).
In Washington, wolves are classified as endangered under Washington state law although packs residing within the NRM in the eastern third of the state are federally delisted (Becker et al. 2015). State recovery objectives include a minimum of 15 successful breeding pairs for three consecutive years, with ≥4 pairs in each of three recovery regions (Wiles et al. 2011). Non-lethal preventative methods (e.g. range riders and shepherds, guard dogs, carcass removal) and lethal removal of wolves are currently used to deter livestock depredation (Becker et al. 2016). Since lethal removal could delay recovery objectives, it is important to implement proactive, science-based preventative methods to mitigate livestock depredation by wolves as recolonization occurs.
Livestock depredation risk by wolves has been addressed at different spatial scales, influencing the interpretation of results. At a regional scale, both higher wolf abundance (Kaartinen et al. 2009, Wielgus and Peebles 2014) and lethal removal of wolves (Harper et al. 2008, Wielgus and Peebles 2014) increased livestock depredation risk in multiple studies. At the scale of a wolf pack territory larger pack size also increased the probability of recurring livestock depredations although, in contradiction to the regional scale, lethally removing wolves decreased the probability of future depredations events because fewer individuals remained to depredate within that local area (Bradley et al. 2015). Land cover types such as percentages of pasture (Kaartinen et al. 2009), coniferous forest (Treves et al. 2004), and tall vegetation (Davie et al. 2014), have also been shown to predict risk at multiple scales. However, the relationship between land cover types and depredation risk varies since larger spatial scales (i.e. townships; Treves et al. 2004) tend to reflect locations where livestock are grazed, whereas smaller scales (i.e. depredation sites; Davie et al. 2014) may represent habitat in which livestock are most vulnerable to predators. Therefore, considering the biology of the species and the spatial (and temporal) scales for which management actions are designed to effect are important components of study design.
Minimizing wolf–livestock conflicts requires identifying conditions that place livestock at risk and focusing outreach and interventions at a local scale (Treves et al. 2011). Risk mapping, commonly used to predict spatial locations where hazards will occur (Treves et al. 2011), has been increasingly implemented in human–carnivore conflict mitigation (Kissling et al. 2009, Marucco and McIntire 2010, Soh et al. 2014, Miller 2015, Miller et al. 2015). In the Great Lakes region, predictive maps were developed to identify areas at risk to livestock depredation by wolves in Wisconsin and Minnesota (Treves et al. 2004, 2011) and later used to predict risk for the Upper Peninsula of Michigan (using models derived in Treves et al. 2004 and Edge et al. 2011). To date no predictive risk maps exist for the NRM and neighboring states, which differ in landscape and livestock management practices from the Great Lakes region (Fritts et al. 1992). In addition, few studies assessing livestock depredation risk by wolves have included wolf pack demographics, prey abundance and landscape characteristics.
We investigated characteristics of wolf pack territories in Idaho and Montana from 1991 through 2008 to predict cattle depredation risk by a recolonizing wolf population in Washington. Although our primary goal was to quantify cattle depredation risk for Washington, there have been too few wolf packs and depredation events in the state to model depredation risk without incorporating data from neighboring states. Our first objective was to develop risk models for Idaho and Montana (i.e. training dataset) to test hypotheses that cattle depredations by wolves were associated with multiple factors including wolf demographics, cattle abundance, prey abundance and land cover types. We predicted that cattle depredation risk would increase for larger wolf packs with a history of cattle depredation in forested areas with higher cattle and wild prey abundance. Our second objective was to evaluate the accuracy of the training models by predicting cattle depredation risk by wolves in Washington from 2008 through 2016 (i.e. testing dataset). Our third objective was to map and forecast cattle depredation risk by wolves throughout potential wolf habitat in Washington.
Study area
Training and testing data were obtained from wolf-occupied areas of Idaho and Montana (1991–2008) and Washington (2008–2016), respectively (Fig. 1). Forecast maps from training models were developed for areas of > 40% probability of wolf habitat occurrence (Maletzke et al. 2016) in Washington. The study area (including Idaho, Montana and Washington) covers approximately 302 850 km2 (–124°50′56.36″– –109°18′19.66″W, 43°29′59.56″– 49°0′9″N) and exhibits a wide range of precipitation, topography and vegetation types typically found in the Cascade and northern Rocky Mountain ranges and the Columbia and Snake River Basins. In wolf-occupied areas of Idaho and Montana, elevation ranged from 224 to 3895 m, 13% between 224 and 1000 m, 56% between 1000 and 2000 m, 30% between 2000 and 3000 m, and 1% more than 3000 m. Forests cover 59% (58% evergreen, <1% deciduous and mixed woods) of wolf-occupied areas in Idaho and Montana, shrub-scrub 20%, farmlands 4% (2% pastures and 2% croplands), and urbanized areas 1% (Multi-resolution Land Characteristics Consortium 2015). Within potential areas of wolf recolonization in Washington, altitude ranges from 0 to 2491 m, 54% between 0 and 1000 m, 44% between 1000 and 2000 m, and 2% between 2000 and 3000 m. Forests cover 70% (67% evergreen, 1% deciduous, and 2% mixed woods) of the Washington study area, shrub-scrub 17%, farmlands <1%, and urbanized areas 3% (derived from Homer et al. 2015). Farmlands and urbanized areas are primarily found on flat terrain across the study area. The climate ranges from temperate rainforest in the west to semi-arid steppe in the east.
Wolves primarily prey on Rocky Mountain elk Cervus elaphas nelson (Hebblewhite 2005, Spence 2017), white-tailed deer Odocoileus virginianus (Kunkel et al. 2004, Spence 2017), mule deer O. hemionus (Kunkel et al. 2004, Spence 2017), moose Alces alces (Kunkel et al. 2004, Spence 2017), and in the Greater Yellowstone Area, American bison Bison bison (Smith et al. 2000). If recolonization in western Washington occurs, wolves may also prey on Roosevelt elk C. e. roosevelti and black-tailed deer O. h. columbianus.
Over the past two decades changes in human land use (Leu et al. 2008), game hunting regulations (Mackie et al. 1998), and management strategies for forest fires (Smith 2000), timber harvest (Wisdom et al. 2005), and large carnivores (e.g. cougars and wolves; Dusek et al. 2006) have affected the distribution and abundance of large ungulate populations in the region. White-tailed and mule deer are the most abundant and widely distributed ungulates in the study area (Montana Fish Wildlife and Parks 2016, Western Association of Fish and Wildlife Agencies 2016), followed by Rocky Mountain elk (Innes 2011).
Figure 1.
Study area including wolf packs in Idaho (2008), Montana (2008), and Washington, USA (2016) used in cattle depredation risk models and prediction maps for Washington.
The wolf population in the study area has steadily increased following recolonization and reintroduction (US Fish and Wildlife Service et al. 2016). In 2011, the NRM wolf population was delisted in all states except Wyoming and is currently managed by state and tribal wildlife agencies (US Fish and Wildlife Service et al. 2012). The Wyoming wolf population was delisted in 2017. As of December 2015, a minimum of 1704 wolves in 282 packs were documented in the Montana, Idaho and Wyoming portions of the NRM and a minimum of 90 wolves in 18 packs were documented in Washington (US Fish and Wildlife Service et al. 2016).
Domestic livestock production, primarily beef cow-calf Bos taurus and ewe-lamb Ovis aries operations, occurs on private farms and public lands outside of national parks and most wilderness areas (US Fish and Wildlife Service 1994). For many livestock operations, calves and lambs are born in late winter or early spring and, while still nursing, released on public and private allotments for the summer grazing season (1 May – 31 October). Although all domestic neonates are vulnerable during this time period, calves are commonly depredated more often than other age groups but lambs are not (Fritts et al. 1992, Oakleaf et al. 2003, Bangs et al. 2005, Bradley and Pletscher 2005).
Methods
Data collection
We compiled data on cattle depredations, wolf demographics, cattle abundance, ungulate harvest statistics, and land cover types for Idaho (2004–2008), Montana (1991–2008), and Washington (2008–2016). Data collected from Idaho and Montana only extended through 2008 because the wolf harvest seasons implemented in these states during 2009 (Carolyn et al. 2010, Mack et al. 2010) decreased the availability of pack demographic information. Due to the size of the study area and changes in wolf management from federal to state authority, we obtained data from multiple sources (Table 1).
Spatial boundaries of wolf pack (defined as ≥ 2 adult wolves traveling together; Wiles et al. 2011) territories were derived from very high frequency and global positioning system collared wolves. Annual territories were calculated by the overseeing wildlife agency using either the 95% kernel density estimation (KDE) method or the minimum convex polygon (MCP) method of home range analysis. We censored lone and dispersing wolves from the dataset.
We defined a depredation as a cattle mortality (calf, mature cow, bull or steer) with wolves as the confirmed or probable cause of death. Wolves were confirmed as the cause of death by the wildlife agency when there was reasonable physical evidence (e.g. fresh tracks, bite marks, feeding pattern) that the animal was killed by wolves (USDA Wildlife Services 2007). Probable depredations had some evidence to suggest depredation by wolves, but lacked sufficient evidence to confirm the cause of death (USDA Wildlife Services 2007). We included both confirmed and probable losses because implementation of non-lethal deterrents, wolf removals and livestock compensation programs were based on these numbers (Carolyn et al. 2009, Nadeau et al. 2009, Becker et al. 2016). A pack was considered a ‘depredating pack’ and included in analyses if it acquired ≥1 confirmed or probable depredations in a (Gregorian calendar) year on public or private land. A pack with no confirmed or probable depredations in a year was defined as a ‘non-depredating pack’ for that year, regardless of whether they depredated the previous or following year. These data were considered estimates since not all depredations were detected or reported.
Table 1.
Variables, descriptions, and sources for all response and predictor variables included in full models to assess cattle depredation risk by wolves in Idaho and Montana from 1991–2008 and develop risk maps for Washington, USA.
Sheep depredations were not included in these analyses because county-wide sheep abundance data were unavailable for multiple years during the study period and sample sizes were too small to analyze. Other stock (e.g. horses Equus caballus, llamas Lama glama and goats Capra aegagrus hircus) were also not included due to small sample sizes since wolves killed these species infrequently in the NRM.
Predictor variables of cattle depredation risk included multiple measures of wolf demographics, cattle and ungulate abundance, and land cover characteristics (Table 1). We defined years of pack residency as the minimum number of years the pack was confirmed prior to 2008 (Idaho and Montana) or 2016 (Washington). Because depredation risk may increase after a pack has learned to prey on livestock (Harper et al. 2005, Sime et al. 2007) we defined a previously depredating pack as having ≥ 1 confirmed or probable depredations in a previous year, regardless of whether a depredation occurred that year. If complete pack removal was implemented following depredations and wolves recolonized the area, we designated it as a new pack.
Multiple studies conducted on farms and private pastures concluded livestock depredation risk was affected by wolf abundance (Mech et al. 2000, Bradley and Pletschet 2005, Kaartinen et al. 2009), and adult wolf mortality (Harper et al. 2008). We calculated estimates of pack size, minimum number of adults, adult mortality and breeder mortality as the minimum number of individuals prior to the first depredation event (depredating pack) or 15 April, the average first date of the denning season in the NRM (Bradley and Pletscher 2005; non-depredating pack), in a calendar year. If the first depredation occurred prior to 15 April, minimum pack size estimates from the previous year were used. Pack reproduction was confirmed by visiting den and rendezvous sites to determine whether or not the pack produced ≥1 pups in a given year. We defined adult wolf mortality as the number of known mortalities (or translocations) prior to the first depredation event that calendar year for depredating packs and prior to 15 April for non-depredating packs. Since the exact age of a dead wolf is difficult and usually an estimate, we considered all dead wolves estimated at ≥ 1 years-old as adults. Breeder mortality was calculated as the number of known breeding wolf (total, male and female) mortalities prior to the first depredation event that calendar year in addition to losses after 15 April the previous calendar year. This definition accounted for breeders that reproduced but may not have remained long enough to teach pups (potential depredating yearlings the following year) to hunt. During the time of this study, most known wolf mortalities occurred post-depredation as the result of lethal control. In the early years of recovery translocation was used to mitigate wolf–livestock conflict prior to control efforts (Bangs et al. 1998). Although translocating wolves is not a form of mortality, it removes the wolf from its natal pack territory with a similar purpose to lethal control. Therefore, we separated adult or breeder mortalities due to lethal control (69% – government lethal control, 10J rule, and legal take) or translocation (<1%) from other causes of death (27% – ‘other human-caused’ such as illegal take, capture-related or vehicle collision; 3% – natural) since the former were directly related to cattle depredation.
We calculated cattle abundance as the estimated number of beef cattle in a pack territory per year. Cattle abundance was chosen over cattle density to provide a relevant unit of measure to wildlife managers and livestock producers. County-wide beef cattle estimates included animals grazed on private pastures in addition to public grazing allotments. Because pack territories may overlap multiple counties, we weighted cattle abundance per county by the area (km2) of wolf pack territory within each county to estimate the relative number of beef cattle available to wolves. When a pack territory overlapped multiple counties, we calculated the weighted sum of cattle abundance for all counties.
Some cattle are registered in one county but grazed on an allotment in a different county for the summer grazing season. To account for these discrepancies in the county-wide cattle data, we also calculated the total number of cattle on grazing allotments overlapping a pack territory per year.
Antlered deer and elk harvest statistics are often positively correlated with abundance (Wood et al. 1989, Dusek et al. 2006), therefore, we calculated catch-per-unit-effort (CPUE; no. of bulls or bucks harvested / hunter days) annually for elk and deer (mule, white-tailed and black-tailed) within each game management unit (GMU). Annual totals of bulls and bucks harvested and hunter days by GMU were obtained from each state wildlife agency and included both general and permitted hunt seasons. For deer species, all bucks were included since the definition of ‘antlered’ bucks differed between states and changed during the time period of the study. All elk bulls were included. If data were unavailable for a GMU in a particular year (e.g. no harvest season) we averaged total bucks, bulls and hunter days from the prior and following years. Because GMUs and pack territories vary in size and may overlap boundaries, we weighted prey CPUE per GMU by the area (km2) of wolf pack territory within each unit. When a pack territory overlapped multiple GMUs, we summed the weighted prey CPUE for all GMUs.
Livestock depredations occur where livestock grazing areas overlap wolf pack territories, which for our study area included both forest (grazing allotments) and pasture (private farms) land cover types. We calculated the proportion of forest, pasture and developed land cover types per pack territory. We used developed land cover from the National Land Cover Dataset (Multi-resolution Land Characteristics Consortium 2015) as a surrogate for housing unit density (United States Census Bureau 2015) since it was available at a finer temporal scale (i.e. every 5 years instead of 10 years).
Density of paved roads can affect wolf movement patterns (Thurber and Peterson 1994) although the predictive power of road density for territory establishment and livestock depredation risk is debated (Thiel 1985, Mech et al. 1988, Mladenoff et al. 1997, Treves et al. 2004, Mech 2006). Since the scale of this study is a wolf pack territory, road density would not likely identify areas of high human or wolf use within pack territories, and was therefore excluded.
Analysis and mapping
To discriminate high-risk pack territories from those with low cattle depredation risk, we included a comparison set of all KDE and MCP wolf pack territories with no depredation history from 1991–2008 (similar to Treves et al. 2011 and Venette et al. 2010). Only pack territories with available data for all predictor variables were maintained for analyses. The final dataset contained 98 depredating pack territory/years (28% of total pack territory/years) and 258 non-depredating pack territory/years (72% of total pack territory/years) in Idaho and Montana.
We developed cattle depredation risk models using generalized linear mixed models (GLMMs). GLMMs combine two statistical frameworks widely used in ecology, linear mixed models (which include both fixed and random effects, such as relative cattle abundance and wolf pack name, respectively) and generalized mixed models (which handle non-normal stochastic distributions; Bolker et al. 2009). A mixed modeling framework was chosen over traditional logistic regression because inter-annual observations within a wolf pack territory were not independent, therefore, wolf pack was included as a random (or grouping) effect in the intercept. Since baseline depredation risk for each wolf pack was unknown, we used marginal likelihood to calculate the random effects. Year was not included as a random effect in the slope because the sample size was not adequate for model convergence (Grueber et al. 2011). We addressed temporal autocorrelation using other methods.
We performed binomial (non-depredating or depredating pack territory/year) GLMMs in the glmmADMB package (Skaug et al. 2016) to assess how predictor variables (Table 1) were associated with cattle depredation risk by wolves. We conducted all statistical analyses in R ver. 3.2.2. (< www.r-project.org>, accessed 6 September 2015).
We included a predictor variable (or fixed effect) in a multivariate model if it was not collinear with a stronger predictor of cattle depredation risk from previous peer-reviewed wolf–livestock depredation studies. We tested remaining predictors for multicollinearity and removed a variable if the variance inflation factor (VIF) > 3.0 (Zuur et al. 2010, Cade 2015). Lastly, we standardized continuous predictor variables using centering procedures from Gelman (2008) to ease interpretation of interactions by placing binary and continuous predictors on a common scale (
= 0, SD = 0.5; Table 2). Spearman's rank correlation revealed high correlation (rs > ±0.7) between total wolf mortality and adult wolf mortality, total breeder mortality and breeder male and female mortality, and deer CPUE and elk CPUE. Breeder and adult male mortality can alter prey selection (Sand et al. 2006, Harper et al. 2008) and pack behavior (Mech 1999, Mech and Boitani 2003), therefore, we maintained total adult and breeder mortality for further analyses (Table 1). We combined deer CPUE and elk CPUE into one prey CPUE variable. We dropped the minimum number of wolf adults, pups and pup mortality because they were already represented in the total pack size (Kaartinen et al. 2009, Bradley et al. 2015) and pack reproduction (i.e. wolves denning; Bradley and Pletscher, 2005) variables. Percent developed cover was not included since the range within wolf pack territories was 0–2% and we removed percent pasture cover because most cattle depredations in Washington from 2008–2016 were verified on forested, public lands. We dropped years of pack residency from analyses due to multicollinearity with all continuous predictors.
Table 2.
Mean () and Standard deviation (SD) of continuous predictor variables used to assess cattle depredation risk by wolves in Idaho, Montana and Washington, USA. Standardized values (using methods in Grueber (2008)) were used for analyses. Sampling unit was an annual wolf pack territory.
The final GLMM model set was based on hypothesized relationships between predictor variables and cattle depredation risk by wolves (Johnson and Omland 2004; Table 3). The model set included ten variables (total pack size, pack reproduction, adult wolf mortality via lethal control, adult wolf mortality via other causes of death, breeder wolf mortality via lethal control, breeder wolf mortality via other causes of death, relative cattle abundance, prey CPUE, depredation the previous year and forest cover), four interactions (relative cattle abundance × prey CPUE, total pack size × adult breeder mortality via lethal control, total pack size × adult breeder mortality via other causes of death, adult wolf mortality via lethal control × depredation the previous year), and the random effect (wolf pack). The interaction between cattle and prey represents the hypothesis that cattle depredations may increase as wolves are attracted to the presence and abundance of prey in an area (Mladenoff et al. 1997, Treves et al. 2004, Bradley and Pletschet 2005). We also hypothesized the loss of breeding wolves prior to the denning season would decrease total pack size that year, and wolf packs which had depredated cattle the previous year may have experienced adult wolf mortalities via lethal control.
We ranked models using Akaike's information criterion corrected for small sample sizes (AICc) and the most likely model was determined if the Akaike weight (wi; the weight of evidence in favor of model i being the most likely model in the set) ≥0.90 (Burnham and Anderson 2002). If no model was the most likely and the top models were not nested (Arnold 2010), we calculated a model averaged prediction for each pack territory/year using the fewest number of models to attain a sum model set weight ≥ 0.90 (i.e. top model set; Burnham and Anderson 2002, Johnson and Omland 2004).
Because nearby pack territories may be similar in prey abundance and landscape characteristics and the distance between pack territories decreased as available habitat was recolonized, we tested the Pearson's residuals from the full model for spatial and temporal autocorrelation using global Moran's I (Getis and Ord 1992) and graphing autocorrelation through time (Venables and Ripley 2002), respectively. We internally validated the top training models using ten repeats of ten-fold cross validation and, for the binomial model set, the subsequent receiver operating characteristic (ROC) curves. Lastly, we examined the relationship between each predictor variable and cattle depredation risk while holding other predictors constant at their mean and assessed whether the model accurately predicted depredating packs in the two highest risk bins (i.e. intermediate risk 61–80% and high risk 81–100%).
To validate risk models for Washington, we used the top training models to calculate the probability of cattle depredation risk for all MCP wolf pack territories in the state from 2008–2016. We assessed the area under the ROC curve for model fit and whether the model accurately predicted depredating packs in the two highest risk bins. In the initial model set, a negative relationship between cattle depredation risk and percent forest cover was indicative of where most confirmed or probable cattle depredations were found on un-forested private pastures in Idaho and Montana (77% private, n = 223). Therefore, those models did not accurately predict risk for public, forested lands in Washington (32% private, n = 11). We removed forest cover to increase predictive accuracy of cattle depredation risk in Washington (see Supplementary material Appendix 1 for a table of the forest cover model set).
To forecast cattle depredation risk in areas of Washington when wolves occur or are predicted to recolonize, we predicted cattle depredation risk for a 9-km2 raster of > 40% probability of wolf habitat occurrence (Maletzke et al. 2016). National Parks and wilderness areas were removed since livestock grazing is unlawful and occurs intermittently in these areas, respectively. Cattle depredation risk was not forecasted for Indian Reservations as cattle and prey abundance data were unavailable. After these exclusions, the total area of forecasted cattle depredation risk by wolves was 51 174 km2 (or 29.7% of the total land area in Washington). Because this analysis was at the scale of a wolf pack territory, we calculated the probability of a cattle depredation for each pixel on a 933-km2 grid of hypothesized wolf pack territories (Maletzke et al. 2016) and used an 8.60-km2 moving window (half the radius of a hypothesized wolf pack territory) to average cattle depredation risk for each 9-km2 pixel.
Results
The full model included nine variables (total pack size, pack reproduction, adult wolf mortality via lethal control, adult wolf mortality via other causes of death, breeder wolf mortality via lethal control, breeder wolf mortality via other causes of death, relative cattle abundance, prey CPUE, and depredation the previous year), four interactions (relative cattle abundance × prey CPUE, total pack size × adult breeder mortality via lethal control, total pack size × adult breeder mortality via other causes of death, adult wolf mortality via lethal control X depredation the previous year), and the random effect (wolf pack). Graphical validation of the Pearson's residuals against fitted values indicated no heteroskedasticity issues in the full model. Since most pack territory/years (73%) contained no depredations, we tested the model for overdispersion (χ2 residual df = 0.377, p = 1.00) which indicated no zero-inflation. Low intercept variation (
= 2.25, SD = 1.50) showed the variation among wolf pack territories was indistinguishable from measurement error. No spatial or temporal autocorrelation was found in the Pearson’s residuals from the full model (global Moran's I z-score = –0.578, p = 0.561). The CATTLE-YR1DEP and CATTLE-MORTALITY models were in the top model set (Table 3). Since the CATTLE-YR1DEP model held the majority of the model set weight (wi = 0.818) and was a nested version of the CATTLE-MORTALITY model we only used the CATTLE-YR1DEP for further analysis. As predicted, relative cattle abundance and previous years' depredation were reliable predictors of a wolf pack depredating cattle. The difference in the log odds of a wolf pack depredating cattle was 1.940 (SE = 0.567, 85% CI = 1.124–2.757) for every additional cow in a wolf pack territory and 1.322 (SE = 0.452, 85% CI = 0.672–1.973) if a pack depredated the previous year.
Table 3.
Generalized linear mixed models used to assess cattle depredation risk by wolves in Idaho and Montana from 1991–2008 and develop risk maps for Washington, USA. Models were evaluated based on Akaike's information criterion for small sample sizes (AICc). Variable code descriptions are located in Table 1. The change in AlCc values between the top model and the proceeding models (ΔAICc), degrees of freedom (df), and model weight (wi) are included.
Ten repeats of ten-fold cross validation produced an average prediction error (the difference between the predicted and actual probability for each observation) of 0.309. However, an average of only 19% (n = 19) of predictions for known depredating packs were in the two highest risk bins, indicating underprediction of the top model. The average area under the ROC curve (AUC) was 0.732 and showed that for a specificity of 50% the average model sensitivity was 78%, whereas at a specificity of 75% the average sensitivity was 65%. Therefore, the predicted probability of cattle depredation in 50% (or 25%) of wolf pack territories would be accurate 78% (or 65%) of the time.
Cattle depredation risk in Washington
We predicted cattle depredation risk with the CATTLE-YR1DEP model for all MCP wolf pack territories in Washington from 2008–2016 (n = 43) to validate the training dataset on independent observations. Cattle depredations occurred in 7 pack territory/years (16% of total pack territory/years).
For the probability of cattle depredation using countywide cattle abundance, external validation identified 3 of 7 (43%) depredating pack territory/years in the two highest risk bins. A mean prediction error of 0.437 and an AUC of 0.583 indicated poor training model fit on the Washington dataset. For a specificity of 50% the model sensitivity was 57%, whereas at a specificity of 25% the model sensitivity was 41%. In contrast, for allotment cattle abundance external validation produced a mean prediction error of 0.321 and an AUC of 0.746, indicating fair training model fit on the Washington dataset. For a specificity of 50% the model sensitivity was 100%, whereas at a specificity of 25% the model sensitivity was 57%. However, only 2 of 7 (29%) depredating pack territory/years were predicted in the two highest risk bins.
Predicted probabilities were graphed to illustrate the relationship between predictor variables from the top model and cattle depredation risk (Fig. 2). Predicted depredation probability increased to an asymptote of 100% at about 6000 cattle per wolf pack territory (Fig. 2). The median cattle depredation probability for wolf packs with no depredation the previous year was around 5% whereas for wolf packs with a depredation the previous year cattle depredation probability was around 50% (Fig. 2).
To map cattle depredation risk for recolonizing wolves in Washington, we assumed no previous year depredations. Although wolf packs with a history of livestock depredation are often identified as higher risk to subsequent depredations (Harper et al. 2005, Sime et al. 2007), for the recolonizing wolf population in Washington predicting cattle depredation risk for newly established wolf packs with no known depredation history was more relevant.
During the summer months when most cattle depredations occur (Becker et al. 2016), cattle in wolf-occupied areas on the east side of the Cascade Mountains are often placed on grazing allotments whereas cattle west of the Cascade Mountains are located on small, private farms. Thus, relative cattle abundance per grazing allotment and county likely represented cattle depredation risk east and west of the Cascade Mountains, respectively. Average number of cattle per county and grazing allotment from 2008–2016 were used to estimate relative cattle abundance in each 9-km2 pixel. County-wide and allotment relative cattle abundance identified 10.3% and 1.4% of hypothesized wolf habitat within the two highest risk bins, respectively (Fig. 3). The average mapping error (mean difference between predicted probabilities for hypothetical pack territories and moving window averaged map pixels) was 13.1% (SD = 17.0) and 4.8% (SD = 8.6) for county-wide and allotment cattle depredation risk maps, respectively.
Discussion
The initial results showed that a decrease in forest cover best predicted cattle depredation risk in Idaho and Montana where most cattle depredations were found on un-forested private pastures, but did not adequately predict risk for forested public lands in Washington. The key findings from the model set with forest cover removed indicated relative cattle abundance in wolf pack territories and cattle depredation the previous year had the greatest effect on cattle depredation risk in Idaho and Montana. These results were similar to other wolf–livestock depredation studies associating increased depredation risk with increased cattle abundance (Mech et al. 2000, Bradley and Pletschet 2005, Hanley 2017), and packs with a history of livestock depredation (Harper et al. 2005). Fenced cattle pastures in Idaho and Montana were predicted to experience depredations if they contained >310 head of cattle (Bradley and Pletschet 2005). However, unlike this analysis, elk presence was the best predictor of pastures with depredations in that study (Bradley and Pletschet 2005). Multiple studies have concluded encounter rates between wolves and cattle likely increase as wolves are attracted to the presence and abundance of prey in an area (Mladenoff et al. 1997, Treves et al. 2004, Bradley and Pletschet 2005) whereas others found wolves switch to livestock when natural prey become scarce (Meriggi and Lovari 1996). In this study, prey CPUE was not a strong predictor of cattle depredation risk potentially indicating both mechanisms were present. For instance, although Bradley and Pletschet (2005) found local elk abundance increased cattle depredation risk, ungulate population declines in Montana following the severe winter of 1996–1997 also led to increased livestock depredations by wolves (Montana Fish Wildlife and Parks 2003). Since the scale of this study was the wolf pack territory and encompassed many pastures and grazing allotments, local prey abundances or declines were not represented.
The relationship between wolf lethal control and livestock depredation risk is currently debated. Studies conducted in North America at a regional scale found that lethal control increased or had little effect on subsequent livestock depredation rates (Musiani et al. 2005, Harper et al. 2008, Wielgus and Peebles 2014, Kompaniyets and Evans 2017). However, Poudyal et al. (2016) noted a decrease in subsequent depredations and concluded the pack scale may be more appropriate to address this question. A study conducted at the scale of a wolf pack territory in the NRM during the same time period by Bradley et al. (2015) indicated lethal control directed towards partial pack removal reduced the probability of subsequent livestock (primarily cattle and sheep) depredations in wolf pack territories by 29% over a five-year period. A Michigan wolf–livestock depredation study concluded lethal control was associated with an insignificant decrease in depredation risk at the control site, but found an insignificant increase in depredation risk that same year at sites within 5.42 km (Santiago-Avila et al. 2017). In this study, we found only a weak relationship between adult wolf mortality and cattle depredation risk at the scale of a wolf pack territory. The second-best model, CATTLE-MORTALITY (Table 3), contained the top predictors (relative cattle abundance and depredation the previous year) in addition to adult wolf mortality from both lethal control and other causes of death (e.g. illegal take, capture-related, natural). However, 1) the adult wolf mortality variables were weak predictors (i.e. 85% C.I.'s overlap zero), 2) the model weight was only 0.144, and 3) the top model was a nested version of CATTLE-MORTALITY with a much stronger model weight (wi = 0.818). Therefore, we used the more parsimonious top model with no spurious predictors (Arnold 2010). Adult wolf mortality from lethal control may have been an uninformative parameter because wolves in some packs were killed for depredating sheep as well as cattle, muddling model inference for depredation risk on cattle alone. The relationship between cattle depredation risk and adult wolf mortality from other causes of death are likely so varied and nuanced that deriving a clear relationship would be difficult and require additional analyses beyond the scope of this paper.
Figure 3.
Predicted cattle depredation risk for > 40% probability of wolf habitat occurrence (Maletzke et al. 2016) in Washington, USA. Colors categorize risk by 9-km2 pixels in five equal-sized bins for (A) county-wide relative cattle abundance and (B) allotment cattle abundance per 933-km2 hypothetical wolf pack territory. The moving window produced an average mapping error of 13.1% and 4.8% for county-wide and allotment depredation probabilities, respectively.
Relative cattle abundance on grazing allotments was a better predictor for cattle depredation risk in Washington than county-wide cattle estimates (AUC of 0.746 and 0.583, respectively). This likely happened because the majority of cattle depredations in Washington from 2008–2016 occurred on grazing allotments, whereas most cattle depredations in Idaho and Montana occurred on private farms. This trend is indicative not only of where cattle are located during the summer grazing season, but of differences in the number of cattle grazed on private farms. In wolf-occupied counties of Montana and Idaho, 5.54% and 3.34% of private farms graze ≥500 cattle, respectively (National Agricultural Statistics Service 2012a). In contrast, only 1.43% and 0.03% of private farms graze ≥500 cattle in probable wolf-occupied counties on the east and west sides of the Cascade Mountains of Washington, respectively (National Agricultural Statistics Service 2012a). Since these risk maps were only forecasted for areas of suitable wolf habitat (> 40% probability of wolf habitat occurrence) in Washington, wolves may overlap more cattle on private farms in marginal habitat as recolonization continues and suitable territories are filled (Brown 2016).
South-central Washington (Kittitas and Yakima counties) and localized areas in the Blue Mountains and east of the Cascade Mountains (Okanogan County) were identified as high risk using county-wide cattle abundance (Fig. 3A). Yakima county, located in south-central Washington east of the Cascade Mountains (in the Southern Cascades and Northwest Coast Recovery Region), is the second largest producer of beef cattle in the state with 15 414 cattle on 675 farms (National Agricultural Statistics Service 2012a). This hot spot may also predict sheep depredation risk since the largest number of sheep (6525 ewes and lambs, number of farms not disclosed; National Agricultural Statistics Service 2012b) and area of active sheep grazing allotments also occur in Yakima county. No cattle depredations occurred within hotspots identified as high risk in part because wolves have not recolonized south-central or western Washington. However, three cattle depredations were confirmed near these areas – one in eastern Okanogan county (North Cascades Recovery Region), two in northern Kittitas county (North Cascades Recovery Region), and one east of the Blue Mountains (Eastern Washington Recovery Region) – in 2012, 2015 and 2016, respectively.
Several hotspots of intermediate-high depredation risk occurred west of the Cascade Mountains using countywide cattle abundance (Fig. 3A). However, 99.17% of private farms in western Washington graze <100 cattle. Therefore, realized depredation risk may be lower than predicted by county-wide cattle abundance if smaller cattle herds closer to human residences experience fewer depredations (Mech et al. 2000, Bradley and Pletscher 2005). No wolf packs currently reside in these areas and we were unable to verify whether depredation risk was lower than predicted.
Active cattle grazing allotments only occur east of the Cascade Mountains, and hotspots in Okanogan, Ferry and southwestern Yakima counties were identified as intermediate depredation risk (Fig. 3B). These areas include densely clustered grazing allotments with >200 cattle per allotment (see Hanley 2017 for analysis of cattle depredation risk in grazing allotments). Wolves currently inhabit Okanogan and Ferry counties, and assuming no depredations the previous year and five wolves per pack territory, wolf-occupied areas in Ferry county overlapped areas identified as intermediate depredation risk.
Model uncertainty and limitations
Model predictions and associated risk maps should be interpreted with care. Top models from hypothesis-based model sets may contain different variables than model sets created using stepwise regression. Wolf packs with no confirmed or probable cattle depredations may have killed cattle but carcasses were not found or reported, contributing to model uncertainty and predictor variability. In addition, wolf packs with or without a history of cattle depredation may have depredated sheep or other livestock. These depredations were not accounted for in these analyses and may not result from the same pack processes.
Predicted depredation risk may differ from realized depredation risk since cattle abundance does not necessarily reflect cattle availability to wolves. Factors such as human presence have been associated with decreased wolf depredations on fenced cattle (Mech et al. 2000, Bradley and Pletscher 2005). In addition, the average relative countywide cattle abundance for wolf pack territories in Idaho and Montana was nearly half that of Washington wolf pack territories (Table 2), likely because the average wolf pack territory size in Idaho and Montana ( = 383.2, SD = 327.4 km2) was half that of recolonizing wolf pack territories in Washington ( 844.2, SD = 469.0 km2). This discrepancy is due to differences in sample size between the training and validation areas, pack territories shrinking as the wolf population increased in Idaho and Montana (Montana Fish Wildlife and Parks 2003), and a switch from very high frequency to geographic positioning system wolf collars in the latter years of the study. Therefore, predicted depredation probabilities using county-wide cattle data in Washington may inflate risk since 1) county-wide cattle abundance likely overestimates cattle availability to wolves and 2) training models were based on fewer average cattle per pack territory. Conversely, cattle abundance on allotments may underestimate cattle availability since cattle grazed on private farms within a wolf pack territory were excluded.
Fine-scale data of wolf movements and activity patterns from GPS locations (which were not available for these analyses) might have increased the predictive power of this study. Miller et al. (2015) found predictive accuracy for modeling cattle depredation risk by tigers Panthera tigris in central India increased with decreasing spatial resolution, as the attack and killing stages of the hunt occurred over just a few square meters. Although wolves are coursing predators and the stalk and chase stages of the hunt can occur over a few meters to several kilometers (Peterson and Ciucci 2003), the scale of this study broadly encompassed the entire pack territory, limiting insights into within-territory processes. In smaller-scale studies from Idaho and Montana, depredations by wolves were correlated with spatial proximity of cattle to wolf den and rendezvous sites (Oakleaf et al. 2003, Bradley and Pletscher 2005). Locations from collared wolves could be used to reduce the spatial scale of future analyses to wolf ‘core areas’ (<50% home ranges, usually including den and rendezvous sites) which overlap cattle to varying degrees throughout the summer grazing season (Oakleaf et al. 2003). Although spatial overlap between cattle and wolf core areas does not always lead to depredation (Wydeven et al. 2004, Brown 2016), increased encounter rates between wolves and cattle increase depredation risk.
The prey CPUE variable produced inconclusive results in this study. We merged elk and deer abundance due to correlation (rs > 0.7) and were unable to include a metric of moose abundance since population estimates were unavailable for the entire study period and harvest was primarily through permitted hunts. In addition, hunter effort can be influenced by factors such as weather conditions and land access (Hewitt 2011), biasing estimates of ungulate catchper-unit effort. Since individual wolf packs specialize on different prey species (Mech and Peterson 2003, Spence 2017) and prey species composition varies by season and abundance (Kunkel and Pletscher 2001, Smith et al. 2010) incorporating accurate prey density estimates by species would be preferred, although not logistically feasible in many cases due to budget constraints or availability of historic data.
Adaptive management and future improvements
Adaptive management allows researchers and managers to develop and improve management practices through an iterative learning process (Rehme et al. 2011). Washington Department of Fish and Wildlife has successfully achieved wildlife management objectives for multiple species using this approach (Beausoleil et al. 2013, Maletzke et al. 2016). Our risk maps provide localized areas where WDFW can focus preventative methods to reduce depredation risk by wolves and federal agencies and livestock producers can select grazing areas. We suggest using the relative countywide cattle abundance risk map (Fig. 3A) to implement wolf–livestock education programs and depredation prevention tools in the South Cascades and Northwest Coast Recovery Region, and using the allotment cattle abundance risk map (Fig. 3B and risk maps in Hanley 2017) to do the same in the North Cascades and Eastern Washington Recovery Regions.
As wolves continue to recolonize Washington, priority areas for cattle depredation risk will change and expand. Secondly, applying interventions such as shifting livestock to low-risk grazing areas could alter the probability of cattle depredation (Miller 2015). Because this is a dynamic system, we developed an interactive web page using the Shiny package in RStudio (< https://shiny.rstudio.com>) to provide management agencies with a means of updating cattle depredation risk, given current conditions within a wolf pack territory, using the top models from this study (Supplementary material Appendix 1 Fig. A1).
In the future, fine-scale analyses could further guide adaptive management regimes in Washington and adjacent wolf-occupied states. Locations from GPS-collared wolves in conjunction with GPS locations, land ownership and allotment information (e.g. name, number of cattle, dates of use, animal husbandry practices) from confirmed and probable livestock depredations would provide a robust – both spatial and temporal – dataset. Documenting animal husbandry practices (e.g. type, intensity, and duration pre- and post-depredation) and prey abundance in wolf-occupied areas also would be useful.
Management implications
Predator–livestock interactions are complex, dynamic and challenging to predict. Variables included in risk models are spatially and temporally dependent, limiting inferences to other scales and locations. Addressing hypotheses at various scales and validating outcomes with new information from monitoring programs are essential to adequately mitigating human–carnivore conflict (Miller 2015, Williams and Brown 2016). Interactive web pages using programs such as Shiny by RStudio are a low-cost, effective means of informing stakeholders and policy-makers with up-to-date results. The risk maps produced in this study may be used by federal and state wildlife agencies as a baseline for determining high risk areas to implement depredation prevention measures and a template for future analyses as wolves continue to recolonize Washington. With patience, tolerance and commitment from all stakeholders it is possible to recover large carnivore populations and sufficiently protect livestock interests in the United States and worldwide.
Acknowledgements
We thank personnel from the government agencies and private companies who provided data for this research. C. Goldberg co-developed statistical methods and edited all drafts.
Funding – Funding for this study was provided by the Washington State Legislature and Washington Department of Fish and Wildlife.
References
Appendices
Supplementary material (available online as Appendix wlb-00419 at < www.wildlifebiology.org/appendix/wlb-00419>). Appendix 1.
