Wildfire and the ecological niche: Diminishing habitat suitability for an indicator species within semi-arid ecosystems.

Globally accelerating frequency and extent of wildfire threatens the persistence of specialist wildlife species through direct loss of habitat and indirect facilitation of exotic invasive species. Habitat specialists may be especially prone to rapidly changing environmental conditions because their ability to adapt lags behind the rate of habitat alteration. As a result, these populations may become increasingly susceptible to ecological traps by returning to suboptimal breeding habitats that were dramatically altered by disturbance. We demonstrate a multistage modeling approach that integrates habitat selection and survival during the key nesting life-stage of a bird species of high conservation concern, the greater sage-grouse (Centrocercus urophasianus; hereafter, sage-grouse). We applied these spatially explicit models to a spatiotemporally robust dataset of sage-grouse nest locations and fates across wildfire-altered sagebrush ecosystems of the Great Basin ecoregion, western United States. Female sage-grouse exhibited intricate habitat selection patterns that varied across regional gradients of ecological productivity among sagebrush communities, but often selected nest sites that disproportionately resulted in nest failure. For example, 23% of nests occurred in wildfire-affected habitats characterized by reduced sagebrush cover and greater composition of invasive annual grasses. We found survival of nests was negatively associated with wildfire-affected areas, but positively associated with higher elevations with increased ruggedness and overall shrub cover. Strong site fidelity likely drove sage-grouse to continue nesting in habitats degraded by wildfire. Hence, increasing frequency and extent of wildfire may contribute disproportionately to reduced reproductive success by creating ecological traps that act as population sinks. Identifying such habitat mismatches between selection and survival facilitates deeper understanding of the mechanisms driving reduced geographic niche space and population decline at broad spatiotemporal scales, while guiding management actions to areas that would be most beneficial to the species.

INTRODUCTION

Global conservation policy is challenged by population decline, extirpation, and species extinction resulting from altered disturbance regimes that are hallmarks of Anthropocene environments. One dominant form of disturbance currently altering resource distribution and ecosystems at global scales is wildfire (Flannigan, Krawchuk, de Groot, Wotton, & Gowman, 2009). Sudden shifts in wildfire regimes have been attributed to a warming and drying climate (Pechony & Shindell, 2010), yet direct human interactions (Syphard, Keeley, Pfaff, & Ferschweiler, 2017) and expanding invasive plant species (D’Antonio & Vitousek, 1992; Pausas & Keeley, 2014) also play a significant role.

Environmental stressors often act in tandem to promote changes in wildfire regimes. During recent decades, a positive feedback cycle between non‐native annual grasses and wildfire has been identified, wherein fire frequency and size are amplified where encroachment of non‐native invasive grasses has occurred (Balch, Bradley, D’Antonio, & Gómez‐Dans, 2013). Exotic grasses facilitate increases in highly flammable fuel biomass and can outcompete native plant communities following disturbance, particularly in increasingly arid environments (Brooks & Chambers, 2011; Brooks et al., 2004; Shriver et al., 2018). In some circumstances, the grass‐fire cycle can be exacerbated by changing climate patterns, especially when seasonally high precipitation is followed by extreme drought (Balch et al., 2013; Chambers et al., 2014). As fire frequency increases relative to historical regimes (Whisenant, 1990), native vegetation communities become more susceptible to permanent state transitions, and wildlife communities are at greater risk of local extirpation due to loss of habitat (Coates, Ricca, et al., 2016; Rockweit, Franklin, & Carlson, 2017).

Habitat specialists are disproportionately sensitive to altered disturbance regimes characteristic of global change (Clavel, Julliard, & Devictor, 2011). According to ecological niche theory (Hutchinson, 1957), specialist species perform exceptionally well under a relatively narrow range of environmental conditions (i.e., niche width or niche breadth), while generalist species inhabit a wider range of conditions with more moderated performance throughout (Futuyma & Moreno, 1988). Because they are less adaptive to rapidly changing conditions and/or increased competition from generalist species, specialists tend to suffer population declines when rapid change occurs (Munday, 2004). In many systems, wildfire disturbance causes abrupt and immediate loss of habitat for specialists (Brooks, Matchett, Shinneman, & Coates, 2015). While prefire communities may eventually recover, short‐term losses in habitat are likely to have population impacts, as rapid environmental change can lead to mismatches between habitat occupancy or selection and the fitness potential of those habitats.

Because individuals within a population vary in their ability to recognize habitat quality, some make “good” decisions, whereas others misinterpret environmental cues and occupy sites that fail to maximize their chances of survival and reproduction (Raynor, Beyer, Briggs, & Joern, 2017). Mismatches between individual habitat selection and individual success can result in the emergence of ecological traps (Battin, 2004; Robertson & Hutto, 2006), that can have adverse population‐level effects when many individuals are affected (Remeš, 2000; Schlaepfer, Runge, & Sherman, 2002). One notable circumstance that can lead to an ecological trap is an abrupt decline in habitat suitability, while selection cues used by individuals are unchanged (Robertson & Hutto, 2006). While adaptive mechanisms such as breeding site fidelity may foster population growth under relatively stable conditions for specialist species (Chalfoun & Schmidt, 2012; Schroeder & Robb, 2003), these behaviors may result in maladaptive habitat selection following disturbance (Piper, Palmer, Banfield, & Meyer, 2013). Rapid environmental change can thereby lead to mismatches between occupancy and performance, with associated population impacts depending on the magnitude of change. Such occurrences may be increasingly common under altered disturbance regimes, possibly resulting in apparent sinks or pseudosinks (Heinrichs, Lawler, & Schumaker, 2016; Remeš, 2000; Watkinson & Sutherland, 1995), where a habitat's suitability is reduced such that it cannot sustain long‐term population persistence without immigration (e.g., true sink; Pulliam, 1988), or its carrying capacity has substantially declined (e.g., pseudosink; Watkinson & Sutherland, 1995). Because spatial variation in reproductive success drives population‐level source–sink dynamics (Heinrichs et al., 2016; Pulliam, 1988), a more comprehensive understanding of variation in habitat contributions to species performance within landscapes fragmented by wildfire is needed to identify the existence of maladaptive habitat selection, ecological traps, and novel sink habitats (Remeš, 2000; Robertson, Rehage, & Sih, 2013).

