Age-specific survival rates, causes of death, and allowable take of golden eagles in the western United States.

In the United States, the Bald and Golden Eagle Protection Act prohibits take of golden eagles (Aquila chrysaetos) unless authorized by permit, and stipulates that all permitted take must be sustainable. Golden eagles are unintentionally killed in conjunction with many lawful activities (e.g., electrocution on power poles, collision with wind turbines). Managers who issue permits for incidental take of golden eagles must determine allowable take levels and manage permitted take accordingly. To aid managers in making these decisions in the western United States, we used an integrated population model to obtain estimates of golden eagle vital rates and population size, and then used those estimates in a prescribed take level (PTL) model to estimate the allowable take level. Estimated mean annual survival rates for golden eagles ranged from 0.70 (95% credible interval = 0.66-0.74) for first-year birds to 0.90 (0.88-0.91) for adults. Models suggested a high proportion of adult female golden eagles attempted to breed and breeding pairs fledged a mean of 0.53 (0.39-0.72) young annually. Population size in the coterminous western United States has averaged ~31,800 individuals for several decades, with λ = 1.0 (0.96-1.05). The PTL model estimated a median allowable take limit of ~2227 (708-4182) individuals annually given a management objective of maintaining a stable population. We estimate that take averaged 2572 out of 4373 (59%) deaths annually, based on a representative sample of transmitter-tagged golden eagles. For the subset of golden eagles that were recovered and a cause of death determined, anthropogenic mortality accounted for an average of 74% of deaths after their first year; leading forms of take over all age classes were shooting (~670 per year), collisions (~611), electrocutions (~506), and poisoning (~427). Although observed take overlapped the credible interval of our allowable take estimate and the population overall has been stable, our findings indicate that additional take, unless mitigated for, may not be sustainable. Our analysis demonstrates the utility of the joint application of integrated population and prescribed take level models to management of incidental take of a protected species.

INTRODUCTION

The prescribed take level (PTL) model offers a scientifically credible approach to combine demographic information with policy objectives to estimate sustainable allowable take rates, and the uncertainty around those rates, for a species or group of species subject to harvest (Runge et al., 2009). The PTL approach has been used recently to estimate allowable take levels for game and nongame birds subject to recreational harvest or lethal removal to abate damage to crops, livestock, or property (Johnson et al., 2012; Koneff et al., 2017; Runge et al., 2004; Zimmerman et al., 2019). Management agencies are also called upon to manage take of birds that is incidental to other activities in some instances, and the PTL approach can theoretically be used to assess allowable take levels in these situations as well (Runge et al., 2009). In this paper, we apply the PTL approach to the management of incidental take of golden eagles (Aquila chrysaetos) in the western United States, and in so doing demonstrate utility of the approach in managing incidental take.

Golden eagles are an iconic, apex predator that occurs throughout the northern hemisphere (Katzner, Kochert, et al., 2020). Golden eagles and humans have a complex history, ranging from the eagle's religious significance in many indigenous cultures (Watson, 2010) to relentless persecution in some developed and developing nations (Watson, 2010; Whitfield et al., 2004). Over the past century, concern has grown over the status of golden eagles in many areas, including the United States, where human persecution and habitat loss are thought to have negatively affected the species (Harlow & Bloom, 1988; Scott, 1985). That concern motivated the U.S. Congress to pass the Bald and Golden Eagle Protection Act (Eagle Act) (Millsap et al., 2007) in 1940. The Eagle Act makes take (killing, injuring, disturbing, or possessing) bald (Haliaeetus leucocephalus) or golden eagles or their parts illegal in most situations in the United States, but it also recognizes several legitimate uses of eagles, and authorizes the U.S. Fish and Wildlife Service (Service) to permit the take of eagles in these cases. The Eagle Act explicitly identifies incidental take of eagles as a prohibited act, and the authorization of incidental take is one of the cases where permits may be issued. However, Congress specified that the Service's authority to permit take is contingent on a determination that the take is compatible with the preservation of the eagle species being taken. Thus, the Service is legislatively mandated to make a scientific determination regarding the sustainability of any take of eagles that it authorizes, both individually and cumulatively (U.S. Fish and Wildlife Service, 2016a). How best to make that determination is an important question because mistakes can result in overharvest and population declines, and the scientific credibility of such decisions are often contentious and subject to litigation (Beans, 1996).

The Service revised regulations governing the take of eagles in 2016 (U.S. Fish and Wildlife Service, 2016a), and as part of that process it used PTL models with preliminary demographic information from 1997 to 2014 to estimate the allowable take rate. This assessment suggested that most mortality of golden eagles in older age classes was from anthropogenic causes, and that the estimated amount of take exceeded the allowable take level indicated by the PTL model (U.S. Fish and Wildlife Service, 2016b). Under the demographic theory of life history evolution (Reznick et al., 2002), golden eagles exhibit characteristics of species under strong selection for high adult survival, delayed maturation, and low reproductive effort. Generally, population models for species with these characteristics suggest that adult survival is the most important demographic parameter relative to population trajectory (Hunt et al., 2017; Monzón & Friedenberg, 2018; Whitfield et al., 2004). Consequently, the finding that anthropogenic mortality was potentially driving survival rates of older age classes of golden eagles to unsustainable levels in the United States raised considerable concern (U.S. Fish and Wildlife Service, 2016a).

Motivated by these preliminary findings, we launched an effort to build on the Service's analysis to better answer two critically important questions in the context of incidental take of golden eagles: (1) What is the allowable take of golden eagles in the coterminous western United States (hereafter western United States)? And (2) is the current level of anthropogenic mortality sustainable? Our assessment, reported in this paper, incorporates four major improvements over the prior unpublished analyses. First, we extended the period of study to include the years from 1997 to 2016, allowing incorporation of data from two additional years of intensive golden eagle study in the western United States (Brown et al., 2017). Second, we compiled data from numerous studies across the western United States that deployed Argos, GPS‐Argos, or GPS/GSM transmitters (hereafter collectively referred to as transmitters) on golden eagles. These data, when combined with information from band recoveries, constitute the most extensive data set ever assembled for estimating survival rates of golden eagles and afforded us a unique opportunity to assess the potential influence of transmitters on survival. Third, we launched an aggressive, coordinated effort to recover transmittered eagles that died, and to have the remains of those eagles evaluated by pathologists to determine the cause of death. Transmittered birds provide less biased data on causes of mortality compared to band recoveries and incidental finds of dead individuals because both of the latter are dominated by birds that die in places where they are more likely to be discovered (Schaub & Pradel, 2004). Deaths of transmittered eagles can often be detected regardless of where the death occurs. Finally, for our data analysis, we combined the information on survival and causes of mortality with available information on fecundity and annual population size in an integrated population model (IPM). Estimates of vital rates from the IPM were then used to inform our PTL model. IPMs take advantage of the direct relationship between demographic rates and population size by formally combining the available information on all relevant parameters to obtain more precise estimates of each (Kéry & Schaub, 2012; Schaub & Abadi, 2011). The use of IPMs to help estimate allowable take rates with the PTL approach is a relatively new concept (Gamelon et al., 2021), thus our analysis represents an important contribution to the joint application of these two methods.

METHODS

Study area and timeframe

