Despite regional variation, Gymnorhinus cyanocephalus (Pinyon Jay) densities generally increase with local pinyon–juniper cover and heterogeneous ground cover

Abstract

Traditionally, local-scale habitat relationship models are developed over small spatial extents, limiting model transferability and inference outside the study area. Thus, habitat managers frequently lack fine-scale information regarding the influence of vegetation composition and structure on site suitability or species abundance. Gymnorhinus cyanocephalus (Pinyon Jay) represents one declining species for which managers have limited information regarding the influence that vegetation composition and structure have on abundance at broad scales. To address this need, we developed a hierarchical Bayesian abundance model using summertime bird and vegetation data collected under the Integrated Monitoring in Bird Conservation Regions program to explain jay abundance as a function of local conditions. Our G. cyanocephalus abundance model allowed abundance relationships with pinyon pine (Pinus edulis and P. monophylla) and juniper (Juniperus spp.) to vary by ecoregion, thereby accounting for potential regional differences in habitat associations. We found G. cyanocephalus abundance was generally positively associated with pinyon pine and juniper cover; however, habitat relationships varied by ecoregion. Additionally, we found positive associations between jay abundance and grass cover, sagebrush cover, and percent bare ground. Our results agree with prior research suggesting mechanical removal of pinyon pine and juniper trees for sagebrush restoration or fuel treatments may negatively affect G. cyanocephalus. Managers wishing to reduce pinyon and juniper tree cover without negatively affecting G. cyanocephalus may benefit from targeting sites where both large-scale distribution models and our local habitat relationships suggest G. cyanocephalus are likely to occur in low numbers. Additionally, our modeled relationships indicate restoration that increases grass cover, sagebrush cover, and bare ground, while maintaining pinyon and (or) juniper cover, may lead to increased local densities of G. cyanocephalus.

RESUMEN

Tradicionalmente, los modelos de relación de hábitat a escala local se desarrollan en extensiones espaciales pequeñas, limitando la transferibilidad del modelo y la inferencia fuera del área de estudio. Por lo tanto, los gestores de hábitats frecuentemente carecen de información a escala fina sobre la influencia de la composición y estructura de la vegetación en la conveniencia del sitio o la abundancia de la especie. Gymnorhinus cyanocephalus representa una especie en declive para la cual los gestores tienen información limitada sobre la influencia que la composición y estructura de la vegetación tienen en la abundancia a gran escala. Para abordar esta necesidad, desarrollamos un modelo jerárquico bayesiano de abundancia utilizando datos de aves y de vegetación de verano recopilados bajo el programa de Monitoreo Integrado en Regiones de Conservación de Aves para explicar la abundancia de G. cyanocephalus en función de las condiciones locales. Nuestro modelo de abundancia de G. cyanocephalus permitió que las relaciones de abundancia con el pino piñonero (Pinus edulis y P. monophylla) y el enebro (Juniperus spp.) variaran según la ecorregión, teniendo en cuenta las posibles diferencias regionales en las asociaciones de hábitat. Encontramos que la abundancia de G. cyanocephalus generalmente se asociaba positivamente con la cobertura de pino piñonero y enebro; sin embargo, las relaciones de hábitat variaban según la ecorregión. Además, encontramos asociaciones positivas entre la abundancia de G. cyanocephalus y la cobertura de hierba, la cobertura de artemisa y el porcentaje de suelo desnudo. Nuestros resultados concuerdan con investigaciones previas que sugieren que la eliminación mecánica de árboles de pino piñonero y enebro para la restauración de artemisa o los tratamientos de combustible pueden afectar negativamente a G. cyanocephalus. Los gestores que deseen reducir la cobertura de árboles de pino piñonero y enebro sin afectar negativamente a G. cyanocephalus pueden beneficiarse de seleccionar sitios donde tanto los modelos de distribución a gran escala como nuestras relaciones de hábitat local sugieran que G. cyanocephalus se presente probablemente en números bajos. Además, nuestras relaciones modeladas indican que la restauración que aumente la cobertura de hierba, la cobertura de artemisa y el suelo desnudo, mientras se mantiene la cobertura de pino piñonero y (o) enebro, puede conducir a un aumento de las densidades locales de G. cyanocephalus.

INTRODUCTION

Biodiversity is declining globally due to over-exploitation and land-use change (Caro et al. 2022). As a result, conservationists frequently model distributions of declining wildlife to help identify regions for protection and (or) restoration (Peterson and Soberón 2012, Porfirio et al. 2014, Harms et al. 2017). However, misclassifications in remote-sensing data (Qin and Liu 2022) and limitations in the transferability of ecological relationships (Yates et al. 2018) may result in a lack of precision to inform fine-scale decisions aimed at population recovery across large spatial extents. Therefore, managers may benefit from incorporating fine-scale habitat descriptions into management decision-making to help recover populations locally, and across regional scales.

