Identification of carnivore kill sites is improved by verified accelerometer data

Quantifying kill rates is central to understanding predation ecology. However, estimating kill rates and prey composition in carnivore diets is challenging due to their low densities and cryptic behaviors limiting direct observations, especially when the prey is small (i.e., < 5 kg). Global positioning system (GPS) collars and use of collar-mounted activity sensors linked with GPS data can provide insights into animal movements, behavior, and activity. We verified activity thresholds for American black bears (Ursus americanus), a bobcat (Lynx rufus), and wolves (Canis spp.) with GPS collars containing on-board accelerometers by visual observations of captive individuals’ behavior. We applied these activity threshold values to GPS location and accelerometer data from free-ranging carnivores at locations identified by a GPS cluster algorithm which we visited and described as kill sites or non-kill sites. We then assessed use of GPS, landscape, and activity data in a predictive model for improving detection of kill sites for free-ranging black bears, bobcats, coyotes (C. latrans), and wolves using logistic regression during May–August 2013–2015. Accelerometer values differed between active and inactive states for black bears (P < 0.01), the bobcat (P < 0.01), and wolves (P < 0.01). Top-performing models of kill site identification for each carnivore species included activity data which improved correct assignment of kill sites by 5–38% above models that did not include activity. Though inclusion of activity data improved model performance, predictive power was less than 45% for all species. Collar-mounted accelerometers can improve identification of predation sites for some carnivores as compared to use of GPS and landscape informed covariates alone and increase our understanding of predator–prey relations.


Background
Kill rate is the central process of predation and estimation of kill rates, and prey composition in carnivore diets, is important for understanding the effects of predators on their prey [1]. However, due to typical low densities and often cryptic behavior, estimating carnivore diet and predation is challenging. Coupling global positioning system (GPS) technology with field investigations can provide estimates of kill rates to address effects of carnivores on prey populations [2][3][4] without the need for direct observation.
Monitoring animal movements from GPS collars is common in studies of carnivores and other species difficult to observe [5][6][7][8], including use of clusters of GPS locations to detect potential predation sites [9][10][11]. Unfortunately, discerning among behaviors using GPS location data alone can be difficult for species that are central place foragers (e.g., ring-billed gulls [Larus delawarensis]; [12]) or species that forage and sleep in close proximity (e.g., American black bear [Ursus americanus]; [13]). In addition, challenges exist in identifying predation sites for species predating on small-and medium-sized prey when handling time is short and prey remains are limited [11,14]. While GPS collars can track spatial movements, obtaining reliable information regarding animal activity can be problematic [5,15,16]. Integrating landscape characteristics with behavior can improve identification of activities at location clusters [17], but multiple behaviors may occur at each cluster and relating behavior to activity may prove more informative. The ability to differentiate between inactive and active states when paired with GPS data could be used to discern rest sites from areas used while active (e.g., foraging [18]).
Accelerometers can lend insight into animal behavior by measuring their movement at short intervals over long periods of time [19]. Accelerometers have been used to monitor and verify activity in mammals including moose (Alces alces [20]), elk (Cervus elaphus [21,22]), bears (U. maritimus and U. arctos [23]), brown hares (Lepus europaeus [24]), African elephants (Loxodonta africana [25]), Japanese macaques (Macaca fuscata [26]), and dingoes (Canis dingo [27]). Identifying active and inactive periods within areas of intensive use can provide information on behaviors associated with habitat use (e.g., foraging, resting, or traveling) and clarify benefits of that area [28]. In addition, knowing when animals are active or inactive is important when investigating travel efficiency or energy expenditure [29] or behaviors such as foraging or predation events [19,30]. Thus, pairing activity data with GPS movement data can improve our understanding of movements, resource selection, and inactivity [31] as it relates to carnivore predation ecology. Linking activity data with GPS data has improved detection of carnivore predation sites, most notably with cougar (Puma concolor) where prey items are typically large and prey remains at clusters are identifiable [30,32].
We build on the GPS cluster model to predict kill sites by incorporating activity data collected by accelerometers on GPS collared American black bears, bobcats (Lynx rufus), coyotes (Canis latrans), and wolves (C. spp.) to identify kill sites during May-August when prey are typically of small or medium size [11,14,33]. We validated accelerometer values to discern active from inactive behaviors then incorporated activity data into a competing model framework to assess performance of cluster characterization as kill sites or non-kill sites. We hypothesized that measured activity from accelerometers could be used to discern active from inactive behaviors and improve detection of kill sites by defining behavior state at each cluster. We predicted clusters with greater numbers of active locations and greater proportion of time spent active would be positively associated with the likelihood of a cluster being a kill site.

