Integrating spatially explicit indices of abundance and habitat quality: an applied example for greater sage-grouse management.

Predictive species distributional models are a cornerstone of wildlife conservation planning. Constructing such models requires robust underpinning science that integrates formerly disparate data types to achieve effective species management.Greater sage-grouse Centrocercus urophasianus, hereafter 'sage-grouse' populations are declining throughout sagebrush-steppe ecosystems in North America, particularly within the Great Basin, which heightens the need for novel management tools that maximize the use of available information.Herein, we improve upon existing species distribution models by combining information about sage-grouse habitat quality, distribution and abundance from multiple data sources. To measure habitat, we created spatially explicit maps depicting habitat selection indices (HSI) informed by >35 500 independent telemetry locations from >1600 sage-grouse collected over 15 years across much of the Great Basin. These indices were derived from models that accounted for selection at different spatial scales and seasons. A region-wide HSI was calculated using the HSI surfaces modelled for 12 independent subregions and then demarcated into distinct habitat quality classes.We also employed a novel index to describe landscape patterns of sage-grouse abundance and space use (AUI). The AUI is a probabilistic composite of the following: (i) breeding density patterns based on the spatial configuration of breeding leks and associated trends in male attendance; and (ii) year-round patterns of space use indexed by the decreasing probability of use with increasing distance to leks. The continuous AUI surface was then reclassified into two classes representing high and low/no use and abundance. Synthesis and applications. Using the example of sage-grouse, we demonstrate how the joint application of indices of habitat selection, abundance and space use derived from multiple data sources yields a composite map that can guide effective allocation of management intensity across multiple spatial scales. As applied to sage-grouse, the composite map identifies spatially explicit management categories within sagebrush steppe that are most critical to sustaining sage-grouse populations as well as those areas where changes in land use would likely have minimal impact. Importantly, collaborative efforts among stakeholders guide which intersections of habitat selection indices and abundance and space use classes are used to define management categories. Because sage-grouse are an umbrella species, our joint-index modelling approach can help target effective conservation for other sagebrush obligate species and can be readily applied to species in other ecosystems with similar life histories, such as central-placed breeding.

Introduction

The ability to predict the occurrence of imperilled species is a critical component of conservation planning (Johnson & Gillingham 2004; Rushton, Ormerod & Kerby 2004; Early, Anderson & Thomas 2008) because of its importance in preserving biodiversity at local, regional and global scales. Predictive distributional models can delineate priority areas for habitat preservation, identify areas where management and restoration activities are likely to be most beneficial, evaluate risks of anthropogenic activities and identify low‐risk areas where human activities can occur with minimal impact. Each of these objectives may require predictive distributional models that operate at different spatial or temporal scales. For example, models at range‐wide or regional extents are important for assessing species range contractions. Landscape‐scale models may be required to identify specific geographical areas with high potential for preservation. Local‐scale models may be best for evaluating effects of site‐specific developments.

All distributional models require some form of spatial data on the target species (Rushton, Ormerod & Kerby 2004). The utility of such data depends on the match between the spatiotemporal resolution of the distributional information and the scale at which the model is to be applied. Discrete breeding locations and survey counts often are used to define species occupancy at broad spatial scales (e.g. regional), particularly for taxa such as birds that are often surveyed during the breeding season. These types of count data can also provide a useful index as a proxy for abundance (Stephens et al. 2015). Breeding location data sometimes can be relatively coarse, however, if they do not provide spatially explicit information on habitat use other than immediate breeding location. In lekking species, for example, leks (traditional breeding grounds) are located at discrete locations where animals readily can be counted but that do not fully represent the broader area of breeding occupancy. For such central‐place breeders, point‐based models, such as lek location models, are less suited for identifying ecological relationships that predict animal space use at local scales and the associated relationships with habitat features. Such data sets may not be suited for fine‐scale predictions of occupancy that often are necessary for conservation planning.

