Age‐specific foraging strategies among pumas, and its implications for aiding ungulate populations through carnivore control

Humans have been controlling carnivore numbers for centuries. Predator hunting, however, may indirectly influence predator‐prey dynamics unintentionally by influencing the age‐ and sex‐structure of predator populations that exhibit intraspecific (IS) variation in prey selection. We tested for IS in a small population of pumas in the Greater Yellowstone Ecosystem, United States, and identified foraging strategies shared by multiple individuals. Further, we tested extrinsic and intrinsic variables that explained differences in foraging strategy. Our top model was composed of a single intrinsic characteristic, Age. In short, the older the animal, the larger the prey it specialized upon. Our provocative results suggest that the current controversial strategy of increasing puma culling to aid mule deer, as currently underway in Colorado, may in fact exacerbate problems for mule deer by changing the age‐structure of the puma population to predominantly younger animals that are more likely to hunt deer over elk.

Humans have been controlling carnivore numbers for centuries. Predator hunting, however, may indirectly influence predator-prey dynamics unintentionally by influencing the age-and sex-structure of predator populations that exhibit intraspecific (IS) variation in prey selection. We tested for IS in a small population of pumas in the Greater Yellowstone Ecosystem, United States, and identified foraging strategies shared by multiple individuals. Further, we tested extrinsic and intrinsic variables that explained differences in foraging strategy. Our top model was composed of a single intrinsic characteristic, Age. In short, the older the animal, the larger the prey it specialized upon. Our provocative results suggest that the current controversial strategy of increasing puma culling to aid mule deer, as currently underway in Colorado, may in fact exacerbate problems for mule deer by changing the age-structure of the puma population to predominantly younger animals that are more likely to hunt deer over elk.

K E Y W O R D S
carnivore hunting, intraspecific prey selection, ungulate management, wildlife management

