Introduction

Golden Eagles (Aquila chrysaetos) are strongly influenced by availability of key prey species, and the availability of prey and foraging habitat may be limiting factors for Golden Eagles throughout their annual cycle (US Fish and Wildlife Service 2016). The distribution and abundance of medium-sized (0.5–4.0 kg) prey species such as hares, rabbits, and ground squirrels can affect breeding success, breeding territory occupancy (Steenhof et al. 1997, Preston et al. 2017), and population-level distribution of Golden Eagles (Schweiger et al. 2015). Habitat modification and changes in land use across large areas of the western United States of America (USA) have raised questions about the status of prey populations (Simes et al. 2015) and the impacts of changes in prey availability to Golden Eagles (Heath and Kochert 2016, Kochert et al. 2020).

Because the Golden Eagle is protected under the Bald and Golden Eagle Protection Act (16 United States Code 668-668d), compensatory mitigation is required to offset authorized (permitted) take (mortality). Golden Eagles are also identified as species of conservation need by various states, tribes, nongovernmental organizations, and other stakeholders interested in conducting habitat-based management (Allison et al. 2017). Authorized mitigation options are currently limited to power pole retrofitting, but compensatory mitigation could potentially be achieved through activities that increase prey availability (US Fish and Wildlife Service 2013). An improved understanding of prey selection by Golden Eagles, the dynamics of prey populations, and the effects of changes in prey communities on eagle populations is needed to support conservation actions and habitat-based compensatory mitigation. However, strategies for prey or habitat-based mitigation are lacking due to the difficulty in quantifying relationships among habitat, prey, and Golden Eagle demography.

Broad patterns in the distribution and abundance of typical Golden Eagle prey species are recognized. A review of Golden Eagle diets in the western United States found that jackrabbits (Lepus spp.) and cottontails (Sylvilagus spp.; Family Leporidae, hereafter leporids) were the primary prey identified at occupied nesting sites in the majority of study areas, but ground squirrels (Otospermophilus spp., Urocitellus spp.), marmots (Marmota spp.), and prairie dogs (Cynomys spp.; Family Sciuridae, hereafter sciurids) were primary prey in study areas where leporids were absent or scarce (Bedrosian et al. 2017). Black-tailed jackrabbits (Lepus californicus) occur in diverse landscapes, including grasslands, shrublands, deserts, agricultural areas, and rangelands across the western USA (Flinders and Chapman 2003, Simes et al. 2015). Vegetation communities that are generally favored by this jackrabbit are dominated by sagebrush (Artemisia spp.), greasewood (Sarcobatus vermiculatus), or shadscale (Atriplex confertifolia; Gross et al. 1974, Flinders and Chapman 2003). Desert cottontails (Sylvilagus audubonii) tend to occur in arid and semi-arid shrublands, shrub-grasslands, and scrub and riparian habitats (Orr 1940, Davis et al. 1975, Chapman and Willner 1978). They also occur in open grasslands, provided scattered shrub cover is available (Lightfoot et al. 2010). Mountain cottontails (S. nuttallii) are often associated with sagebrush habitats, particularly in the northern part of their range (Orr 1940, Powers and Verts 1971, Sullivan et al. 1989). Rock squirrels (Otospermophilus variegatus), the predominant sciurid prey of Golden Eagles in Utah (Bedrosian et al. 2017), are residents of rocky habitats including large boulders, talus slopes, rocky hillsides, and canyons (Oaks et al. 1987, Ortega 1987) and also adapt to disturbed environments, including areas along stone walls and roadside irrigation ditches where they feed from nearby cultivated fields (Marsh 1994). Marmots typically live on vegetated talus slopes or in well-drained rock outcrops in meadows (Frase and Hoffman 1980). Other suitable marmot habitat occurs within a wide variety of rocky areas, such as lava fields, rimrock, cliffs, and canyon walls, within many vegetation types (Bailey 1936). Despite these well-documented qualitative observations among Golden Eagle mammalian prey species and vegetative communities, habitat preferences are poorly understood from a quantitative habitat management perspective.

