Recent Trends in Survival and Mortality of Wolves in Minnesota, United States

Survival is a key determinant of population growth and persistence; computation and understanding of this metric is key to successful population management, especially for recovering populations of large carnivores such as wolves. Using a Bayesian frailty analytical approach, we evaluated information from 150 radio-tagged wolves over a 16-year time period to determine temporal trends and age- and sex-specific survival rates of wolves in Minnesota, United States. Based on our analyses, overall annual survival of wolves during the study was 0.67, with no clear evidence for age- or sex-specific differences in the population. Our model demonstrated statistical support for a temporal trend in annual survival; the highest survival was predicted at the beginning of the time series (0.87), with lowest survival (0.55) during 2018. We did not observe evidence that survival was markedly reduced during years when a regulated hunting and trapping season was implemented for wolves (years 2012–2014). However, cause-specific mortality analysis indicated that most mortality was human-caused. While the estimate for increasing human-caused mortality over time was positive, the evidence was not statistically significant. Anthropogenic causes resulted in ∼66% of known mortalities, including legal and illegal killing, and vehicular collisions. Trends in wolf survival in Minnesota may reflect an expanding distribution; wolf range has spread to areas with more human development during the study, presumably leading to increased hazard and reduced survival. Our results provide foundational information for evaluating and guiding future policy decisions pertaining to the Great Lakes wolf population.


INTRODUCTION
Human expansion and persecution have endangered terrestrial carnivores worldwide (Ripple et al., 2014). Some carnivore species or populations however, have recovered to varying degrees from past exploitation, expanding from "refuges" back into parts of their former range. As their populations expand or human populations encroach new areas, large carnivores compete with humans for space or wild prey and may kill or injure pets, livestock and humans, while humans often respond by removing individuals, groups or entire populations (Cardillo et al., 2004;Treves and Bruskotter, 2014). Carnivore conservation in the Anthropocene is thus a challenge, and characterization of key vital rates such as survival is crucial for aiding recovery and understanding population dynamics.
Although the ecological roles of carnivores in their respective ecosystems lack generalities and fuel discourse (reviewed in Ford and Goheen, 2015), large carnivores indisputably remain as flagbearers for conservation of large landscapes (Ritchie et al., 2012). For their conservation roles and their connection to the human psyche, restoration efforts continue to be implemented and evaluated, and have often resulted in optimistic recoveries of carnivore species such as the tiger Panthera tigris (Jhala et al., 2021), puma Puma concolor (Jansen and Logan, 2002), Asian lion P. leo leo (Jhala et al., 2019), brown bear Ursus arctos (Lamb et al., 2018), gray wolf Canis lupus (Wydeven et al., 2009;Mech and Boitani, 2010;Ripple and Beschta, 2012;Nowak and Mysłajek, 2016), African wild dog Lyacon Pictus (Gusset et al., 2010), and black-footed ferret Mustela nigripes (Jachowski et al., 2011).
However, the Anthropocene has disparately affected the persistence and distribution of carnivore populations (Ripple et al., 2014), rendering certain sub-populations more vulnerable than other populations of the same species . Thus, it is imperative to prioritize management actions for effective conservation. An understanding of population dynamics is necessary for successful management (Caughley, 1994), whether for: (i) increased protection, (ii) supplementation, (iii) reintroduction, and/or (iv) population control. For example, wolf populations in Yellowstone and Isle Royale have been restored through reintroductions, recolonization, and supplementation (Ripple and Beschta, 2012;Hervey et al., 2021) and the Gir lions in India have increased from less than 50 animals to ∼700 in the past 100 years because of committed protection and conservation (Jhala et al., 2019). Some carnivore populations have also been controlled to maintain ecological carrying capacity and to reduce human-carnivore conflict (e.g., American black bears Ursus americanus, Garshelis et al., 2020).
Wolves are iconic predators across their extant range, but perhaps nowhere are their contemporary fates so entwined with human actions as they are in North America. Wolves evoke strong and polarizing reactions of support and persecution, and are thus involved in intense conservation conflicts. Such conflicts are often aggravated because people connect with wolves as symbols of pristine wilderness, reconciliation, invasion, disease, and government overreach. As a consequence, the status of wolves still ranges from complete to no legal protection, resulting in a mosaic of management emphasis across regional to national scales (Bruskotter, 2013). This was recently exemplified in the dialog around the delisting of wolves from the federal Endangered Species Act (ESA) 1 in the United States. Historically, the last remaining wolves in the "Lower 48 states" were protected in Minnesota and gradually expanded to repopulate northern Wisconsin and Michigan's upper peninsula 2 . Consequently, Minnesota has a long history of successful wolf management aided by scientific monitoring, including intense 1 https://www.fws.gov/home/wolfrecovery/ regional study (e.g., Olson, 1938;Van Ballenberghe et al., 1975;Fuller, 1991;Mech et al., 2000;Erb and Benson, 2004;. However, wolves in Minnesota still face threats from habitat alteration, and mortalities from escalating linear infrastructure, roadkill, depredation management, and illegal killing 2 . Wolves in North America, with the exception of the Northern Rocky Mountain population, have been re-listed on the Endangered Species Act; the Minnesota population considered as 'threatened' 1 . It is therefore imperative to monitor the population for long-term viability. Annual survival rate is a key parameter that informs our understanding of the ecological dynamics and persistence potential of a population. Herein, we determine survival rates for wolves with data from 150 radio-or GPS-collared wolves spanning 16 years (2004-2019) across Minnesota. Although estimates of annual or seasonal survival rate can provide important information for management, awareness of any temporal trend in survival can be crucial for policy formulation. Such a temporal analysis of survival becomes even more necessary when populations are affected by environmental stochasticity, human impacts that vary over time, and socio-biological factors (such as territoriality, competition and density dependence), which often affect groupliving carnivores (Cubaynes et al., 2014;O'Neil et al., 2017). Using long-term known-fate data on individuals, we estimate human versus natural mortality while testing for the effect of time on wolf survival with four major questions: (1) how does annual survival vary between years, (2) does survival show a trend over time, (3) is there a particular inflection point or period where survival rates or trends change, and (4) did recreational wolf hunting and trapping (years 2012-2014) affect survival rates of wolves in Minnesota.