One declining species of growing conservation interest is the Gymnorhinus cyanocephalus (Pinyon Jay), a colonial and cooperative songbird that demonstrates a close mutualism with pinyon and juniper woodlands of southwestern North America. Gymnorhinus cyanocephalus populations have demonstrated both long- (Sauer et al. 2020) and short-term (Van Lanen et al. 2023) population declines. These range-wide population declines have resulted in the proposed listing of the G. cyanocephalus (Defenders of Wildlife 2022) under the Endangered Species Act (United States 1983) and increased concern regarding the conservation of this species (Somershoe et al. 2020). Declining G. cyanocephalus populations are puzzling, given pinyon pines (Pinus edulis and P. monophylla) and juniper (Juniperus spp.; hereafter, “pinyon–juniper”) cover is increasing at broad spatial extents across the western United States (Sankey and Germino 2008, Filippelli et al. 2020).

Researchers and conservationists have posited several hypotheses regarding potential drivers of G. cyanocephalus population declines. Some evidence suggests mechanical removal of pinyon and juniper woodlands to recover sagebrush ecosystems or reduce wildfire fuel may be reducing G. cyanocephalus habitat (Johnson and Sadoti 2019, Zeller et al. 2021, Van Lanen et al. 2023). Other studies have demonstrated reduced mast production of pinyon pines (Redmond et al. 2012, Wion et al. 2020) and local tree mortality has increased (McDowell et al., 2015), both due to warming and drying climates and drought. Climate-driven reductions in pinyon nut production have subsequently been suggested as a potential driver of G. cyanocephalus population declines (Boone et al. 2018). Other research suggests changes in climate or land use may be altering woodland structure (including tree age and density) and health, potentially negatively affecting G. cyanocephalus habitat (Johnson et al. 2017a, Boone et al. 2018). Still other research has proposed jay populations may be declining due to anthropogenic development (Kleist et al. 2018, Somershoe et al. 2020, Van Lanen et al. 2023) or declining insect populations (Collen et al. 2012; Dirzo et al., 2014). Of the above hypotheses, only mechanical removal of pinyon–juniper woodlands has been empirically tested and supported by experimental findings (Johnson et al. 2018, Magee et al. 2019). Thus, there remains little consensus or scientific support for the potential drivers associated with these long-term G. cyanocephalus population declines.

Regardless of the cause(s) of G. cyanocephalus declines, jay recovery may benefit from studies evaluating fine-scale habitat characteristics for the species throughout much of its range. Such information can inform site-specific decision-making, including guiding G. cyanocephalus habitat management and steering pinyon–juniper reduction projects away from areas expected to support large numbers of jays. Although several studies have previously evaluated G. cyanocephalus habitat characteristics, these studies either leverage remotely sensed data for habitat associations (Johnson et al. 2016, Zeller et al. 2021, Van Lanen et al. 2023) or develop fine-scale habitat associations at small geographic extents within the jay range (Johnson et al. 2013, Novak et al. 2021). In both cases, the consideration of scale in ecological modeling comes into play (Levin 1992), in which regional-scale efforts may fail to account for local differences in ecological processes and local-scale efforts may suffer from poor transferability. Furthermore, most habitat associations have been largely restricted to the southern portion of the species’ range (Boone et al. 2018, Johnson and Balda 2020) and prior work suggests G. cyanocephalus habitat relationships may vary considerably across space (Johnson and Sadoti 2019).

Issues of model transferability when evaluating G. cyanocephalus habitat associations are underscored by the immense variability in pinyon–juniper woodland characteristics. Communities broadly characterized as “pinyon–juniper woodland,” which support G. cyanocephalus, are extremely diverse in both structure and species composition. In total, 5 distinct species of juniper (Juniperus deppeana, J. monosperma, J. occidentalis, J. osteosperma, and J. scopulorum) and 2 species of pinyon pine (P. monophylla and P. edulis) occur within areas occupied by G. cyanocephalus (Johnson and Sadoti 2023) and stands may be comprised of single species or mixtures of species (Romme et al. 2009). As a result, 77 distinct vegetation classes of pinyon-dominated and another 71 classes of juniper-dominated communities have been recognized (Romme et al. 2009). Stand structure within pinyon–juniper communities also varies, with sites supporting persistent woodlands (characterized by variable densities of pinyon and or juniper trees, coarse soils, and infrequent disturbance), pinyon–juniper savannas (characterized by mixtures of trees and grasses, likely shaped by low-intensity fires), and wooded shrublands (characterized as shrublands with intermittent pinyon and [or] juniper trees; Romme et al. 2009). The tremendous diversity in pinyon–juniper woodland characteristics is also coupled with site-specific variation including soil conditions, understory structure, and associated long and short-term weather. Collectively, local-scale variability is almost certain to result in significant variability in resources available to G. cyanocephalus, including pinyon pine mast, juniper berry production, arthropod communities, and nesting cover.

