Multiscale occupancy of the Lesser Prairie-Chicken: the role of private lands in conservation of an imperiled bird
 
Présence multiéchelle du Tétras pâle : rôle des terres privées dans la conservation d'un oiseau en péril

• Abstract
• Introduction
• Methods
  • Study area
  • Covariate development
  • Model justification and hypotheses
  • Sampling design and field surveys
  • Statistical analysis
• Results
  • Annual variation in site occupancy
  • Multiscale covariate relationships
  • Exploratory model selection
  • Mapping the range-wide occupancy distribution
  • Simulated effect of changes to CRP-enrolled land
• Discussion
  • Annual variation in site occupancy
  • Multiscale covariate relationships
  • Exploratory model selection
  • Simulated effect of changes to CRP-enrolled land
• Responses to this article
• Acknowledgments
• Literature Cited

INTRODUCTION

Avian species are facing an unprecedented rate of decline globally. In particular, grassland avifauna of North America and Europe have experienced more dramatic declines than forest dwelling species (BLI 2018). The causal factors for these changes in abundance and distribution have varied but largely resulted from conversion of grassland to other uses, e.g., intensive agriculture or energy development (BLI 2018). In North America, grassland-obligate birds are sensitive to the loss and fragmentation of native prairie (Brennan and Kuvlesky 2005), and tend to select the largest patches of grassland available (Ribic et al. 2009). Thus, conservation efforts have focused primarily on retention of native grasslands, and second on restoration (or reclamation) to improve connectivity and quantity of grasslands. Most remnant grasslands are privately owned and used to produce forage for livestock (Samson et al. 2004).

Historically, primary causes of disturbance in prairies were periodic fire and herbivory from large herds of nomadic grazers, which promulgated spatial heterogeneity in vegetative cover and composition. However, under current land uses, most prairie remnants are infrequently (or more likely never) burned and repeatedly experience uniform grazing, which facilitates vegetative homogeneity across the prairies (Samson et al. 2004, Fuhlendorf et al. 2006). Loss of heterogeneity is spatially hierarchical; at the broadest scales agriculture conversion of prairie produces monocultures of crops, and redundant grazing and lack of fire in uncropped land can lead to woody encroachment or monotypic stand of herbaceous cover (Fuhlendorf et al. 2006). Evidence is mounting that these landscape changes are detrimental to grassland-obligate species and ecosystem function (Engle et al. 2008, Hovick et al. 2015, Pavlacky et al. 2017).

Species presence at fine, localized spatial scales, e.g., an individual home range or habitat patch, is conditional upon life-history needs at broader landscape scales (Johnson 1980). These concepts of hierarchal habitat selection are not new (Johnson 1980, Cody 1985), but analytical approaches are more prevalent now to rigorously evaluate species occurrence and habitat use across multiple scales. Such investigations are critical as anthropogenic disturbances continue to remove and fragment available prairie landscapes. Multiscale occupancy models enable us to explicitly examine the hierarchy of habitat selection across distributional limits while simultaneously accounting for imperfect detection (Nichols et al. 2008, Pavlacky et al. 2012). In the case of prairies, we can begin to characterize both broader landscapes that are suitable for occupancy by avian species (first-order selection; Johnson 1980), and conditions therein that are suitable for occupancy at finer spatial scales (second-order selection; Johnson 1980). These models allow us to test hypotheses for habitat use conditional on hierarchical landscape structure for obligate grassland species.

Reversing loss, fragmentation, and degradation of prairie ecosystems may require coordinating bird conservation with large-scale programs to stabilize upland bird populations associated with grasslands (Brennan and Kuvlesky 2005). One such approach is to use a flagship species to focus restoration or conservation efforts (Caro 2010, Miller et al. 2017). In the Southern Great Plains, the Lesser Prairie-Chicken (LEPC, Tympanuchus pallidicinctus) is an iconic species that requires habitat heterogeneity at multiple scales to meet life history needs (Hagen and Elmore 2016). LEPC is a species of conservation concern because it only occupies approximately 16% of its historical distribution and is currently under consideration for protections under the Endangered Species Act. Causes behind loss, fragmentation, and degradation of native prairie are threats to long-term persistence of LEPC (USFWS 2016). Thus, understanding biotic and abiotic factors driving LEPC occupancy across the Southern Great Plains may assist in delivering effective conservation and restoration of prairie ecosystems beyond LEPC.

Over 95% of land in LEPC distribution is privately owned (Van Pelt et al. 2013), used mostly for livestock grazing and agricultural cropping. These patterns of ownership and land use may require unique solutions to restore habitat and conserve imperiled species (Briske et al. 2017). Technical and financial resources through the U.S. Department of Agriculture (USDA) can restore prairie (e.g., Conservation Reserve Program [CRP]) and improve habitat quality and ecosystem services (e.g., Environmental Quality Improvement Program [EQIP]) on private lands.

CRP provides rental payments to landowners to convert marginal farmland back to grasslands for a contractual period of 10 years. EQIP offers technical and/or financial assistance to landowners to develop conservation plans for their land to improve their resiliency to drought and increase forage production for livestock. Funding and prioritization of these conservation actions are guided by the “Farm Bill,” i.e., Agricultural Improvement Act of 2018, which establishes payment schedules, enrollment goals, and quotas. Although not designed to deliver species conservation per se, direct and indirect benefits from the land conservation practices therein have been demonstrated to benefit prairie-obligate species like LEPC. In particular, CRP can increase LEPC food availability, nesting cover, and population resilience to drought (Fields et al. 2006, Ross et al. 2016a, Sullins et al. 2018). Alternatively, EQIP investments have been spatially targeted through USDA’s Lesser Prairie-Chicken Initiative (LPCI) with the goal of maximizing conservation benefit to LEPC. There are 27 LPCI conservation practices to improve or maintain rangeland health, but prescribed grazing and brush management, i.e., removal of invasive trees, are most important regionally (Bartuszevige and Daniels 2013). Additionally, as part of a range-wide LEPC monitoring and habitat improvement program administered by The Western Association of Fish and Wildlife Agencies (WAFWA), financial assistance was provided to landowners to develop conservation easements as offsite mitigation for energy developments (Van Pelt et al. 2013). Here, we developed predictive models of LEPC multiscale occupancy to investigate covariate relationships over time, expanding upon previous efforts that were temporally limited (Hagen et al. 2016). Predictive multiscale models were used to investigate second-order habitat relationships (Johnson 1980) using the theory of hierarchical habitat use (Cody 1985), where habitat use at the small scale (56.25 km²) is conditional on habitat use at the large scale (225 km²). Relationships between occupancy and covariates of interest have implications for landscape conservation at multiple scales (George and Zack 2001), perhaps suggesting management actions that could maintain or increase the range-wide extent of occurrence of LEPC.

We included covariates that described habitat composition and configuration, anthropogenic development, drought-related climatic conditions, and conservation efforts; and we examined the effects of covariates from 2012 to 2016 (Table 1). Our objectives were to (1) quantify annual variation in the probability of LEPC occupancy at two spatial scales over five years of study throughout the entire species range, (2) identify the most important predictors of LEPC occupancy range-wide and within each of four ecoregions, (3) map LEPC occupancy probability range-wide as a function of the most important predictor variables, and (4) project how changes in restoration programs, i.e., CRP, may affect species distribution. We predicted fragmentation effects, i.e., mean patch size, of grassland, shrubland, or native habitat would be more important than amounts of habitat, i.e., land cover, at the large scale (Hagen et al. 2016). We expected occupancy to decline with anthropogenic development, cumulatively or by individual type of development, at the small scale (Hagen et al. 2011). Because at certain life stages LEPC are sensitive to drought (Grisham et al. 2016), we predicted that range expansion or contraction at the small scale would be correlated with spatial and temporal variation in the duration of drought during spring and summer months. We predicted that CRP and woodland landcover would influence site occupancy at the large scale through processes of habitat loss and fragmentation, and we expected prescribed grazing and conservation easements would influence site occupancy at the small scale by improving habitat condition. Finally, because ecoregions differ with respect to biophysical and anthropogenic processes, we predicted the above responses would vary by ecoregion (Hagen and Elmore 2016).

METHODS

Study area

