Intrinsic traits of woodland caribou Rangifer tarandus caribou calves depredated by black bears Ursus americanus and coyotes Canis latrans

Individuals in substandard physical condition are predicted to be more vulnerable to predation. Support for this prediction is inconsistent partly as a result of differences across systems in the life histories of predator and prey species. Our objective was to examine the physical condition of woodland caribou Rangifer tarandus caribou calves depredated by two predators with different life histories in Newfoundland, Canada. Black bears Ursus americanus are capable of chasing calves at high speeds over short distances and primarily prey on calves <1 month of age. Coyotes Canis latrans are cursorial predators that pursue prey over longer distances, which is expected to result in the selection of substandard individuals. We hypothesized that 1) black bears will kill calves in substandard physical condition, while 2) coyotes will kill calves from across the distribution of individual conditions. We used mitochondrial DNA species identification tests to assign predator species to calf mortalities. We then used molecular identifications and field observations to build a predictive model using generalized boosted trees to predict the predator species where a molecular identification was unavailable. We tested our hypotheses using Cox proportional hazards models under a competing risks framework. Bears killed younger calves and lighter calves, while coyotes killed heavier calves. Coyotes also killed more late-born calves, which might suggest prey switching as calves become more abundant later in the season. Our findings suggest that the physical constraints of predators play a greater role than predator hunting strategies in this system, but other processes are likely influential. The tendency for coyotes to kill heavier calves might result from sustained coyote predation over time, following the removal by black bears of lighter calves during their first month of age. This research illuminates the complexity of predator–prey interactions in Newfoundland and highlights an important source of variability for predator–prey systems.

For prey species, individuals in substandard physical condition are assumed to be more vulnerable to predation (Errington 1946) potentially resulting from behaviours that increase risk or via a lesser capacity to escape predators (Lima andDill 1990, Curio 1993). Although numerous studies demonstrate increased predation on individuals that are perceived to be in worse physical condition (Kunkel et al. 1999, Wright et al. 2006, Krumm et al. 2009), others fail to support this trend (Karanth andSundquist 1995, Anderson et al. 2007). The life history characteristics of prey and predator explain some of the variability underlying the distribution of depredated prey with regards to individual physical condition. For example, disproportionate predation on substandard individuals is theorized to be more likely for prey species that are difficult to catch (Temple 1987, Wirsing 2002. Predator hunting strategies, however, are also thought to play a role in determining which individuals are depredated from a prey population (Michalko and Pekár 2016).
Cursorial predators, such as most canids (family Canidae), are expected to pursue prey over long distances in more open habitats (Schaller 1972), which is predicted to result in the more frequent capture of substandard individuals, since these individuals are more likely to be identified during a long chase and fall behind their stronger conspecifics (Kunkel et al. 1999). In contrast, stalking and ambushing predators are not expected to demonstrate strong selection for substandard individuals, because chases are generally short (Caro and Fitzgibbon 1992). Similarly, incidental and opportunistic predators are predicted to kill prey from across the distribution of individual conditions. Incidental predation is the killing of an unexpected prey following a random encounter (Vickery et al. 1992). Incidental predation is opportunistic, but opportunistic predation also can include the purposeful selection for habitats where prey occur, followed by random encounters within those habitats (Vickery et al. 1992).
In Newfoundland, black bears Ursus americanus and coyotes Canis latrans are the primary predators of woodland caribou Rangifer tarandus caribou calves (Mumma et al. 2014, Bastille-Rousseau et al. 2016a). Both black bears (Kunkel et al. 1999, Zager andBeecham 2006) and coyotes (Bastille-Rousseau et al. 2016a) often prey heavily on neonates (<1 month of age) where ungulates are present, but coyotes are more capable of capturing older calves and fawns (Bishop et al. 2005). Recent research in Newfoundland suggests that bears might be actively hunting caribou calves in spring (Rayl et al. 2018), although black bear predation has been shown to be opportunistic in other locales (Bastille-Rousseau et al. 2011). Whether opportunistic or more directed (stalking), bears are only capable of high speed chases over short distances (Kolenosky and Strathearn 1987). In contrast, coyotes pursue prey over longer distances, although they do show some flexibility in hunting strategy dependent on habitat structure (Murray et al. 1995).
Our objective was to evaluate differences in the physical condition of depredated calves by black bears and coyotes. A common source of error when evaluating cause-specific mortality is the reliable identification of predator species at kill sites. To decrease this uncertainty, we used molecular techniques to identify the predator species at calf mortality locations (Mumma et al. 2014), and then built a predictive model using molecular predator identifications and field data to determine the predator species for mortalities where a molecular identification was unavailable.
We used this information and that of surviving calves in cumulative incidence functions to evaluate changes in predation risk from black bears and coyotes, conditional on calf age (0-60 days), and in a competing risks analysis to evaluate the probability of being killed by a black bear or coyote as a function of physical condition. Because black bears have difficulty capturing older calves (>1 month, Zager and Beecham 2006), we anticipated that they would also have difficulty capturing younger, stronger calves. We, therefore, hypothesized that 1) black bears would disproportionately prey upon substandard individuals despite their hunting strategy (short chases). In contrast, the cursorial hunting strategy (longer chases) of coyotes led us to predict that 2) coyotes also would select for substandard individuals. Through our examination of the physical condition of individual prey depredated by two predators with different life histories, we highlight an important of source of variability underlying many predator-prey systems.