Advances in radio‐ and satellite‐telemetry allow the location of individual animals to be measured with great precision (Cagnacci et al. 2010). Analytical tools can allow individual‐based location data to be expressed as probabilistic space use by animals at a range of temporal and spatial scales in order to address a variety of applied and ecological questions (Hebblewhite & Merrill 2008). Advances in geographical information systems (GIS) allow predictions of species occurrence based on estimated relationships with underlying ecological features that influence their distribution. In particular, habitat selection indices (HSIs) generated from resource selection functions (RSFs) yield a powerful empirical approach to link animal spatial distribution to environmental characteristics (Manly et al. 2002; Gillies et al. 2006), which can then be combined with geospatial data to predict the relative probability of occurrence across adjacent unsampled areas (Boyce & McDonald 1999; Manly et al. 2002; Johnson et al. 2006). However, conducting telemetry studies with a sufficient spatial extent to avoid over‐extrapolation to unsampled areas is often challenging, especially across remote areas. Modelling approaches that rely on ecological relationships derived from multiple sampled populations within a region could be useful for region‐wide habitat inferences. This could be especially effective in a hierarchical framework, whereby inference is made to the ‘average’ population in a region.

Given the limitations in spatial and temporal data in most ecological studies, the use of indices as proxies for space use, abundance and measures of habitat quality are widespread but often used separately to inform environmental management and policy decisions (Stephens et al. 2015). However, where sufficient data exist, course (e.g. point count) and fine‐scale (e.g. probabilistic) data could be integrated to create an analytical tool that reflects known relationships and projects across broad geographical areas to guide landscape‐scale decisions. Such an integrated tool also could be scaled down to guide decisions at a local scale. For example, the integration of point‐based lek‐survey data and telemetry‐based locational data provide both coarse and fine resolution inputs, respectively, that can be integrated to create an analytical tool across a range of spatial scales.

Here, we describe a process for categorizing region‐wide habitat quality and prioritizing areas for conservation and management of greater sage‐grouse Centrocercus urophasianus; hereinafter ‘sage‐grouse’ across a relatively large portion of their range. Sage‐grouse are a North American lek‐breeding grouse in decline largely as a result of loss, degradation and fragmentation of sagebrush Artemisia spp. ecosystems (Knick & Connelly 2011) and currently occupy approximately one‐half of their historic range (Schroeder et al. 2004). This species has been a candidate for listing under the US Endangered Species Act (ESA) of 1973 (CFR 2010, 2015) and is listed as an Endangered Species in Canada (Stiver 2011). Sage‐grouse are considered an umbrella species (Rowland et al. 2006) based on their dependence on large, contiguous expanses of sagebrush habitat to meet multiple life‐history needs (Knick & Connelly 2011) that concomitantly encompass other sagebrush‐dependent wildlife. In addition, although protection for sage‐grouse and their habitat under the US ESA was recently ruled to be unwarranted, that decision was based in part on enactment of effective conservation measures (CFR 2015). Hence, spatially explicit decision‐support tools, operating at multiple spatial scales, are needed largely because of recent intensification in sage‐grouse population management and policymaking. The process described here incorporates two indices: one derived from large‐scale distributional data based on known lek locations and their counts and the other derived from fine‐scale information of individual sage‐grouse based on telemetry data across 12 subregions. This joint distributional index model can be a powerful decision‐support tool for land and wildlife managers and policymakers.

Materials and methods

Geographical extent of analysis

The geographical extent of our study was defined by the outer perimeter of all sage‐grouse population management units (PMU; NDOW 2014) in Nevada and north‐eastern California. We included a 10‐km buffer beyond the PMU's outer perimeter to ensure adequate representation of available sage‐grouse habitats (Fig. 1). This represented an area of 21·5 million hectares that approximated the total potential sage‐grouse range in Nevada and north‐eastern California, excluding those sage‐grouse associated with the bi‐state distinct population segment on the eastern side of the central Sierra Nevada Mountains (Benedict et al. 2003). Floristically, the region was typical of sagebrush Artemesia spp. shrub‐steppe communities in the Great Basin (see Appendix S1, Supporting Information).

Figure 1

Region‐wide extent and subregion boundaries used in resource selection function analyses for greater sage‐grouse habitat and management category mapping in Nevada and north‐eastern California. MCP, minimum convex polygon.

Sage‐grouse telemetry data

