Metrics for conservation success: using the ‘Bird-Friendliness Index’ to evaluate grassland and aridland bird community resilience across the Northern Great Plains ecosystem

Aim Evaluating conservation effectiveness is essential to protect at-risk species and to maximize the limited resources available to land managers. Over 60% of North American grassland and aridlands have been lost since the 1800s. Birds in these habitats are among the most imperiled in North America, yet most remaining habitats are unprotected. Despite the need to measure impact, conservation efforts on private and working lands are rarely evaluated, due in part to limited availability of suitable methods. Location Northern Great Plains Methods We developed a novel metric to evaluate grassland and aridland bird community response to habitat management practices, the Bird-Friendliness Index (BFI), consisting of density estimates of grassland and aridland birds weighted by conservation need and a functional diversity metric to incorporate resiliency. We used the BFI to inform three assessments: 1) a spatial prioritization to identify ecologically significant areas for grassland and aridland birds, 2) estimation of temporal trends in grassland and aridland bird community resilience, and 3) evaluation of the effects of land management practices on grassland and aridland bird communities. Results The most resilient bird communities were found in the Prairie Potholes region of Alberta, Saskatchewan, northern Montana, and North Dakota, and the lowest BFI values in the southern and western regions of the Northern Great Plains. BFI values varied little over time on average, but trends varied regionally, largely in response to interannual relative variability in grassland and aridland bird densities. Main conclusions BFI values increased in response to simulated habitat management, suggesting that practices recommended for use in bird-friendly grassland habitat management plans will increase the abundance and resilience of the grassland and aridland bird community, and will be detected using the BFI. The BFI is a tool by which conservationists and managers can carry out accountable conservation now and into the future.

Additionally, species-specific abundance alone is rarely a reliable predictor of habitat quality (Johnson, 2008). Instead, metrics incorporating distributions, habitat relationships, and diversity of multiple species better encapsulate the response of ecological communities to conservation actions (Goyert et al., 2016;Nuttle, Leidolf, Burger, & Loiselle, 2003).
By incorporating suites of species using habitat in different ways in multi-species metrics, we can incorporate the variability among species in their habitat use, foraging preferences, and functional contributions to ecosystems (e.g., pest control, seed dispersal; Şekercioğlu, Wenny, & Whelan, 2016). Together, these suites of species and their associated functional traits and services influence the integrity (i.e., ecosystem structure, composition and function; Karr, 1981) and resilience (i.e., ability to resist or recover from disturbances; Scheffer, Carpenter, Dakos, & van Nes, 2015) of the larger ecological community (Fischer et al., 2007).
Functional traits such as diet, body size and habitat use influence species' ecological impacts, and communities composed of species with complementary suites of traits will have greater ecological integrity and resilience (Cardinale et al., 2012;McGill, Enquist, Weiher, & Westoby, 2006). Consequently, functional diversity measures incorporating functional traits and species abundance provide an indirect way to measure resilience and integrity (Standish et al., 2014), and as such are increasingly used in large-scale assessments of North American bird and ecological communities (Schipper et al., 2016;Schleuter, Daufresne, Massol, & Argillier, 2010).
In order to evaluate the success of habitat enhancement in conserving birds and their habitats on private lands in North American grasslands and aridlands, we developed a metric 6 evaluating bird community response to habitat management practices. This metric, the Bird-Friendliness Index (BFI), uses avian count, functional trait, and conservation status data together with remotely sensed environmental data to evaluate the capacity of a landscape to support an abundant, diverse, and resilient bird community. The BFI was designed to enable inference at multiple spatial and temporal scales, and incorporates standardization methods that facilitate spatial comparisons among North American grasslands and aridlands spanning expansive climatic and community gradients. Here we demonstrate the use of the BFI to inform three assessments: 1) spatial prioritization to identify ecologically significant areas for grassland and aridland birds, 2) estimation of temporal trends in relative grassland and aridland bird community resilience, and 3) evaluation of the effects of land management practices on grassland bird communities. We use the Northern Great Plains of North America as a case study, and highlight the insights the BFI provides into the plight of grassland and aridland birds and the integrity and resilience of the larger ecosystems in which they reside.

