Multiple anthropogenic interventions drive puma survival following wolf recovery in the Greater Yellowstone Ecosystem

Abstract Humans are primary drivers of declining abundances and extirpation of large carnivores worldwide. Management interventions to restore biodiversity patterns, however, include carnivore reintroductions, despite the many unresolved ecological consequences associated with such efforts. Using multistate capture–mark–recapture models, we explored age‐specific survival and cause‐specific mortality rates for 134 pumas (Puma concolor) monitored in the Greater Yellowstone Ecosystem during gray wolf (Canis lupus) recovery. We identified two top models explaining differences in puma survivorship, and our results suggested three management interventions (unsustainable puma hunting, reduction in a primary prey, and reintroduction of a dominant competitor) have unintentionally impacted puma survival. Specifically, puma survival across age classes was lower in the 6‐month hunting season than the 6‐month nonhunting season; human‐caused mortality rates for juveniles and adults, and predation rates on puma kittens, were higher in the hunting season. Predation on puma kittens, and starvation rates for all pumas, also increased as managers reduced elk (Cervus elaphus) abundance in the system, highlighting direct and indirect effects of competition between recovering wolves and pumas over prey. Our results emphasize the importance of understanding the synergistic effects of existing management strategies and the recovery of large, dominant carnivores to effectively conserve subordinate, hunted carnivores in human‐dominated landscapes.

2017; Treves & Bruskotter, 2014). Regardless of human motivations, large carnivores are particularly susceptible to human-caused mortalities due to slow life histories that include a reliance on adult female longevity and long interbirth intervals (Darimont, Fox, Bryan & Reimchen, 2015). Without management interventions that aid carnivores, many apex predators are predicted to decline to extinction (Ripple et al., 2014;Wolf & Ripple, 2016).
Management strategies to conserve and restore large carnivores include reintroduction and translocation efforts (Seddon, Griffiths, Soorae & Armstrong, 2014). These strategies are mostly implemented in Europe and North America, where increasing tolerance may allow viable populations of large predators to exist in humandominated landscapes (Chapron et al., 2014). Studies following the recovery of large carnivores have also made important contributions to our current understanding of the direct and indirect effects of predators upon ecological communities, including predator-prey interactions and trophic cascades affecting community assemblages (e.g., Ripple & Beschta, 2011). The effect of recovering carnivores on other, subordinate, carnivores has received much less attention. This is surprising given known effects of competition within carnivore guilds, including reduced abundance and survivorship among subordinate carnivores (Harihar, Pandav & Goyal, 2011;Levi & Wilmers, 2012;Roemer, Donlan & Courchamp, 2002) and shifts in space use or prey selection by subordinate carnivores (Harihar et al., 2011;Lendrum et al., 2014;Ruth et al., 2011).
The reintroduction of wolves (Canis lupus) to Yellowstone National Park, USA, began in 1995 and is touted as one of the most successful conservation stories of all time (Smith & Ferguson, 2012).
They are also a conservation success story, in that puma populations rebounded in the west of the United States and Canada after 1965, when wildlife managers in nearly every western state stopped paying state bounties for killing pumas and introduced managed puma hunting with limits in restricted seasons (Mattson & Clark, 2010).
Previous research has suggested that pumas are subordinate competitors in the presence of wolves (Kortello, Hurd & Murray, 2007;Ruth, 2004), but research has failed to demonstrate clear fitness consequences for pumas competing with wolves despite numerous recorded observations of wolves killing pumas (e.g., Kunkel, Ruth, Pletscher & Hornocker, 1999;Ruth et al., 2011). Pumas are also a game species legally hunted throughout much of their range. High hunting pressure can severely reduce the abundance of pumas and destabilize established population structures (Robinson, Wielgus, Cooley & Cooley, 2008;Stoner et al., 2013). Thus, while direct fitness consequences of recovering wolf populations on pumas are likely, demonstrating them requires data and analytical approaches capable of separating their potential effects from those attributable to other limiting factors, such as hunting.
Our objective was to test for the effects of recovering wolves on puma survival in the southern Greater Yellowstone Ecosystem (GYE), where existing management objectives also supported the legal hunting of pumas and decreasing the local elk population through "liberal" hunting seasons (WGFD, 2014). To achieve our objective, we used 14 years of monitoring data from 134 individually marked pumas collected concurrent with the recolonization of wolves and elk reduction to determine drivers of puma survival and cause-specific mortality rates.

