DOI: 10.5061/dryad.ksn02v757
title: Avian species richness and abundance shows stronger response to bison grazing intensity than to ecosystem productivity

description: 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.
description: <em>Study
 sites</em> 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<sup>2</sup> NBR
  boundary. YNP is managed by the U.S. National Park Service. The YNP study
 site took place in 890 km<sup>2</sup> 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. <em>Study
 species</em> 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 (<em>Spizella
 pallida</em>), grasshopper sparrow (<em>Ammodramus
 savannarum</em>), vesper sparrow (<em>Pooecetes
 gramineus</em>), and western meadowlark (<em>Sturnella
 neglecta</em>). The sagebrush-steppe species were Brewer’s sparrow
 (<em>Spizella breweri</em>), green-tailed towhee
 (<em>Pipilo chlorurus</em>), and sage thrasher
 (<em>Oreoscoptes montanus</em>). 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. <em>Sampling
 frame</em> <em>                </em>We
 selected survey plots of 250 x 250 m<sup>2</sup> 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. <em>Avian
 surveys</em> <em> </em>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. <em>Explanatory
 Variables</em> <em>Bison grazing intensity data
 </em>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). <em>Productivity data</em> 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). <em>Covariate
 correlation and collinearity </em>We tested for correlation between
 bison grazing intensity and productivity by calculating the
 R<sup>2</sup> 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. <em>Data
 processing</em> 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. <em>Multi-species occupancy
 model</em> 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 (<em>NBR</em> or
 <em>YNP</em>), bison grazing intensity
 (<em>patties</em>), productivity (<em>NDVI</em>),
 an interaction between bison grazing intensity and productivity, and year,
 for species <em>i</em>, plot <em>j</em>, year
 <em>k</em> and site <em>l</em>, 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
 (<em>p</em>) 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 <em>j</em>, year <em>k</em>, site
 <em>l</em>, is a derived parameter resulting from the sum of
 occupancy for each species, including augmented species. <em>Abundance
 model</em> 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 (<em>NBR</em> or
 <em>YNP</em>), bison grazing intensity
 (<em>patties</em>), productivity (<em>NDVI</em>),
 an interaction between bison grazing intensity and productivity, and year,
 for species <em>i</em>, plot <em>j</em>, year
 <em>k</em> and site <em>l</em>, 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 (<em>Charadrius montanus</em>) and burrowing owl (<em>Athene cunicularia</em>) 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.
description: Microsoft Excel is required to open the data files.