| Bird-Friendliness Index overview
The BFI takes into account species density, conservation status, and diversity of the entire community of grassland birds at a site. Bird abundance, and changes in relative abundance over time, are presently the most common tool for monitoring and evaluating landbird abundance (e.g., Rosenberg et al., 2016;Sauer et al., 2017), and form the foundation of the BFI. Including conservation status ensures that the BFI does not underestimate conservation need by diluting the contribution of rare species (Beissinger, 2000). Lastly, the BFI also includes a measure of functional diversity as a measure of the intactness and resilience of both the grassland and aridland bird community and their larger ecological community (Fischer et al., 2007).
The BFI is the sum product of estimated avian density and conservation status times a measure of functional diversity based on species traits (i.e., diet, foraging strata, and body mass), as follows: It is calculated in six steps, all of which are calculated at the resolution of 1 km 2 grid cells. This resolution was chosen because it represents the smallest scale at which the avian survey data used here are summarized (see Avian data), and it is large enough to contain entire home ranges of most study species. First, densities for each species were estimated for surveyed grid cells using point counts that include one or more auxiliary methods enabling estimation of detection-corrected density (e.g., distance or time of first detection; Matsuoka et al., 2014).
Second, individual species models relating density to remotely-sensed environmental covariates were built and used to predict densities across the study area. Third, species densities were multiplied by a conservation score for each species, so that threatened and rare species carry a greater weight. Fourth, conservation status-corrected densities for each species were summed within each grid cell. Fifth, summed conservation-corrected densities for each grid cell were multiplied by an index of functional diversity calculated using the predicted densities ( Step 2) and corresponding functional traits for the species present within a grid cell. weather, density-dependent cycles, and background population-level trends.

| Study area
Our study area encompasses the Northern Great Plains ecosystem, as defined by the West-  (Comer et al., 2018). The study area also includes extensive arid sagebrush steppe habitat intermixed with grasslands, ranging from 1-25% of regional land cover in Montana and the Dakotas, to >75% cover in northeastern Wyoming (Connelly, Knick, Schroeder, & Stiver, 2004).

| Study species
We used the 2016 State of North America's Birds species assessment to identify 107 bird species that use grasslands or aridlands as their primary or secondary breeding habitat (North American Bird Conservation Initiative, 2016). We selected grassland and aridland birds because the two habitats are intermixed at relatively fine spatial scales in this region, and private landholdings on which management practices are implemented often include both habitat types and, consequently, communities. This list formed the pool of species evaluated for inclusion in the BFI contingent on data availability in the study area (see Avian data).

