Projecting current and future location , quality , and connectivity of habitat for breeding birds in the Great Basin

We estimated the current location, quality, and connectivity of habitat for 50 species of breeding birds in four mountain ranges in the central Great Basin (Lander, Nye, and Eureka Counties, Nevada) and projected the future location, quality, and connectivity of habitat for these species given different scenarios of climate-induced land-cover change. In the United States, such models are relevant to federally mandated management of wild animals by state-level agencies. We sampled birds during the breeding seasons of 2001–2009 with fixed-radius point counts. For each species, we used boosted regression trees to model incidence (proportion of years a location was surveyed in which the species was present) as a function of topography and current land cover and climate. To assess model fit, we calculated the proportion of binomial deviance explained. We used cross-validation to estimate the predictive accuracy of the models. We applied the conservation planning program Zonation to identify locations where incidences of multiple species were maximized through time given current land cover and two scenarios of land-cover change, expansion of pinyon–juniper woodland into sagebrush shrubsteppe and contraction of riparian woodland. Models based on a set of 13 covariates derived from remotely sensed data had some predictive capacity for 41 of 50 species. Model outputs suggested substantial changes in amount of habitat for many species following projected expansion of pinyon–juniper woodland, but less pronounced changes following projected contraction of riparian woodland. Zonation analyses indicated that the spatial distribution of the highest-quality habitat for the avian assemblage was relatively consistent through time under both scenarios. Breeding birds in the Great Basin commonly are grouped in management plans on the basis of their general association with land-cover classes such as pinyon–juniper woodland, sagebrush shrubsteppe, and riparian woodland. However, even within these groups, the environmental attributes that explained a high proportion of variation in species’ incidences and the projected responses to different scenarios of land-cover change varied considerably among species.


