Objectives

Nevada and northeastern California contain a substantial proportion of the present-day range-wide distribution of sage-grouse (Coates and others, 2021b). To provide information necessary for management decisions, including BLM’s land use planning amendment and associated environmental impact statements (86 FR 66331), our primary objective was to map sage-grouse habitat and produce example habitat management areas in these states. To produce this map, we investigated habitat responses by sage-grouse in Nevada and northeastern California, using a comprehensive research study design, wherein habitat selection and survival patterns were considered and mapped across six distinct annual life stages and seasons: (1) nesting, (2) early brood-rearing, (3) late brood-rearing, (4) spring, (5) summer, and (6) winter. The overarching goal was to establish example habitat management areas that may be used by the BLM and other land managers to guide conservation and management decision-making in conjunction with on-the-ground field data and local expertise. The specific objectives were:

Environmental Landscape Covariates

We quantified landscape conditions potentially associated with sage-grouse habitat using remotely sensed geographic information system (GIS) products. Landscape conditions included vegetation type, hydrologic features, and topography. We considered all landscape variables at multiple spatial scales (Aldridge and others, 2008; Casazza and others, 2011) calculated from sage-grouse movement patterns. We calculated landscape characteristics using a circular moving window (neighborhood analysis tool, ArcGIS Spatial Analyst) within radii representing the minimum (167.9 m; 8.7 hectare [ha]), mean (439.5 m; 61.5 ha), and maximum (1,451.7 m; 661.4 ha), daily distances traveled by sage-grouse primarily during the spring and summer seasons (Coates and others, 2016b). We also included neighborhoods corresponding to 75, 260 and 370 m for nest and brood analyses based on the assumption that movement patterns would be different during these life stages because a female is tied to either a single nest site or a brood that has limited mobility (Dudko and others, 2018). The smallest radius (75 m; 1.8 ha) was intended to accommodate landscape attributes more directly associated with cover characteristics at sage-grouse nest locations.