| Avian data
We used avian point count data collected under Bird Conservancy of the Rockies's Integrated Bird Monitoring in Bird Conservation Regions (IMBCR) program (Woiderski et al., 2018). The IMBCR uses a hierarchically nested sampling design in which spatially balanced random stratification is used to select 1 km 2 grid cells within BCRs. Each grid cell includes up to 16 point count locations evenly spaced at distances of 250 m. At each point, five-(2009) or six-minute (2010)(2011)(2012)(2013)(2014) unlimited-distance point count surveys were conducted between one-half hour before and five hours after sunrise during the late spring and early summer. Observers recorded species, age (where possible), first minute of observation, type of detection (call, song, visual), day, time of day (2010-2014 only), and primary habitat type at the point, and measured distance to all individuals and flocks using a laser rangefinder. We excluded flyover observations, juvenile birds, and observations with missing minutes or distances.
We assigned IMBCR survey points to spatial strata defined as unique combinations of BCRs and states or provinces, as used in hierarchical modeling of Audubon Christmas Bird Count data ( Figure 1; Soykan et al., 2016). We used a subset of the IMBCR data for strata in which Of these, 34 species had sufficient data from which to estimate density (see Density estimation) and were included in the BFI (Table 1). These included 26 species that use grassland habitats and 12 species that use aridlands as a major breeding habitat. Four species (Bell's Vireo, Lark Sparrow, Loggerhead Shrike, and Western Kingbird) use both grasslands and aridlands as major breeding habitats (Table 1), and many species use both habitats facultatively. The final dataset for these 34 species included 45,152 observations of 99,842 birds.

| Environmental data
We used 17 environmental variables in individual species models to predict avian densities across the study area (see Density prediction; Table S1). All environmental variables were sampled for 799,015 grid cells of 1 km 2 , covering the study area and aligning with the avian density grid cells. We used statistically downscaled climate normals from 1981-2010 derived from the Climatic Research Unit Timeseries 3.22 dataset to represent long-term climatic conditions (Wang, Hamann, Spittlehouse, & Carroll, 2016). We obtained annually-derived ecosystem metrics (e.g., net primary productivity, soil moisture, evaporation) from the NASA CASA terrestrial ecosystem model (Potter et al., 1993;Potter, Klooster, Huete, & Genovese, 2007 Renfrew et al., 2013;Wiens, 1974).

| Density estimation
As the avian dataset included both time of first detection and distance of observations, we used a formulation of time removal and distance sampling models developed by Sólymos et al., (2013) and implemented in R package detect (Sólymos, Moreno, & Lele, 2016). We estimated point-level density corrected by two components of detection probability, availability and perceptibility, using conditional multinomial maximum likelihood estimation. Availability -the probability that a bird provided a visual or auditory cue during sampling and was thus available to be detected -was estimated using the minute of first observation. Each individual was counted only once, thus individuals were 'removed' once detected. Perceptibility -the conditional probability that birds available for detection were actually detected -was estimated as a function of distance from observer (Sólymos et al., 2013).
We combined data from 2009-2014 in a single model per species to maximize the number of species with sufficient data for modeling, but allowed availability and perceptibility to vary among years. We obtained model convergence for species with 23 or more observations during all years. The number of observations per species ranged from 23 (Canyon Wren) to 34,857 (Western Meadowlark), with a median value of 544 ± 6,533 SD individuals. To account for species absence or non-detection, we created zero observation records for all surveys at which each species was not observed.
We grouped observations into three time periods: 0-2 minute, 2-4 minute, 4-6 minute (4-5 minutes in 2009), and three distance bins: 0-50 m, 51-100 m, and >100 m. Availability was estimated as the singing rate, defined as the probability a bird sang or gave a visual cue during a survey. Singing rate was collapsed across distance bins and modeled as a log-linear function.
We controlled for temporal and environmental variables affecting singing rate by including linear terms for year, primary habitat, and spatial strata, and both linear and quadratic terms for day of year, in addition to the intercept. We calculated perceptibility using the effective detection radius, defined as the distance beyond which as many birds were detected as remained undetected within. We controlled for temporal and environmental variables affecting perceptibility by including linear terms for year, primary habitat, spatial strata, and surveyor identity, and both linear and quadratic terms for day of year, in addition to the intercept.
Models including all combinations of covariates, and no covariates, were built and model selection was conducted using AIC (Burnham & Anderson, 2002).
Singing rate and effective area sampled (in ha), defined as the area within the effective detection radius, were calculated at the point level then used to estimate a correction factor for each survey accounting for imperfect detection. Correction factors were log transformed and entered as offsets in generalized linear mixed models calculated using R package lme4 (Bates, Mächler, Bolker, & Walker, 2015). We modeled detection-corrected annual density at the point level using a Poisson distribution and included primary habitat as a fixed effect, and grid cell and year as random effects. We left-truncated density estimates at 0.0008 birds/ha, equivalent to the density estimate for a point with a single bird observed at the maximum observed distance in the dataset (2000 m). Annual densities at the 1 km 2 grid cell level were estimated by taking the median of annual point-level densities (birds/ha) within each surveyed grid cell multiplied by 100.

| Species habitat modeling
We built species habitat models, with density as the response variable, for each species using boosted regression trees (BRTs), a machine learning approach that is ideal for modeling complex species-environment relationships with multiple, and often highly correlated, predictors (Elith, Leathwick, & Hastie, 2008). As many species had skewed density distributions with absences in many grid cells, we implemented a hurdle model approach in which we separately modeled occurrence and density. The response variables were presence/absence for the presence model, and detection-corrected abundance scaled to 1 km 2 for the density model.
For the density model, we first rounded densities and used a Poisson distribution. If these models failed to converge, we log-transformed estimated densities to improve fit and used a Gaussian distribution. The predictor variables included the 17 environmental variables described above, plus year (Table S1). Models were fit using R packages dismo (Hijmans, Phillips, Leathwick, & Elith, 2015) and gbm (Ridgeway, 2017).
BRTs use three regularization parameters to shrink the number of terms in the final model and thus avoid overfitting. Learning rate shrinks the contribution of each tree to the full model, bag fraction specifies the proportion of data selected from the training set at each iteration, and tree complexity sets the number of nodes and, thus, level of interactions between predictors (Elith et al., 2008). We iteratively tuned these parameters to optimize model fit while ensuring a minimum of 1,000 trees using default parameter ranges recommended by Elith et al., (2008): learning rate 1.0e -10 -0.1, bag fraction 0.5-0.75, and tree complexity 1-3. At each step we used 10-fold cross-validated area under the curve (AUC; occurrence model only) and residual deviance (both models) to select the optimal parameter value.
Species distribution models such as these are susceptible to spatial autocorrelation, which can result in biased presence/density-environment relationships (Boria, Olson, Goodman, & Anderson, 2014). We took two steps to reduce residual spatial autocorrelation in the model.
First, we conducted geographic filtering by overlaying a 10 km 2 grid and randomly selecting only one grid cell from each 10 km 2 grid for inclusion in the analyses (Boria, Olson, Goodman, & Anderson, 2014). We chose this scale based on Wilsey et al., (2019), which tested three geographic resolutions for filtering and found that 10 km produced the greatest model performance. Second, we used spatially stratified cross-validation by dividing the data sets into 12 gridded bins by latitude and longitude, and withholding one latitudinal bin for testing at each fold (Roberts et al., 2017). We tested for residual spatial autocorrelation in the final models using Moran's I, calculated in R package ape (Paradis, Claude, & Strimmer, 2004), and evaluated model performance using cross-validated AUC and TSS (presence only), deviance explained, and correlation. We evaluated the influence of each predictor variable by calculating mean relative importance values for all 34 species.
Because of the randomness introduced by the geographic filtering, as well as the iterative BRT estimation process itself, we produced 25 bootstrapped geographically filtered datasets. We selected 25 bootstrap iterations as there was a maximum of 24 surveyed grid cells within a single 10 km 2 filtering grid. Each of the 25 presence/absence and density models were used to generate predicted probability of occurrence or density across the entire study area using the 18 predictor variables. For each species, we calculated median predicted probability of occurrence and density for each grid cell and across the 25 bootstraps. We calculated a minimum probability of occurrence threshold using maximum sensitivity + specificity in R package SDMTools (VanDerWal et al., 2014), and used this to mask median predicted density such that densities were retained only at grid cells that surpassed the threshold hurdle.

| Conservation scores
We incorporated the breeding season Combined Conservation Score from the State of North America's Birds (North American Bird Conservation Initiative 2016), as these scores were calculated consistently for all native bird species across North America trends using a process developed over many years and updated in response to peer-review. Combined Conservation Sscores for ranged from 7 -20 based on population size, distribution, threats and trends (Panjabi, Blancher, Dettmers, & Rosenberg, 2012). Species-specific breeding season Combined Conservation Scores were multiplied by predicted densities for each 1 km 2 grid cell to produce conservation-corrected densities, which were then summed across all 34 species per grid cell.
We explored scaling the Combined Conservation Scores by dividing each score by the root mean square (range: 0.60 -1.37). BFI values calculated with scaled Combined Conservation Scores were directly correlated with BFI values calculated with raw Combined Conservation Scores in each year (ρ = 0.98 -0.99, p < 0.0001). Therefore, we used the raw scores in BFI estimation.

| Functional Diversity
Functional diversity indices take a variety of forms, differing in the resolution of the trait information they include and the use of presence versus abundance data (Schipper et al., 2016;Schleuter et al., 2010). We implemented a two-step approach that used environmental filtering to group species into functional groups (Table 1) using both high-resolution continuous traits (i.e., body mass, proportion dietary composition) and discrete categorical traits (i.e., foraging strata) obtained from EltonTraits 1.0 (Wilman et al., 2014). We converted foraging strata from categorical to integer form, subdivided each functional trait into two bins based upon the distribution of values for that trait, and constructed a series of functional groups ). These functional groups were then used in place of species identity to calculate a functional Shannon's diversity index, as both abundance and evenness have been shown to influence ecosystem level processes independent of richness (Petchey & Gaston, 2006).

