Temperate grassland ecosystems are one of the most threatened ecosystems worldwide, and their loss endangers the grassland songbirds that rely upon them. This guild of birds has shown long-term declines in North America. At the same time, American bison (Bison bison) are becoming more common through reintroductions, and they may make significant modifications to grassland songbird habitat. To support conservation for this guild, we sought to understand the importance of bison grazing and ecosystem productivity to the species richness, occupancy, and abundance of this avian community. We conducted dependent double-observer bird counts, measured bison grazing intensity with patty counts, and used remote-sensed Normalized Difference Vegetation Index (NDVI) data to measure ecosystem productivity. Our work took place in the National Bison Range near Moiese, Montana and in Yellowstone National Park in Wyoming. We found that species richness was positively correlated with patty counts, and had a weak negative correlation with NDVI. Occupancy probability for six of seven grassland songbird species was positively correlated with patty counts, and for six of seven species was negatively correlated with NDVI. Abundance of vesper sparrow (Pooecetes graminueus) and western meadowlark (Sturnella neglecta) were positively correlated with patty counts, although for western meadowlark, this trend became less positive with increasing patty counts. Our work suggests that managers may want to encourage a broad range of bison grazing intensities to ensure that vegetative conditions related to bison grazing are present for all species.

Study sites

We conducted this study in two intermountain valleys of the North American Rocky Mountain Range -  National Bison Range (NBR), MT, USA and Yellowstone National Park (YNP), WY, USA (Fig. 1). Bison have been present in both locations for over 100 years, and therefore we expect the longterm effects of grazing, and any potential effects on songbird communities, are well established. Differences between the study sites in bison management, landscape size, and vegetation composition provide a useful representation of conditions that grassland songbirds may encounter in other intermountain valleys of the area. The NBR was managed by the U.S. Fish and Wildlife Service and maintained approximately ~350 bison at the time of this study. NBR had approximately eight individual pastures in which bison were allowed to move freely amongst pastures (e.g., open gates, removal of some fencing) within the 76 km² NBR  boundary. YNP is managed by the U.S. National Park Service. The YNP study site took place in 890 km² area of the Northern Range where approximately 3,500-4,000 bison are present (Mosley and Mundinger 2018). Bison movements are unmanaged within YNP, and in contrast with NBR, bison can move out of the YNP study area. While both study areas are intermountain valleys, NBR vegetation is generally described as Palouse Prairie and YNP as sagebrush-steppe.

Study species

We included grassland songbird species associated with grassland and sagebrush-steppe ecosystems for which we made at least 10 observations. The grassland species were clay-colored sparrow (Spizella pallida), grasshopper sparrow (Ammodramus savannarum), vesper sparrow (Pooecetes gramineus), and western meadowlark (Sturnella neglecta). The sagebrush-steppe species were Brewer’s sparrow (Spizella breweri), green-tailed towhee (Pipilo chlorurus), and sage thrasher (Oreoscoptes montanus). All species were included in estimating species richness and occupancy. We estimated abundance for vesper sparrow and western meadowlark, species observed on both NBR and YNP, to further explore how bison grazing intensity influences individual species.

Sampling frame

                We selected survey plots of 250 x 250 m² and stratified plot selection across three levels of expected bison grazing using a habitat suitability index (HSI) for bison (Steenweg et al. 2016) that incorporated slope (USGS 2015), distance to water (USGS 2013), and vegetation type (USGS 2012). HSI values range from 0 to 1, with 0 representing the least suitable habitat, and 1 representing the most suitable habitat for bison. We binned the HSI values from both study sites into three strata of expected bison grazing: Low (0.17–0.46), Medium (0.46–0.61), and High (0.61–0.80). We randomly selected survey plots in each stratum that contained ≥ 75% grassland or shrub steppe vegetation. On NBR, we sampled 10 plots in the Low stratum, 30 plots in the Medium stratum, and 30 plots in the High stratum each year. Compared to Medium and High stratum (43% and 32% of NBR, respectively), there was less Low stratum available for sampling on the NBR (approximately 11% of NBR), because much of this area was also forested. To compensate, we sampled a larger number of Low stratum plots at the YNP site. In YNP, we sampled 30 plots in the Low stratum, 13 plots in the Medium stratum, and 12 in the High stratum in 2016. In 2017, we sampled 25 plots in the Low stratum, 8 in the Medium stratum, and 11 in the High stratum.