Data used in this study were collected independently and for various purposes by biologists and researchers across the coterminous United States west of the 100th meridian (Figure 1). Golden eagle population estimates were derived from the Breeding Bird Survey (BBS) and a west‐wide aerial transect survey of golden eagles (WGES). The BBS has been conducted annually in late spring since 1966 throughout the study area (Sauer & Link, 2011). The WGES has been conducted annually during late summer 2006–2016, and covers Bird Conservation Regions (BCR) 9 (Great Basin), 10 (Northern Rockies), 16 (Southern Rockies/Colorado Plateau), and 17 (Badlands and Prairies); collectively, these BCRs represent about 85% of the golden eagle's range in the western United States (Nielson et al., 2014). Data on survival came from uniquely numbered leg bands provided by the U.S. Geological Survey Bird Banding Laboratory (BBL; U.S. Geological Survey Bird Banding Laboratory, 2020) and transmitters that we attached to golden eagles in the study area during 1997–2016. We also incorporated data from golden eagles that were leg‐banded prior to 1997 but found dead and reported (i.e., recovered) during 1997–2016. Band and transmitter deployment did not follow any overall design, but most regions of the western United States were represented. Data on fecundity (i.e., females produced annually per breeding pair) were provided by various monitoring programs across the study area from 1996 to 2014 (U.S. Fish and Wildlife Service, 2016b).

FIGURE 1

Map of (a) golden eagle count locations (dots are Breeding Bird Survey routes, blue lines are aerial transects that were surveyed), (b) transmitter deployment locations, (c) fecundity study areas, and (d) banding locations that provided information or data that we used in the integrated population model for golden eagles in the coterminous western United States, 1997–2016

Data sources

Population information

Population size information for the IPM came from the WGES and the BBS (Millsap et al., 2013) (all data used in the analyses for this article are available online; see Data Availability). The WGES provided late‐summer golden eagle population size estimates for the four BCRs covered by the survey. The WGES was timed to precede most post‐fledging dispersal and large‐scale fall migratory movements of golden eagles (Bedrosian et al., 2018; Katzner, Kochert, et al., 2020; McIntyre et al., 2008; Murphy et al., 2017). Good et al. (2007) and Nielson et al. (2014) provide details on the design and methods used to collect data for the WGES. The BBS was conducted annually by biologists and volunteers at established survey points along standardized roadside routes during the breeding season for most bird species (April–June; Robbins et al., 1986). The BBS provided an index of golden eagle abundance for each route each year, and unlike the WGES data, were available for the entire study area and period (Robbins et al., 1986, Sauer & Link, 2011).

Banding and tagging data

We filtered banding data provided by the BBL to eliminate individuals marked prior to 1997, those marked with auxiliary devices other than leg bands, and those discovered by means other than having been found dead. We supplemented this banding data set with records from golden eagles that had been banded prior to 1997 but were recovered during our study period. We also incorporated data from golden eagles that were transmittered within the study area during the study period. The majority of these were GPS‐Argos or GPS/GSM tags, but 6% were Argos‐only transmitters and thus not GPS capable. These eagles were transmittered by the authors and our collaborators as part of 14 different studies with varying objectives (list of transmittered eagles is available online; see Data Availability), but we collaborated on the recovery and post‐mortem examination of many of the eagles that died. Collectively, these eagles were transmittered widely across the study area, and subsequent movements exposed these eagles to nearly all regions of the western United States where golden eagles occur (Brown et al., 2017). Transmittered golden eagles included in the study were equipped with transmitters weighing from 22 to 116 g. Transmitters were manufactured by Cellular Tracking Technologies (Rio Grande, NJ, USA), Microwave Telemetry, Inc. (Columbia, MD, USA), North Star Science and Technology (Oakton, VA, USA), Sirtrack (Hawkes Bay, NZ), and Telonics, Inc. (Mesa, AZ, USA). Transmitters were attached with Teflon ribbon (Bally Ribbon Mills, Bally, PA, USA) backpack mount configurations (Buehler et al., 1995), except for 24 units that were mounted on tails (Harmata et al., 2018).

Transmitters typically provided multiple daily locations for each golden eagle. When transmitters failed to provide data when expected, we examined the most recently transmitted data to determine if a transmitter might have failed, or if the lack of transmission could have been caused by a mortality. If a transmitter failure was suspected, we usually conducted a search to attempt to recover the tag and confirm the cause of the failure (e.g., removed transmitter). If recovery was not possible immediately, no further data were transmitted, and there was no evidence of a mortality in the available data, we considered the transmitter to have failed. Many of the failures in this category were expected given the age of the transmitters and their expected life of 3–4 years. When we suspected a transmittered eagle had died, we immediately initiated a recovery effort, usually in coordination with a federal, state, or tribal wildlife law enforcement officer. Varying location download frequencies and uncertainty regarding status usually led to a delay of 2–7 days between the time of death and recovery of a transmittered eagle. If an eagle was suspected of having died due to illegal activity, its remains were typically seized as evidence, and a necropsy was later performed at the Service's Forensics Laboratory. If no illegal activity was suspected and the proximate cause of death was not obvious, the remains were typically sent to the Southeastern Cooperative Wildlife Disease Study or the U.S. Geological Survey National Wildlife Health Center for necropsy. Some dead eagles were too decomposed or scavenged for meaningful post‐mortem evaluation; in these cases, we recorded the cause of death as uncertain. For the others, we based our categorization of a given eagle's cause of death on findings provided by the institution that performed the necropsy, or on evidence at the scene when the proximate cause of death was obvious (e.g., vehicle collision). We acknowledge that in some cases where the cause of death was assigned in the field, underlying factors that might have contributed to the death would have been overlooked (e.g., a sick eagle may be more likely to collide with a vehicle). Additionally, some eagles that were necropsied in the early years of our study were not subjected to full toxicology tests. The net result of these two factors is that our cause‐of‐death findings may under‐represent the degree to which poisons and other contaminants (e.g., rodenticides, lead) contributed to golden eagle deaths.

Fecundity information

The Service conducted a meta‐analysis of published data on the number of young fledged per breeding pair of golden eagles in the study area as part of the regulatory rule revision in 2016 (U.S. Fish and Wildlife Service, 2016b). We used the results of this meta‐analysis, adjusted to reflect the number of female young fledged per breeding female per year (assuming a 1:1 sex ratio in fledglings), for the fecundity inputs for our IPM.

Data analysis