| BFI calculation and standardization
BFI values were calculated by multiplying the summed conservation-weighted densities by functional diversity for each 1 km 2 grid cell. Raw values were standardized by scaling from zero to one using a logistic distribution, to accommodate the lower limit of zero and the occasional, exceptionally high BFI value. For the spatial prioritization and evaluation of temporal trends we standardized using all grid cells to enable comparison across the study area as a whole.
However, for the assessment of bird community response to simulated management we standardized using only grid cell values from the surrounding spatial strata (South Dakota -Badlands and Prairies; Figure 1) to highlight effects of management relative to nearby areas with similar climate and landcover.

| Sensitivity analyses
We conducted a sensitivity analysis to assess the relative influence of species density, conservation score, and functional diversity on BFI values. For each of 1,000 bootstrap iterations we sampled each species' annual densities from a log-normal distribution using estimated means and standard deviations increased by 10% with a total sample size equivalent to the number of 1 km 2 grid cells (n = 799,015). For each of 1,000 bootstrap samples of conservation scores we drew species values from a uniform distribution with the minimum (4) and maximum (20) conservation weight serving as the bounds of the distribution. For functional diversity we calculated the mean and standard deviation of the functional evenness values across the entire NGPE in 2011 and again increased the standard deviation by 10%. For functional diversity we sampled values 1,000 times from a beta distribution (which produces values ranging from 0-1) using estimated annual means and estimated standard deviations increased by 10%. We calculated a resampled BFI for each iteration, and estimated effect size as the degree of overlap between the distribution of BFI values calculated from the actual data and the distribution of BFI values for each bootstrap sample to determine effect size. Degree of overlap was calculated by using integration of the area of overlap between the two curves. We subtracted the mean overlap from one so that the greatest effect size would correspond to the component with the smallest degree of overlap.