Complex drivers of increasing wildfire, subsequent habitat losses, and decline of specialist species necessitate tractable tools to help identify and ameliorate impacts on the available geographic niche space of sensitive wildlife species using reliable measures of habitat suitability (Fordham, Akçakaya, Araújo, Keith, & Brook, 2013). Science‐based conservation action thereby requires identifying geographic areas in altered environments with characteristics that meet specific life‐history requirements of vulnerable populations. In practice, habitat suitability models characterize the ecological niche, or the “n‐dimensional hypervolume” of environmental conditions that results in long‐term population persistence (Hirzel & Le Lay, 2008; Holt, 2009; Hutchinson, 1957; Pulliam, 2000). Projecting geographic realizations of the ecological niche across broad spatial expanses has become integral to allocating conservation resources across species’ actual or potential range (Guisan & Thuiller, 2005; Ocampo‐Peñuela, Jenkins, Vijay, Li, & Pimm, 2016; Panzacchi, Van Moorter, Strand, Loe, & Reimers, 2015; Santika, McAlpine, Lunney, Wilson, & Rhodes, 2015), especially under rapidly changing environmental conditions driven by disturbance.

However, habitat suitability modeling often focuses on species distribution (Gaillard et al., 2010; Guisan & Zimmermann, 2000; Hirzel & Le Lay, 2008), without providing information about species performance. In rapidly changing environments, species performance may be more indicative of the fundamental ecological niche (that is, where population growth rate is positive on average; Panzacchi et al., 2015; Pulliam, 2000), providing early warning signs of potential extirpation in poorly performing habitats, whereas species distribution may not reflect such losses in habitat suitability until an area is no longer occupied. To better meet conservation challenges, habitat suitability models would benefit from moving beyond measures of apparent habitat suitability that are based primarily on species’ presence and may not reliably predict demographic responses to underlying environmental characteristics (Heinrichs, Bender, Gummer, & Schumaker, 2010; Matthiopoulos et al., 2015; Thuiller et al., 2014).

Our objectives were twofold. First, we demonstrate a spatially explicit modeling approach to address multiscale habitat selection and variability in demographic responses within a wildfire‐impacted ecoregion. Our approach was designed in part to identify areas where habitat selection patterns deviated from reproductive success, thereby implying maladaptive habitat selection. Second, using the results of our models, we show how a shift in habitat associated with increasing wildfire frequency and severity contributes to shrinking geographic niche space for a specialist species. We applied our models to a robust spatiotemporal dataset of nest locations and fates of greater sage‐grouse (Centrocercus urophasianus; hereafter, sage‐grouse), a ground‐nesting sagebrush‐obligate species that acts as an indicator for the ecological integrity of the sagebrush biome (Hanser & Knick, 2011) and occupies semi‐arid ecoregions characteristic of accelerated wildfire regimes (Scholze, Knorr, Arnell, & Prentice, 2006). While sagebrush cover and the predominance of drought drives range‐wide sage‐grouse distributions (Aldridge et al., 2008; Schroeder et al., 2004), reproduction within their population range depends on predator composition, structure of habitat patches, and localized disturbances (Aldridge & Boyce, 2007; Conover & Roberts, 2017; Heinrichs et al., 2018; Walker, Naugle, & Doherty, 2007). Under these conditions, the increasing frequency and extent of localized disturbance may act to decouple individual habitat selection from reproductive success, likely cascading to more widespread effects on species’ population viability.

A key example of such local‐to‐regional scaling patterns can be found in the Great Basin ecoregion of the western United States, a region larger than 75% of countries worldwide. Here declining sage‐grouse populations over 30 years are linked to cumulative impacts of an accelerated cycle of wildfire and Bromus tectorum (i.e., cheatgrass) invasion that destroys thousands of hectares of sagebrush annually and mimics effects of drought (Balch et al., 2013; Coates, Ricca, et al., 2016). While a few local‐level studies have identified specific reproductive life‐stages impacted by wildfire (Foster, Dugger, Hagen, & Budeau, 2019), the demographic mechanism driving patterns of population decline at large ecoregional scales is uncertain. Hence, large gaps remain in understanding how sage‐grouse habitat suitability has been altered by wildfire. For example, the nesting life‐stage may be especially vulnerable to maladaptive habitat selection (Cutting et al., 2019) resulting in high incidence of reproductive failure because sage‐grouse exhibit strong interannual fidelity to nesting sites (Fischer, Apa, Wakkinen, Reese, & Connelly, 1993; Holloran & Anderson, 2005). As sage‐grouse population growth is sensitive to variation in nesting survival (Taylor, Walker, Naugle, & Mills, 2012), identifying such patterns of spatial heterogeneity, particularly in relation to widespread wildfire and annual grass invasion, could indicate possible ecological traps and help guide the prioritization of habitats for strategic conservation actions.

MATERIALS AND METHODS

Overview

Using sage‐grouse as a focal species, we quantified scale‐dependent factors driving nest site selection and nest survival across broad spatial scales in order to identify wildfire impacts and other environmental influences on variation in nesting productivity across a broad ecoregion spanning mesic and xeric shrub communities. To investigate the consequences of habitat selection and explore the potential for a source–sink reproductive landscape, we sought to classify nesting habitat on a scale ranging from adaptive (high selection, high survival) to maladaptive (high selection, low survival). By incorporating regional heterogeneity in vegetation communities through differences in soil moisture and temperature that index ecological productivity (Chambers et al., 2014), the resulting model would be transferable across regions, providing a scientific mapping tool to identify variation in potential reproduction based on likelihood of nest site selection and subsequent likelihood of nest survival across a wildfire‐impacted ecoregion. By spatially integrating habitat selection with reproductive success, this modeling approach is a broadly applicable framework for any species or community that can be monitored across expansive landscapes. We synthesized the resulting predictive mapping products to demonstrate the scale and prevalence of maladaptive breeding habitat selection for sage‐grouse within an ecosystem currently being impacted by rapid environmental change, primarily from wildfire.

Data collection and materials

We collected data at 12 field sites within the Great Basin ecoregion during 2009–2018, where data from 2009 to 2017 were used as training data and 2018 was withheld for validation (Appendix S1‐Figure S1). Vegetation was dominated by Wyoming big sagebrush (Artemisia tridentata wyomingensis) and black (Artemisia nova) and low (Artemisia arbuscula) sagebrush, depending on elevation (Appendix S1). At high elevations, mountain big sagebrush (Artemisia tridentata spp. vaseyana) was more frequent. Non‐sagebrush shrubs included rabbitbrush (Chrysothamnus ssp.), Mormon tea (Ephedra viridis), snowberry (Symphoricarpos ssp.), western serviceberry (Amelanchier alnifolia), and antelope bitterbrush (Purshia tridentata). Conifer forests most frequently comprised single‐leaf pinyon pine (Pinus monophylla) and Utah juniper (Juniperus osteosperma; hereinafter, “pinyon‐juniper”). Non‐native annual grasses included cheatgrass (B. tectorum) and medusahead rye (Taeniatherum caput‐medusae). Native perennial grasses included needle and thread (Hesperostipa comata), Indian ricegrass (Achnatherum hymenoides), and squirreltail (Elymus elymoides).