Avian surveys

 Observers walked a transect through the middle of each sample plot (Fig. 2) surveying birds up to 125 m to each side. Beyond 125 m, human detection of birds declines dramatically (Ralph et al. 1995). We used a dependent double-observer method described by Nichols et al. (2000), that has been employed in similar types of land cover (e.g., Tipton et al. 2008, 2009; Golding et al. 2017). To execute the dependent double observer method, two observers walk transects through the middle of a plot. The primary observer walks first and communicates all adult birds they see within the plot to the secondary observer. The secondary observer walks 3-5 meters behind the primary observer on the transect, records the birds the primary observer detects, and records any additional birds the primary observer fails to detect.Visual confirmation was required, thus auditory-only detections were not recorded. Observers used rangefinders to confirm each bird was within plot boundaries. Observers alternated who was the primary and secondary observer consecutively. This field protocol was established to account for imperfect detection when estimating abundance and we converted abundance data into occupancy data to estimate occupancy. Avian surveys were conducted in the early morning hours, from dawn until 1000 Mountain Daylight Time except in weather conditions (e.g., rain or winds ≥24 km/hr) that impaired detection of birds. In NBR, bird surveys were conducted from May 20 to July 7 in 2016 and from May 19 to June 30 in 2017. In YNP, bird surveys were conducted from June 1 to July 1 in 2016 and from May 31 to July 9 in 2017. Each plot was sampled twice within a year.

Explanatory Variables

Bison grazing intensity data We measured bison grazing intensity by counting bison patties in each plot (Fig. 2). Density of patties estimates bison grazing intensity at patch-level spatial scales and reflects vegetative responses to grazing intensity (Tastad 2013). Following Sliwinski (2011), patties were either individual, well-formed piles, or several closely associated piles. Observers counted patties within 1 m to each side of four transect lines per plot. Observers walked two transects across the entire plot after each avian survey. The transects ran in an East-West direction, and a North-South direction (Fig. 2).

Productivity data We measured productivity using the mean of the cumulative Normalized Difference Vegetation Index (NDVI) summarized at the plot level from April 1 to June 30 of each year. We chose this timeframe to capture the beginning of the vegetation growing season and the end of the bird breeding season. NDVI is a remote-sense measurement of photosynthetic pigments in vegetation and is considered a proxy for primary production of vegetation (Tucker and Sellers 1986, Paruelo et al. 1997). NDVI values range from -1 to 1, with 1 indicating the highest measurement of vegetative productivity. We used an established NDVI product that uses gap-filling and smoothing techniques with a resolution of 30 m (Robinson et al. 2017). All spatial data was processesed in the NAD 1983 datum using the Geodetic Reference System (GRS) 1980 spheroid with ArcGIS v10.6 (ESRI 2017).

Covariate correlation and collinearity We tested for correlation between bison grazing intensity and productivity by calculating the R² value and associated p-value, using a linear regression model. We tested for collinearity between bison grazing intensity and productivity by calculating the variance inflation factor between the two covariates.

Data processing

The dataset includes the raw data used for analyses in the paper. Formatting guidelines for the multi-species occupancy model can be found in Zipkin et al. 2010, and formatting guidelines for the multi species dependent double observer model can be found in Golding et al. 2017.

Multi-species occupancy model

To analyze occupancy across covariates, we used a Bayesian multi-species occupancy model (MSOM, Zipkin et al. 2010). This hierarchical model uses parameters estimated from the observation process to describe parameters in the ecological process of interest (occupancy). The observation process is almost always imperfect, but failure to detect a species can be distinguished from true absence of a species through repeated sampling of a plot. The field data is therefore corrected for detection, and then used to model the ecological relationship. The use of community-level hyper-parameters adds another level to the hierarchy because species-level parameters are assumed to be drawn from a common distribution (Zipkin et al. 2010). Hyper-parameters represent a mean response of all species, enable better estimates of species-specific estimates, and allow estimation for species that are rare in the data (Zipkin et al. 2010).

We modeled occupancy as a function of site (NBR or YNP), bison grazing intensity (patties), productivity (NDVI), an interaction between bison grazing intensity and productivity, and year, for species i, plot j, year k and site l, using the logit-link function. We included an intercept for each site to account for differences between NBR and YNP sites. We included the square of patty counts to test for our hypothesis that intermediate grazing would support highest species richness, and we included a covariate for year to account for other interannual variation in occupancy. 

We allowed detection (p) to vary by species and Julian date. We augmented the number of species in the community as described by Royle et al. (2007), to estimate the total number of species in the community, including species not observed during avian surveys. Species richness at plot j, year k, site l, is a derived parameter resulting from the sum of occupancy for each species, including augmented species.

Abundance model

To examine individual songbird species responses to bison grazing intensity, we used a multispecies dependent double-observer abundance model (MDAM; Golding et al. 2017). The MDAM incorporates the survey process of the dependent double-observer method into the likelihood statement, accounting for imperfect detection and creating detection-adjusted abundance estimates. 