| Temporal trends and management case study
We analyzed relative temporal trends in BFI values by spatial strata (i.e., unique combinations of states/provinces and Bird Conservation Regions; Figure 1) and year using general linear models calculated in R package nlme ( (Lenth, 2018) using Tukey adjustments for multiple comparisons, and plotted using ggplot2 (Wickham, 2016).
We simulated an evaluation of the effects of management actions on a theoretical private property. In these simulations, we modeled the effects of common bird-friendly habitat management actions including reducing cropland cover, CO2, evapotranspiration, and N2O, and increasing litter biomass, net primary productivity, soil moisture, proportion grassland, grassland cohesion, and grassland patch area. These covariates were selected for the simulation modeling based on known relationships with bird occurrence and abundance, or because they respond to -and therefore serve as proxies for -reduction of cattle rotation rates and consequently grazing intensity, a common recommendation of bird-friendly habitat management plans (Brennan & Kuvlesky, 2005;Johnson et al., 2019;North American Bird Conservation Initiative, 2016;Undersander, Temple, Bartelt, Sample, & Paine, 2000). Grassland cover, patch area, and connectivity; primary productivity; and litter biomass are important predictors of grassland birds, and associated with higher occurrence rates and densities Fisher & Davis, 2010;Harrower et al., 2017;Renfrew et al., 2013). Similarly, many grassland birds have positive associations with soil moisture and negative associations with drought (evapotranspiration), particularly in the drier western Great Plains (George et al., 1992;Niemuth et al., 2008;Wiens, 1974). These environmental variables may be influenced by management practices that increase cover, connectivity, and density of grasses. Conversely, many grassland birds have lower densities in croplands or avoid them altogether (Best et al., 1997;Pool et al., 2014). Cattle grazing increases nitrous oxide (N2O) and carbon dioxide (CO2) emissions, therefore these variables serve as proxies for cattle grazing intensity (Allard et al., 2007;Oenema, Velthof, Yamulki, & Jarvis, 1997).
Covariate changes were simulated to accumulate incrementally over time at a rate of 10% per year with some random variation included. These simulated environmental conditions were then used to re-calculate bird densities and functional diversity, and re-estimate the BFI each year within that property. We simulated habitat improvement beginning in 2011 to simulate a before-after control-impact design with monitoring for two years prior to management. We then compared changes in simulated BFI values post-management to actual, estimated BFI from 2009 to 2014. 0.00 -0.75; Table S2). All Moran's I values were ≤ 0.05, indicating that the predictors fully accounted for spatial autocorrelation in the data (Table S2).