Golden Eagles in North America demonstrate a preference for jackrabbits, when available (Steenhof and Kochert 1988), as well as the capacity for dietary shifts in response to changing abundance of primary prey species (MacLaren et al. 1988, Steenhof and Kochert 1988, Keller 2015, Watson and Davies 2015, Preston et al. 2017). At a continental scale, Golden Eagles exhibit more diverse diets in systems lacking in leporids, as the frequency of leporids in the diet is negatively related to both overall dietary breadth and the frequencies of sciurids, other mammals, and birds (Bedrosian et al. 2017). Long-term studies of Golden Eagle diet may provide deeper insights into Golden Eagle responses to variability in prey, habitat, and climate. Here, we examine a data set that presents a unique opportunity to explore long-term patterns in Golden Eagle diets. This collection effort represents the largest and longest running data set of its kind, representing 2445 nesting site (i.e., location of the nest substrate on the landscape) visits from 1970 to 2014. The diversity of habitats represented by the 254 nesting sites (Newton and Marquiss 1982, Steenhof and Newton 2007, Steenhof et al. 2017) are likely to exhibit land cover characteristics that could influence Golden Eagle dietary composition (Keller 2015).

Our primary objectives were to (1) identify differences among diets of eagles in different areas and time periods, and (2) identify environmental variables associated with assemblages of Golden Eagle prey. We investigated three spatial questions: (1) do observed prey communities correspond to a priori ecoregion assignments, (2) are environmental variables associated with common prey species, and (3) how does prey diversity vary among nesting sites, ecoregions, and environmental variables? We also explored three temporal questions: (1) what species increase as prey during years when eagles prey less on leporids, (2) does rainfall influence prey use, and (3) does prey use change over time, independently of other time-varying conditions? The answers to these questions should be useful for informing conservation management strategies (including Eagle Conservation Plans and Eagle Take Permits) for eagles in the western USA, particularly in the context of wind development. A better understanding of the relationships between habitat and prey use by Golden Eagles will also be critical to the development of habitat management strategies as a mechanism of eagle mitigation.

Methods

Study Site and Field Methods. Field surveys for breeding Golden Eagles were conducted from 1970–2014 in a large area (approximately 14,500 km²) of north-central Utah, USA, that included portions of Box Elder, Juab, Tooele, Salt Lake, Wasatch, Duchesne, and Sevier Counties. According to the Commission for Environmental Cooperation Ecoregion Level III classifications (Commission for Environmental Cooperation [CEC] 1997), nesting sites were situated in the Central Basin and Range, Wasatch and Uinta Mountains, Colorado Plateaus, and Wyoming Basin ecoregions (Table 1, Fig. 1). The study area is predominantly high desert shrub-steppe containing sagebrush, greasewood, shadscale, and invasive cheatgrass (Bromus tectorum). Areas of the Wasatch front included chaparral and forested habitats. Elevations ranged from 1280 to 2600 masl and the adjacent mountains were >3600 masl.

Table 1.

Summaries of individual prey remains data identified in Golden Eagle nesting sites in north-central Utah, USA, by the Commission for Environmental Cooperation (CEC) Level III Ecoregion classifications (CEC 1997).

Figure 1.

Golden Eagle nesting sites (n = 254) in the long-term study area in Utah, USA, with CEC Level III Ecoregions (CEC 1997).

We haphazardly sampled eagle nesting sites for prey remains (i.e., entered nests either while banding young eagles, or following the completion of breeding activities for that year). The mean number of visits to each nesting site per year was 1.52 visits. We identified prey remains in situ at nesting sites. The minimum number of individual animals represented by the remains was estimated for each species identified. Hereafter, we refer to observed prey remains identified at nesting sites simply as “prey,” recognizing they are likely biased samples of the actual prey captured and consumed. All prey was removed from the nesting site at each visit so that the same prey could not be counted again during subsequent nest visits.

Data Preparation. Prior to analysis, we pooled ecologically similar prey that were rare in our sample (<0.1%), with exceptions for species that were of particular interest. For birds, we retained all distinct species of raptors, owls, gamebirds, waterfowl, and the American Crow (Corvus brachyrhynchos), American White Pelican (Pelecanus erythrorhynchos), Common Loon (Gavia immer), Great Blue Heron (Ardea herodias), and Double-crested Cormorant (Phalacrocorax auritus). We pooled other water-associated birds in two size classes, small (270–416 g) and medium (487–810 g). We binned the remaining rarely observed birds according to average body mass (<86 g, 104–200 g; Dunning 2007). For mammals, we retained badger (Taxidea taxus), coyote (Canis latrans), domesticated mammal species, Rocky Mountain elk (Cervus elaphus nelsoni), gray fox (Urocyon cinereoargenteus), kit fox (Vulpes macrotis), long-tailed weasel (Mustela frenata), mink (Neovison vison), porcupine (Erethizon dorsatum), raccoon (Procyon lotor), ringtail (Bassariscus astutus), and white-tailed jackrabbit (Lepus townsendii). We combined the remaining mammals into the following groups: rare sciurids; woodrats; pocket gophers and voles; and kangaroo rats, pocket mice, and grasshopper mice. For reptiles, only gopher snakes (Pituophis catenifer) were retained and the remaining snake species were pooled. All fish were pooled.

