Crossing boundaries in conservation: land ownership and habitat in ﬂ uence the occupancy of an at-risk small mammal

. Ecological research is critical for informing management of at-risk species, for example, by identifying what drives occupancy of species across landscapes. However, restricting research to public land and omitting private land, as commonly occurs in ecological research, can bias inferences because important drivers of population and community patterns may vary with land ownership. We conducted a landscape-scale study of a species of conservation concern, white-tailed prairie dog (WTPD; Cynomys leucurus ), across multiple land ownerships in Wyoming. We quanti ﬁ ed how WTPD occupancy varied with both land ownership and biotic and abiotic factors. We established a baseline occupancy rate for WTPD in Wyoming, and quanti ﬁ ed how this baseline would have been biased if we had restricted our study to public land. We surveyed 440 sites throughout the Wyoming range of WTPD, which included sites on public (275), private (80), checkerboard (55), and Wind River Indian Reservation (30) land. We found that WTPD occupancy varied signi ﬁ cantly with land ownership and elevation, with WTPD being four times as likely to occupy private as public land at the median study area elevation. This difference in occupancy may have been driven by variation in habitat quality between ownership types, but we could not de ﬁ nitively deter-mine the underlying cause. Regardless of land ownership, WTPD occupancy increased with bare ground, but also when recent plant biomass (as estimated by NDVI) was higher than a site ’ s long-term average biomass. In other words, WTPD tended to occupy sparsely vegetated sites, but occupancy increased with short-term increases in biomass. The strong land ownership effect illustrates how study area delineation in relation to land ownership can in ﬂ uence research inferences and management actions, and why it is coun-terproductive to presuppose that certain land ownerships contain low-quality habitat. Ecologically, the variation in occupancy resulting from long-term site conditions, such as overall density of vegetation, relative to short-term site conditions, such as changes in plant biomass, can be used to improve long-term population monitoring for WTPD.