| Estimated density
Climate and landcover variables explained the most variation in presence/absence and density of the grassland and aridland birds studied here. Climatic moisture deficit explained the most variation in occurrence (mean ± SE: 18.51 ± 2.80%), followed by terrain ruggedness index (mean ± SE: 12.39 ± 2.52%) and proportion cropland (mean ± SE: 11.22 ± 2.55%; Figure 2A).
Proportion cropland cover explained the most variation in density (mean ± SE: 13.07 ± 3.37%), followed by climatic moisture deficit (mean ± SE: 11.42 ± 3.26%) and terrain ruggedness index (mean ± SE: 11.13 ± 3.24%; Figure 2B). Year explained the least variation in both occurrence and density. Though densities and occurrence frequencies varied among years, environmental predictors that varied across both time and space explained more variation than year alone.
Relative variable importance varied among species consistent with species-specific habitat preferences, e.g., White-throated Swift density was explained most by terrain ruggedness, while occurrence and density of grassland birds like Upland Sandpiper were explained most by grassland area, patch area, and cohesion ( Figures S1-S34).

| Bird-Friendliness Index
BFI values showed large-and fine-scale regional variation. The most resilient grassland and aridland bird communities were found in the Prairie Potholes region of Alberta, Saskatchewan, northern Montana, and North Dakota ( Figure 3). Conversely, the lowest BFI values were found in the southwestern regions of the study area, notably southern Montana, Wyoming, and southern South Dakota. The regional patterns were generally consistent over time, though high BFI values were more concentrated across the Prairie Potholes in 2013-2014 than earlier years.

| Sensitivity analyses
Sensitivity analyses revealed that BFI values were most sensitive to variation in bird densities ( Figure 4). BFI calculated with resampled bird densities overlapped with original BFI values by 87.11 ± 0.03%, while resampling conservation scores produced 96.43 ± 0.02% overlap and resampling functional diversity produced 94.41 ± 0.02% overlap (Figure 4).

| Spatial prioritization
We identified top-priority areas supporting resilient grassland and aridland bird communities as grid cells within the top 20% of BFI values in all years (2009)(2010)(2011)(2012)(2013)(2014). A total of 427,070 km 2 fell within the top 20%, concentrated in the Prairie Potholes, North Dakota, northern South Dakota, and central Nebraska ( Figure 5).