Identifying activity thresholds
We observed captive black bears during August-September 2011 and a captive bobcat and wolves during August-September 2013 at the DeYoung Family Zoo, Menominee County, Upper Peninsula of Michigan, USA. Carnivores were housed in open-air enclosures and fed twice daily at prescheduled times by zoo personnel. Three male and 2 female black bears shared a 0.3-ha enclosure, 8 male and 5 female wolves shared a 0.2-ha enclosure, and the male bobcat occupied a 36-m 2 2-story enclosure by itself. Zoo personnel conducted immobilizations in conjunction with each individual's biannual health check and fitted 3 black bears (2 male, 1 female), 1 male bobcat, and 2 wolves (1 male, 1 female) with a Lotek 7000 GPS collar (Lotek Wireless, Newmarket, Ontario, Canada) equipped with an on-board tri-axial accelerometer. In the default pre-programmed mode, collar-mounted accelerometers measured gravitational acceleration along two axes (X [anterior-to-posterior] and Y [side-to-side]) and did not report measurements about the third axis (i.e., Z [up-and-down]). Accelerometers measured activity four times per second simultaneously on each axis, where the difference in acceleration between two consecutive measurements was given a value from 0 to 255. Collars averaged and stored activity data over the default pre-programmed sampling interval of 300 s with associated date and time. As accelerometers functioned on 8-s loops, the true measure of mean activity was 296 s. We chose this interval to match free-ranging collared carnivores where collars had storage limitations for activity data due to deployment length.
We had 4 observers directly recorded black bear, bobcat, and wolf behaviors (e.g., inactive, walking, running) during observation segments conducted in 300-s intervals (up to 4 h in duration) and recorded the time to the nearest second at each behavior change. We did not use a double observer study design to verify categorized behavior groups between individuals. We initially categorized behavior as resting (i.e., inactive), loafing (e.g., minor head movement while lying, sitting, or standing), feeding/drinking, walking/trotting/running (i.e., locomotion), and other movements (e.g., play, dominance interactions, etc.). We summed the duration of each behavior observed for each 300-s interval and assigned the dominant behavior to that interval. As there were no accelerometer values provided for the Z axis, we could not directly replicate common metrics such as overall dynamic body acceleration [34] and instead used the sum of the accelerometer X and Y axis values recorded during each interval to quantify activity and used this as the sampling unit to test for differences among behaviors. Although an analysis of variance [35] identified differences for some behaviors, there was significant overlap within inactive and active behaviors when we used the Tukey-Kramer method [36] for unequal sample sizes as a post-hoc test (see Additional file 1) so we collapsed behaviors into inactive (i.e., resting or loafing) and active (all other behaviors) states. We used Welch's t test for unequal variances to identify differences between active and inactive states' means with α = 0.05. We used the upper confidence limits of the means for loafing behaviors (i.e., head movement only) as determined by accelerometer data from captive carnivores as the breaking value between active and inactive states. Though we did not collar captive coyotes we used the wolf activity threshold for coyotes as locomotion and morphology are similar for these canids, though body mass of wolves is 2-3 times greater than that of coyotes [37,38]. Mississippi State University's Institutional Animal Care and Use Committee approved collaring and observation of captive carnivores (protocols 09-004, 12-012).