The study area spanned the entire estimated occupied range of LEPC in 2011 (~80,000 km²), including portions of five states: Colorado, Kansas, New Mexico, Oklahoma, and Texas (McDonald et al. 2014; Fig. A1.1). Because of expected geographic variation in LEPC habitats, distribution, and abundance, we stratified the study area into four ecoregions for ecoregion-level analyses: Shinnery Oak Prairie (SOPR), Sand Sagebrush Prairie (SSPR), Mixed Grass Prairie (MGPR), and Shortgrass/CRP Mosaic (SGPR; McDonald et al. 2014; Fig. A1.1).

Covariate development

We derived covariates that described anthropogenic land uses, drought-related climatic conditions, conservation actions, and vegetative landcover at two spatial scales (225 km² grid cells and 56.25 km² quadrants) within a Geographic Information System (Table 1). As possible, the value of a grid- or quadrant-level covariate could change annually as updated source datasets were available. However, covariates representing primary road density, transmission line density, and landcover types sourced from the National Land Cover Database (Tables 1, A1.1) were assumed to be constant through time.

Model justification and hypotheses

To identify the most important predictors of LEPC occupancy, we used predictive models and the method of multiple working hypotheses (Chamberlin 1965); we used model-based tests of hypotheses for effects of landscape structure, anthropogenic development, conservation practices, and drought-related climatic conditions on site occupancy. We did this using covariate relationships at two spatial scales (56.25 km² quadrants and 225 km² grid cells) for which LEPC may respond differently (Fuhlendorf et al. 2002). The multiscale occupancy model provides inference to relationships between occupancy patterns and covariates of interest at these spatial scales. Animals select habitat hierarchically (Hutto 1985), and understanding occupancy patterns at multiple spatial scales is imperative for successful management of wildlife and their habitats (Chalfoun and Martin 2007).

The primary management question for landscape structure involved addressing the importance of habitat loss and fragmentation on range contraction and expansion. We drew inferences about habitat loss and fragmentation from patterns of landscape composition (landcover type) and configuration (patch size), respectively.

We predicted that landscape configuration and mean patch size of grassland, shrubland, or native habitat would be important for site occupancy of LEPC (Hagen et al. 2016). However, we were uncertain whether lands enrolled in CRP would contribute to core habitat patches or function as landscape context of between-patch matrix habitat. We hypothesized that LEPC would respond negatively to increases in landcover and patch size of cropland (Haukos and Zavaletta 2016). We also considered an alternate hypothesis that LEPC would respond positively to landscape heterogeneity (Fahrig et al. 2011), wherein probabilities of occupancy would be highest at intermediate values of cropland landcover or patch size (Ross et al. 2016a). We also investigated quadratic responses for grassland, shrubland, and native habitat to represent hypotheses for landscape heterogeneity involving nonlinear responses to suitable habitat at the landscape scale. Finally, we investigated possible interactions between ecoregion (as a factor) and continuous landscape composition and configuration covariates because we hypothesized that habitat-occupancy relationships likely varied by ecoregion.

We developed hypotheses for anthropogenic disturbance using covariates for vertical structures, oil and gas wells, primary road density, transmission lines, and landcover associated with anthropogenic development (Table 1). We predicted LEPC occupancy would decline with increasing anthropogenic development cumulatively and by individual type of development (Bartuszevige and Daniels 2016).

We tested hypotheses for conservation efforts in the ecoregions using covariates for CRP-enrolled land, LPCI-prescribed grazing, and WAFWA conservation easements. We predicted that LEPC occupancy would increase with increasing landcover of LPCI core conservation practices, including prescribed grazing and CRP (USFWS 2011, Hagen et al. 2016). We tested whether LEPC responded to CRP in relation to (1) its contribution to land cover and configuration of general habitat, (2) its additive area as between-patch matrix habitat, and (3) the additive effect of its mean patch size. We evaluated quadratic relationships for percent landcover and patch size of CRP to investigate whether occupancy was highest at intermediate amounts of these covariates.

Climate change in the Southern Plains is expected to influence the population viability of LEPC. Interactions between spring precipitation and vegetation cover may potentially influence key population vital rates, such as nest survival and recruitment (Grisham et al. 2016). We predicted that LEPC range expansion would be correlated with spatiotemporal variation in the duration of drought conditions during spring and summer (Table 1). Extreme temperatures during summer drought can negatively influence abundance of invertebrate prey, with consequences for LEPC recruitment and survival (Grisham et al. 2016).

Sampling design and field surveys

The sampling design and field methodology were detailed by McDonald et al. (2014), and we offer only a brief summary here. As part of a range-wide LEPC monitoring program administered by WAFWA, a spatially balanced sampling procedure was used to select 15 km × 15 km grid cells to survey for LEPC. Survey efforts varied annually and by ecoregion, but approximately 250–300 total grid cells were surveyed each year, with each grid cell surveyed once per year (McDonald et al. 2014). We subdivided each grid cell into four quadrants (7.5 km × 7.5 km each; Hagen et al. 2016). During 2012–2016, surveyors flew a single 7.5-km line transect through each quadrant during March, April, or May and recorded detections of prairie-chickens within 300 m of the line using a double-observer method and four observers in the helicopter. LEPCs, Greater Prairie-Chickens (T. cupido), and their hybrids co-occur in portions of the SGPR, but are not reliably distinguishable during aerial surveys. Therefore, on-the-ground visits were conducted to verify species identification in areas where mixed-species groups were possible (McDonald et al. 2014, Hagen et al. 2016).

Statistical analysis

Sampling framework for multiscale occupancy

We aggregated and summarized LEPC detection data such that large-scale occupancy (ψ) corresponded to LEPC presence on 15-km × 15-km grid cells, small-scale occupancy (θ) corresponded to species presence in four quarters of the grid cell (7.5-km × 7.5-km quadrants) given presence within grid cells, and detection probability (p) corresponded to detections by observers given presence in quadrants. The encounter history was arranged by treating independent observers in the helicopter as independent sampling occasions to estimate the probability of detection (Hagen et al. 2016). We pooled encounters of LEPC across observer in the front-left seat and pilot in the front-right seat (first occasion or search). Similarly, we pooled encounters across observers in the back-left seat and back-right seat (second occasion or search). This yielded an encounter history for each quadrant with two occasions (or searches) each time a survey flight was conducted. The multiscale model can be interpreted as having two components of the observation process; the probability of detection by observers (p) given presence in the quadrants and the probability of availability within quadrants (θ) given presence in the grid cells (Hagen et al. 2016, Pavlacky et al. 2012). To account for heterogeneity in detection probability from multiple observers on a single day (MacKenzie et al. 2018), we considered covariates for ordinal date and time since sunrise to account for within-season variation due to timing of surveys. Although, we did not repeatedly sample line transects within seasons, estimates of small-scale occupancy can be interpreted as the probability of spatial availability, which accounts for the probability that LEPC may be unavailable for sampling at line transects across days and years of the study (Pavlacky et al. 2012, Hagen et al. 2016).

Implicit dynamics multiscale occupancy

We estimated detection and occupancy probabilities of LEPC using implicit dynamics (MacKenzie et al. 2018) version of multiscale occupancy model (Nichols et al. 2008, Pavlacky et al. 2012). The multiscale occupancy model provides inference to relationships between occupancy patterns and covariates of interest at two spatial scales. Parameters of the model were (1) detection probability p_ijot for observer o, quadrant j, grid cell i and year t given the quadrant and grid cell were occupied; (2) small-scale occupancy θ_ijt for quadrant j, grid cell i and year t given the grid cell was occupied; and (3) large-scale occupancy ψ_it for grid cell i and year t.

We ranked the candidate set of models using Akaike’s Information Criterion adjusted for sample size (AIC_c; Burnham and Anderson 2002), with sample size defined by the number of 15-km × 15-km grid cells surveyed across years. We evaluated multimodel support for hypotheses using evidence ratios and cumulative AIC_c weights for balanced model sets ([w₊(j)]; Burnham and Anderson 2002). We determined variable support for unbalanced model sets according to γ_i = [w_i/(1-w_i)]/[f_i/(1-f_i)], where w_i is the cumulative AIC_c weight and f_i is the frequency of the covariate i in the model set (Doherty et al. 2012). We evaluated effect sizes and conditional 90% confidence intervals (CIs) for β parameters with respect to 0. We model-averaged estimates of large-scale or small-scale occupancy for models with ΔAIC_c < 4 (Burnham and Anderson 2002).

Annual variation in site occupancy