INTRODUCTION
Effective management and conservation of atrisk flora and fauna requires an understanding of population trends and drivers of population dynamics (Williams et al. 2002, Morrison et al. 2006. Research conducted to assess population parameters and inform management decisions is necessarily restricted to a study area; however, the delineation of the area itself can strongly influence research conclusions (Wiens 1989, Knight 1999, Spies et al. 2007). Wildlife research is often limited to public land due to challenges accessing private land (Haufler and Kernohan 2009), but critical drivers of population patterns, such as habitat quality, ecosystem management, or land-use, often vary with land ownership (Scott et al. 2001, Spies et al. 2007, Donnelly et al. 2016. Studies restricted to a single ownership may therefore yield incomplete and potentially misleading results if a substantial portion of the population or habitat is omitted, which in turn may influence management actions. Research that includes multiple land ownerships can thus more effectively inform conservation and management and, potentially, result in new ecological insights (Knight 1999, Haufler and Kernohan 2009, Kremen and Merenlender 2018. Predicting species' responses to potential stressors, such as from development or climate change, is also critical for successful conservation. If we understand the threats to a species within the context of the species' natural history, we can predict likely population responses based on the expected increase or decrease in threats overtime Murray 2011, Sih et al. 2011). Species monitoring should therefore incorporate environmental and anthropogenic factors expected to influence population metrics, such as abundance or occupancy.
Understanding how species respond to stressors is particularly important in grassland and shrubland ecosystems, which are underrepresented in habitat conservation and overrepresented in habitat loss worldwide (Hoekstra et al. 2005). These ecosystems have been extensively degraded and are continually threatened by factors such as oil and gas development and agriculture (Ellis et al. 2010, WGFD 2017. These factors often negatively impact native flora and fauna, but species responses vary (Vila et al. 2011, Northrup and Wittemyer 2013, Ehlers Smith et al. 2015, so it is critical to understand the relative importance of stressors for a particular species. In the temperate grasslands and shrublands of North America, prairie dogs (Cynomys spp.) are potential keystone species that improve habitat condition for other species in the community and are a vital prey source for many predators (Seglund et al. 2006, Delibes-Mateos et al. 2011).
The endangered black-footed ferret (Mustela nigripes) is an obligate prairie dog predator and thus highly dependent on healthy prairie dog colonies for survival (Campbell III et al. 1987). Prairie dogs are therefore an important component of ecosystem management especially because they are subject to numerous direct and indirect stressors (Delibes-Mateos et al. 2011).
The white-tailed prairie dog (WTPD; Cynomys leucurus) is a species of conservation concern that inhabits the grasslands and shrublands of the intermountain west, USA (Clark et al. 1971, Seglund et al. 2006, USFWS 2017. In addition to indirect impacts from habitat loss and degradation, WTPD has experienced direct losses from poisoning and especially introduction of nonnative sylvatic plague (Yersinia pestis; Seglund et al. 2006). This constellation of stressors resulted in WTPD being petitioned in 2002 for listing under the United States' Endangered Species Act (ESA; 16 U.S.C. § 1531); litigation ultimately extended the U.S. Fish and Wildlife Service's (USFWS) review of the petition until the end of 2017, when USFWS concluded that WTPD was not warranted for listing under the ESA (USFWS 2017, 82 FR 57562).
In response to the 2002 ESA petition and known threats to WTPD, the Western Association of Fish and Wildlife Agencies (WAFWA) formed a working group to conduct a WTPD conservation assessment (Seglund et al. 2006) and design a standardized rangewide monitoring protocol (Andelt et al. 2009, USFWS 2017. The protocol is based on WTPD occupancy because an occupancy approach is repeatable, statistically robust, and avoids the pitfalls associated with mapping colony boundaries or counting active burrows (MacKenzie et al. 2006, Andelt et al. 2009, USFWS 2017. To assess the status of WTPD, Wyoming, Colorado, and Utah, which comprise >98% of the predicted WTPD range, conduct standardized WTPD occupancy surveys (USFWS 2017).
We conducted a landscape-scale study across multiple land ownerships that assessed the drivers of WTPD occupancy in Wyoming, which contains approximately 77% of predicted WTPD range (USFWS 2017). WTPD habitat in Wyoming and rangewide encompasses multiple land ownerships, and habitat quality, management, and accessibility for surveys differ by these ownerships (Scott et al. 2001, Pocewicz et al. 2009). We stratified our study area by land ownership and addressed four primary objectives: (1) quantify how WTPD occupancy varies with land ownership while controlling for biotic and abiotic variables, (2) evaluate potential mechanistic drivers of occupancy, such as biotic, abiotic, and anthropogenic factors, (3) establish a Wyoming occupancy rate to provide a baseline for continued WTPD monitoring, and (4) compare conservation inferences from an analysis based on all land ownerships to an analysis restricted to public land. More broadly, we evaluated how inferences from wildlife research can be altered by study area delineation, such as limiting surveys to public land.

Study design
To select survey sites, we first developed a habitat suitability model for WTPD within their known or suspected range in Wyoming (Keinath et al. 2010; see Appendix S1: Table S1 for suitability model criteria). We overlaid a 500 × 500 m grid of survey quadrats over the suitability model and defined our sampling frame as all quadrats that were comprised of 50% or more suitable habitat. Since factors that affect WTPD ecology, such as habitat quality and land use, vary by land ownership in Wyoming (WY; Scott et al. 2001, Pocewicz et al. 2009, WGFD 2017, we divided potential sites into four distinct ownership strata: public, private, Wind River Indian Reservation (WRIR) of the Eastern Shoshone and Northern Arapaho tribes, and checkerboard ( Fig. 1). Public land (51.0% of suitable habitat in WY) is open to the public and is predominantly managed by the U.S. Bureau of Land Management (BLM) for multiple uses such as commercial (e.g., energy development), recreation, and conservation (BLM 2007). Private land (22.9% of suitable habitat in WY) is generally closed to public access and is managed by independent landowners predominately as rangelands for livestock grazing, although a variety of uses occur (WGFD 2017). The WRIR (4.5% of suitable habitat in WY) has limited public access and is managed by the Eastern Shoshone and Northern Arapaho Tribes for protection of Tribal resources, fish and wildlife conservation and recreation, mineral and energy production, and livestock grazing (Smith and Baldes 1982). The checkerboard (21.6% of suitable habitat in WY) is a swath of land along the Union Pacific Railroad/ Interstate highway 80 corridor that alternates between private and public ownership for adjacent square mile sections. The complex ownership pattern of the checkerboard restricts public access and results in a mix of public (BLM) and private land management (Haufler and Kernohan 2009). Finally, from the population of suitable sites, we selected sites with a random spatially balanced sample (balanced acceptance sampling [BAS]; Robertson et al. 2013). Sites were selected independently for each ownership stratum, and number of sites per stratum was proportional to stratum area.
To better assess the effect of energy development, a possible stressor for prairie dogs Collinge 2004, Sackett et al. 2012), we selected additional sites with high oil and gas development based on point locations of active well pads (Wyoming Oil and Gas Conservation Commission, http://wogcc.state.wy.us), which increased variation in energy development between sites. For analysis, the additional oil and gas sites were included in their corresponding ownership stratum, and effects were modeled with oil and gas predictors, rather than creating an additional energy development stratum (see Appendix S1 for additional details).

Field surveys
To survey WTPD across all-ownership stratum, we surveyed 500 × 500 m sites using a mix of aerial and ground surveys. Sites were surveyed either twice on the ground and once aerially (three occasions), twice on the ground (two occasions), or once aerially (one occasion). Surveys were conducted between 5/23/2016 and 8/ 29/2016 between approximately 07:00 and 17:00 h, except when wind speed exceeded approximately 23 mph (Beaufort scale level 5; World Meteorological Organization 1970) or rain was moderate or greater. Ground surveys used a double-observer design: two technicians conducted simultaneous, independent surveys, resulting in two survey occasions. Technicians searched for WTPD while walking the perimeter of each site and scanned the site's interior for WTPD for five minutes from each corner using v www.esajournals.org  10 × 42 power binoculars; technicians always surveyed the same site simultaneously, but, in order to maintain independence, walked the perimeter of each site in opposite directions and never surveyed from the same corner simultaneously (Andelt et al. 2009). We surveyed all private land sites aerially due to challenges acquiring private land access on a statewide scale. Aerial surveys were conducted on 5/19, 5/20, 7/26, and 7/27/2016, using Cessna 182 planes. Biologists and pilots both scanned for WTPD during three passes over each site, which comprised a single survey occasion (see Appendix S1 for additional flight methods). For all surveys, a site was only classified as occupied if ≥1 WTPD was visually observed within the site. We excluded other methods, like aural detection of calls or signs of presence, due to difficulties determining the location of a call or age of sign such as scat and burrows.

Occupancy modeling
For detection and occupancy probability modeling, we chose predictors a priori based on prairie dog ecology and a literature review; all predictor interactions represented a priori hypotheses. We focused on factors likely to directly or indirectly affect food availability or the costs of above-ground activity, such as predation risk, which are well-established controls on small mammal behavior and fitness (Longland and Price 1991, Orrock et al. 2004, Facka et al. 2010, Avila-Flores et al. 2012. We estimated detection probability as a function of survey method using 96 sites that were surveyed both aerially and on the ground (Table 1). This approach assumed that survey method did not interact with land ownership to influence detection probability, since the 96 sites that were surveyed both aerially and on the ground excluded private land (Table 1). This assumption is reasonable because (1) all sites were filtered based on our habitat suitability model, which reduced variation between ownership strata in variables expected to influence detection, such as slope, and (2) site-and time-specific variables, such as temperature, were included to model additional variation in detection (Table 2; Andelt et al. 2009).
To identify drivers of WTPD occupancy and reduce the likelihood of a spurious ownership effect, we assessed a suite of remotely sensed predictors summarized at multiple temporal and spatial scales (Table 3). Climate predictors (Facka et al. 2010, Savage et al. 2011, Avila-Flores et al. 2012, Eads and Hoogland 2017 were summarized from monthly Daymet rasters (1-km resolution, accessed May 2017). Daymet data are interpolated from daily weather observations from ground-based meteorological stations (Thornton et al. 2014). Topography predictors (Andelt et al. 2009, Magle andCrooks 2009) and substrate predictors (Savage et al. 2011) were derived from a Digital Elevation Model (30-m resolution) and Natural Resources Conservation Service STATSGO data (30-m resolution), respectively. For vegetation predictors (Avila-Flores et al. 2012, Baker et al. 2013, Olson et al. 2017, bare ground and shrub cover were summarized from rasters developed by Homer et al. 2012 (30m resolution) and NDVI (Normalized Difference Vegetation Index) was summarized from rasters provided by Merkle et al. 2016 (250-m resolution, 8-day interval). NDVI is an index of plant greenness that is correlated with plant biomass in arid and semi-arid ecosystems (Hobbs 1995). For anthropogenic development predictors (Johnson and Collinge 2004, Magle and Crooks 2009, Sackett et al. 2012, Olson et al. 2017, overall disturbance was based on an index developed by the University of Wyoming's Geographic Information Science Center (30-m resolution, http://nrex. wyo.gov) that attempts to integrate development from multiple sources, such as agriculture, energy development, and roads. Oil and gas summaries were based on point locations of active well pads (Wyoming Oil and Gas Conservation Commission, http://wogcc.state.wy.us). In addition to including these predictors in occupancy models, we summarized predictor values across sites and grouped by ownership stratum to aid with interpretation of results.
Spatial null predictors assessed factors that may be partially confounded with land ownership and create spurious ownership effects. For example, there is more private land in the eastern portion of the WTPD range in Wyoming (Fig. 1), so we included the coordinates for each site's centroid to assess whether occupancy varied along an east-west or north-south gradient that was not an effect of land ownership per se. Similarly, an ownership stratum effect may depend v www.esajournals.org on whether a site is surrounded by the same stratum or a mix of strata. To account for site context, we included the percent of a 2.5-km site buffer that was the same ownership stratum as the site (stratum.evenness; Table 3), based on potential distances of prairie dog dispersal and sensitivity to development (Johnson and Collinge 2004).
We used single-season occupancy models to estimate probabilities of detection (p) and occupancy (Ψ) for WTPD (MacKenzie et al. 2006) as functions of environmental and anthropogenic variables. Multi-model inference was conducted with Akaike's information criterion corrected for small sample sizes (AIC c ), and ΔAIC c and AIC c weight were used to compare the relative support for each model (Burnham and Anderson 2002, Arnold 2010, Cade 2015. When comparing nested models with AIC, it is important to account for models with uninformative predictors that are carried by informative predictors (Arnold 2010, Murtaugh 2014. Models with uninformative predictors take AIC weight that would otherwise go toward higher ranked models, distorting model set interpretation. When comparing nested models that differed by a single predictor, where a predictor could be continuous, a factor, or an interaction, the additional predictor was only considered informative if it resulted in a lower AIC than the smaller model. Models with uninformative predictors were excluded when calculating AIC c weights (see repositories at https://github.com/jceradini for a custom Program R script [AIC_Uninforma-tiveParameters_RMarkOcc] for identifying uninformative predictors and recalculating AIC weights for RMark occupancy models). To minimize the effects of multicollinearity, we used a cutoff of 2.5 for the variance inflation factor (Zuur et al. 2010).
We used a multistage modeling approach (Doherty et al. 2012), which enabled the comparison of a suite of predictor variables. Predictors advanced to the next stage if they were, (1) within two ΔAIC c of the best model, (2) ranked higher than the constant model, and (3) the highest ranked among predictors summarizing the same information but for different time intervals or spatial scales (e.g., spring vs. fall precipitation, or development within a 1-km vs.

2-km buffer).
Stage 1: detection models.-Occupancy was held as the largest additive model that would  Intercept-only (non-spatial) Notes: Predictors were spatially explicit summaries for each site unless noted otherwise. Data sources and literature justifications are listed under occupancy methods.
† Predictor was included as a linear model and a quadratic model.
‡ Temperature during aerial surveys was taken from nearby weather stations.
v www.esajournals.org converge, while detection models competed that were all additive combinations of predictors likely to influence WTPD detection (Table 2); see repositories at https://github.com/jceradini for a custom Program R script [AllAdditiveCom-bos_RMark] for creating RMark formulas for all additive combinations of predictors.
Stage 2: predictor-category occupancy models.-Detection was held as the best model from stage 1 while models within each occupancy predictor category competed (Table 3). Excluding null variables, each occupancy predictor was included twice: once as an interaction with ownership strata and once with an additive strata term.
Stage 3: cross-category occupancy models.-Detection was held as the best model from stage 1 while all additive combinations of the models that advanced from stage 2 competed. Stage 4: model evaluation and final model set.-We conducted four sequential assessments of the models that comprised 90% of the AIC c weight in stage 3 (i.e., the 90% AIC confidence set) in order to identify the best models: Notes: Predictors were spatially explicit summaries for each site. Precipitation and NDVI were summarized for spring (April-June), fall (August-October), and/or the growing season (grow; April-October). Abbreviations: HUC, Hydrologic Unit Code; NDVI, Normalized Difference Vegetation Index; SD, standard deviation. Data sources and references are listed under occupancy methods.
§ Predictor was tested as linear and quadratic only at the 2-km scale, otherwise was linear.
v www.esajournals.org 1. Occupancy effects can be spurious if detection is not properly modeled (MacKenzie et al. 2006). To account for this, we reran the confidence set models with two variations to detection models: (1) using alternative competitive detection models from stage 1, and (2) using the occupancy predictors from the confidence set for both the occupancy and detection components of the model, so detection predictors were replaced with the best occupancy predictors. 2. To retain only necessary interaction terms, we compared confidence set models that had interactions to reduced models excluding those interactions. 3. To test the effect of stratum after accounting for the best predictors, we compared models in the confidence set with and without a stratum term. 4. Finally, we quantified the predictive accuracy of the final confidence set using the area under the curve (AUC) statistic adapted for occupancy models.
Assessing the predictive accuracy of a model is critical when results are intended to inform management decisions. AIC provides a relative measure of prediction error that is useful for model comparison, but does not estimate a model's absolute predictive accuracy Anderson 2002, Aho et al. 2014). Predictive accuracy for presence-absence models is commonly assessed using the area under the curve (AUC) statistic (Hosmer et al. 2013). However, it is problematic to estimate AUC for occupancy models because the true occupancy state is unknown for sites with zero detections, and unoccupied sites can thus not be differentiated from occupied sites where the species was not detected (i.e., false negatives; MacKenzie et al. 2006). Assessing predictive accuracy using the naïve occupancy data (i.e., observed 1s and 0s) is thus counterproductive because it does not account for bias due to imperfect detection and therefore undermines the rationale behind occupancy models (Mackenzie et al. 2006). In order to estimate AUC, we adapted the method of Zipkin et al. 2012 to estimate the latent true occupancy state (z) for sites with zero detections using a parametric bootstrap approach (see Appendix S2 for AUC methods). AUC estimates for occupancy models should be interpreted cautiously since they are dependent on predicting the unknown occupancy state for sites with no detections.
We assessed statistical significance of predictor estimates (β) and odds ratios (exp[β]) based on their 95% confidence intervals (CI) and z-values (β/standard error[β]), with a z-value of approximately ≥|2| indicating a significant effect. Goodness of fit was assessed on the global model using the parametric bootstrap in program PRE-SENCE Bailey 2004, Hines 2006).
To estimate a statewide occupancy rate, we first model-averaged the final stage 4 all-ownerships model set to estimate occupancy probability for each ownership stratum Anderson 2002, Cade 2015). We then multiplied each stratum's occupancy estimate by the number of available sites in the stratum, resulting in the number of sites predicted to be occupied. The predicted number of occupied sites was then summed across strata and divided by the total number of available sites in the sampling frame, resulting in an overall occupancy rate. This process was bootstrapped, resampling sites by stratum, to estimate bias-corrected CI (Manly 1998).
To estimate cumulative detection probabilities, which help to inform future WTPD survey designs, we estimated P* for each site as the probability of detecting at least one WTPD across all survey occasions (Kéry and Schaub 2012). We then summarized P* across sites to estimate overall cumulative detection probability (see Appendix S2 for P* formula).
To quantify how inferences would change if our study were restricted to public land rather than conducted across multiple ownerships, we repeated analysis stages 1-4 with a subset of the complete dataset that was only public land sites. All analyses and plotting were done in program R (Wickham 2009, R Core Team 2017).

Model building and evaluation
We surveyed 440 sites across all-ownership strata (Fig. 1, Table 1). Naïve WTPD occupancy across ownership strata was 0.14; that is, WTPD was observed on 63 of the 440 sites surveyed. Naïve occupancy does not account for detection probability or false negatives, which is why we v www.esajournals.org estimated occupancy rate based on our occupancy models. The additive global model for the all-ownerships analysis passed the goodness-offit test (ĉ = 0.84, P = 0.51).
Stage 2: predictor-category occupancy models.-Stage 2 was used to select the top predictors within each predictor category ( Table 3) that would advance to the additive models in stage 3. For the all-ownerships analysis, the only climate variable that advanced to stage 3 was ppt.spring. 2016, which is the annual mean cumulative precipitation for spring in the year we conducted surveys. Among topography models, the elevation × stratum interaction advanced to stage 3, while the sand × stratum interaction advanced from the substrate predictors. Two vegetation variables, NDVI.grow.ratio 2 and bare, advanced to stage 3. NDVI.grow.ratio 2 is a quadratic version of a ratio between short-term and long-term average NDVI during the growing season, so it represents sitespecific change in NDVI; bare is mean percent bare ground cover. The best model in the anthropogenic development set was stratum-only, so no predictors were sufficiently informative to advance to stage 3 (Appendix S4: Tables S2-S6).
For  : Tables S2-S6). Stage 3: cross-category occupancy models.-Stage 3 was used to compare models that combined the best predictors from each category in stage 2. For the all-ownerships analysis, the top five informative cross-category models comprised 90% of AIC c weight (i.e., the 90% confidence set; Table 4, Appendix S4: Table S7). Every model in the 90% confidence set contained NDVI.grow. ratio 2 and all except one contained bare. The best two models comprised 75% of AIC c weight and were nearly identically ranked.
In the public stratum analysis, the top three models comprised the 90% confidence set (Table 5, Appendix S5: Table S7). Every model in the confidence set contained percent bare ground (bare) and the standard deviation of 2001-2015 growing season NDVI (NDVI.grow.sd.2001(NDVI.grow.sd. .2015. Stage 4: final model set and model evaluation.-Stage 4 evaluated the stage 3 90% AIC confidence set models in order to establish a final model set and assess model predictive accuracy. Model selection results for the all-ownerships analysis did not change with any of the post hoc occupancy and detection model variations. Additionally, confidence set models that included a  Table S7 for complete model set). The italicized models are the final models used for interpretation and prediction. Models that met the uninformative predictor criteria are not shown and were excluded when calculating w. Abbreviations are K, number of parameters; AIC c , Akaike's information criterion corrected for small sample sizes; ΔAIC c , AIC ciminimum AIC c ; w, AIC c model weight; NLL, negative log-likelihood; and AUC, area under the curve (95% CI) based on a parametric bootstrap. Interactions included main effects, but those terms are not shown in the table. For all models, detection probability (p) = survey method + DOY 2 + slope.
v www.esajournals.org stratum term had moderate to strong support relative to reduced models excluding stratum (ΔAIC c for reduced models ≥4.56).
Besides stratum × elevation, no interactions in the 90% confidence set improved model fit relative to a nested model excluding the interaction (Table 4). Additionally, ppt.spring.2016 was weakly supported in one model in the confidence set but was uninformative after controlling for bare ground (Appendix S4: Table S7). The best two stage 3 models were therefore the final allownerships models and were used for interpretation of effects and prediction of occupancy probabilities (Table 4).
In the stage 4 public-only model evaluation, model selection results did not change with any of the alternative post hoc occupancy and detection models. The stage 3 confidence set included different additive combinations of the same four predictors, and the direction and strength of the predictor estimates did not vary meaningfully between models, so we interpreted effects from the top ranked cross-category model; however, we model-averaged all the confidence set models for predicting occupancy probabilities (Table 5).

Model interpretation and prediction
Drivers of WTPD occupancy based on all land ownerships.-Occupancy varied with ownership stratum in the final all-ownerships models; however, the strength and shape of the strata effects differed between the models (Table 4; see Appendix S6 for effect size descriptions for the publiconly analysis). The best model contained a marginally significant interaction between strata and elevation. Predicted occupancy probabilities and odds ratios from the best model both suggested that this interaction modeled important ecological variation (Figs. 2, 3; Appendix S3: Table S1). Specifically, the odds of WTPD occupying private land were significantly higher than public land when elevation was approximately >2040 m; otherwise, the effect was insignificant (Figs. 3, 4; see Appendix S3: Table S1 for pairwise odds ratio comparisons of ownership strata). For example, WTPD was 4.33 times (95% CI = 1.04, 18.04, z = 2.01) more likely to occupy private than public sites at the median study area elevation (2051 m). There was also a narrow range of elevations where WTPD was more likely to occupy checkerboard than public land (approximately 1980-2080 m) and checkerboard than WRIR land (~1490-1780 m; Fig. 4, Appendix S3: Table S1). Finally, although elevation significantly influenced the effect of strata on occupancy, strata did not significantly alter the effect of elevation, which was not significant within any stratum (Fig. 2). The elevation effect was strongest for private sites, where WTPD occupancy was 1.50 times (95% CI = 0.98, 2.27, z = 1.89) greater with every 100 m increase in elevation.
In the second-best model, the effect of ownership stratum did not vary with elevation but rather generally followed a similar pattern as the stratum effects at median elevation in the best model ( Fig. 4; Appendix S3: Table S1). However, in the second-best model both private and checkerboard sites were more likely to be occupied than WRIR sites (OR = 14.07, 95% CI = 1.64, 120.78, z = 2.41 and OR = 5.42, 95% CI = 1.11, 26.35, z = 2.09, respectively), whereas Notes: Model set contained all additive combinations of the models that advanced from each category in stage 2; only the 90% AIC c confidence set is shown (see Appendix S5: Table S7 for complete model set). The top ranked model was used for interpretation, and all three models were used for prediction. Models that met the uninformative predictor criteria are not shown and were excluded when calculating w. Abbreviations are K, number of parameters; AIC c , Akaike's information criterion corrected for small sample sizes; ΔAIC c , AIC ciminimum AIC c ; w, AIC c model weight; NLL, negative log-likelihood; and AUC, area under the curve (95% CI) based on a parametric bootstrap. For all models, detection probability (p) = shrub15. binary + temp + slope.sd 2 .
v www.esajournals.org those effects only occurred within a narrow elevation range in the best model. Occupancy increased significantly with bare ground, and the effects were nearly identical in the best two models. In the best model, the odds of occupancy increased by 1.18 (95% CI = 1.05, 1.32, z = 2.83) for every 5% increase in bare ground (Fig. 5). NDVI.grow.ratio 2 is a quadratic version of average 2015 NDVI as a percentage of average 2001-2014 NDVI, summarized across growing seasons (April-October). This NDVI ratio represents if 2015 average NDVI was less than (<100%), equal to (100%), or greater than (>100%) the 2001-2014 NDVI average, summarized across growing seasons. The odds of occupancy increased as 2015 greenness (NDVI) increased relative to long-term (2001-2014) mean greenness, and occupancy probability was highest where this ratio was highest (NDVI.grow.ratio 2 ; Figs. 5, 6). The effect was similar in the top two models. In the best model, the effect was significantly positive when 2015 greenness was approximately >112% of the long-term mean. For example, when 2015 greenness was 120% of the long-term mean, the odds of occupancy increased by 1.93 (95% CI = 1.36, 2.74, z = 3.68) for every 5% increase in greenness. The effect was insignificant when 2015 greenness was approximately 98-112% of the long-term mean, and negative when 2015 greenness was <98% of the long-term mean. Finally, summaries of site predictor values predominantly varied with stratum and partially with stratum and elevation ( Fig. 7; Appendix S7). Comparison of all-ownerships and public land only results. -There was partial congruence between the best models for the public-only and all-ownerships analyses (Fig. 5,Tables 4,5). Both models contained positive effects of bare ground (bare) and varying effects of different NDVI summaries on WTPD occupancy. In addition, the best public-only occupancy model contained slope, which was not in the best all-ownerships occupancy model; however, slope was in the all-ownerships detection model and the standard deviation of slope (slope.sd 2 ) was in the public-only detection model, so slope likely affected both WTPD occupancy and detection ( Fig. 5; Appendix S3: Fig. S1). In the best public-only model, air temperature during the survey (temp) had a weak negative effect on detection and the best all-ownerships detection model contained a quadratic day-of-year (DOY 2 ) effect. DOY and temp are highly correlated (r = 0.94 for the public-only dataset), so they likely represented similar processes affecting detection (Appendix S3: Fig. S1).
By contrast, the best two public-only occupancy models contained the standard deviation of 2001-2015 fall precipitation (ppt.fall.sd.2001.2015), whereas neither of the two final all-ownerships models contained a precipitation predictor (Tables 4,5). In the best public-only model, WTPD detection was higher on sites with lower shrub cover, although the effect was not significant (shrub.15.binary), but shrub cover did not affect detection in the all-ownerships analysis (Appendix S3: Fig. S1). Finally, the best all-ownerships detection model contained a strong survey method effect, which was absent from the best public-only detection model. The public-only dataset contained ground-only, and aerial and ground surveys, but no aerial-only surveys Fig. 3. Odds of WTPD occupying private land relative to public land sites estimated from the best occupancy model in the final model set, which contained a stratum × elevation interaction. The horizontal gray line is at one, so estimates above the line represent an increase in odds, and vice versa below the line. WTPD was significantly more likely to occupy private than public land when elevation was approximately >2040 m. Solid and dotted lines are the predicted mean and 95% CI, respectively. To facilitate interpretation, UCL estimates above 30 were removed when plotting; however, all mean and LCL estimates are shown. All pairwise comparisons of occupancy odds ratios for ownership strata are in Appendix S3: Table S1. Fig. 4. WTPD occupancy probability by ownership stratum and analysis. Probabilities were estimated by model-averaging predictions from the final models for the all-ownerships analysis and the public stratumonly analysis, respectively. Estimates for the all-ownerships analysis are for the median study area elevation, and all model predictors other than stratum were held at their mean or median for prediction. Points and lines are means and 95% CI, respectively. WRIR is Wind River Indian Reservation. (Table 1), which likely influenced the survey method effect.
Overall occupancy rates and detection probability. -We estimated an overall occupancy rate to provide a baseline for continued WTPD monitoring by model-averaging predicted occupancy across the final models (Table 4). Although not significantly different based on 95% CI, the mean occupancy rate from the all-ownerships analysis was 43% higher than the mean rate from the publiconly analysis. The all-ownerships occupancy rate was 0.20 (95% CI = 0.12, 0.29) across all strata or 70,621 (95% CI = 41,377, 103,512) of the 352,561 available sites predicted occupied. The overall occupancy rate in the public-only analysis was 0.14 (95% CI = 0.09, 0.25) or 29,792 (95% CI = 17,951, 51,590) of the 207,337 available public land sites predicted occupied.
For the all-ownerships analysis, median cumulative detection probability summarized across sites and survey occasions (P*) was 0.90 (95% CI for median = 0.88, 0.92) for public land, 0.42 (0.33, 0.46) for private land, 0.80 (0.50, 0.89) for checkerboard land, 0.87 (0.68, 0.93) for WRIR land, and 0.83 (0.79, 0.87) across all strata. Median cumulative detection Predictions were restricted to observed covariate ranges, and all model variables were held at their mean or median for prediction except for the predictor on the x-axis. For the all-ownerships analysis, predictions were averaged across ownership strata and CI was estimated via bootstrapping; solid and dotted lines are the predicted mean and 95% CI, respectively.

DISCUSSION
Estimates of population vital rates, such as occupancy, and subsequent management decisions can be biased by how study areas are delineated such as limiting surveys to public lands (Wiens 1989, Knight 1999, Spies et al. 2007. We found that WTPD occupancy varied with land ownership and elevation in Wyoming, illustrating the importance of surveying across multiple land ownerships. WTPD occupancy increased with the amount of bare ground on a site, but also when plant biomass (as estimated by NDVI) was higher than a site's long-term average plant biomass. Climate, substrate, anthropogenic development, and other vegetation predictors did not affect occupancy directly, but predictor summaries suggested that differences in occupancy between ownership strata were in part caused by variation in climate and vegetation factors. Finally, we found meaningful differences in occupancy probability and occupancy rate between an analysis based on all land ownership types versus an analysis based only on public lands, enabling us to quantify the bias that would have resulted from restricting our study to public land. Despite the net occupancy differences resulting from land ownership, the underlying habitat factors driving WTPD occupancy were similar across analyses.

Effects of land ownership on white-tailed prairie dog occupancy
We found that WTPD was more likely to occupy sites on private land than public land, with the land ownership effect interacting with elevation in the best model but not in the closely ranked second-best model. Land ownership may be an important driver of wildlife population and community dynamics if factors such as habitat quality, ecosystem management, or land-use vary with ownership (Scott et al. 2001, Spies et al. 2007, Donnelly et al. 2016). We could not directly test what factors caused WTPD occupancy to vary with ownership, but multiple lines of evidence suggest that land ownership effects Fig. 6. Odds ratios of WTPD occupancy for a 5% increase in the quadratic NDVI ratio predictor (NDVI.2015.ratio 2 ; observed range = 87, 166%) estimated from the best model in the final all-ownerships model set. NDVI ratio is average 2015 NDVI as a percentage of average 2001-2014 NDVI, summarized across growing seasons (April-October); the NDVI ratio represents if 2015 average NDVI was less than (<100%), equal to (100%), or greater than (>100%) the 2001-2014 NDVI average, summarized across growing seasons. The horizontal gray line is at one, so estimates above the line represent an increase in odds, and vice versa below the line. The NDVI effect was significantly negative when 2015 NDVI was approximately <98% of the long-term mean, and significantly positive when 2015 NDVI was approximately >112% of the longterm mean. Solid and dotted lines are the predicted mean and 95% CI, respectively. To facilitate interpretation, UCL estimates above 35 were removed when plotting, however, all mean and LCL estimates are shown.
were not spurious and were likely driven by differences in habitat quality. In addition to controlling for detection probability and biotic and abiotic factors, we used spatial null variables (Table 3) to test whether ownership effects were a spurious sampling artifact. We had no support for any of the spatial null models, suggesting that differences between ownerships represented an ecologically meaningful pattern.
Habitat quality and land ownership.-Summaries of our occupancy predictors showed that factors likely to influence WTPD habitat quality, such as precipitation and plant cover, varied significantly with ownership ( Fig. 7; Appendix S7). WTPD is primarily herbivorous, relying on a variety of plants for food and as a main water source (Seglund et al. 2006). Plant quality and quantity, and precipitation through its regulation of plant production (Chapin III et al. 2011), affect prairie dog body condition, overwinter survival, reproductive success, and population dynamics (Hoogland 2001, Yeaton andFlores-Flores 2006, Fig. 7. Growing season NDVI (NDVI.grow) is average NDVI across the 2001-2015 growing seasons (April--October). Growing season precipitation (ppt.grow) is average cumulative precipitation across the 2001-2015 growing seasons. Both variables were predictors in the occupancy analyses. Points in the left column are individual site values plotted against elevation; fitted lines and 95% error bands are from a linear regression model. In the right column, site values were averaged by stratum and error bars are 95% CI. Only the public and private strata are shown in the left column; all strata are in the right column. WRIR is Wind River Indian Reservation. See Appendix S7 for additional predictor summary plots. Facka et al. 2010, Avila-Flores et al. 2012, Grassel et al. 2016.
Summaries of climate and vegetation predictors suggested that WTPD habitat quality was higher on private than public sites, which may partly explain the land ownership effect. Private sites received more precipitation and had more vegetation (less bare ground and higher NDVI; Fig. 7; Appendix S7). Moreover, as elevation increased, the differences between climate and vegetation variables for private and public sites tended to increase, which may explain why private sites were only more likely to be occupied than public sites at higher elevations in the best model (Fig. 7). Greater precipitation and plant biomass on private sites could result in increased forage and plant moisture content, which in turn could increase WTPD fitness (Hoogland 2001, Yeaton and Flores-Flores 2006, Facka et al. 2010) and thus increase the likelihood of private sites being occupied. However, most prairie dog species, including WTPD, require habitat that is a mix of open space and vegetation cover, likely to facilitate predator detection and social interactions (open habitat), and to ensure vegetation for forage and cover (closed habitat; Seglund et al. 2006, Avila-Flores et al. 2012, Baker et al. 2013. So, the effect of increased vegetation on prairie dogs can vary with ecological context (Avila-Flores et al. 2012).
Higher habitat quality on private lands due to greater precipitation and productivity likely reflects broader socio-political factors, because in Wyoming and globally more productive and profitable land tends to be privately owned and underrepresented in conservation (Scott et al. 2001, Hoekstra et al. 2005, Pocewicz et al. 2009). This has ramifications for species conservation, as demonstrated in Oregon, California, and Nevada, where Donnelly et al. 2016 found that private land was <20% of their study area but contained 75% of the mesic resources (e.g., wetlands and riparian areas), which were important for greater sage-grouse (Centrocercus urophasianus) breeding and abundance. If private lands contain a disproportionate amount of high-quality habitat then collaboration with landowners is essential to species conservation (Knight 1999, Scott et al. 2001, Haufler and Kernohan 2009, and our results suggest that such collaboration is important for conserving WTPD. Although conditions may be more suitable for WTPD on private than public sites due to inherent environmental differences, differences in WTPD occupancy could also be due to variation in land management. For example, variation in livestock grazing, agriculture, and prairie dog control between ownership strata could influence WTPD habitat and population dynamics (Pauli andBuskirk 2007, WGFD 2017). Unfortunately, data on private lands in Wyoming and worldwide are generally lacking (Haufler andKernohan 2009, Pocewicz et al. 2009), so we were unable to test these effects. Quantifying these land-use factors in future research could improve our understanding of WTPD population dynamics and better inform conservations efforts.
Plague dynamics and land ownership.-Environmental differences between ownership strata may also differentially influence Sylvatic plague dynamics and thus WTPD occupancy. Sylvatic plague, which is transmitted to prairie dogs by fleas carrying the plague bacterium (Yersinia pestis), is a threat to prairie dogs and a major driver of their population dynamics (George et al. 2013, USFWS 2017. Plague is endemic throughout the WTPD range, with periodic and spatially heterogeneous epizootic outbreaks resulting in largescale die-offs (USFWS 2017), so it is likely present across all our ownership strata. Plague dynamics can be influenced by precipitation, with increased precipitation resulting in decreased prairie dog flea load (Eads et al. 2016). Specifically, increased precipitation can indirectly improve prairie dog body condition by increasing forage availability (Facka et al. 2010, Eads et al. 2016. Prairie dogs in better condition, in turn, have more effective defenses against fleas, such as devoting more time to grooming and exhibiting stronger immune responses. This could be an additional explanation for why WTPD occupancy would be higher on the wetter private lands. However, other studies have shown contradictory effects of precipitation on prairie dog-plague dynamics because increased precipitation can also increase flea abundance, making the net effect of precipitation unclear (Snäll et al. 2008).

Effects of habitat factors on white-tailed prairie dog occupancy
In order to understand what factors drive WTPD population dynamics so we can better inform management (Williams et al. 2002, v www.esajournals.org Morrison et al. 2006), we included a variety of relevant habitat and anthropogenic development predictors in our analysis. The important predictors of WTPD occupancy were similar for all ownerships. This suggests that fundamental aspects of WTPD ecology did not vary with stratum and that our analyses identified biologically meaningful drivers of WTPD populations. Overall, our results indicate that WTPD occupancy varied significantly with two vegetation metrics, bare ground (bare) and NDVI (NDVI.grow.ratio 2 ). We found limited support for climate predictors and no support for substrate or anthropogenic development predictors.
Bare ground and NDVI ratio. -WTPD requires a mix of open habitat, which likely enables predator detection and social interactions, and closed habitat with vegetation for forage and cover (Seglund et al. 2006, Avila-Flores et al. 2012, Baker et al. 2013. We thus expected occupancy models to indicate the importance of both open and closed habitats, for example, through quadratic relationships with vegetation predictors, which is what we found. WTPD occupancy increased with bare ground regardless of ownership (Fig. 5), reflecting WTPD's association with relatively open habitats. WTPD occupancy should not increase indefinitely with bare ground since habitat quality should decline as bare ground approaches 100% due to a lack of food and cover. However, a quadratic bare ground effect was not supported, potentially because our sites only contained a range of bare ground values that have a positive effect on occupancy (Appendix S7).
Regardless of ownership, we found a quadratic effect of the NDVI ratio (NDVI.grow.ratio 2 ); occupancy increased when 2015 plant biomass (NDVI) was approximately >112% of the longterm average and decreased when it was <98% of the long-term average (Fig. 6). NDVI.grow.ratio is a ratio between short-term and long-term average NDVI, so it represents short-term, sitespecific change in NDVI. It captures spatial and temporal variation in NDVI, whereas absolute NDVI only captures spatial variation. Therefore, NDVI.grow.ratio is not correlated with bare ground (r = 0.08) whereas bare ground is negatively correlated with absolute values of NDVI (r range = −0.85, −0.75). A site can have relatively open habitat as evidenced by high percent bare ground, but also have a high NDVI ratio that indicates plant biomass increased in 2015 relative to previous years. We were therefore able to capture WTPD habitat requirements for a mix of open habitat and plant cover, resulting in occupancy models that better reflect WTPD natural history. Olson et al. (2017) found a negative effect of shrub cover and a positive effect of herbaceous cover on WTPD density in Wyoming, reflecting similar WTPD habitat requirements. We did not include herbaceous cover in our analysis because remotely derived herbaceous estimates performed poorly in validation tests (Homer et al. 2012). However, we agree that herbaceous cover is an important component of WTPD habitat (Seglund et al. 2006, Baker et al. 2013).
An NDVI ratio could be valuable for assessing trends in WTPD populations because it synthesizes spatial and temporal variation in plant biomass. First, it can help control for stochastic, inter-annual variation in plant biomass among sites and survey periods. For example, if WTPD occupancy fluctuated across sampling periods and/or sites in response to short-term changes in NDVI, then the unexplained variation may obscure long-term population trends if an NDVI ratio was not included. Second, it can provide evidence for whether long-term trends are being driven by large-scale shifts in climate that are resulting in temporal changes in local habitat conditions. For example, climate change forecasts generally predict warming and drying conditions in much of WTPD range (Dai 2011, USFWS 2017. Based on our analysis, this could result in directional shifts in biomass, which would translate into temporally decreasing NDVI ratios across large portions of WTPD range with a commensurate drop in WTPD occupancy. Such a trend and associated causal link would not be identifiable without inclusion of NDVI ratios.
Effects of anthropogenic development on whitetailed prairie dog occupancy WTPD occupancy did not vary with anthropogenic development predictors at any of the three spatial scales tested. The lack of a development effect suggests that WTPD is partly tolerant of disturbance from development, which has also been shown for black-tailed prairie dogs (BTPD; Cynomys ludovicianus). For example, BTPD colonies maintained genetic connectivity in a landscape fragmented by urbanization and agriculture (Sackett et al. 2012), and BTPD densities increased with urbanization and road density (Johnson and Collinge 2004). By contrast, Olson et al. (2017) found that WTPD density increased as the distance to oil and gas wells increased. Olson et al. (2017) may have captured a more subtle development effect because they assessed WTPD density rather than occupancy (Gaston et al. 2000), and they paired WTPD locations with oil and gas locations, whereas our predictors were summarized to the site scale or greater.
Predictor summaries (Appendix S7) indicated that our sites captured substantial variation in development at the scale of the 1-km and 2-km site buffers, but only marginal development variation at the site scale. Therefore, our results may suggest that WTPD tolerate development at broader scales (1 and 2 km), while no effect at the site scale could be due to lack of power resulting from insufficient variation in development between sites. Future studies should attempt to include greater variation in site-level development to further address this question.

Monitoring and management recommendations
One obstacle to landscape-scale species monitoring is the need to survey large regions across multiple land ownerships. We found that a combination of ground-based and aerial occupancy surveys was effective for surveying WTPD across multiple ownerships. However, ground and aerial detection of WTPD were significantly different, so it is essential to design a study where detection probability can be estimated, for example, through repeat or independent dual-observer surveys, otherwise results will be biased (MacKenzie et al. 2006). This approach also enabled us to estimate detection probability for all-ownership strata. However, private stratum sites were only surveyed for one occasion (Table 1), resulting in substantially lower median cumulative detection probability (median P* for private = 0.42) compared to the public, checkerboard, and WRIR strata sites (P* = 0.90, 0.80, and 0.87, respectively), which were surveyed one, two, or three occasions. Future WTPD monitoring should attempt to conduct repeat surveys on a subset of private stratum sites in order to improve detection and occupancy probability estimates (MacKenzie et al. 2006).
Climate and vegetation factors, primarily recent plant biomass relative to long-term averages as indexed by NDVI ratios, open habitat as indexed by bare ground, and potentially precipitation, were important predictors of WTPD occupancy. Incorporating these factors into the WTPD monitoring framework can help refine trend analyses in the face of stochastic climate variability and directional shifts in climate. Further, in addition to directly monitoring WTPD populations, managers can monitor changes in these temporally and spatially explicit climate and vegetation variables that may index WTPD habitat quality and, in turn, WTPD populations. This will help managers anticipate population changes and prioritize at-risk regions within the WTPD range. This is particularly important for WTPD because rangewide monitoring will likely only occur every 3 or 6 years (USFWS 2017), while remotely sensed indices of habitat quality can be estimated annually. Additional research should be conducted to improve the mechanistic understanding of how these factors are driving WTPD occupancy.
Our full analysis based on all land ownerships and our analysis restricted to only public land both suggest meaningful differences in occupancy probabilities and rates based on land ownership. Since these occupancy metrics are central to WTPD monitoring, we would have mischaracterized the status of WTPD in Wyoming if we had restricted our study to public lands, resulting in inappropriate management recommendations and biasing baseline estimates. For example, our research revealed that effective WTPD conservation will require collaborating with private landowners to prioritize WTPD habitat on higher elevation private land, such as within 2000-2500 m where the odds of occupancy were highest. Finally, our study suggests that landscape-scale flora and fauna monitoring and research should account for land ownership effects and consider how results and management recommendations may be influenced by restricting surveys to public land (Knight 1999, Scott et al. 2001, Haufler and Kernohan 2009). Department (WGFD). Bob Lanka, Zack Walker, and Nichole Bjornlie of the WGFD were essential partners on this project, helping to secure funding and providing valuable feedback on study design and this manuscript. We are grateful for the dedication of many field technicians, especially Jessica Sellers, who conducted WTPD surveys throughout Wyoming, and for WGFD and Bureau of Land Management personnel who assisted with field logistics. We thank Wyoming Natural Diversity Database staff for support throughout the project, and Dr. Trent McDonald for statistical advice. We thank the anonymous reviewers for providing valuable feedback on this manuscript. There are no conflicts of interest between this research, publication, and the authors.