Study Area
We evaluated wolf survival and mortality within northern Minnesota, United States. Wolf distribution and abundance has expanded south and west within the state since the 1970s 2 (Figure 1). The study area was primarily characterized by mixed northern hardwood forest, bog, wetland, and agricultural land cover types (Erb and Sampson, 2013;Erb and Humpal, 2021). Human population and road densities were generally low (typically < 5 humans / km 2 and < 2 km of roads / km 2 ; MN DNR unpublished data), with primary land uses being recreation, logging, and some mining, with increasing agriculture in the south and west of wolf range. Prey and food sources for wolves included white-tailed deer (Odocoileus virginianus), beaver (Castor canadensis), moose (Alces alces), and sometimes fish and berries (Gable et al., 2016Homkes et al., 2020).
After ESA protections were established in the 1970s, public harvest of wolves was prohibited until their delisting in 2012. However, lethal control of wolves in response to livestock depredation was permitted throughout the study period and region (available reporting indicates that the annual number of wolves killed ranged from 95 to 216 for 2011-2020). During 2012-2014, state law permitted a regulated, public harvest of wolves. These hunting and trapping seasons occurred in the late fall and early winter months, with a total of 915 wolves killed (Erb and Sampson, 2013). Federal (ESA) protection of wolves has varied during the period of our analysis, but public harvest has not been allowed since 2014 2 .

Data Collection
For radiotelemetry, most wolves were captured using foothold traps (Erb and Humpal, 2021). Some wolves were also captured with the use of live-restraining neck snares (Gese et al., 2019), and a few by aerial darting from a helicopter. Wolves captured with foothold traps were captured between May-October and livesnared wolves were caught between December-March, however, winter live-snaring only started from 2013. Upon capture, each individual was weighed, sexed, and aged prior to release. Based on tooth-wear and color, coupled with appearance and behavior, individuals were aged into three categories: pups (<1 year), juveniles (1-2 years) and adults (>2 years). Post-mortem aging of some of the tagged wolves from tooth cross sections aligned with our initial assessments. Telemetry equipment ranged from VHF-only (20%) to VHF/GPS collars (remaining 80%). Most GPS radio-collars were programmed to take 3-6 locations per day, and wolves fitted with VHF-only radio-collars were relocated at approximately 7 to 10-day intervals throughout the year, and in some cases, primarily from early winter through spring (Erb and Humpal, 2021). All captures were done as per regulations and guidelines from the Minnesota Department of Natural Resources.