We captured sage‐grouse using nighttime spotlighting (Giesen, Schoenberg, & Braun, 1982; Wakkinen, Reese, Connelly, & Fischer, 1992) and deployed very high frequency (VHF), necklace‐style radiotransmitters (Kolada, Sedinger, & Casazza, 2009). A subset of individuals was fitted with combined GPS Platform Transmitter Terminals (North Star Science and Technology, LLC, King George, Virginia) and VHF radiotransmitter systems. We verified locations of nests visually when telemetry data indicated two consecutive observations at the same location for a female sage‐grouse. We took extra care to avoid flushing the female. We then attempted to relocate the female twice weekly during nesting (March–June) to monitor status and nest fate. Nests were considered successful if ≥1 egg hatched.

Cumulative burned area

We implemented a temporally dynamic spatial data model to represent annual “cumulative burned area” surfaces (hereafter, CBA) across sagebrush communities within our study area, adapted from Coates, Ricca, et al. (2016). Importantly, this model accounted for time lag effects in wildfire by allowing recovery rates to vary across different sagebrush communities within burned areas, dependent on underlying conditions of elevation, soil moisture, and soil temperature gradients. First, we identified wildfire perimeters that had burned intensely enough to remove underlying vegetation from national burn severity mosaics spanning 1984–2016. Using an iterative approach, we then updated burned area on an annual basis by adding new fires that occurred while allowing sagebrush to initiate recovery in burned areas following a fire. The spatial data model was initially defined at a 30 m pixel resolution, where each cell took on the value of 0 (unburned or recovered), or 1 (indicating burned and/or not recovered), with initial settings based on first year of data availability (1984). Sagebrush recovery rates were based on sagebrush community type (mountain big sagebrush, Wyoming big sagebrush) and classification indices of resilience to disturbance and resistance to invasion (RR), which was characterized by gradients of elevation, soil temperature, and soil moisture (Brooks et al., 2015; Chambers et al., 2014; Maestas, Campbell, Chambers, Pellant, & Miller, 2016). Sagebrush recovery was reflected in the CBA model once 20% of prefire cover was reached (Connelly, Schroeder, Sands, & Braun, 2000), with recovery projected to occur after 9 years in mountain sagebrush vegetation communities with cool, moist soils (high RR), and after 16 years in moderate RR classes. In contrast, permanent state transitions to annual grass were predicted in Wyoming big sagebrush communities with warmest and driest soil class (low RR; Table S2; Coates, Ricca, et al., 2016; Miller et al., 2011). Recovery values were reset to 0 when new burns occurred. Once CBA models were generated across all years, we updated sage‐grouse nest locations with underlying CBA values associated with the appropriate year, to be evaluated alongside other landscape variables described in the following section. Further detail of CBA methodology is reported in Appendix S2.

Landscape spatial data

Using Geographic Information Systems (GIS), we quantified a suite of candidate predictors to explain variation in nest site selection and survival. Predictors represented shrubland and grassland vegetation, anthropogenic infrastructure, water sources, topography, and soil‐based measures of ecological productivity, along with CBA. Vegetation datasets characterizing percentages of shrubland and grassland cover types were provided by US Geological Survey via the National Land Cover Database (Xian, Homer, Rigge, Shi, & Meyer, 2015). These data included layers depicting sagebrush and shrub cover and height components (e.g., mountain big sagebrush, Wyoming big sagebrush, low sagebrush, rabbitbrush, bitterbrush), bare ground, residual litter, annual grasses, and herbaceous vegetation. We considered agricultural cropland, woodlands, and wet meadows using Landsat‐based mapping products, and included a high‐resolution mapping product to represent pinyon‐juniper conifer cover (Gustafson et al., 2018). We considered the Normalized Difference Vegetation Index (NDVI) to further represent vegetation productivity. To represent topography, we mapped elevation, slope, transformed aspect, topographic roughness, and several additional indices to characterize heat load, curvature, and potential wetness (Appendix S3). Hydrographic features, including streams, springs, and other water bodies were gathered from the National Hydrography Dataset (US Geological Survey, 2017). Anthropogenic infrastructure was represented in the form of roads and transmission lines, while CBA was estimated for each year of the study (see Section 2.2; Appendix S2).

Because species relationships with environmental attributes are hierarchical and scale‐dependent (McGarigal, Wan, Zeller, Timm, & Cushman, 2016), we evaluated predictors at four spatial scales based on estimated movements of telemetered sage‐grouse (Coates, Casazza, et al., 2016). Using focal neighborhood analyses (Hagen‐Zanker, 2016) with ArcGIS® Spatial Analyst (ArcGIS 10.4: Environmental Systems Research Institute, Inc.), we summarized data at four spatial scales by evaluating each predictor within circular moving windows with radii of 75 m (1.8 ha), 167.9 m (8.7 ha), 439.5 m (61.5 ha), or 1,451.7 m (661.4 ha), where each 30 m pixel represented the contribution of all surrounding pixels occurring within the moving window (see Appendix S3 for details). We also estimated distances to features such as streams, wet meadows, roads or transmission lines. All GIS‐based spatial variables and scales are described with complete details in Appendix S3‐Table S3.

We sought to quantify habitat using covariates with predictive properties that were transferable to other locations and time periods (Paton & Matthiopoulos, 2016) for the purposes of large‐scale mapping of nesting habitat. One approach for addressing regional variability in species‐habitat relationships involves modeling “functional responses in habitat selection,” where local habitat preferences depend on habitat availability at some larger spatial scale (Aarts, Fieberg, Brasseur, & Matthiopoulos, 2013; Mysterud & Ims, 1998). We therefore hypothesized that sage‐grouse respond to regional variation in ecological productivity and shrubland community composition along gradients of soil temperature and moisture that correspond with ecological resistance and resilience (RR; see Section 2.2). Because the index of RR characterizes productivity of sagebrush communities and has been mapped across western ecosystems (Chambers et al., 2014; Maestas et al., 2016), we applied it as a proxy for regional habitat availability. That is, we considered habitat functional response interactions (Aarts et al., 2013; Matthiopoulos, Hebblewhite, Aarts, & Fieberg, 2011; Mysterud & Ims, 1998) of all vegetation covariates (e.g., shrub, grass and conifer cover) with a regional index of RR (adapted from Maestas et al., 2016) that distinguished mesic areas with cooler temperatures and greater soil moisture from drier and more xeric areas. In contrast to other landscape predictors, the RR index was quantified and summarized at the scale of local sage‐grouse populations (moving window radius = 21.58 km; Appendix S3), to allow other habitat effects to vary depending on generalized regional conditions affecting the population (Aarts et al., 2013; Paton & Matthiopoulos, 2016; Appendix S3). Functional response interactions were specified as first‐order expectations (Matthiopoulos et al., 2011; Paton & Matthiopoulos, 2016), where RR was assumed to represent regional availabilities with respect to habitat predictors in the linear predictor l(x), which took the form:where each x was a habitat covariate representing local conditions and R was the regional RR index.