| Study area and wolf reintroductions
Our study area encompassed approximately 2,300 km 2 of the GYE in southern Teton County, Wyoming (Appendix S1, Figure A1.1). The area was constrained by Yellowstone National Park to the north, Grand Teton National Park to the west, and the National Elk Refuge to the south and supported an abundance of ungulate prey, both in terms of species (n = 7) and biomass. Detailed descriptions of climate, topography, habitat, and community composition including ungulate prey species are presented in Elbroch, Lendrum, Newby, .
Wolves were first reintroduced north of our study area in Yellowstone National Park in 1995 (Smith & Ferguson, 2012 Elk in our study area were part of the migratory Jackson herd and cooperatively managed by the Wyoming Game and Fish Department (WGFD), National Park Service, and the USFW's National Elk Refuge (NER). The Jackson elk herd typically travels long distances and congregates adjacent the town of Jackson, WY, where they receive supplemental feeding in winter on the NER and adjacent public lands on feed lots managed by the WGFD. Our study covered the time period in which managers implemented liberal hunting quotas across jurisdictional boundaries to reduce the Jackson herd from 16,000 in 2000 to 11,000 animals, of which managers wanted 5,000 to win-

| Puma captures, monitoring, and age classifications
Puma monitoring began in 2001 and we followed puma capture and immobilization protocols described in Elbroch et al. (2013) and approved by the Jackson Institutional Animal Care and Use Committee (Protocol 027-10EGDBS-060210). We fit pumas with a VHF (Telonics, Mesa, AZ) or GPS (Telonics, Mesa, AZ; Televilt, Lindesberg, Sweden; Vectronics, Berlin, Germany; Lotek Wireless, Ontario, Canada) collar. We hand-captured kittens between 5 and 7 weeks old without the aid of immobilization drugs and fit them with custom-made, lightweight, expandable VHF collars (Telonics, Mesa, AZ, USA).
Attempts to locate kittens wearing VHF collars were made every 2 days. All other pumas wearing VHF collars were located at minimum weekly from the ground and monthly from aircraft. Location data were acquired by GPS collars 4-12 times per day. All collars were equipped with mortality sensors, which alerted researchers when an individual had not moved for ≥8 hr. We investigated mortality sites and determined the cause of death through interpreting field signs (e.g., bite marks, footprints), necropsies conducted with a veterinarian, and based on blood and tissue samples analyzed by the Wyoming State Veterinary Laboratory.

| Determining minimum annual puma densities
Each year, we determined minimum puma density in our study area based on overlapping home ranges (Rinehart, Elbroch & Wittmer, 2014). Annual home ranges for adult pumas were determined using fixed-kernel density estimators (Worton, 1989) in ArcGIS 10, and isopleth calculations in the Geospatial Modeling Environment (Beyer, 2012); methods are further described in Lendrum et al. (2014).
We determined the boundaries of the area in which we consistently searched for pumas each winter, and in which we believed we had captured all resident pumas. In ArcGIS 10, we created a polygon of our capture area and quantified each puma's residency within this polygon (Rinehart et al., 2014). "Minimum puma densities" for our 890 km 2 capture area were then determined by summing the residency estimates for all adult pumas with overlapping home ranges for each year. We also scaled density estimates to pumas per 100 km 2 for comparisons with previous research, aware that scaling introduces extrapolation bias impacting the precision of estimates (Schonewald-Cox, Azari & Blume, 1991;Smallwood & Schonewald, 1998).

| Puma survival analyses
We estimated puma survival probabilities and cause-specific mortality rates (m i ; Schaub & Pradel, 2004;Marescot, Forrester, Casady & Wittmer, 2015) with multistate capture-mark-recapture (CMR) models in E-SURGE (Choquet, Rouan & Pradel, 2009;Lebreton, Nichols, Barker, Pradel & Spendelow, 2009). CMR models were best suited for analyzing our data because of their ability to account for both right and left censored data and to accommodate encounter histories based on different data sources (e.g., kittens observed in dens, individuals monitored with VHF vs. GPS technology) and sampling intervals (Cubaynes et al., 2014;Devineau et al., 2010). However, standard CMR models have to meet multiple assumptions to avoid biasing parameter estimates and model overparameterization (Choquet, Lebreton, Gimenez, Reboulet & Pradel, 2009;Fletcher et al., 2012). Therefore, prior to final model selection, we employed a range of goodness-of-fit (GOF) tests to determine if our data met the assumptions for CMR models, to test whether our methods may have biased estimates of survival or recapture probabilities, and to determine whether our models fit the data and explained variation in our selection parameter, puma survival (Appendix S3). We also tested for potential overdispersion in our data due to siblings in the same litter dying from the same cause more than expected under the assumptions of independence (Ruth et al., 2011), and subsequently adjusted AIC values using the variance inflation factor (ĉ; Appendix S3).
For our analyses, we categorized pumas into three age classes based on differences in life histories and survival reported in the literature.
We defined kittens (n = 75) as individuals <6 months old. Kittens are completely dependent on their mothers and experience high mortality from both predation and starvation (Logan & Sweanor, 2001;Ruth et al., 2011). We defined juveniles (n = 22) as individuals ≥6 and <18 months old. Juveniles remain dependent on their mothers, but are less susceptible to predation (they better avoid predators by climbing or running); (Logan & Sweanor, 2001)). Juveniles experience higher risks of starvation, can be legally hunted once they are one year old and separate from their mothers (WGFD, 2006), and experience risks associated with dispersal (Quigley & Hornocker, 2010;Stoner et al., 2013). We pooled all individuals ≥18 months old into an adult age class (n = 37) when pumas are expected to establish stable territories and become reproductively active (Logan & Sweanor, 2001;Quigley & Hornocker, 2010).
At each time step of the model, individuals (i.e., kittens, juveniles, adults) occupied one of the following seven states: "alive" (survival rate in matrix below); "recently dead from hunting, poaching, or management action" (i.e., human causes; cause-specific mortality rate m h ); "recently dead from predation by wolves, bears (Ursus spp.), or other pumas" (cause-specific mortality rate m p ; the small number of predation events limited our ability to further differentiate between predator species); "recently dead from starvation" (cause-specific mortality rate m s ); "recently dead from other natural causes including disease and exposure during cold weather" (cause-specific mortality rate m o ); or "recently dead from unknown causes" (cause-specific mortality rate m u ). All dead individuals were eventually assigned to a permanent "dead" state independent of their actual cause of mortality (Lebreton, Almeras & Pradel, 1999). We also censored emigrating pumas (i.e., "alive outside the study area") from these analyses, meaning that they were only included until their departure from the study site, ignoring their subsequent fates. Censoring dispersers in this way mitigated an inflation in mortality estimates due to a reduced sample size, as well as allowed us to quantify mortality rates specific to our study population.
The survival transition matrix S from time t to time t + 1 was written as The matrix could further be decomposed into two equivalent biological processes, survival and its associated probability , and the probability of mortality by a specific cause m i , as is equal to the complementary of the sum of mortality rates (Schaub & Pradel, 2004).
We began model building by identifying biologically relevant covariates that might explain change in puma survival and causespecific mortality rates, from which we built competing a priori models to test in an information theoretic framework (Burnham & Anderson, 2002). We also tested whether recapture/detection probabilities were time-dependent (Culina, Lachish, Pradel, Choquet & Sheldon, 2013), before building candidate models from puma covariates, including age (i.e., kittens, juveniles, adults), and additional ecological and time-based covariates described below.
We included the following numeric covariates in our models: (a) annual wolf counts for our study area as reported by the USFW, annual puma harvest numbers for Unit 2, the hunting unit in which we studied pumas, and (e) annual minimum puma densities, calculated as described above (resident adults/890 km 2 ; Appendix S4, Table A4.1). We used two elk metrics as a measure of prey availability and bottom-up effects, as elk were the primary prey for pumas in our study area (Elbroch et al., 2013): first, we employed total elk numbers in the Jackson herd, and second, the portion of the Jackson elk herd off-Refuge, to highlight those elk more likely truly available to pumas in our study (Elbroch, Lendrum, Robinson & Quigley, 2016).
We excluded density/abundance estimates of wolves, pumas, and elk for 2001 and 2015 from the analysis because puma capture efforts did not occur throughout these years.
We also included two time-based covariates that captured greater ecological complexity than numeric covariates without overparameterizing models; time-based covariates, however, also introduced complexity when interpreting results. First, we tested for potential seasonal variation in survival probabilities (ϕ) and cause-specific mortality rates to account for environmental variability due to weather, prey availability, puma foraging behaviors (Elbroch et al., 2013), and anthropogenic top-down effects, primarily human hunting (Wyoming Game and Fish Department, 2006).
We split each year into two 6-month seasons, following legal hunting periods for pumas, which we expected to yield differences in age-specific survival and mortality rates because human-caused mortality is generally the driver of puma population dynamics across hunted populations (Quigley & Hornocker, 2010;Stoner et al., 2013): (a) we defined the "hunting season" as 1 October to 31 March of the following year, during which pumas were legally hunted. The hunting season also captured the following additional ecological variation: elk returned to low-elevation winter ranges in November and aggregated in large herds near supplementary feeding stations, deer migrated out of the study area, competition between wolves and pumas likely increased near shared prey, and deep snows and cold temperatures increased the risk of starvation (Elbroch, Lendrum, Newby, Quigley & Thompson, 2015;Elbroch et al., 2013); (b) We defined the "nonhunting season" as 1 April-30 September, during which puma hunting was closed, and during which elk migrated to summer ranges at higher elevations and TA B L E 1 Model selection results for our a priori CMR models estimating survival (ϕ) and cause-specific mortality rates (m) for humancaused, predation, starvation, other, and unknown mortalities of pumas in the Southern Yellowstone Ecosystem. Recapture probabilities (p) were constant (i) given that models accounting for temporal variation in detection were nonidentifiable. Covariates included annual wolf counts (Nwolf), annual counts of the Jackson elk herd (Nelk), annual counts of the Jackson herd that wintered off the National Elk Refuge (NelkOff), annual puma abundances (Npuma), annual pumas harvested in Hunting Unit 2 (Nharvest), and two time-dependent covariates: 2001-2005 versus 2006-2015, reflecting changes in wolf density (low/high), and a seasonal comparison based upon the 6-month hunting season for pumas (hunt/no-hunt). All models estimated age-specific differences (kitten, juvenile, adult) in survival probabilities and causespecific mortality rates of resident pumas across years, and here we report model descriptions: the number of parameters (n), deviance, QAIC c , ∆QAIC c , and QAIC w initiating a steep period of wolf population growth (Appendix S2, Table A2.1). Elk numbers were highest in the "low-wolf" period, and more variable, but generally declining during the "high-wolf" period.
The low-wolf period was also characterized by high human harvest of pumas, and the high-wolf period by substantially lower but variable human harvest of pumas.
To mitigate issues to do with overparametizing models and including uninformative parameters (Arnold, 2010;Burnham & Anderson, 2002), and to more directly test our specific hypotheses, we devised a short list of competing a priori models (Table 1).
All models included one parameter for detection probability, five parameters representing the interaction effects of different cause-specific mortalities, and five parameters for initial state probabilities (Marescot et al., 2015;Schaub & Pradel, 2004); each model also included the parameters associated with additional numeric (e.g., annual elk counts) or categorical (e.g., seasons, age) covariates. Then we tested which model best fit our selection parameter, puma survival, through a comparison of QAIC c corrected for slightly overdispersed data. We considered the top model and any subsequent model differing by <4 QAIC c units to have produced substantial empirical support for explaining variation in puma survival and cause-specific mortality rates (Burnham & Anderson, 2002).