Potential covariates were generated from remote-sensed or climate data. Environmental variables were derived by sampling remote-sensed land cover (National Land Cover Database [NLCD]; US Geological Survey 2020) and elevation (National Elevation Dataset 2015) databases. We calculated the means of elevation and percent cover by land cover class in a 6.4-km-radius (128.7-km² area) neighborhood around nesting sites to approximate the scale of likely foraging habitats (metadata available in  S1 Table (rapt-55-01-06_s02.pdf) in Dunk et al. 2019). Golden Eagles could have used more than one nest structure within a nesting site: here, we selected a single geographic point that represented the focal point of eagle nesting activities for that particular territory. In addition to the more generalized land cover categories of cultivated cropland, forest, grass, shrub, and wetland (Homer et al. 2015), we also considered land cover of sagebrush (hereafter “sage”; LANDFIRE 2010), alfalfa (Medicago sativa; CropScape 2013), and pinyon pine (Pinus spp., hereafter “pinyon”; Fire and Invasive Assessment Team 2015) that may be more relevant to potential prey populations in this region. Preliminary investigations suggested that land cover within 6.4 km of eagle nesting sites did not change much between 2001 and 2016, and that the prey remains found at sites with greatest land cover changes did not differ much from sites with less land cover change ( Supplemental Materials 1 (rapt-55-01-06_s01.pdf)). Hence, we used a single year's land cover estimate for each nesting site to represent conditions for the entire study period. We characterized annual precipitation by finding the mean of annual rainfall across Box Elder, Juab, Tooele, Salt Lake, Wasatch, Duchesne, and Sevier Counties (PRISM Climate Group 2015). We then classified each year as “average,” “dry,” or “wet” depending on whether the average precipitation in that year was within the 98.5% CI of the overall mean or not. We performed analyses in software environment R (R Core Team 2015), and processed geospatial data in ArcGIS (ESRI 2017).

Multivariate Analysis: Nonmetric Multidimensional Scaling (NMDS). We investigated broad patterns in prey species without considering environmental information by assessing community structure via nonparametric ordination. In the multivariate analyses (because we expected to examine the association with ecoregional locations post hoc) we did not consider prey from the three eagle nesting sites in the Wyoming Basin ecoregion because of the small sample size. We summed the prey collected from every visit to a territory because of the many zero values in the full data set (final n=251). We selected nonmetric multidimensional scaling (NMDS) as our method for unconstrained ordination because the technique will map observed community dissimilarities nonlinearly into ordination space, and will accommodate nonlinear species responses of any shape (Oksanen et al. 2015). We opted to transform our counts of prey using two different dissimilarity/distance indices, the Chao index and Anderson's log base 10 version of the Gower index (Anderson et al. 2006), and perform ordination on both resulting distance matrices because of the characteristics of our data, rather than rely on the more typical Bray-Curtis dissimilarity transformation. The Chao index should be robust to differing sample sizes, which was desirable given that we knew nesting sites were visited a variable number of times across years (1 to 68 visits per nesting site; Anderson and Millar 2004, Chao et al. 2005). We also noted that the sample was numerically dominated by black-tailed jackrabbit remains, so we also performed ordination using the modified Gower dissimilarity measure parameterized by log base 10, which should reduce the influence of the relatively high abundances of black-tailed jackrabbits and instead emphasize species composition of the prey at the territory (Anderson et al. 2006).

We noted the stress returned by each ordination and examined Shepard plots of the ordination fits to guide our choice of specifying ordination along two or three axes. The NMDS ordination results depended only on the dissimilarity matrices of the prey collections, but we explored any patterns that might correlate with ecoregion or environmental variables by fitting the corresponding factors and vectors to the ordination results. The correlation between these fits and the ordination results, as well as the visualization of the centroids of ecoregional groups or surfaces of the environmental variables, illustrated possible ecoregional or habitat influences on the prey communities. We also compared the strength of the relative potential explanatory power of each type of variable using this approach.