The lack of fine-scale G. cyanocephalus habitat associations across larger spatial extents which accounts for variability in pinyon–juniper woodlands therefore represents a knowledge gap, potentially inhibiting managers’ ability to implement appropriate management actions to stabilize or reverse G. cyanocephalus population declines. Such information can guide recovery efforts, identify potential refugia for G. cyanocephalus, help managers engaging in conifer removal practices avoid areas with characteristics associated with high G. cyanocephalus densities, and provide insight into environmental impact assessments. To address this information gap and help inform G. cyanocephalus recovery efforts, we developed a hierarchical abundance model identifying fine-scale habitat characteristics associated with high summertime G. cyanocephalus densities. We then evaluated (1) the strength and direction of relationships between G. cyanocephalus abundance and vegetation structure and composition; and (2) the extent to which relationships between pinyon and juniper cover and jay abundance vary across portions of the G. cyanocephalus range. We expected G. cyanocephalus abundance would be positively associated with herbaceous, juniper, pinyon, sagebrush, and woody cover and negatively associated with bare ground and grass cover. We expected intermediate values of mean canopy height and elevation would be associated with high G. cyanocephalus abundance. Our abundance–habitat relationship model, which allows for regionally varying relationships across the species’ range, has broad applications for other species and systems and addresses previous researcher concerns (Mannocci et al. 2020). In particular, our model formulation may be appropriate when species occur in variable habitats or across large geographic extents.

METHODS

Study Area

We modeled G. cyanocephalus abundance across 13 level III ecoregions (U.S. Environmental Protection Agency 2011), spanning portions of 16 states (Arizona, California, Colorado, Idaho, Kansas, Montana, North Dakota, Nebraska, New Mexico, Nevada, Oklahoma, Oregon, South Dakota, Texas, Utah, and Wyoming) within the western United States (Figure 1). Ecoregions within our study area encompassed grassland, sagebrush, riparian, and forested ecosystems, with elevations ranging from –85 to 4,397 m.

Study Data

Gymnorhinus cyanocephalus count data

We modeled fine-scale G. cyanocephalus abundance using data associated with summertime point-count surveys collected under the Integrated Monitoring in Bird Conservation Regions (IMBCR) program from 2008 to 2020 (Pavlacky et al. 2017). Raw data were provided by the Bird Conservancy of the Rockies and may be made available upon reasonable request (see Data Availability section). Sample units under the IMBCR program are represented by 1-km² grid cells, each containing 16 point-count stations. IMBCR point-count stations are arranged uniformly throughout the grid cell and separated by 250 m. Grid cells to be sampled are selected using a spatially balanced random algorithm (Stevens and Olsen 2004, Pavlacky et al. 2017). Seasonal timing of surveys was dependent upon elevation, latitude, survey logistics, expert local knowledge, and the ecosystem represented within the grid cell (Hanni et al. 2012), with all surveys occurring between April 26 and July 31 each year.

Trained observers conducted avian point-count surveys in the early morning (Buckland 2006, Hanni et al. 2012), with counts occurring between 30 min prior to local sunrise and finishing no later than 5 hr after sunrise (mean survey duration was 3.38 hr). During point-count surveys, observers recorded all birds detected, the radial distance to each detected individual, mode of detection, and the minute interval within which the detection occurred. Observers were instructed to avoid double-counting of individuals within a survey (Hanni et al. 2012). We removed all G. cyanocephalus records demarcated in the raw data as potential migrants and/or individuals flying over, but not using, resources within the grid (Hanni et al. 2012). Point counts lasted 5 min in 2008 and 2009 and 6 min in all subsequent years.

We originally wished to model G. cyanocephalus abundance at the subsample level of the point-count station. However, the large data set and model complexity resulted in prohibitively long model runs and large computational memory requirements. We therefore summed the number of independent G. cyanocephalus detections for all sampled point-count stations within a grid each year. For analyses, each independent detection was assigned to a point-count minute interval and one of 10 equal-sized distance bins to facilitate estimation of availability and detection (see Model Structure and Modeling Procedures section).

We also excluded all data within ecoregions for which there were fewer than 20 detections during the 13 years of our study to facilitate model fitting.

Site visit data

To account for diurnal rhythms in G. cyanocephalus activity, we calculated the mean start time of all point counts conducted within a grid each year. We then subtracted the mean point-count start time from the local sunrise time on the date of the survey, derived using the sunrise function in the maptools package (Bivand and Lewin-Koh 2020) in R (R Core Team 2020). To model behavioral changes in availability related to G. cyanocephalus detection (likelihood a jay was vocalizing or perched in a visible location) throughout the summer sampling season, we calculated the ordinal date each grid cell was sampled during each year of the study. Finally, we accounted for variation in observer ability by deriving a binary variable, “observer experience,” indicating if the point-count surveyor had conducted IMBCR surveys in a prior breeding season.

Vegetation data