INTRODUCTION
Governmental and nongovernmental entities worldwide charged with managing wild animals (henceforth, wildlife) seek to understand how projected changes in land use, land cover, and climate may affect species' probabilities of persistence.Decision-making about the locations and intensity of human activities (e.g., agriculture, livestock grazing, timber harvest, recreation) or settlements may consider where habitat quality for one or more species is likely to be high or low given multiple scenarios of environmental change.The gradient of habitat quality for one or more species sometimes is framed as a gradient of relative conservation value (e.g., Forest Stewardship Council 2006, Barlow et al. 2010).
In the United States, individual states have primary responsibility for management of wildlife.In 2000, the U.S. Congress enacted the State Wildlife Grants Program [U. S. Code 16 (2000) 669 (c)].To receive funds from this program, in 2005 the 50 states and five U.S. territories each developed a state wildlife action plan to maintain or restore populations of rare or endangered species and to prevent endangerment of species that currently are common (Association of Fish and Wildlife Agencies 2006).Some states are revising their wildlife action plans to more comprehensively respond to climate change and other stressors, such as development of renewable and non-renewable sources of energy, changes in disturbance regimes, and expansion of non-native invasive species.The highestpriority taxa vary among regions and jurisdictions.For example, the Nevada Department of Wildlife wishes to explore how changes in climate may affect the distribution and quality of habitat for breeding birds (Nevada Department of Wildlife 2012).
It is well established that efforts to conserve a species are more likely to be effective when patches of its habitat are not highly isolated from each other (Hanski 1999).The Great Basin and Mojave Desert are among the few ecosystems in the conterminous United States in which movement of wildlife is relatively unimpeded by human infrastructure (Theobald et al. 2012).For this reason, and because the ecosystems have considerable topographic and microclimatic heterogeneity, species in the Great Basin and Mojave Desert may have relatively high potential to persist despite environmental change.Advances in technology and computing power have facilitated development of quantitative methods that can examine the quantity, quality, and configuration of habitat at spatial extents that are most relevant to ecosystem-level management (e.g., in the western United States, Canada, and Australia, hundreds to thousands of km 2 ) (e.g., Urban and Keitt 2001, McRae 2006, McRae et al. 2008).
However, as in virtually all ecosystems, changes in land use and land cover are gradually changing the locations and attributes of habitat for many species in the Great Basin and Mojave Desert.For instance, the Lincoln County [Nevada] Conservation, Recreation, and Development Act of 2004 (H.R. 4593; http://www.gpo.gov/fdsys/pkg/PLAW-108publ424/content-detail.html) provided authority to divert groundwater from the central and eastern Great Basin to the Las Vegas, Nevada metropolitan area.Such diversion also may affect surface flows and the species associated with seeps, springs, and riparian areas, including many endemic taxa.Development of renewable energy and changes in fire dynamics are among the numerous environmental shifts affecting regional habitat quantity, configuration, and quality (e.g., Balch et al. 2013).Of about 240 species of birds that breed in Nevada, 72 were designated as species of conservation priority in the state's 2005 wildlife action plan and about 40 others are identified as priorities by regional or continental partnerships (Wildlife Action Plan Team 2006).
Our objectives were developed in consultation with the Nevada Department of Wildlife and philanthropic organizations that support diverse environmental conservation programs in the Great Basin and Mojave Desert.We aimed to estimate the current location, quality, and connectivity of habitat for 50 species of breeding birds in the central Great Basin, which falls within central Nevada, and to project the future location, quality, and connectivity of habitat given different scenarios of land-cover change (Thomson et al. 2009).As described below, these 50 species were detected at a sufficient proportion of our sampling locations for modeling to be feasible.We focused on species associated with three land-cover types that are high priorities for Nevada Department of Wildlife: woodland dom-inated by single-leaf pinyon (Pinus monophylla) and juniper (Juniperus osteosperma, J. occidentalis), shrubsteppe dominated by sagebrush (Artemisia spp.), and riparian woodland dominated by deciduous trees (e.g., aspen [Populus tremuloides]) and shrubs (e.g., willow [Salix spp.]).We defined habitat for each species by identifying functional relations between empirical data on species occurrence and data on geophysical attributes, land cover, and the composition and structure of vegetation.
We examined two scenarios of land-cover change between 2010 and 2100: expansion of pinyon and juniper woodland and contraction of riparian woodland.Across the Intermountain West, expansion and increases in density of pinyon and juniper, especially into sagebrush steppe, has been observed for decades (Burkhardt andTisdale 1976, Tausch et al. 1981).Multiple natural and anthropogenic factors, not mutually exclusive, have been suggested as mechanisms.For example, beginning in the 1880s, a 60-year period of relatively warm temperatures and high precipitation in the Great Basin may have facilitated woodland expansion.This hypothesis is consistent with pollen records spanning the Holocene (Miller andWigand 1994, Gray et al. 2006).Additionally, intensive and extensive grazing by cattle and sheep than began in the 1870s removed herbaceous vegetation that might compete with pinyon and juniper for limiting resources (Miller and Rose 1995).Furthermore, episodic expansion of pinyon and juniper is consistent with data on historical fire frequencies and dynamics (Bukowski and Baker 2013).
A widespread program to remove pinyon and juniper through prescribed fire and other fuels management treatments has been implemented on the basis of assumptions that doing so will increase habitat quality and quantity for Greater Sage-Grouse (Centrocercus urophasianus), a candidate for listing under the U.S. Endangered Species Act.Occurrence of Greater Sage-Grouse is associated most closely with presence of relatively extensive, contiguous stands of mature sagebrush (Aldridge and Boyce 2007, Hanser et al. 2011, Knick et al. 2013).Preventing the species' listing is among the highest priorities of resource agencies in the Great Basin.There is, however, a temporal gap of 25 to .75 years, depending on elevation and subspecies of sagebrush, between treatment and reestablishment of mature sagebrush (Welch and Criddle 2003, Baker 2011, Bukowski and Baker 2013) that provides the habitat structure needed by Greater Sage-Grouse.Moreover, fire can fragment highquality habitat for other species of breeding birds that are associated with pinyon and juniper woodlands (Gillihan 2006, Fleishman andDobkin 2009).
Dynamics of woodland contraction in areas near permanent or ephemeral water bodies are difficult to project.Nevertheless, observations and projections across the southwestern United States suggest that precipitation is decreasing in frequency and increasing in intensity and variability (Garfin et al. 2013), and the probability of multidecadal drought now exceeds 15% per 50 years (J.T. Overpeck, personal communication).Even if total precipitation changes little, increases in temperature will decrease water availability and likely will intensify human appropriation of surface water and ground water.
We used the conservation planning program Zonation (Moilanen and Kujala 2008), an analytical method and software for spatial prioritization, to quantify the projected probability that a given location will retain connected, high-quality habitat for multiple species given current conditions and the two future scenarios.Our work recognizes that species differ in their environmental associations and responses to environmental change, yet most are unlikely to be managed individually.Projections of species' responses to current and future configurations of land cover across the Intermountain West may inform decisions about management of specialstatus species, land-use permitting or zoning, and where to focus, restoration, fuels treatments, fire suppression, and post-fire rehabilitation.

Field methods
Our study area (about 8,000 km 2 above 1800 m and .28slope) includes much of the adjacent Shoshone Mountains and Toiyabe, Toquima, and Monitor ranges (Lander, Nye, and Eureka Counties).Almost all of the land is managed by federal agencies that have multiple-use mandates, and the state of Nevada has jurisdiction over man-agement of the animals that inhabit the area.
We sampled birds during the breeding seasons (late May through June) of 2001-2009 with 75-m (2001-2004) or 100-m (2005-2009) fixed-radius point counts in canyons that drain the east and west slopes of the mountain ranges.Most point centers were .350m apart.We located points along the full elevational gradient of each canyon, typically with two or three points per 100 vertical m.Points were positioned to sample the dominant land-cover types throughout the canyons.
We recorded coordinates of each point with a global positioning system (GPS).During each visit, we recorded by sound or sight all birds using terrestrial habitat within the point.Visits were conducted in calm weather and within 3.5 hours after dawn.We surveyed each point three times per year for 5 min per visit.All points were surveyed in !2 years, and 50% were surveyed for !5 years; 47 points (13%) were surveyed in 2 years, 18 (5%) in 3 years, 115 (32%) in 4 years, 41 (12%) in 5 years, 21 (6%) in 6 years, 56 (16%) in 7 years, 44 (12%) in 8 years, and 14 (4%) in 9 years.We surveyed a total of 356 points: 60 in the Shoshone Mountains, 152 in the Toiyabe Range, 75 in the Toquima Range, and 69 in the Monitor Range.We considered a species to be present at a point in a given year if it was detected during one or more of the three visits.
In addition to the data from the above 356 points, we used data from 15 variable-width point-count locations (Bibby et al. 2000) within 7.5 km 2 of the Seven Mile study site (Monitor Range, Eureka and Nye Counties, Nevada), part of the Sagebrush Steppe Treatment Evaluation Project (SageSTEP) project.These 15 points were visited in May and June of 2006and 2007(McIver et al. 2010).Each point was surveyed twice per year for 10 min per visit, during which all nonflying individuals were recorded.Points were !600 m apart.Surveys were conducted between 15 min after sunrise and 10:15 A.M. on days without precipitation or high wind.For the analyses presented here, we excluded detections that were .100m from an observer.We considered a species to be present at a point in a given year if it was detected during one or both of the visits.
To characterize vegetation composition and structure (which commonly are associated with habitat quality for birds) in the field, we established plots centered on the 356 long-term bird-survey points (Appendix).

Species distribution models
Our response variable was the incidence, I ij , for each species at each point-count location, where I ij ¼ (number of years species i was detected at point-location j )/(number of years point-location j was surveyed).The incidence may be interpreted as the probability that a species will be detected at a given location in any one breeding season, or as a probability of occurrence during a single breeding season (Thomson et al. 2005).Summing the predicted incidences across all 30 m 3 30 m grid cells in the study area, including locations for which field data do not exist, yields the expected number of occupied cells per year (also referenced as occupied area).For example, if the predicted incidences for each of four cells are 0.5, the expected number of occupied cells is two.We built models for 39 species of breeding birds that occurred in at least 5% of points and had mean I !0.05, plus 11 species of interest that occurred in fewer than 5% of points or had mean I , 0.05 [American Kestrel (Falco sparverius), Rednaped Sapsucker (Sphyrapicus nuchalis), Western Wood-Pewee (Contopus sordidulus), Ash-throated Flycatcher (Myiarchus cinerascens), Loggerhead Shrike (Lanius ludovicianus), Tree Swallow (Tachycineta bicolor), Juniper Titmouse (Baeolophus ridgwayi ), Orange-crowned Warbler (Oreothlypis celata), Virginia's Warbler (Oreothlypis virginiae), Sagebrush Sparrow (Artemisiospiza nevadensis), and Bullock's Oriole (Icterus bullockii )] (Table 1).An additional 55 species that we recorded from 2001 to 2009 were detected too infrequently to model.
We used boosted regression trees (BRTs) (Friedman 2001, Elith et al. 2008), an ensemble modeling method based on classification and regression trees, to model the incidences of each species.The BRTs combine numerous classification and regression-tree models for inference and prediction.Ensemble methods account for uncertainty in model structure, reduce the probability of overfitting, and increase the similarity between model projections and observations (Friedman 2001, Wintle et al. 2003).
We first predicted incidences across the study area as a function of 13 covariates related to topography, current land cover, and current climate that one reasonably might expect to be associated with occurrence of breeding birds (Table 2).We derived these data in a geographic information system (GIS) and converted them to 30-m rasters.The advantage of deriving covariates from remotely sensed data as opposed to field measurements is that values typically can be calculated for the full extent of a given landscape and values of the response variable projected for locations that have not been visited.We resampled data on precipitation and minimum temperature from 800-m to 30-m resolution in ArcGIS (version 9.3,ESRI,Redlands,California).We assigned each 30-m cell the mean values of the 49 cells in the surrounding 210-m neighbor- Note: Mean incidence calculated as [(number of years species i was detected at point-location j )/(number of years pointlocation j was surveyed) averaged among point-count locations] (in brackets).
v www.esajournals.orghood.We chose this neighborhood size because it is similar in extent to the radius sampled by the point counts.
To map riparian woodland, we acquired and post-processed two sets of 21 60-km SPOT-5 multispectral satellite images (four spectral bands, 10-m pixels) that covered the study area.The first set was captured from July through August 2005, when deciduous vegetation was photosynthetically active.The second set was captured during October and November 2007, when deciduous vegetation was dormant.The difference in false-color infrared (composite of near infrared, red, and green color bands) between the two sets of images allowed us to better discriminate between deciduous and evergreen vegetation and thus to identify riparian woodland (Xu et al. 2010).We used a set of five vegetation-structure covariates that were measured in the field (Table 2) along with the 13 topography, land cover, and climate covariates to build a second set of models of incidences.This second set of models allowed us to examine the proportion of deviance explained by vegetation structure and to compare the accuracy of models built with and without field measurements of vegetation.One might expect field measurements to be more strongly associated with distributions of species that remotely sensed covariates, but again values of these covariates cannot be projected as easily for locations that have not been visited.The SageSTEP points were omitted from these analyses because vegetation structure was not measured at those points.
For the 13-and 18-covariate models, we fit 3000 initial models for each species with a boosting algorithm (gbm package in R) (R Development Core Team 2012).We assumed a Gaussian distribution for model errors and weighted points by numbers of years surveyed.Thus, the information each point contributed to the model was a function of the number of years it was surveyed.The weighting reflects an expectation that consistent site-level presence or absence (e.g., present in 4/4 years, I ¼ 1) is more reliably associated with habitat quality than presence or absence in one year (I ¼ 1).We allowed two-way and three-way interactions among covariates (i.e., a maximum interaction depth of the 3 in the gbm function).The optimal number of trees was determined by an out-of-bag estimator, which is based on the reduction in deviance in observations not used to select each new regression tree (Freidman 2001).Final inference (association between covariates and incidences) and prediction were based on the optimal number of trees.
We used the relative-importance measure in the gbm package to evaluate the relative strengths of association between each covariate and each species.Relative-importance measures are based on the number of times a covariate is selected for splitting in a tree, weighted by the improvement in model fit as a result of each split, and averaged over all trees in the ensemble (Friedman 2001, Elith et al. 2008).Values were scaled to sum to 100.
To assess model fit and predictive accuracy, we calculated the proportion of binomial deviance explained, I 2 dev ¼ 1 À deviance (I predicted )/deviance (I constant ), where I predicted is the predicted incidence, I constant represents a null model with constant incidences for all points, and deviance is negative twice the binomial likelihood (Yen et al. 2011).
We used cross-validation to estimate the predictive accuracy of the models.Data from each discrete spatial unit (i.e., an individual canyon or the one SageSTEP unit) were used in turn as test data for models built with data from the remaining units.We combined the test predictions from each cross-validation to produce a single vector of test predictions, which we used to calculate a combined I 2 dev(cv) .Our calculation of incidence assumes that if a species was detected at a point on any visit within a year, the location was occupied in that year.The calculation also assumes no false negatives within years.Because most species are detected imperfectly, the latter assumption is unrealistic.Consequently, there may be some bias in incidences, and in estimated speciesenvironment relations, if detection probability ( p; the probability of detecting the species at a site given the species is present) is related to habitat attributes (MacKenzie et al. 2006).
For six species with relatively high numbers of detections we compared the results of the BRT models with hierarchical Bayes models that incorporated imperfect detection within years (MacKenzie et al. 2006;Appendix).We selected species that were associated with land-cover types that are high priorities for federal and state management agencies: riparian woodland [Warbling Vireo (Vireo gilvus), MacGillivray's Warbler (Geothlypis tolmiei )], pinyon-juniper woodland [Blue-gray Gnatcatcher (Polioptila caerulea), Black-throated Gray Warbler (Setophaga nigrescens)], ecotone between pinyon-juniper woodland and sagebrush shrubsteppe [Green-tailed Towhee (Pipilo chlorurus)], and ecotone between pinyon-juniper woodland and riparian woodland [Spotted Towhee (Pipilo maculatus)].These models estimated detection probability and incidence.The inclusion of imperfect detection in hierarchical models had little effect on estimated species-environment relations, and did not increase explained variation in crossvalidation tests (Appendix).Furthermore, in cross-validation tests the BRT models with Gaussian errors explained more variation than all hierarchical models (regardless of whether detection probability was incorporated; Appendix) and BRT models with binomial errors.Therefore, we report only the results from the BRT models with Gaussian errors.

Occurrence given scenarios of land-cover change
We initially calculated values of all covariates at 30-m or finer resolution.For computational tractability, we averaged values to 150-m resolution (5 3 5, 30-m cells) and used the fitted BRT models to project species' incidences across the study region at that resolution.
We estimated incidences for each bird species given current values of covariates and values associated with expansion of pinyon-juniper woodland and contraction of riparian woodland.We based the first scenario on Bradley's (2010) estimates, which assumed that relative probability of expansion of pinyon-juniper woodland decreases as distance from contemporary woodland increases.The latter assumption was based on observed expansion of woodland from 1986 through 2005.We assumed that in 2100, pinyonjuniper woodland will be present in all 30-m cells in which it is currently present and in all 30-m cells with relative probability of expansion of !20% (Bradley 2010).We calculated the proportion of pinyon-juniper woodland in the 210-m neighborhood to include both the current and expanded cover of woodland.
Because we classified dominant land cover of each pixel, we assumed that if pinyon-juniper woodland expanded into a 30-m cell, it replaced any sagebrush present in that cell.We projected values of the four sagebrush covariates (Table 2) in 2100 by subtracting the proportion of expanded pinyon-juniper from the current proportion of sagebrush.We replaced total canopy cover (Table 2) values in cells to which pinyon-juniper expanded with the canopy-cover values of their nearest neighbors in which pinyon-juniper currently is present.We did not change values for riparian woodland, deciduous canopy in riparian woodland, topography, or climate.Downscaled climate projections for the Great Basin at the resolution of our analyses that were available when the analyses were being conducted (Coupled Model Intercomparison Project Phase 3 [CMIP3] and the Fourth Assessment Report [AR4] of the Intergovernmental Panel on Climate Change) were unlikely to be accurate given the dearth of weather stations in the region and the complexity of the terrain.This scenario effectively doubled the number of cells within the study area in which the dominant land-cover type was pinyon-juniper woodland (Fig. 1).
For the second scenario, we reduced the extent of riparian woodland.We changed 10-m grid cells currently classified as riparian woodland to non-riparian vegetation if current vegetation in any of the four abutting cells was classified as non-riparian (we ignored the four diagonal neighbors in this reclassification).This resulted in a relatively uniform narrowing of riparian corridors throughout the study region.The total area of riparian woodland was approximately halved (Fig. 2), but the total length of riparian woodland changed relatively little (20%).This change in extent was not derived from a rigorous model of land-cover change but is plausible and useful for exploring the feasibility of projecting the future location and quality of habitat given multiple scenarios of land-cover change.We assumed sagebrush and pinyon-juniper woodland replaced the lost riparian woodland.First, for each 210-m neighborhood, we calculated the proportion of 10-m cells from which riparian woodland was lost.Then, to obtain projected values for cover of all sagebrush (Artemisia tridentata tridentata, A. t. wyomingensis, A. t. vaseyana, and A. nova), Wyoming and basin big sagebrush (Artemisia t. wyomingensis, A. t. tridentata), mountain big sagebrush (A.t. vaseyana), black sagebrush (A.nova), and pinyon-juniper woodland in 2100, we increased the current proportion of each land-cover type within each 210-m neighborhood by the proportion of lost riparian woodland (e.g., pinyon-juniper in 2100

Location, quality, and connectivity of habitat
For each species given each land-cover change scenario, we calculated the proportional change in total occupied area (i.e., the sum of predicted incidences for all grid cells in the study area) and the proportional change in landscape contagion.Landscape contagion is a pixel-based measure of the extent to which habitat patches (defined here v www.esajournals.orgv www.esajournals.orgas grid cells with incidences above a given threshold) are aggregated.We used contagion as a spatially extensive measure of structural connectivity.We assessed contagion of grid cells with high predicted incidences (i.e., probable habitat) for each species given each scenario with the contagion index p 2 log( p)/( p À 1), where p is the probability that two adjacent cells have high incidences (Parresol 2011).We defined predicted incidence as high if it exceeded the 75th percentile of the current values (results with a 50% threshold were similar).
Zonation software (Moilanen 2007) is designed to identify areas that are associated with a high incidence for multiple species, accounting for each species' resource requirements and connectivity of those resources.Thus, Zonation seeks connected cells that collectively are associated with occurrence of many species.We applied Zonation to identify the set of cells that maximized summed incidences for multiple species through time given current and future scenarios of land cover (Carroll et al. 2010).
We used the basic ''core-area'' grid-cell removal rule, which ranks cells according to their proximity and contributions to the summed predicted incidences of all species.The Zonation algorithm iteratively removes cells while minimizing reductions in summed incidences (occupied area).For example, if cell A represents 80% of the occupied area for species 1 and 5% of the occupied area for species 2, and cell B represents 50% of the occupied area for species 1 and 50% of the occupied area for species 2, then cell B will be removed before cell A because cell A represents a greater percentage of the putative habitat for one species.Thus, if a location represents a high proportion of the occupied area for any species then the location will receive a high rank, even if incidences of other species are low (assuming species are weighted equally).
We conducted Zonation analyses for the full set of 50 breeding birds and for subsets of species associated with pinyon-juniper woodland, sagebrush shrubsteppe, and riparian woodland (Table 1).Some species occupy more than one land-cover type, but for computational tractability we selected the association we thought to be strongest (Dobkin andWilcox 1986, Mac Nally et al. 2008).
First, we ranked cells on the basis of their proportional contributions to modeled current species distributions and identified cells where breeding birds are expected to have highest incidences given current values of covariates (core locations).We loaded spatial data on modeled incidences for the given set of species (whether all 50 species or the species associated with a given land-cover type) into Zonation.We weighted all species for which models had some predictive capacity in cross-validation tests (I 2 dev(cv ) .0) equally.The species with I 2 dev(cv) ¼ 0 (nine of the 50; Table 3) were given zero weight.Thus, these species did not affect cell rankings, but Zonation still provided information on projected change in occupied area for each species given different scenarios (to the extent models for these species are accurate).We smoothed distributions with the dispersal kernel set to 200 m for all species.We assumed that all species use resources in areas within a 200-m radius of the focal point (i.e., of the midpoint of a 150-m grid cell).Therefore, incidences within each cell were partially adjusted according to values for adjacent cells.As a result, clusters of cells with high incidences were retained and relatively isolated cells were removed.To increase the speed of calculations, we removed 100 cells at each iteration.
We generated three sets of Zonation solutions for each scenario for the full set of 50 breeding birds and for the three subsets of species associated with different land-cover types.For the first solution, we ranked cells on the basis of projected future distributions only.The methods for this analysis were the same as those for the Zonation ranking on the basis of current incidences.
For the second and third solutions, we ranked cells on the basis of both current and future incidences.For the second solution, we weighted current and future incidences equally.This analysis identified locations in which the incidences of multiple species are projected to be high given current or projected future land cover.The analysis assumes that species' dispersal distances are unlimited and that species recognize attributes associated with high habitat quality.For the third solution, we accounted for the proximity of areas with high current and future incidences.We used the species interactions option in Zonation to link current and future incidences for each species.We created two additional data layers for each species: first, current incidences modified by the expected future incidences in cells within 100 m (centroid to centroid) (i.e., connectivity of present habitat to future habitat), and second, future incidences modified by current incidences (connectivity of future habitat to present habitat).

Model building and bootstrap evaluation
Mean incidences for the 50 species of breeding birds ranged from 0.01 to 0.55, with a median of 0.09 and mean of 0.14 (0.14 SD) (Table 1).Models based on the set of 13 covariates derived from remotely sensed data had some predictive capacity (I 2 dev(cv) .0) in cross-validation tests for 41 of 50 species (Table 3).Including field measurements of vegetation structure as covariates did not change model fits (I 2 ) substantially, but reduced predictive accuracy in cross-validation tests for many species (Table 3).Among covariates, proportion of riparian woodland and elevation explained the greatest relative proportion of variation in incidences (i.e., were most strongly associated with occurrence) when field measures were not included as covariates (13.2 6 2.3 [mean 6 SE] and 12.4 6 2.2 for riparian woodland and elevation, respectively) (Table 4).Precipitation (10.9 6 2.9) and elevation (10.3 6 2.3) explained the greatest proportion of variation in incidences when field measures were included as covariates.
There were considerable differences among species with respect to variables that explained a high proportion of variation in probabilities of incidences (Appendix: Table A1, A2, and A3).For example, among the 44 species for which sample sizes were relatively robust, strengths of associ-ation of Mourning Dove (Zenaida macroura), Spotted Towhee, Lark Sparrow (Chondestes grammacus), Dark-eyed Junco (Junco hyemalis), Lazuli Bunting (Passerina amoena), and House Finch (Carpodacus mexicanus) with elevation were greater than the mean plus one standard deviation (remotely sensed covariates only; Appendix: Table A1).Strengths of association of Gray Flycatcher (Empidonax wrightii ), Western Scrub-Jay (Aphelocoma californica), Mountain Chickadee (Poecile gambeli ), Juniper Titmouse, Black-throated Gray Warbler, and Chipping Sparrow (Spizella passerina) with cover of pinyon-juniper woodland also were greater than the mean plus one standard deviation (Appendix: Table A1).

Location, quality, and connectivity of habitat
Zonation analyses indicated that the areas of highest-quality pinyon-juniper woodland habitat (i.e., core areas, or areas associated with high probability of incidence of many species) currently are located in the Toiyabe Range and at higher elevations in the Toquima and Monitor Ranges.Because a high proportion of species associated with riparian woodland occur only in Notes: Measures are based on the number of times a covariate is selected for splitting in a tree, weighted by the improvement in model fit as a result of each split, and averaged over all trees in the ensemble (Friedman 2001, Elith et al. 2008).Values were scaled to sum to 100.Left column, values for models that included 13 covariates derived from remotely sensed data; right column, values for models that included the latter 13 covariates and 5 covariates measured in the field (Table 2).Maximum strengths of association appear in parentheses.Species-specific strengths of association are reported in Supplemental Material.
v www.esajournals.orgthose woodlands, riparian woodlands tended to receive higher ranks than uplands.Riparian woodlands within or adjacent to patches (aggregations of cells) with high incidences for nonriparian associated species had particularly high ranks because Zonation seeks connected cells that collectively are occupied by a high proportion of the regional pool of species.
Projected incidences given a scenario of expansion of pinyon-juniper woodland suggested that area occupied would be substantially reduced (10-30% reduction in number of cells expected to be occupied) for nearly all species associated primarily with sagebrush shrubsteppe (Fig. 3).Area occupied was projected to decrease by nearly 30% for Black-throated Sparrow (Amphispiza bilineata) and Sage Thrasher (Oreoscoptes montanus), and by 10-20% for Rock Wren (Salpinctes obsoletus), Sagebrush Sparrow (Artemisiospiza nevadensis), Vesper Sparrow (Pooecetes gramineus), and Brewer's Sparrow (Spizella breweri ).In contrast, projected area occupied increased from 10% to more than 30% for most species associated primarily with pinyon-juniper woodland (Fig. 3).Area occupied was projected to increase by .25% for Gray Flycatcher, Juniper Titmouse, Plumbeous Vireo (Vireo plumbeus), Chipping Sparrow, Mountain Chickadee, Hairy Woodpecker (Picoides villosus), and Lark Sparrow.Changes in the landscape contagion index given the woodland expansion scenario largely mirrored changes in total occupied area (the correlation between change in occupied area and change in contagion was 0.85).Proportional changes in contagion were greater than proportional changes in occupied area.For example, all species for which we projected .10%decline in occupied area had projected declines in contagion .50%(Appendix: Fig. A1).
Zonation ranks based on the full set of 50 species given expansion of pinyon-juniper woodland were similar to Zonation ranks based on current distributions (Fig. 4), suggesting little change in the spatial distribution of high-quality habitat following woodland expansion.Woodland expansion increased the ranks of some lowelevation areas, especially in the Toquima Range, but decreased the ranks of other low-elevation areas (Fig. 4C).
The riparian-contraction scenario reduced the area projected to be occupied for nearly all riparian-associated species (Fig. 5), although the estimated reductions in total area occupied were ,4% for most species.The largest reductions in area occupied were for Warbling Vireo (10%), House Wren (Troglodytes aedon; 8%), and Rednaped Sapsucker (Sphyrapicus nuchalis; 5%).In these riparian woodlands, Red-naped Sapsucker provides nest cavities for an array of nonexcavating, obligate cavity-nesters (Dobkin andWilcox 1986, Dobkin et al. 1995).Projected area occupied did not increase by more than 1% for any species.The landscape contagion index (connectivity) did not change substantially for any species given the riparian contraction scenario; proportional changes ranged from À0.05 through þ 0.04.
Zonation ranks of habitat quality given contraction of riparian woodland generally were similar to ranks based on current distributions (Fig. 6).Ranks of some areas along the east slope of the Toiyabe Range decreased, but not substantially.Similarly, ranks of some areas in the southern Monitor Range increased slightly, but remained relatively low.

DISCUSSION
Many states currently base wildlife management plans on land cover.That is, they assume suites of species are associated with each given land-cover type and will have essentially identical responses to management of that land-cover type.Our models suggested that despite wellsupported general associations between species of breeding birds and land-cover types, land cover often was not the covariate most strongly associated with occurrence.We projected that expansion of pinyon-juniper woodland would reduce the amount of habitat for species associated primarily with sagebrush shrubsteppe by as much as 30%.Yet there was little change in the spatial distribution of the highest-quality habitat for the full set of 50 species following expansion of pinyon-juniper woodland.Because some of these patches of high-quality habitat fall within designated wilderness areas (Arc Dome in the Toiyabe Range, Mt.Jefferson in the Toquima Range, and Table Mountain in the Monitor Range), there likely is high potential for conservation of these patches to remain a management priority.
A species' probability of persistence is strongly related to its local abundance (Gaston et al. 2000, Zuckerberg et al. 2009), the number of populations of the species, and the geographic distribution of those populations (Lande 1993, Foley 1994).Traditional viability analyses are based on time-series estimates of demographic parameters such as abundance, survival, and reproduction (Beissinger and McCullough 2002).Estimating these parameters requires extensive field surveys and following individual animals.Indirect methods of viability analysis that use area occupied (estimated from presence/absence data) as a measure of a species' geographic distribution within the survey area are a tenable surrogate when intensive demographic analyses are impractical or impossible.Area occupied, the viability state variable, serves as a surrogate measure of the species' abundance in the survey area (Gaston andBlackburn 2003, Borregaard andRahbek 2010).Area occupied is much more feasible to estimate than abundance (Robbins et al. 1989, Zielinski et al. 2005, Pollock 2006).
Differences in the resolution at which vegetation was measured in the field (within, at maximum, 0.53 ha; Appendix) or derived from remotely sensed data (within 4.4 ha) may explain why the accuracy of cross-validation tests decreased when field measurements of vegetation structure were included as covariates.The larger v www.esajournals.orgarea may coincide better than the smaller area with territory sizes or the area within which birds typically seek resources.Reliance on naı ¨ve model fits (I 2 dev ) alone would have led to erroneous inferences about the ability of field-measured variables to explain incidences of most species.Independent cross-validations are essential for estimating predictive accuracy, especially when data are spatially structured.
The apparently small effect of contraction of riparian woodland on birds associated with that land-cover type is inconsistent with our understanding of riparian bird species' strong dependence on these woodlands in the arid western United States (Dobkin and Wilcox 1986, Krueper et al. 2003, Earnst et al. 2012), and may largely relate to the fact that we calculated the proportion of 10-m cells from which riparian woodland was lost within a 210-m neighborhood (about the same spatial extent sampled by our point counts).Although the number of riparian cells was reduced by approximately 50%, the loss was modeled as a uniform contraction in the size of a riparian-woodland patch rather than as loss or fragmentation of patches.The length of the riparian-woodland corridor was not reduced and decreases in connectivity (as measured by contagion) among riparian cells were small.In most cases, lost riparian-woodland cover generally represented ,25% of the total area of 210 m- neighborhoods, and was never .50% of the total area.Our data and models indicate whether species are present in a neighborhood, but do not address abundance, which likely is more sensitive than presence to changes in percentage of riparian cover.
Results for some individual species at least partly reflect habitat associations that are relatively complex and difficult to capture in models that classified pixels on the basis of the dominant cover type rather than percent cover.For example, modeled associations between Greentailed Towhee and pinyon-juniper woodland were relatively weak.Green-tailed Towhee nest in or under sagebrush within open pinyonjuniper woodland (males frequently use trees for song posts), especially on and adjacent to slopes.The species typically does not occur in expansive sagebrush without emergent trees, and it does not occur in mature pinyon-juniper woodland because it depends on shrubs that are absent in denser woodlands (Fleishman andDobkin 2009, Dobbs et al. 2012).Thus, the projected 6% decrease in occupied area (summed incidence) in the pinyon-juniper expansion scenario makes biological sense if expansion increases woodland density and decreases the extent of sagebrush.Additionally, woodland may expand to areas with relatively rough topography in addition to areas with relatively low elevation or slope (Bradley and Fleishman 2008).
Models that associate species distributions with environmental covariates are relatively common, but rarely are based on extensive field data for a full assemblage as opposed to empirical data on a small number of species or bioclimatic envelope models for an assemblage.We used multiple years of empirical data on the presence and absence of breeding birds to generate species-specific estimates of occurrence, relate those estimates to remotely sensed and field measures of topography and land cover, and project how occurrence will shift given alternative scenarios of land-cover change.The methods we applied could accommodate alternative scenarios of change developed by any group of researchers, managers, or stakeholders; are statistically rigorous; and are transferable among taxonomic groups and in space and time.Evaluation of spatial and temporal trade-offs of changes in land cover, whether natural or human-caused, can inform decisions about passive or active management strategies that likely are compatible with achieving multiple biological objectives.

Vegetation sampling
We measured three radial 30-m lines from the center of the point.Lines were separated from each other by 1208.The distal end of each line was the center of a circle with 11.3-m radius (0.04 ha).Within each circle, we recorded identities and sizes (either diameter at breast height or basal diameter, depending on plant morphology) of all live trees and standing dead trees.Vegetation at each plot was measured once from 2001 through 2006; most woody vegetation (as measured here) changed little from 2001 through 2009.
We used a concave spherical densiometer to estimate proportion of canopy cover.To estimate frequency of shrubs and ground vegetation, we used an ocular tube with measurements taken at a 458 angle downward from the line of sight (Noon 1981).We recorded occurrence of approximately 30 dominant tree, shrub, and herbaceous taxa.We collected 21 densiometer and ocular tube readings at each plot: one each at 8 m, 16 m, and 24 m along the 30-m line from the center of the plot toward the perimeter of each circle, and one while facing in each of the four cardinal directions from the center of each circle.We averaged cover values for the canopy for each bird-survey point.We aggregated frequency values for shrubs and ground vegetation, and occurrence data for dominant plant species, into a relative measure of frequency at each point.

Hierarchical Bayes models
Incidences are similar to so-called naı ¨ve estimates of occupancy, which can be calculated for a given site by dividing the total number of detections at that site by the total number of surveys.However, because most species are detected imperfectly, naı ¨ve estimates can be less accurate than those that account for detection probability, the probability of detecting the species at a site given the species is present (MacKenzie et al. 2006).For six species with relatively high numbers of detections, we compared the proportion of deviance explained by the boosted regression tree (BRT) models with the proportion explained by hierarchical Bayes models of Incidence that assumed either perfect or imperfect detection within years (MacKenzie et al. 2006).We selected species that were associated with high priority land-cover types: riparian woodland [Warbling Vireo (Vireo gilvus), Table A1.Strengths of association (relative proportion of deviance explained) between occurrence (measured as incidence, see text) of breeding birds and 13 covariates, derived from remotely sensed data (Table 2), included in boosted regression tree models.

Methods
We used a hierarchical Bayesian model to assess whether failure to account for imperfect detection would bias estimated effects of covariates on incidences.The full model for the detection history of species q at point i in year j, O ij , was In the above model, q ij is the probability of detecting the species in a single visit to site i in year j, which depends on the occupancy status of site i in year j (Y ij ¼ 1 if species present, 0 if species absent) and the single-visit detection probability (p i , the probability of detecting the species in a single visit given it is present).The true occupancy status within a year, which is unknown if O ij ¼ 0, is modeled as a Bernoulli variable with probability p ij , which is equivalent to the Incidence (site-specific probability of occupancy in a single year), and is modeled as a function of point-level covariates X and spatial random effects (canyon and point).
The canyon-level and point-level random effects account for spatial variation not accounted for by covariates, including expected spatial correlations among points within the same canyon.The covariates included linear and quadratic transformations of 11 covariates.We used Bayesian model averaging, implemented via reversible jump Markov chain Monte Carlo (Lunn et al. 2009), to calculate posterior probabilities that linear and quadratic terms had nonzero effects [Pr(b 6 ¼ 0)] and to estimate modelaveraged coefficients, which account for uncertainty in model structure (Wintle et al. 2003).We fitted two models for each of six species.The first model allowed for imperfect detection (p assigned uninformative prior distributions for each site and estimated as part of the model) whereas the second model assumed perfect detection (p ¼ 1 for all sites).We then compared the estimated coefficients, and probabilities of non-zero effects, to determine whether estimated effects of covariates on Incidences were sensitive to assumptions about detection probability.
For each of the six species, we compared the proportion of variation in incidences explained by four models: occupancy model with model averaging, assuming imperfect detection probabilities; occupancy model with model averaging, assuming detection probabilities of 1.0; boosted regression tree (BRT) models that treated incidences as a Gaussian response variable [similar to the models described in the main text, but with a slightly different set of covariates]; and BRT models that treated incidences as a binomial variable.We included the same covariates in all models.The occupancy models did not allow interactions among covariates and used quadrat-ic functions to accommodate nonlinear relations between the response variable and covariates.The BRT models allowed interactions among covariates and higher-order nonlinear relations between the response variable and covariates.
We measured the extent to which each model explained variation in incidences with both a standard r 2 , which treated incidence as a continuous, normally distributed variable; and deviance r 2 , which treated incidence as a binomial variable, thus weighting sites on the basis of the number of years in which they were visited.
We calculated the maximum absolute regression coefficients (maximum from either linear or quadratic terms) for occupancy models for each variable, the relative strengths of association (relative amounts of variation explained) in BRT models for each species, and the posterior probability that variables have non-zero coefficients (i.e., explain substantial variation in the response variable) in occupancy models.Because absolute regression coefficients and relative strengths of association are calculated on different scales, relative values or ranks rather than absolute values should be compared.Posterior probabilities .0.75 generally are considered to be substantial evidence of an association between the covariate and the response variable.However, it is possible to have a high posterior probability for a covariate that is weakly associated with the response variable (i.e., strong evidence of a weak association).Therefore, comparisons between results of occupancy and BRT models are more reliable when coefficients are compared to strengths of association.

Results
Among the four models (hierarchical Bayes model with model averaging, assuming imperfect detection probabilities; hierarchical Bayes model with model averaging, assuming detection probabilities of 1.0; BRT models that treated incidences as a Gaussian response variable; and BRT models that treated incidences as a binomial variable), the BRT models for five of six species explained about equal or greater proportion of deviance in the response variable than the hierarchical models (Table A4).Whether the BRT model with incidence treated as a Gaussian or binomial response explained greater deviance depended on whether the R 2 or deviance R 2 was examined.With the exception of Warbling Vireo, accounting for imperfect detection probabilities generally did not increase the proportion of variance explained, nor the set of covariates strongly associated with occurrence.
The hierarchical and BRT models generally suggested that the same covariates had the greatest strengths of association with the response variable (Table A5).Given that the BRT models allow interactions among covariates and nonlinear relations between response variables and covariates, the results of the hierarchical and BRT models will never be identical.Results from the hierarchical models that assumed imperfect and perfect detection were similar, although the models that assumed imperfect detection often had larger effect sizes (larger coefficients).Accordingly, assumptions about detection probability did not affect which covariates were identified as strongly associated with incidences, but allowing for imperfect detection often led to larger estimates of effect size.Because species distribution maps standardized, Zonation § Occupancy model assuming perfect detection probability.
} Boosted regression tree model that treats incidence as a Gaussian response variable.# Relative strength of association (amount of variation explained).jj Boosted regression tree model that treats incidence as a binomial variable.
Posterior probability that covariates have non-zero coefficients (i.e., explain substantial variance).
v www.esajournals.orgresults are affected by relative effect sizes rather than the absolute effect sizes.Therefore, not accounting for within-year detection probabilities should not have a substantial effect on Zonation results.

Estimation of contagion
We calculated the contagion index as p 2 log( p)/ ( p À 1), where p is the probability that two randomly selected adjacent cells have occupancy probabilities above the 75th percentile of the baseline occupancy probabilities for a given species.p depends on the proportion of cells with high probabilities of occupancy (about 0.25 for observed data; varies for future scenarios) and the conditional probability that an adjacent cell has a high probability of occupancy given that the focal cell has a high probability of occupancy.Results were largely driven by changes in the number of cells with a high probability of occupancy.The degree of aggregation of remaining cells with high probabilities of occupancy was relatively similar given both scenarios of land-cover change (expansion of pinyon-juniper woodland and contraction of riparian woodland) (Fig. A1).

Fig. 1 .
Fig. 1.Shaded relief of study region given a scenario in which pinyon-juniper woodland expands from 2010 through 2100 (Bradley 2010).Mountain ranges from left to right are Shoshone Mountains, Toiyabe Range, Toquima Range, and Monitor Range.Dots represent point-count locations.Uniform light color indicates areas both below 1800 m and with slope ,28.Blue, current distribution of woodland; red, expanded distribution.

Fig. 2 .
Fig. 2. Scenario of contraction of riparian woodland.Areas where riparian woodland was projected to contract are indicated in red.Color intensity is proportional to the proportion of riparian woodland lost within a 210-m neighborhood (brighter colors correspond to loss of a greater proportion of riparian woodland).Dots represent point-count locations.

Fig. 3 .
Fig. 3. Proportional change in summed incidences (occupied area) for species of breeding birds according to a scenario in which pinyon-juniper woodland expands by 2100.Probabilities were summed after predicted incidences were smoothed according to the 100-m dispersal kernel.White bars: species primarily associated with pinyon-juniper woodland; black bars: species primarily associated with sagebrush shrubsteppe; grey bars: species primarily associated with riparian woodland.

Fig. 4 .
Fig. 4. Zonation rankings on the basis of (A) present and (B) future incidences of breeding birds under a scenario of expansion of pinyon-juniper woodland by 2100.Blue areas have the highest predicted incidences for at least one species.(C) change in Zonation ranking from present to future.The relative contribution of blue areas to overall quality and connectivity of habitat is expected to increase given the scenario, whereas the relative contribution of red areas is expected to decrease given the scenario.(D) Zonation rankings on the basis of both current and future probabilities of occurrence and accounting for proximity of future areas with high incidences to areas that currently have high incidences for each species.Blue areas are projected to have consistently high incidences over time.Dots represent point-count locations.

Fig. 5 .
Fig. 5. Proportional change in summed incidences for species of breeding birds according to a scenario in which riparian zones contract by 10 m on all sides.Probabilities were summed after predicted incidences were smoothed according to the 100-m dispersal kernel.White bars: species primarily associated with pinyon-juniper woodland; black bars: species primarily associated with sagebrush shrubsteppe; grey bars: species primarily associated with riparian woodland.

Fig
Fig. A1.Species-specific proportional changes in landscape contagion (degree of aggregation of grid cells with occupancy probabilities .75thpercentile of baseline values) given proportional changes in total occupied area (summed occupancy probabilities) for scenarios of (A) expansion of pinyon-juniper woodland and (B) contraction of riparian woodland.

Table 1 .
Species of breeding birds recorded in the Shoshone Mountains and Toiyabe, Toquima, and Monitor Ranges (Lander, Nye, and Eureka Counties, Nevada, USA), their primary land-cover associations, proportion of point-count locations at which the species was detected during one or more years, and mean incidence.

Table 2 .
Covariates included in models of incidences of individual species of breeding birds.Values for covariates derived from remotely sensed data are means of the 49 30 m 3 30 m cells (210 m 3 210 m neighborhood) surrounding the centroid of the point-count location.

Table 3 .
Model fit (I 2 dev ) and cross-validation (I 2 dev(cv) ) (proportion of deviance explained) values for models with 13 covariates derived from remotely sensed data (n ¼ 356 points) and with the latter covariates plus five covariates measured in the field (n ¼ 371 points) (Table2).The additional 15 points included in the latter models had minimal effect on I 2 dev and I 2 dev(cv) values for those models (i.e., differences in cross-validation results cannot be attributed to the additional points).

Table 4 .
Relative mean strengths of association (relative proportion of deviance in incidences explained by each covariate as compared to other covariates; not absolute proportion of deviance in incidence explained) between occurrence of breeding birds and covariates in all boosted regression tree models.

Table A4 .
Proportion of deviance in incidence explained by models with different assumptions about probability of perfect detection and distribution of the response variable.

Table A5 .
Continued.Maximum absolute regression coefficient (from both linear and quadratic terms).