Region‐wide modelling and mapping of sage‐grouse habitat and management scenarios relied on two primary sources of sage‐grouse data: lek surveys and telemetry locations (illustrated in a conceptual model in Fig. 2). A summary of the overall methodology is located in Appendix S2. We used data from 19 sage‐grouse telemetry studies located across the Great Basin (Nevada and California) over 15 years (1998–2013), provided by Nevada Department of Wildlife (S. P. Espinosa) and California Department of Fish and Game (S. C. Gardner). In total, 35 883 telemetry locations of 1612 individual sage‐grouse were compiled into a region‐wide data base. We split location data into three independent data sets: (i) an RSF model training subset that contained 80% of location data; (ii) a classification subset that contained 10%; and (iii) a validation subset that contained 10% (Fig 2). Individual sage‐grouse were randomly assigned to the three data sets at the given proportions, with no individual sage‐grouse occurring across data sets. We assigned individual sage‐grouse to data sets rather than grouse locations to prevent autocorrelation among training, validation and classification data sets attributable to random effects associated with individual grouse. Further details on capture and radiotelemetry techniques can be found in Appendix S3.

Figure 2

Diagram showing conceptual model for greater sage‐grouse habitat selection model and habitat management category map for Nevada (NV) and north‐eastern California (CA). Input data sets (blue boxes) were subjected to a series of processing steps (black boxes) to produce interim and final spatially explicit maps (red parallelograms). RSF, resource selection function.

Delineating subregions

Spatial associations between sage‐grouse and existing PMU boundaries (NDOW 2014) were used as a starting point for delineating subregions for habitat selection analyses (see Fig. S1). We calculated the spatial extent of each subregion (n = 19) from a minimum convex polygon (MCP) that encompassed all telemetry locations, then buffered each MCP by the averaged maximum daily sage‐grouse movement (1451 m; Fig. 1). Seven subregions had insufficient sample size to be included in the training data set (criteria for inclusion were >20 radiomarked grouse and >100 telemetry locations). As another validation set of telemetry data, we used data from these seven subregions to evaluate model performance from sage‐grouse locations in areas not used to train the model. After data screening, 12 subregions were represented in the model training data set.

Classification of landscape habitat features

We used GIS to quantify a broad suite of biotic and abiotic variables that represented hypothesized ecological relationships with sage‐grouse populations. Land‐cover types representing the dominant vegetation within 900 m2 pixels derived from Landsat imagery were classified into binary raster layers (see Appendix S2). The non‐sagebrush land‐cover classes used in our analysis were annual grass, perennial grass, lowland non‐sage shrub, upland non‐sage shrub, wet meadow, riparian, pinyon Pinus spp. and juniper Juniperus spp. conifer (hereafter, pinyon‐juniper conifer), non‐pinyon‐juniper conifer, agricultural cropland, open water and bare ground. All sagebrush taxa were condensed into a single ‘sagebrush’ land‐cover class for analysis.

Because variation in sage‐grouse habitat selection can be strongly scale‐dependent (Doherty et al. 2008; Casazza, Coates & Overton 2011; Aldridge et al. 2012), we measured proportions of each land‐cover type at three different spatial scales relevant to sage‐grouse movement patterns. Specifically, we used a circular moving window in ArcGIS™ Spatial Analyst (Environmental Systems Research Institute, Redlands, CA) with a search radius of 167·9 m (8·7 ha), 439·5 m (61·5 ha) or 1451·7 m (661·4 ha) to classify the proportion of cells within each buffered distance that were classified as each land‐cover type. The search radius lengths are biologically relevant in that they represented the averaged minimum, mean and maximum daily distance travelled by sage‐grouse in this study. Within each spatial extent, we also measured the variety of land‐cover types (i.e. the number of unique land‐cover types), the variety of edge types (i.e. the number of unique combinations of adjacent land‐cover types) and the total length of interface of two adjacent cover types (hereafter, ‘interspersion’).

