Multiscale occupancy of the Lesser Prairie-Chicken: the role of private lands in conservation of an imperiled bird

Grasslands are one of the most imperiled ecosystems globally, and the Lesser Prairie-Chicken (Tympanuchus pallidicinctus) is an iconic grassland-obligate species with high conservation priority in the USA. Lesser Prairie-Chicken conservation is compounded by its requirement for a spatial hierarchy of heterogeneous habitats, coupled with nearly all (> 95%) of its range being privately owned. The U.S. Department of Agriculture currently offers technical and financial resources that facilitate prairie restoration, e.g., Conservation Reserve Program (CRP), and improve habitat quality and ecosystem services, e.g., Environmental Quality Improvement Program, on private lands. We modeled Lesser Prairie-Chicken occupancy at two scales relative to covariates that described landscape composition and configuration, anthropogenic development, drought-related climatic conditions, and conservation efforts from 2012 to 2016. Large-scale (225 km2) occupancy was most associated with shrubland, grassland patch size, and CRP range-wide. Patterns of small-scale (56.25 km2) occupancy varied regionally, but key covariates included shrubland, grassland, and CRP landcover. These covariate relationships may be useful for identifying conservation practices at different spatial scales and habitat factors that influence Lesser Prairie-Chicken distributions ecoregionally. Notably, CRP-enrolled lands appear to serve as a surrogate for prairie habitat in some ecoregions, especially in conjunction with larger extant patches of native habitat. Although not as influential as CRP at large scales, every 2.25 km2 of prescribed grazing increased the odds of site occupancy by 11%. In addition to supported covariates, we found that for every 0.56 km2 of industrial development at small scales and 2.25 km2 of woodland cover (10%-canopy) at large scales, odds of occupancy decreased by 22% and 13%, respectively. Our results suggest that increased amounts of native grassland and shrubland, and in particular higher levels of CRP enrollment could expand LEPC distribution by as much as 17% (1418–1744 km2). Moreover, our findings illustrate the potential for federal conservation policies to benefit the distribution of an imperiled species. Présence multiéchelle du Tétras pâle : rôle des terres privées dans la conservation d'un oiseau en péril RÉSUMÉ. Les prairies sont l'un des écosystèmes les plus menacés au monde, et la conservation du Tétras pâle (Tympanuchus pallidicinctus), espèce emblématique qui dépend des prairies, est prioritaire aux États-Unis. La conservation du Tétras pâle représente un énorme défi en raison du fait qu'il dépend d'une hiérarchie spatiale d'habitats hétérogènes et que la quasi-totalité (> 95 %) de son aire de répartition se trouve en terres privées. Le ministère étatsunien de l'agriculture offre actuellement des ressources techniques et financières qui encouragent la restauration des prairies, par exemple le Conservation Reserve Program (CRP), et améliorent la qualité des habitats et les services écosystémiques, par exemple l'Environmental Quality Improvement Program, sur les terres privées. Nous avons modélisé la présence du Tétras pâle à deux échelles en fonction de covariables décrivant la composition et la configuration du paysage, l'aménagement du territoire, les conditions climatiques liées à la sécheresse et les efforts de conservation de 2012 à 2016. La présence du tétras à grande échelle (225 km2) était surtout associée aux arbustes, à la taille des parcelles de prairie et au CRP dans l'ensemble de son aire de répartition. La tendance de la présence à petite échelle (56,25 km2) a varié selon les régions, mais les principales covariables comprenaient les arbustes, les prairies et l'étendue géographique du CRP. Ces relations entre covariables peuvent permettre aux gestionnaires d'identifier les pratiques de conservation à différentes échelles spatiales et les facteurs d'habitat qui influent sur la répartition écorégionale du Tétras pâle. En particulier, les terres inscrites au CRP semblent servir d'indicateurs d'habitat de prairies dans certaines écorégions, surtout en combinaison avec de grandes parcelles d'habitat naturel. Bien qu'ils ne soient pas aussi influents que les CRP à grandes échelles, chaque 2,25 km2 de pâturage dirigé augmentait de 11 % les chances d'occupation d'un site. En plus des covariables indicatrices, nous avons constaté que pour chaque 0,56 km2 de développement industriel à petites échelles et 2,25 km2 de couvert forestier (10 % de couvert) à grandes échelles, les chances d'occupation diminuaient de 22 % et 13 %, respectivement. Nos résultats indiquent qu'une augmentation des superficies de prairies et d'arbustes indigènes, et en particulier une adhésion accrue au CRP, favoriseraient une expansion de la répartition du Tétras pâle de 17 % (1 418-1 744 km2). De plus, nos résultats illustrent le potentiel des politiques fédérales de conservation à augmenter la répartition d'une espèce en péril.


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, grasslandobligate 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 , 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. 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 . Here, we developed predictive models of LEPC multiscale occupancy to investigate covariate relationships over time, expanding upon previous efforts that were temporally limited . 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 . 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).

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

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