| Temporal trends in relative grassland and aridland bird community resilience
BFI values were consistent over time across the Northern Great Plains study area, with mean values of 0.50 ± 2.0 e -4 SE in each year (slope = 8.9e -5 ± 5.4e -5 , p = 0.00; Table 2). However, the best-fitting model had an interaction between year and strata, indicating that both mean BFI values and temporal trends varied between strata ( showed slight declines ( Figure 6).

| Case study: effects of land management practices on grassland and aridland bird
communities BFI values representing grassland bird community response to simulated management were 11% higher in 2014 than BFIs estimated from observed data ( Figure 7). Additionally, BFI values with simulated bird-friendly habitat management significantly increased over the six-year period (slope = 0.036 ± 0.005 SE, p = 0.002), while estimated BFI values (without bird-friendly habitat management) did not change during the same time period (slope = 0.014 ± 0.007 SE, p = 0.125). This suggests that practices recommended for use in bird-friendly grassland and aridland habitat management plans will increase the abundance and resilience of the bird community, and will be detected using the BFI.

| DISCUSSION
The ability to quantify the impacts and success of conservation and management actions is crucial to the adaptive management cycle (Nichols & Williams, 2006). Conservation Area (Comer et al., 2018). Similarly, many Grassland Priority Conservation Areas (Pool & Panjabi, 2011) were identified across the Prairie Potholes, Dakotas, as well as in a few isolated pockets in northeastern Wyoming coinciding with top-priority areas identified by the BFI. A recent spatial prioritization analyzing current and future climate suitability also identified the Prairie Potholes and North Dakota as high-priority areas, similar to our findings (Grand, Wilsey, Wu, & Michel, 2019).
Conversely, the lowest BFI values were found in the southwestern regions of the study area, notably southern Montana, Wyoming, and southern South Dakota. This finding is consistent with a recent grasslands prioritization based primarily on climate suitability (Grand et al., 2019). These regions tend to be drier with more sparse vegetation, as well as lower soil moisture and productivity. While these regions provide critical habitat for aridland bird species such as Greater Sage-grouse (Burkhalter et al., 2018;Connelly, Hagen, & Schröder, 2011), they are less suitable for grassland birds that prefer lusher and more productive habitats Fisher & Davis, 2010;Harrower et al., 2017;Renfrew et al., 2013). The BFI was calculated using data for 26 predominantly grassland-breeding birds and 12 predominantly aridlandbreeding birds. Thus BFI values were lower in this region overall as many grassland birds have higher mean densities than aridland birds (Table 1), and the BFI is most sensitive to changes in bird densities (Figure 4).
Evaluation of the drivers of occurrence and abundance revealed that moisture availability, terrain ruggedness, and availability of suitable habitat (grassland, shrubland) were the most important predictors. This conclusion fits well in line with what we already know of the Northern Great Plains ecosystem and grassland bird biology. The Northern Great Plains are characterized by high variability in terms of moisture, undergoing periodic bouts of large-scale drought or decreased precipitation particularly in southern and western regions, which in turn has large-scale impacts on vegetation growth (Laird, Fritz, Maasch, & Cumming, 1996). Both short-term (~2-3 years) and long-term (24 years) studies have shown large impacts on many grassland bird species, with most responding negatively to decreases in annual precipitation and moisture levels, indicating their sensitivity to drought (George et al., 1992;Gorzo et al., 2016;Niemuth et al., 2008).
Grassland birds are also associated with low elevation sites with limited topographical variation (Wiens, 1973;Niemuth et al., 2017). In addition to being dry, the western region of the Northern Great Plains study area has greater variability in terrain and more high elevation sites, contributing to the lower BFI values. Finally, the importance of grassland, shrubland, and cropland cover in explaining variation in bird distribution and abundance is consistent with grassland and aridland bird ecology. Grassland cover is an important predictor of grassland birds associated with higher occurrence rates and densities, while shrubland and cropland cover are associated with lower occurrence rates and densities Fisher & Davis, 2010). Conversely, many aridland species such as Greater Sage-grouse, Bell's Vireo, and Loggerhead Shrike are positively associated with shrubland habitats, increasing the importance of this predictor (Budnik, Thompson, & Ryan, 2002;Connelly et al., 2011). We were unable to incorporate annual landcover in this case study due the lack of robust annual landcover classifications that span the US / Canada border. However, future applications of the BFI should take advantage of annual landcover data, where available, because of these inherent sensitivities.
The BFI was designed for spatial comparisons within a single time period, or evaluating site-level trends in resiliency relative to other sites or the region as a whole, but not for evaluating population trends at the species or region level. As such, the standardization procedure was implemented, which resulted in relatively stable BFI values over time across strata or the study area as a whole, despite interannual variation in grassland and aridland bird density. Standardization enables isolation of local-scale effects, such as bird response to habitat management or regional differences in land use, from large-scale processes such as long-term population-level declines due to shared drivers such as climate oscillations or wintering ground effects (Gorzo et al., 2016;Macías-Duarte et al., 2017). The effects of this standardization procedure were apparent in the habitat management case study. Though neither mean BFI across the study area nor mean observed BFI at the case study site increased over six years, simulated habitat management increased grassland and aridland bird densities and diversity relative to the surrounding landscape, producing a significant increasing trend in BFI. However, the standardization step could be removed for studies with an objective of tracking long-term trends in ecological resilience across large spatial scales, rather than evaluating local response to conservation actions.
Conservation efforts should aim to do more than prevent the extinction of species, but rather should be aimed at preventing species from becoming threatened in the first place, as well as provide conditions that enhance ecosystem processes (Rodrigues, 2006). By identifying early on which species and communities are doing well or poorly, and where, the BFI can pinpoint priority strongholds for conservation, or opportunities for restoration, both of which may contribute to population stabilization or even growth. To be able to do this effectively requires indicators that are rigorous, repeatable, and easily understood (Balmford et al., 2005).
National Audubon Society deploys the BFI as an accountability metric for the Conservation Ranching Program (https://www.audubon.org/conservation/ranching), a market-based conservation solution that aims to conserve imperiled grassland and aridland bird species and habitats through partnership with private landowners. Other public/private grassland and aridland conservation efforts, e.g., World Wildlife Fund's Grasslands program (https://www.worldwildlife.org/habitats/grasslands), Sodsaver (enacted as part of the Agricultural Act of 2014) and Natural Resources Conservation Service's EQIP (https://www.nrcs.usda.gov/wps/portal/nrcs/main/national/programs/financial/eqip/) may also benefit from the accountability provided by this community ecological resilience metric.
The BFI is a tool by which conservationists and managers can carry out accountable conservation now and into the future.      Locations of strata are shown in Figure 1.