Statistical Analysis
We used a Bayesian shared frailty model (Halstead et al., 2012;Heisey, 2012) to capture variation in annual survival of wolves from 2004 to 2019, and as a function of sex and estimated age at capture. In addition, we partitioned the hazard rate from the frailty model into cause-specific mortality rates over the same time period (Heisey and Patterson, 2006). For each individual wolf, we created encounter histories that reflected the time period (number of days) between initial date of capture and date last known alive. Individuals that were characterized by "loss to follow-up" had undetermined fates and were right censored with an encounter history endpoint being the last alive date (DeCesare et al., 2016;Moore, 2016). Individual collars that were detected in mortality mode were located in the field, where fate was determined. The time period between last known alive and date confirmed dead was interval censored, such that date of death was considered unknown and was imputed by the frailty model (i.e., fate could have occurred any time between last date known alive and date confirmed death). Known fate mortality was classified as human (e.g., legally or illegally killed, vehicle collision), natural (e.g., disease, starvation, intraspecific strife), or unknown (cause undetermined). We generated time-varying covariates to estimate temporal effects and age. Specifically, we set the initial day of season, year, and age for each individual encounter history. As encounter histories progressed, the corresponding covariates for year and age were updated in alignment with each individual encounter history across time. Age was classified as either adult or non-adult (pup or juvenile) due to uncertainty in age estimation. As such, pups and juveniles graduated to adults after 2 and 1 years, respectively, with the date of graduation set to the 15th of April to correspond with the approximate wolf biological year.
We identified five a priori candidate models to test relative support for each of our research questions, specifically whether survival changed over time, whether a change in survival (e.g., an inflection point) was evident during the study, and whether years with legal hunting and trapping resulted in lower survival. The shared frailty model infers survival across a specified encounter history interval via a hazard function, expressed as where UH represents the unit hazard, with the unit defined as length of interval (daily, weekly, monthly). In this model, γ 0 is the intercept, providing a constant baseline hazard that can be offset by any number of fixed covariates X, or random effects κ. We specified the following models to evaluate evidence for temporal and/or hunting effects on survival: For models including a harvest effect, the years associated with public harvest were represented by a binary indicator covariate for harvest (x hunt = 1) vs. no harvest (x hunt = 0). We included sex and age (time-varying) in all models, which were also treated as binary fixed covariates, where the respective coefficients, β 1 and β 2 , modeled the influence of being male (x male = 1) vs. female (x male = 0), and adult (x adult = 1) vs. pup or juvenile (x adult = 0), respectively. We did not include a separate parameter for pups due to sample size, as very few pups with known fates were included in the dataset. We ranked models using the Bayesian widely applicable information criterion (WAIC; Vehtari et al., 2017). We also included results of leave-one-out (LOO) cross-validation and deviance information criterion (DIC) to check for consistency across evaluation metrics (Spiegelhalter et al., 2014;Vehtari et al., 2017). Following identification of the top performing model, we refit the model with an additional parameter to capture residual variation across years (e.g., year effects are not assumed independent of one another) and added a component to evaluate competing mortality sources. For each known mortality, cause was assigned as k ∈ {1 = human, 2 = natural, 3 = unknown}, using a categorical distribution with cause probabilities (p 1 , p 2 , p 3 ), where K k=1 p k = 1 (Cross et al., 2015;Stenglein et al., 2018). We tested evidence for the following mortality cause effects by relating covariates (α) to cause probabilities via the logit link function (Stenglein et al., 2018): first, human cause changes over time (relative to other causes), second, a difference in human vs. natural cause for adults relative to sub-adults, and third, a difference in unknown cause vs. other causes for adults relative to sub-adults.
Following estimation of the unit hazard, we calculated annual survival as where CH represents the annual, cumulative hazard and S represents the annual, cumulative survival rate. For reporting, we derived estimates of annual hazard and survival from the final model under the condition that ratios of male to female and adult to non-adult were equivalent to those observed during the study. In addition, we derived estimates of annual survival for each of the four age and sex classes (adultmale, adult-female, juvenile-male, juvenile-female) assuming an average year during the study.
We specified uninformative Normal prior distributions (µ = 0, τ = 1/σ 2 = 0.01) for all fixed covariate effects, as  Parameters were estimated with respect to the daily unit hazard (UH), and estimates were obtained from Gibbs MCMC sampling.
Frontiers in Ecology and Evolution | www.frontiersin.org well as the residual year effects. For the baseline hazard, we specified a weakly informative prior, γ 0 ∼ N(µ, σ 2 ), based on a comprehensive literature review of wolf survival (annual estimates ranging from 0.24-0.91), and used moment matching to appropriately place the hyperparameters on the UH scale (µ = −7.47, σ = 1.12). We estimated all model parameters for the frailty model using the Gibbs MCMC sampler in JAGS 4.2.0 (Plummer, 2003), by way of the jagsUI package (Kellner, 2019) implemented in R 4.0.3 (R Core Team, 2020). We ran three chains of 50,000 iterations and retained every 10th sample following a burn-in phase of 20,000. We calculated DIC using the jagsUI package and WAIC and LOOIC from the model's posterior predictive distribution using the loo package (Vehtari et al., 2020). We visually examined chains and calculated Gelman-Rubin statistics to verify chain convergence (r < 1.05). We report median values of posterior distributions and 95% credible intervals (CRI) for each parameter, unless otherwise stated.

