The Golden Eagle (Aquila chrysaetos) has great religious importance to many indigenous North American peoples, including the Hopi Tribe and Navajo Nation of the southwestern United States. Hopi oral traditions indicate their ancestors harvested nestling Golden Eagles prior to the arrival of Europeans to the region, and this religious practice continues today. Despite contemporary conservation concern for Golden Eagles, no studies have evaluated potential negative effects of religious harvest on populations of this species. We conducted aerial and ground searches for Golden Eagle nesting territories on the Navajo Nation from 1996–2005, and monitored occupancy and reproductive rates of territories in three study areas: one area where Hopi annually harvested eaglets, and two areas without harvest. We analyzed 9 yr of data (1997–2005) using multi-season occupancy models and generalized linear mixed models to test for differences among study areas in occupancy dynamics, and production of early-season and fledging-age nestlings. We found no significant differences in probabilities of occupancy, persistence, or colonization of territories between study areas. Territories in harvest and control areas produced similar numbers of nestlings early in the season; however, significantly fewer (53%) reached fledging age in the harvested area, suggesting collection of nestlings led to locally depressed fledgling production. Given possible declining trends of Golden Eagle populations in the southwestern U.S., we recommend continued monitoring and more intensive demographic studies to better understand the effects of religious harvest on the population of Golden Eagles nesting on the Navajo Nation.
The need to understand effects of harvest on animal populations gave rise to the science of wildlife management in North America (Leopold 1933), and subsequent research across diverse taxa has informed sustainable harvest of many species (Krausman and Leopold 2013). However, effects of harvest on nongame species (including migratory birds other than upland gamebirds and waterfowl) have rarely been studied because legal harvest of these species is generally small and allowed only in special circumstances. This includes raptors, for which limited harvest of some species is permitted for use in falconry (Millsap and Allen 2006) and Native American religious ceremonies (Davidson 2009). Eagles figure prominently in the beliefs and traditions of many Native American tribes (Waldman 2006), and ceremonial harvest of Golden Eagles (Aquila chrysaetos) and Bald Eagles (Haliaeetus leucocephalus) is protected by the American Indian Religious Freedom Act (AIRFA) of 1978 (42 United States Code § 1996). Despite widespread concern for conservation of Golden Eagles in North America (Kochert et al. 2002, U.S. Fish and Wildlife Service [U.S.F.W.S.] 2009), research on falconry take of this species is limited (Millsap and Allen 2006), and no peer-reviewed studies have evaluated effects of religious take on Golden Eagle populations.
Golden Eagles are protected by the Migratory Bird Treaty Act (16 United States Code § 703–712), and receive additional protections with explicit penalties for illegal take or disturbance under the Bald and Golden Eagle Protection Act (16 United States Code § 668–668d). These laws and the published regulations they have promulgated empower the U.S.F.W.S. to permit accidental or intentional lethal take of Golden Eagles provided there is no net loss to Golden Eagle populations (U.S.F.W.S. 2009). Although one assessment suggested Golden Eagles were declining across the western U.S. (Kochert and Steenhof 2002), more recent estimates have shown a slight increasing trend for populations in the northwestern U.S. and a slight decreasing trend in the southwestern U.S. (Millsap et al. 2013). The U.S.F.W.S. recommends management to maintain stability of Golden Eagle populations at the level of Bird Conservation Regions (BCRs). Relevant to this study, the 471,600-km2 Southern Rockies Colorado Plateau BCR (of which our 67,873-km2 Navajo Nation study area makes up 14.4%) was the only BCR with a negative population trend for Golden Eagles during 2007–2012 (Nielson et al. 2014). Apparent negative trends of populations in the southwestern U.S. create a pressing need to understand potential factors affecting survival and reproduction of Golden Eagles in this region.
The Hopi Tribe of northeastern Arizona have practiced ceremonial collection and sacrifice of nestling Golden Eagles for several hundred years (Fewkes 1900). Since 1986, the Hopi Tribe has received U.S.F.W.S. permits to take ≤40 eaglets annually from nests on the Hopi Reservation and adjacent public, private, and sovereign lands in northeastern Arizona, including the western third of the Navajo Nation that surrounds the Hopi Reservation (Fig. 1). The Navajo, who venerate live eagles and use eagle feathers in traditional ceremonies, were required by a 1996 court order to provide permits for the Hopi to take up to 12 eaglets per year from nests on the Navajo Nation. As part of land-settlement negotiations in 2006, the Navajo Nation agreed to increase the permitted take by the Hopi Tribe from Navajo lands by 50% to 18 eaglets (U.S.F.W.S. 2009).
Collection of eaglets directly reduces the reproductive output of affected nesting territories, but the population-level effects of this practice are unknown. Likewise, little is known about the relative importance of religious harvest relative to the other factors that may negatively affect Golden Eagle populations, including natural disturbance (Steenhof et al. 1997), human modification of habitat (Watson 2010), electrocution by power lines (Lehman et al. 2007), collision with wind turbines (Smallwood and Thelander 2008) and vehicles (Franson et al. 1995), environmental contaminants (Langner et al. 2015), shooting, poisoning, and trapping (Kochert et al. 2002). The only published study to analyze effects of falconry harvest on Golden Eagles suggested maximum sustainable yield could be achieved with 31% annual take of nestlings and juveniles, but recommended allowed harvest be capped at 5% (Millsap and Allen 2006). Other studies are equivocal on effects of nestling harvest on raptor populations: marked declines in occupied nesting sites of Golden Eagles in western China were attributed in part to overharvest for falconry (Ma 2013), whereas a study of Prairie Falcons (Falco mexicanus) in Wyoming documented negative effects of nestling harvest on productivity and return rates of breeders, but no likely effect on population size (Conway et al. 1995). Removal of potential future breeders from the population is expected to affect Golden Eagles primarily at local and regional scales, rather than range-wide, because of high natal philopatry in this species. For example, 90% of all banded nestling Golden Eagles encountered as adults were found within 175 km of their natal nest (Millsap et al. 2014). Local reductions in the number of breeding-age Golden Eagles could lead to lower colonization rates of available territories and habitat, with the potential to reduce local population density.
In addition to loss of nestlings, collection of eaglets causes disturbance at nests. Human disturbance of nesting Golden Eagles has been associated with reduced occupancy (Whitfield et al. 2004, Kaisanlahti-Jokimaki et al. 2008) and productivity (Boeker and Ray 1971, Watson 2010, Steenhof et al. 2014), although results vary by type of disturbance and response. Further research is, therefore, necessary to evaluate potential negative effects of nestling collection and associated disturbance on Golden Eagles. The desire to better understand effects of this annual take on the Golden Eagle population nesting on Navajo lands led the Navajo Nation leaders to initiate the study described herein.
Our primary objective was to evaluate effects of religious harvest of eaglets on the Golden Eagle population nesting within the Navajo Nation. We analyzed 9 yr (1997–2005) of data on nesting territory occupancy dynamics, and production of early-season and fledging-age nestlings from one study area where eaglets were harvested annually prior to fledging and two adjacent nearby study areas without harvest. We predicted territories in the harvested and control areas would produce similar numbers of nestlings, but fewer young would reach fledging age in the harvest area because eaglets would be removed from nests. We predicted lower rates of fledgling production would result in lower colonization rates for territories in the harvested study area, with the possibility that occupancy rates would also be locally depressed. Additionally, we predicted territories in the harvested study area would have lower probabilities of occupancy and persistence because breeding adults would be disturbed by collection activities.
Study Area. We conducted our study within the Navajo Nation, an approximately 71,000-km2 semi-autonomous Native American territory in the Four Corners region of the southwestern U.S., of which 60% is in northeastern Arizona, 32% is in northwestern New Mexico, and 8% is in southeastern Utah. Our 67,873-km2 study area comprised the contiguous Navajo Nation, excluding disjunct Navajo land holdings (3797 km2) and the Reservation of the Hopi Tribe (6555 km2), which is completely surrounded by the Navajo Nation in Arizona. Our study area has diverse topography characteristic of Golden Eagle habitat, including mesas, buttes, canyons, forested mountains, and expansive desert shrublands and grasslands. Elevation ranges from >3000 masl in the Chuska Mountains on the Arizona–New Mexico border to 830 masl at the confluence of the Little Colorado and Colorado Rivers on the western boundary of the Navajo Nation with Grand Canyon National Park. Cliffs, many as tall as 200 m, occur extensively in some areas and provided all known Golden Eagle nest sites. Annual precipitation for the region averages 25.4 cm (Western Regional Climate Center 2016).
For this study, we divided the Navajo Nation into three study areas: East (E), Central (C), and West (W; Fig. 1). The W study area (22,241 km2), in which Hopi practitioners collected nestling Golden Eagles annually, encompassed the western third of the Navajo Nation that surrounds the Hopi Reservation, with its eastern boundary defined by U.S. Highway 191, Navajo Route 59, U.S. Highway 160, and Arizona Highway 98. We divided the remaining two-thirds of the Navajo Nation into the C (24,034 km2) and E (21,599 km2) control areas along the natural barrier of the Chuska Mountains on the Arizona–New Mexico border, and excluded outlying lands that were not within the contiguous area of the Navajo Nation. The three study areas were broadly comparable in size, habitat composition, and density of known Golden Eagle nesting territories.
Nest Surveys and Measurement of Reproductive Rates. We located nests of Golden Eagles with fixed-wing aerial surveys of potential Golden Eagle nesting habitat during March of 1996–2005. Most nesting territories (n = 124) were located during initial years of nest surveys (1996–1999), and fewer nesting territories (n = 26) located during later years. A small number of nesting territories that were in the Navajo Natural Heritage Program database prior to the study were found opportunistically (n = 10) or located during rotor-winged aerial surveys of New Mexico in the 1980s (n = 5).
Our terminology followed the definitions of Steenhof and Newton (2007) and McIntyre and Schmidt (2012), and recent recommendations of Steenhof et al. (2017). We defined a Golden Eagle nesting territory as a 3.2-km-radius circular buffer around a nest or nest cluster located during the nest survey, where only one pair of eagles nested in a year. We monitored the territories during 1997–2005, and confirmed territory occupancy from air by presence of (1) incubating Golden Eagle, (2) adult pair on nest or with view of nest cliff, or (3) recent nest repair with or without adult present, indicated only by fresh greenery and excluding evergreen desert plants (e.g., Mormon tea [Ephedra viridis]), and (4) eggs and/or nestlings in a nest; and from ground by any of the above criteria, or observation of territorial behavior (i.e., undulating flight). Territory occupancy surveys used a removal design with ≤2 annual visits. We conducted the first survey each year in March when eagles were incubating (Stahlecker et al. 2010) using fixed-wing aircraft (Cessna 206), and made second visits from the ground in March and April. Collection of eaglets by Hopi practitioners occurred mainly during May and the first two weeks of June. All aerial surveys were conducted by DWS and DGM, and most ground visits by DWS. As in other studies of raptor territory occupancy (e.g., Kochert et al. 1999, Olson et al. 2005, Jiménez-Franco et al. 2011), our sample consisted of nesting territories with a known history of use; “occupancy,” therefore, represented the probability that historical Golden Eagle territories were reoccupied (Wallace et al. 2016), as distinct from its more common definition in the ecological literature as probability of occupancy for randomly selected patches (i.e., “patch occupancy”; MacKenzie et al. 2006).
We assessed productivity at nests where incubation was observed with a fixed-wing aerial survey in early May, followed by a second survey in early June of nests with nestlings observed in May. If status could not be determined from aircraft, we made ground visits to accessible nests. We aged nestlings using plumage descriptions (Watson 2010). To compare the number of nestlings before and after the harvest period, we distinguished two stages of productivity: we defined early-season nestling production as the number of young per occupied territory counted during the early-May survey, and fledgling production (i.e., “productivity” as defined by Steenhof et al. 2017) as number of young per occupied territory that reached approximately 51 d (80% of fledging age) counted during the early-June survey. We used nestling counts instead of binary response variables (e.g., nesting success or nest survival) because harvest may entail removal of only part of a brood, resulting in fewer fledglings, but not nest failure.
Model Sets and Selection. We estimated probabilities of detection (p), initial occupancy (ψ1), extinction (ϵ), and colonization (γ) of Golden Eagle territories using multi-season “dynamic” occupancy models (MacKenzie et al. 2006) in Program MARK (White and Burnham 1999) with the RMark interface (Laake 2013) for Program R (R Core Team 2015). For our study, p represented the probability of successfully detecting occupancy on a single survey of a Golden Eagle territory, given the territory was occupied; ψ1 represented the probability that a territory was occupied during the first year of the study in 1997; ϵ represented the probability that a territory that was occupied during the previous season was not occupied during the current season; and γ represented the probability that a territory that was not occupied during the previous season became occupied in the current season. We derived estimates of annual occupancy (ψt), which represented the probability that a territory was occupied at year t, using the equation of MacKenzie et al. (2006): . We converted ϵ to probability of persistence (ϕ; ϕ = 1 − ϵ) because we thought ϕ, which represented probability that a territory occupied during the previous season remained occupied in the current season, was more biologically appropriate given high territory fidelity in Golden Eagles (Kochert et al. 2002). We excluded newly located territories from the sample for occupancy modeling until the year after they were found to avoid introducing a positive bias in occupancy probability, and included productivity data for all available territories and years.
To avoid fitting large numbers of models for all possible combinations of covariates, we developed models using a stepwise approach (e.g., Olson et al. 2005, Dugger et al. 2011): (1) we considered a suite of models for p, while maintaining a general structure on other parameters (i.e., group + linear time trend). We considered models for p in which detection among years was constant, year-specific, and a function of additive and quadratic trends; and detection within seasons was constant and a function of additive trends. We did not evaluate occasion-specific detection probabilities because they cannot be estimated under a removal sampling design (MacKenzie et al. 2006). We tested for differences in detection between the E, C, and W study areas by comparing models with no effects of study areas, additive effects of study areas, and interactions between study areas and time trends. We ranked models using the Akaike Information Criterion adjusted for small sample size (AICc; Burnham and Anderson 2002), and retained forms from each modeling step with ΔAICc < 2. (2) We used the top model for p to select forms separately for ψ1, ϵ, and γ, while maintaining a general structure on other occupancy parameters. We considered the same among-year time structures and group effects for ϵ and γ as for p, and only group effects for ψ1 because it does not vary through time. (3) We constructed final model sets consisting of all combinations of forms retained for each parameter. (4) We selected best-approximating models that ranked within 2 AICc of the top model and did not contain uninformative parameters (i.e., nested versions of other competitive models with ≥1 additional covariate; Arnold 2010). This modeling approach allowed us to account for the possibility that observational and ecological processes would vary over time within study areas, while including simpler reduced-parameter models that we expected to be more parsimonious in the absence of strong differences between groups.
We modeled counts of early-season and fledging-age nestlings using Poisson generalized linear mixed models with package lme4 (Bates et al. 2015) in Program R. We evaluated models of nestling counts as constant, year-specific, and a function of additive, curvilinear, and quadratic trends; and models with no effects of study areas, additive effects of study areas, and interactions between study areas and time trends. We included a random effect of nesting territory to account for repeated measures of nesting territories among years. We assessed goodness of fit for Poisson regression models by comparing deviance residuals to model degrees of freedom with χ2 tests.
For all analyses, we interpreted models within 2 AICc of the top model and without uninformative parameters (Arnold 2010) to have strong support. We evaluated covariate point estimates (β̂) and their 95% confidence intervals (CI) as indicators of the direction and strength of relationships. We interpreted lack of overlap between coefficient CI and zero as evidence of significant relationships, and when coefficient CI overlapped zero, we interpreted the degree of asymmetrical overlap to indicate the relative strength of weaker relationships (Ramsey and Schafer 2002).
Occupancy Dynamics. From 1997–2005, we documented 150 unique Golden Eagle nesting territories, and monitored occupancy of an average of 121.11 territories per year (range = 57–149, SD = 32.47). Sample size and density of territories was similar among study areas: 50 territories in E study area (2.31 territories/1000 km2), 52 territories in C study area (2.16 territories/1000 km2), and 48 territories in W study area (2.16 territories/1000 km2).
There were two competitive models of occupancy with ΔAICc < 2 (Table 1). We selected the top-ranked model (Table 2) as best-approximating, and dismissed the other model because it was a nested version of the top model with one additional uninformative covariate (Arnold 2010). Although the top model apparently had only moderate weight (wi = 0.48), deleting all models with uninformative parameters from the model set would have resulted in full support for the top model (wi = 1.00). The best-approximating model predicted initial occupancy was high (mean = 0.98, CI = 0.88–1.00) and did not vary by study area (Table 2). Extinction probability (mean = 0.20, range = 0.02–0.32, SD = 0.12) did not differ between study areas and followed a significant quadratic time trend (β̂T = 1.40, CI = 0.33–2.48; β̂TT = −0.11, CI = −0.21 to −0.01; Fig. 2, Table 2). Accordingly, persistence probability (mean = 0.80, range = 0.68–0.98, SD = 0.12) decreased from 2000–2003, then increased slightly (Fig. 2). Colonization probability (mean = 0.26, range = 0.06–0.54, SD = 0.17) did not differ between study areas, and increased from 1997–2005 with a significant additive time trend (β̂T = 0.42, CI = 0.07–0.77; Fig. 2). Annual occupancy probability (mean = 0.74, range = 0.54–0.98, SD = 0.09) did not differ among study areas and declined from 1997–2003, then returned to 2001–2002-levels by the end of the study in 2005 (Fig. 3).
Multi-season models of occupancy for Golden Eagle nesting territories in the Navajo Nation study area, 1997–2005. Models included parameters for probabilities of initial occupancy (ψ1), extinction (ϵ), colonization (γ), and detection (p), as functions of intercepts-only (.), three study areas (study_area), years (year), and linear (T) and quadratic (TT) time trends. Bold text indicates best-approximating model without uninformative parameters. Number of parameters (K), values of the Akaike Information Criterion for small sample sizes (AICc), ΔAICc, and Akaike weight (wi) are shown for each model.
Beta-coefficients from the best-approximating model for occupancy dynamics of Golden Eagle nesting territories in the Navajo Nation study area, 1997–2005. Initial occupancy (ψ1) did not vary by study area, extinction probability (ϵ) followed a quadratic trend (T + TT), and colonization probability (γ) followed an additive trend (T). Detection probability (p) varied by year (year). Shown are model parameter, structure, beta coefficient (β), SE, and lower (LCI) and upper (UCI) 95% CI. Model selection results are included in Table 1.
Detection of Territory Occupancy. Detection probability (mean = 0.61, range = 0.29–0.79, SD = 0.18) varied significantly between years and did not differ among study areas (Fig. 3, Table 2). Although we could not model detection probabilities separately for first and second visits within years, we most often detected Golden Eagles that were associated with occupied nests during the first (aerial) survey visit: 535 (87.3%) of 613 total documented territory occupancies were incubating birds, including 457 (74.6%) detected during the first (aerial) survey occasion, and 78 (12.7%) detected during the second (ground) survey occasion. Only 78 (12.7%) detections were sightings of two eagles or repaired nests, with 32 (5.2%) detected on the first survey, and 46 (7.5%) on the second survey.
Productivity. We monitored productivity of an average of 70.2 territories per year (range = 31–94, SD = 22.1). We found strong support (wi = 0.88) for the top-ranked model of early-season nestling production (Table 3). This model (Table 4) suggested the number of nestlings per occupied territory during the early-May survey (mean = 0.58, range = 0.24–0.80, SD = 0.20) varied significantly between years, but did not differ among study areas (Fig. 4). Although the second-ranked model included an effect of study area, we disregarded this model because it was not competitive (ΔAICc = 4.07), had low support (wi = 0.12), and coefficient CI of study area effects were centered on zero. Non-zero standard deviation for the random effect of site in the top-ranked model (σ = 0.37) suggested its inclusion was warranted, and Pearson's χ2 test provided no evidence of lack of fit (P = 1.00).
Models of early-season nestling production for Golden Eagle territories in the Navajo Nation study area during early-May, 1997–2005. Poisson generalized linear mixed models of nestling counts as functions of intercepts-only (.), fixed-effects of study areas (study_area), years (year), and linear (T), curvilinear (lnT), and quadratic (TT) time trends. All models included a random effect of territory. Interactions are denoted with colons. Bold text indicates best-approximating model without uninformative parameters. Number of parameters (K), values of the Akaike Information Criterion for small sample sizes (AICc), ΔAICc, and Akaike weight (wi) are shown for each model.
Beta-coefficients from the best-approximating model of early-season nestling production for Golden Eagle territories in the Navajo Nation study area, 1997–2005. This Poisson generalized linear mixed model suggested nestling counts varied by year and included a random-effect of nesting territory. Shown are model parameter, beta coefficient (β), SE, and lower (LCI) and upper (UCI) 95% CI. Model selection results are included in Table 3.
There was strong support (wi = 0.85) for the top-ranked model of fledgling production (Table 5). This model (Table 6) suggested the number of nestlings that reached 80% of fledging age per occupied territory varied significantly by year, and was significantly higher in the C (β̂C = 0.76, CI = 0.49–1.12) and E (β̂E = 0.73, CI = 0.46–1.11) study areas than in the W study area. The number of fledglings per occupied territory predicted by this model was >50% lower in the W study area (mean = 0.26, range = 0.12–0.39, SD = 0.09), than the E (mean = 0.54, range = 0.24–0.82, SD = 0.20) and C (mean = 0.56, range = 0.25–0.85, SD = 0.20) study areas, and differences were strongest in 5 of 9 yr when CI of group means did not overlap (Fig. 4). The second-ranked model included an effect of study area and a quadratic time trend that interacted with study area; we disregarded this model because it was not competitive (ΔAICc = 4.10) and had low support (wi = 0.11). Non-zero standard deviation for the random effect of site in the top-ranked model (σ = 0.37) suggested its inclusion was warranted, and Pearson's χ2 test provided no evidence of lack of fit (P = 1.00).
Models of fledgling production for Golden Eagle nesting territories in the Navajo Nation study area, 1997–2005. Poisson generalized linear mixed models of number of nestlings reaching 80% of fledging age as functions of intercepts-only (.), fixed-effects of study areas (study_area), years (year), and linear (T), curvilinear (lnT), and quadratic (TT) time trends. All models included a random effect of territory. Interactions are denoted with colons. Bold text indicates best-approximating model without uninformative parameters. Number of parameters (K), values of the Akaike Information Criterion for small sample sizes (AICc), ΔAICc, Akaike weight (wi), and deviance are shown for each model.
Beta-coefficients from the best-approximating model of fledgling production for Golden Eagle territories in the Navajo Nation study area, 1997–2005. This Poisson generalized linear mixed model suggested the number of nestlings reaching 80% of fledging age varied by year and was higher in the Central and Eastern study areas than the Western study area, and included a random-effect of nesting territory. Shown are model parameter, beta coefficient (β), SE, and lower (LCI) and upper (UCI) 95% CI. Model selection results are included in Table 5.
Territory Occupancy Dynamics. Few other studies have accounted for imperfect detection in models of Golden Eagle territory occupancy. Occupancy probabilities from our study (mean = 0.74, range = 0.54–0.98) were similar to those from a study in Wyoming, which reported mean annual occupancy probabilities of 0.69 and 0.80 that ranged from 0.28–0.99 as a function of nest height (Wallace 2014). Although estimates from studies that did not account for imperfect detection are expected to be biased low (MacKenzie et al. 2006), unadjusted occupancy rates from other long-term studies of Golden Eagles were nonetheless higher than detection-adjusted occupancy probabilities from our study. For example, territory occupancy over 22 yr in Denali National Park, Alaska (mean = 0.87, range = 0.80–0.94; McIntyre and Schmidt 2012), 23 yr in Idaho (mean = 0.90, range = 0.83–1.00; calculated from Table 1 in Steenhof et al. 1997 and M. Kochert pers. comm.), and 37 yr in the Italian Alps (mean = 0.85, range = 0.60–1.00; Fasce et al. 2011) all averaged higher than our study. Although lower territory occupancy rates from our study could reflect resource limitation in a desert environment, comparison of our results with other studies of Golden Eagles is complicated by the substantially larger size of our study area. The studies described above focused on relatively smaller areas chosen for high densities of nesting Golden Eagles, whereas we may have observed lower occupancy rates, in part, because our larger study area incorporated more marginally suitable habitat in which some territories were only sporadically occupied.
This is the first study, to our knowledge, to estimate the occupancy dynamics of persistence and colonization for Golden Eagle territories, although dynamic occupancy models have been used to analyze long-term datasets on other raptors, including Peregrine Falcons (Falco peregrinus; Bruggeman et al. 2016) and Spotted Owls (Strix occidentalis; Olson et al. 2005, Dugger et al. 2011). Persistence and colonization rates did not differ between study areas, suggesting habitat suitability was comparable across the study area. Colonization probabilities that increase from an initial low point can be expected in studies like ours that begin with a sample of recently occupied territories because high initial occupancy leaves few territories available for colonization (Dugger et al. 2011). Apparent declines in occupancy and persistence over the study period may also be artifacts of our sampling design: because we started with a nearly saturated sample, the initial decline and subsequent increase in these rates may actually reflect the return of the sample to an unknown average level of occupancy. This is a common issue that complicates interpretation of results in territory-based occupancy studies of raptors (Dugger et al. 2011, Wallace et al. 2016). Longer-term monitoring is, therefore, necessary to draw conclusions about trends in occupancy dynamics in our study area.
Detectability. Detection probabilities during our study (mean = 0.61, range = 0.29–0.79) were within the range of other fixed-wing aerial surveys: Booms et al. (2010) reported detection probability of 0.68 for at least one Golden Eagle adult at historic nest cliffs in Alaska, and Wallace (2014) reported detection probability of 0.95 for territories with nests in lone trees and 0.38 for territories with nests on cliffs and rock outcrops in Wyoming. Detection rates were higher in years of better reproductive success, which reflected the fact that nesting Golden Eagle pairs were more easily located when they spent more time at known nests within territories, especially when incubating. Significant variation in detection probability among years of our study underscored the importance of accounting for imperfect detection in studies of Golden Eagle occupancy.
Reproductive Rates. The number of nestlings counted early in the brooding period did not vary significantly among study areas, suggesting similar numbers of young were hatched and/or survived until early-May each year at territories across the Navajo Nation. Early-season nestling counts have not been reported by other researchers, so comparisons were not possible. However, the number of young that reached minimum fledging age (the standard measure of productivity for raptors, Steenhof et al. 2017) was significantly lower in the W study area. Average fledging production for the C (mean = 0.56, range = 0.25–0.85) and E (mean = 0.54, range = 0.24–0.82) study areas were within the range reported in long-term studies from Alaska (mean = 0.62, range = 0.05–1.21; McIntyre and Schmidt 2012), Idaho (mean = 0.79, range = 0.34–1.38; Steenhof et al. 1997), and Italy (mean = 0.55, range = 0.27–2.00; Fasce et al. 2011), whereas average fledgling production in the W study area (mean = 0.26, range = 0.12–0.39) was similar to the lowest levels documented in other studies.
Effects of Religious Harvest. Although we cannot directly attribute the entirety of depressed fledgling production in the W study area to religious harvest, the similarity of occupancy rates and early-season nestling production among all three study areas suggests environmental conditions to support nesting by Golden Eagles were comparable across the Navajo Nation during our study. The number of young reaching fledging age in the W study area was consistently lower, but tracked annual changes in the C and E study areas, suggesting fledgling production across the Navajo Nation was influenced by fluctuations in unknown environmental or anthropogenic factors. However, during the four of 9 yr when the difference among study areas was less pronounced, fledging rates in the control areas (C and E) declined toward those of the harvest area (W), rather than productivity in the harvest area increasing into the average range of the control areas. Taken together, these relationships suggest that an additive factor was consistently suppressing fledgling production in the W study area and dampening the wider annual variation seen in the control areas. We therefore suggest that our results provide evidence of religious harvest as the most likely cause of low fledgling production in the W study area.
We used study areas as a proxy for nestling harvest because the secrecy surrounding this religious practice prevented us from reliably documenting the nesting territories and times at which eaglets were collected. Accounts from Navajo wildlife officers and Hopi collection reports indicate a minimum total harvest during our study of 87 of 178 nestlings (49%) in W study area, and no confirmed harvest in the C and E study areas. Known loss of 49% of nestlings to religious harvest is similar to our model-based predictions, which estimated a statistically significant difference of 53% fewer fledglings in the W compared to the C and E study areas. Moreover, estimates from our study represent minimum harvest levels because collection of additional eaglets was known to occur after our final survey flight in some years.
We found no evidence to support our predictions that lower fledging rates and/or disturbance of nests reduced persistence and occupancy of territories in the W study area. This lack of relationships could be explained by immigration of breeding-age adults into the W study area, acclimation of this population to disturbance over hundreds of years of harvest, enhanced subsequent survival, reduced age at first breeding, or other factors beyond the scope of this study.
Study Limitations and Suggestions for Future Research. Our goal was to improve understanding of the effects of an understudied threat on Golden Eagles in a region where populations of this species may be declining. Accordingly, we used the best available methods to analyze a legacy dataset on religious harvest of nestling Golden Eagles within the Navajo Nation. However, the methods used to collect these data resulted in some important limitations on the conclusions that can be drawn from our study. We acknowledge that our scope of inference is limited to the nesting territories in our sample because nest surveys were based on a presumed census of high-quality Golden Eagle habitat, not a design-based sampling scheme. We do not, however, expect sampling bias affected our conclusions on differences among study areas because we used consistent methods to locate nesting territories in all areas. We were not able to collect any data on environmental factors (e.g., prey abundance, weather, vegetation), but acknowledge other researchers have found Golden Eagle reproductive success was correlated with environmental factors, particularly prey abundance (Steenhof et al. 1997, McIntyre and Schmidt 2012). Limited funding and personnel restricted us to relatively few survey occasions per year to document occupancy (≤2 visits), and production of early-season (one visit) and fledging-age nestlings (one visit). More survey occasions would likely have improved detection rates and precision of estimates from occupancy models, and accuracy of production estimates.
Our study was relatively short (9 yr) and changes in demographic rates of long-lived species, like Golden Eagles, may occur over longer time periods. The two longest studies of Golden Eagle populations in North America documented declines in demographic rates: in Idaho, a breeding population studied for four decades declined from 35 to 25 territorial pairs due to fire-induced habitat changes and increased human activity (Kochert et al. 1999, M. Kochert pers. comm.). Declines in breeding success over 23 yr in Alaska could not be tied to conditions on breeding grounds, suggesting that conditions on wintering grounds, which include the southwestern U.S., may have negatively affected breeding condition of returning migrants (McIntyre and Schmidt 2012). Additionally, our ability to interpret trends in occupancy may have been limited by adding territories to our sample over the course of the study. Added territories could have introduced either a positive or negative bias in occupancy rates, depending on favorability of environmental conditions for breeding in years when occupied nests were found. Thus, we suggest future long-term monitoring should focus on a fixed sample of nesting territories.
Although we found no differences in occupancy or early-season nestling production between harvest and control areas, other studies suggest the potential for time-lags in effects of human disturbance on Golden Eagles. For example, Gregory et al. (2003) found that nest interference suppressed fledgling production as many as 5 yr later in Scotland and Whitfield et al. (2004) showed that persecution affected the age-structure of a Golden Eagle population, leading to increased occurrence of mixed-age (i.e., adult and subadult) pairs or subadult and subadult pairs. Lag-effects may be unlikely in our study population, given that Hopi oral history indicates religious harvest has occurred for hundreds of years. However, ongoing harvest pressure may be of greater conservation concern in the context of cumulative effects from other contemporary stressors to Golden Eagle populations (e.g., electrocution, collision, lead poisoning, habitat loss). Although protection of Native American religious rights is an important part of Golden Eagle management objectives of the U.S.F.W.S., regulation of religious harvest is a tool by which managers could affect immediate reductions in mortality of the fledgling age class.
Our study was limited by the data collected to evaluate effects of nestling harvest on territory occupancy dynamics, and nestling and fledging production rates. Future studies should use data from marked individuals and probabilistic sampling to investigate other demographic variables, such as survival of nestlings from harvested and non-harvested nesting territories on the Navajo Nation, effects of harvest on age structure of the regional population, and density of nesting pairs and individuals. Although our sampling methods did not allow us to confirm ages for the majority of territory occupants, aging of nesting eagles should be a part of future monitoring efforts on the Navajo Nation because increasing rates of breeding subadults would indicate a declining floater pool (Hunt 1998, Ferrer et al. 2003, Penteriani et al. 2006). We suggest future studies of long-lived raptors on the Navajo Nation and elsewhere use design-based sampling, survey methods and models to account for imperfect detection, and incorporate environmental factors known to affect demographic rates of raptor populations. Given potential declines of Golden Eagle populations in the Southern Rockies and Colorado Plateau region (Nielson et al. 2014), we suggest continued monitoring of the population nesting on the Navajo Nation is justified and necessary.
This effort was initially funded by the Navajo Nation Council, then by the Navajo Nation Department of Justice. Salaries and expenses of DGM and CSS were also paid by Navajo Natural Heritage Program (NNHP) and Navajo Nation Department of Fish and Wildlife (NNDFW). Funding for analysis and manuscript writing were provided by NNDFW. Gloria Tom was NNDFW Director throughout most of the study and provided steadfast support. NNDFW Wildlife Manager, Jeff Cole, also provided necessary assistance throughout the project. Pilot Lee Rhodes of Gallup Flying Service flew all fixed-wing surveys with safety always the priority. We thank Pat Kennedy of Oregon State University for recommending ZPW join our team. Finally, we thank the Academy of the Environment, University of Nevada, Reno, for paying publication costs.