A large-scale perspective for managing prairie avifauna assemblages across the western US: influences of habitat, land ownership and latitude

Study area

The study area is approximately 140,000,000 ha within the prairie landscape of the western United States defined by the Great Plains and North American Deserts ecoregions (Commission for Environmental Cooperation (CEC), 1997; Fig. 1). This landscape is comprised of seven western US states and six Bird Conservation Regions (BCRs)—distinct ecological regions with similar biotic and abiotic characteristics and natural resource management issues (US North American Bird Conservation Initiative (US NACBI), 2000). Management of the lands within the study area is by private landowners (51.6%), US Department of Interior—Bureau of Land Management (BLM; 22.2%) and other federal agencies (e.g., US Department of Agriculture—Forest Service (USFS; 16.5%)), Native American tribes (2.8%), and US states’ management agencies (3.1%; Table 1).

Sampling: design-based

The bird data consisted of point count surveys collected from 2009–2012 as part of a regional, multi-agency landbird monitoring collaboration, the Integrated Monitoring in Bird Conservation Regions (IMBCR) project. This design-based program covers all land area within the monitoring landscape, regardless of land ownership, habitat, or other factors (e.g., roads, terrain). Bird surveys were conducted across areas designated by fixed attributes such as ownership boundaries, state boundaries or BCR boundaries (White et al., 2012), hereafter referred to as strata. Within each stratum, grid-scale sampling plots were selected using a generalized random tessellation stratified sample, a spatially-balanced sampling algorithm (Stevens & Olsen, 2004; Theobald et al., 2007). Each sampling plot measured 1 km² in area (Fig. 1) and contained 16 bird point count stations evenly spaced at 250 m to form a 4 × 4 grid of points. Each 1 km² sampling plot represented the grid-scale, while a single point count station within a given plot represented the point-scale. Observers conducted a 6-min survey at each point count beginning 30 min prior to sunrise and concluding prior to 11 AM (local time) once per year. For every bird visually detected, the species, horizontal distance to the detected bird and minute interval were recorded; however, we truncated data into detection/non-detection data. We overlaid land cover data from the USGS National Gap Analysis Program (US Geological Survey, Gap Analysis Program (US GAP), 2011) in ArcGIS 10 to select sampling plots containing prairie habitats (<30% Forested & Woodland; n = 413 sampling plots) for use in the analysis.

Guild delineation

Patterns of prairie bird occurrence and habitat associations can be examined at various levels, ranging from the individual species to the species assemblage. The guild concept is one way of classifying species assemblages; however, “guild” is a term that lacks consistent definition in the literature (Verner, 1984; Szaro, 1986; Simberloff & Dayan, 1991). We defined a guild specifically as either grassland (n = 17 species; Vickery et al., 1999) or sagebrush (n = 3 species; Paige & Ritter, 1999; Hanser & Knick, 2011) obligates (Table 2). We explicitly chose obligates to reflect their specific habitat preferences as opposed to generalists which reflect a different set of processes. The 20 species included in our analyses are designated as a conservation priority (e.g., species petitioned for US federal listing; US state-level species of concern; USFS Management Indicator Species). We limited our analyses to passerines and upland shorebirds (i.e., Charadridae) detected ⩾ 10 times in ⩾ 1 year.

Environmental covariates

We examined the relationship of land ownership, habitat and latitude to regional species richness of grassland and sagebrush obligate avifauna. For each of the 413 sampling plots, we calculated the percent area covered by each type of land ownership (private and public) and habitat (grassland, sagebrush, and other). Public ownership was defined as lands managed by US federal or state agencies. We also included lands managed by North American tribes in the public ownership category, given collective management strategies for many stakeholders. Private ownership included all lands not in public ownership. Habitat was assessed by calculating the percentage of grassland and sagebrush within each sampling plot was based on the National Vegetation Classification (NVC) (i.e., Class/Subclass, Form) and the Ecological Systems Classification (US Geological Survey, Gap Analysis Program (US GAP), 2011). Grassland was defined by four NVC Ecological Systems: (1) mixed grass prairie, (2) shortgrass prairie, (3) sand prairie/sandhill steppe and, (4) montane and foothill shrubland and grassland. Sagebrush was defined by two NVC Ecological Systems: (1) sagebrush shrubland and (2) steppe and semi-desert scrub shrub. Additional habitats in the sampling plot—cultivated cropland and pastures, wetland and open water, forest and woodlands, and developed—were grouped as ‘other’ habitat. Lastly, we used the centroid of each sampling plot to calculate latitude. The latitudinal gradient theory of poleward biodiversity decline is one of the longest recognized (MacArthur & Wilson, 1963), albeit not completely understood, patterns in ecology (Hillebrand, 2004; Mittelbach et al., 2007; Salisbury et al., 2012). Studies conducted at multiple scales, on many taxa and across different biogeographic regions have proposed numerous explanations for this gradient (Hillebrand, 2004; Mittelbach et al., 2007; Salisbury et al., 2012; Karanth et al., 2014). Latitude functions as a surrogate for a suite of environmental factors, thus, interlinking and confounding the distinct underlying hypotheses (Rahbek, 2005). Therefore, using latitude as a predictor allowed us to make inference on covariates of interest, while reducing the parameter space necessary to explain variation in avian occurrence.