| INTRODUCTION
Humans have been killing predators for centuries to augment ungulate populations and other game species we value (Reynolds & Tapper, 1996). In theory (e.g., Hairston, Smith, & Slobodkin, 1960), controlling predators should allow prey species like elk (Cervus canadensis) and deer (Odocoileus spp.) to grow in number, while simultaneously serving the dual purpose of removing our competitors for shared resources. Predator control, however, is also driven by special interest groups that provide fiscal donations to wildlife agencies and lobby for supportive policy and legislation that prioritizes specific game species (Hurley et al., 2011;Todd, 2002). For example, managing for larger prey populations allows for increased hunting opportunities.
Hunting license sales for ungulates, rather than carnivores, generate substantial revenue for local communities and the institutions managing wildlife for the public trust (Anderson, Lindzey, Knopff, Jalkotsky, & Boyce, 2010;Mattson, 2014). In part, the financial benefits associated with prioritizing ungulates over carnivores, and the incentives to reduce competition for shared prey we value, has contributed to the over-hunting of large carnivores around the world (Darimont, Fox, Bryan, & Reimchen, 2015;Ordiz, Bischof, & Swenson, 2013;Treves, 2009).
Food webs and predator-prey dynamics are complex, however, and the assumption that limiting predators results in increased prey does not always play out in natural systems. For example, the number and abundance of prey types in community assemblages (e.g., primary vs. secondary prey; Holt, 1977), predator morphology and energetics (Carbone, Mace, Roberts, & Macdonald, 1999), and intraspecific (IS) variation in prey selection exhibited by individual predators (Araújo, Bolnick, & Layman, 2011) influence predator-prey dynamics.
Uniquely, research on IS among large carnivore species has been dominated by the discussion about how to manage specialists that preferentially select livestock and rare prey of special concern (e.g., the "problem" animal; Linnell, Odden, Smith, Aanes, & Swenson, 1999), rather than the intrinsic and extrinsic mechanisms that explain IS. A common strategy to reduce the number of specialist predators in a population is to increase predator hunting (e.g., Bergstrom et al., 2014;Rominger, 2017). Hunting predators, however, often skews the sex-and age-structure of predator populations (Milner, Nilsen, & Andreassen, 2007), which may hold consequences for prey selection as exhibited by different sex-or age-classes of top predators. For example, younger ageclasses of African lions (Panthera leo) and pumas (Puma concolor) select for smaller suboptimal prey, as they learn the requisite skills to hunt more dangerous large prey (Elbroch et al., 2017;Hayward, Hofmeyr, O'Brien, & Kerley, 2007).
The puma is a highly adaptable solitary carnivore distributed widely across the Americas. Across its range, pumas eat a great diversity of prey, but in any locale they tend to hunt the common ungulate species. Pumas are heavily managed across their range, mostly through legal sport hunting that is in part driven by a desire for increased ungulate populations for human hunters (Mattson, 2014;Todd, 2002). For example, Colorado Parks and Wildlife (CPW) recently launched a $4.5 million effort to cull pumas and American black bears (Ursus americanus) to bolster local mule deer (O. hemionus) populations (CPW, 2016), under heavy criticisms from the public and the scientific community (Ishikawa, 2017a(Ishikawa, , 2017b. There is conflicting evidence that controlling pumas aids ungulate populations (Forrester & Wittmer, 2013;Hurley et al., 2011), and circumstances appear to be case-specific. Pumas exhibit IS, and their IS can negatively impact rare prey populations in multiprey systems (e.g., bighorn sheep, Festa-Bianchet, Coulson, Gaillard, Hogg, & Pelletier, 2006; Rominger, 2017). The occurrence of pumas that specialize on rare prey is generally stochastic, however, and therefore small prey populations can persist longer than expected under predation pressure (Wittmer, Hasenbank, Elbroch, & Marshall, 2014). With regards to deer specifically, one longterm study showed that puma control increased fawn survival and winter doe-fawn ratios, but that these changes did not reliably translate into increased mule deer population growth (Hurley et al., 2011). Hurley et al. (2011) concluded that predator control was not an effective means of increasing mule deer population growth in the intermountain west.
Here, we tested for IS in a small population of hunted pumas in the Greater Yellowstone Ecosystem, United States, and more specifically identified foraging strategies shared by multiple individuals in a system with diverse ungulate and other prey. Our objective in this paper was to explore potential links between IS and hunting pressure, given current management strategies and the potential implications of skewing the sex-or age-structure of predator populations on predator-prey dynamics. Based upon previous research, we hypothesized that males and females would forage differently (Elbroch, Lendrum, Newby, Quigley, & Craighead, 2013). We predicted that males would exhibit a more specialized diet focused upon elk and deer, both abundant primary puma prey, because of their wider movements and larger territories (Elliot-Smith et al., 2015). We predicted females would exhibit a more generalist diet, as the energetic demands of maternal care would require them to hunt primary and secondary prey. We also hypothesized that resource distributions, more generally, would explain differences in any foraging strategies we detected among pumas (Robertson et al., 2015). This is timely research to contribute to the ongoing discussion over whether predator control aids ungulate populations, or instead continues without scientific support under external pressure from special interests.

| Study area
Our puma study spanned approximately 2,000 km 2 of the GYE north of Jackson, Wyoming (Figure 1). Elevations ranged from 1,800 to >3,600 m. The area was characterized by short, cool summers during which prey were widely dispersed and long winters with frequent snowstorms in which prey were aggregated (Elbroch et al., 2013). The study area supported a diverse community of large mammals. Carnivores included gray wolves, black bears (Ursus americanus), grizzly bears (U. arctos), coyotes (Canis latrans), and red foxes (Vulpes vulpes). Ungulates included elk, mule deer (O. hemionus), moose, bison (Bison bison), pronghorn (Antilocapra americana), bighorn sheep, and a very small number of white-tailed deer (O. virginianus). Additional small puma prey included North American porcupines (Erethizon dorsatum), American beavers (Castor canadensis), Northern raccoons (Procyon lotor), grouse, and diverse waterfowl.
2.2 | Puma captures, GPS collar programming, and prey sampling Following capture protocols approved by an Institutional Animal Care and Use Committee, pumas were captured and fitted with global positioning system (GPS) collars (Lotek Iridium and Globalstar S collars, Newmarket, Ontario; Vectronics Globalstar Plus, Berlin, Germany). GPS collars acquired location data at 2-hr intervals. We identified GPS clusters visually in ArcGIS, which we defined as any ≥2 locations occurring within 150 m, and spanning 4 hr to 2 weeks in duration (Elbroch et al., 2017). We transferred puma location data to handheld GPS units to locate clusters in the field. Prey remains, including hair and bone fragments, were used to identify prey species, and signs of struggle around the carcass, the state of prey remains, presence and location of bite marks, and body parts consumed were used to determine whether the puma had killed the animal or was scavenging.

| IS analysis and identifying foraging strategies
For each puma, we binned prey selection into seven categories (elk, deer, bighorn sheep, moose, pronghorn, carnivores, and small prey). We maintained so many categories to capture potential specialists on rare ungulates of special concern (moose, pronghorn, and bighorn sheep). We lumped white-tailed deer and mule deer because we documented very few instances of puma predation on white-tailed deer (n = 5). Carnivores included all members of the Order Carnivora, regardless of their size. Small prey included all smaller vertebrate prey that was not carnivores, including birds, beavers and porcupines. We then standardized prey selection for individual pumas by dividing the number of prey they killed in each category by the total prey we documented for that puma, to account for unequal sampling across pumas and ensuring all prey classes were represented as a proportion of 1.
Based on a rarefaction analysis (Chao, 1987;Krebs, 2014), we only included pumas for which we had documented ≥10 ungulate prey in our analyses, and so we excluded two pumas before we began. Then we created a Bray-Curtis dissimilarity distance matrix to represent each puma's prey selection in space. We applied nonmetric multidimensional scaling (MDS) to assess the spatial relations between individuals and to determine whether there existed foraging strategies shared among groups of pumas (Dickman & Newsome, 2015).
In MDS, multiple parameters are fit to a single axis, much like principal components analysis (PCA), however, unlike PCA, MDS does not assume a linear relationship among the data being fitted, nor does it employ eigenvalues to assign priority to any axis, influencing its explanatory weight (Holland, 2008). Determining the correct number of dimensions in which to fit the data is somewhat arbitrary, but it is generally recommended that researchers use ≤3, otherwise interpreting the results becomes unwieldy. We ran the analysis with both two and three dimensions (3D), given that three significantly reduced the "stress" as determined with a Kruskal's Stress test (formula 1) (Kruskal & Wish, 1978), indicating a better goodness of fit for our data (Holland, 2008). The results resulting from two and 3D were essentially identical, with regards to identifying pumas with similar foraging strategies, and thus we presented a graphical depiction of our results in two dimensions (2D) for ease of interpretation. Ultimately, individuals were assigned to foraging strategies based upon distance correlations outputted from the 3D MDS analysis.

| Explanatory drivers of foraging clusters
We first applied a power analysis to determine the probability of rejecting a false null hypothesis in an analysis with five foraging categories and 13 individuals was 0.52. This is less than ideal power, however, we believed it was sufficient to provide ecological insights as long as our analyses identified a top model with high fit to our data. We tested both intrinsic and extrinsic covariates that might explain puma foraging strategies. Specifically, we tested whether a puma's location in our study system influenced prey selection due to potential differences in prey availability. We used the Slate-Crystal line, a north-south water course to define a binary FIGURE 1 Study area in the Greater Yellowstone Ecosystem, located northeast of the town of Jackson, Wyoming. The rectangle indicates the area where we caught and monitored pumas intensively covariate east-west, as it was a natural boundary to which puma home ranges appeared to naturally fall to one side or the other over the course of the project. We also used the mean elevation of puma location data as a potential indicator of elk abundance in our system, as elk tended to be found higher than deer throughout the year (Sawyer, Lindzey, & McWhirter, 2005;Smith, 2007). We also included puma sex and age as intrinsic covariates.
Given our small sample of pumas (n = 13), we built seven a priori models of 1 or 2 covariates maximum (Table 1), and employed generalized linear models and a Poisson distribution, to test what biological factors best fit our selection parameter, foraging strategy. Before model building, we tested for multicollinearity among covariates; all covariates were correlated at r < |0.5|, and all variables were included in our analyses. We compared model performance based upon Akaike's Information Criterion adjusted for small sample size (AICc), ΔAICc, and Akaike weights (w i ) (Burnham & Anderson, 2002), and assessed top model performance with R 2 values.

| Post hoc analysis of drivers of prey weight
Based upon the results of the MDS analysis, we conducted a post hoc analysis to further test whether pumas select larger prey with age. There are many factors that influence prey selection, including season. Pumas in general eat smaller prey in summer, regardless of their age (Knopff, Knopff, Kortello, & Boyce, 2010). Therefore we tested three simple a priori models (Table 3) with generalized linear mixed models and a Poisson distribution, to see which best fit our selection parameter, prey weight, as determined from the published literature for sex and ageclasses of each prey type (see Elbroch, Allen, Lowrey, & Wittmer, 2014 for more on determining prey weight). For this analysis, we only included kills (n = 445) from pumas age 0-70 months, during which our MDS analysis showed transitions in foraging strategy. We tested for the effects of puma age in months, two seasons (summer from 1 April-31 October and winter from 1 November-31 March), and included individual pumas as a random effect to account for unequal sampling. We compared model performance based upon AICc, ΔAICc, and w i (Burnham & Anderson, 2002). We assessed top model performance with R 2 values and reported parameter estimates for covariates in top models (β ± SE).

| Foraging strategies
We identified five distinct foraging strategies among pumas ( Figure 2, Table 2). The most common foraging strategy was an elk-specialist diet, in which animals killed elk more than any other prey. The second strategy, which we called "mixed," ate elk, deer, and small prey in similar proportions. The third strategy, exemplified by two independent subadult females, specialized on deer. And the last two foraging strategies each described one puma.
We identified one top model that best fit our data to explain ecological drivers of foraging strategies that captured 96% of the Akaike weight. Our top model performed very well (R 2 = 0.78) and included the single covariate, Age. In short, the older the animal, the larger the prey it specialized upon ( Figure 3).

| Posthoc analysis of drivers of prey weight
Our top model included all variables (R 2 = 0.32) and performed significantly better than the remaining two models in terms of AICc and w i (Table 3). Season (β = −33.6; F 1,12.75 = 12.07, p < 0.001) was a much stronger predictor of prey weight than puma age (β = 1.6; F 1,424 = 81.75, p = 0.004), but both proved significant.

| DISCUSSION
Predator hunting may indirectly and unintentionally influence predator-prey dynamics by changing the age-and sexstructure of predator populations that exhibit IS variation in prey selection. Our top model was composed of a single TABLE 2 Foraging strategies exhibited by individual pumas, and the number of kills included in the analysis (n) for each puma, the mean elevation of their territory, a binary covariate for whether their territory was east or west of the water body Slate/Crystal, the puma's age (months) and puma sex intrinsic characteristic, Age, and explained a large portion of the variation in foraging strategies detected in this small puma population. Our sample of pumas was small, even if the number of prey sampled was large, raising legitimate questions as to whether this pattern will hold true with larger samples, or across puma populations. Nevertheless, our top model performed very well and the ramifications of our results are provocative to contemplate as they pertain to carnivore management, predator culling, and sex-specific predator hunting (e.g., different limits on male vs. female harvest). For future research on this topic, we encourage a repeated measures approach in which researchers collect sufficient time series data on individual puma prey selection to better assess whether pumas indeed transition prey selection as they age as predicted by our conclusions. Our results suggest that the current controversial strategy of puma culling to aid mule deer in multiprey systems, and more broadly, the strategy of culling any predator to aid one prey type in multiprey systems, may result in unexpected consequences. Previous research suggested that heavy puma harvest may increase rather than decrease puma-human conflict because of changes in puma sex-and age-classes in the population (Teichman, Cristescu, & Darimont, 2016), and here we show that heavy puma harvest and puma culling may also exacerbate problems for mule deer if the agestructure of the puma population changes to predominantly younger animals more likely to hunt deer over elk. In fact, little science supports puma culling to aid deer populations (e.g., Hurley et al., 2011), although there is some evidence that it may aid rare ungulate prey, especially where pumas are subsidized by an abundant primary prey (e.g., translocations of desert bighorn sheep, Rominger, 2017;Cain & Mitchell, 2018).
Based on our research and that of Teichman et al. (2016), we argue that any predator control be conducted within an experimental design, to determine whether predator prey selection changes as a result of hunting pressure, and whether prey populations grow as a direct consequence of predator removal. We need rigorous research to scrutinize the influence of predator removal on ungulate population dynamics in natural systems to assess its implementation; we need defensible science to articulate to the public why at times, predator control is necessary (e.g., exotic predator control to aid declining native prey), to limit the abuse of predator control occurring under pressure from special interest groups, and to promote greater inclusivity of all people in wildlife conservation and management (e.g., nonconsumptive constituents).