Metrics for conservation success: Using the “Bird‐Friendliness Index” to evaluate grassland and aridland bird community resilience across the Northern Great Plains ecosystem

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 imperilled 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.


| INTRODUC TI ON
Grasslands and aridlands, and the birds that rely on them, are under threat. Over 60% of North America's grasslands (i.e., grass-dominated habitats) have been lost since the 1800s, with some of the greatest losses in the Northern Great Plains Tallgrass and Mixedgrass Prairie regions (Comer et al., 2018). Similarly, aridlands (i.e., arid shrub-dominated communities) have experienced habitat loss and degradation due to urban sprawl and energy development, drought, desertification and invasive species (North American Bird Conservation Initiative & U.S. Committee, 2014; Ryan et al., 2008).
Though much of the grassland and aridland loss and degradation occurred in the nineteenth and twentieth centuries, agricultural conversion, urban sprawl and energy development continue today (Samson et al., 2004;Turnbull et al., 2014). For example, 3 million acres of cropland-a less-suitable habitat for grassland birds-were created during 2008-2012, 77% of which were converted from grassland and 8% from shrubland (Lark et al., 2015). Additionally, grassland and aridland birds face other ongoing threats including agricultural intensification and pesticides (Stanton et al., 2018). Consequently, populations of grassland-dependent birds have declined by an estimated 53.3% during 1970-2017, while aridland birds declined 17%-46% during 1968, (North American Bird Conservation Initiative, 2014Rosenberg et al., 2019). These extensive declines highlight the need for effective grassland and aridland bird conservation efforts.
Conservation metrics often rely on evaluations of habitat quality for single species that are treated as indicators (Simberloff, 1998). Yet single species rarely serve as effective proxies for larger communities, due to limited overlap in habitat requirements, functional roles and responses to management actions among species (Cushman et al., 2010;Johnson et al., 2019;Larsen et al., 2009;Lindenmayer et al., 2015;Winter et al., 2005). Single-species indicator approaches are particularly ill-suited to communities with threatened species, as more common species often perform poorly as surrogates for species of conservation concern (Stephens et al., 2019). 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 quantify the response of ecological communities to conservation actions (Kapos et al., 2008(Kapos et al., , 2009Nuttle et al., 2003).
By using multispecies 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 et al., 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 et al., 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 than communities with greater redundancy (Cardinale et al., 2012;McGill et al., 2006). Consequently, functional diversity metrics incorporating functional traits and species abundance provide an indirect way to measure resilience and integrity (Standish et al., 2014), and are increasingly used in large-scale assessments of North American bird and ecological communities (Schipper et al., 2016;Schleuter et al., 2010).
Private lands offer a potential source of hope for grassland and aridland birds, as 84% of grasslands and 44% of aridlands are held in private ownership; therefore, private land conservation efforts could have substantial population-level impacts (Askins et al., 2007;North  Through this programme, Audubon rangeland biologists partner with ranchers to co-develop bird-friendly Habitat Management Plans. In order to evaluate the success of habitat enhancement initiatives such as the Conservation Ranching programme in conserving birds and their habitats on private lands in North American grasslands and aridlands, we developed a metric 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 grassland and aridland bird communities spanning broad climatic and community gradients. 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. Here, we address two objectives 1) to demonstrate the use of the BFI to map grassland bird resilience K E Y W O R D S accountable conservation, aridland birds, bird-friendly, conservation, grassland birds, grazing, indices, management, ranching, working lands across the Northern Great Plains and identify ecologically significant areas for grassland and aridland birds and 2) to evaluate the effects of bird-friendly habitat management practices on grassland and aridland bird communities.