RESULTS
We evaluated annual hazard and subsequent survival for 150 individual wolves during 2004-2019. Of these, 84 were females (56.0%) and 66 were males (44.0%). At the time of capture, 91 wolves were estimated to be adults (≥2 years; 60.6%), 44 were yearlings (∼1 year old; 29.3%), and 15 were pups (<1 year old; 10.0%). Mortality information (known fate) was obtained for 59 individuals. Of these mortalities, 39 were attributed to anthropogenic causes (66.1%), 14 were natural deaths (23.7%), and 6 had undetermined causes (10.2%). Our top frailty model was based on Eq. 4, with temporal change in the hazard represented by a log-transformation on the continuous year variable. This model indicated a curvilinear relationship between the hazard and time (Table 1 and Figure 2). The linear time trend model was also competitive (e.g., within 2 WAIC units of top model; Table 1), with less support for models representing an effect of regulated public harvest on the hazard rate. Our model demonstrated statistical support for a temporal trend in annual survival [β 3 = 0.54, 95% CRI = (0.07, 1.04), p > 0 = 0.99], with greatest survival generally occurring in the early years of the time series and lower survival during the later years (Table 2 and Figure 3), though the rate of change slowed over time (Figure 3). The highest survival was predicted at the beginning of the time series [0.87, 95% CRI = (0.68,0.97)], with lowest survival during 2018 [0.55, 95% CRI = (0.39,0.68)]. From all Bayesian information criteria rankings, we did not observe evidence that survival was markedly reduced during years when a public harvest was allowed on wolves (years 2012-2014; Table 1).