We also considered non‐linear functional relationships for each non‐distance variable by including second order polynomial terms (quadratic effect; Appendix S4). We modeled distance variables using non‐linear exponential decay functions (Appendix S3), to account for declining influence as distance from the feature increased (Coates, Brussee, et al., 2016; Nielsen, Cranston, & Stenhouse, 2009).

Variable reduction analysis

We implemented an iterative covariate reduction process to identify an appropriate spatial scale and functional form for each candidate predictor of nest site selection and survival. This was done because applying conventional variable selection approaches on the complete dataset would have been inappropriate and computationally infeasible due to many correlated predictors (|r| > .6; Dormann et al., 2013), and because Bayesian variable selection procedures would require unreasonably long simulations to adequately cover all possible models (Converse, Royle, Adler, Urbanek, & Barzen, 2013; Tenan, O’Hara, Hendriks, & Tavecchia, 2014).

The purpose of variable reduction analysis was to reduce the initial number of model covariates by first ranking and selecting the most influential functional form of each candidate predictor (linear, non‐linear, or interaction), and then removing candidate predictors that were highly correlated with other, higher ranked predictors from consideration (Appendix S4). Importantly, our approach accounted for partial correlations by repeatedly fitting multivariate models with subsets of independent predictors. We did this within a hierarchical simulation framework by specifying generalized linear mixed models (GLMM) using a binomial distribution (logit link for selection model, and log‐exposure link for nest survival; see Appendix S4) under the condition that between 5 and 15 fixed effects predictors were present in each model run, and that pairs of predictors could not have correlation greater than |0.6|. To rank importance for each predictor, we compared the AIC value of each model with and without that predictor, as well as models with the predictor represented as a quadratic term or an interaction with the RR metric. We stored the improvement in AIC (dAIC) for each candidate model, for each predictor. We repeated this process 10,000 times, and ranked predictors based on the average dAIC for their top functional form (linear, quadratic, or interaction). For correlated predictors, only the predictor with the greater dAIC was considered in subsequent models. The simulation thus facilitated selection of spatial scale as well as initial functional form for each predictor. We provide complete details in Appendix S4.

Model for nest site selection

We used a GLMM framework within a Bayesian modeling environment to estimate resource selection functions (RSF; Johnson, Nielsen, Merrill, McDonald, & Boyce, 2006; Northrup, Hooten, Anderson, & Wittemyer, 2013) for nesting sage‐grouse. The RSF approach allowed us to spatially infer the relative probability of nest site selection across our study extent. Our base model took the form:where β0 represents the baseline intercept, Xβ is a vector of selection coefficients β multiplied by fixed covariates X, and γ and η are random intercepts for individual and year, respectively (Gillies et al., 2006). The observations Y followed a Bernoulli distribution, with y = 1 indicating a nest, and y = 0 indicating a random background location.

Background locations were used to characterize available nest site locations (Northrup et al., 2013). For each nest, we conditioned the available distribution of background locations on proximity to nearby leks because females often nest within ~10 km of their affiliated lek (Coates et al., 2013; Holloran et al., 2005). At each field site, we identified active leks within 20 km of any nest. For each nest, we then calculated distance to nearest lek. We used the maximum resulting distance (~17 km) to buffer around leks, creating a regional boundary for random sampling. We generated 10 random locations per nest, per region. We also included a distance to lek predictor in our models to account for potentially confounding clustering of nest locations relative to lek locations. At all locations, we used landscape habitat predictors identified from the variable reduction process to represent habitat, so that strongly correlated variables (r > .6) could not co‐occur in models.

We used Integrated Nested Laplace Approximations (www.r‐inla.org; Rue, Martino, & Chopin, 2009) in R v3.4.3 (R Core Team, 2019) to fit Bayesian hierarchical models. We specified uninformative priors for nest selection model coefficients, , following Rue et al. (2009), where is the precision parameter (). Random intercepts were normally distributed (iid) random variables, with precision hyperparameter ~ logGamma (a = 1, b = 0.01), where a and b were the shape and scale, respectively. After fitting a full model with parameters identified by our initial variable reduction procedure, we reduced the model by sequentially culling parameters with approximate Z < 1.44 (e.g., , where is the mean of the posterior distribution). The value of 1.44 corresponds to a parameter with 85% credible interval (CRI) not overlapping zero, which was useful to identify parameters that did not contribute meaningfully to the model (Arnold, 2010). We computed the fitted RSF by discarding the “nuisance” intercept terms (Johnson et al., 2006; McDonald, 2013) and including all selected covariate parameters in the exponential function:

Nest site selection model validation

After selecting a final model, we used independent testing data collected across seven field sites during the field season of 2018 (n = 144 nests) to validate the model's ability to predict nest locations, using cross‐validation methods for RSFs (Johnson et al., 2006) and used‐habitat calibration (UHC) plots (Fieberg et al., 2018). Using the final selected RSF model fit to training data (2009–2017 nest locations), we generated predictions of selection for the testing data. We saved calibration statistics for the set of testing/training data (Spearman's rank ρ, R 2, β predict) which summarized the association between expected habitat use across 10 habitat classes binned by RSF quantiles (from ) and the number of observed locations occurring within each class (Boyce, Vernier, Nielsen, & Schmiegelow, 2002; Fieberg et al., 2018; Johnson et al., 2006). If the model makes valid out‐of‐sample predictions, then Spearman's ρ and R 2 will be relatively high (>0.5). For additional details on generating calibration statistics, see Johnson et al. (2006) and Fieberg et al. (2018).

We generated UHC plots to compare observed distributions of each habitat variable (i.e., the distribution of values at all nest locations) in the model to its expected distribution (Fieberg et al., 2018). UHC plots are created by using the fitted RSF to assign weights to new locations, thus generating a predictive distribution of used habitat. This predictive distribution is then compared to the observed values (i.e., in the independent testing data) for each habitat variable. If the model predicts well, the observed used habitat in the independent testing dataset will fall primarily within the confidence bounds of the predicted used distribution (Fieberg et al., 2018). UHC plots effectively describe used distributions of habitat proportional to their availabilities, while also constituting a statistically rigorous model validation test. We summarized relative importance of each environmental predictor (main effect, plus additive or interaction effects), by computing the difference between the model‐predicted median used habitat distribution and the median of available habitat (Appendix S5). These differences were based on rescaled standardized distributions of each predictor. Because model coefficients reflect partial correlations conditional on other predictors in the model (i.e., not causal effects; Fieberg et al., 2018), we ranked and presented predictors based on distributions of disproportionate use to availability.

Model for nest survival