IMBCR surveyors recorded ocular estimates of vegetation characteristics within a 50-m radius surrounding each point-count station, immediately upon arrival at the station. The surveyors collected data summarizing the ground cover, understory, and overstory vegetation while remaining stationary. Thus, this period of data collection also served as “cool down” time for birds to revert to behaviors they were engaged in upon the surveyors’ approach. Surveyors recorded woody vegetation greater than 3 m in height within the “overstory” layer; woody vegetation between 0.25 and 3 m in height within the “understory” layer; and estimated % cover amounts of bare or litter covered ground, grass, herbaceous forbs, snow, water, and short woody vegetation <0.25 m in height within the “ground cover” layer (see details below). Beginning in 2012, surveyors estimated residual (from the previous growing season) and live (from the current growing season) grass cover percentages separately. Prior to the 2012 sampling season, observers estimated a single cover value for both live and residual grass. Beginning in the 2018 breeding season, IMBCR surveyors also separated litter cover from bare ground. To maintain consistency across the data set, we summed these split ground cover categories (bare ground and litter cover became “bare/litter”; live and residual grass cover became “grass”) to make all data equivalent throughout the study duration.

To speed ocular estimation, surveyors assigned cover amounts of 0, 1, 5, or multiples of 10% to each cover category. Within both the overstory and understory layers, surveyors recorded the 5 most abundant vegetation species and the associated proportion of that cover each species represented, with summed cover values across these classes equaling 100%. Lastly, surveyors estimated the mean height of the entire overstory and understory layers to the nearest 1 and 0.25 m, respectively. IMBCR surveyors estimated vegetation cover proportions of the understory layer for regions not occupied by vegetation classified in the overstory layer. To account for this, we transformed the observer-collected values to measures of absolute vegetation cover (Supplementary Material).

Pinyon juniper trees and shrubs in our study area frequently reached heights close to the 3-m IMBCR overstory and understory cutoff values. To address this, we calculated a single total percent cover value for these cover types representing cover within both the overstory and understory layers. We did this by multiplying the overstory and understory cover values by the species composition pertaining to these species and summing these values. We calculated mean values within a survey for each vegetation covariate, consistent with the G. cyanocephalus count and site visit data.

Model Structure and Modeling Procedures

We developed a hierarchical abundance model, incorporating distance-sampling (Buckland et al. 2001) and time removal procedures (Amundson et al. 2014), to estimate G. cyanocephalus abundance within grid cells (density) as a function of fine-scale vegetation variables closely adhering to procedures described by Van Lanen et al. (2023). Following similar work, we modeled site suitability for jays using a zero-inflation term informed by a generalized additive model (Van Lanen et al. 2023). Due to concerns regarding model transferability across space for G. cyanocephalus (Johnson and Sadoti 2019), we allowed the intercept of our abundance model to vary by ecoregion. Similarly, we allowed the influence of both pinyon pine cover and juniper overstory cover on abundance to vary by ecoregion to better estimate regional differences in habitat associations with these 2 tree types. To account for potential variance-covariance among these intercept and slope parameters, we included an inverse Wishart distribution prior (Van Lanen et al. 2023). We also incorporated linear effects of bare ground, overstory cover, herbaceous ground cover, grass ground cover, sagebrush cover, understory cover, understory mean height, and woody ground cover on abundance. We suspected G. cyanocephalus may occur at higher densities when elevation and canopy heights were intermediate and therefore incorporated quadratic forms for the influence of these variables (Table 1). To control for annual variation, we included a random year effect on G. cyanocephalus abundance. Finally, we included a survey-specific error term within our G. cyanocephalus abundance model. When plotting associations between abundance and habitat covariates we calculated jay density using the predicted abundance values and maximum detection distance following truncation (Kéry and Royle 2016).

We used a half-normal distance function to estimate detectability, including a scale parameter for the rate at which detection probability declined across distance bins and which we allowed to vary with our derived “observer experience” variable. We modeled G. cyanocephalus availability using a removal model, informed by the minute intervals during points counts (Amundson et al. 2014), and we allowed availability to vary as a function of ordinal date and mean minutes since sunrise when point counts were conducted within a survey. We included linear and quadratic terms for ordinal date to account for potential nonlinearity in jay availability during the summer sampling season resulting from changes in jay behavior (e.g., incubating, feeding nestlings, and feeding fledglings).

We assessed potential correlations among our covariates using Pearson’s pairwise correlation values. When pairs of variables were correlated (r ≥ 0.6), we excluded one of the variables from the abundance model. After removing highly correlated covariates, we assessed multicollinearity among the remaining variables with variance inflation factors (VIF). We calculated VIF from the vif function in the car package (Fox and Weisberg 2019), using a formula in which the number of detected individuals during a survey was a function of all included covariate values, to ensure all values were ≤3.0 (Zuur et al. 2010).