We modeled abundance as a function of site (NBR or YNP), bison grazing intensity (patties), productivity (NDVI), an interaction between bison grazing intensity and productivity, and year, for species i, plot j, year k and site l, using the log function. We included an intercept for each site to account for differences between NBR and YNP sites. We included the square of patty counts to test for the possibility that a species has highest abundance at intermediate grazing intensities, and we included a covariate for year to account for other interannual variation in abundance. 

We allowed detection to vary by species and observer. We included observer in the detection model because observer detection rates will likely have a larger influence on abundance data than on occupancy data. We had six observers who completed a small number of bird surveys compared to the other 14 observers. These observers all received equivalent training before completing surveys however, because of the small number of surveys they completed, we combined the six observers into one observer identity to estimate detection.

For both the occupancy and abundance models, we assessed model convergence using the Gelman-Rubin statistic (Gelman and Rubin 1992) and visual inspection of trace plots, autocorrelation, estimate density, and the running mean. We ran the models in program R, version 3.4.3 (R Core Team 2017), using the package R2jags (Su and Yajima 2015) for the abundance model and the package jagsUI (Kellner 2017) for the occupancy model. We ran the occupancy model with 3 chains, for 60,000 iterations, and an adaptation period of 30,000 iterations. We ran the abundance model  with 3 chains, for 35,000 iterations, and a burn in period of 15,000 iterations.

Literature Cited

ESRI 2017. ArcGIS Desktop: Release 10.6. Redlands, CA: Environmental Systems Research Institute.

Gelman, A., and D. B. Rubin. 1992. Inference from iterative simulation using multiple sequences. Statistical Science 7:457-511.

Golding, J. D., J. J. Nowak, and V. J. Dreitz. 2017. A multispecies dependent double-observer model: a new method for estimating multispecies abundance. Ecology and Evolution 7:3425-3435. 

Kellner, K. 2017. jagsUI: A Wrapper Around 'rjags' to Streamline 'JAGS' Analyses. R package version 1.4.9. https://CRAN.R-project.org/package=jagsUI

Nichols, J. D., J. E. Hines, J. R. Sauer, F. W. Fallon, J. E. Fallon, and P. J. Heglund. 2000. A double-observer approach for estimating detection probability and abundance from point counts. Auk 117:393-408.

Paruelo, J. M., H. E. Epstein, W. K. Lauenroth, and I. C. Burke. 1997. ANPP estimates from NDVI for the central grassland region of the United States. Ecology 78:953–958.

Ralph, C. J., J. R. Sauer, and S. Droege. 1995. Monitoring bird populations by point counts. PSW GTR-149 :1–181. USDA Forest Service, Pacific Southwest Research Station, Albany, California, USA.

R Core Team 2017. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 

Robinson, N. P., B. W. Allred, M. O. Jones, A. Moreno, J. S. Kimball, D. E. Naugle, T. A. Erickson, and A. D. Richardson. 2017. A dynamic landsat derived normalized difference vegetation index (NDVI) product for the conterminous United States. Remote Sensing doi:10.3390/rs9080863

Royle, J. A., R. M. Dorazio, and W. A. Link. 2007. Analysis of multinomial models with unknown index using data augmentation. Journal of Computational and Graphical Statistics 16:67-85.

Sliwinski, M. S. 2011. Changes in grassland songbird abundance and diversity in response to grazing by bison and cattle in the northern mixed-grass prairie. Thesis. The University of Manitoba, Winnipeg, Manitoba, Canada.

Su, Yu-Sung and M. Yajima. 2015. R2jags: Using R to Run 'JAGS'. R package version 0.5-7. https://CRAN.R-project.org/package=R2jags

Tastad, A. C. 2013. The relative effects of grazing by bison and cattle on plant community heterogeneity in northern mixed prairie. Thesis. University of Manitoba, Winnipeg, Manitoba, Canada.

Tipton HC, Doherty PF, Dreitz VJ. 2009. Abundance and density of mountain plover (Charadrius montanus) and burrowing owl (Athene cunicularia) in eastern Colorado. Auk 126:493–499.

Tipton HC, Dreitz VJ, Doherty PF. 2008. Occupancy of mountain plover and burrowing owl in Colorado. Journal of Wildlife Management 72:1001–1006. 

Tucker, C. J., and P. J. Sellers. 1986. Satellite remote-sensing of primary production. International Journal of Remote Sensing 7:1395–1416.

Zipkin, E. F., A. J. Royle, D. K. Dawson, and S. Bates. 2010. Multi-species occurrence models to evaluate the effects of conservation and management actions. Biological Conservation 143:479–484.

Microsoft Excel is required to open the data files.