We used a Bayesian shared frailty model (Halstead, Wylie, Coates, Valcarcel, & Casazza, 2012) to model habitat covariate effects on sage‐grouse nest survival. This model accounts for exposure time by treating each day within each nest's exposure period as a Bernoulli trial (survived or not; Converse et al., 2013; Schmidt, Walker, Lindberg, Johnson, & Stephens, 2010). We created encounter histories for each nest that consisted of known alive and censored days throughout the nesting season, ending with a successful or failed nest. Dates of failure for failed nests were imputed between their last active date and the failure discovery date. Daily unit hazard (UH) and covariate effects on UH were the parameters of interest. This model took the formwhere the subscripts h, i, j, and k refer to individual nest, day, year and female, respectively, and Xβ is a vector of selection coefficients β multiplied by fixed covariates X. Covariates included all landscape predictors identified by the variable reduction procedure, age of female (yearling/adult), and estimated nest initiation day (relative to the mean for year j). We included random effects for year (), and female () to account for within‐individual autocorrelation and unexplained interannual variation. As was done with nest site selection, we culled parameters with Z < 1.44.

We specified uninformative prior distributions for all fixed effects, , where 0.001 was the precision (1/σ 2). Random effects received uninformative uniform priors (Halstead et al., 2012). We estimated nest survival expressed as: where CH represents the cumulative hazard across the complete incubation period (38 days) and s is the cumulative survival rate. Using MCMC, we ran three chains of 50,000 iterations, following a burn‐in of 10,000 iterations, and retained every 10th sample. We implemented frailty models using JAGS 4.2.0 and the package rjags (Plummer, 2003, 2016) for implementation in R (R Core Team, 2019). We visually examined chains and calculated Gelman–Rubin statistics to verify chain convergence ( < 1.05). We report median values of the posterior distribution and 95% CRI for parameter estimates, unless otherwise stated.

Nest survival model validation

We validated the frailty model by comparing fates and encounter histories from the 2018 testing dataset to “replicate” encounter histories generated from the model fit to the training data. Specifically, we calculated a Bayesian predictive p value (Gelman et al., 2013), a measure of fit that compares the empirical survival curve to the predicted survival curve by repeatedly generating data from the model's posterior predictive distribution (Schmidt et al., 2010). From each draw (1:1,000), we calculated the fit statistic which took the form:where is the cumulative survival probability for nest i through time t and are the observations (1 = active, 0 = failed). We calculated D for the observed data and for replicate data and reported the proportion of iterations where (e.g., Schmidt et al., 2010). Contrary to conventional p values in frequentist statistics, this p value represents how unlikely the observed data are compared to simulated data, where values near 0 or 1 indicate that the model rarely produces data resembling the real data (Gelman et al., 2013; Schmidt et al., 2010). To evaluate whether the model was consistently overestimating or underestimating survival of 2018 nests, we also recorded the error between the observed daily survival rate according to known encounter histories and the daily survival rate from predicted encounter histories. Ideally, these values would be centered around 0.

Selection and survival mapping

To map the predicted contribution of pixels on the landscape to sage‐grouse nesting productivity, we generated spatially explicit maps of predicted RSF and predicted cumulative (38 day) nest survival. We reclassified the nest selection map from a continuous surface representing relative nest site selection into categories based on habitat used by sage‐grouse: non‐habitat, low, moderate, high, and highest relative selection. These classes were determined using a percent isopleth approach adapted from Doherty, Evans, Coates, Juliusson, and Fedy (2016). First, nest RSF values were extracted at all training nest locations. Then, nests were sorted by RSF value and partitioned at 25% intervals with the first cutoff value indicating the top 25% RSF values (25% of nests had RSF > x 1). The second and third partitions occurred at the 50th and 75th percentiles, respectively. Values below the 95th percentile cutoff were considered non‐habitat.

To reclassify the survival map, we extracted model‐predicted survival values to the training nest locations. We used a kernel density estimator to approximate the probability distribution of predicted survival at failed and successful nests. We then used the mean of the joint distribution of all nests to distinguish between nests more likely to be successful (predicted survival above the mean) and those less likely to be successful (predicted survival below the mean). We created additional partitions at the 25th percentile of the distribution of failed nests (very low predicted survival), and the 75th percentile of the distribution of successful nests (very high predicted survival), resulting in four classes overall: very low (values <25th percentile of failed), low (values ≥25th percentile of failed but less than mean overall predicted), high (values greater than or equal to mean overall predicted but less than 75th percentile of successful), and very high (values ≥75th percentile of successful nests).