We also calculated distances to point and linear landscape features and specific topographic characteristics that may affect the probability of sage‐grouse occurrence (Aldridge et al. 2008; see Appendix S2). Landscape distance features included perennial streams, springs, wet meadows and open water bodies. For all features, linear distance was calculated as the Euclidean distance from a used or available point. Nonlinear distance metrics were calculated using exponential decay functions (Nielsen, Cranston & Stenhouse 2009). Topographic characteristics included elevation, topographic roughness (e.g. variance in elevation change; (Riley, DeGloria & Elliott 1999) and topographic position indices (TPI; e.g. finer scale depressions or ridges; De Reu et al. 2013).

Subregional RSF modelling

We estimated population‐level resource selection functions (RSFs) using logistic regression (Boyce & McDonald 1999; Manly et al. 2002; Johnson et al. 2006) in a mixed effects model framework (Gillies et al. 2006), where the landscape habitat features (described above) were modelled as fixed effect covariates. Individual sage‐grouse were fit as a random effect (i.e. random intercept) to balance the unequal sampling effort and account for potential autocorrelation among locations within an individual (Gillies et al. 2006). Year was also included as a random intercept for subregions with >1 year of telemetry data to account for temporal intraclass correlation. Five random locations within each subregion MCP were generated for every used location to measure availability of habitat features (Aldridge et al. 2012). Random locations received less weight to allow for equal influence of used (weight = 1·0) and random (weight = 0·2) points. To allow all seasonal use areas to be represented equally, thereby representing the entire annual cycle of sage‐grouse (i.e. breeding, late‐summer/autumn, wintering), we added an additional weight to each location based on the proportion of use occurring during spring/summer (March–August), autumn (September–November) and winter (December–February). Thus, each season had equal representation for estimation of model parameters. We fit all models using the lme4 package (Bates, Maechler & Bolker 2012) in Program R (R‐Core‐Team 2014).

For each subregional map, we employed a multistep model evaluation approach that has been described in (Coates et al. 2014). The first step reduced the number of variables by identifying the most appropriate spatial scale, distance function or topographic covariate that best approximated the probability of occurrence for each corresponding environmental factor. Covariates that represented the best performing scale function met two criteria: (i) model consisted of the lowest bias‐corrected Akaike Information Criterion (AIC) value (Burnham & Anderson 2002); and (ii) null model (i.e. random effects only model) was >2 AIC relative to the single covariate model. Additionally, covariates that represented the best performing distance and topography function were <2 AIC relative to null model. Variables meeting these criteria were carried forward to step 2. During the second step, we constructed a series of additive models containing all possible two‐covariate combinations of those covariates carried forward from the first step. We sought to reduce multicollinearity by constructing correlation matrices and removing models with evidence of correlated effects (r ≥ ¦0·65¦). We then calculated model‐averaged parameter estimates (βs) for each covariate across the set of additive models to account for model selection uncertainty (Burnham & Anderson 2002). The purpose of this stage was to estimate the effect of each covariate, while accounting for the presence of all other covariates (2‐factor models), and use the model‐averaged parameter estimates to calculate an RSF, rather than developing the most parsimonious additive model with multiple covariates (Coates et al. 2014). The RSF took the form:where w(x) is the RSF as a function of model‐averaged coefficients (β1,…, β) for each covariate (X 1,…, X ; (Manly et al. 2002). Although the RSF is not an absolute probability in our study because unused areas were not known, the RSF is useful as an approximation of the probability of selection (Manly et al. 2002). Covariates were excluded from the RSF if their 95% unconditional confidence intervals overlapped zero.

Region‐wide HSI, classification and validation

The subregional RSF consisted of values spanning orders of magnitude. Therefore, we transformed each RSF to an HSI as follows:

The HSI is equivalent to a logistic transformation on the β X , for each covariate X , and was used only to express relative differences in habitat quality, expressed as a 0‐to‐1 index (low to high) for each grid cell (mapped resolution 900 m2). Although we do not interpret HSI values as absolute probabilities, an increase in HSI corresponds to a relative increase in the probability of selection based on environmental covariates. The subregional HSIs were applied across the region and then averaged across each pixel to calculate a single continuous surface representing the region‐wide HSI (Fig 2).

An index of habitat quality may be obtained directly from the region‐wide surface, which likely provides the most flexible and informative estimates for habitat conditions. However, a continuous 0‐to‐1index at the highest possible mapping resolution may yield greater ambiguity for decision‐making when compared with larger areas classified into discrete habitat quality classes. Therefore, we extracted HSI values from the region‐wide model using the independent classification data set of telemetry locations (e.g. 10% of randomly selected sage‐grouse held aside) and formed four discrete categories of habitat selection (i.e. high, moderate, low and non‐habitat) using classification values corresponding to standard deviations (SD) below the mean () HSI as demarcations.

We validated the accuracy of the categorized HSI map with three independent data sets, two comprised of telemetry locations and one of lek locations (Fig. 2). The first set of telemetry locations consisted of 10% of the data held aside for validation within RSF regions (n = 3124). The second set was comprised of the locations from non‐RSF subregions that had insufficient sample size for HSI calculation and were excluded from the model training data set (n = 609, subregions = 7) and provided evaluation of the habitat classification where they were extended outside the area represented by training data. The third evaluation data set was comprised of all known active leks region‐wide. Locations from each validation set were overlaid onto the classified HSI map then evaluated for agreement between percentages of locations falling within each habitat class and SD percentile class used for habitat classification. We also used Cohen's Kappa coefficient (κ) to assess agreement between the frequencies of observed validated HSI values vs. expected values based on SD percentile bins. Values of κ > 0·75 constitute excellent agreement, 0·40–0·75 are acceptable and <0·40 are poor (Fleiss 1981).

Region‐wide abundance and space use index

We developed a region‐wide index that served as a proxy for abundance and space use (AUI) for sage‐grouse using lek data following procedures described in Farzan et al. (2015). To develop this index, we first obtained spatial coordinates for leks, as well as the number of males occupying each lek over the past 5 years, from data bases compiled by the Nevada Department of Wildlife and California Department of Fish and Wildlife. Using lek data to provide information about contemporary abundance, we combined indices describing both density and distance to account for lek location configuration and variation in counts (density index), as well as space use from individual utilization distributions in relation to lek sites (distance index; Fig 2.).

As described in Farzan et al. (2015), the density index was generated using kernel density estimation (Silverman 1986) on point‐based lek location data to create a utilization distribution weighted by the average maximum number of sage‐grouse counted, following similar procedures described in (Doherty et al. 2011). For further detail on this procedure see Appendix S4. The distance index was developed from the results of Coates et al. (2013) and expressed a nonlinear relationship between distance from leks and the relative probability of occurrence across the annual life cycle of sage‐grouse. For further detail, see Appendix S4. To create the AUI, values for density and distance indices at each grid cell (900 m2) were normalized by dividing by the maximum of their respective index then averaged across the region‐wide extent.

For the purpose of creating a more usable decision‐support map for management and policymakers, we generated two classes from estimated AUI values: high use and low/no use areas. High use areas consisted of regions within the cumulative 85% isopleth of AUI values. (Doherty et al. 2011) delineated high abundance population centres that contain 75% of the known breeding population. However, a more conservative demarcation of 85% was used here for greater spatial connectivity among areas of likely sage‐grouse use and is comparable with previously used standards for sage‐grouse breeding density. Also, the 85% value is appropriate for the western portion of their range because sage‐grouse seasonal use areas often are relatively far apart. Low/no use areas of the landscape consisted of areas with <15% of the cumulative AUI density.

Management scenario categorization

Subsequent intersections between categorized HSI and AUI can provide spatially explicit information to managers and policymakers (Fig. 2). As an example, we developed a rubric, in consultation with a stakeholder team composed of experts from state and federal resource agencies and academia, to derive four management scenario categories from these intersections. A brief rationale follows:

Core Areas: Defined as the intersection between all HSI classes (high, moderate and low) and the high use AUI class. This management scenario is intended to incorporate areas where environmental conditions are relatively good with a high certainty of current sage‐grouse occupancy.

Priority Areas for Management Action: Defined as the intersection between high HSI class within the low/no AUI class, as well as those areas that scored as non‐habitat but occurred within high AUI areas. These are priority areas for management attention and potential action, where high‐quality, but unoccupied, habitat exists, and management intervention might be especially fruitful because those areas have potential for occupancy. Where habitat quality is low, but is sometimes occupied by sage‐grouse, as occurs when sage‐grouse move through low‐quality habitat during seasonal migration, management attention to these counterintuitive locations may be a high priority. Specifically, the priority area scenario encompasses the following: (i) high‐quality habitats based on environmental covariates with a lower potential for occupancy given the current distribution of sage‐grouse; and (ii) sage‐grouse incursion into areas of low‐quality habitat that is potentially important for local populations such as corridors of non‐habitat connecting disjunct higher quality habitats.

General Areas: Defined as moderate and low HSI classes present within the low/no AUI class. This scenario represents areas used less frequently by sage‐grouse yet environmental conditions are likely conducive for grouse.

Non‐habitat Areas: Defined as non‐habitat HSI class within the low/no AUI class. This scenario represents habitat of marginal value to sage‐grouse populations.

Results

Although effects of environmental covariates varied across subregions (see Tables S1 and S2), sage‐grouse consistently selected sagebrush at the largest spatial scale, as expected, and generally avoided pinyon juniper at the intermediate scale across all subregions. A covariate representing distance to water was inversely related to sage‐grouse habitat selection in all subregional RSFs, while the effects of other covariates were variable across subregions, particularly cropland. Habitat selection indices were generally highest in the north‐central part of the region and more variable across the south‐eastern part where alternating mountain ranges and valleys occurred (Fig. 3). When categorized using SDs from the classification data set, high value of selected habitat comprised all HSI values >0·5 SD below the (percentile rank range: 30·9–100%; Fig. 4a). Moderate value comprised HSI values between 1·0 and 0·5 SD below the (percentile rank range of 15·0–30·9%; Fig. 4b). Low value comprised HSI values between 1·5 and 1·0 SD below the (percentile rank range: of 6·7 – 15·0%; Fig. 4c). Non‐habitat was comprised of HSI values ≤1·5 SD below the (<6·7%; Fig. 4d). The ratio of the proportion of habitat area gained to the proportion of RSF telemetry points added, which indicated that demarcations beyond 1·5 SD incorporated disproportionately fewer telemetry points per unit area (Fig. S2). As further support, this demarcation value likely corresponds to the appropriate portion of time (<7%) that sage‐grouse might spend in non‐habitat areas while moving between seasonal habitat or in exploration.

Figure 3

Spatially explicit habitat selection indices for greater sage‐grouse in Nevada and north‐eastern California.

Figure 4

Classification of habitat selection index generated from generalized linear mixed effects models using environmental covariates for greater sage‐grouse in Nevada and north‐eastern California. Index was demarcated into four quality classes: high (a), moderate (b), low (c) and non‐habitat (d).

In overlaying a validation data set on habitat classes (Fig. S3), validation of the habitat classes resulted in relatively good agreement based on both percentage and κ (Table 1). Agreement across all classes was exceptionally strong within the RSF subregion validation set and acceptable for the non‐RSF subregional set. On a cumulative basis, 79%, 94% and 97% of leks occurred within the classes: (i) high only; (ii) high and moderate; and (iii) high, moderate and low habitat, respectively. Acceptable agreement occurred for the lek validation set although more leks occurred in high class habitat and fewer leks occurred in low and non‐habitat than expected (Table 1).

Table 1

Summary of Cohen's Kappa coefficient (κ) to assess agreement between the frequencies of observed validated habitat selection index classes vs. expected values based on standard deviation percentile bins for greater sage‐grouse in Nevada and north‐eastern California. RSF, resource selection function

Habitat selection class	Expected %	Validation set
		RSF subregions % (κ)	Non‐RSF subregions % (κ)	Active Leks % (κ)
High	69	68 (0·97)	56 (0·50)	79 (0·73)
Moderate	15	20 (0·83)	34 (0·37)	15 (0·98)
Low	9	7 (0·89)	3 (0·61)	3 (0·50)
Non‐habitat	7	5 (0·81)	7 (0·85)	3 (0·57)

The highest relative abundances occurred within the north‐eastern portion of Nevada (Fig. S4) and leks tended to be more isolated with lower numbers of males counted at leks in the south‐eastern portion of our study area. High and low/no use areas were demarcated by calculating the 85% isopleth of the AUI (Fig. 5). As a result of allowing the distance index to be a component of the AUI, all active leks were surrounded by some amount of high use area. Management categories were established using a stakeholder‐driven rubric, resulting in a spatially explicit categorical management scenario map (Fig. 6). All active leks were surrounded by some amount of core management, based on overlap between high use and the presence of habitat.

Figure 5

Classification of abundance and use index generated from probabilistic estimates of lek density and nonlinear space use relative to distance to leks for greater sage‐grouse in Nevada and north‐eastern California. Index was demarcated by the 85% isopleth into two classes: high and low to no use.

Figure 6

An example of habitat management categories based on the intersection of objectively classified habitat selection and abundance and use indices for greater sage‐grouse in Nevada and north‐eastern California.

Discussion

We have presented a ‘first of its kind’ spatially explicit model of sage‐grouse distribution that integrates indices of habitat quality (i.e. HSI) with information related to sage‐grouse occupancy and abundance based on sage‐grouse lek locations, population size and movement patterns (i.e. AUI). By combining multiple data types, our joint distribution model incorporated large‐scale distributional information with finer scale probabilistic models of telemetry‐based ecological relationships. Ultimately, this allows for greater flexibility when designing conservation strategies for greater sage‐grouse because decisions can be informed by the same predictive tool at local, landscape and regional scales.

The increased availability of remotely sensed imagery used to classify resources across regional extents, coupled with large telemetry data sets, has led to various types of modelling that improves understanding spatial variation in habitat selection for sage‐grouse (Aldridge et al. 2012) and many other taxa. Previous large‐scale mapping efforts that quantified sage‐grouse habitat have either relied solely upon the spatial distribution of leks and their associated habitat characteristics (Knick, Hanser & Preston 2013) or contrasted historic vs. contemporary correlates of occupancy (Aldridge et al. 2008). Typically, detailed studies of habitat selection using RSF analyses are derived from data on radiomarked sage‐grouse at local scales (e.g. Doherty et al. 2008; Doherty, Naugle & Walker 2010). Only recently a few studies have relied on large‐scale collaborative efforts to develop detailed habitat selection models using multiple telemetry data sets across relatively large spatial extents (e.g. Rice et al. 2013; Fedy et al. 2014). In our study, HSIs were derived from such high resolution data across multiple local sites and then scaled up by averaging predicted values across the region‐wide extents. Furthermore, our management scenarios incorporate contemporary sage‐grouse abundance associated with the distribution and relative density of breeding leks. Because sage‐grouse occupancy of areas is closely associated with the distribution of breeding lek sites (Fedy et al. 2012; Coates et al. 2013), output can then be downscaled more robustly to inform local habitat management decisions, for example by targeting actions in specific areas with high‐quality habitat and high likelihood of sage‐grouse occupancy.

The results of the resource selection models aligned closely with well‐known aspects of sage‐grouse ecology. For example, large contiguous expanses of sagebrush are a critical factor that influence sage‐grouse persistence (Aldridge et al. 2008) and consistent selection of sagebrush at the largest spatial scale (667 ha) across subregions in our study fits this pattern. Avoidance of pinyon‐juniper woodlands, typically at more intermediate (61 ha) scales, corroborates avoidance of this land‐cover type by sage‐grouse even when it is more patchily distributed across the landscape (Casazza, Coates & Overton 2011; Baruch‐Mordo et al. 2013; Knick, Hanser & Preston 2013). Selection for habitats in close proximity to some water sources likely reflects high use of forage habitats with higher soil moisture levels that enhance availability of forbs and arthropods critical to the growth and survival of young sage‐grouse (Connelly et al. 2004; Casazza, Coates & Overton 2011). Selection of cropland in some subregions may indicate similar patterns of resource use, although ecological benefits may be limited since cropland exceeding approximately 25 ha km−2 is associated with decreased sage‐grouse persistence range‐wide (Aldridge et al. 2008).

Validation tests indicated exceptional agreement for HSIs derived within RSF subregions. However, our averaging technique may have resulted in the loss of some resolution in partitioning habitat selection classes across some unsampled areas of the region, as suggested by the imperfect validation agreement within non‐RSF subregions where high and moderate habitat classes were under‐ and over‐represented, respectively. The subsequent categorization of priority and general management scenarios may have some similar bias. However, the extent of the core area scenario is unaffected by this possible bias because the category includes both high and moderate habitat classes and exceptional agreement existed for all habitat vs. non‐habitat in the non‐RSF subregions.

Four potential limitations of our modelling approach merit consideration. First, it is important to recognize that GIS‐based habitat models, such as ours, can miss microscale features that are necessary life‐history‐specific resources (e.g. nesting microsite concealment cover and functional forage types). For example, unmapped microhabitat characteristics likely influence variation in resource selection within high‐quality habitat modelled at broader GIS‐based landscape scales. For management purposes, additional microhabitat information would be beneficial. A second limitation is that using RSFs alone as a proxy for habitat quality does not incorporate potential source–sink dynamics that can result in higher occupancy rates of low‐quality habitats as fewer high‐quality habitats become saturated (Johnson, Seip & Boyce 2004). In other words, the relationships between the quality of habitats, demographic performance and population abundance are often decoupled by multiple ecological processes (Stephens et al. 2015). This limitation is, in part, ameliorated by our assignment of habitat management categories that account for environmental conditions conducive to sage‐grouse with the composite index of abundance and space use derived from lek surveys. For example, sage‐grouse population dynamics are driven, in part, by variation in nest survival (Taylor et al. 2012), and the spatial distribution of sage‐grouse leks correlates strongly to the presence of nesting habitat (Connelly et al. 2004). Importantly, a low AUI does not intrinsically indicate an absence of adequate nesting habitat because factors other than nesting habitat also influence sage‐grouse occupancy. Nevertheless, a relatively high AUI value can serve as a surrogate indicator for presence of nesting habitat. Hence, prioritization rubrics based on indices of habitat coupled with abundance likely can help managers identify areas most relevant to sage‐grouse population performance.

Our method of developing a single composite map that represented the annual life cycle of sage‐grouse did not include separate maps reflecting seasonal habitat (Fedy et al. 2014). Nevertheless, we controlled for season by weighting each location based on the proportion of use occurring during the three seasons (spring/summer, March–August; autumn, September–November and winter, December–February). This technique allowed all seasons to be represented equally and provided a single predictive surface representing annual patterns in selection for use by managers and policymakers.

Lastly, we recognize that effective conservation planning involves stakeholder involvement. We stress that the HSI and AUI demarcations and resulting habitat management categories presented here serve simply as an example of the type of information output that can be created with this empirical framework. Nevertheless, this process provides substantial opportunity for stakeholder perspectives to be modelled for predictive purposes or to pursue customized stakeholder goals, and demarcations can be reclassified readily. For example, SD demarcations can be relaxed or constrained to modify classes of habitat quality, yet objectivity remains as long as new probabilities and rationale are reported.

Our approach represents the only quantitative decision‐support tool currently available to inform region‐ and state‐wide sage‐grouse planning decisions in the south‐western portion of their range. For example, the resulting empirically derived surfaces might be used to help inform decisions related to mitigation, evaluate ecological impacts of management actions and delineate areas of highest importance for protection. Additionally, these maps could be updated and the model re‐employed as new data became available. For example, not all lek locations are known, new leks might arise while others cease to exist, and some leks shift locations. Other updates include variance in resource selection and abundance associated with life‐history‐specific habitat requirements (e.g. nesting, brood rearing and overwinter), anthropogenic disturbances (e.g. power lines, energy development) and wildfire. Adaptive management approaches should be employed to allow for such revisions. Also, our approach is not limited to sage‐grouse and can be applied to other taxa with sufficient data for developing indices as proxies for habitat quality and population abundance, particularly germane for taxa with central‐place breeding or foraging strategies.

Fig. S1. Map showing telemetry points (coloured dots) comprising greater sage‐grouse locations available for use in resource selection function modelling, Nevada and north‐eastern California. Names refer to locations associated with Population Management Units designated by the Nevada Department of Wildlife.

Click here for additional data file.

Fig. S2. Relationship between ratio of added potential habitat to added telemetry points validated and standard deviation (SD) cut‐off.

Click here for additional data file.

Fig. S3. Map showing overlay of radio‐telemetry data and lek locations used to validate habitat selection classes for sage‐grouse in Nevada and north‐eastern California.

Click here for additional data file.

Fig. S4. Abundance and space use index based on lek configuration, lek counts, and seasonal movements in relation to lek sites for sage‐grouse in Nevada and north‐eastern California.

Click here for additional data file.

Appendix S1. Summary description of different imagery source and spectral compositing techniques used to create land‐cover classifications for Nevada and north‐eastern California.

Appendix S2. Summary description for methodology used in the conceptual model illustrated in Fig. 2.

Appendix S3. Summary description of sage‐grouse capture and radio‐telemetry techniques.

Appendix S4. Developing a region‐wide abundance and space use index for sage‐grouse.

Table S1. Summary of effects of environmental covariates in resource selection function models for 12 subregions in Nevada and north‐eastern California.

Table S2. Summary for model averaged parameter estimates of covariates associated best approximating spatial scale or distance metric used to create habitat selection indices across 12 subregions in Nevada and north‐eastern California.

Click here for additional data file.