| Puma survival analyses
We identified the cause of mortality for 80 of 108 pumas that died during our study; seven died outside the study area and their fates were censored from the cause-specific mortality analysis. Predominant causes of mortality varied among age classes (Table 2). Predation was the predominant cause of mortality for kittens (9 were killed by wolves), starvation was the predominant cause of mortality for juveniles, and adults most frequently died from human causes.
Our CMR analysis inclusive of cause-specific mortalities resulted in two top models (Table 1). Our top-ranked model highlighted differences in age-specific (kittens, juveniles, adults) variation in survival and cause-specific mortality rates between hunting and nonhunting seasons (Table 3). This model received 2.8 times more empirical support than our second model based on an evidence ratio (Burnham, Anderson & Huyvaert, 2011). We also assessed age-specific survivorship among pumas, defined as the complementary of the sum of mortality rates estimated in our analyses. This method assumes mortality rates are additive, which we could not confirm for certain, however, previous research has shown that pumas do not reliably exhibit compensatory mortality under hunting pressure, as would be expected for a territorial species (Cooley et al., 2009) (Table 3).
Our second top model emphasized additional insights missed in the first model, and age-specific variation in survival and causespecific mortality rates correlated with changes in the number of Jackson elk wintering off the NER. All age classes of puma decreased in survival with decreasing elk availability off the NER, but juveniles were impacted most ( Figure 1). As elk wintering off the NER decreased, starvation across age classes, other natural causes of death such as exposure and disease across age classes, and predation on kittens increased (Figure 2).
TA B L E 2 Distribution of cause-specific mortalities among age classes, and across the hunting and nonhunting seasons. "Other causes" included natural mortality, primarily disease and exposure. Causes of mortality of pumas that died outside the study area are shown in parentheses "()" but were censored from analyses