We combined the reclassified selection and survival maps into a ranked index that represented potential habitat contributions to nesting productivity at each pixel. This map had 16 classes ranking from high (1) to low (16). The highest habitat rank occurred where the highest habitat selection classification intersected the highest survival classification. The lowest rank occurred where the highest habitat selection classification intersected the lowest survival classification. Although ranks that included low predicted survival independent of selection can indicate low nesting productivity, we considered high relative selection coupled with lowest ranked survival to represent maladaptive habitat selection and potential ecological traps, where high relative selection coupled with low survival would result in a high incidence of nest failures and low subsequent reproduction. All spatially explicit mapping results are available in the USGS ScienceBase digital repository (https://doi.org/10.5066/P9TE06L4; Coates et al., 2020).

RESULTS

We located and determined fate of 786 sage‐grouse nests across 12 field sites (Appendix S1‐Table S1). Of all nests, 364 hatched successfully, and 422 failed. Females with nests were primarily adults (age ≥ 2 years at time of nest initiation = 80%). Most nests were initiated between late March and early June (95% CI = 6 April–1 June).

Nest site selection model

Variable reduction analysis resulted in 30 initial environmental metrics for the nest RSF (Appendix S4‐Table S4). The initial nest RSF model included each of these metrics, with associated parameters for interaction and quadratic effects (Table S4). Parameters that were considered uninformative (|Z| < 1.44 after inclusion in the full model) included interaction terms for distance to wet meadow, distance to cropland, elevation, non‐sagebrush shrub, and roughness; non‐linear, second‐order term for density of springs; and main effects for distance to cropland, NDVI, distance to water body, surface curvature, and density of small roads. The final reduced model had strong out‐of‐sample predictive capabilities (Spearman's ρ = 0.997, R 2 = .95, β predict = 0.72). The strongest predictors of habitat selection were indicated by clear shifts in the distribution of used nest locations relative to availability for a given predictor (Figure 1; Table S5‐1). The UHC plots (Figure 1), illustrating observed habitat falling within model‐predicted 95% confidence envelopes, indicate model‐predicted habitat was consistent with observed nest locations across multiple habitat predictors (Fieberg et al., 2018).

FIGURE 1

Used habitat calibration plots based on strongest effects estimated from a hierarchical spatial model of greater sage‐grouse nest site selection in California and Nevada, United States, fit to data from 2009 to 2017, with data from 2018 used for validation. Solid green lines represent empirical kernel density distributions of used habitat from sage‐grouse nest locations. Blue dotted lines represent density distributions of availability derived from random location measurements. Gray shaded areas represent 95% confidence intervals of model predictions of used habitat; empirical distributions falling within the shaded region indicate a well‐calibrated model (Fieberg et al., 2018). Predicted median used and available values are shown by green and blue dotted vertical lines, respectively. Results of the most influential metric and spatial scale are demonstrated for the most important landscape predictors of nest site selection, including (a) herbaceous cover, (b) big sagebrush cover, (c) residual litter, (d) sagebrush height, (e) density of springs, (f) density of intermittent streams, (g) low sagebrush cover, and (h) elevation

Female sage‐grouse selected nest locations with taller sagebrush plants (scale 2) and increased cover of big sagebrush (scale 3), herbaceous plants (scale 4), and litter (scale 1), while accounting for proximity to leks (Figure 1; Tables S5‐1 and S5‐2). Greater spring densities (scale 4), intermittent stream densities (scale 4), proportions of low sagebrush (scale 4), and greater elevations (scale 2) also had positive relationships with nest site selection (Figure 1; Tables S5‐1 and S5‐2). However, proportion of pinyon‐juniper (scale 4), non‐sagebrush shrub (scale 3), bare ground (scale 3), and compound topographic index (scale 3), negatively influenced nest site selection (Tables S5‐1 and S5‐2; Figure S5). Additional predictors retained in the model had less influence (Table S5‐1; Figure S5). The final model supported non‐linear relationships between nest site selection and intermittent stream densities, low sagebrush, and sagebrush height (Table S5‐2). Habitat functional response interactions with RR were supported for cover of big sagebrush, herbaceous plants, and bare ground (Table S5‐2), with weaker effects for agricultural cropland and transformed aspect (Tables S5‐1 and S5‐2).

Nest survival model

We estimated a 38 day nest survival probability of 0.25 (0.16, 0.35), which varied annually (Appendix S6‐Table S6). Variable reduction analysis indicated that fewer predictors (n = 10) influenced nest survival (Appendix S4‐Table S4). The final frailty model indicated strong increased hazard to sage‐grouse nests (reduced survival probability) with earlier nest initiation dates, increased proportion of CBA, and increased litter cover (Table 1; Figure 2). Predicted hazard was reduced with increased shrub cover immediately surrounding nests, greater landscape ruggedness (scale 4; Table 1; Figure 2), and reduced distance to wet meadows (Table 1; Figure 2).

TABLE 1

Posterior summaries of fixed effects estimated from a Bayesian hierarchical shared frailty model of greater sage‐grouse nest survival in California and Nevada, United States, fit to data from 2009 to 2017. Negative coefficients indicate reduced hazard (i.e., increased survival probability) , while positive values indicate increased hazard (reduced survival probability). Effects are sorted by relative influence (|mean[β]|/SD[β])

Variable	Mean	2.5th	Median	97.5th	|Mean|/SD	Influence on survival
% Total shrub cover (scale 2)	−0.185	−0.310	−0.184	−0.061	2.898	+
Roughness (scale 4) a	0.600	0.195	0.597	1.023	2.852	+ a
Initiation date	−0.149	−0.259	−0.149	−0.040	2.691	+
Cumulative burned area (scale 2)	0.138	0.016	0.138	0.261	2.212	−
% Litter (scale 3)	0.115	0.006	0.114	0.226	2.040	−
Proximity to wet meadow	−0.924	−1.950	−0.903	0.003	1.852	+

Scale 1: radius = 75 m (area = 1.8 ha).

Scale 2: radius = 167.9 m (area = 8.7 ha).

Scale 3: radius = 439.5 m (area = 61.5 ha).

Scale 4: radius = 1,451.7 m (area = 661.4 ha).

RR = Index of resistance and resilience, regional scale = 21,584 m (area = 146,574 ha).

Roughness was fit as an exponential decay function to allow its effect to subside at high values, rather than increase linearly. The coefficient of this effect is thereby reversed, indicating reduced hazard (increased survival) with increasing roughness.

FIGURE 2

Modeled effects for important seasonal and landscape predictors on cumulative 38 day nest survival for sage‐grouse nests in Nevada and California, United States during 2009–2017. Nest survival was modeled using a Bayesian hierarchical shared frailty model, and each predictor's effect was plotted conditional on all other predictors’ existence at their mean values at sage‐grouse nests

Predicted survival from the model was consistent with independent testing data (Bayesian p = .62). The deviation between observed and predicted cumulative survival rates for the testing data was centered around 0 (median = 0, interquartile range = [−0.004, 0.009]), indicating minimal bias. The model predicted greater cumulative survival at nests that were successful in the independent testing data (0.31) as opposed to nests that failed (0.25).

Nesting selection and survival mapping

Spatial variation in nest site selection and risk of nest failure indicated that nest site selection had fitness consequences for nesting sage‐grouse. Selection of greater herbaceous cover resulted in increased risk of nest failure in wildfire‐affected areas where annual grass was the dominant cover type (effect of % CBA; Figure 2). In contrast, selection of greater sagebrush cover and height, especially at higher elevations, resulted in increased nest survival, based on effects of shrub cover (correlated with sagebrush cover) and topographic roughness (correlated with elevation; Figure 2). Sage‐grouse habitat selection predictions indicated substantial selected habitat in the northwest, northeast, and south‐central sections of the study area (Figure 3a). Survival predictions varied within habitats selected by sage‐grouse, and strongest selection did not always coincide with high survival (Figures 3b and 4a). The combination of reclassified surfaces for nest selection and survival demonstrates the modeled gradient of nesting productivity (Figure 4a), with highest habitat class ranks (shades of blue; adaptive selection) comprising ~15% of mapped habitat, middle classes (shades of yellow; neutral) comprising ~69%, and lowest classes (shades of red; maladaptive selection) comprising 16%. Below‐average survival was predicted for 58.1% of habitat occurring within 17 km of active leks (Figure 4a). An evaluation of lowest productivity classes demonstrated the extent of burned area and invasion of annual grass influencing areas occupied by sage‐grouse (Figure 4b), where ~51% of the area in these classes was associated with recent burns and/or ≥10% annual grass cover (high/highest selection coupled with low/lowest survival; classes 12–13; 15–16 per Figure 4a).

FIGURE 3

(a) Map of nesting habitat selection scores predicted from a resource selection function (RSF) developed from sage‐grouse nest locations. Nest site selection was modeled using a generalized linear mixed model of used and random locations in a Bayesian modeling environment, and the midpoint of coefficient conditional posterior distributions were used for prediction. Continuous values were reclassified and ranked using a percent isopleth approach with respect to observed nest locations. (b) Map of cumulative 38 day nest survival predicted from a Bayesian hierarchical shared frailty model of sage‐grouse nest fates. The midpoint of coefficient conditional posterior distributions of 38 day nest survival were used for prediction at each 30 m pixel across the landscape. The map shows predicted survival reclassified based on the 25th percentile of all values at failed nests (lowest class), mean value of all nests, and 75th percentile of successful nests. Fate and location data were gathered from ground radio telemetry monitoring of sage‐grouse females in Nevada and California, United States, 2009–2017

FIGURE 4

(a) Ranked index of model‐projected nest site selection integrated with nesting productivity (i.e., nest survival), demonstrating the spatial distribution of adaptive versus maladaptive habitat selection at each 30 m pixel. Hierarchical models of nest selection and survival were fit to landscape covariates within a Bayesian modeling framework in Nevada and California from 2009 to 2017 to develop spatially explicit information about nest site selection and survival consequences across the landscape. Habitat was separated into 16 classes ranking from high (1) to low (16). Habitat ranked highest where the top nest selection and survival classes intersected (adaptive selection), whereas the lowest rank occurred where the top nest selection class intersected with the lowest survival class (maladaptive selection, or potential ecological trap). Areas of greatest sage‐grouse occupancy are highlighted within 17 km of all active leks in the study area. (b) Map of burned areas from a cumulative burned area (CBA) model representing burn scars and annual grass cover >10% in the same study region. Burned areas were represented based on wildfires occurring between 1984 and 2017 that had not recovered to 20% sagebrush cover, or had otherwise experienced a permanent state transition to annual grass based on simulated sagebrush regrowth across varying soil temperature and moisture regimes

DISCUSSION

Semi‐arid shrubland ecosystems are undergoing accelerated rates of change driven by anthropogenic land use and disturbance (Coates, Ricca, et al., 2016; Knick et al., 2003), and are vulnerable to desertification resulting from changing climate, intensifying wildfire cycles, and invasion of exotic annual grasses (Pellant, Abbey, & Karl, 2004; Rigge et al., 2019). Habitat alteration and degradation over relatively short time periods can lead to population decline and increased extinction risk, particularly for specialist wildlife species that are less able to adapt to rapid change (Clavel et al., 2011). In this respect, geographic availability of the ecological niche space occupied by certain species is likely diminishing in response to intensifying wildfire regimes and the spread of exotic grass species, as measured by declines in habitat that contribute positively to reproductive success and population growth (Coates, Ricca, et al., 2016; Rockweit et al., 2017). Thus, altered wildfire regimes fueled by invasive annual grasses may contribute to large‐scale functional homogenization (Clavel et al., 2011) in affected environments over relatively broad time scales.

We used a novel methodological approach combined with a long‐term, spatially expansive dataset to integrate hierarchical resource selection and frailty modeling with landscape predictors at multiple spatial scales. The results delineated a complex mosaic of adaptive and maladaptive habitat selection, implicating a source–sink reproductive landscape for imperiled sage‐grouse breeding in Great Basin shrubland ecosystems. Specifically, we found that sage‐grouse frequently establish nests in areas that were highly selected for, but their probability of nest success was relatively poor. In other words, this misalignment in selection and survival patterns delineate areas on the landscape that represent ecological traps (Battin, 2004; Delibes, Gaona, & Ferreras, 2001), which were largely a result of past wildfires that have reshaped vegetation communities within sagebrush ecosystems (Blomberg, Sedinger, Atamian, & Nonne, 2012; Chambers et al., 2014; Reisner, Grace, Pyke, & Doescher, 2013). These changes appear to be occurring at a pace that exceeds that of sage‐grouse’ ability to adapt by modifying their habitat selection behaviors. For example, selection for herbaceous plant characteristics (i.e., herbaceous and residual cover) did not consistently result in improved nest survival, as sage‐grouse preferences for greater herbaceous cover often led to nests being established in previously burned areas that were associated with reduced nest survival. Within our study, 23% of nest locations (n = 180) occurred in areas that had recently burned without adequate recovery time These sites were typically characterized by a significant annual grass component coupled with reduced shrub cover, resulting in negative fitness consequences. In fact, our models predicted that 58.1% of classified habitat within 17 km of active leks (~53,440 km2) would yield below‐average nest survival. Furthermore, 51% of habitat classes characterized by high selection but lowest survival in this region were associated with burned areas that had not recovered and/or comprised ≥10% annual grass cover. This evidence of maladaptive habitat selection (Remeš, 2000) may be explained in part by rigid patterns of breeding and nest site fidelity behavior (Fischer et al., 1993; Holloran & Anderson, 2005; Schroeder & Robb, 2003), wherein sage‐grouse demonstrate strong philopatry with respect to previous breeding and nesting areas that no longer support high rates of successful nesting following habitat alterations.

Breeding site fidelity can influence female sage‐grouse habitat selection patterns on long‐ and short time scales. Individual nest site fidelity has been shown to have a strong short‐term influence on nest site selection, with females most frequently selecting nests within 400–800 m of prior year nest locations (medians reported, overall range of distances = 0.005–33 km; Fischer et al., 1993; Holloran & Anderson, 2005; Schroeder & Robb, 2003). Sage‐grouse also demonstrate natal lek fidelity (Dunn & Braun, 1985) and select nest locations near breeding leks (Coates et al., 2013; Holloran & Anderson, 2005), suggesting that fidelity can be passed through generations, leading to more long‐term influences. Site fidelity can be advantageous in relatively static environments, where experienced females enhance their knowledge of the environment, such as locations of water, prevailing predator communities, and other features that previously resulted in a successful nest (Schroeder & Robb, 2003). However, novel threats in habitats altered by wildfire may outweigh advantages of site fidelity (Knick & Rotenberry, 2000), potentially creating ecological traps (Battin, 2004; Robertson et al., 2013) where individuals either fail to adapt to a changing environment or are attracted to postburn pulses in grassy vegetation (Foster, 2016).

Areas impacted by wildfire that lose shrub communities and become dominated by annual grasses (Lockyer, Coates, Casazza, Espinosa, & Delehanty, 2015) may not provide adequate concealment cover to protect nests from visual predators such as ravens (Coates & Delehanty, 2010; O’Neil et al., 2018), yet behaviorally or genetically entrained nest site fidelity may lead female sage‐grouse to nest in these depauperate sites, where nest survival is apparently reduced (Foster et al., 2019). Continued increases in wildfire frequency and severity, leading to further encroachment of annual grass cover throughout the Great Basin, likely will result in reduced population recruitment in heavily impacted areas by way of reduced nest survival (Blomberg et al., 2012; Coates, Ricca, et al., 2016; Foster et al., 2019). To further understand patterns in site fidelity post hoc, we evaluated data from females tracked across multiple breeding (n = 150) and observed substantial nest site fidelity (median distance between consecutive year nests = 604 m; 95% CI = 55–6,076 m). Furthermore, the range of distances did not appear to be influenced by previous burns as indicated by our CBA model (two‐sided non‐parametric permutation test, Z = −0.031, p = .614). In addition, we documented two females that nested in year 2017 prior to a significant wildfire and returned to nest within the burned area the following year, 230 and 1,820 m from their previous nest locations, respectively. These observations combined with our nest survival modeling results imply little evidence of adaptive nest site selection in response to wildfire (see also Foster et al., 2019). Instead, sage‐grouse demonstrate rigid philopatric behaviors conducive to a system of which interannual variation in shrub cover is typically low, and such behaviors of site fidelity appear to be detrimental to sage‐grouse reproduction under accelerated fire regimes. It is also important to recognize that in the absence of wildfire, sage‐grouse prefer contiguous tracts of sagebrush that offer overhead nesting cover (Connelly et al., 2000). Thus, our findings largely indicate that philopatric behaviors for areas that offered overhead concealment from shrubs prior to wildfire events is maladaptive and reduces reproductive success because burned areas frequently transition to invasive annual grass with reduced shrub cover. This process provides one possible mechanism for previously documented long‐term patterns of population declines of sage‐grouse associated with increased burned areas (Coates, Ricca, et al., 2016).

Negative impacts of accelerated wildfire on wildlife species and their habitats will likely vary depending on biophysical factors, such as localized climate and soil conditions that influence an ecosystem's resilience to disturbance and resistance to invasion from exotic species (RR; Chambers et al., 2014). Our results supported functional responses in resource selection by sage‐grouse that varied across an RR gradient that we quantified at the scale of sage‐grouse subpopulations in our study. For example, selection for big sagebrush and herbaceous cover declined with increasing RR, selection for south‐facing slopes increased with RR, and avoidance of agricultural fields decreased (Table 1). Habitat functional responses suggest that sagebrush and herbaceous cover may be more limiting in areas of low RR. As RR increases, typically with increases in elevation and average precipitation, the availability of big sagebrush and herbaceous cover subsequently increases, resulting in selection patterns more proportional to regional availabilities. The combined effects of these interactions and vegetation community effects reveals an overall preference for higher regional RR, comprising higher elevation, mesic landscapes with greater shrub cover and wet meadow availability that contribute positively to sage‐grouse nest selection and survival at broad landscape scales. Even in these more resilient sites, however, the accelerated pace of the wildfire‐grass cycle creates spatiotemporal lags between relatively slower shrubland recovery and faster demographic responses of nesting sage‐grouse to disturbance (Coates, Ricca, et al., 2016; Ricca & Coates, 2020), which likely result in increasing suboptimal habitat use by sage‐grouse and other sagebrush‐obligate species. Shrubland recovery may be further hindered by warmer temperatures, severe drought and intensifying El Niño events that contribute to positive grass‐fire feedback loops (Balch et al., 2013). Although warmer and drier sites may support sagebrush specialist species such as sage‐grouse during favorable precipitation years (Coates et al., 2018) these habitats generally characterize low overall productivity and are at risk of unsustainable rates of transition to annual grass states (Coates, Ricca, et al., 2016). Current and projected climate trends indicate that conditions will continue to be favorable for such state transitions (Cai et al., 2014; Wang et al., 2019), suggesting that sagebrush‐shrub communities and the species that depend on them face a challenging future. These circumstances underscore the need to quantify sensitive habitats in terms of their potential to promote reproductive success, population growth, and persistence (e.g., approximation of the fundamental niche), thereby facilitating evaluation of range contraction rates under varying conservation and land management scenarios.

We identified nesting productivity consequences of habitat selection that have important implications for conservation and management of habitat in shrubland ecosystems. Sage‐grouse populations use high productivity habitat types (i.e., “source” habitats) that compensate for reduced performance in marginal habitat, indicating critical value in preserving source habitats to conserve populations over the long term. Second, identifying sink habitat types, as predicted by maladaptive selection patterns, allows for additional exploration and analysis in these areas to determine: (a) the extent to which these habitats are being occupied, (b) the specific mechanisms leading to poor performance in these areas, and (c) the potential to remediate. For example, according to our classification, a substantial proportion (~58%) of nest locations apparently consists of environmental attributes that predict moderate to high probability of nest selection with low expected survival rates and probably represent a form of ecological trap. In some cases, modifications or restoration, such as sagebrush planting projects in burned areas, may be possible to improve marginal habitats to be more consistent with the higher productivity habitats, leading to improved population performance. However, external factors that were not considered in our models such as predators, disease, or genetic limitation may interact with those factors that we considered in explaining the reduced nest performance in these areas. Our analytical approach can identify ecological traps and possible sink habitats across other life‐stages (e.g., adult and juvenile selection and survival), and may be extended to any species where data on selection and survival can be collected. Finally, parameter estimates from selection and survival models, coupled with spatially explicit estimates of sagebrush recovery (Ricca & Coates, 2020), inform generalizable planning tools (Ricca et al., 2018) for efficiently prioritizing localized postfire restoration efforts across large spatial extents to recover losses in breeding habitat suitability.

Our study stresses the importance of understanding relationships between species’ distribution and their fundamental ecological niche in the context of altered disturbance regimes and rapid global change. Landscape‐level species conservation efforts will increasingly require knowledge of fitness or population performance within occupied distributions to identify areas at risk of range contraction. We demonstrated spatially expansive mismatches between the current distribution of sage‐grouse and population performance, where former habitats fail to support high nesting productivity as wildfire becomes more frequent and severe. Adaptive mechanisms such as dispersal may be outpaced by current rates of habitat alteration for sage‐grouse, and this may be true of many other species dependent on sagebrush within these ecosystems. The implications of maladaptive behaviors among species within environments modified by anthropogenic disturbance and development likely also lead to poor demographic performance and perhaps population declines. We caution that representations of the ecological niche may be misleading when measures of habitat selection or occupancy are not coupled with demographic responses in species distribution modeling, especially for species prone to ecological traps. We encourage future investigations spanning multiple life stages to measure the impacts of marginal habitats on population growth rates or reproductive success, thereby facilitating characterizations of gains or losses to a species geographic realization of the fundamental niche under dynamic, altered disturbance regimes.

PERMITS

Sage‐grouse captured and marked in this study were done so with permission from the Nevada Department of Wildlife and California Department of Fish and Wildlife, and in accordance with USGS Western Ecological Research Center Animal Care and Use Committee permit WERC 2015‐002, and Idaho State University Institutional Animal Care and Use Committee Protocol 770.

CONFLICT OF INTEREST

The authors declare no conflicts of interest. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the US Government.

AUTHOR CONTRIBUTION

S.T.O., P.S.C., B.E.B., M.A.R., S.P.E., S.C.G., and D.J.D. conceived the study and contributed critically to interpretation of results and the writing process; S.T.O., P.S.C., B.E.B., and M.A.R. analyzed the data and wrote the manuscript; all co‐authors gave final approval.

Appendix S1

Click here for additional data file.

Appendix S2

Click here for additional data file.

Appendix S3

Click here for additional data file.

Appendix S4

Click here for additional data file.

Appendix S5

Click here for additional data file.

Appendix S6

Click here for additional data file.