Statistical analysis Sampling framework for multiscale occupancy
We aggregated and summarized LEPC detection data such that large-scale occupancy (ψ) corresponded to LEPC presence on 15km × 15-km grid cells, small-scale occupancy (θ) corresponded to species presence in four quarters of the grid cell (7.5-km × 7.5km 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 . 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 , 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.

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) smallscale 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 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 We approximated the variance of λ t using the delta method (Powell 2007) and calculated 90% log-normal CIs for λ t (Burnham et al. 1987). Table 2. Summary of the number of covariates allowed during the plausible-combinations stage of model selection for multiscale occupancy models fit to ecoregional subsets of the Lesser Prairie-Chicken (Tympanuchus pallidicinctus) data from the range-wide monitoring program, 2012-2016. The multiscale occupancy model included submodels for large-scale occupancy (ψ), small-scale occupancy (θ), and detection probability (p). When a quadratic effect was included, the main effect for that covariate was also included, resulting in two covariates in the model for each quadratic term. When an interaction effect was included, the main effect for each interacting covariate was also included, resulting in three covariates in the model for each interaction term. Ecoregion acronyms are defined in text.

Range-wide relationships between covariates and multiscale occupancy.
We used all range-wide data in the four ecoregions from 2012 to (McDonald et al. 2016), but did not use auxiliary data collected within SGPR and SOPR during 2015 ). 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. 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 highlikelihood 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  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 largescale 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

Simulated effect of changes to amount of CRPenrolled 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 CRPenrolled lands in increments of 5%. For example, to simulate the amount and distribution of CRP enrollment under an 80%-ofbaseline 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 ). 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 with an estimated variance of .
The change in area occupied (ΔAO k ) was calculated as , with an estimated variance of

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).  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 plausiblecombination 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   (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   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  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 smallscale occupancy for the 2016 sampling frame within the occupied

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 secondorder 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 , 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 ). 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 droughtrelated 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 longterm 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 .
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.  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 secondorder 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 .
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 largescale 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. 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), 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 . 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   Limits (LCL and UCL, respectively) for the detection (p) of the lesser prairie-chicken from the best and second best model and range-wide monitoring program, 2012-2016. Table A1.3. Plausible combinations model selection for large-scale occupancy (ψ) of the lesser prairie-chicken from the range-wide monitoring program, 2012-2016. The model-selection metrics are the value of the minimized -2 log-likelihood function -2loge(ℒ)], parameter number (K), Akaike's Information Criterion adjusted for sample size (AICc), difference between model and minimum AICc value (ΔAICc) and AICc weight (wi). High-weight submodels with wi < 0.001 and high likelihood submodels with -2log(ℒ) < maximum -2log(ℒ) of high weight models are shown.  Table A1.4. Plausible combinations model selection for small-scale occupancy (θ) of the lesser prairie-chicken from the range-wide monitoring program, 2012-2016. The model-selection metrics are the value of the minimized -2 log-likelihood function -2loge(ℒ)], parameter number (K), Akaike's Information Criterion adjusted for sample size (AICc), difference between model and minimum AICc value (ΔAICc) and AICc weight (wi). High-weight submodels with wi < 0.001 and high-likelihood submodels with -2log(ℒ) < maximum -2log(ℒ) of high weight models are shown.