| D ISCUSS I ON
The reintroduction and recovery of large carnivores inside and outside of protected areas are as much a social as an ecological triumph (Chapron et al., 2014 TA B L E 3 Cause-specific mortality rates attributable to human, starvation, predation, other, and unknown causes. Estimates and standard error in brackets were adjusted using the Heisey and Fuller (1985) method for three age classes of pumas resident in the study area, during the hunting and nonhunting seasons F I G U R E 1 The relationship between elk numbers off the National Elk Refuge (NER) and age-specific puma survival of mortality for pumas throughout the western United States and Canada (Quigley & Hornocker, 2010), and thus these results were not unexpected. Higher kitten and juvenile starvation in the hunting season may also have been influenced by the harvest of females that orphaned kittens too young to forage and defend themselves (Wyoming Game and Fish Department, 2006). Hunting limits for pumas were reduced over the course of our study, however, further reductions in Unit 2 or additional reductions in adjacent hunting units to increase immigration rates (sensu Robinson et al., 2008) may be necessary to facilitate stability in the local puma population.
The second anthropogenic influence on puma survival was the reduced availability of elk, their primary prey (Elbroch et al., 2013).
Top predators are generally regulated by prey availability (Wallach, Izhaki, Toms, Ripple & Shanas, 2015), and in our system it was juvenile survival that was impacted most, followed by adult puma survival ( Figure 1). As elk decreased, starvation across age classes and predation on kittens both increased (Figure 2 WGFD, 2016). Through these varied competition mechanisms, wolves likely contributed to increased puma starvation across age classes, as has been observed following wolf recolonization in other areas (Kortello et al., 2007;Ruth, 2004). Over the same time period, elk availability off the NER dropped by 70% and wolf numbers increased >600% (Figure 3). These results suggest that humans-through their influence on top-down and bottom-up forces-have successfully facilitated a change from a system dominated by pumas since 1926, when wolves became extirpated (Haines, 1996) to one dominated by wolves, as in historic times. Given the differences between pumas and wolves in social organizations (solitary felid vs. social canid), hunting behaviors (ambush vs. cursorial), and species-specific carrying capacities (lower vs. higher), changing relative predator densities is likely to result in wide-reaching ecosystem changes. Predicting and addressing the additive effects of multiple, management actions are difficult. Once they occur, however, they require flexible conservation strategies that encourage the coexistence of people and predators (Chapron et al., 2014) and better link the comanagement of predators and prey.

ACK N OWLED G M ENTS
We thank the Summerlee Foundation, Charles Engelhard

CO N FLI C T O F I NTE R E S T
None declared.
F I G U R E 3 Annual wolf counts, elk numbers off the National Elk Refuge (NER), and estimated number of pumas for our study across our approximately 2,300 km 2 study area (see Appendix S4)

AUTH O R CO NTR I B UTI O N S
LME, HQ, and DC conceived the project and conducted the research. LM conducted the CMR analyses, with input from HW and LME. LME, LM, and HW drafted the manuscript, and HQ and DC contributed feedback. All authors approved the final version.

DATA ACCE SS I B I LIT Y
Data will be archived in Dryad.