Study area
The island of Newfoundland (111 390 km 2 ) is located off Canada's eastern coast and has a cool, maritime climate with interspersed coniferous forest, windswept barrens and peatlands (McManus and Wood 1991). Woodland caribou are widely distributed and the only native ungulate on Newfoundland. Since 1998, the Newfoundland caribou population has decreased by >66% (Bastille-Rousseau et al. 2013, Weir et al. 2014, largely as a result of increased calf predation driven by changes in climate and densitydependence (Bastille-Rousseau et al. 2016a, Mahoney et al. 2016. Historically, gray wolves, black bears and Canada lynx Lynx canadensis were considered the major predators of caribou calves in Newfoundland (Weir et al. 2014). Wolves, however, were extirpated in the early 1900s, and despite the recent discovery of several wolves in Newfoundland, are not thought to have re-established a breeding population (CBC 2012). In contrast, non-native coyotes are thought to have immigrated across sea ice and are now widely distributed across the island (Mumma et al. 2015). In addition to coyotes and black bears, current potential predators include lynx, red fox Vulpes vulpes and bald eagles Haliaeetus leucocephalus (Weir et al. 2014). We selected three study areas (La Poile, Middle Ridge and Northern Peninsula) that correspond to the calving grounds of four of Newfoundland's main caribou herds (La Poile, Middle Ridge and the adjacent Northern Peninsula and St. Anthony herds) (Fig. 1, Rayl et al. 2014).