Field study site
We conducted our field study in portions of Iron and Baraga counties in Michigan's Upper Peninsula, USA (46.27ºN, 88.23ºW) spanning about 1000 km 2 . Forested lands (86%) dominated the landscape and consisted of northern hardwood and boreal forest types with wetlands and limited agricultural lands interspersed [39]. May-August temperatures ranged from average highs of 25.8 °C during July to lows of 1.7 °C during May. Average rainfall during May-August was 36.7 cm [40]. Carnivore densities within the study area were 25.9/100 km 2 for black bears, 3.8/100 km 2 for bobcats, 23.8/100 km 2 for coyotes, and 2.8/100 km 2 for wolves and density for white-tailed deer adults and fawns (Odocoileous virginianus), a dominant prey during the study period, was 571/100 km 2 during 2013-2015 [41].

Cluster investigations
During February-June of 2013-2015, we captured freeranging black bears, bobcats, coyotes, and wolves using culvert traps, #3 Victor Soft-Catch padded foothold traps (Oneida Victor, Euclid, Ohio, USA) or modified cage traps, #3 Victor Soft-Catch padded foothold traps, and modified MB 750 foothold traps (Minnesota Brand, Pennock, Minnesota, USA), respectively. We used an intramuscular injection of ketamine hydrochloride (Ketaset ® , Fort Dodge Laboratories, Inc., Fort Dodge, Iowa, USA) and xylazine hydrochloride (X-Ject E ™ , Butler Schein Animal Health, Dublin, Ohio, USA) to immobilize bobcats (10:1.5 mg/kg), coyotes (4:2 mg/kg), and wolves (10:2 mg/kg) and used an intramuscular injection of tiletimine/zolazopam (Telazol ® ) to immobilize black bears (7 mg/kg [42]). We fitted each carnivore with a Lotek 7000MU (black bears) or 7000SU (bobcats, coyotes, and wolves) GPS collar (Lotek Wireless, Newmarket, Ontario, Canada) programmed to attempt a GPS-fix every 15 min and collect activity data averaged and stored over the default pre-programmed interval of 300 s. We administered yohimbine hydrochloride (0.15 mg/kg; Hospira©, Forest Lake, Illinois, USA) to antagonize the xylazine hydrochloride before we released bobcats, coyotes, and wolves at the capture site to facilitate recovery. Mississippi State University's Institutional Animal Care and Use Committee approved all capture and handling procedures (protocols 09-004, 12-012). We downloaded GPS data using ultra-high frequency communication with a handheld unit twice weekly by airplane.
We identified potential kill sites using a rule-based algorithm [11] developed in R (version 3.0.2; R Core Team 2018) that delineated clusters of ≥ 4 GPS locations within 50 m of each other and occurring within 24 h, then calculated a geometric mean for each cluster center. Two personnel, with or without a trained detection dog, investigated cluster sites within 10 days of cluster formation to increase detection of kills [11]. Detection dogs were trained to locate hair, feathers, bones, and other prey remains. In all cases we searched the area systematically throughout the 50 m radius that included the locations used to identify the cluster. At each cluster we recorded carnivore tracks, scat, prey remains, and other animal sign to determine site use (i.e., kill or non-kill site). We confirmed presence of a predation (i.e., kill site) by detection of prey remains with the presence of fresh blood, minor decomposition of carcass remains, and hemorrhaging on the hide or muscle tissue of remains [9,10]. If prey remains were present but lacked evidence of predation we considered the cluster a non-kill site. Where we identified bed or den sites at clusters with no evidence of kill site, we considered these non-kill sites. We confirmed den sites where pups, kittens, or dens were present. We confirmed a bed site when matted vegetation or soil with hair was present to confirm species use. If we identified evidence of a kill site and a bed site at a cluster, we recorded the site as a kill site. We did not investigate clusters formed within 5 days of predator capture to reduce potential effects of their recovery from capture and handling.