We modeled fine-scale jay abundance in a Bayesian framework using NIMBLE (de Valpine et al. 2017) version 0.13.1 (NIMBLE Development Team 2022) in Program R (version 4.1.1; R Development Core Team 2020). Following convention, we removed the 10% furthest observations from the data set (Buckland et al. 2001) and centered and scaled all covariate values prior to analysis to facilitate MCMC chain convergence and improve model fit. We imputed missing covariate data within our data set by specifying a normal probability distribution for each variable (Royle 2009) and priors for the mean (N [0, 10,000]) and standard deviation (Unif [0, 3]). In total, 3,580 of a possible 208,260 covariate-grid-year combinations contained missing data, representing 1.72% of covariate values. We specified priors for parameters associated with detection, availability, and abundance as N(0, 10,000). We further specified priors for abundance intercepts and juniper and pinyon juniper effects associated with the inverse Wishart distributions as:

with 4 degrees for freedom.

We fit the model with four parallel chains using the parallel (R Development Core Team 2020), dclone (Solymos 2010), and snow (Tierney et al. 2018) packages. Model code and associated data are available via the accompanying U.S. Geological Survey (USGS) data release (Van Lanen et al. 2024). We ran the model for 1,000,000 iterations, discarding the first 200,000 iterations and thinning the remaining samples by 100, resulting in 32,000 retained samples for assessing chain convergence, model fit, and biological inference. To evaluate chain convergence, we visually inspected traceplots with the MCMCvis package (Youngflesh 2018) and ensured Gelman–Rubin potential scale reduction factors were <1.1 (R-hat; Gelman and Rubin 1992). We assessed model fit by comparing observed and predicted counts using a chi-square discrepancy posterior predictive check (i.e., Bayesian P-value; Gelman et al. 1996, Kéry and Royle 2016).

RESULTS

Our dataset contained avian and vegetation data from 4,720 independent IMBCR grid cells, 13,884 IMBCR grid cell visits (Supplementary Material Table S1; combination of grid cell and year), resulting in a total of 162,467 individual point-count surveys across ecoregions (Supplementary Material Table S2). Our dataset of G. cyanocephalus detections contained 4,467 independent observations within 12 distinct ecoregions after truncating the furthest 10% of detections (Table 2). Following the truncation of the data, the maximum radial distance of remaining detections was 370 m, resulting in 10 distance bins 37 m in width. Pearson’s pairwise correlation assessment identified 2 sets of highly correlated covariates: mean overstory height and percent overstory cover (r = 0.67), and percent understory cover and percent sagebrush cover (r = 0.62). We removed percent overstory cover from our final model (Table 1) because we assumed overstory cover was at least partially addressed via pinyon pine and juniper cover amounts. Similarly, we retained sagebrush cover instead of understory cover (Table 1), because sagebrush cover has previously been shown to influence G. cyanocephalus densities (Van Lanen et al. 2023). After removing these variables, we calculated a maximum VIF of 2.53, indicating acceptable levels of multicollinearity within our dataset. Gelman–Rubin potential scale reduction factors were between 1.00 and 1.06 for all parameters of interest. The posterior predictive check (Bayesian P-value) was 0.30, indicating satisfactory model fit.

Gymnorhinus cyanocephalus availability was largely stable throughout the summer sampling period (Table 3; Figure 2) and the probability of detecting jays did not vary with mean minutes since sunrise when point counts were conducted. Observers who had conducted IMBCR surveys in a previous breeding season detected G. cyanocephalus at slightly higher rates within a given distance bin compared to inexperienced observers (Table 3; Figure 2).

Within regions where they occur, G. cyanocephalus abundance was positively associated with the percentage of bare ground, grass, and woody vegetation within the ground cover layer. Grass cover had the strongest influence on G. cyanocephalus abundance (Table 3; Figure 3). The percentage of absolute sagebrush cover in the understory layer was positively and moderately associated with increasing G. cyanocephalus densities. Mean heights of both understory and overstory vegetation had little influence on G. cyanocephalus densities (Table 3; Figure 3).

Habitat associations within ecoregions supporting the different species of juniper occurring within our study area were either positive or neutral (Table 3; Figure 4). We found 4 positive relationships (Arizona/New Mexico Mountains, Central Basin and Range, Colorado Plateaus, and Wasatch and Uinta Mountains) and 3 neutral relationships (Arizona/New Mexico Plateau, Mojave Basin and Range, and Northern Basin and Range) among ecoregions supporting Utah juniper (J. osteosperma); 3 positive (Middle Rockies, Wasatch and Uinta Mountains, and Wyoming Basin) and 3 neutral (Arizona/New Mexico Plateau, Northwestern Great Plains, and Southern Rockies) relationships in ecoregions supporting Rocky Mountain juniper (J. scopulorum); 1 neutral association in the Northern Basin and Range ecoregion supporting Western juniper (J. occidentalis); and 1 positive (Southwestern Tablelands) and 2 neutral (Arizona/New Mexico Plateau and Southern Rockies) relationships in ecoregions supporting oneseed juniper (J. monosperma; Shinneman et al. 2023).