Multivariate Analysis: Constrained Ordination and Permutational MANOVA. We took two approaches to more directly assess whether the relationship between prey species collections at nesting sites was more strongly influenced by either ecoregional classification or environmental variables, and whether prey differed with precipitation. Our first method was to perform canonical correspondence analyses (CCA), and the second was to examine the variance/dispersion and location within dissimilarity space (by testing the homogeneity of multivariate dispersion and partitioning the sums of squares of dissimilarity matrices via permutational MANOVA, details below; McArdle and Anderson 2001, Oksanen 2005, Anderson et al. 2006). As with the NMDS, we excluded the prey from our three Wyoming Basin ecoregion nesting sites. We also either summed prey for all visits to a particular nesting site (as for the NMDS), or summed prey by both site and annual precipitation (final n = 560).

We compared the results of separate CCA models that assessed either the influence of ecoregion, or the influences of our continuous environmental variables. Although CCA can be performed using collinear environmental variables without drastic effects on the final ordination results, we chose to find parsimonius models by performing forward stepwise selection of our continuous environmental variables based on the CCA analogs of deviance and Aikaike Information Criterion (AIC) (functions step, add1.default, and drop1default in R package vegan). We additionally explored the influence of precipitation on prey community structure by treating annual precipitation (average, dry, or wet) as a categorical variable, as well as the number of nest visits, by assessing the influence of using each of those variables as a conditioning variable (partial canonical correspondence analysis or pCCA).

We were interested in quantifying any heterogeneity in beta diversity (considered here as the variability in species composition among sampling units for a given area) between groups in general, and also note that heterogeneity between groups can complicate the interpretation of permutational MANOVA, which tests for differences in group means (i.e., the location in multivariate space) or the effect of continuous variables. Therefore, we first determined whether the groups of prey as defined by ecoregions or precipitation year differed in beta diversity (as represented by heterogeneity in multivariate space) using a multivariate extension of Levene's test of the equality of variances (i.e., function betadisper in R package vegan) on the dissimilarity matrices created earlier (see NMDS section above; Anderson 2006, Anderson et al. 2006). Significant results were explored by post hoc Tukey honest significant difference (HSD) tests, which created a set of confidence intervals on the differences between the means of the groups (Yandell 1997). Then we performed permutational multivariate analysis of variance on the distance matrices, testing the possible influences of ecoregion, type of precipitation year, and continuous habitat variables on prey community structure (i.e., function adonis in R package vegan). We assessed those habitat variables that had been identified as important by the CCA ordinations, and further examined their marginal significance via permutation testing.

Univariate Regression Models. We explored potential influences on the abundance of black-tailed jackrabbits, cottontails, rock squirrels, and yellow-bellied marmots (Marmota flaviventris), i.e., the four most numerous prey species in the prey data set, using Bayesian hierarchical Poisson generalized linear mixed models accounting for overdispersion (Kéry 2010). We assessed whether habitat characteristics of Golden Eagle nesting sites, climate and temporal trends, or the diversity of prey or co-occurrence of jackrabbits (for the three non-jackrabbit species) affected the number of important prey. Interannual trends were explored by testing linear and quadratic effects of time. Climate was assessed by testing for influence of the precipitation status of the year for that nesting attempt (average, dry, or wet, based on state mean precipitation, see above). We additionally considered the possibility of a lag effect of precipitation, by also considering the previous year's precipitation status. We quantified the diversity of prey by calculating the Hill number (a linear assessment of diversity, calculated here as the exponent of the Shannon-Weaver diversity index; Jost 2007, Marion et al. 2015) considering all species of prey for that prey record. With two different methods, we tested whether the abundance of jackrabbits in each nesting site's prey influenced the abundance of the other three species (either assessing the influence of the number of jackrabbit remains found at that particular nesting site or by classifying each year into “average,” “fewer,” or “more” jackrabbit remains recovered based on deviations from the 0.985 confidence interval of the overall mean number of jackrabbit remains recovered). Our intent in classifying years by jackrabbit abundance was to assess whether prey composition switched in years when the overall use of jackrabbits was lower or higher than average. The overall use and availability of jackrabbits may not adequately represent changes in counts at individual nesting sites, due to parental behavior and individual prey preferences. Because our two jackrabbit variables were related to each other, we did not include both terms in the same model, and instead separately analyzed two series of models for each species and then compared the best models of each series afterwards.

We also assessed the possible influence of habitat variables for each site, after checking for collinearity between habitat measures and excluding categories with Pearson's pairwise correlation coefficients > 0.5. Thus, the remaining habitat variables included proportion land cover of crop, forest, grass, and sage. Grass cover and forest were moderately correlated (Pearson's correlation coefficient r = –0.50), so models were initially built with either grass or forest variables, and then compared, and the most influential of those two variables was retained for further analyses. All continuous predictor variables were Z-standardized (means transformed to zero with unit variance). An offset term for the number of nest visits was included to account for sampling effort, as well as a random effect of site identification to control for repeated measures at the same nesting site. Global models, with all covariates of interest, were fitted first, then uninformative parameters were removed and models were refit. The best models were determined by lowest deviance information criteria and with 95% Bayesian confidence intervals of all beta coefficient estimates not overlapping zero.