Modeling kill sites
We coded clusters as kill sites (1) or non-kill sites (0) as the response variable for logistic regression models. For each cluster we recorded covariates derived from GPS data (number of GPS locations, time of day of formation, number of days at cluster) and landscape data (distance to nearest road, hydrologic feature, and edge, and patch size) thought to improve identification of kill sites [10,11]. We used the number of 24-h periods at a cluster as the number of days and the total number of GPS locations the algorithm used to create each cluster as the number of GPS locations. We calculated time of day as a binary variable with 0 forming at night (between local sunset and sunrise) and 1 forming during the day (between local sunrise and sunset). We calculated distance to nearest road (Michigan Geographic Framework, all roads v17a), hydrologic feature (Michigan Geographic Framework, hydrography lines v17a), and nearest edge (National Land Cover Database [39]) from the geometric mean of the cluster in ArcMap 10.3 (Environmental Systems Research Institute, Redmond, California, USA). We calculated patch size as the area of the land cover classification [39] that contained the geometric mean of the cluster in Geospatial Modeling Environment [43].
We also identified covariates of activity we predicted would improve identification of kill sites. We matched each 15-min GPS relocation in the cluster with the overlapping 300-s accelerometer interval and used that value to determine the activity state of that individual using the captive carnivore active and inactive state thresholds. We did not use all 300-s intervals that occurred throughout the duration of the cluster as individuals would leave the cluster and return throughout the start and end periods as defined by the algorithm (i.e., ≥ 4 locations within 50 m within a 24-h period). We counted the number of active 300-s intervals within each cluster (NACT) and calculated proportion of time active (i.e., number of active intervals/total number of intervals; PACT) within each cluster as covariates of cluster activity. We also summed the accelerometer values for all 300-s intervals that overlapped with GPS relocations within the cluster for the sum of cluster activity (SACT) and calculated the mean of those overlapping 300-s interval accelerometer values for the mean of cluster activity (MACT). We included the quadratic term for MACT and SACT as individuals may rest where kills occur, and we therefore did not expect the relationship of activity and kill sites to be linear. For example, if kill clusters are associated with resting at the same site then kill probability would peak at middle values for MACT or SACT where lesser values would be non-kill sites where only resting or loafing occurred and greater values would be non-kill sites where only play or continuous foraging (e.g., black bears foraging on herbaceous vegetation) occurred.
We built three base models organized by the covariate data type used to predict kill sites (i.e., GPS model, landscape model, and activity model). We excluded any covariates that were highly correlated (|r| > 0.8) to reduce model complexity. The GPS model included number of cluster locations, number of days cluster was visited, and time of day the cluster was initiated. The landscape model included distance to nearest edge, hydrologic feature, and road, and patch size at the geometric mean of the cluster. The activity model included the SACT, MACT, and the quadratic term for SACT and MACT. We used these three base models to construct seven models a priori as our model set for each species which included all combinations of the GPS, landscape, and activity models. We used a likelihood ratio test to assess the use of a global fixed-or mixed-effect generalized linear model with a random effect for individual [17,44]. As these models did not differ for any species (black bear, χ 2 = 3.30, P = 0.07; bobcat, χ 2 = 0.89, P = 0.35; coyote, χ 2 = 0.00, P = 1.00; wolf, χ 2 = 0.55, P = 0.50), we used the simpler model structure of the fixed-effect generalized linear model with the "logit" link for all analyses using package 'stats' within program R [45]. We ran all model set combinations and identified the top-ranked model for each species with Akaike Information Criterion for small sample sizes (AIC c ) and calculated Akaike weights to measure model support and aid in model selection [46]. We considered the model with the lowest AIC c score as the most supported model and used that model in assessing performance. As models within 2 AIC c of the top model are often considered equally supported [46], we report model-averaged parameter estimates where applicable to facilitate predictive ability. We calculated the probability output (i.e., the cut-off value) for each species top-ranked model as the minimum difference among classification rate (proportion of correctly classified kill and non-kill sites), sensitivity (proportion of true kill sites correctly classified), and specificity (proportion of true non-kill sites correctly classified) to optimize the correct assignment of kill sites for all classification metrics. We validated model accuracy by bootstrapping with 999 iterations of the top model for each species and randomly selecting 70% of the data for each run and withheld 30% of the data to test the model upon and report the means and 95% confidence intervals for the iterations.