Gymnorhinus cyanocephalus abundance associations with juniper and pinyon pine were relatively consistent across the 12 ecoregions we investigated. Overall, G. cyanocephalus abundance was significantly (95% credible intervals did not overlap zero) and positively associated with juniper cover for 7 of the 12 ecoregions investigated (Table 3; Figure 4). Similarly, G. cyanocephalus densities were significantly and positively associated with pinyon pine cover within 5 ecoregions (Arizona/New Mexico Plateau, Arizona/New Mexico Mountains, Colorado Plateaus, Southern Rockies, and Wasatch and Uinta Mountains; Table 3; Figure 4), all of which overlapped the range of P. edulis. We found no evidence of significant negative influences of pinyon pine or juniper cover on G. cyanocephalus abundance within any of the investigated ecoregions (Table 3; Figure 4).

DISCUSSION

This study provides quantitative descriptions of fine-scale summertime (primarily, postbreeding) habitat characteristics associated with G. cyanocephalus densities throughout most of its range. Our findings validate previous concerns regarding model transferability (Johnson and Sadoti 2019), as we found relationships regarding the influence of pinyon pine and juniper tree cover on G. cyanocephalus densities varied throughout our study area, though were generally neutral or positive. The positive associations we found with contrasting ground cover conditions (e.g., bare ground, grass, and woody cover) suggest G. cyanocephalus may occur at higher densities in regions with heterogeneous ground cover conditions. The habitat associations we developed can be used by managers to (1) target potential fuels and conifer removal treatments within local areas which are unlikely to support high jay densities, and (2) inform restoration efforts designed to bolster G. cyanocephalus populations.

Habitat Relationships

The scale at which we evaluated habitat relationships (1 km²) most closely approximates the “colony scale” as detailed in prior work (Johnson and Sadoti 2023); however, we recognize a seasonal mismatch between our study, which used summertime detections, and prior studies characterizing colony habitat characteristics during the G. cyanocephalus breeding season (in spring and early summer). Nevertheless, the general pattern of positive associations with pinyon pine and/or juniper tree cover agrees with prior work (Johnson and Sadoti 2023), including local-scale studies within the southern and western portions of the G. cyanocephalus US range (Johnson et al. 2016, 2017b, Boone et al. 2021, Novak et al. 2021) and a landscape-scale study throughout much of the range of G. cyanocephalus (Van Lanen et al. 2023).

We found evidence that associations between jay abundance and pinyon pine and juniper cover varied throughout our study area, as has been suggested by other researchers (Johnson and Sadoti 2019) and which is apparent when viewing studies from disparate regions (Johnson and Sadoti 2023). Interestingly, we found G. cyanocephalus densities were positively associated with both pinyon and juniper cover in the south-central ecoregions (Wasatch/Uinta [19], Colorado Plateau [20], and Arizona/New Mexico Mountains [21]) of our study area (Figure 1), including central and eastern Utah and central Arizona. These ecoregions represent areas where oneseed juniper (J. monosperma), Utah juniper (J. osteosperma), and 2-needle pinyon (P. edulis) overlap (Shinneman et al. 2023). In other ecoregions we investigated, G. cyanocephalus densities were positively associated with either juniper or pinyon pine cover, but not both. We did not find any clear distinction among habitat associations within ecoregions supporting various juniper species in our study area. Jay densities were generally positively associated with pinyon pine cover within ecoregions overlapping the distribution of P. edulis (5 positive relationships [ecoregions 19, 20, 21, 22, and 23] and 2 neutral relationships [ecoregions 14 and 26]), but found no positive associations between G. cyanocephalus abundance and pinyon pine cover within the range of single-leaf pinyon (3 neutral relationships [ecoregions 13, 14, and 80]; Shinneman et al. 2023). Prior research found P. edulis seeds have higher fat and protein content compared to those produced by P. monophylla (Sutton 1984). Although it is possible the higher fat and protein content of P. edulis mast drive the positive association that we observed between G. cyanocephalus and pinyon cover, direct evaluations of relationships between G. cyanocephalus abundance and demography and pinyon mast quantity, nutritional quality, and processing times, remain areas for further investigation.

We found neutral associations between G. cyanocephalus abundance and pinyon pine or juniper cover within 3 ecoregions representing the boundaries of our study area (Mojave Basin and Range [14] in the southwest, Northwestern Great Plains [43] in the northeast, and Northern Basin and Range [80] in the northwest of our study area). These 3 ecoregions were among the 4 ecoregions with the fewest number of G. cyanocephalus detections in our dataset. We therefore acknowledge a lack of statistical power to detect a strong association between jay abundance and pinyon and juniper cover along the boundary of their range (i.e., Northwestern Great Plains) and (or) due to reduced IMBCR sampling effort (i.e., Mojave Basin and Range and Northern Basin and Range) in these areas. Additionally, the Northwestern Great Plains and Northern Basin and Range are largely devoid of pinyon pine, as reflected by the rug plots in Figure 4. There are reports of G. cyanocephalus using ponderosa (P. ponderosa) and Jeffrey (P. jeffreyi) pine (Johnson and Balda 2020), a habitat variable which we did not specifically evaluate. Thus, G. cyanocephalus may be more closely associated with ponderosa and Jeffrey pine cover within the Northwestern Great Plains and Northern Basin Range regions, where pinyon pine is absent. We suggest assessment of habitat associations near the boundary of the G. cyanocephalus range, particularly in the western portion (supporting single-leaf pinyon), northwestern region (supporting ponderosa and Jeffrey pine, without pinyon pine), and northeastern portion (supporting ponderosa pine, without pinyon pine), may represent an area for additional research. Additionally, we caution managers in the use of these modeled relationships within ecoregions lacking ample samples across the range of covariate values (Figure 4).