All prior probabilities were selected to be uninformative. In most cases, for each model we ran two chains for 14,000 iterations, discarding the first 12,000 runs as adaptation and burn-in, and then drawing 2000 samples for parameter estimates. Convergence was assessed by visual inspection of the model trace, and the Raftery-Lewis and Gelman-Rubin diagnostic tests implemented in the R package coda (Gelman and Rubin 1992, Plummer et al. 2006). If the Gelman-Rubin statistic for any parameter exceeded 1.1, then the model was run for an additional 20,000 iterations and convergence reassessed. We assessed model fit by posterior predictive checks, or Bayesian P-values (Hobbs and Hooten 2015). Briefly, we simulated a new data set at each iteration of the converged chains. We then computed the proportion of times that the estimate of mean response or the sum of squares of model fit of the simulated data set exceeded the corresponding estimate from the observed data set. If this statistic was close to 0 or 1, we concluded that the model did not represent the distribution of the data well.

Results

We sampled 254 nesting sites with 2445 nesting site visits (range = 1–68 visits/site). We identified a minimum of 26,734 individual prey animals representing 147 prey species (Table 1; summarized by species in  Supplemental Materials 2 (rapt-55-01-06_s02.pdf):  Table S1 (rapt-55-01-06_s02.pdf)). The top four most frequently observed prey were black-tailed jackrabbits (67.3%), cottontail spp. (9.2%), rock squirrels (3.5%), and yellow-bellied marmots (2.7%). Prey varied from 1 to 134 individual remains per nesting site, per year (mean = 16.7 and median = 11). The range of values representing land cover overlapped among ecoregions (Fig. 2). Therefore, it is possible that CEC ecoregions are inadequate boundaries as a priori ecosystem assignments for eagle nesting sites.

Figure 2.

The distribution of environmental variables of land cover in a 6.4-km neighborhood surrounding Golden Eagle nesting sites (n = 254) by ecoregion in Utah, USA. Excludes data from three nesting sites in the Wyoming Basin ecoregion.

Multivariate Analysis: Nonmetric Multidimensional Scaling (NMDS). We found that the collections of prey were best described by NMDS when three axes were used. The modified Gower index resulted in slightly lower stress and ordination fit, but the final stress values and fit were still somewhat comparable (for Chao index, stress = 0.166, nonmetric fit r² = 0.972, linear fit r² = 0.833; for modified Gower index, stress = 0.139, nonmetric fit r² = 0.981, linear fit r² = 0.877). Prey clustering towards the extremes of axes represented suites of nesting sites with similar prey (Fig. 3A). These prey groups roughly correspond to species typical of sagebrush-steppe, wetlands/agriculture, and montane ecosystems (Fig. 3B).

Figure 3.