| Bird-Friendliness Index overview
We calculated the BFI using 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 populations (e.g., Rosenberg et al., 2016;Sauer et al., 2017) and form the foundation of the BFI. We also include conservation status to ensure that the BFI does not underestimate conservation need by diluting the contribution of vulnerable 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: First, densities were estimated on a species-specific basis (here: Baird's Sparrow) for surveyed grid cells using distance and time removal information to correct for imperfect detection (1). Estimated densities were combined with remotely sensed environmental covariates in bird-habitat models to predict densities across the Northern Great Plains (NGP) study area (2; here Swainson's Hawk, Baird's Sparrow and Western Meadowlark). Predicted densities across the NGP were multiplied by species-specific breeding season Combined Conservation Scores (CCSb; 3). Predicted densities and species' functional traits (e.g., diet composition) were used to calculate functional diversity for each grid cell (4). Conservation-weighted densities were summed for each grid cell and multiplied that cell's functional diversity (5). Finally, values were scaled from zero to one using a logistic transformation (6) We calculated the index in six steps, all of which are calculated at the resolution of 1 km 2 grid cells ( Figure 1). We chose this resolution because it represents the smallest scale at which the avian survey data used are summarized (see Section 2.4), and is large enough to contain entire home ranges of most study species.
Calculation of BFI scores at this resolution also enables BFI values to be easily aggregated to larger scales (e.g., ranch, state) by calculating the mean or median of grid cells included within the area of interest.

| Study area
Our study area encompasses the Northern Great Plains ecosystem, as defined by the West-Central Semi-arid Prairies level II ecoregion Sand Prairie to the south (Comer et al., 2018). The study area also includes extensive arid sagebrush steppe habitat (i.e., aridlands) intermixed with grasslands (Connelly et al., 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. These species were evaluated for inclusion in the BFI (see Section 2.4). primary habitat type at the point and measured distance using a laser rangefinder. We excluded flyover observations, juvenile birds and observations with missing minutes or distances.

| Avian data
We assigned IMBCR survey points to spatial strata defined as unique combinations of BCRs and states or provinces ( Figure S1a).
We defined strata this way because BCRs delineate ecoregions with similar bird communities and habitats; combining with states/ provinces ensures that each stratum faces similar policy and management regimes. These strata are also used in hierarchical modelling of Audubon Christmas Bird Count data (Soykan et al., 2016) and North American Breeding Bird Survey data (Sauer et al., 2017)

| Environmental data
We used 12 environmental variables in individual species models to project avian densities across the study area (see Section 2.7).
We sampled all environmental variables for 799,015 1 km 2 grid cells, covering the study area and aligned with the avian density grid cells.
The 12 variables comprised five classes of environmental predictors (climate, ecosystem, land cover, temporal and topographical) and were collected or calculated either once (e.g., land cover from 2010) or annually (e.g., ecosystem variables; Table 2). We included three land cover types-grassland, shrubland and cropland-as they represent habitat frequently used by grassland and aridland birds; other potential land cover types (e.g., developed) were too sparse to explain much variation in bird occurrence or abundance. We cal- We selected environmental variables based on known relationships with distribution and abundance of grassland and aridland birds (Table S1). Correlations among predictors were minimal (all r ≤ |.70|).

| Density estimation
As the avian dataset included both time of first detection and distance of observations, we used a conditional multinomial maximum-likelihood formulation of time removal and distance sampling models developed by Sólymos et al. (2013) and implemented using command cmulti in R package detect (Solymos et al., 2016). This method improves upon distance estimation alone, which allows estimation of just one component of detection probability (perceptibility), by utilizing minute of first observation to TA B L E 1 Grassland and aridland bird species included in Bird-Friendliness Index estimation for the NGP, their major breeding habitats, functional species grouping, breeding season Combined Conservation Score (CCSb) and mean density (birds/km 2 ) at surveyed grid cells estimate availability, defined here as the singing rate (Matsuoka et al., 2014). Availability and perceptibility are further defined in the Supporting Information.
We combined data from 2009 to 2014 in a single model per species to maximize the number of species with sufficient data for modelling, but allowed availability and perceptibility to vary among years.
We modelled species with a minimum of 60 detections (Buckland et al., 1993). To estimate point-level density, we entered the correction factors as offsets in generalized linear mixed models calculated using R package lme4 (Bates et al., 2015). We used a Poisson distribution and included primary habitat as a fixed effect and the unique combination of grid cell and year as a random effect. We left-truncated density estimates at 0.08 birds/km 2 , equivalent to the density estimate for a point with a single bird observed at the maximum observed distance in the dataset (2,000 m). We extrapolated annual densities across the entire 1 km 2 area of all surveyed grid cells by taking the median of annual point-level densities (birds/ha) within each surveyed grid cell multiplied by 100.

| Species habitat modelling
We built species habitat models, with density as the response variable, for each species using boosted regression trees (BRTs), and used these to estimate densities of our 22 focal species in all 799,015 grid cells across the study region. BRTs are a machine-learning approach ideal for modelling complex species-environment relationships with TA B L E 2 Environmental variables used as predictors in the species distribution models, with their sources, frequency of collection or calculation and citations Annual CASA ecosystem model (Potter et al., 1993(Potter et al., , 2007 Ecosystem Net primary productivity (g C/ m 2 ) Annual CASA ecosystem model (Potter et al., 1993(Potter et al., , 2007 Ecosystem Nitrous oxide flux (g N 2 O/m 2 -day) Annual CASA ecosystem model (Potter et al., 1993(Potter et al., , 2007 Ecosystem Soil moisture (cm; 0-10 cm depth) Annual CASA ecosystem model (Potter et al., 1993(Potter et al., , 2007 Land Derived from digital elevation model (Riley et al., 1999) multiple predictors and are robust to correlations among predictors (Elith et al., 2008). Many species had skewed density distributions with absences in many grid cells (i.e., zero inflation), which violates the Poisson model assumption that the mean equals the variance.
Therefore, we implemented a hurdle model approach in which we separately modelled occurrence and density, then combined the models to estimate density only at grid cells that met a threshold occurrence level. For the occurrence model, we used presence/absence as the response variable and for the density model detectioncorrected abundance scaled to 1 km 2 . 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.
We used the 12 environmental variables described above plus year as predictor variables (

| Conservation scores
We incorporated the breeding season Combined Conservation Score

| 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 (Varela et al., 2014) to group species into functional groups (Table 1) based on dietary, foraging and size-based traits (Table 3). We then calculated Shannon's diversity using functional group identities in place of species identity. The methods for calculating functional diversity are further detailed in the Supporting Information.

| BFI calculation and evaluation
We calculated BFI values by multiplying the summed conservationweighted densities by functional diversity for each 1-km 2 grid cell.
We chose to multiply rather than add densities and functional diversity, as the summed conservation-weighted densities are orders of magnitude larger than functional diversity indices and vary widely among grid cells. As a result, functional diversity information would be swamped if simply added to densities. Alternatively, raising densities to the power of functional diversity would greatly increase the range and skewness of BFI scores.
Finally, we standardized raw BFI values by annually scaling from zero to one using a logistic distribution, to accommodate the lower limit of zero and the occasional, exceptionally high BFI value.
Standardization expresses each cell's bird-friendliness relative to the study area, producing a ranked index scaled from zero to one that is easily interpreted and compared across an extensive, climatically and topographically variable range. It also controls for large-scale factors such as climatic fluctuations or population-level density-dependent processes that influence population trends, thus isolating the effects of local management actions. For study area mapping, we standardized using all grid cells to enable comparison across the NGP as a whole. However, for the assessment of bird community response to simulated habitat management, we standardized using only grid cell scores from the surrounding spatial strata (e.g., South Dakota-Badlands and Prairies; Figure S1a) to highlight effects of management relative to nearby areas with similar climate and land cover.
We quantitatively evaluated the BFI in two ways. First, we quantified spatial variation by analysing differences in mean BFI values among strata across the study period using a generalized linear mixed model with strata as a fixed effect and year as a random effect. We explored the fit of four distributions (gamma, beta, lognormal and Gaussian) using R package fitdistrplus (Delignette-Muller & Dutang, 2015). As the beta distribution best approximated the distribution of BFI values, we fit the model in R package glm-mTMB (Brooks et al., 2017) and conducted Tukey's post hoc evaluations of differences among strata using emmeans (Lenth, 2020).
Second, we evaluated the BFI's ability to identify ecologically significant areas for grassland and aridland birds by comparing BFI values among areas previously identified as priorities using other criteria.
We used the Consensus Priorities identified by Grand

| Sensitivity analyses
We conducted a sensitivity analysis to assess the relative influence of species density, conservation score and functional diversity on BFI values. We used bootstrapping (1000 iterations) to evaluate the influence of increasing the variability (here, standard deviation) of each component of the BFI on the index value. Details of the sensitivity analysis methods are provided in the Supporting Information.

| Habitat management case study
We simulated an evaluation of the effects of management actions on a theoretical private property, for example an Audubon Conservation Ranch. In these simulations, we modelled the effects of common bird-friendly habitat management actions including reducing cropland cover and N 2 O and increasing litter biomass, net primary productivity, soil moisture, proportion grassland, grassland cohesion and grassland patch area. Our selection criteria for using these covariates are explained in the Supporting Information.
Climate and land cover variables explained the most variation in occurrence and density of the grassland and aridland birds studied here. Climatic moisture deficit explained the most variation in occurrence (mean ± SE: 17.10 ± 3.54%), followed by proportion cropland (mean ± SE: 13.60 ± 3.28%) and proportion grassland (mean ± SE: 12.51 ± 2.47%; Figure 2a). Proportion cropland cover explained the most variation in density (mean ± SE: 17.58 ± 4.92%), followed by climatic moisture deficit (mean ± SE: 13.42 ± 3.24%) and proportion grassland (mean ± SE: 9.18 ± 1.94%; 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 space and, in some cases, time explained more variation than year alone. Relative variable importance varied among species consistent with species-specific habitat preferences; for example, Rock Wren density was explained most by terrain ruggedness, while occurrence and density of grassland birds like Sprague's Pipit and Western Meadowlark were explained most by grassland proportion cover and patch area ( Figures S2-S23). Drivers of grassland and aridland bird occurrence and abundance are further discussed in the Supporting Information.
BFI values were higher in areas identified by one or more previous grassland conservation prioritizations and increased with the number of prioritizations the area was included within (Figure 4b).
F I G U R E 2 Mean variable importance and 95% confidence intervals for 12 variables used as predictors in presence/absence (a) and density (b) species habitat models for 22 grassland and aridland bird species across the Northern Great Plains during 2009-2014. Climate and land cover variables explained the most variation in occurrence and density

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

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

| D ISCUSS I ON
The ability to quantify the impacts and success of conservation and management actions is crucial to the adaptive management cycle (Nichols & Williams, 2006). Here, we show that the Bird-Friendliness Index (BFI) served as a proxy for identifying ecologically significant areas for grassland and aridland birds, with mean BFI values increasing with the degree of consensus across published grassland conservation prioritizations. Moreover, the BFI detected regional variation in community richness and functional diversity, further supporting its utility as a metric for identifying ecologically significant areas for grassland and aridland birds. The BFI also effectively detected South Dakota had a more balanced and diverse representation of functional traits but lower densities in its grassland and aridland bird community, whereas Alberta had higher densities overall-or higher densities of species of conservation concern-but lower FD.
Conversely, the lowest BFI values were found in the south-western regions of the study area, notably southern Montana, Wyoming and southern South Dakota. 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 Brewer's Sparrow, they are less suitable for grassland birds that prefer lusher and more productive habitats (Fedy et al., 2018;Fisher & Davis, 2010;Harrower et al., 2017;Renfrew et al., 2013).
As a result, they have fewer bird species, many of which share similar functional traits (e.g., ground-dwelling insect-and granivores; Parks . There are many opportunities for further refinement and application of the BFI. Additional bird survey data, through inclusion of additional years and/or an expanded survey area, would enable estimation of robust density estimates for additional grassland and aridland bird species, enabling the BFI to represent a broader slice of the community. Further, we were unable to incorporate annual land cover in this case study due to the lack of robust annual land cover classifications that span the US/Canada border. Future applications of the BFI should take advantage of annual land cover data, where available, because of these inherent sensitivities. We designed the BFI for the purpose of conducting 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. Our standardization approach enabled us to isolate local-scale effects, such as bird response to habitat management or regional differences in land use, from largescale 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). 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.
The BFI serves as a tool for quantifying grassland and aridland bird community response to management, which enables its use to inform a robust adaptive management process (Lancia et al., 1996).
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 providing conditions that enhance ecosystem processes (Rodrigues, 2006). Under this framework, scientists and land managers would work together to refine habitat management protocols to increase abundance and functional diversity-and consequently resilience-of grassland and aridland bird communities on privately managed lands. To be able to do this effectively requires indicators that are rigorous, repeatable and easily understood (Balmford et al., 2005).
Recent grassland and aridland bird population declines call for implementation of new conservation and restoration efforts that demonstrably improve habitat conditions and slow or reverse declines. 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. In short, the BFI is a tool that can be used by conservationists and managers to develop and measure accountable conservation now and into the future.

ACK N OWLED G EM ENTS
We thank Max Alleger (Missouri Department of Conservation) and

CO N FLI C T O F I NTE R E S T
None.

PE E R R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/ddi.13163.

DATA AVA I L A B I L I T Y S TAT E M E N T
The R scripts used to conduct the analyses and simulations are avail-