We found several positive associations with resource conditions which, based on previous studies, likely provide food or caching sites for G. cyanocephalus. We found a positive association between G. cyanocephalus densities and bare ground, which is seldom investigated explicitly in scientific studies. Our findings do agree, however, with anecdotal observations by Marzluff and Balda (1992) and Johnson and Balda (2020), who detailed G. cyanocephalus use of bare areas for caching pinyon seeds. Their observations, coupled with our findings, suggest explicit evaluation of bare ground in future G. cyanocephalus habitat relationship modeling may improve model predictions. Similarly, we found a strong positive association with grass cover, which also has not been specifically investigated in most prior G. cyanocephalus home-range, colony, or nest site selection habitat association work (Johnson and Sadoti 2023). However, researchers also have documented increased use of mixed lower-elevation grasslands and pinyon and/or juniper woodlands during the nonbreeding season (Johnson et al. 2016) and while caching food (Boone et al. 2021). Jays are also known to feed on grasshoppers (Johnson and Balda 2020), which are likely more abundant in areas with higher grass cover. Thus, cover for caching seeds and invertebrates may drive the positive association we observed between grass cover and G. cyanocephalus abundance. Finally, the positive association between G. cyanocephalus abundance and woody ground cover we found also aligns with Boone et al. (2021) who found G. cyanocephalus cached seeds in areas with higher woody ground cover compared to adjacent sites. We recognize managing for increased bare ground is likely to result in reductions in grass and woody ground cover (and vice versa) and thus represents a management conundrum. Our results suggest G. cyanocephalus have heterogeneous resource needs within their home range to provide conditions suitable for their full suite of daily behaviors (e.g., caching and foraging) and throughout their annual cycle (e.g., nesting and brood-rearing). Thus, management designed to produce a mosaic containing patches of bare ground, grass, and woody ground cover, within or surrounded by pinyon and juniper canopy cover, may best support this species throughout its daily and annual activity patterns; however, we were unable to evaluate this explicitly.

The positive association we observed between G. cyanocephalus densities and sagebrush cover adheres to modeled landscape-scale habitat associations (Van Lanen et al. 2023), evidence from nest site selection investigations (Novak et al. 2021), and associations developed from a foraging and caching study (Boone et al. 2021). Together, these findings suggest G. cyanocephalus may frequently occur within wooded shrublands, particularly outside of the breeding season (Johnson and Sadoti 2023). Wooded shrublands, in these systems, are often found at pinyon–juniper and sagebrush ecotones, which represent sites frequently targeted for conifer removal projects to restore sagebrush within the northern portion of the G. cyanocephalus range (Boone et al. 2018). These conifer removal projects may cause undesired trade-offs for G. cyanocephalus resulting from mechanical conifer removal (Boone et al. 2018, Johnson et al. 2018, Johnson and Sadoti 2019, Van Lanen et al. 2023). Findings from 2 experimental studies evaluating occupancy rates and site use of G. cyanocephalus before and after conifer removal further support the negative consequences of conifer removal treatments for G. cyanocephalus (Johnson et al. 2018, Magee et al. 2019).

Study Limitations

Although we suggest sites capable of supporting high summertime G. cyanocephalus densities may help indicate important areas for G. cyanocephalus conservation, we recognize abundance may be a misleading indicator of habitat quality for many species and systems (Van Horne 1983) and jay habitat use may vary seasonally. Habitat relationships based upon demographic rates conducted throughout the full annual cycle therefore may provide better inference regarding high-quality habitat characteristics. Additionally, observational studies represent a weak form of inference regarding wildlife response to management, and before-after-control-impact experiments, beyond the single study conducted within the G. cyanocephalus range (Magee et al. 2019), may provide stronger inference regarding G. cyanocephalus response to conifer removal efforts.