Data collection
From late May though early June, 1-4 day old calves were located via helicopter and hand-captured in 2010-2012. Calves were aged based on hoof wear and affixed with very high frequency telemetry collars containing a motion sensitive transmitter. When field conditions permitted, calf sex, mass (kg) and hindfoot length (cm) was recorded. Calves were monitored every 1-3 days via helicopter and fixed-wing aircraft until August, when mortalities became rare. Mortality investigations were conducted on the day of detection, to reduce the probability of scavenging, using standard protocols (Supplementary material Appendix 1 Fig. A1) and included the sampling of carcasses for predator saliva using sterile, cotton swabs (Mumma et al. 2014).
Field observations included a categorical assessment of the state of the carcass (hereafter carcass treatment). Categories for carcass treatment included buried, decapitated, dismembered (limbs separated from body), intact, and sparse (few remains). The percentage of the carcass remaining (0-100) and the approximate distance to cover (m) were also recorded, along with whether or not the hide (skinning) and skull cap where separated from the carcass, and if there was visible throat trauma. Skinning, skull cap removal and throat trauma were categorized as yes, no or unknown, since a yes or no determination was not always possible when the majority of a carcass was consumed or when it appeared that a part of the carcass (e.g. head) and collar had been moved from the mortality site and buried. We also assigned an approximate age for each depredated calf using the estimated age at time of capture and mortality date.
Swabs were dipped in ethanol before swabbing up to four locations on each calf carcass or up to four tissues when the majority of the carcass was consumed. We collected multiple swabs from each location as back-up. All swabs were stored in paper envelopes placed in sealed plastic bags containing silica desiccant. Swabs were extracted in a laboratory dedicated to low quantity DNA samples and extracted using Qiagen DNeasy tissue kits. We used a mitochondrial DNA (mtDNA) control region fragment analysis method to identify all samples to species (De Barba et al. 2014). All failed samples were analyzed with an additional test to identify Canada lynx samples using mtDNA primers developed for Iberian lynx Lynx pardinus by Palomares et al. (2002). Additional details can be found in Mumma et al. (2014). We determined allele sizes using an Applied Biosystems 3130xl Genetic Analyzer and associated GeneMapper ver. 3.7 software.
Samples that failed fragment analysis testing were amplified and sequenced using mtDNA cytochrome B primers that amplify genetic material for most carnivores (Farrell et al. 2000) using conditions described by Onorato et al. (2006). These primers can identify black bear, Canada lynx and red fox, but not coyote. Any remaining failed samples were amplified and sequenced using the canid-specific mtDNA control region ScatID primers (Adams et al. 2003) to identify coyote samples that failed initial testing. We did not use molecular tools to test for the presence of bald eagles, because they frequently scavenge carcasses and are generally infrequent predators of caribou calves (O'Gara 1994).

Predictive model
We used predator identifications from swab samples (n = 69) and the corresponding field observations (carcass treatment, percentage consumed, distance to cover, skinning, skull cap removed, throat trauma and approximate age) in generalized boosted tree models (Freund and Schapire 1997, Friedman et al. 2000, Friedman 2001) to predict the predator species at mortality locations where the molecular testing failed (n = 43). We selected generalized boosted tree models because of their predictive power (Elith et al. 2006). Boosting is a step-wise process, first building an initial tree and then sequentially building trees on the residuals from the prior model, each time using a randomly selected proportion (bag fraction) of the observations and selecting a limited number of the most informative covariates. We set the bag fraction to 0.5 (Elith et al. 2008) and the contribution of each tree to 0.001 (learning rate, De'ath 2007). We modelled error using a Bernoulli distribution for our binary response variable attributing mortality locations to a black bear or coyote. We used five-fold cross validation to avoid overfitting and evaluated model fit using the area under the receiver operating curve (Mason and Graham, 2002) for the cross-validated data (CV AUC). We also calculated the relative influence (RI) of each predictor. RI (%) is a measure based on the improvement contributed to a model each time a predictor is selected, and then averaged across all models (Friedman 2001). Analyses were implemented in program R (<www.rproject.org>) using the packages gbm (Ridgeway 2017) and dismo (Hijmans et al. 2017).

Cause-specific mortality analyses
As a first step to exploring differences between calf mortality from black bears and coyotes, we built cumulative incidence functions under a competing risks framework (Fine and Gray 1999). This allowed us to visually assess differences in the cumulative probability (Dignam et al. 2012) of being killed by a bear or coyote as a function of calf age. Given that age was a covariate in our predictive model, we built two sets of cumulative incidence functions to evaluate the robustness of our results. First, we built cumulative incidence functions using only the molecular predator identifications, and then built a second set of functions using both the molecular and predicted predator identifications.
We then used data augmentation to evaluate the effect of intrinsic calf traits on cause-specific hazards (Lunn and McNeil 1995) in order to test our hypotheses regarding the physical condition of depredated calves. Data augmentation creates a replicate of the dataset corresponding to each outcome. For each replicate, the corresponding outcome is preserved while alternative outcomes are censored. The augmented dataset was modelled using Cox proportional hazards models (Cox 1972). As with other survival models, Cox proportional hazards models are concerned with modelling the time until mortality (outcome) or survival time T. Specific to Cox models is the underlying assumption of proportionality to the baseline hazard function h 0 (t) -with no assumptions in regards to its shape -that increases or decreases for an individual i dependent on a vector of covariates (Murray et al. 2010). Under a competing risks framework, the effects (β k1 , β k2 , …, β kp ) of each covariate (x 1 , x 2 , …, x p ) are estimated for each outcome K. We used calf mass at time of capture, hindfoot length and mass/hindfoot length as indices of physical condition. Calf mass has been associated with lower survival for neonates in ungulates (Couturier et al. 2009). We assumed that calf mass and hindfoot length at time of capture corresponded to mass and size at time of birth. We calculated mass/hindfoot length to account for the possibility that a large, lean calf might be in worse physical condition than a small, stout calf. Timing of ungulate parturition is theorized to match spring greenup, but might also serve to overwhelm predator populations (predator swamping), and thereby reduce risk (Ims 1990). In Newfoundland, most calves are born within a 12-day window around 30 May (Bergerud 1975). To account for potential changes in risk through time, we estimated calf birth date (date of capture -estimated age). We also anticipated that male calves might be larger than female calves and incorporated sex to account for this additional source of variability. We then evaluated our three condition indices in relation to birth date (t-tests of Pearson correlation coefficients, α = 0.05) and sex (t-tests, α = 0.05) to evaluate potential relationships.
We used our condition indices, sex and birth date to build competing, cause-specific, Cox-proportional hazards models and selected the best supported models using Akaike's information criteria (Akaike 1998) for small sample sizes (AIC c , Burnham and Anderson 2002). Our outcomes included mortality from black bears, coyotes and other causes (or unknown). We stratified each outcome to allow different baseline hazard functions (Kleinbaum and Klein 2012). Surviving individuals were right-censored at 60 days of age. We included a random intercept for study area and year to account for spatial and temporal variability. To aid with model convergence, we multiplied mass/ hindfoot length by 10. We only included a single condition index (mass, hindfoot length or mass/hindfoot length) in a given model to avoid multicollinearity and included interaction terms to account for potential relationships between calf sex and our condition indices. We also included a quadratic term for birth date to evaluate the potential for a non-linear relationship between birth date and the probability of mortality. We evaluated the assumption of proportionality using metrics based on Schoenfeld residuals (Fox 2002) and determined pseudo R 2 (Cox and Snell 1989) values corrected for censored data (O'Quigley et al. 2005) for our best supported models. We also assessed multicollinearity using the variance inflation factor (Graham 2003). All analyses were conducted in program R using the packages Hmisc (Harrell 2018a), survival (Therneau 2018b), coxme (Therneau 2018a), rms (Harrell 2018b) and MuMIn (Bartón 2018

Results
We hand-captured 328 calves in 2010, 2011 and 2012. From late May through July of each year, we investigated 129 mortality locations of collared calves. Collars containing blood and bite marks were the only evidence discovered at 17 of these locations. Only 105 out of the remaining 112 mortalities were swabbed because of a shortfall in swab collection supplies in 2010. We identified a single predator species at 67.6% (35 black bear, 34 coyote and 2 red fox) and multiple predator species at 5.7% (n = 6) of the mortality locations where a swab was collected. Molecular identification failed to identify any predator species at the remaining 26.7% (n = 28) of mortalities (Supplementary material Appendix 1 Table A1).

Predicting predators
We disregarded mortalities where red fox DNA was discovered, since they were rare (n = 2) and suspected to result from scavenging, and built our predictive model using the 69 mortalities attributed to black bear or coyote. We limited each tree to two splits (predictors) given our sample size (Elith et al. 2008). Cross-validation identified 4150 trees as the optimal number of models to maximize prediction and avoid overfitting. Our final model was highly predictive (CV AUC = 0.903).
The most informative variable in our model was age (RI = 25.9%), followed by skull cap removal (RI = 24.7%), carcass treatment (RI = 23.6%), distance to cover (RI = 9.9%), throat trauma (RI = 8.7%), percent consumed (RI = 6.9%), and hide removal (skinning, RI = 0.4%) (Fig. 2). Bear mortalities were characterized by younger calves, sparse remains or intact carcasses, skull cap removal, close to cover, no or unknown throat trauma, and higher consumption (Fig. 2). Coyote kills were characterized by older calves, buried, decapitated or dismembered carcasses, no skull cap removal, farther from cover, throat trauma and lower consumption (Fig. 2). There was, however, a minimal difference in the model between black bears and coyotes at the highest levels of consumption (>90%), because a high percentage of consumption was recorded for coyote kills when carcasses were decapitated and the head and collar were moved, even though the actual kill site location and percent consumption was unknown. Hide removal (skinning) was an uninformative variable with low influence (Fig. 2) and little discernment between black bears and coyotes, which was related to two different skinning techniques. At some black bear mortalities, the hide was entirely removed from the carcass, while coyotes tended to leave the hide attached, but skinned the hide away from the body and down the limbs. The final model was used to predict the predator species (20 black bear and 23 coyote) at the 43 remaining kill sites with an unknown predator (Supplementary material Appendix 1  Table A1).

Cause-specific mortality
We excluded 70 calves from our cause-specific mortality analyses, because they were collared in an area undergoing coyote removals (n = 37) (Lewis et al. 2017), were missing sex, mass or hindfoot length information (n = 31), or died from capture-related causes (n = 2). Five additional calves were excluded because they died prior to verifying that they had re-bonded with the female after collaring and one calf was excluded because of collar failure. One-hundred and forty-three of the remaining 247 calves were alive at the end of their first 60 days and 43, 44 and 16 of the deceased calves were attributed to black bears, coyotes or 'other causes', respectively. 'Other causes' of calf mortality included 13 calves where only a bloody or chewed collar was found, two calves that died of natural causes, and one calf that drowned (Supplementary material Appendix 1 Table A2).
Cumulative incidence functions demonstrated a high probability of mortality for caribou calves within the first two weeks of life (Fig. 3). The highest level of early (<2 weeks of age) mortality was attributed to black bears, after which the probability of being killed by a bear declined and appeared to stabilize (Fig. 3). The probability of being killed by a coyote was also high within the first two weeks of life, before declining, but demonstrated a secondary increase between 40 and 55 days. The shape of the cumulative incidence functions for bears and coyotes were similar when modelled using molecular and predicted predator species identifications (Fig. 3) versus the molecular identifications alone (Supplementary material Appendix 1 Fig. A2). The cumulative risk of bear and coyote predation declined when only molecular identifications were used as a result of the higher number of mortalities attributed to unknown (other) causes ( Fig. 3 and Supplementary material Appendix 1 Fig. A2).
Several models generated support in our competing risks analysis (Table 1) (∆AIC c ≤ 2, Burnham and Anderson 2002). Our best supported model included mass (pseudo R 2 = 0.320), but a model including mass, sex and birthdate (∆AIC c = 0.19, pseudo R 2 = 0.396) was also supported, as were models including the mass/hindfoot length, sex and birth date (∆AIC c = 1.49, pseudo R 2 = 0.389) and mass/ hindfoot length alone (∆AIC c = 1.43, pseudo R 2 = 0.389). Calves with lesser mass and mass/hindfoot length were more frequently killed by black bears as indicated by hazard ratios and corresponding upper 95% confidence bounds that were <1 (Fig. 5). Calves with greater mass and mass/hindfoot length appeared to be more frequently killed by coyotes as indicated by hazard ratios >1, although lower 95% confidence bounds slightly overlapped 1 (Fig. 5). There did not appear to be any influence of sex or birth date on predation by black bears or sex on predation by coyotes (Fig. 5). Coyotes did, however, appear more likely to kill late-born calves (Fig. 5).

Discussion
Black bear predation was higher for younger calves, and in alignment with our first hypothesis, bears were also more likely to kill lighter calves and calves with a lower mass/ hindfoot length. Coyotes, despite their cursorial hunting strategy, did not disproportionately kill calves in substandard physical condition. Previous research demonstrated that bear Ursus spp. predation on neonate ungulates declines with calf age (Adams et al. 1995, Singer et al. 1997, likely because of difficulty capturing older and faster calves and fawns. Our results indicate that this constraint also limits predation by black bears on heavier, younger calves. In contrast, these potentially less-vulnerable calves did not experience a reprieve from coyote predation. Temple (1987) suggested that selection for substandard individuals is predicated upon the difficulty a predator has in catching prey. Although black bears are capable of running at high speeds, they lack the endurance of coyotes to pursue prey over longer distances. In Newfoundland, one explanation for the patterns observed is that the physical abilities of bears and coyotes have a greater impact on prey selection than hunting strategies and that capturing calves proves a more difficult task for black bears than coyotes. Indeed, our results suggest a greater probability of being killed by a coyote for calves with a higher mass and mass/ hindfoot length. This potentially results from competition with black bears. High predation by bears on younger, lighter calves might skew the distribution of surviving calves towards heavier individuals. These older calves might remain susceptible to coyote predation, but would be rarely killed by black bears, thus resulting in coyote-killed calves being in better physical condition on average, than bear-killed calves. Other explanations are also possible, such as differences in female caribou behaviour and space-use, dependent on maternal condition (and correlated with calf condition, Singer et al. 1997), which might increase co-occurrence and risk from one predator, while decreasing risk from the other (Duquette et al. 2014). For example, if female caribou in poor condition select for areas with greater forage that are more frequently used by black bears, it could increase risk from bear predation.
A more definitive explanation (temporal changes in spaceuse) is available for the secondary spike in coyote predation on older calves. Research in British Columbia demonstrated increased caribou calf mortality from wolves, coincident with females leaving their calving grounds (Gustine et al. 2006). In Newfoundland, females in the La Poile and Middle Ridge study areas move away from their calving grounds in mid-July (Rayl et al. 2014). This migration and concurrent increase in risk from coyotes, likely results from greater habitat overlap with coyotes and potentially increased calf detection during migration (Bastille-Rousseau et al. 2016b).
Apart from age, we observed a negative relationship between estimated calf birth date and risk from coyotes, but observed no discernable trend for bears. This might suggest that coyotes key-in on calves (prey switching) once calves become readily abundant across the landscape. Increased predation for late-born caribou and elk calves by predators has been demonstrated in other carnivore-ungulate systems (Adams et al. 1995, Singer et al. 1997), but not consistently across studies (Capreolus capreolus - Gaillard et al. 1993, Antilocapra americana -Jacques et al. 2015. Differences in the life histories of the predator species in these systems   might explain the inconsistent findings between birth date and mortality. Alternative mechanisms have been proposed to explain the lower survival of late-born calves, such as correlations between female or calf condition and birthing date (Fairbanks 1993) and birthing asynchrony with spring green-up (Plard et al. 2015); adult females of late-born calves might be limited by forage quality, thereby inhibiting their ability to provide nourishment and slowing calf growth. We did not have information on maternal condition, but did not detect a relationship between calf condition and birth date. Further, if late-born calves were in worse condition (lesser mass), we would anticipate bears to be selecting for late-born calves, since they demonstrated a propensity for killing lighter individuals.
Our research demonstrates the novel amalgamation of molecular tools and a modelling approach to reduce the uncertainty of predator identifications at prey mortality locations. Our model was highly predictive and revealed key differences in observations at calf mortalities attributed to bears (younger, skull cap removal, high or no consumption, closer to cover) versus coyotes (older, moderate consumption, further from cover, throat trauma). These observations would likely be informative for other bear-coyote-caribou systems, but more broadly, we would recommend the application of molecular and modelling approaches to reduce the uncertainty of predator identification in other multi-predator systems.

Conclusion
The physical constraints of black bears is a probable explanation for the tendency of bears to kill calves in worse physical condition. In contrast, the speed and endurance of coyotes would seem to allow them unfettered access across the distribution of available calves, thus suggesting that hunting strategy is of secondary importance, as it pertains to the selection of substandard individuals, in this system. Other factors, however, such as competition between bears and coyotes for calves and changes in predator and prey space-use, likely contribute additional variability. This study highlights the utility of molecular tools and a modelling approach to reduce the uncertainty of predator identifications at prey mortality sites and provides insights into mechanisms that determine who is depredated by whom.