We investigated annual variation in large-scale (ψ) and small-scale (θ) occupancy of LEPC to evaluate temporal patterns within four ecoregions. Candidate sets for ψ and θ were each composed of five models, including the full model (ecoregion + year + ecoregion * year) and reduced models (ecoregion + year), (ecoregion), (year), and intercept (.) only. We modeled detection (p) according to three continuous covariates for ordinal date, time after sunrise, and annual trend, and three factor covariates for ecoregion, observer, and year (Table 1). We excluded detection models containing both a continuous covariate for annual trend and factor covariate for year. The candidate model set for detection included all subsets of five covariates and intercept only model p(.), for a total of 61 models. We fit all subsets of covariates and parameters (Doherty et al. 2012) for a total of 1200 models using the RMark interface (Version 2.2.4, Laake 2013, R Core Team 2017) for program MARK (Version 8; White and Burnham 1999). Although Burnham and Anderson (2002) cautioned against all subsets model selection, applying model-averaging to all subsets of models may reduce bias in parameters and model weights (Doherty et al. 2012). It is recommended to consider fewer models than the number of samples, as in our analysis (n = 1404), and this is expected to reduce the risk of spurious results (Burnham and Anderson 2002). We estimated rates of change in occupancy (λ_t) between the t temporal oscillations in small-scale occupancy (θ_t) according to λ_t = θ_t+1/ θ_t (MacKenzie et al. 2018). We approximated the variance of λ_t using the delta method (Powell 2007) and calculated 90% log-normal CIs for λ_t (Burnham et al. 1987).

Range-wide relationships between covariates and multiscale occupancy.