Identification of behavior groups
We directly observed and recorded the behaviors of the 3 captive black bears for 4835 min ( x = 1578.0, SD = 1622.0) on 21 occasions, the bobcat for 1440 min on 7 occasions, and 2 wolves for 5245 min ( x = 2582.5, SD = 272.2) on 26 occasions. Black bears, the bobcat, and wolves were inactive 45.1%, 68.2%, and 72.4% of the time, respectively, and active 54.8%, 29.3%, and 18.9% of the time, respectively. Active and inactive states were differentiated by accelerometer values (Fig. 1) for black bears (active x = 91.5, inactive x = 22.9, P < 0.01), the bobcat (active x = 73.1, inactive x = 17.5, P < 0.01), and wolves (active x = 83.3, inactive x = 12.7, P < 0.01). The breaking values for active and inactive behaviors were 35.9 for black bears, 36.8 for bobcats, and 30.7 for coyotes and wolves (see Additional file 2).
We investigated 2591 clusters with an average of 4.6 days following formation (range 1-11 days) with (n = 1143) or without (n = 1448) detection dogs. Though detection dogs were generally better at finding kill sites (Additional file 3) we pooled cluster investigations as the trends of cluster determination were similar and the number of known kill sites in the sample was small. Of the investigated clusters 2110 had all necessary covariates for analysis of visited black bear (n = 755), bobcat (n = 350), coyote (n = 549), and wolf (n = 456) clusters. Average number of locations at black bear, bobcat, coyote, and wolf clusters was 21.4 (SD = 24.7, range = 5-300), 20.6 (SD = 24.7, range = 5-246), 21.3 (SD = 16.58, range = 5-110), and 26.9 (SD = 23.3, range = 5-247), respectively. We found kill sites at 2.1% of black bear (n = 16), 16.9% of bobcat (n = 59), 8.0% of coyote (n = 44), and 11.0% of wolf (n = 50) clusters. Prey types were dominated by white-tailed deer fawns (34% of kill sites) as well as other small-bodied species including snowshoe hares (Lepus americanus), beavers (Castor canadensis), muskrats (Ondontra zibethicus), red squirrels (Tamiasciurus hudsonicus), raccoons (Procyon lotor), ruffed grouse (Bonasa umbellus), sandhill crane colts (Grus canadensis), turkeys (Meleagris gallopavo), Canada geese (Branta canadensis), and various songbirds. Adult or yearling white-tailed deer were the largest bodied prey, found only at wolf clusters, and represented 16% of kill sites. We excluded NACT and PACT as covariates for activity data as these were highly correlated (|r| ≥ 0.8) with SACT and MACT of clusters, respectively. The top supported models for each species included activity data ( Table 1) and each had one equally supported model for which we calculated model-averaged parameter estimates ( Table 2). The top supported model for bobcats, coyotes, and wolves included MACT as a significant term. In addition, MACT 2 and LOCS were significant terms in the bobcat and wolf top models. Kill site probability was greater for ROAD in the top coyote model and HYDRO in the top wolf model. No model terms were significant for black bear. Classification rate was above 53% for all carnivore models, the greatest for coyote clusters at 68% (Table 3). Predictive power (sensitivity) was relatively low across species (28-44%) but greatest for wolves. Top models were generally better at classifying non-kill clusters (56-72%; specificity) than kill clusters.

Discussion
We demonstrated substantial improvements in identification of kill sites in 4 species of carnivores by considering behavior state (active or inactive). Activity was included in the top model for each species and identification of

Table 1 Competing model rankings by Akaike Information Criterion for small sample sizes (AIC c ) for factors influencing prediction of carnivore kill clusters
'GPS' includes number of cluster locations, number of days cluster was visited, and time of day the cluster was initiated. 'Activity' includes the sum of cluster activity, mean cluster activity, and the quadratic term for mean cluster activity and sum of cluster activity. 'Landscape' includes distance to nearest edge, hydrologic feature, and road, and patch size at the geometric mean of the cluster. 'Global' model includes all previous covariate terms. clusters as kill or non-kill sites was improved 5-38% across species compared to uninformed GPS-only cluster models for prey of black bears, bobcats, coyotes, and wolves during summer. Activity data from on-board accelerometer sensors have been used to identify differences in mobile and immobile behaviors for species such