We leveraged ocular estimates of vegetation characteristics surrounding point-count stations in our habitat relationship modeling. Although doing so provided on-the-ground information, avoiding misclassification and biases sometimes associated with remotely sensed data (Li et al. 2021), ocular vegetation estimates may suffer from biases, imprecision, and interobserver variability (Block et al. 1987). Furthermore, our use of vegetation information within a 50-m radius surrounding the point-count station as an index of vegetation characteristics within areas used by G. cyanocephalus, despite including G. cyanocephalus detections nearly 370 m from the point-count location in our data set. Thus, there was a mismatch in the spatial scales between the bird and vegetation data we used in our modeling effort. Furthermore, we did not consider regional landscape characteristics which undoubtedly influence G. cyanocephalus densities and habitat use. We also note the vegetation data we used did not differentiate between species of pinyon pine or juniper. Species-specific differences within pinyon and juniper trees may convey differential benefits to G. cyanocephalus, which we were unable to explicitly investigate in areas where distributions of multiple species overlap. Finally, the IMBCR program does not sample across the entire occupied range of the G. cyanocephalus and sampling within some areas where G. cyanocephalus occur (Johnson and Balda 2020), including Arizona, Nevada, New Mexico, and Oregon was relatively sparse (Figure 1). Therefore, caution should be considered when applying our modeled habitat relationships within regions with sparse IMBCR sampling (Figure 1; Figure 4).

Management Applications

Our findings provide information relevant to the management and recovery of G. cyanocephalus in 2 primary ways. First, the habitat associations can be used to inform habitat restoration efforts designed to bolster G. cyanocephalus populations. Our results indicate maintaining a sagebrush understory, grass cover, and patches of bare ground, in addition to pinyon and juniper cover, may lead to sites supporting high G. cyanocephalus densities. Therefore, management practices that include seeding and planting of sagebrush and grasses, as well as grazing and invasive forb control to increase patches of bare ground, may benefit G. cyanocephalus populations. Secondly, the habitat associations we developed here can be used in conjunction with regional habitat suitability (Zeller et al. 2021) or density distribution maps (Van Lanen et al. 2022, 2023) to assess a site’s potential to support G. cyanocephalus. Such an exercise could be used to target areas for future G. cyanocephalus research, including studies that may benefit from efficiently locating colonies and/or individuals. Similarly, managers seeking to reduce pinyon and juniper tree cover, such as efforts to restore sagebrush habitat for sage-grouse and/or to manage fuel loads (Boone et al. 2018), can assess local site conditions to determine whether G. cyanocephalus are likely to occur at high densities. Specifically, habitat managers and fuel crews can use Figures 3 and 4 to estimate an index of local site suitability, based upon current vegetation cover, and evaluate the relative influence of vegetation manipulations on the G. cyanocephalus. In such applications, current expected G. cyanocephalus index values could be used as target thresholds for evaluated sites. Deriving G. cyanocephalus population indexes may also be particularly useful for evaluating multiple proposed project areas to determine where potential negative effects can be minimized, or potential positive effects can be maximized, with respect to G. cyanocephalus. Our results support a growing body of empirical (Johnson et al. 2018, Magee et al. 2019) and model-derived evidence (Zeller et al. 2021, Van Lanen et al. 2023), that suggests management which reduces pinyon–juniper cover is likely to negatively affect G. cyanocephalus habitat. Therefore, an assessment using both regional maps of G. cyanocephalus density or habitat use, coupled with our fine-scale vegetation associations, can inform environmental impact assessments, and ultimately help managers determine if pinyon–juniper removal at a particular site is consistent with management and conservation objectives.

Finally, we suggest our modeling approach serves as a valuable framework for other species and systems in which abundance–habitat relationships are liable to vary across a studied geographic footprint and/or when species are associated with variable habitats. By allowing abundance–habitat relationships to vary with ecoregion we were partially able to address the persistent issues of scale and transferability in ecology. As a result, our model estimated individual abundance–habitat relationships for relatively small ecological regions while producing model results across a large geography. We provide our code via the associated data release (see Data Availability section) to assist any similar future efforts.

Acknowledgments

We are grateful to the Bird Conservancy of the Rockies and the IMBCR partnership for sharing their bird and vegetation data with us. Our model was developed with assistance from David Koons and Adam Green. Austin Nash and Giancarlo Sadoti kindly provided suggested edits, greatly enhancing this manuscript. Our manuscript was further improved by helpful suggestions from associate editor, Desirée Narango. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government.

Funding statement

We appreciate the Bureau of Land Management and the U.S. Geological Survey for funding this work.

Ethics statement

This study used data collected previously by a number of other organizations. Data were passively collected and no animals were handled in the gathering of these data. Data in analyses were used with written permission from the data owner.

Conflict of interest statement

The authors declare that they have no conflicts of interest.

Author contributions

N.V.L. and C.A. conceived the idea, design, experiment. N.V.L. performed the experiments. N.V.L. wrote the paper. N.V.L., A.M., and C. A. developed and designed methods. N.V.L. and A.M. analyzed the data. C.A. contributed substantial materials, resources, or funding.

Data availability

Data and R code used to model density–habitat relationships in this study are available via a USGS data release (Van Lanen et al. 2024). The raw data underlying this article cannot be shared publicly to protect the rights of the data owners and the anonymity of private landowners who participated in the IMBCR program. Researchers interested in obtaining raw IMBCR data should contact Bird Conservancy of the Rockies and submit a data-sharing request through the online portal at https://forms.monday.com/forms/80217ca03b7de4ed366dce0a608172d8?r=use1.

LITERATURE CITED