Nielson et al. (2014) used the WGES count data and a combination of distance sampling (Thomas et al., 2010) and mark–recapture methods (Borchers et al., 2006) to estimate population sizes and detection functions for golden eagles in each of the four BCRs covered by the WGES. The BBS data were analyzed using a log‐linear hierarchical model that included fixed effects for trend, strata (the intersection of states and BCRs), and year, and random effects of observer and survey route (Sauer et al., 2013; Sauer & Link, 2011). The WGES provided estimates of golden eagle density for each BCR surveyed as well as the total WGES area, whereas the BBS provided an index to golden eagle population size for each route over the entire western United States Millsap et al. (2013); see supplement therein for code for this analysis) integrated the WGES and BBS analyses to extend the temporal and geographic scope of inference possible from the WGES beyond the four BCRs covered by the survey to the entire western United States covered by the BBS. A key component of this integration was the use of data from both surveys from the four BCRs where both surveys were performed to derive a function to scale BBS indices to golden eagle density as estimated by the WGES. We followed the approach described in Millsap et al. (2013) by using the overall median scaling factor for BBS counts from the jointly sampled BCRs to derive annual population size and variance estimates for the entire study area during 1967–2016. Population size estimates for golden eagles for the BCRs sampled by both surveys during 2006–2016 were derived by integrating both data sets, whereas elsewhere and in other years the estimates were derived only from the BBS indices adjusted by scaling factors estimated from overlap BCRs. Because we used the overall median scaling factor derived from the joint surveys across a large portion of the golden eagle's range in the western United States to adjust BBS counts, we propagated all components of the adjustment and the uncertainty into resulting population size estimates through Markov chain Monte Carlo simulations. We used the full BBS time series (1967–2016) to derive population size (N) estimates but truncated the composite N to 1997–2016 in the IPM to conform to the timeframe for that model. The means and standard errors of N for each year were supplied as the observed population information for the IPM. We elected to supply the parameters of the annual population size estimates from the WGES‐BBS model as data rather than run that model within the IPM because each model alone required many hours to converge, and together convergence time was prohibitive. By reading in the population estimates and their standard errors we incorporated all the uncertainty in the population estimates in the IPM more efficiently.

Integrated population model