Table 2 Model averaged parameter estimates and 95% confidence intervals for prediction of kill sites from clusters of black bear, bobcats, coyotes, and wolves
Parameters include mean activity (MACT), sum of activity (SACT), number of global positioning system locations (LOCS), number of days spent at the cluster (DAYS), time of day the cluster was initiated (TOD), land cover patch size the cluster was within (PATCH), distance to nearest hydrologic feature from cluster center (HYDRO), distance to nearest land cover edge from cluster center (EDGE), and distance to nearest road from cluster center (ROAD as puma [30] and improved kill site detection by 10% from GPS data alone [32]. The accelerometers we used to collect activity data are commercially available alone or in combination with tracking devices, and are consistent in collection and inexpensive. Kill sites were best identified for black bear using activity and GPS covariates. As black bears often forage and sleep in the same locations [13], identifying clusters with greater proportions of active locations can help differentiate sleeping sites from mixed-use sites and reduce the need for field investigation of up to 64% of clusters. Though activity was included in the top model for bobcats and improved model performance over GPS covariates alone, the predictive power (34%) and overall classification rate (60%) was similar to Svoboda et al. [11] where only GPS-derived covariates were in their top model. This could be a result of a large percentage of bobcat clusters being classified as unknown (> 70.0%; see Additional file 3), where we did not find evidence of a kill site or resting site and may have misclassified predation events as non-kill sites which confounded results. Coyote classification of kill sites (28%) was less than for other species; which may be due to their omnivorous diet where foraging on fruits, insects, or small mammals [47] could result in undetected predation events by leaving no remains or their short duration. In contrast, the classification of non-kill sites (72%) for coyotes was the greatest for any species due to high levels of overlap between kill and non-kill cluster covariates which may be better differentiated with more refined behavioral classifications. The greater identification of wolf kill sites (44%) may be attributed to wolves being the only carnivore that killed large-bodied prey. Killing larger prey increases energy expenditure and handling time, making detection of kill sites easier as observed with other carnivores taking large prey [30,32].
We emphasized maximizing the overall classification of clusters for accuracy, sensitivity, and specificity. Predictive power can be improved by changing the cut-off value but at the cost of increasing the number of false-positives (i.e., incorrect prediction of kill sites), which would increase the number of incorrect sites investigated by researchers [11]. Although including activity data offers improved detection of predations over GPS and landscape covariates alone, further refinement is necessary as our predictive power was < 45% for all species, which points to the challenges of identifying kill sites when prey are small.
Though we distinguished only active from inactive behaviors, using visual observations of collared animals can improve our ability to detect differences in activity values for similar behaviors (e.g., capture of different prey [48]). Because individuals may perform several behaviors with a similar accelerometer response, applying and verifying accelerometer values to varying behavioral characteristics is important to improve interpretation of activity. The sampling interval we used for the accelerometer to define activity could have influenced our ability to accurately estimate behavior [30]; we recognize that activity data averaged across a 300-s sampling interval may be less informative when bursts of activity and periods of inactivity occur within the sampling interval. Examining activity through validation of several sampling intervals for analyses may improve the categorical assignment of behaviors identified beyond inactive and active states. In addition, as model accuracy and predictive ability can be further reduced for clusters where both kill sites and resting sites occur, we recommend consideration of multinomial response models that incorporate activity to further refine behavior at clusters. Alternatively, using activity thresholds to remove all inactive GPS locations before creating clusters may reduce problems with confounding influence of inactive locations. However, if inactive handling time is important to reach the defined cluster interval (i.e., number of locations) or the location interval is too great, some predations may be missed.

Table 3 Assessment of top model performance (selected by the lowest Akaike Information Criterion score) by classification rate (proportion of correctly classified kill and non-kill sites), sensitivity (proportion of true kill sites correctly classified), and specificity (proportion of true non-kill sites correctly classified)
The cut-off value was determined by the minimum distance between the three metrics for model assessment to maximize correct classification of kill and non-kill sites. Standard error (SE) for each metric is also displayed.