To characterize shrubland cover types, we used Rangeland Condition Monitoring Assessment and Projection base layers from the National Land Cover Database (RCMAP; Rigge and others, 2020), where each 900-square meter (m²) pixel represented a continuous percentage of cover. Vegetation cover types included annual percent cover of sagebrush, non-sagebrush shrub, and total shrub cover. We used Rangeland Analysis Platform data (RAP; Jones and others, 2018; Allred and others, 2021) to characterize annual estimates of percent cover of bare ground, annual grass, perennial grass, litter, total shrub, and tree cover. We used data from the Sage Grouse Initiative representing mesic areas (Donnelly and others, 2016) to extract alfalfa fields. To characterize burned areas, we created spatiotemporal surfaces of cumulative burned area using Monitoring Trends in Burn Severity (U.S. Geological Survey, 2020) data combined with information on sagebrush recovery rates (Coates and others, 2016c; O’Neil and others, 2020). We included an interaction between annual grass and burned areas in summer selection, brood selection, and brood survival models to account for varying responses to annual grass and recent burns previously documented during key brood-rearing periods (Brussee and others, 2022). Hydrologic feature representations included combined intermittent and perennial stream types, and perennial streams only, which were gathered from the National Hydrography Dataset (U.S. Geological Survey, 2017). We obtained data representing wet meadow and riparian features from the U.S. Fish and Wildlife Service National Wetlands Inventory (https://www.fws.gov/program/national-wetlands-inventory/wetlands-data). We measured distances to hydrologic, wetland, and alfalfa features using Spatial Analyst (ArcGIS 10.6.1). We transformed all distance-based variables so their effects would decay with increasing distances using the formula exp(−d/α), where α represented either the median distance value measured from all locations, or 1,451 m (maximum daily distance traveled), whichever was smaller. For alfalfa, we also considered log-transformed distance to account for the potential for human-modified landscapes to affect sage-grouse beyond 1,451 m, primarily by affecting the predator community (Coates and others, 2016a, 2020b). We also calculated density of linear and polygon features (streams and alfalfa fields) within each moving window scale. Topographic variables derived from digital elevation models (900-m² pixels) included elevation, topographic roughness representing variance in elevation (Riley and others, 1999; Evans and others, 2014), topographic position index (TPI; De Reu and others, 2013), heat load index (HLI; McCune and Keon, 2002), and transformed aspect (TRASP; Roberts and Cooper, 1989). We averaged each of these indices across each circular moving neighborhood to capture scale-dependent responses.

Study Design

To evaluate sage-grouse use of the landscape relative to all accessible habitat, we characterized available habitat in GIS by generating five random locations for each used location for each life stage and season (Northrup and others, 2013). For nests and broods, we calculated the 95th percentile of distance from used locations to the nearest lek (for analysis of nest data) and to nests (for analysis of brood data) separately for each study site. Using these distances, we created buffers around leks and nests within which we then conditioned the available distribution of background locations (Holloran and Anderson, 2005; Coates and others, 2013). We defined spring as March 16–June 30, summer as July 1–October 15, and winter as October 16–March 15. For seasonal models, we restricted available locations to within the 99-percent kernel utilization distribution (UD) using all locations from each study site for a given season to delineate each boundary. We generated five random locations for every used location with the boundary. We extracted values of landscape vegetation characteristics, distance metrics, and topographic indices at used and random locations for all analyses. We centered and standardized continuous landscape variables that represented index values or percent cover and centered exponential decay variables.

Sage-grouse require specific vegetation components during brood-rearing which are functionally different from those required for nesting. When the two habitat types are not adjacent, females must move their broods to meet the requirements of this later reproductive life stage. For these reasons, we classified vegetation characteristics used during the transitional period (defined as 20 days after hatch) as early brood-rearing habitat, whereas those used from days 21 to 51 post hatch were classified as late brood-rearing habitat. We performed our analyses separately for early and late brood-rearing stages (Blomberg and others, 2013; Brussee and others 2022).

For all analyses, we withheld 15 percent of all birds and reserved it to use as testing data for model validation. The data from the remaining 85 percent of birds was used to train models and predict habitat. Model validation techniques are described in the following sections and were used to check model fit and predictive capacity of all final models across seasons and life stages.

Analysis: Variable Selection

In total, we considered 90 landscape variables as candidates in our models for selection and survival. However, we could not fit models to the complete dataset because we measured similar habitat characteristics across multiple spatial scales; strong pairwise correlations existed between some of these covariates (|r|>0.65), and model performance suffers when large numbers of covariates are used to predict a response, resulting in overfitting (Tibshirani and others, 2012; Dormann and others, 2013). To consider all candidate variables equally while accounting for these potential issues, we implemented Bayesian variable and scale selection techniques to evaluate statistical support for each variable and scale and subsequently identify a subset of the most predictive variables to be used for habitat mapping.

Our variable and scale selection models combined two methodological approaches to first identify the most effective variable among a group of similar and typically correlated variables and second to test statistical support for inclusion of those variables when included in models with subsets of other potentially important variables. The first method, known as Bayesian latent indicator scale selection (BLISS; Stuber and others, 2017), was specified to identify the most appropriate scale within a group of candidate scales. Briefly, the BLISS method facilitates estimation of a latent categorical probability for each variable within a defined group. For example, given six measurements of sagebrush cover at six different scales, BLISS estimates a posterior probability (0–1) for each scale. The highest probability can then be interpreted as the scale with the most statistical support (Brussee and others, 2022; O’Neil and others, 2023). The second method, known as Gibbs or indicator variable selection (Kuo and Mallick, 1998; O’Hara and Sillanpää, 2009), involves specifying an indicator dummy variable, termed “w,” to be assigned to each candidate variable (or each group of candidate variables). When fit using Markov Chain Monte Carlo (MCMC) sampling methods, the inclusion probability for each variable is updated, and the support for inclusion can be evaluated from the proportion of MCMC iterations where w=1. Hierarchical modeling using MCMC sampling methods allows BLISS and indicator variable selection to be incorporated within the same model framework, such that the probability of inclusion for all candidate variables could be evaluated from Bayesian posterior distributions. We assigned flat priors for sub-variables within groups (as in, equal prior probability of selection for each scale) and assigned a prior probability of 0.5 for each w, representing a 50/50 initial probability of inclusion. We fit models using MCMC procedures described in proceeding sections (“Analysis: Computation”) and carried forward the variables receiving substantial statistical support when fitting a final predictive model for each season and life stage. Final variables were included if w was greater than 0.50 from the highest supported scale for each candidate group.

Analysis: Models for Habitat Selection

For each life stage and season, we analyzed habitat selection by contrasting used locations with random available locations (Johnson and others, 2006; Northrup and others, 2013) to infer sage-grouse spatial habitat patterns. Specifically, the base model took the form shown in equation 1:

The response observations Y were assumed to follow a Bernoulli distribution, with y=1 indicating a sage-grouse location, and y=0 indicating a randomly sampled background location. After estimation of all model parameters, we computed the resource selection function, $\widehat{w}\left(x\right)$, by discarding the intercepts and applying the exponential link: $\widehat{w}\left(x\right)=\text{exp}\left(X\beta \right)$ (Johnson and others, 2006; McDonald, 2013; Northrup and others, 2013). We validated our models by applying resource selection function (RSF) cross-validation methods to the independent testing data subset. We interpreted each model’s ability to make predictions proportional to the probability of selection from linear model fit statistics (β, R², Spearman’s rank coefficient ρ), comparing the proportion of observations in the testing dataset to those expected based on ordinal RSF bins from the fitted model (Johnson and others, 2006; Fieberg and others, 2018). Full model specification is reported in appendix 1.

Analysis: Models for Survival

We used a hierarchical logistic exposure model (Rotella and others, 2004; Shaffer, 2004; Catlin and others, 2015) to estimate survival across all life stages and evaluate habitat covariate effects. The logistic exposure model assumes a daily survival rate (DSR) that follows a Bernoulli distribution, wherein each day the event of interest represented whether the nest or brood survived. Exposure time is accounted for through the length of each encounter history. For each life stage, we created encounter histories for each bird that consisted of known alive and censored intervals throughout the season. For nests and broods, encounter histories ended with either a successful or failed nest (or brood), and exposure time was calculated as the length of time between visits. To model environmental effects, we considered only the landscape covariates that were determined to be important in the variable and scale selection process. For brood models, we created two time-varying categorical brood stage variables representing early versus late. Because it was not always possible to perform counts exactly 10 days apart, when time intervals overlapped both periods, we split the encounter history and divided the exposure time to correspond to each brood age. When broods failed during an interval that overlapped 21 days for brood age, we allowed the covariate information to inform effects for early and late broods (Brussee and others, 2022). For example, for early broods, the categorical brood stage multiplier variable was 0 for late brood intervals, thus only covariates determined to be important for early brood survival varied, whereas variables determined to be important for survival of late broods were held constant at 0. Similar to the base model for habitat selection, DSR was modeled through the logit link function shown in equation 2:

Observation data were indexed for nest or brood (h) with a row for each interval (i) where $y_{h,i}$ ~ $Bernoulli\left({\theta }_{h,i}\right)$, and where $y_{h,i}=1$ if nest/brood h survived interval i, $y_{h,i}=$ 0 if the nest/brood failed, and ${\theta }_{h,i}$ is the probability of nest/brood h surviving interval i: ${\theta }_{h,i}=DS{R}^{t_{h,i}}$. In this expression, ${\theta }_{h,i}$ represents the daily survival probability, DSR, raised to the length t of interval i for nest/brood h.

We validated our models by comparing observed encounter histories and fates in the testing data to simulate “replicate” encounter histories generated from final model posterior predictive distributions (Schmidt and others, 2010). We did post-hoc 1,000 simulations to calculate a Bayesian predictive P-value (Gelman and others, 2013), where values approaching 0 or 1 indicate poor model fit, and values closer to 0.5 indicate good fit. Full model specifications are reported in appendix 2.

Analysis: Computation

We used NIMBLE (version 1.0.1; de Valpine and others, 2017, 2022) in R (version 4.1.3; R Development Core Team, 2022) to estimate all models within a Bayesian hierarchical framework using MCMC sampling. To prevent overestimation of habitat effects and optimize predictions from each model, we implemented L-1 regularization (Tibshirani and others, 2012; Gerber and Northrup, 2020) by specifying Lasso (as in, Laplace, or double-exponential) prior distributions for all coefficients β with an uninformative hyperprior to be estimated for the tuning parameter λ (Park and Casella, 2008; Hooten and Hobbs, 2015). We included the best scale for all important variables identified from Bayesian variable selection process such that strongly correlated variables (r>0.65) did not co-occur in models. If strongly correlated variables were selected during the Bayesian variable selection process, we included the variable with the highest inclusion probability. All random intercepts were assumed to follow a Gaussian distribution centered on zero, such that their effects could be interpreted as deviations away from the population’s grand mean (Gelman and Hill, 2007; Kéry, 2010; Kéry and Royle, 2015). Using MCMC methods, we ran 3 chains of 30,000 iterations, after a burn-in of 20,000 iterations, and retained every fifth sample. We verified chain convergence visually and using the Gelman-Rubin statistic ($\widehat{r}$<1.1). For parameter estimates, we reported median values and 95-percent credible intervals of the posterior distribution. We used high-performance computing resources (Falgout and Gordon, 2022) when necessary to facilitate model runs of adequate length with large datasets.

Space Use Index

We created an abundance and space use index (ASUI) using peak male sage-grouse abundances at leks from the period between the most recent population nadir for the Great Basin region (2013) and the most recent counts (2021) in combination with a surface that represented distance to leks. The lek abundances were used as the input for kernel point density models (Doherty and others, 2016) using the Kernel Density (KDE) tool in ArcMap. We created a surface representing distance to lek to a maximum of 30 kilometers (km) and reclassified that surface using non-linear curves previously derived for sage-grouse, where values of pixels proximal to leks are weighted much more heavily than remote pixels (Coates and others, 2013, 2016b). The KDE and exponential decay rasters were then normalized to have values between 0 and 1. The normalized surfaces were averaged together to produce the ASUI. From the ASUI, we calculated isopleths at every 5-percent probability contour and considered the region contained within the 85-percent isopleth to be high use and areas outside of that isopleth to be low use (Coates and others, 2016b). Lastly, previous research indicates that sage-grouse use is concentrated within 5 km of leks (Coates and others, 2013), so we buffered all leks by 5 km and merged those buffers with the high ASUI so that remote leks were not underrepresented.

Habitat Mapping

We developed spatially explicit index maps of habitat to guide conservation decision-making throughout the study region, using methods adapted from previous research (Coates and others, 2016b, 2020a; Doherty and others, 2016; O’Neil and others, 2020; Brussee and others, 2022). For habitat mapping, we used the median value of the posterior distribution of slope parameters to represent each covariate effect and applied these in the model equations, where the matrix X included the relevant standardized raster values for each covariate representation at each 900-m² pixel across the study area. For selection maps, we transformed estimates using a habitat selection index (HSI) as $HSI=\frac{w\left(x\right)}{1+w\left(x\right)}$, which indicates relative habitat use proportional to availability on a scale of 0–1 (Coates and others, 2016b). For survival maps, and to obtain cumulative survival, we exponentiated values of DSR based on the number of days in the period (nesting=38 days, early brood-rearing=21 days, late brood-rearing=28 days). We produced nine initial maps, each representing predictions of habitat selection and expected survival on a continuous scale for an associated life stage (nest, early brood, late brood) or season (spring, summer, and winter).

For management purposes, we used the continuous selection and survival layers to calculate a habitat suitability index which we combined with the abundance and space use index to develop a single metric that characterizes habitat selection and survival with sage-grouse occupancy. To create the habitat suitability index, we first calculated seasonal selection and survival indices, where the spring index was a combination of seasonal spring, nest, and early brood-rearing layers, the summer index included seasonal summer and late brood-rearing layers, and the winter selection index only included the seasonal winter selection layer. To calculate the seasonal selection indices, we relativized each layer by dividing by its maximum value, added the individual selection components together and then scaled the seasonal selection indices between 0 and 1. To create seasonal survival indices, we followed the same procedure but only created seasonal survival indices for spring, which included nest and early brood-rearing survival layers, and summer, which only included the late brood-rearing survival layer. To create the habitat suitability index, we added the seasonal selection and survival indices together and relativized them to create seasonal suitability indices for spring, summer, and winter. We then added the three seasonal suitability layers together and scaled the resulting layer between 0 and 1 to calculate the annual habitat suitability index. Finally, we added the annual habitat suitability index to the abundance and space use index, which was multiplied by 0.5 to give greater weight to habitat suitability than occupancy.

To aid interpretation for conservation and management purposes, we categorized individual habitat selection and survival maps. Because output from RSFs represent the relative probability of selection, we first categorized each HSI within mesic and xeric hydrographic regions (Brussee and others, 2023) into four categories: (1) non-habitat, (2) low, (3) moderate, and (4) high (Coates and others, 2016b, 2020a). Categorizing by hydrographic region prevented areal predictions of habitat from being disproportionately skewed in one region relative to the other. To delineate the habitat selection categories, we set three cut-points from the observed location data by updating locations spatially with underlying RSF scores from each selection map and identifying the 5th, 25th, and 50th RSF percentiles. Resource selection function values lower than the 5th percentile were considered non-habitat, those between the 5th and 25th percentile were assigned low, values between the 25th and 50th values were assigned moderate, and values above the 50th percentile were assigned high.

To delineate the survival categories (very low, low, moderate, and high) for each life stage within hydrographic regions, we first updated locations with underlying predicted cumulative survival values from the logistic exposure model for each life stage by hydrographic region. We then used a kernel density estimator to approximate the probability distribution of survival at all failed and successful locations. We assigned the highest category as values exceeding the 75th percentile of the distribution of successful nests or broods. We assigned the lowest category as values less than the 25th percentile of the distribution of failed nests or broods. The final cut-point was calculated from the mean of the joint distribution of failed and successful nests or broods to distinguish between high versus moderate (above the mean value but lower than the 75th percentile of successful) and low (below the mean value but higher than the 25th percentile of failed) versus very low. Importantly, delineating categories separately by hydrographic region does not imply differences in survival responses to environmental covariates. The purpose of this procedure was to make sure that priority habitat areas were not disproportionately mapped within sage-grouse distributional extent. For example, in the Great Basin region, sagebrush plant community type, herbaceous cover, and grass and forb height depend largely on soil type, moisture, and temperature characteristics. In more mesic regions, vegetation characteristics on average are more conducive to sage-grouse habitat suitability, relative to more xeric regions (Brussee and others, 2023), yet sage-grouse occupy substantial areas in both hydrographic types. Delineating habitat categories by hydrographic region helped to balance the amount of area assigned to highly ranked habitat categories more proportionally across occupied areas.

We then created a source-sink index map by combining the categories from the individual habitat selection and survival maps for each reproductive life stage. Specifically, we categorized all habitat such that the highest rank was assigned to pixels with the highest selection and survival, representing high quality source habitats. In contrast, the lowest ranks were assigned to pixels that received high selection but very low survival because these pixels represented maladaptive habitat selection and were most likely to contribute to ecological traps or habitat sinks. Source habitats were defined as any pixel that supported high selection and high survival for a given life stage. Importantly, a given pixel was only considered source habitat if it was not sink habitat in another life stage to avoid prioritizing selected areas supporting very low survival.

We combined all six individual habitat selection maps (nesting, early brood-rearing, late brood-rearing, spring, summer, and winter) to create a composite map. Because habitat requirements can differ across seasons and life stages, we prioritized all periods equally and for each pixel we selected the highest value (high, moderate, low, or non-habitat) across seasons and life stages. This allowed important habitat within any period to be prioritized regardless of its value during other seasons or life stages. Then, we used the composite habitat selection map and the identified source habitat in combination with the abundance and space use index described earlier to delineate specific example habitat management categories, with the following procedure and criteria:

To avoid prioritizing reproductive sink habitat outside of occupied areas, we reduced a pixel’s value by one level if that pixel was considered sink habitat (high selection coupled with very low survival) in any of the three reproductive life stages within the low ASUI class (that is, if a pixel was considered general habitat it was reduced to other habitat or if it was considered other habitat it was reduced to non-habitat). We did not apply this adjustment within the high ASUI class because all selected habitat within occupied areas was determined to merit designation as priority.

Habitat Mapping: Adjustments and Modifications

Finally, we made additional adjustments to the example habitat management categories to incorporate two phenomena that could not be integrated directly into habitat models and source-sink mapping, yet had clear implications for the landscape’s potential to support current and future sage-grouse population stability and persistence: (1) pre-fire conditions in fire scars that recently burned and still had the potential to recover (less than 16 years old; Coates and others, 2016c) and (2) connectivity between important nesting and brood-rearing habitats. We describe these modifications in the next sections.

Habitat Recovery Potential in Recently Burned Areas

Modern remotely sensed land cover data products now represent vegetation cover estimates on an annual time step, so wildfires are represented in these data through their effect on vegetation cover from one year to the next. Although a vast improvement over static land cover data, an implication is that habitat mapping efforts based on the most current conditions do not inherently account for vegetation regrowth after fire and may lead areas with high potential for recovery to be defined as poor quality habitat or not habitat at all from lack of sagebrush or shrub cover. Therefore, to capture pre-fire conditions, we recreated seasonal and life stage selection maps within the perimeter of fires that occurred in the last 16 years (Coates and others, 2016c) using temporally varying spatial layers (RAP or RCMAP) from the year before the fire. We then followed the steps outlined earlier, by taking the highest value for each pixel across all seasons and life stages and combining that layer with the ASUI to delineate pre-fire priority habitat. To allow for recovery potential within occupied habitat, we then included any pixels that were considered priority (PHMA or PHMA+) within the high ASUI class based on the pre-fire layers as priority in the final map of example habitat management categories. This ensured burned areas that originally had high-quality habitat and are currently occupied but still have the potential to recover (burned less than 16 years ago) were prioritized.

Movement and Corridor Analysis

A second potential limitation of the methods described to delineate example habitat management categories is that although they may successfully identify areas that are important for sage-grouse for a given life stage, movement among areas is likely underrepresented. As such, it is possible that the areas categorized as PHMA or PHMA+ may form a fragmented pattern on the landscape. Although the areas between fragments may not be used frequently for specific life stages, they may need to be traversed to access and move between patches of high-quality habitat. Connectivity is particularly important when the fragmented patches are used at different life stages (for example, nesting to brood-rearing) because it may be necessary for a sage-grouse to travel between patches that are used in different ways in different seasons. Therefore, there may be areas between priority habitat that, although not classified as priority habitat by the above methodology, are corridors for traveling between habitat patches, with important conservation value for sage-grouse. To incorporate connectivity, we used a least-cost paths (LCP) approach to identify areas that may serve as corridors between areas of high abundance and areas that were identified as source areas. The LCP algorithm originated in graph theory for finding the LCP across weighted edges in a graph (Dijkstra and others, 1959) and is easily extended to finding LCPs on a raster representing movement resistance; that is, a raster where each cell contains a value indicating the difficulty of movement through that cell (Huber and Church, 1985; Van Bemmelen and others, 1993; Etherington, 2016). Given a start and end point, the LCP algorithm finds the path that accumulates the lowest total cost. This is a frequently used method in ecology for evaluating connectivity on landscapes (Adriaensen and others, 2003; Sawyer and others, 2011; Etherington, 2016).

The movement resistance surface used in connectivity analyses is critical; a resistance surface that does not accurately represent a species’ movement behavior will likely yield incorrect results. Previous studies have used a variety of data sources to develop resistance surfaces, including expert opinion, presence and absence data, and movement path data. Of these data sources, movement path data is generally preferred because it explicitly represents an animal’s movement (Zeller and others, 2012). Therefore, we created a resistance surface using movement data from GPS-tracked sage grouse. To capture movement between nesting and brood-rearing habitat, we first fit an integrated step-selection model using the GPS data from female sage-grouse that had successful nests and subsequently moved up in elevation to access important brood-rearing habitat (Avgar and others, 2016). We used the “amt” package (Signer and others, 2019) in R to format location data into steps and randomly sample three available steps for each used step, with step lengths and turning angles for available steps drawn from a Gamma and a von Mise’s distribution, respectively (Avgar and others, 2016; Signer and others, 2019). We included the log-transformed step length and the cosine of the turning angle in the model to allow for unbiased inferences regarding habitat selection and movement and evaluated habitat variables at the end of each step (Avgar and others, 2016). We included variables outlined earlier that represented vegetation type, hydrologic features, and topography. Because movement models evaluate selection at a smaller scale than resource selection functions, we measured all variables only within the 75-m radial buffer to balance small-scale movement with accuracy of the underlying remotely sensed layers. We used the “inla” package (Rue and others, 2009) in R to fit conditional Poisson models with stratum-specific intercepts (for example, step selection functions), where strata consisted of matched used and available steps and were modeled as random effects with large, fixed variance (Muff and others, 2020).

We then used the model output to create a predictive surface across the study area following the methods outlined earlier for creating a map from a habitat selection model. We scaled this surface to be between 0.01 and 1, where 1 represents the highest possible movement resistance. We used 0.01 as the lowest value, because using zero as the lowest value in LCP analysis can result in unrealistic paths that wander through a patch of zero-resistance cells.

Because the LCP algorithm operates on pairs of points, it was necessary to create sets of points to use as start and end locations. Sage-grouse move from nesting habitat around leks into higher elevation wet meadows to support their broods. To create a set of starting locations, we multiplied the previously created nest selection surface by the ASUI to create a probability surface that captured areas of high abundance and use and high nesting potential. We then sampled 10,000 points from this probability surface with the restriction that all sampled points must fall in a polygon created by taking the top 85th percentile of the ASUI surface and then buffering the result by 5 km. Doing this ensured that all start points would fall within areas classified as suitable sage grouse habitat. For the destination points, we used areas identified as brood source habitat. To avoid finding paths to very small patches of brood source habitat, we removed all patches that were made up of less than 10 cells (9,000 m²). This decision was more based on optimization of model output than on biological mechanisms. Having many very small patches was not considered biologically realistic for corridor analysis; ultimately this decision point only led to the removal of less than 1 percent of all patches and, as such, the decision turned out to be fairly inconsequential. For each remaining patch, we identified the centroid of each 900-m² cell that formed the border of the patch and used these points as the destination points.

To identify paths between the start and end points, we iterated over the set of 10,000 start points. For each point, we first identified all destination points (brood source locations) within 15.3 km, which was the maximum distance from leks for all used locations during the late brood-rearing period. Ideal brood-rearing habitat is generally at higher elevations where moisture is retained into the summer. Therefore, we removed any destination points that were at a lower elevation than the start point. We then determined all LCPs between the start point and the destination points and counted the number of paths that passed through each cell of the resistance surface for each of the 10,000 start points. The results from each iteration were then summed to create a single raster with 30- by 30-m resolution where the cell values represent the total number of LCPs that passed through that cell. In general, the resulting raster had very narrow corridors because LCPs between different pairs of points may coalesce to a single, low-cost path. However, true corridors will be wider than the narrow paths identified by the algorithm. To account for this issue, we smoothed the raster using a circular moving average with a radius of 439 m, which represented the mean daily distance traveled by sage-grouse (Coates and others, 2016b).

We then converted this raster to binary values to incorporate the results into the previously described map of example habitat management categories. Because we were interested in finding corridors that were not already included in the priority habitat classification, we first masked out all areas of the raster that fell into the high-use ASUI class. We then converted all the cells of this raster that were in the 95th percentile or above to 1, whereas all other cells were converted to 0 (after visual inspection, we determined that using the 95th percentile as the cutoff point produced reasonable results). Although this process captured the areas with the most LCPs, in some cases this resulted in patches of cells that did not connect to a high ASUI area or source area. Therefore, we removed contiguous areas that did not connect to a high-use ASUI area and source area. However, this does not remove offshoots that may not reach either a high ASUI area or a source area. To identify and remove these areas, we added an additional processing step. We first masked the resistance surface and brood source areas with the binary raster such that all areas where the binary raster had “0” were removed from the resistance and source rasters. These two rasters were then used to rerun the LCP analysis, using the original 10,000 starting points. This restricted all paths to fall within the already-identified corridors and meant that no paths would travel through offshoots that did not lead to either a high ASUI area or a source area. Using these results, we removed all areas from the original binary raster that did not have a path within 439 m to match the moving window distance used previously. The resulting raster represented the final corridors, which were incorporated as general habitat in the map of example habitat management categories.

Finally, we masked out perennial waterbodies and towns, which do not represent habitat for sage-grouse. To capture recreational activity around urban centers, we buffered town boundaries based on population size using an asymptotic regression curve with a maximum of 1,000 m, which was chosen because all leks were a minimum of 1,100 m from towns, and we did not want to remove any lek locations from designated habitat. We also masked all buildings regardless of location using the Esri building footprints layer (https://hub.arcgis.com/maps/esri::microsoft-building-footprints-tiles/about), which we buffered by 100 m to capture the land parcel and infrastructure associated with each building. Urbanization can also spread beyond town boundaries, however, so we hand digitized all areas within 10 km of town buffers with permanent structures based on both the Esri building footprints layer and the National Land Cover Database’s imperviousness layer (Dewitz and U.S. Geological Survey, 2021). Digitized areas were meant to capture regions around neighborhoods or concentrated structures, not areas around isolated houses or ranches. We then reduced these digitized areas by one level (that is, from priority+ to priority, priority to general, general to other, or other to non-habitat). To ensure the adoptability and effectiveness of the final output, we employed a months-long and iterative stewardship-based process of coproduction to create the final map of example habitat management categories, the process of which is described in appendix 3. We provide all raster surfaces, including selection maps, survival maps, example habitat management categories, source habitat, sink habitat, pre-fire priority habitat, and corridors to support land management and conservation efforts (Coates and others, 2024).

Nest Selection and Survival

We located and monitored 1,220 nests from 885 individual hens across 13 years and 15 study sites. The 95th percentile of distances between lek locations and nearest nest varied by study site, with a median of 6,332 m (range: 4,260–13,873 m). From variable selection analyses (for full results, see appendix 4), the most effective landscape predictor variables from correlated groups of variables for the nest selection model were bare ground (r=75 m), litter cover (r=439 m), perennial grass cover (r=1,451 m), total shrub cover (RAP; r=75 m), sagebrush cover (r=1,451 m), tree cover (r=167 m), elevation (r=167 m), HLI (r=439 m), TRASP (r=1,451 m), topographic roughness (r=439 m), TPI (r=167 m), total stream density (r=1,451 m), perennial stream density (r=1,451 m), alfalfa (r=1,451 m), and proximity to wetland (α=0.21). The most effective predictors for the nest survival model were annual grass cover (r=167 m), bare ground (r=1,451 m), total shrub cover (RCMAP; r=75 m), tree cover (r=1,451 m), HLI (r=439 m), and topographic roughness (r=370 m).

Nesting sage-grouse demonstrated apparent selection of greater values for litter cover, total shrub cover, sagebrush cover, elevation, TRASP (more south-facing slopes), and total stream density, while demonstrating apparent avoidance of greater values of bare ground, perennial grass cover, tree cover, HLI, topographic roughness, TPI, perennial stream density, alfalfa, and proximity to wetland (table 1; fig. 2). All coefficients were interpreted as strongly informative, with all effect probabilities exceeding 0.95 (table 1). Our nest selection model had reasonable out-of-sample predictive capabilities (Spearman’s rank ρ=0.96; R²=0.85; β_predict=0.86)

Habitat effects on nest survival were positive for total shrub cover, tree cover, HLI, and topographic roughness and were negative for annual grass cover and bare ground (table 2; fig. 2). All six habitat effects were interpreted as at least moderately informative (pd>0.9; table 2). In addition, there was little evidence for an effect of hen age on nest survival. Throughout the study area, cumulative 38-day nest survival was estimated to be 0.25 (95-percent credible intervals [CRI]=0.20–0.30). Predicted survival estimated from our final model of nest survival was consistent when compared with testing data (Bayesian P-value=0.59) and training data (Bayesian P-value=0.60) encounter histories, evidence that our model had good model fit and strong out-of-sample predictive capabilities.

Brood Selection and Survival

We monitored 480 broods from 437 individual hens across all years and study sites collecting a total of 1,155 locations during the early brood-rearing period and 1,469 locations during the late brood-rearing period. The 95th percentile of distances between an individual’s brood location and its nest varied by study site with a median of 4,945 m during the early brood-rearing period (range: 1,123–11,758 m) and 11,631 m during the late brood-rearing period (range: 2,394–54,036 m). From variable selection analysis (for full results, see appendix 4), the most effective landscape predictor variables from correlated groups for the early brood habitat selection model component were bare ground (r=75 m), total shrub cover (RAP; r=75 m), tree cover (r=167 m), elevation (r=167 m), HLI (r=75 m), topographic roughness (r=75 m), TPI (r=1,451 m), total stream density (r=370 m), proximity to alfalfa (α=1,451), and proximity to wetland (α=0.24). For the late brood habitat selection component, the most effective predictor variables were annual grass cover (r=439 m), bare ground (r=75 m), litter cover (r=1451 m), perennial grass cover (r=75 m), total shrub cover (RCMAP; r=75 m), tree cover (r=439 m), elevation (r=260 m), HLI (r=370 m), topographic roughness (r=260 m), TPI (r=1,451 m), proximity to stream (α=1,451), proximity to alfalfa (α=1,451), whether or not an area was burned, and an interaction between annual grass and whether or not an area was burned. For brood survival models, during the early brood component the most effective predictors from variable selection were bare ground (r=439 m), litter cover (r=1,451 m), non-sagebrush shrub cover (r=1,451 m), tree cover (r=370 m), TRASP (r=1,451 m), topographic roughness (r=167 m), total stream density (r=167 m), proximity to alfalfa (α=1,451), and whether or not an area was burned. During the late brood component, the most effective survival predictors were annual grass cover (r=75 m), bare ground (r=1,451 m), TPI (r=75 m), total stream density (r=1,451 m), and proximity to perennial streams (α=1,451 m).

Broods demonstrated apparent selection of greater values for total shrub cover, elevation, TPI, and areas closer to alfalfa, while demonstrating apparent avoidance of greater values of bare ground, tree cover, HLI, and topographic roughness, during the early and late periods (table 1; fig. 3). Early broods demonstrated apparent avoidance of areas with higher total stream densities and closer proximities to wetlands, whereas late broods demonstrated apparent selection of areas closer to streams (table 1; fig. 3). Late broods also demonstrated apparent selection of annual grass, litter cover, perennial grass cover, and burned areas but demonstrated decreased selection for annual grass within burned areas. All habitat effects for the early and late brood-rearing periods were interpreted as moderately or strongly informative (table 1). Selection models had adequate out-of-sample predictive capabilities for the early (Spearman’s rank ρ=0.90; R²=0.73; β_predict=0.90) and late (Spearman’s rank ρ=0.95; R²=0.82; β_predict=0.68) brood-rearing periods.

Habitat effects on brood survival were negative for bare ground for early and late periods. Early broods also demonstrated positive responses to litter cover, non-sagebrush shrub cover, tree cover, TRASP, and burned areas while demonstrating negative responses to topographic roughness, proximity to alfalfa, and total stream density (table 2; fig. 3). Late brood survival was lower with greater values of TPI and higher with more annual grass cover, greater total stream densities, and closer to perennial streams. Of the nine habitat effects estimated for early brood survival, seven were interpreted as moderately or strongly informative, with burned area indicating weaker evidence (table 2). For late brood survival, effects of annual grass cover, bare ground, and TPI were interpreted as moderately informative, with the remaining two habitat effects showing weaker evidence of effects (table 2). In addition, there was not moderate or strong evidence of an effect of hen age on brood survival, but there was evidence for a positive effect of brood age on survival (table 2). The estimated daily survival rate for broods was 0.985 (95-percent CRI=0.980–0.990). Predicted survival estimated from our final model of brood survival was consistent when compared with testing data (Bayesian P-value=0.64) and training data (Bayesian P-value=0.58) encounter histories, evidence that our model had good model fit and strong out-of-sample predictive capabilities.

Spring Selection

We monitored 2,171 individuals during the spring across years and study sites, collecting a total of 22,735 locations. The median 99-percent kernel UD for all birds at a study site was 1,431 square kilometers (km²; range: 720–3,743 km²). From variable selection analysis (for full results, see appendix 4), the most effective landscape predictor variables from correlated groups for the spring habitat selection model component were annual grass cover (r=75 m), bare ground (r=75 m), litter cover (r=439 m), perennial grass cover (r=1,451 m), total shrub cover (RAP; r=439 m), sagebrush cover (r=1,451 m), tree cover (r=439 m), elevation (r=167 m), HLI (r=439 m), TRASP (r=75 m), topographic roughness (r=167 m), TPI (r=439 m), total stream density (r=1,451 m), perennial stream density (r=1,451 m), log-transformed distance to alfalfa, and proximity to wetland (α=0.22). In the spring, sage-grouse demonstrated apparent selection for greater values of litter cover, sagebrush cover, elevation, TPI, total stream densities, and closer to alfalfa, while demonstrating apparent avoidance of areas with greater values for annual grass cover, bare ground, perennial grass cover, total shrub cover, tree cover, HLI, TRASP, topographic roughness, perennial stream densities, and closer to wetland (table 1; fig. 4). All 16 habitat effects were interpreted as strongly informative (table 1). The spring selection model had strong out-of-sample predictive capabilities (Spearman’s rank ρ=0.99; R²=0.94; β_predict=0.91).

Summer Selection

We monitored 1,778 individuals during the summer across years and study sites, collecting a total of 5,919 locations. The median 99-percent-kernel UD for all birds at a study site was 1,362 km² (range: 493–3,141 km²). From variable selection analysis (for full results, see appendix 4), the most effective landscape predictors from correlated groups for the summer habitat selection model component were bare ground (r=1,451 m), litter cover (r=75 m), perennial grass cover (r=1,451 m), total shrub cover (RAP; r=1,451 m), total shrub cover (RCMAP; 75 m), tree cover (r=1,451 m), elevation (r=167 m), HLI (r=1,451 m), TRASP (r=1,451 m), topographic roughness (r=167 m), TPI (r=439 m), proximity to stream (α=540), proximity to perennial stream (α=1,451), log-transformed distance to alfalfa, proximity to wetland (α=0.18), and whether or not an area was burned. In the summer, sage-grouse demonstrated apparent selection of areas with greater values of total shrub cover within 75 m, elevation, TRASP, TPI, closer to total and perennial streams, and closer to alfalfa (table 1; fig. 4). Sage-grouse demonstrated apparent avoidance of areas with greater values of bare ground, litter cover, perennial grass cover, total shrub cover within 1,451 m, tree cover, HLI, topographic roughness, closer to wetlands, and burned areas (table 1; fig. 4). All 16 habitat effects were interpreted as strongly informative (table 1). The summer selection model had strong out-of-sample predictive capabilities (Spearman’s rank ρ=1.00; R²=0.92; β_predict=0.86).

Winter Selection

We monitored 264 individuals during winter across years and study sites, collecting a total of 5,508 locations. The median 99-percent kernel UD for all birds at a study site was 2,058 km² (range: 946–4,110 km²). From variable selection analysis (for full results, see appendix 4), the most effective landscape predictor variables from correlated groups for the winter habitat selection model component were annual grass cover (r=370 m), bare ground (r=1,451 m), litter cover (r=439 m), total shrub cover (RCMAP; r=1,451 m), tree cover (r=1,451 m), elevation (r=75 m), HLI (r=1,451 m), TRASP (r=1,451 m), topographic roughness (r=439 m), TPI (r=439 m), proximity to stream (α=660), perennial stream density (r=439 m), proximity to alfalfa (α=1,451), and proximity to wetland (α=0.21). Wintering sage-grouse demonstrated apparent selection of areas with greater values for annual grass cover, litter cover, total shrub cover, elevation, TRASP, TPI, and closer to streams and alfalfa, while demonstrating apparent avoidance of areas with greater values for bare ground, tree cover, HLI, topographic roughness, perennial stream densities, and closer to wetland (table 1; fig. 4). Of the 14 habitat effects included in the final model, all were interpreted as strongly informative (table 1). The winter selection model had strong out-of-sample predictive capabilities (Spearman’s rank ρ=0.99; R²=0.97; β_predict=0.92).

Movement Model

We monitored 18 GPS-marked individuals to inform the resistance surface that was used to quantify connectivity between nesting and brood-rearing habitat. Based on an integrated step-selection model, females moving from nesting to brood-rearing habitat demonstrated apparent selection for greater total shrub cover, litter cover, higher densities of streams, areas closer to wetlands, higher elevations, and greater values of TPI (table 3). In contrast, females demonstrated apparent avoidance of higher tree cover, more bare ground, greater terrain roughness, and greater values of HLI (table 3).

Results: Habitat Mapping

Maps demonstrating seasonal sage-grouse selection patterns indicated that habitat potential (estimated from habitat selection models) was similar during the spring (fig. 5) and summer (fig. 6) but more widespread in the winter (fig. 7) when considering the overall area falling into the moderate or high categories. For reproductive life stages, habitat potential based on habitat selection models was more widespread across the study area during the nesting (fig. 8) and early brood-rearing periods (fig. 9) compared to the late brood-rearing period (fig. 10). Habitat suitability from survival models seemed possibly most limited during the nesting period (fig. 11) compared to the early brood-rearing (fig. 12) and late brood-rearing periods (fig. 13).

Areas of high-use based on the ASUI (fig. 14) generally aligned with areas of high use based on the habitat selection models (figs. 5–10) and the combined metric of habitat suitability and occupancy (appendix 5). Priority habitat seemed to be concentrated in the northern parts of California and Nevada and more fragmented in the southern parts of the study area (fig. 15). A total of 1.5 percent of the study area was classified as priority+ (5,112 km²), 15.5 percent as priority (51,683 km²), 6.4 percent as general (21,359 km²), and 9.1 percent as other (30,410 km²) with the remainder classified as non-habitat. Including corridors that connect occupied nesting habitat to brood source habitat resulted in small increases in general habitat (0.18 percent), particularly in the southern part of the study area (fig. 16).