First two axes of nonmetric multidimensional scaling (NMDS) ordination of prey remains collected from Golden Eagle nesting sites (dissimilarities calculated using Anderson's 2006 modification of the Gower index) in Utah, USA, with gray circles indicating nesting sites in multivariate space and (A) notable prey species labeled, or (B) environmental gradients overlaid, but not informing the ordination. Environmental covariates are proportion of land cover in a 6.4-km radius around each eagle nesting site (centered to 0 and standardized to unit variance).

When we fitted either the ecoregional classification of each nesting site (as a factor) or the land cover habitat variables (continuous variables of forest, crop, sage, pinyon, and wetland land cover, as suggested by CCA results below), we found that both approaches correlated with the NMDS ordinations. However, the ecoregion factors had half or less explanatory value than did the most highly correlated continuous variables (e.g., r² for forest land cover of NMDS axes 1 and 2, Gower dissimilarity matrix = 0.5335, vs. r² for ecoregion = 0.2192). Crop cover, which did not vary between ecoregions, correlated with all three NMDS axes.

Multivariate Analysis: Constrained Ordination. We found that continuous environmental variables could explain 18% of observed trends in the prey data sets when constrained ordinations were performed using CCA, vs. 10% explained by ecoregion (Fig. 4A). Six continuous variables, including the proportion of forest, crop, pinyon, wetland, alfalfa, and sage land cover, along with elevation, contributed to the ordination (Fig. 4B). Permutation tests suggested that the first four constrained axes had significant explanatory power. Most of the pattern on the first axis was driven by differences between sage and pinyon vs. forest and elevation. The second axis was differentiated by crops, wetlands, and alfalfa. Our explorations with adding the number of nest visits as a conditioning covariate in a partial CCA yielded similar results with very little explanatory power coming from the conditioning variable (<1%). Likewise, we found that precipitation year did not inform the ordination (<1%).

Figure 4.

Relationships between ecoregion (A) or environmental gradients (B) and prey remains assemblages at Golden Eagle nesting sites in Utah, USA, as illustrated by canonical correspondence analysis (CCA). Eagle nesting sites in multivariate space shown as gray open circles. (A) Centroids of ecoregions (CBR = Central Basin and Range, WUM = Wasatch and Uinta Mountains, CP = Colorado Plateaus) in black. (B) Environmental gradient vectors illustrated as proportionally scaled black arrows. Environmental covariates are proportion of land cover in a 6.4-km radius around each eagle nesting site (centered to 0 and standardized to unit variance).

Multivariate Analysis: Multivariate Levene's test. Beta diversity, quantified as differences in dispersion/variance among groups in dissimilarity space, varied in relation to ecoregion but not precipitation. The pattern of heterogeneity was different depending on how dissimilarity was calculated. Prey at nesting sites within the Wasatch and Uinta Mountains were more heterogeneous than the other two ecoregions when dissimilarity was estimated by the Chao index (overall ANOVA, F = 19.946, P < 0.001; P-adjusted from Tukey test for mountains vs. basin and range < 0.001, for mountains vs. Colorado plateau = 0.012). In contrast, the Wasatch and Uinta Mountains were significantly more heterogeneous than only the Colorado plateau when the modified Gower dissimilarity index was used, which should be less influenced by the abundance of common species such as jackrabbits (overall ANOVA, F = 4.116, P = 0.017; P-adjusted for mountains vs. Colorado plateau = 0.013).

Multivariate Analysis: Permutational MANOVA. Because the dispersion of prey communities did differ between ecoregions, our analyses with permutational MANOVA must be interpreted with caution. However, permutational MANOVA is relatively robust to mild violations of the assumption of homogeneity of variance among groups at larger sample sizes, as we analyze here (Anderson and Walsh 2013). Again, results varied slightly depending on dissimilarity index used. For Chao index, ecoregion by itself influenced prey community composition when tested by permutational MANOVA (F = 51.543, r² = 0.29362, P = 0.001). Alternatively, the continuous variables of forest, pinyon, wetland, and alfalfa all influenced prey community composition (Table 2). Precipitation year did not affect prey community composition. For modified Gower index, ecoregion by itself influenced prey community composition when tested by permutational MANOVA (F= 13.412, r² = 0.0976, P= 0.001). With the relative influence of jackrabbit abundance reduced by the choice of diversity index, a larger collection of continuous variables influenced prey community composition (forest, crop, alfalfa, pinyon, sage, and wetland, Table 2). Also, precipitation year influenced prey composition, although the grouping variable had very little explanatory power (F = 2.0101, r² = 0.00712, P = 0.006). Post hoc comparisons suggested that the wet years showed prey communities that differed from average and dry years.

Table 2.

Land cover differences related to contrasting prey remains assemblages at Golden Eagle nesting sites in Utah, USA, as indicated by permutational multiple analysis of variance (PERMANOVA) analysis. Comparisons were based on either Chao or the modified log base-10 Gower indices. Shown are degrees of freedom (df), sums of squares (SS), F-ratios (F) and the Monte Carlo asymptotic P values (P [MC]).

Univariate Regression Models. Habitat and diversity covariates generally influenced the abundance of individuals in the top four prey species more strongly (as evidenced by larger standardized beta coefficients) than climate or temporal trends (Fig. 5). Both jackrabbits and cottontails were more abundant at low levels of forest land cover ( Supplemental Materials 2 (rapt-55-01-06_s02.pdf):  Fig. S1 (rapt-55-01-06_s02.pdf) and  S2 (rapt-55-01-06_s02.pdf)). In contrast, the abundance of rock squirrels and marmots increased as forest cover increased ( Fig. S3 (rapt-55-01-06_s02.pdf) and  S4 (rapt-55-01-06_s02.pdf)). Jackrabbit and rock squirrel abundances were also influenced by grass cover, with jackrabbit abundance maximized at an intermediate level of grass cover, and rock squirrel abundance maximized at low levels of grass cover. Cottontails increased with sage cover, whereas rock squirrels decreased. The only species that responded to crop cover was rock squirrels, which increased in abundance as crop cover increased.

Figure 5.

Comparison of standardized beta parameter estimates (means ± 95 Bayesian credible intervals) for Poisson models explaining the abundance of black-tailed jackrabbits (Jackrabbits), cottontail spp. (Cottontails), rock squirrels (Rock squirrels), and yellow-bellied marmots (Marmots) in prey remains recovered from Golden Eagle nesting sites in Utah, USA, from 1970–2014. Parameter estimates are shown only if that variable is included in the top model for that species. “Crop,” “grass,” “sage,” and “forest” refer to the proportion of land cover dominated by that vegetation type within 6.4 km of the nesting site, “diversity” is the Hill number estimated to describe the overall diversity of species found within the nesting site's prey remains (larger numbers indicate higher diversity), “jackrabbits” is the minimum number of black-tailed jackrabbit remains found at each nest, “year” is the overall time trend, and “dry year,” “wet year,” “previous year dry,” and “previous year wet” all indicate whether the annual precipitation averaged across weather stations within the study area was outside the 0.985 CI of the mean for that year or the previous one.

We found mixed evidence supporting a feeding strategy in which other primary prey species were the most important when jackrabbits were not heavily used (Fig. 6). Contrary to our expectations, none of the other three major species were less abundant when jackrabbit abundance increased. We found no relationship between the abundance of non-jackrabbit prey and whether average numbers, few, or many jackrabbits were found in that year. Additionally, cottontails and rock squirrels were more abundant in a nesting site's prey remains when jackrabbits were also abundant. However, overall diversity in prey influenced abundance, with jackrabbit abundance decreasing when the overall prey diversity increased, whereas numbers of the other three species increased up to intermediate levels of overall prey diversity.

Figure 6.

Annual frequencies of black-tailed jackrabbit (circles) and cottontails (crosses) found in the prey remains (mean of minimum number of individuals) in Golden Eagle nesting sites in Utah, USA.

The abundance of all four species was influenced by overall precipitation patterns, albeit weakly. More and larger effects were seen for precipitation patterns in the same year as the nesting attempt vs. precipitation in the previous year. The two leporids, jackrabbits and cottontails, were more abundant in dry years, and jackrabbits were also more abundant when the previous year was dry. Rock squirrels and marmots were less abundant in wet years. Time also influenced the abundance of two species in the prey remains, with cottontail numbers declining over time, and marmot numbers at a maximum during the middle of the study, but effect sizes were weak.

Discussion

The abundance of key prey species for breeding Golden Eagles was influenced by habitat within nesting sites, and we successfully predicted broad patterns of prey at eagle nesting sites using land cover. Prey communities were better predicted by environmental variables than by ecoregional location, and our multivariate analyses provided evidence for three groupings of prey. Our results highlighted the likely importance of grass to black-tailed jackrabbits, sagebrush to cottontails, crops to rock squirrels, and forest to yellow-bellied marmots. Given the importance of jackrabbits and cottontails as a breeding season food source for Golden Eagles in the Great Basin (Bedrosian et al. 2017), our results suggest a possible opportunity to effect change in this ecosystem to the benefit of eagles through the conservation of sagebrush-steppe landscapes.

Environmental covariates that were top predictors of prey species abundance were concordant with published descriptions of prey species' habitat associations, suggesting our quantifications of these relationships were grounded in real ecological processes and were not statistical artifacts. For example, the proportion of land cover dominated by grass that maximized leporid numbers in the prey records, around 30%, was substantially higher than the observed mean value of grass land cover, 10%, that surrounded nesting sites (Fig. 2,  Supplemental Materials 1 (rapt-55-01-06_s01.pdf)). Our data set cannot tell us whether those higher levels of grass land cover contributed directly to higher jackrabbit populations, or whether the higher proportion of grass increased the abundance of leporids in prey for different reasons, such as increasing eagle foraging efficiency. Importantly, the NLCD grass category does not distinguish between exotic annual grasses, such as cheatgrass (which may also be present in the shrub land cover category), and native grasses (Boyte et al. 2016). Our analyses also do not address whether the spatial configuration of land cover patches near nesting sites was related to prey. But our results do suggest that surveys of prey populations may be useful to test hypotheses as to the relationship between leporid abundance and grass/sagebrush cover throughout Utah and the western USA, including extent, composition, and configuration.

Another striking pattern in the data was that occurrence of a primary prey species, black-tailed jackrabbits, was a notable predictor of both cottontails and rock squirrels, suggesting there may be synchronicity among prey population fluctuations. Jackrabbit populations can experience geographically broad multi-annual fluctuations (reviewed in Flinders and Chapman 2003, and Simes et al. 2015). Cottontail populations may also exhibit multi-annual cyclic fluctuations in abundance (Chapman and Litvaitis 2003). Cottontail population fluctuations were highly correlated with Greater Sage-Grouse (Centrocercus urophasianus) populations in Wyoming (Fedy and Doherty 2011), and it is possible that population fluctuations in our focal prey species were also correlated in Utah during the study period.

We note that visible prey remains were catalogued in the field as opposed to collecting all nest material for identification in a laboratory through comparison to museum samples. We suggest caution when interpreting such data as they do not represent all possible prey. However, this sampling bias was consistent among sites and years, as the same investigator identified and recorded all prey remains data. Similarly, prey remains gathered from raptor nesting sites are known to be a biased sample of prey in general (Collopy 1983, Redpath et al. 2001). Prey remains from Golden Eagles are likely particularly biased since this long-lived species will often retain a breeding territory but not attempt breeding in every year, presumably somewhat in response to perceived poor food resources in that year (Slater et al. 2017, Kochert et al. 2020). In this study, the number of nesting sites sampled each year varied due to the investigators' schedule as well as variability in egg-laying and territory occupancy ( Fig. S5 (rapt-55-01-06_s02.pdf)). Further information on territory occupancy and reproduction attempts by these eagles would have improved our ability to infer underlying trends in prey populations.

An understanding of prey-habitat relationships is particularly urgent given the extent and severity of threats to the Great Basin ecosystem, and the potentially negative effects of prey habitat loss or modification to resident and migratory Golden Eagles. Exotic annual grasses, such as cheatgrass, have also replaced or degraded shrublands and grasslands within large portions of the western USA (Pellant and Hall 1994, Reid et al. 2008, Brooks et al. 2016, Pyke et al. 2016). These grasses can provide highly combustible and continuous fuel for fires, which can degrade jackrabbit habitat and negatively impact Golden Eagles (Knick and Dyer 1997, Kochert et al. 1999, Heath and Kochert 2016). Long-term conservation of Golden Eagles in the western USA may hinge in large part on prevention of the loss of prey habitat to the exotic grass/fire cycle (US Fish and Wildlife Service 2016). Conservation of jackrabbits and other prey is largely concordant with existing strategies for preserving or restoring habitat for other wildlife species of concern, such as Greater Sage-Grouse and Brewer's Sparrows (Spizella breweri; Range-wide Interagency Sage-grouse Conservation Team 2012, Bureau of Land Management 2015). However, management prescriptions to benefit these species may well differ from those that are optimal for jackrabbits and other Golden Eagle prey. Thus, we encourage landscape-level efforts (to maintain or enhance habitat for the benefit of other sagebrush-steppe wildlife species) that include a component to monitor the population-level responses of leporids and other key prey species.

Information on spatial and temporal patterns in Golden Eagle diets, along with species-specific information on important prey, can be used to inform prey management with an emphasis on supporting the development of compensatory mitigation, conservation banking, and other resource management programs related to Golden Eagles. Despite the long duration of this study, we found little evidence that relative amounts of rain or unmeasured, long-term temporal trends strongly influenced the prey remains found at used Golden Eagle nesting sites. In contrast, land cover variables that are relatively straightforward to quantify were more strongly associated with the overall patterns of prey remains. Management actions that conserve, improve, or restore sagebrush-steppe habitat are likely to positively influence prey abundance and subsequently food availability, productivity, and overall fitness for Golden Eagles.

Supplemental Materials (available online).  Supplemental Materials 1 (rapt-55-01-06_s01.pdf). Investigation into land cover change at Golden Eagle nesting sites in Utah.  Supplemental Materials 2 (rapt-55-01-06_s02.pdf).  Table S1 (rapt-55-01-06_s02.pdf): Summary of Golden Eagle prey remains.  Figure S1 (rapt-55-01-06_s02.pdf): Predicted number of black-tailed jackrabbits within prey remains for Golden Eagle nesting sites.  Figure S2 (rapt-55-01-06_s02.pdf): Predicted number of cottontails (all species) within prey remains for Golden Eagle nesting sites.  Figure S3 (rapt-55-01-06_s02.pdf): Predicted number of rock squirrels within prey remains for Golden Eagle nesting sites.  Figure S4 (rapt-55-01-06_s02.pdf): Predicted number of yellow-bellied marmots within prey remains for Golden Eagle nesting sites.  Figure S5 (rapt-55-01-06_s02.pdf): Number of Golden Eagle nesting sites visited by year.