The IPM we developed for golden eagles in the western United States integrated information from survival analyses with the fecundity and population size information described above to produce joint estimates of latent and unobservable true (as opposed to observed from sampling) annual population size, as well as estimates of all relevant demographic and detection parameters (Figures 2 and 3; Data S1). We used a state‐space modeling approach that combined a process model describing the true, but unobservable, population process with an observation model that linked the true population process with our sample data (Kéry & Schaub, 2012). We used the estimates of population size (N_OBS ) from the integrated WGES‐BBS model as our observed population information to input into the IPM. The core component of our IPM was an underlying process model that linked changes in population size and demographic rates. For this study we used a density‐dependent, post‐breeding, stage‐structured Lefkovitch matrix model (Caswell, 2001) with four stages: first year (Y1), second year (Y2), third year (Y3), and after third year (AY3): We summed age specific estimates and assumed the likelihood for the process model was Parameters with the subscript t indicate demographic rates assumed to vary annually, whereas those without the subscript t were assumed to be constant. The parameter p.breed was the proportion of AY3 female golden eagles that attempted to breed each year. We assumed . The WGES and BBS did not provide information on population age structure, so we partitioned the into the four age classes by estimating the stable stage distribution, using the mean survival (S) and fecundity (F) values from previous studies and then using them as probabilities for a Dirichlet distribution After the initial year, we allowed to be estimated as a parameter in the IPM. We computed estimates of the initial stable stage distribution, as well as parameter elasticity values (the effect of a proportional change in a parameter on the population's finite rate of increase, λ) using the popbio package (Stubben & Milligan, 2007) in Program R (R Core Team, 2015). We specified the observation model as: where N_OBS_SE was the standard error for the population size estimates from the integrated WGES‐BBS model.

FIGURE 2

Graphical representation of the integrated population model for golden eagles in the coterminous western United States

FIGURE 3

Graphical representation of a cause‐of‐death submodel implemented within an integrated population model for golden eagles in the coterminous western United States. The ψ parameters represent transition probabilities for individuals that died, initially into bins for known and unknown causes of death, and for those for which the cause of death could be determined, into subsequent bins representing major mortality factors. These transition probabilities were combined with age‐specific estimates of population size (N) and survival rates (S) to obtain age‐specific estimates of the number of individuals that died annually (D) from each mortality factor

We had annual estimates of golden eagle population size throughout our study area and time frame, but data to estimate survival and fecundity were patchy in time and space. The survival data were too sparse to fit models that considered temporal trends in rates, and preliminary analyses suggested there was more support for a spatially constant model. For fecundity, prior analyses indicated more support for a spatially and temporally constant model (U.S. Fish and Wildlife Service, 2016b). Considering their life history and strong natural selection for high and stable survival, we assumed that annual changes in golden eagle population size were likely mainly due to annual variation in recruitment. Accordingly, we structured the IPM to sample from the mean and variance of the predictive distribution for fecundity as the primary mechanism to reconcile changes in population size from year to year. In our model we allowed only AY3 birds to reproduce, consistent with evidence that most breeding golden eagles at least four years of age (Bloom & Clark, 2001; Steenhof et al., 1983; Whitfield et al., 2004). We considered incorporating some form of reduced fecundity for first‐time breeders (i.e., newly transitioned Y3 to AY3 individuals) in the IPM because evidence suggests such individuals have lower fecundity (Whitfield et al., 2004). We elected not to incorporate this feature directly because doing so would have required speculation as to the degree to which to discount fecundity. Moreover, the distribution we sampled for fecundity already incorporated this information indirectly in proportion to the relative frequency of occurrence of young breeders in the populations in the fecundity meta‐analysis.

Because we had estimates for all the necessary parameters in the IPM, we were able to estimate not only the proportion of deaths attributable to different causes, but the annual number of deaths (N_COD) by cause (i) for each age class (a) using the formula where is the mean annual number of individuals in age‐group a from the IPM, is the mean annual survival rate for age‐group a, and is the probability a dead eagle from age‐group a died from mortality factor i. For the cause‐of‐death analysis, we collapsed the age classes into Y1 and after first year (AY1) given the relatively small sample sizes and finding in U.S. Fish and Wildlife Service (2016b) that most age‐related variation in cause of death was between these age groups.

We used Bayesian models to estimate unknown parameters in all analyses. We used a Markov chain Monte Carlo method implemented in the jagsUI package (Kellner, 2018) in Program R to estimate the posterior distributions of parameters. Unless otherwise noted, we used uninformative priors in our Bayesian models, and employed three chains of 60,000 iterations each, with the first 50,000 discarded as burn‐in. Throughout the paper, we present parameter estimates as the median values of distributions and the 95% credible intervals (CI) unless otherwise noted.

We used three models to estimate S in the IPM: (1) a dead‐recovery model to analyze data from eagles banded and recovered between 1997 and 2016, (2) a known‐fate model to analyze the data from eagles banded prior to 1997 but recovered during 1997–2016, and (3) a combination known‐fate and cause of death model to analyze data from transmittered eagles. For the latter model, we first estimated survival allowing for different age‐specific survival rates for transmittered and band‐only golden eagles, and we compared these estimates within age classes for evidence of a tag‐effect, as recommended by Cooch and White (2014). Following that, we implemented a version of the model that assumed age‐specific survival rates were the same between transmittered and band‐only eagles. We assumed an annual time step for our survival models. Capture and marking, for transmitters and leg bands, occurred throughout the year rather than during short distinct time periods each year. Consequently, annual time steps did not overlap completely among individuals, nor did they follow calendar years. For example, the annual time step for a bird banded on 1 January 2000 ended on 31 December 2000, whereas the annual time step for an eagle banded on 1 December 2000 ended on 30 November 2001. This precluded us from being able to estimate annual survival rates, so we assumed that survival was constant in all models. We believe this is a reasonable assumption for two reasons. First, previous analyses of banding data for golden eagles showed that survival models that allowed for temporal variation were not as competitive as those with constant survival rates (U.S. Fish and Wildlife Service, 2016b). Second, relative stability in survival rates is consistent with observations in trends (i.e., stable population with little year‐to‐year fluctuation; Millsap et al., 2013) and the high sensitivity of population growth rate to survival (i.e., because population growth rate is highly sensitive to survival, we would have expected greater annual variability in population indices if survival varied annually at the scale of our analysis). We also assumed that reporting rates were constant through time because no concerted programs exist to increase or decrease locating and reporting bands through time (e.g., waterfowl harvest and recovery rates vary through time with changes in regulatory seasons for achieving management objectives; U.S. Fish and Wildlife Service, 2019). Therefore, we did not consider temporal variability in survival or reporting rates in this analysis but acknowledge that this assumption is important to assess if more intensive sampling for survival is implemented.

We used the dead recovery model described by Seber (1970) to estimate S , and the pooled reporting rate (r_BAND; the probability that a leg‐banded dead golden eagle in any age class was encountered, and its band number reported to the BBL) for the band recovery data for eagles leg‐banded after 1996. We analyzed the data by creating a matrix where each row represented a cohort of eagles banded and released in a year and the columns represented the number of recoveries from eagles in that cohort for each year of the study. We assumed each row of the matrix was a multinomial random variable: The recovery matrix was indexed by the number of years of banding (t = 1–20 [1997-2016]) and number of years of recoveries plus an extra column for birds never recovered (j = 1–21; the “.” in the above equation indicates the data range is the entire row). The represents the number of eagles banded in year t and Pr represents the cell probabilities (i.e., combinations of and r_BAND representing the distribution of recoveries among years). We assumed uniform priors in the interval of 0 to 1 for all and r_BAND.

For golden eagles banded and released prior to our study but recovered after 1996 we assumed that survival and reporting rates were constant, and that cell probabilities were conditional on birds being recovered during the study. We simulated typical band and dead recovery data and analyzed those data as we did for this paper and found that the approach was unbiased if survival and recovery rates were constant (Appendix S1). Given this model is unbiased with constant survival and recovery probabilities, the assumption of constant survival and recovery are reasonable for this population as described previously, and the additional information provided by the traditional band recovery and telemetry data, we felt confident that these additional band recoveries could help inform survival without biasing estimates. For the purposes of this model, we specified that eagles banded prior to 1997 and recovered subsequently were “first” banded and released in 1997 in whatever age class they were in at that time. In this format, only the survival probabilities of these eagles after 1997 were considered in our analysis. We used a multinomial model to analyze these data, which had to be scaled by the probability of the final cell in the matrix (i.e., the probability an individual is not recovered), because no individuals in this data set were never recovered. The recovery vector, represented the annual recoveries of eagles in age class a that were released prior to 1997 and recovered sometime between 1997 and 2016 (). The cell probabilities were calculated from the same survival probabilities used for the band recovery and transmitter models (i.e., birds banded from 1997 to 2016); we fixed r_BAND = 1 because the sample was limited to known recoveries. We assumed that leg‐band only eagles alive in 1997 had the same survival probabilities whether they were recovered or not. That is, there was some unknown number of banded eagles that were banded prior to 1997 that were still alive in 1997, but which were never recovered and reported subsequently. We assumed the eagles banded prior to 1997 had the same age‐class‐dependent survival probability as birds banded from 1997 to 2016.

For analysis of the known‐fate data from transmittered eagles, we used a multistate model expanded from model 9.5 in Kéry and Schaub (2012) and as described in detail in Millsap et al. (2019). Our model included 20 states: four survival states (Y1, Y2, Y3, and AY3) for working transmitters; four survival states for failed transmitters; two states for eagles recovered dead from unknown causes (one for working transmitters and one for failed transmitters); nine recovery states for each of the cause of death sources described below; and one absorbing state for eagles that had died. We assumed the probability of resighting a live (p_PTT) or recovering a dead (r_PTT) golden eagle with a transmitter was 1, and we considered each transmittered eagle to be in this state unless its transmitter failed. Transmittered golden eagles with failed transmitters transitioned to resight and dead recovery probabilities equal to those of band‐only golden eagles (p_BAND and r_BAND, respectively); we assumed that p_BAND = 0 and r_BAND was the same for eagles with only bands. We adapted the general approach for estimating the relative importance of different sources of mortality developed by Schaub and Pradel (2004) to our transmitter known‐fate data set by incorporating several additional transition parameters: (1) a parameter for the probability cause of death could be determined (); (2) a parameter for the probability that a transmitter failed (); and (3) nine mortality‐specific transition parameters (). Specific causes of death (superscript i) for which we estimated probabilities were: (1) collision (vehicle, power line, wind turbine), (2) electrocution (power line), (3) shooting, (4) poisoning (including lead toxicosis), (5) trap (typically, caught in leg‐hold trap intended for other wildlife), (6) fighting (likely intraspecific aggression), (7) disease, (8) accident (e.g., trauma from a collision not associated with an anthropogenic feature), and (9) starvation (possibly exacerbated by exposure to harsh weather; Figure 3). We were unable to analyze causes of death at a finer scale (e.g., collisions with vehicles versus powerlines) because we did not have sufficient information in most cases to make that determination. As noted earlier, we separately estimated for the two age groups, so a represented either Y1 or AY1. Therefore, the cell probabilities for the transition matrix where transmittered eagles with working transmitters died and causes of death were determined were and for eagles that died for which a cause of death could not be determined, the cell probabilities were Thus, to assign a mortality into one of the nine categories, an eagle had to die () from one of the nine mortality causes (), with a working transmitter (), and the cause of death was determined ().

We considered 20 observation states: four for being detected alive with a working transmitter (Y1, Y2, Y3, or AY3); four for being detected alive with a failed transmitter, which we fixed to a probability of zero for all four age classes because we assumed p_BAND = 0; two for observed dead with mortality source unknown (one for working transmitters and one for failed transmitters); nine for observed dead from one of the nine cause of death factors we considered; and one for not seen alive or found dead. We assumed Dirichlet priors for the to ensure that the probabilities for the nine possible causes of death, given an eagle died and its fate could be determined, summed to 1 (Kéry & Schaub, 2012). We assumed uniform priors between 0 and 1 for and .

For fecundity, we assumed that the production of Y1 eagles was the product of two processes: (1) the proportion of AY3 eagles that reproduced (p.breed) and (2) the number of female juveniles fledged per breeding female (). We assumed that the proportion of adults that bred each year was a binomial random variable and that density dependence largely acted through this parameter in two ways: (1) territoriality and limitations on the number of suitable breeding territories (locations suitable for a pair of golden eagles to build a nest and attempt to reproduce) placed an upper limit on the number of adults that could attempt to breed (the remainder of which were adult floaters) and (2) prey availability determined the proportion of pairs that occupied suitable breeding territories that attempted to breed each year (Ferrer et al., 2004; Hunt, 1998; Whitfield et al., 2004). We used a logit link to correlate population size with the proportion of AY3 birds that bred: To improve convergence, we divided by 1000 for this analysis. Probability of nesting for paired golden eagles has been observed to correlate closely with prey abundance, particularly Lepus spp. and Sylvilagus spp. (Katzner, Kochert, et al., 2020). Previous studies estimated that when prey density is high, between 93% and 97% of golden eagle pairs attempt to breed (Steenhof et al., 1997). Therefore, we assumed that at low population sizes when prey density was not limited approximately 95% of AY3 birds would attempt to breed each year ( = 3 on the logit scale) and treated this parameter as a constant in the model. Although golden eagles can defer breeding when prey populations are low (McIntyre & Schmidt, 2012; Tjemberg, 1983; Watson et al., 1992), for the PTL analysis we were interested in vital rates at low density under average environmental conditions (Runge et al., 2009). We believed a high rate of breeding was likely under these conditions based on empirical data for large species of Accipitriformes (Hunt, 1998; Tjemberg, 1983; Turrin, 2014). We assumed that the probability of breeding would decrease as the density of breeders increased and nesting territories became limiting (Hunt, 1998; i.e.,  < 0). We do not know of any studies that have measured the proportion of golden eagles nesting as a function of golden eagle density. Therefore, we plotted the proportion of birds breeding as a function of observed breeding density using a range of possible values to specify an informative prior (Appendix S2). Based on this assessment, we assumed

The U.S. Fish and Wildlife Service (2016b) conducted a meta‐analysis to estimate the predictive distribution for the number of young fledged per breeding pair of golden eagles using fecundity data from 1997 through 2014 from across North America. Although this meta‐analysis included data on golden eagle fecundity in northeastern Canada and Alaska, these data were sparse compared to those available for the western United States. Additionally, the analysis revealed no evidence of large‐scale geographic differences in the number of young fledged per breeding pair (U.S. Fish and Wildlife Service, 2016b). Thus, we felt comfortable using the mean and variance from this model to specify a fecundity distribution for the western United States, from which we sampled in the IPM. We had to divide the predictive distribution by 2 to convert the number of young fledged per breeding pair to F, the number of female young fledged per breeding pair. The fecundity distribution used in the IPM was which corresponds to an expected fecundity of 0.27 (SD = 0.087).

Estimating allowable take

We used updated parameters from the IPM and the PTL approach applied to other migratory bird species (U.S. Fish and Wildlife Service, 2016b; Johnson et al., 2012; Runge et al., 2009; Zimmerman et al., 2019) to estimate allowable take for golden eagles. The PTL approach requires a management objective (F o), an estimate of the maximum possible growth rate (r max), and an estimate of population size (N) to estimate an allowable PTL with the formula In this formulation, r max/2 is the harvest rate at maximum sustained yield (h max), F o × r max/2 is the desired harvest rate (h), and h N is the desired take limit (H). We selected F o that was consistent with the objective of maintaining a stable or increasing population size relative to the estimated population size in 2009 (U.S. Fish and Wildlife Service, 2016a).

The r max used in PTL analyses should represent the expected growth rate in the absence of density dependence and anthropogenic take under average environmental conditions (Runge et al., 2009). Therefore, r max should be higher than observed growth but realistic in the current environment. We used a four‐stage, post‐breeding projection matrix model to estimate r max based on the estimated survival and recruitment parameters from the IPM: We used simulation to incorporate uncertainty in matrix parameters and estimated r max as the dominant eigenvalue (λ) from the projection matrix minus 1, for each iteration of the simulation. Because r max should represent growth rates at low population density and without prey limitations, we assumed p.breed = 1 (i.e., all AY3 golden eagles attempt to breed), and we sampled fecundity values from a uniform distribution from the median to the maximum observed fecundity reported in U.S. Fish and Wildlife Service (2016b). For survival, we sampled from uniform distributions between the median and upper 95% credible intervals from the IPM results for each age class. We used the most recent (2016) estimate of the latent population size in the IPM state‐space model as N in the PTL. We transformed the mean and standard deviation of N 2016 to the log‐normal scale and sampled values from that distribution to include in the simulations. Thus, we included uncertainty in N in the overall uncertainty in allowable take. Our simulation included 10,000 iterations.

We conducted PTL analyses using two variations on the underlying discrete logistic model, one assuming a linear relationship in density dependence (Runge et al., 2009), and the other allowing for nonlinear density dependence, the theta‐logistic model (Gilpin et al., 1976; Johnson et al., 2012; Zimmerman et al., 2019). The theta‐logistic model is specified as where θ is a shape parameter that describes the form of density dependence. The linear discrete logistic model is a special case of the theta‐logistic where θ = 1.0. We used the regression model fit by Johnson et al. (2012) and the approach described in Zimmerman et al. (2019) to estimate θ for our golden eagle study population.

Determining in our case requires an estimate of carrying capacity. We estimated carrying capacity by assuming population regulation in golden eagles is consistent with Moffat's equilibrium (Hunt, 1998), whereby the primary population regulation mechanism is an upper limit on the number of suitable breeding territories. To apply this concept for golden eagles in our study area, we (1) estimated the maximum number of simultaneously used breeding sites over the period 1997 to 2016 (N sites) and then (2) simulated population growth using our matrix model with N sites limiting the number of AY3 females that could breed in a given year. We simulated population growth at r max but with the constraint imposed by N sites on the proportion of AY3 female breeders until the population reached carrying capacity and leveled (i.e., λ ~ 1.0). Under this approach, population equilibrium is reached due to a declining proportion of eagles breeding offsetting the growth rate from other demographic parameters. We estimated N sites as the maximum number of female fledglings estimated from the IPM in any year divided by mean fecundity (i.e., total number of juvenile females divided by 0.27, the mean number of juvenile females fledged per year per breeding pair).

We calculated the difference between the observed take from the IPM and the prescribed allowable take limits under the linear and nonlinear discrete logistic models to assess how current anthropogenic mortality related to the prescribed take level. For this assessment, we drew a sample from a normal distribution with the mean and standard deviation for annual observed take from the IPM at each iteration of the PTL simulation. We then subtracted the PTL allowable take estimates from the observed take estimate at each iteration. This allowed us to calculate the percentage of times in the simulation that observed take was less than or equal to the prescribed allowable take level under assumptions of linear and nonlinear density dependence.

RESULTS

Survival

We included data from 512 transmittered and 3128 banded golden eagles in our survival analyses (Table 1). Transmitters used in this study were relatively reliable, with an overall 13% (11%–15%) failure rate per year; ~11% of the 13% of transmitters that failed were transmitters that had reached the end of their expected life (i.e., expected failures). Among transmittered eagles that survived the period over which they were monitored, the mean number of days tracked was 585 (n = 336, SD = 487).

TABLE 1

Sample sizes of banded and transmittered golden eagles for survival analyses

Marker type	Deployed (age at deployment)				Recovered a
	First year	Second year	Third year	After third year
Banded after 1997	2656	88	11	254	199
Banded before 1997		14	16	89	119
Transmitters	292	23	4	193	175
Total	2948	125	31	536	493

Note: All bands and transmitters were deployed in the coterminous western United States from 1997 to 2016 (see Figure 1 for a map of banding and tagging locations), except for eagles banded before 1997. Golden eagles in this category were banded prior to 1997 but were recovered between 1997 and 2016 (see Methods for a description of how these band recoveries were analyzed).

Recovered banded eagles were found dead and reported to the U.S. Geological Survey bird Banding Laboratory. We used transmitter data to identify dead transmitter‐tagged eagles, which we recovered.

Annual survival rate estimates for golden eagles tagged with GPS transmitters and golden eagles with only leg bands were nearly identical for all age classes, although estimates from bands were slightly less precise for all but the AY3 age class (Figure 4). The greater precision was not unexpected, as the recovery rate for transmitters was 1.0 and we estimated the dead band recovery/reporting rate as 0.09 (95% CI = 0.08–0.10).

FIGURE 4

Annual age‐specific survival rates of golden eagles in the western United States, 1997–2016, based on band recoveries from dead eagles, recoveries dead transmitter‐tagged eagles, and with both marker types combined. Error bars span the 95% credible intervals

Given the lack of evidence of a transmitter effect on survival, we pooled data from both marker types to obtain annual survival rate estimates in the final IPM. The pooled model‐estimated mean annual survival rates ranged from 0.70 for Y1 birds to 0.90 for birds in the AY3 stage. Annual survival steadily increased with age among the four age classes.

Causes of death

We recovered 175 dead golden eagles that were tagged with transmitters and determined the cause of death for 126 (72%). The observed causes of death were starvation and emaciation (N = 37), shooting (N = 16), collisions (N = 16; five with vehicles, two with wind turbines, two with power lines, one with a train, and six undetermined), accidents (N = 15; four predation, two impacts with natural features, one drowning, one burned in a wildfire, and seven trauma but in natural settings), electrocution (N = 13), poisoning (N = 10; four lead, three multiple substances including lead, two Aldicarb, and one anticoagulant rodenticide [Cholorophacinone and Diphacinone]), disease (N = 8; four West Nile virus, two septicemia, one complications from knemidocoptiasis, and one with multiple issues), intraspecific fighting (N = 7, 6 of which were AY3 individuals), and trapping (N = 4; one Conibear, one snare, and two unspecified traps).

Our model indicated that leading causes of death differed between Y1 and AY1 golden eagles, with natural mortality factors predominant for Y1 (75%) and anthropogenic factors predominant for AY1 (74%) individuals (Table 2). The model‐estimated primary cause of death for Y1 golden eagles was starvation (50%), most cases of which occurred prior to or shortly after dispersal from the natal territory. For AY1 golden eagles, the model indicated most deaths were from shooting (20%), collision (18%), electrocution (14%), and poisoning (13%; Figure 5). We extrapolated from these proportions to estimate that anthropogenic factors (collision, electrocution, shooting, poisoning, and trapping) accounted for approximately 2572 deaths of golden eagles annually in the western United States.

TABLE 2

Estimated annual number of first year (Y1) and older (AY1) golden eagles in the coterminous western United States, and the average number of deaths by cause per year during 1997–2016

	Median	Lower 95% credible interval	Upper 95% credible interval
Y1
N alive at start	4386	3132	6038
Deaths per year
Collision	51	11	143
Electrocution	69	20	174
Shot	69	20	174
Poisoned	32	4	109
Caught in trap	88	30	203
Fight	32	4	109
Disease	88	30	204
Accident	182	86	346
Starvation	656	416	1001
Subtotal take a	333	187	559
Subtotal natural a	997	652	1432
AY1
N alive at start	27,281	23,374	31,779
Deaths per year
Collision	560	322	877
Electrocution	437	231	731
Shot	601	354	926
Poisoned	395	201	675
Caught in trap	191	67	409
Fight	191	68	408
Disease	150	45	351
Accident	274	118	523
Starvation	150	45	348
Subtotal take a	2239	1819	2670
Subtotal natural a	804	520	1160

Note: Cause‐of‐death estimates were derived from transmittered eagles with functioning transmitters (n = 512) that died and were recovered (n = 175), and for which the cause‐of‐death could be confidently determined (n = 126).

Subtotals were estimated as derived parameters in the integrated population model separately from the causes of mortality, so the column totals do not equal the subtotal values and the credible intervals for the subtotals are smaller than the sum of the credible intervals for the individual causes of mortality.

FIGURE 5

Comparison of the probability of different causes of death between first‐year and older golden eagles in the coterminous western United States, 1997–2016. Error bars span the 95% credible intervals. Categories of causes of death are Acci, accident; Coll, collision; Dise, disease; Elect, electrocuted; Fight, intraspecific fight; Pois, poisoned; Shot, shot; Starv, starvation. See Methods for additional details

Recruitment

The model‐estimated proportion of AY3 golden eagles that attempted to breed each year was high and relatively constant (range = 0.85–0.86) and changed little from the prior ( = −0.55 [−0.96–0.04] on the logit scale). During the 20 years of our study, fecundity averaged 0.27 (95% CI = 0.20–0.36), the range in annual median estimates was 0.26–0.27, and there was no evidence of a temporal trend (see IPM summary output; Data Availability).

Population size and projection matrix

Annual IPM estimates of golden eagle population size in the western United States ranged from 31,229–32,257 (range in 95% CIs = 29,254–35,324) during the study period. The IPM annual estimates of the underlying true population size (mean coefficient of variation [] = 4%) were more precise than the observed population size ( = 11%) and both indicated very low temporal variability in population size during the study period (Figure 6). The observed mean estimate of λ during 1997–2016 was 1.00 (95% CI = 0.96–1.05); λ over the most recent 5 years was also 1.00 (95% CI = 0.96–1.05). Elasticity analysis suggested that growth rate was more sensitive to small proportional changes in survival of older birds than of younger birds or of fecundity (elasticity S AY3 = 0.76 vs. elasticity = 0.08 for all other parameters). The stable‐stage distribution based on mean parameter estimates suggested that Y1 birds made up 0.14 of the population, Y2 = 0.10, Y3 = 0.08, and AY3 = 0.68, although we allowed these values to vary annually in the IPM. Given these vital rates, the model estimated a generation time for golden eagles of 12.6 years.

FIGURE 6

Integrated population model estimates of golden eagle population size in the western United States, 1997–2016. The observed counts are annual estimates of population size after breeding in late summer, derived from aerial transect surveys and Breeding Bird Survey counts (see Methods for more details). Shaded polygons are the upper and lower bounds for the 95% credible intervals

Potential take limit

We estimated there was a maximum of 8602 simultaneously used golden eagle breeding sites in the coterminous western United States during our study, and that carrying capacity was ~51,000 individuals. Our PTL simulations for the discrete logistic model with linear density dependence indicated that F o = 0.75 was consistent with the Service's goal of maintaining a stable or increasing population given our estimates of r max (median = 0.11) and N2009 (median = 31,900). With this model we estimated that 1280 golden eagles from our study population could be taken annually while meeting the management objective (Table 3). For the θ‐logistic model, we estimated a median value of θ = 2.50, which implies the effect of density dependence on golden eagles is nonlinear and greatest as populations approach carrying capacity. With the θ‐logistic model, we estimated that 2227 golden eagles could be taken annually while meeting the objective (Figure 7). The percentage of simulations where observed take exceeded allowable take was 99% under the assumption of linear density dependence and 63% with the θ‐logistic model.

TABLE 3

Results of a potential take limit analysis for golden eagles in the coterminous western United States, 1997–2016

Variable	Mean	SD	Lower 95% CL	Median	Upper 95% CL
Population size (2016)	32,269	1543	29,349	32,256	35,415
AY3 survival	0.94	0.02	0.90	0.94	0.98
Y3 survival	0.92	0.02	0.88	0.92	0.96
Y2 survival	0.87	0.02	0.83	0.87	0.91
Y1 survival	0.73	0.02	0.70	0.73	0.77
Fecundity	0.45	0.10	0.28	0.45	0.61
r max	0.11	0.03	0.05	0.11	0.17
h linear	0.04	0.01	0.02	0.04	0.06
H linear	1276	385	522	1280	1994
θ	4.04	4.91	0.35	2.50	17.24
h nonlinear	0.07	0.03	0.02	0.07	0.13
H nonlinear	2283	920	707	2227	4182

Note: Demographic rates are values expected in the absence of density dependence and without anthropogenic take; these were used to estimate r max, the maximum growth rate possible under average environmental conditions. Estimates with the subscript linear were generated by using a discrete logistic model and assuming a linear density dependence; those with the subscript nonlinear were generated by using a discrete logistic model and assuming nonlinear density dependence, with θ = 2.5. Parameters reported are r max, the maximum growth rate; h, the allowable harvest rate; θ, a shape parameter that describes the form of density dependence; and H, the allowable take limit.

FIGURE 7

Estimates of the prescribed take limits (H) for golden eagles in the western United States under assumptions of linear (θ = 1.0) and nonlinear (θ = 2.5) density dependence; error bars are 95% credible intervals. The horizontal line is the estimated median annual actual take; the shaded polygon is the 95% credible interval

DISCUSSION

Population trend

The IPM suggests golden eagle populations are, and have for some time been, stable in the western United States. This finding is important because, unlike recent previous analyses of count data alone (Millsap et al., 2013; Nielson et al., 2014), it takes into account new demographic data, especially the survival information that confirms high rates of anthropogenic AY1 mortality. To some degree, this stability likely reflects shortcomings of the survey methods, particularly the BBS, to capture annual variation in population size. For example, the population estimates based on the BBS data alone (i.e., 1968–2005) show less annual variation than the subsequent observations that include the WGES data. The insensitivity to short‐term temporal variation would not affect the long‐term trend, however. Thus, one key finding from this study is that current levels of mortality averaged across the western United States appear to have so far been compatible with the maintenance of stable golden eagle populations at that large scale.

This conclusion, however, should not be interpreted to mean that golden eagle populations are stable everywhere in the western United States. For example, Scott (1985) documented declines in golden eagle populations in southern California and Watson et al. (2020) provide evidence that declines may be occurring in eastern Washington. Moreover, new sources of mortality would likely be additive, potentially causing a change to the population's trajectory (Katzner et al., 2017).

Survival rates

One concern with survival analyses based on transmitter‐tagged individuals is the possibility that transmitters negatively affect survival. Like Crandall et al. (2019), we found no evidence of a negative effect on golden eagles with transmitters of the variable weights and shapes included in this study. Another potential source of bias is the destruction of transmitters by humans when eagles are illegally killed, but unlike studies of golden eagles in other countries (Whitfield & Fielding, 2017), the transmitter failure rates we observed were consistent with expected normal rates of failure in such devices. We also acknowledge that some forms of mortality, such as vehicle collision or electrocution, could be associated with transmitter failure, resulting in these causes of death being under‐represented in the data. These were among the most frequent causes of mortality that we documented, however, and we do not have any evidence of this occurring widely in our sample. Given this, we believe our data provide relatively unbiased estimates of golden eagle survival and causes of death across the western United States.

The mean AY3 survival rate we estimated for golden eagles in the western United States is lower than that reported by other researchers working with this species in south‐central Montana (0.93), roughly comparable to rates in west‐central California (0.83–0.90, depending on breeding status), but slightly higher than that reported for southwestern Montana (0.86; Crandall et al., 2019, Harmata et al., 2018, Hunt et al., 2017; many of the eagles in these individual studies were also included in the present study). Fewer studies have reported comparable survival rates for other golden eagle age classes, but the rate we estimated for Y1 golden eagles is lower than that reported from west‐central California (0.84) and the southern Colorado Plateau (0.79; Hunt et al., 2017; Murphy et al., 2017). McIntyre et al. (2006) reported FY survival rates for migratory Alaskan golden eagles of 0.19–0.34, but the migratory nature of that population likely contributed to the relatively low survival. The majority of golden eagles in this study were either non‐migratory or relatively short‐distance migrants.

Because we constrained our survival rate estimates to be constant over time, the estimates and uncertainty reported here represent measures of survival pooled over the broad range of environmental and anthropogenic conditions experienced by golden eagles in the western United States during a 20‐year study period. Thus, differences between our results and others’ findings in local areas over shorter time periods are unsurprising. Undoubtedly, there are fine‐scale spatial and temporal differences in survival that reflect the relative preponderance of anthropogenic risks, the abundance of prey, the propensity to undertake long‐distance dispersal, and precipitation cycles that are smoothed over in our analysis (Hunt et al., 2017; Murphy et al., 2017; Crandall et al., 2019; Murphy et al., unpublished manuscript). Our objective in this analysis was to obtain estimates of survival representative of the study area as a whole and that incorporated all of the uncertainty associated with this variability, rather than to specifically model the local or regional variation in survival, because the U.S. Fish and Wildlife Service manages permits for incidental take of golden eagles at the larger scale (U.S. Fish and Wildlife Service, 2016a).

Causes of golden eagle mortality

Humans caused nearly 60% of all golden eagle mortality in the coterminous western United States, and over 70% of AY1 mortality. High rates of anthropogenic mortality have been reported in most other studies of survival in this species (Crandall et al., 2019; Hunt et al., 2017; McIntyre, 2012; Whitfield & Fielding, 2017). In this study, however, we were able to quantitatively discriminate between the relative importance of anthropogenic and natural mortality factors at the population level. Elasticity analysis of the vital rates in the population projection matrix confirmed that adult survival is the most important demographic parameter relative to λ, thus the rate of anthropogenic mortality among adult golden eagles is profoundly relevant to population status.

Russell and Franson (2014) reported on causes of death of 1427 incidentally found golden eagles that were turned into the U.S. Geological Survey National Wildlife Health Center from 1975–2013. Their analysis suggested leading causes of death were trauma and electrocution (each accounting for 27% of cases examined), followed by shooting (14%); emaciation accounted for only 6% of deaths. Our results indicate starvation is a much higher source of golden eagle mortality, particularly for Y1 eagles, than previously reported. We suspect this reflects the recovery bias described by Schaub and Pradel (2004), given that eagles that die both from natural causes in remote areas and from illegal activity are much less likely to be opportunistically encountered and submitted for necropsy than eagles that die accidentally in areas frequented by humans. Because we detected these deaths via transmitters, our study likely provides a more accurate depiction of the relative importance of starvation and some other mortality factors.

Population‐level effect of anthropogenic mortality

The U.S. Fish and Wildlife Service's (2016b) previous allowable take analysis for golden eagles assumed a linear relationship in the effect of density dependence on population growth. Our study shows that the current level of take does not appear to be resulting in population declines. This would be an unlikely outcome under the assumption of linear density dependence, where 99% of simulations of model‐estimated take exceeded the allowable take limit. Results for the θ‐logistic model were more consistent with a stable population because the allowable take limit was only exceeded in 59% of the simulations. A review by Williams (2013) concluded that species with life history strategies favoring high adult survival and relatively low reproductive potential (i.e., K‐selected species like the golden eagle) are particularly likely to demonstrate non‐linear density dependence, with values of θ > 1.0. Thus, based on both empirical and theoretical evidence, we conclude that the discrete θ‐logistic framework likely models allowable take for the golden eagle more reliably than does the linear.

If anthropogenic mortality was entirely additive, the take rate we observed, which may slightly exceed maximum sustainable yield as estimated with the θ‐logistic model, could cause golden eagle populations to decline. In some other long‐lived species of predators for which anthropogenic mortality matches or exceeds natural mortality, take appears to be compensated for, to some degree, by decreases in density‐dependent feedback that lead to increases in other vital rates (Creel & Rotella, 2010; Gantchoff et al., 2020; Lebreton, 2005). Péron (2013) evaluated the capacity for demographic compensation according to life history attributes across a range of taxa and found that such capacity is inversely related to generation time. Our estimate of generation time for golden eagles in the coterminous western United States is close to the maximum considered in Péron's (2013) analysis. Thus, we would expect that the capacity to absorb added mortality, particularly take of adults, would be limited.

To some degree, the discrete θ‐logistic model accounts for this density‐dependent feedback. However, flexibility in age at first breeding is a density dependent response that is not fully accounted for in our IPM and thus is not in our PTL analysis. Plasticity in the age at first breeding is thought to contribute to resilience to high levels of mortality, and it is also established that variation in age at first breeding is a density‐dependent response for many members of the order Accipitriformes (Ferrer et al., 2004, Millsap et al., 2019, Morandini et al., 2019; but see also Monzón & Friedenberg, 2018). As such, this response can provide a buffering mechanism when populations are in decline and competition for nesting territories is relaxed.

Our IPM results suggest the proportion of AY3 golden eagles breeding each year in the coterminous western United States was high and temporally stable, but we did not allow for breeding by younger individuals in our IPM. Yet, golden eagles in subadult plumage (primarily Y3 in this context, Katzner, Kochert, et al., 2020) can and do successfully breed (Hunt et al., 2017; Murphy et al., 2019; Steenhof et al., 1983), likely most often in situations where there are insufficient numbers of AY3 individuals to fill all available breeding opportunities. As relevant examples, Whitfield et al. (2004) demonstrated that the proportion of breeding golden eagles in subadult plumage in Scotland increased in areas where eagles were being persecuted and adult mortality was high. Similarly, Watson et al. (2020) documented high levels of subadult breeding in eastern Washington over a period of declining reproductive success. Thus, one possible means by which golden eagles compensate for high take rates or declining fecundity, at least locally, is through increased rates of subadult breeding. If true, it is unclear whether this is a sustainable strategy given that young accipitriform breeders may have lower fecundity and experience higher post‐breeding mortality (Millsap et al., 2019; Watson et al., 2020; Whitfield et al., 2004).

Conservation implications

Although our findings suggest that golden eagle populations may be resilient to the current level of mortality in the western United States, our model outputs contain high levels of uncertainty in estimates of population size and trend. This uncertainty is important when considering management implications of this study. Additionally, some human activities that contribute to mortality of golden eagles are increasing in frequency across the western United States, such as wind energy development (Beston et al., 2016; Diffendorfer et al., 2019). Thus, the overall rate of anthropogenic golden eagle mortality will likely increase in the future.

The Service has adopted a risk‐averse management strategy for golden eagles that specifies the 20th quantile of the distribution of H (the allowable sustainable take limit given the management objective) as the annual take limit for permitting purposes (U.S. Fish and Wildlife Service, 2016a). This decision acknowledged considerable uncertainty in the relevant data and was intended to reduce the risk of failing to meet the agency's management objective of maintaining a stable population. Based on our PTL analysis, this policy would result in an annual prescribed take limit of 1441 golden eagles in the western United States. The observed annual take that we report exceeds this prescribed take limit by over 1100 individuals annually. Because most of the current take of approximately 2500 individuals annually is illegal, it is not directly under the Service's control and therefore the Service is compelled to give it allocation precedence over new take that would be authorized in permits (U.S. Fish and Wildlife Service, 2016a). For this reason, current Service regulations require compensatory mitigation actions that will “replace” golden eagles authorized to be killed under most new incidental take permits (U.S. Fish and Wildlife Service, 2016a). Our findings indicate this management approach reduces the risk of over‐harvesting and thus is likely to meet the Service's objective (U.S. Fish and Wildlife Service, 2016a). Mitigation to offset golden eagle take authorized under permits has focused on mortality factors for which management solutions exist, such as electrocution abatement (Dwyer et al., 2016; Mojica et al., 2018). Notable work has also occurred on other mortality factors with mitigation potential, such as use of non‐toxic ammunition for hunting and removal of animal carcasses from highways, which could reduce golden eagle mortality from secondary lead poisoning and vehicle collisions, respectively (Allison et al., 2017; Cochrane et al., 2015; Lonsdorf et al., 2018).

A particularly troubling finding from our study is that the apparent leading cause of AY3 golden eagle mortality in the western United States is shooting, accounting for nearly 700 golden eagle deaths annually. This finding is consistent with other recent work documenting high rates of illegal shooting of raptors and other species in Idaho (Katzner, Carlisle, et al., 2020). Purposeful killing of golden eagles is a criminal infraction under the Eagle Act, yet prosecution and associated penalties apparently have been insufficient to curtail widespread shooting. In a recent economic valuation study, the estimated cost of effectively replacing a golden eagle by using the most common mitigation method currently employed under the Service's permit program, retrofitting of power lines to abate electrocution risk, ranges from US$15,200 to US$38,000 per eagle (Hosterman & Lane, 2017). Thus, the economic cost to offset the biological effects of shooting on the golden eagle population in the western United States is theoretically between US$10,275,200 and US$25,688,000 annually. The high biological and socioeconomic cost of shooting of golden eagles suggests that biologists, wildlife law enforcement agents, and prosecutors may wish to focus on actions and initiatives that might better address this mortality factor.

Our modeling exercise identifies three priority areas of uncertainty relevant to management of golden eagles in the coterminous western United States. First, better information on the process of recruitment in golden eagles and the propensity of breeding by AY3 individuals would identify if and how demographic compensation may be occurring in this species. Ultimately, such research could lead to more accurate estimates of allowable take and less need for stringent risk management policies. Second, a better understanding of relationships between golden eagle demography and prey and habitat conditions could expand the suite of offsetting mitigation measures to include habitat enhancement and landscape conservation options. These currently are not practical mitigation alternatives under the Service's permits for take of golden eagles because their demographic benefits cannot yet be quantified. Finally, continuous monitoring of a representative sample of individual golden eagles from annual cohorts with technologies such as GPS tagging would improve identification of temporal trends in anthropogenic mortality and critical vital rates such as survival, age at first breeding, and probability of breeding. Current guidance suggests necessary sample sizes for such a program for golden eagles in the western United States could be as high as 1500 marked individuals annually (Lindberg & Walker, 2007). This number, however, does not take into account the potential for leveraging information from other available data to help inform the estimates of these vital rates. Such data leveraging is a strength of IPMs, and even modest sample sizes, as we have shown here, can provide important insights into demographic processes that can inform management programs.

CONFLICT OF INTEREST

The authors do not have any conflicts of interest.

AUTHOR CONTRIBUTIONS

Brian Millsap, Guthrie Zimmerman—methodology, formal analysis, writing (equal); William Kendall—methodology, formal analysis (equal); all others—capture and tagging, editing (equal).

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Data S1

Click here for additional data file.