Analysis: model-based

We used a multi-species occupancy/richness hierarchical model as described by Dewan & Zipkin (2010), Zipkin, Dewan & Royle (2009) and Zipkin et al. (2010) to estimate the number of species within each guild present at each sampling plot. We combined individual species occurrence models into a single model by assuming that species covariate effects arise from a common distribution, allowing for more precise estimates of occupancy (Zipkin, Dewan & Royle, 2009; Kéry, 2010). In the context of estimating occupancy/richness, hierarchical models can distinguish absence from non-detection by incorporating models that specify presence versus absence as one process and detection versus non-detection as another process (Zipkin, Dewan & Royle, 2009).

For species-level models, we assumed that the the latent state of occurrence for a given species z_i, is a Bernoulli process where the probability that species i is present at sampling plot j(z_i,j = 1) is ψ_i,j (MacKenzie & Kendall, 2002). We modeled species (i) and site (j) level heterogeneity in occurrence probabilities using a logit link function with relevant covariates such that: ${\displaystyle logit\left({\psi }_{ij}\right)={\alpha }_{0,i}+{\alpha }_{1i}{PUBLIC}_j+{\alpha }_{2i}{GRASS}_j+{\alpha }_{3i}{SAGE}_j+{\alpha }_{4i}{LAT}_j+{\tau }_i,}$ where α₀ is the intercept and α₁ − α₄ are the effects of the environmental covariates on species i and τ_i is a random effect for species i. We used covariates on the percent of grassland (GRASS) and sagebrush (SAGE) habitat along with public land ownership (PUBLIC) and latitude (LAT).

Because species are imperfectly detected during sampling (MacKenzie & Kendall, 2002), we assumed that true occurrence, (z_i,j = 1), is a latent process that is only partially observable. If an observer detected species i at sampling plot j during point count k, denoted x_i,j,k = 1, then it can be determined that z_i,j = 1. However, if a species is not detected it could be that the species was either absent or undetected during sampling.

We assumed detection probability (p_i) of species i was unlikely to vary by land ownership. However, landcover is likely to influence variation in species detectability through vegetative structure and available habitat, so we included covariates for grassland and sagebrush cover and included latitude to incorporate other unmeasured sources of variation in species specific models: ${\displaystyle logit\left(p_{i,j}\right)={\beta }_{0,i}+b_{1,i}GRAS{S}_j+b_{2,i}SAG{E}_j+b_{3,i}LATITUD{E}_j.}$ The species-specific occurrence and detection probabilities are related by a community-level hierarchical component that assumes each probability parameter (e.g., α_iPUBLIC_j) arises from a common distribution. Model parameters were estimated using a Bayesian analysis of the model with naive prior distributions for coefficients (Normal[μ = 0, σ² = 1, 000]) and variance of random effects (Gamma[α = 0.1, β = 0.1]) as there was no previous information on grassland bird occurrence at this spatial scale. We used JAGS (Plummer, 2003) using the R2jags package in R (R Development Core Team, 2013). We ran 4 Markov chain Monte Carlo chains of 20,000 iterations each with the first 1,000 discarded for burn-in. We used the Gelman–Rubin R statistic to assess chain convergence (Gelman & Rubin, 1992).

To assess the fit of data to our model we conducted a posterior predictive check by calculating a Bayesian p-value. Specifically, we calculated deviance for each MCMC iteration (s) following methodology from Broms, Hooten & Fitzpatrick (2016). However, we calculated deviance for each species using the observed, ${\displaystyle D_i^{\left(s\right)}=-2\sum _{j=1}^{J}log\left(x_{i,j}|{\psi }_j^{\left(s\right)},p_j^{\left(s\right)}\right)}$ and predicted data (${\tilde{x}}_{i,j,k}\sim Bern\left[z_{i,j}{p}_{i,j,k}\right];where{z}_{i,j}\sim Bern\left[{\psi }_{i,j}\right]$). ${\displaystyle {\tilde{D}}_i^{\left(s\right)}=-2\sum _{j=1}^{J}log\left({\tilde{x}}_{i,j}|{\psi }_j^{\left(s\right)},p_j^{\left(s\right)}\right).}$ We calculated a Bayesian p-value as the proportion of MCMC samples in which observed deviance was greater than that calculated using predicted data for each species. Generally, a model has poor fit to the data if the p-value is >0.95 or <0.05 (Broms, Hooten & Fitzpatrick, 2015).

Lastly, we developed a predictive map of occurrence across the study area by applying species-specific coefficient estimates to spatial covariate layers. Individual species occurrence predictions were summed across cells to generate single a spatial prediction of species richness.