DISCUSSION
Survival is a key determinant of population growth and persistence, and computation and understanding of this metric is key to successful population management (e.g.,  Mech et al., 1998O'Neil et al., 2017 and understanding our relationship with wolves . Foremost, our results indicate a gradual increase in hazard and associated decline in median annual wolf survival in Minnesota over 16 years (Figure 2). While our results likely hint toward lower overall survival in Minnesota's wolves during the study period than the contiguous populations in Wisconsin and Michigan (Table 3), populations in those two states over our study timeframe were growing rapidly as wolves recolonized and recovered (Stenglein et al., 2015a;O'Neil et al., 2020). In contrast, the Minnesota wolf population was not in a consistent growth phase during most of the period of our analysis. Carnivore species often show higher survival in years when the population expands into new range or is recovering, while the rate diminishes as it nears carrying capacity because of density dependent factors such as competition and territoriality (Banerjee and Jhala, 2012;Cubaynes et al., 2014;O'Neil et al., 2017). The temporal trajectory of our study possibly mirrors a similar trend where higher survival rates were documented in the earlier periods of the time series (first 8-9 years) than the later (last 6-7 years) (Figure 2). The higher mortality in the later part of the time series might be because of intra-specific competition for food and space, density dependent population regulation, as well as human-induced causes. Wolf survey counts indicate that the population in Minnesota increased until 2007Minnesota increased until -2008 to about 2,500-3,000 individuals, reaching relative stability thereafter (Erb and Humpal, 2021), potentially corroborating density dependence as a mechanism that likely dampened annual survival and population growth. In addition, this population has expanded its distribution into areas with greater human population density and development (Figure 1). Although the most recent years when legal public harvest was allowed (2012-2014) were not characterized by significant alteration in survival rates, human-induced mortalities were the most common reason for wolf deaths in Minnesota, as was the case in adjoining populations in Michigan's Upper Peninsula (O'Neil et al., 2017) and Wisconsin (Stenglein et al., 2015a). Human-caused mortality also showed a positive temporal trend, but we recommend caution in interpreting this owing to the lack of statistical strength. Anthropogenic causes resulted in ∼66% of known mortalities, which includes legal and illegal killing, and vehicular collisions. The road and highway network have expanded in Minnesota over the past 50 years and the state governance has committed to further increase this network till 2040 (Minnesota Transport Alliance, 2011). Vehicular collisions continue to cause wolf deaths in Minnesota, and they likely can increase because of increase in road networks.
Although we did not detect any significant age-or sexspecific differences in annual survival rates of wolves, wolf behavior suggests younger dispersing wolves are often more vulnerable to mortality (Mech and Boitani, 2010). Dispersal events often compel young wolves to navigate around or across risky unfamiliar wolf territories and human-dominated landscapes and roads, thereby reducing their survival chances (Barber-Meyer et al., 2021). More information from younger/dispersing wolves is necessary to confirm such patterns.
While we evaluated extensive and long-term data on collared wolves, we acknowledge that fate could only be confirmed for ∼40% of the collared wolves. Data from the remaining individuals had to be censored owing to the lack of information following endpoints when either the collar failed or contact was lost. It is possible that censoring and mortality could be confounded in some cases, i.e., if censored wolves died and collars failed in such a way that the mortality went undetected (Liberg et al., 2012;Stenglein et al., 2018), but we found no evidence or reason to indicate that to be of concern here. First, the author team included those who work directly with collared wolves within the study area, with experiences indicating that known collar loss and failure has been extremely common among marked individuals. Further, there is little reason to believe that the circumstances (e.g., illegal killing followed by destruction of the transmitter) causing such misclassification are common in our study region; illegal kills are commonly detected as a mortality cause. Second, most prevailing literature from neighboring regions with similar socio-political environments suggests that the percentage of wolves with lost collars that may have confounded the detection of dead wolves has been quite low ( < 1%, Stenglein et al., 2015b), thereby reducing any significant positive survival biases in our interpretations. Finally, in the Bayesian shared frailty model that we have used to analyze survival, the censored endpoint for a given wolf does not depend on a survival assumption beyond the endpoint. Instead, the assumption is that the wolf continues to exhibit the expected survival rate (given the model) beyond the point of censoring. While misclassification of a large number of dead wolves as censored could ultimately result in optimistic survival estimates, all indications suggest this type of error is rare. Having incorporated these checks and balances, we are confident in the survival estimates from our analysis.
Our results provide baseline information on the recent trend in annual survival rate of and cause specific mortality of wolves in Minnesota. These demographic parameters would be helpful to inform policy decisions for wolves in the Great Lakes population. Future research exploring site-specific variability in these demographic parameters can provide spatial contexts to the trends that we have reported here and augment our current understanding of the Great Lakes wolf population. Studies have revealed that spatially varying survival rates can be crucial for prioritizing management actions within a landscape, wherein certain areas can be "riskier" for carnivores (characterized by lower survival) and act as population sinks (Robinson et al., 2008;Stenglein et al., 2015a;O'Neil et al., 2017;Barber-Meyer et al., 2021). Identification of such source-sink dynamics through the characterization of area-specific variation in population parameters is important to successful monitoring and management, especially where the population is or has been increasing in number and distribution, and human-caused mortality is a significant contributor to annual survival. Such analyses would also provide fine-scale patterns of survival in protected and nonprotected areas (e.g., Barber-Meyer et al., 2021), within core wolf areas versus expanding frontiers, and on tribal lands where cultural differences toward wolf recovery manifest. Additionally, an increase in number and distribution of wolves expands the human-wolf interface, thereby increasing risk of conflicts. Hence, we recommend continued monitoring of collared wolves to further investigate temporal and spatial patterns of mortality and survival amidst shifting management authority (e.g., state, federal, tribal), landscape conditions, and public attitudes.

DATA AVAILABILITY STATEMENT
The original contributions presented in this study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
All research protocols were undertaken by the Minnesota Department of Natural Resources (MN DNR) as per their guidelines and regulations.