We used all range-wide data in the four ecoregions from 2012 to 2016 (McDonald et al. 2016), but did not use auxiliary data collected within SGPR and SOPR during 2015 (Hagen et al. 2016). Prior to model selection, we used a variable screening step to identify potential curvilinear quadratic relationships for continuous covariates, and two-way interactions between covariates for each parameter. We fit univariate and quadratic models for each covariate, and to evaluate two-way interactions, we fit additive and multiplicative models for covariates. We selected quadratic relationships for entry into the analysis when AIC_c was lower than a univariate relationship. In a similar fashion, we selected a two-way interaction for entry into the analysis when AIC_c was lower for the multiplicative model than the additive model. For (ψ, we fit quadratic and interaction models while holding constant θ(Ecoregion) and p(Observer + Ecoregion + Year). For (θ), we fit quadratic and interaction models while holding constant ψ(Year) and p(Observer + Ecoregion + Year). For (p), we fit quadratic and interaction models while holding constant ψ(Year) and at θ(Ecoregion).

Model selection procedures increase in complexity when models contain multiple submodels and when research objectives require modeling the effect of multiple, potential predictor variables (Bromaghin et al. 2013). The multiscale occupancy model we used was composed of three separate submodels: large-scale occupancy (ψ), small-scale occupancy (θ), and detection probability (p; Nichols et al. 2008, Pavlacky et al. 2012, Hagen et al. 2016). Moreover, our objectives necessitated modeling ψ, θ, and p as functions of multiple predictor variables. We therefore adopted a two-staged model-selection approach to first select plausible structures for each submodel, then to consider all possible combinations of plausible submodel structures (Table 2).

The plausible-combinations approach proceeded in two steps. First, we identified plausible covariate relationships for each parameter independently, and second we combined all-model subsets of the submodels across parameters to identify parsimonious full-models (Bromaghin et al. 2013). For each parameter, we selected high-weight submodels with AIC_c weight w_i > 0.01 and high-likelihood submodels with -2log(L) < maximum [-2log(L) of high-weight models] for entry into the second step of the plausible-combinations model selection analysis (Bromaghin et al. 2013).

Prior to plausible combinations model selection, we constrained the candidate set of models by omitting submodels with correlated covariates (Pearson’s ρ > 0.6). We flagged submodels with diminutive (< 0.00001) standard errors (SE), and submodels with small (< 0.5) or large (> 5) t-ratios (β/SE) for inspection. Additionally, we constrained the candidate set of models by omitting submodels with redundant covariates, defined as those representing similar biological hypotheses with nonexhaustive and nonexclusive classification, e.g., grassland landcover and native habitat were not allowed in the same model. Finally, we expanded the candidate set of models by appending submodels that replaced main effects by supported quadratic relationships for each covariate, and two-way interactions.

First, we ran all subsets of 29 covariates for ψ with a maximum of three covariates per models while holding constant θ(Ecoregion) and p(Observer + Ecoregion + Year), resulting in a candidate set of 8249 plausible models. For θ, we ran all subsets of 19 covariates with a maximum of four covariates per model while holding constant ψ(Year) and p(Observer + Ecoregion + Year), resulting in a candidate set of 16,218 models. For p, we ran all subsets of six covariates with a maximum of four covariates per model while holding constant ψ(Year) and θ(Ecoregion), resulting in a candidate set of 46 models. Although, considering a large number of candidate models may increase the risk of spurious results (Burnham and Anderson 2002), the plausible combinations approach requires orders of magnitude fewer models than the all-subsets approach while providing a more comprehensive search of parameter space than other sequential model building approaches (Bromaghin et al. 2013). Although it may be questionable to consider sets of thousands of candidate models due to the risk of spurious results, model complexity and the goal to predict occupancy at multiple scales in our study required this approach.

Second, we combined all subsets of the high-weight and high-likelihood submodels across parameters (Bromaghin et al. 2013), for a total of 40 plausible model combinations. We considered additional models in an exploratory a posteriori to check whether data supported models with greater complexity than those with conditions imposed in the a priori analyses. We added each nonredundant candidate covariates one at a time to the top a priori selected models, and we evaluated all subsets of amended submodels, for a total of 8526 models. In addition, we ran the highest ranking covariate model from Hagen et al. (2016) to investigate previously identified effect size for LPCI prescribed grazing [p(Ecoregion + observer) θ(Ecoregion + CRP) ψ(NativeHabitatPatch + Grazing)]. We evaluated effect sizes and conditional 90% CI for the covariate β coefficients with respect to 0, including odds ratio interpretations for the β coefficients (MacKenzie et al. 2018). We model-averaged estimates of large-scale or small-scale occupancy for all models within ΔAIC_c < 4 (Burnham and Anderson 2002).

Mapping range-wide occupancy

We model-averaged predictions of ψ and θ according to covariate values in the sampling frame. We multiplied conditional estimates of _ij for each of the j quadrants in grid cell i by the corresponding estimate i _i to arrive at unconditional estimates of small-scale occupancy (_ij * _i) for all quadrants in the sampling frame (Nichols et al. 2008, Pavlacky et al. 2012). We approximated SEs for model-averaged unconditional estimates of small-scale occupancy using the delta method (Powell 2007). We estimated the coefficient of variation (CV) for unconditional estimates of small-scale occupancy to quantify uncertainty around the predicted occupancy distribution. We mapped predicted occupancy for 477 grid cells that intersected the occupied range of LEPC plus a 16 km buffer (Van Pelt et al. 2013). Because some covariates change over time, we used covariate values for 2016 to map predicted occupancy presented herein.

Simulated effect of changes to amount of CRP-enrolled land

Modeling range-wide relationships between covariates and multiscale occupancy revealed that CRP was a primary driver of LEPC occupancy. Thus, we conducted a simulation to demonstrate the changes in predicted occupancy as a function of fluctuations in CRP. We used the final set of 40 models to the range-wide dataset to estimate model-averaged occupancy rates assuming nine scenarios for CRP enrollment, the baseline scenario and eight alternatives. The baseline scenario was defined as the amount of CRP enrollment observed in 2016. Alternative scenarios included 80–120% of this baseline amount of CRP-enrolled lands in increments of 5%. For example, to simulate the amount and distribution of CRP enrollment under an 80%-of-baseline scenario, we multiplied the baseline amount of CRP in each quadrant and grid cell by 0.80. All other covariates were held constant. We then generated a model-averaged prediction of the unconditional probability of small-scale occupancy (_ij * _i) for each of the 1908 quadrants in the 477 grid cells that intersected the occupied range of LEPC plus a 16 km buffer (Van Pelt et al. 2013). The area occupied by LEPC was estimated as the mean value of
(_ij * _i), multiplied by the area of the ecoregion. We summarized the effect that changes in CRP levels had on occupancy by calculating the predicted change in total area occupied for each scenario relative to the estimated baseline area. We used the delta method (Powell 2007) to calculate 90% CI for change in area occupied based on quadrant-specific variances of unconditional small-scale occupancy estimated for each quadrant. For simplicity, we represent unconditional small-scale occupancy (_ij * _i) as simply θ_j in the equations below, N_k is the total amount of quadrants within ecoregion, A_k is the area of the ecoregion (km²), and quadrants are again indexed by j. The average change in unconditional small-scale occupancy from the baseline to an alternate CRP scenario
(_k) was calculated as

(1)

with an estimated variance of

(2)

The change in area occupied (ΔAO_k) was calculated as

(3)

with an estimated variance of

(4)

RESULTS

Annual variation in site occupancy

We found almost no evidence for annual variation in ψ and only limited evidence in θ (Table A1.1, Fig. 1). However, range dynamics oscillated temporally in θ (Fig. 1). Within the SGPR, LEPC experienced a 27% range contraction between 2012 and 2013 ( = 0.73; 90% CI = 0.53, 1.01), 15% range contraction between 2013 and 2014 ( = 1.15; 90% CI = 0.83, 1.60), 29% range expansion between 2014 and 2015 ( = 1.29; 90% CI = 0.97, 1.72) and 19% range contraction between 2015 and 2016 ( = 0.81; 90% CI = 0.62, 1.07; Fig. 1). Considering the change between low and high occupancy years in the SGPR, LEPC range expanded by 49% between 2013 and 2015 ( = 1.49; 90% CI = 1.09, 2.02; Fig. 1). The CI for between 2012 and 2013 narrowly covered 1, and CIs between 2013 and 2015 excluded 1, indicating a measurable range contraction and expansion, respectively over these time frames. Magnitude of temporal patterns of change in range was similar in regions (Fig. 1), because spatial and temporal effects of occupancy were additive. In addition, θ was greater in the SGPR than in the other ecoregions across years (Table 3). Detection was greater for back-seat observers than front-seat observers, and p increased over survey years, with increasing time after sunrise, and ordinal date (Table A1.2).

Multiscale covariate relationships

The first step of plausible-combinations model selection identified two plausible submodels for large-scale occupancy, ψ(CRP + GrassPatch + Shrub) and ψ(CRP + CRP² + GrassPatch + Shrub). Because adding CRP² did not appreciably decrease the -2log(L), it was not considered a competing model (Arnold 2010). Thus, we considered only ψ(CRP + GrassPatch + Shrub) as the single plausible submodel for ψ in the second step of plausible-combination models selection (Table A1.3).

We identified three plausible submodels for θ (Table A1.4). However, we did not consider the third-ranked model (ΔAIC_c = 6.00), containing Grass² as a competing model, because it did not appreciably decrease the -2log (L) relative to the second ranked model (Arnold 2010). Thus, we considered the top two models in the second step of the plausible-combinations. We identified 20 plausible submodels for p (Table A1.5), and we considered these models in the second step of the plausible-combinations.

In the second step, we ran all subsets (n = 40) of plausible submodels across parameters. Overall, we found considerable model selection uncertainty and 7 candidate models with ΔAIC_c < 2 (Table 4). The top-ranked model for multi-scale occupancy contained effects of Shrub, GrassPatch size, and CRP on ψ; CRP, Grass, Shrub, Ecoregion, CRP*Ecoregion, and Grass*Ecoregion on θ; and effects of Observer, Trend, and Time on p (Table 4). The second ranked model included the additional effect of Ecoregion on p, but was otherwise equivocal to the highest ranking model (Table 4).

Large-scale occupancy increased with increasing Shrub, GrassPatch, and amount of CRP (90% CIs ≠ 0; Fig. 2, Table 5). Small-scale occupancy varied ecoregionally showing large positive effects with increasing amounts of CRP in the SGPR, smaller positive effects in the SSPR and MGPR, and a much smaller positive effect of CRP in the SOPR (Fig. 3; 90% CIs ≠ 0 for interaction terms).

Small-scale occupancy increased with increasing Shrub, and the slope of the positive effect of Shrub was identical in all ecoregions (Fig. 4). However, the effect was most pronounced in the SOPR because site occupancy and the landcover of shrubland was greater in this ecoregion (Fig. 4). Small-scale occupancy showed a large increase with increasing CRP in the SGPR, followed by much smaller effects in the SSPR, MGPR, and SOPR (Fig. 4, Table 5). The interaction between CRP*Ecoregion indicated the slope of the CRP was much less in the SSPR, MGPR and SOPR than the slope of the CRP effect in the SGPR (90% CIs ≠ 0; Fig. 4, Table 5).

Multimodel inference using cumulative AIC_c weights indicated models with the Grass*Ecoregion interaction showed an 80% probability of occurring in the top model, whereas models including Grass² and interactions with Ecoregion showed a 20% probability of occurring in the top model. In the top-ranked model, small-scale occupancy showed a large linear increase with increasing Grass in the SGPR, a smaller positive effect in the MGPR, and diminutive effects of Grass in the SOPR and SSPR (Fig. 5, Table 5). The interaction between Grass*Ecoregion indicated the slope of Grass was much less in the SSPR, SOPR, and MGPR than in the SGPR (Fig. 5, Table 5). The 90% CIs for interaction terms of Grass in the best model did not include 0, but the CI for the interaction term for Grass and MGPR narrowly included zero in the second best model, suggesting a marginal effect for this Ecoregion (Table 5).

Exploratory model selection

The top-ranked model for the exploratory analysis of LEPC occupancy contained the effects of Woodland-10 on ψ and Development on θ (Table 6). The evidence ratios indicated the effect of Woodland-10 [w₊(j) = 0.89] was 11 times more likely than that Woodland-5 [w₊(j) = 0.08], and 220 times more likely than that with Woodland-1 [w₊(j) < 0.01], and ψdeclined with increasing Woodland-10 (Fig. 6). The cumulative AIC_c weights indicated the effect of Development on θ had an 88% probability of being selected in the top model. The effect of Road [w₊(j) = 0.07] was 3 times more probable than Well [w₊(j) = 0.02], 74 times more probable than Vertical [w₊(j) < 0.01], and 135 times more probable than Transmission [w₊(j) < 0.01], and θ declined with increasing Development (90% CIs ≠ 0, Fig. 7, Table 7). The highest ranking model from Hagen et al. (2016) showed much less support than the highest ranking range-wide model (ΔAIC_c = 168). Beta coefficients for the Grazing (β = 10.5; SE = 5.7; CI = 1.1, 20.0) and NativeHabitatPatch (β = 3.4; SE = 1.4; CI = 1.0, 5.7) covariates demonstrated large and precise effect sizes. Odds ratios indicated occupancy increased by 11% for every 1% increase in the landcover of prescribed grazing (1.11), and 3% for every 100 ha increase in the mean patch size of native habitat (1.03).

Mapping the range-wide occupancy distribution

We predicted range-wide covariate relationships for ψ and θ using unconditional estimates of small-scale occupancy in 2016 (_ij * _i; Fig. 8). The model-based estimate of unconditional small-scale occupancy for the 2016 sampling frame within the occupied range was 0.073 (CI = 0.057, 0.091). Estimates of mean unconditional small-scale occupancy by ecoregion were 0.049 (CI = 0.047, 0.052) for the SOPR, 0.031 (CI = 0.029, 0.033) for the SSPR, 0.073 (CI = 0.071, 0.076) for the MGPR, and 0.123 (CI = 0.121, 0.126) for the SGPR (Fig. 8). The overall model-based estimate of area occupied within the estimated occupied range of LEPC was 7953 km² (CI = 7847, 8059). The area occupied by ecoregion was SOPR = 1370 km² (CI = 1309, 1432), SSPR = 490 km² (CI = 463, 517), MGPR = 2576 km² (CI = 2516, 2636), and SGPR = 3517 km² (CI = 3461, 3574).

Simulated effect of changes to CRP-enrolled land

The amount of land enrolled in CRP varied across the Southern Great Plains, and the SSPR and SOPR had the highest average percentage (Table 8). However, the relationship between CRP and LEPC occupancy was most pronounced in the SGPR (Figs. 3, 9). In this ecoregion, the models predicted that a 5% increase in CRP relative to 2016 (an additional 114 km² enrolled) would result in an additional 228 km² (90% CI = 147 - 310 km²) of LEPC occupied area, i.e., 6.1% increase. Changes in CRP had less-pronounced effects in the remaining ecoregions (Figs. 3, 9). We estimated that the most CRP-liberal scenario (20% more CRP than 2016), would increase LEPC occupied area range-wide by 1581 km² (90% CI = 1418–1744 km²), i.e., a 16.6% increase. The gains in occupied area under this most CRP-liberal scenario were predicted to be within the SGPR (952 km², 90% CI = 864–1039 km²), followed by the MGPR (375 km², 90% CI = 282–468 km²), SSPR (157 km², 90% CI = 110–205 km²), and SOPR (97 km², 90% CI = 7–187 km²) ecoregions (Fig. 9). Increases in CRP associated with this most CRP-liberal scenario was predicted to increase occupancy probabilities range-wide; however, predicted gains in occupancy were most pronounced in portions of west-central Kansas (Fig. 9).

DISCUSSION

Our work provided a multiscaled portrait of biotic and abiotic factors that explained occupancy rates of LEPC across their entire distribution. It revealed conditions necessary for a landscape to be suitable for occupancy and conditions therein suitable for occupancy at local scales. Our results were surprising in the sense that occupancy was insensitive to large-scale weather patterns and exhibited little annual variation for the duration of our study. However, we were able to identify important thresholds of landscape composition linked to occupancy at both scales. We found landcover of CRP, Grass, and Shrub influenced second-order habitat use (Johnson 1980) at both spatial scales, but the effect of grassland fragmentation operated at the large scale, with greater variation among ecoregions at the small scale. Covariate relationships may be useful for informing conservation practices at different spatial scales and for identifying habitat factors that influence species distributions (Hagen et al. 2016, Pavlacky et al. 2017). These findings highlight potential benefits to LEPC of restoring prairie landscapes through CRP or other conservation programs that convert tilled agricultural land back to grassland. Additionally, our findings elucidate the importance of landscape scale management, e.g., tree removal and mitigation of development, to increase the likelihood of occupancy.

Annual variation in site occupancy

The pattern of implicit dynamics suggested that the geographic range of LEPC, as estimated by the probability of occupancy, remained stable at the large scale but fluctuated temporally at the small scale. Stable range dynamics at the large scale with fluctuating range dynamics at the small scale may be expected because site occupancy exhibits greater correlation with abundance at smaller spatial scales (Noon et al. 2012). We found little evidence of annual variation in large-scale occupancy (ψ) of LEPC between 2012 and 2016, but discovered moderate evidence for annual variation in the probability of small-scale occupancy (θ; Fig. 1). However, range contraction and expansion oscillated temporally in θ, with a 27% contraction between 2012 and 2013 (λ = 0.73; 90% CI = 0.53, 1.01) and a 49% expansion between 2013 and 2015 (λ = 1.49; 90% CI = 1.21, 1.82), and marginal temporal changes between other years (Fig. 1). Best models for annual variation in θ included the additive effect of ecoregion, suggesting θ was greater in the SGPR ecoregion than in others, but with the same pattern of dynamics among all ecoregions (Fig. 1). Annual estimates of θ in the SGPR and MGPR regions paralleled regional abundance estimates, but θ for the SOPR and SSPR regions did not (McDonald et al. 2016). Nevertheless, observed annual variation in site occupancy was relatively small compared to landscape effects observed in the range-wide covariate analysis.

Multiscale covariate relationships

Overall, we found strong support for hypotheses that LEPC occupancy was correlated with landscape composition and configuration, and conservation programs involving CRP at both spatial scales. We found some support for hypotheses that LEPC occupancy was correlated with Development at the small scale, and the effect was several orders of magnitude greater than at the large scale. The exploratory analysis showed a strong negative effect of Development with trivial probabilities of θ when landscapes contained > 10% of this landcover. Our results did not support hypotheses for shifts in extent of occurrence in response to drought-related covariates. These patterns are commensurate with other work, which has demonstrated reductions in abundance as it related to drought but site occupancy was largely unchanged especially at these spatial scales (Ross et al. 2016a, Fritts et al. 2018). Therefore, management of landscape mosaics (Haukos and Zavaletta 2016) may be more influential than drought-related climatic patterns (Grisham et al. 2016) at both spatial scales.

Lack of drought-related effects may not be surprising, given we observed low annual variation in ψ, and only moderate annual variation in θ. Although interactive models of landcover and drought were not well supported, overall, occupancy shifted from CRP to Shrub during drought conditions. This pattern was contradictory to regional variation, i.e., SGPR, in CRP use by radiomarked LEPCs during drought (Sullins et al. 2018). However, the general importance of CRP to occupancy was evident but we were unable to detect strong relationships with drought. LEPC abundance was correlated with lagged climatic effects (Ross et al. 2016a, b), and we expected that LEPC abundance may have been more sensitive to climatic conditions than occupancy.

Landscape composition and configuration, and conservation efforts were the most important predictors of ψ when modeling pooled data over all ecoregions. However, we found little support for relative effects of anthropogenic disturbance or drought-related patterns at this scale. It is possible that our modeling was unable to detect correlations between occupancy and Development because of the legacy effects of habitat conversion that occurred decades prior to our work. Such large-scale conversion may have previously displaced local populations and those areas may have been excluded from our sampling frame. Shrub, GrassPatch, and CRP were the most important correlates of ψ in these pooled data for the range-wide analysis. This confirms previous assertions that landscape-level habitat loss and fragmentation are among the most important factors for long-term population persistence of LEPC (Haukos and Zavaletta 2016). The most suitable landcover mosaics within 225 km² landscapes were composed of average grassland patches > 19 ha embedded in a matrix with Shrub > 22% or CRP > 12%. These landscape relationships were more important than patch sizes of GeneralHabitat or NativeHabitat, suggesting that CRP did not contribute to the patch sizes of native habitat. Because such a large proportion of grassland patches were small (75% of grassland patches ≤ 30 ha; 90% ≤ 66 ha), mean patch size may not reflect the most likely landscape metric of interest when evaluating LEPC habitat, however it appears to provide an adequate measure of fragmentation (Spencer et al. 2017).

This mosaic pattern of landscape configuration and composition suggested that habitat fragmentation was a more important determinant for LEPC presence in Grass, whereas total area of habitat, regardless of fragmentation, was more important in CRP. Because patch size was more important than landcover for Grass, landscape management to maintain large patches of grassland may be an important conservation strategy in modified agricultural landscapes (Kareiva and Wennergen 1995). In contrast, the effects of CRP on ψ suggested CRP functioned primarily as between-patch matrix habitat and did not contribute to the patch sizes of grassland vegetation. Landscape management to implement CRP and maintain high landcover of sand sagebrush and shinnery oak shrubland may be important conservation strategies, and losses of these habitat types in any configuration is expected to reduce ψ. We speculate that the amount of CRP or Shrub may increase landscape permeability and facilitate LEPC dispersal between patches of native grassland (Ricketts 2001, Niemuth 2011).

Models with CRP consistently outperformed those containing LPCI-prescribed grazing on large-scale occupancy. However, we observed large additive effects of prescribed grazing and patch size of native vegetation on ψ in less supported models with predicted increases (11%) in occupancy as the proportion of prescribed grazing increased (1% of landcover enrolled), although these were ~50% less than previously observed (Hagen et al. 2016). Because of close connections between habit loss and local extinction, landscape management of CRP may have greater potential to increase the extent of occurrence in agricultural landscapes, compared to prescribed grazing to address habitat degradation in remnants. However, because habitat degradation is expected to become increasingly important above a threshold of available habitat (Fischer and Lindenmayer 2007), future research to investigate the role of prescribed grazing along gradients of landscape composition may provide important insights. Hagen et al. (2016) predicted the effect of prescribed grazing would operate at the small-scale, and it is possible that the effects of grazing management operate on third-order habitat use, i.e., use within an organisms home-range, rather than second-order habitat use, i.e., home ranges selected within geographic distribution, at the landscape scale (Haukos and Zavaletta 2016). Lesser Prairie-Chicken third-order space use has been empirically linked and positively correlated with LPCI-prescribed grazing plans, and threshold effects have been demonstrated for stocking density and patch-burn frequency (Lautenbach 2017, Kraft et al. 2020). However, the sparseness of our data may have reduced the ability to detect these effects a priori among the large suite of covariates (n = 31) especially those that affect connectivity and fragmentation (Hagen et al. 2016).

Landscape composition and conservation efforts were the most important predictors of θ in these pooled data. There were ecoregional interactions with CRP, Shrub, and Grass when data were pooled across ecoregions. The odds ratios indicated that θ increased by 11%, 5%, 5%, and 0.8% for every 1% increase in of the amount of CRP in the SGPR, MGPR, SSPR, and SOPR, respectively (Fig. 4). The interaction between ecoregion and CRP, suggested the landscape management of CRP in 56.25 km² landscapes will be most effective in the SGPR, moderately effective in the MGPR and SSPR, and least effective in the SOPR. The positive and additive effects of Shrub were parallel, and the odds of θ increased by 5% for every 1% increase in Shrub (Fig. 5), in the SOPR. In the SGPR, θ exceeded the inflection point of the logistic relationship when Grass was > 58%. The odds ratio indicated that θ increased by 2% for every 1% increase of Grass in the MGPR. The interaction effects suggest that managing for greater grassland landcover may be less important in the SOPR and SSPR than in the SGPR and MGPR. This finding may not be surprising because sand sagebrush and sand shinnery oak are codominant with grasses in SSPR and SOPR, respectively, and these shrubs are crucial to long-term resilience of vegetation communities in these more drought susceptible ecoregions (Hagen and Elmore 2016). Large-scale loss of shrub landcover to grass has not been conducive to LEPCs in these regions (Haukos and Zavaleta 2016). The small effect of Grass on small-scale occupancy suggested the amount of grassland habitat at the local scale was not a limiting factor in the SOPR and SSPR, and the large effect of Shrub suggested a greater reliance on shrubland habitat in these regions.

Exploratory model selection

Our a priori model selection indicated that large patches of native grassland, CRP, and shrubland cover were associated with large-scale occupancy. Exploring additional covariates suggested that tree cover was also highly associated with LEPC occupancy at the large scale, likely influencing distributional limits. In our post-hoc exploratory analyses of ψ, we discovered that occupancy decreased 13% for every 1% increase in Woodland-10. This confirms the findings of field studies of LEPC space use and movement where individual birds exhibited strong avoidance to tall woody cover (Boggie et al. 2017, Lautenbach et al. 2017). These field studies focused on space use relative to honey mesquite (Prosopris glandulosa) and eastern red cedar (Juniperus virginiana) in the SOPR and MGPR, respectively. Interestingly, our finding that occupancy was negatively correlated with the presence of relatively low levels (< 5 and < 10%) of woodland cover was commensurate with fine-scale avoidance of other prairie-grouse (McNew et al. 2012, Severson et al. 2017). At broader scales the increase in woodland cover has been linked to declines in lek attendance of LEPC and other prairie-grouse species (Fuhlendorf et al. 2002, McNew et al. 2012, Baruch-Mordo et al. 2013, Hagen et al. 2019). Thus, while increasing available prairie can be achieved through CRP (or other surrogate grassland restoration) tree removal may offer additional avenues for restoration (Fuhlendorf et al. 2002, Hagen et al. 2019).

In our post-hoc exploratory analyses of θ, we found limited evidence that Development was negatively influencing LEPC occupancy and it varied regionally. The larger absolute effects were observed in the SGPR and MGPR, but in either case with every 1% (0.56 km²) increase in Development the odds of occupancy decreased by 22%. Various studies have examined space use and movements of individual LEPCs (Pitman et al. 2005, Pruett et al. 2009, Hagen et al. 2011, Plumb et al. 2018), but ours is the first to rigorously estimate occupancy as a function of human development at broad scales. Growing evidence suggests that individuals tend to exhibit a tolerance threshold measured in minimum distances, and that may translate to reduced occupancy as the proportion of the landscape is occupied by human developments (Bartuszevige and Daniels 2016).

Simulated effect of changes to CRP-enrolled land

Given the positive relationships found between LEPC occupancy and amounts of land enrolled in CRP, we described the effects of changes in CRP were expected to have on occupied area of LEPC across their range. LEPC occupancy was most sensitive to changes in CRP levels in the SGPR. Disproportionately large gains in occupancy associated with increases in CRP enrollment in the SGPR were exemplified in the CRP-occupancy relationship for ψ and θ (Figs. 3C, 4D). In both cases, when CRP was rare on the landscape (approximately 0–15% coverage), as little as 5–10% additional percent coverage of CRP was estimated to produce an approximate doubling in the probability of occupancy in the ecoregion. This strong relationship between occupancy and CRP may be partially explained by the relative scarcity of CRP in this ecoregion (Table 8).

We caution that the distribution of LEPC at the range-wide scale may be constrained by factors that we did not consider. In particular, this analysis assumes that changes in the probability of occupancy translate directly to changes in the area occupied by LEPC, with no assumed barriers to dispersal and colonization of previously unoccupied areas. We also made the simplifying assumptions that CRP expansion could only occur within quadrants where some CRP was already present and that changes in CRP enrollment did not affect the values of other environmental covariates also used as predictor variables in the occupancy models. A spatially explicit simulation could be conceived to add (or remove) CRP according to some process (random or supervised) and modify the values of other environmental covariates where CRP was added (or removed). Despite these assumptions, our results suggested that higher levels of CRP enrollment would expand LEPC occupied area and illustrated potential impacts that conservation policies have on the distribution of a species of high conservation concern.

Our work confirms prior assertions that broad-scale habitat loss and fragmentation are primary factors eroding the long-term population persistence of LEPC (Van Pelt et al. 2013). Although threats to functional prairie ecosystems persist through loss and degradation of habitat, our work provides insights that may prove useful in maintaining, improving, or restoring these systems through conservation actions taken on private lands. Notably, CRP appears to serve as a surrogate for prairie habitat in some ecoregions, especially in conjunction with larger extant patches of native habitat. These findings support previous research demonstrating the long-term effectiveness of CRP as a tool to reduce grassland habitat fragmentation at landscape scales (Spencer et al. 2017). The results suggest managers may be able to increase the current distribution of LEPC through increasing the amount of available grasslands. If CRP is combined with strategies to stem the tide of woody encroachment, it should yield increases in patch size and usable space for LEPCs. Similarly, targeted approaches to siting of energy developments and restoration of abandoned sites could prove useful in shifting the distribution of LEPCs.

LITERATURE CITED

Arnold, T. W. 2010. Uninformative parameters and model selection using Akaike’s information criterion. Journal of Wildlife Management 74:1175-1178. https://doi.org/10.1111/j.1937-2817.2010.tb01236.x

Bartuszevige, A., and A. Daniels. 2013. Evaluation of Farm Bill Programs to benefit Lesser Prairie-Chicken (Tympanuchus pallidicinctus). U.S. Department of Agriculture, Natural Resources Conservation Service, Washington, D.C., USA.

Bartuszevige, A. M., and A. Daniels. 2016. Impacts of energy development, anthropogenic structures, and land use change on Lesser Prairie-Chickens. Pages 205-220 in D. A. Haukos, and C. W. Boal, editors. Ecology and conservation of Lesser Prairie-Chickens. Studies in Avian Biology (no. 48), CRC, Boca Raton, Florida, USA. https://doi.org/10.1201/b19678-15

Baruch-Mordo, S., J. S. Evans, J. P. Severson, D. E. Naugle, J. D. Maestas, J. M. Kiesecker, M. J. Falkowski, C. A. Hagen, and K. P. Reese. 2013. Saving Sage-Grouse from the trees: a proactive solution to reducing a key threat to a candidate species. Biological Conservation 167:233-241. https://doi.org/10.1016/j.biocon.2013.08.017

BirdLife International (BLI). 2018. State of the world’s birds: taking the pulse of the planet. BLI, Cambridge, UK.

Boggie, M. A., C. R. Strong, D. Lusk, S. A. Carleton, W. R. Gould, R. L. Howard, C. Nichols, M. Falkowski, and C. Hagen. 2017. Impacts of mesquite distribution on seasonal space use of Lesser Prairie-Chickens. Rangeland Ecology and Management 70:68-77. https://doi.org/10.1016/j.rama.2016.09.006

Brennan, L. A., and W. P. Kuvlesky. 2005. North American grassland birds: an unfolding conservation crisis? Journal of Wildlife Management 69:1-13. https://doi.org/10.2193/0022-541x(2005)069<0001:nagbau>2.0.co;2

Briske, D. D., B. T. Bestelmeyer, J. R. Brown, M. W. Brunson, T. L. Thurow, and J. A. Tanaka. 2017. Assessment of USDA-NRCS rangeland conservation programs: recommendation for an evidence-based conservation platform. Ecological Applications 27:94-104. https://doi.org/10.1002/eap.1414

Bromaghin, J. F., T. L. McDonald, and S. C. Amstrup. 2013. Plausible combinations: an improved method to evaluate the covariate structure of Cormack-Jolly-Seber mark-recapture models. Open Journal of Ecology 3:11-22. https://doi.org/10.4236/oje.2013.31002

Burnham, K. P., and D. R. Anderson. 2002. Model selection and multimodel inference: a practical information-theoretic approach. Springer-Verlag, New York, New York, USA.

Burnham, K. P., D. R. Anderson, G. C. White, C. Brownie, and K. H. Pollock. 1987. Design and analysis methods for fish survival experiments based on release-recapture. American Fisheries Society Monograph 5. American Fisheries Society, Bethesda, Maryland, USA.

Caro, T. M. 2010. Conservation by proxy: indicator, umbrella, keystone, flagship, and other surrogate species. Island Press, Washington, D.C., USA.

Chalfoun, A. D., and T. E. Martin. 2007. Assessments of habitat preferences and quality depend on spatial scale and metrics of fitness. Journal of Applied Ecology 44:983-992. https://doi.org/10.1111/j.1365-2664.2007.01352.x

Chamberlin, T. C. 1965. The method of multiple working hypotheses. Science 148:754-759. https://doi.org/10.1126/science.148.3671.754

Cody, M. L. 1985. Habitat selection in birds. Academic, New York, New York, USA.

Doherty, P. F., G. C. White, and K. P. Burnham. 2012. Comparison of model building and selection strategies. Journal of Ornithology 152:317-323. https://doi.org/10.1007/s10336-010-0598-5

Engle, D. M., B. R. Coppedge, and S. D. Fuhlendorf. 2008. From dust bowl to green glacier: human activity and environmental change in the Great Plains grasslands. Pages 253-271 in O. W. Van Auken, editor. Western North American juniperus communities: a dynamic vegetation type. Springer, New York, New York, USA. https://doi.org/10.1007/978-0-387-34003-6_14

Fahrig, L., J. Baudry, L. Brotons, F. G. Burel, T. O. Crist, R. J. Fuller, C. Sirami, G. M. Siriwardena, and J.-L. Martin. 2011. Functional landscape heterogeneity and animal biodiversity in agricultural landscapes. Ecology Letters 14:101-112. https://doi.org/10.1111/j.1461-0248.2010.01559.x

Fields, T. L., G. C. White, W. C. Gilgert, and R. D. Rodgers. 2006. Nest and brood survival of Lesser Prairie-Chickens in west central Kansas. Journal of Wildlife Management 70:931-938. https://doi.org/10.2193/0022-541X(2006)70[931:NABSOL]2.0.CO;2

Fischer, J., and D. B. Lindenmayer. 2007. Landscape modification and habitat fragmentation: a synthesis. Global Ecology and Biogeography 16:265-280. https://doi.org/10.1111/j.1466-8238.2007.00287.x

Fritts, S. R., B. A. Grisham, R. D. Cox, C. W. Boal, D. A. Haukos, P. McDaniel, C. A. Hagen, and D. U. Greene. 2018. Interactive effects of severe drought and grazing on the life history cycle of a bioindicator species. Ecology and Evolution 8:9550-9562. https://doi.org/10.1002/ece3.4432

Fuhlendorf, S. D., W. C. Harrell, D. M. Engle, R. G. Hamilton, C. A. Davis, and D. M. Leslie. 2006. Should heterogeneity be the basis for conservation? Grassland bird response to fire and grazing. Ecological Applications 16:1706-1716. https://doi.org/10.1890/1051-0761(2006)016[1706:SHBTBF]2.0.CO;2

Fuhlendorf, S. D., A. J. W. Woodward, D. M. Leslie, and J. S. Shackford. 2002. Multi-scale effects of habitat loss and fragmentation on Lesser Prairie-Chicken populations of the US Southern Great Plains. Landscape Ecology 17:617-628. https://doi.org/10.1023/A:1021592817039

George, T. L., and S. Zack. 2001. Spatial and temporal considerations in restoring habitat for wildlife. Restoration Ecology 9:272-279. https://doi.org/10.1046/j.1526-100x.2001.009003272.x

Grisham, B. A., A. J. Godar, and C. P. Griffin. 2016. Climate change. Pages 221-242 in D. A. Haukos, and C. W. Boal, editors. Ecology and conservation of Lesser Prairie-Chickens. Studies in Avian Biology (no. 48), CRC, Boca Raton, Florida, USA.

Hagen, C. A., and D. Elmore. 2016. Conclusions, synthesis, and future directions. Pages 345- 351 in D. A. Haukos and C. W. Boal, editors. Ecology and conservation of Lesser Prairie-Chickens. Studies in Avian Biology (no. 48) CRC, Boca Raton, Florida, USA.

Hagen, C. A., D. C. Pavlacky Jr, K. Adachi, F. E. Hornsby, T. J. Rintz, and L. L. McDonald. 2016. Multiscale occupancy modeling provides insights into range-wide conservation needs of Lesser Prairie-Chicken (Tympanuchus pallidicinctus). Condor 118:597-612. https://doi.org/10.1650/CONDOR-16-14.1

Hagen, C. A., J. C. Pitman, T. M. Loughin, B. K. Sandercock, R. J. Robel, and R. D. Applegate. 2011. Impacts of anthropogenic features on habitat use by Lesser Prairie-Chickens. Pages 63-75 in B. K. Sandercock, K. Martin, and G. Segelbacher, editors. Ecology, conservation, and management of grouse. Studies in Avian Biology (no. 39), University of California Press, Berkeley, California, USA. https://doi.org/10.1525/9780520950573-007

Hagen, C. A., K. A. Taylor, A. M. Bartuszevige, A. B. Daniels, M. T. DeLeon, and M. Bogaerts. 2019. Slowing the tide of mesquite invasion: using a bioindicator species to deliver conservation triage. Journal of Arid Environments 168:46-55. https://doi.org/10.1016/j.jaridenv.2019.05.004

Haukos, D. A., and J. C. Zavaletta. 2016. Habitat. Pages 99-132 in D. A. Haukos and C. W. Boal, editors. Ecology and conservation of Lesser Prairie-Chickens. Studies in Avian Biology (no. 48), CRC, Boca Raton, Florida, USA. https://doi.org/10.1201/b19678

Hovick, T. J., R. D. Elmore, S. D. Fuhlendorf, D. M. Engle, and R. G. Hamilton. 2015. Spatial heterogeneity increases diversity and stability in grassland bird communities. Ecological Applications 25:662-672. https://doi.org/10.1890/14-1067.1

Hutto, R. L. 1985. Habitat selection by nonbreeding, migratory land birds. Pages 455-476 in M. L. Cody, editor. Habitat selection in birds. Academic, Orlando, Florida, USA.

Johnson, D. H. 1980. The comparison of usage and availability measurements for evaluating resource preference. Ecology 61:65-71. https://doi.org/10.2307/1937156

Kareiva, P., and U. Wennergen. 1995. Connecting landscape patterns to ecosystem and population process. Nature 373:299-302. https://doi.org/10.1038/373299a0

Kraft, J. D., D. A. Haukos, M. R. Bain, M. B. Rice, S. Robinson, D. S. Sullins, C. A. Hagen, and J. C. Pitman. 2020. Using grazing to manage herbaceous structure for a heterogeneity-dependent bird. Journal of Wildlife Management, in press.

Laake, J. L. 2013. RMark: an R interface for analysis of capture-recapture data with MARK. Alaska Fisheries Science Center Processed Report 2013-01. Alaska Fisheries Science Center, National Oceanic and Atmospheric Administration, National Marine Fisheries Service, Seattle, Washington, USA.

Lautenbach, J. D. 2017. The role of fire, microclimate, and vegetation in Lesser Prairie-Chicken habitat selection. Thesis. Kansas State University, Manhattan, Kansas, USA.

Lautenbach, J. M., R. T. Plumb, S. G. Robinson, C. A. Hagen, D. A. Haukos, and J. C. Pitman. 2017. Lesser Prairie-Chicken avoidance of trees in a grassland landscape. Rangeland Ecology and Management 70:78-86 https://doi.org/10.1016/j.rama.2016.07.008

MacKenzie, D. I., J. D. Nichols, J. A. Royle, K. H. Pollock, L. L. Bailey, and J. E. Hines. 2018. Occupancy estimation and modeling: inferring patterns and dynamics of species occurrence. Second edition. Academic, London, UK. https://doi.org/10.1016/C2012-0-01164-7

McDonald, L., G. Beauprez, G. Gardner, J. Griswold, C. Hagen, F. Hornsby, D. Klute, S. Kyle, J. Pitman, T. Rintz, D. Schoeling, and B. Van Pelt. 2014. Range-wide population size of the Lesser Prairie-Chicken: 2012 and 2013. Wildlife Society Bulletin 38:536-546. https://doi.org/10.1002/wsb.417

McDonald, L., K. Nasman, T. Rintz, F. Hornsby, and G. Gardner. 2016. Range-wide population size of the Lesser Prairie-Chicken: 2012, 2013, 2014, 2015, and 2016. Western EcoSystems Technology, Inc., Laramie, Wyoming, USA.

McNew, L. B., T. J. Prebyl, and B. J. Sandercock. 2012. Effects of rangeland management on the site occupancy dynamics of prairie-chickens in a protected prairie preserve. Journal of Wildlife Management 76:38-47. https://doi.org/10.1002/jwmg.237

Miller, R. F., D. E. Naugle, J. D. Maestas, C. A. Hagen, and G. Hall. 2017. Special Issue: Targeted woodland removal to recover at-risk grouse and their sagebrush steppe and prairie ecosystems. Rangeland Ecology & Management 70:1-8. https://doi.org/10.1016/j.rama.2016.10.004

Nichols, J. D., L. L. Bailey, A. F. O'Connell, N. W. Talancy, E. H. C. Grant, A. T. Gilbert, E. M. Annand, T. P. Husband, and J. E. Hines. 2008. Multi-scale occupancy estimation and modelling using multiple detection methods. Journal of Applied Ecology 45:1321-1329. https://doi.org/10.1111/j.1365-2664.2008.01509.x

Niemuth, N. D. 2011. Spatially explicit habitat models for prairie grouse. Pages 3-20 in B. K. Sandercock, K. Martin, and G. Segelbacher, editors. Ecology, conservation, and management of grouse. Studies in Avian Biology (no. 39), University of California Press, Berkeley, California, USA. https://doi.org/10.1525/9780520950573-003

Noon, B. R., L. L. Bailey, T. D. Sisk, and K. S. McKelvy. 2012. Efficient species-level monitoring at the landscape scale. Conservation Biology 26:432-441. https://doi.org/10.1111/j.1523-1739.2012.01855.x

Pavlacky, D. C., Jr., J. A. Blakesley, G. C. White, D. J. Hanni, and P. M. Lukacs. 2012. Hierarchical multi-scale occupancy estimation for monitoring wildlife populations. Journal of Wildlife Management 76:154-162. https://doi.org/10.1002/jwmg.245

Pavlacky, D. C., Jr., P. M. Lukacs, J. A. Blakesley, R. C. Skorkowsky, D. S. Klute, B. A. Hahn, V. J. Dreitz, T. L. George, and D. J. Hanni. 2017. A statistically rigorous sampling design to integrate avian monitoring and management within bird conservation regions. PLoS ONE 12:e0185924. https://doi.org/10.1371/journal.pone.0185924

Pitman, J. C., C. A. Hagen, R. J. Robel, T. M. Loughin, and R. D. Applegate. 2005. Location and success of Lesser Prairie-Chicken nests in relation to vegetation and human disturbance. Journal of Wildlife Management 69:1259-1269. https://doi.org/10.2193/0022-541X(2005)069[1259:LASOLP]2.0.CO;2

Plumb, R. T., J. M. Lautenbach, S. Robinson, D. A. Haukos, V. L. Winder, C. A. Hagen, D. S. Sullins, J. C. Pitman, and D. K. Dahlgren. 2018. Lesser Prairie-Chicken space use in relation to anthropogenic structures. Journal of Wildlife Management 83:216-230. https://doi.org/10.1002/jwmg.21561

Powell, L. A. 2007. Approximating variance of demographic parameters using the delta method: a reference for avian biologists. Condor 109:949-954. https://doi.org/10.1093/condor/109.4.949

Pruett, C. L., M. A. Patten, and D. A. Wolfe. 2009. Avoidance behavior by prairie-grouse: implications for development of wind energy. Conservation Biology 23:1253-1259. https://doi.org/10.1111/j.1523-1739.2009.01254.x

R Core Team. 2017. R: a language and environment for statistical computing. Version 3.3.0. R Foundation for Statistical Computing, Vienna, Austria. [online] URL: https://www.r-project.org

Ribic, C. A., R. R. Koford, J. R. Herkert, D. H. Johnson, N. Niemuth, D. E. Naugle, K. K. Bakker, D. W. Sample, and R. B. Renfrew. 2009. Area sensitivity in North American grassland birds: patterns and processes. Auk 126:233-244. https://doi.org/10.1525/auk.2009.1409

Ricketts, T. H. 2001. The matrix matters: effective isolation in fragmented landscapes. American Naturalist 158:87-99. https://doi.org/10.1086/320863

Ross, B. E., D. A. Haukos, C. A. Hagen, and J. C. Pitman. 2016a. Landscape composition creates a threshold influencing Lesser Prairie-Chicken population resilience to extreme drought. Global Ecology and Conservation 6:179-188. https://doi.org/10.1016/j.gecco.2016.03.003

Ross, B. E., D. Haukos, C. Hagen, and J. Pitman. 2016b. The relative contribution of climate to changes in Lesser Prairie-Chicken abundance. Ecosphere 7:e01323. https://doi.org/10.1002/ecs2.1323

Samson, F. B., F. L. Knopf, and W. R. Ostlie. 2004. Great Plains ecosystems: past, present, and future. Wildlife Society Bulletin 32:6-15. https://doi.org/10.2193/0091-7648(2004)32[6:GPEPPA]2.0.CO;2

Severson, J. P., C. A. Hagen, J. D. Tack, J. D. Maestas, D. E. Naugle, J. T. Forbes, and K. P. Reese. 2017. Better living through conifer removal: a demographic analysis of sage- grouse vital rates. PLoS ONE 12:e0174347. https://doi.org/10.1371/journal.pone.0174347

Spencer, D., D. Haukos, C. Hagen, M. Daniels, and D. Goodin. 2017. Conservation Reserve Program mitigates grassland loss in the Lesser Prairie-Chicken range of Kansas. Global Ecology and Conservation 9:21-38. https://doi.org/10.1016/j.gecco.2016.11.004

Sullins, D. S., J. D. Kraft, D. A. Haukos, S. G. Robinson, J. H. Rietz, R. T. Plumb, J. M. Lautenbach, J. D. Lautenbach, B. K. Sandercock, and C. A. Hagen. 2018. Demographic consequences of Conservation Reserve Program grasslands for Lesser Prairie-Chickens. Journal of Wildlife Management 82:1617-1632. https://doi.org/10.1002/jwmg.21553

United States Department of Agriculture (USDA). 2014. Common land unit geospatial data. Memorandum of understanding between the USDA and the Western Association of Fish and Wildlife Agencies, signed 21 February 2014. USDA, Farm Service Agency, Economic and Policy Analysis, and Commodity Credit Corporation, Washington, D.C., USA.

United States Fish and Wildlife Service (USFWS). 2011. Conference report for the Natural Resources Conservation Service Lesser Prairie-Chicken Initiative (LPCI). U.S. Department of Agriculture, Natural Resources Conservation Service and U.S. Department of Interior, Fish and Wildlife Service, Washington D.C., USA.

United States Fish and Wildlife Service (USFWS). 2016. U.S. Fish and Wildlife Service removes Lesser Prairie-Chicken from List of Threatened and Endangered Species in accordance with Court Order. USFWS, Washington, D.C., USA. [online] URL: https://www.fws.gov/news/ShowNews.cfm?ref=u.s.-fish-and-wildlife-service-removes-lesser-prairie-chicken-from-list-o&_ID=35739

Van Pelt, W. E., S. Kyle, J. Pitman, D. Klute, G. Beauprez, D. Schoeling, A. Janus, and J. B. Haufler. 2013. The Lesser Prairie-Chicken range-wide conservation plan. Western Association of Fish and Wildlife Agencies, Cheyenne, Wyoming, USA.

White, G. C., and K. P. Burnham. 1999. Program MARK: survival estimation from populations of marked animals. Bird Study 46:S120-S139. https://doi.org/10.1080/00063659909477239