Introduction

The degradation, fragmentation, and loss of native shrub-steppe ecosystems, as well as the concurrent decline of wildlife populations that depend on them, are among the most pressing issues facing land managers across western North America (Davies et al., 2011; CFR, 2015a). In sagebrush ecosystems of the Great Basin, distribution and abundance of pinyon (primarily Pinusmonophylla) and juniper (primarily Juniperus osteosperma) woodlands (hereafter, “pinyon-junipeR&Rdquo;) has increased dramatically (i.e., >150%) since European settlement (Miller et al., 2008) owing to changes in land-use practices (Romme et al., 2009), climate (Miller and Wigand, 1994; Romme et al., 2009), and disturbance regimes (Miller and Rose, 1999). For the purposes of this paper, we define encroachment to include both expansion (establishment of pinyon-juniper into areas previously devoid of trees) and infill (increasing closure of previously sparse pinyon-juniper canopies), as modified from Miller et al. (2013). Although pinyon-juniper is a native component contributing to landscape heterogeneity in the Great Basin and some encroachment may stem from natural recovery of pinyon-juniper woodlands previously cleared by European settlers (Romme et al., 2009), the overall current rate of encroachment is profoundly influencing contemporary sagebrush ecosystem processes (Miller et al., 2005; Davies et al., 2011). Accordingly, a variety of management actions have been directed toward decreasing the rate of pinyon-juniper expansion into sagebrush plant communities (Tausch et al., 2009), and many of these actions are the focus of studies presented in this volume of Rangeland Ecology & Management.

Habitat selection theory generally predicts that animals will occupy areas that optimize their fitness (i.e., survival and reproduction; Rosenzweig, 1981; Morris, 1989), and identifying mechanisms that link selection of environmental features with fitness is a lynchpin for implementing management actions that improve habitat suitability (Van Horne, 1983; Aldridge and Boyce, 2007; Casazza et al., 2011). Yet individuals do not always make decisions that maximize their fitness as they move through their environment, perhaps owing to variability in howthey perceive environmental cues (Matthiopoulos et al., 2015) or respond to intraspecific competition for limited resources (Fretwell and Lucas, 1970). Moreover, maladaptive selection can lead to the formation of ecological traps, which can be defined as environments that provide attractive cues that yield lower survival and/or reproduction output, thereby decoupling an individual's perception of the habitat's fitness-harming traits (Robertson and Hutto, 2006). Directly linking an organism's fitness to its decisions becomes even more daunting and complex when these sources of variation are coupled with logistical difficulties in obtaining sufficient data on movement processes, environmental features, and fitness repercussions in a time-dependent fashion for multiple individuals of a species that occupy remote and vast landscapes (Morales et al., 2010).

Greater sage-grouse (Centrocercus urophasianus, hereafter sage-grouse) require large continuous areas of sagebrush-dominated ecosystems for population persistence (Knick et al., 2013), and this species is considered an indicator species for the health of sagebrush ecosystems because they require distinct ecological states to fulfill their diverse life history requirements at large spatial scales (Rowland et al., 2006; Hanser and Knick, 2011). Populations of sage-grouse have declined concomitantly with the loss and fragmentation of sagebrush ecosystems that now occupy slightly more than half of their former range (Schroeder et al., 2004; Miller et al., 2011). In large parts of the Great Basin, encroachment of pinyon-juniper has been identified as a primary threat to sage-grouse populations (CFR, 2015a) by contributing to fragmentation of continuous expanses of sagebrush and accelerating a positive feedback between wildfire and invasive annual grass (the other primary threat in the Great Basin) that often eliminates and replaces sagebrush (Brooks et al., 2004; Balch et al., 2013; Chambers et al., 2014a).

Several studies have documented strong avoidance of pinyon-juniper by sage-grouse at multiple spatial scales and across different grouse life history stages (Doherty et al., 2008; Atamian et al., 2010; Casazza et al., 2011; Knick et al., 2013) even at relatively low density (e.g., b4% canopy cover; Baruch-Mordo et al., 2013). Importantly, avoidance of pinyon-juniper by sage-grouse can have population-level consequences to brood survival (Casazza et al., 2011) and lek persistence (Baruch-Mordo et al., 2013), and can lead to genetic isolation (Oyler-McCance et al., 2005; Oyler-McCance et al., 2014). Different levels of pinyon-juniper cover (e.g., sagebrush dominant to pinyon-juniper woodland) may vary in their effects on sage-grouse behavior and population dynamics. For example, important resources to sage-grouse such as food and concealment cover decrease disproportionately as the percent of pinyon-juniper overstory increases (Bates et al., 2005; Miller et al., 2005; Miller et al., 2011). Additional tall vertical structures (such as trees) that provide perching and nesting habitat in an otherwise flat landscape can increase risk of avian predation (Coates et al., 2014a; Howe et al., 2014), which sage-grouse may perceive as a threat that changes with the density of trees on the landscape.

For management purposes, continuous encroachment of pinyon-juniper is often categorized into three transitional phases (i.e., Phase I, II, and III) indicating the dominant vegetation influencing ecological processes (Miller et al., 2005; Miller et al., 2013). For plant community structure, Phase I is characterized by relatively low pinyon-juniper canopy cover and overall dominance of sagebrush and associated perennial grasses. During Phase II, herbaceous sagebrush understory begins to thin significantly and become codominant with pinyon-juniper, while Phase III is characterized by dominance of pinyon-juniper and little to no herbaceous sagebrush understory. For sage-grouse, these transitional phases may elicit different demographic and behavioral responses that have important implications for management of pinyon-juniper and the mechanisms underlying degradation of sage-grouse habitat suitability.

A deeper understanding of relationships between evolved environmental cues that influence sage-grouse behavioral choices (e.g., to avoid or select an area) and the demographic consequences (e.g., mortality) of those choices could benefit conservation actions for sage-grouse populations. In the case of restoring sagebrush ecosystems by removing recently established pinyon-juniper trees, identifying these complex patterns and processes in relation to transitional phase is especially important because specific management actions designed to improve habitat suitability (e.g., removal of all trees in lower-density stands vs. thinning of higher-density stands) may elicit different behavioral responses from sage-grouse, which in turn might yield unique consequences for sage-grouse fitness components and management efficacy. Furthermore, quantifying thresholds within phases where management actions achieve desired goals (e.g., maintain or increase fitness) is paramount for effective conservation planning (Baruch-Mordo et al., 2013). Identifying linkages between ecological mechanisms driving both the “how” and “why” (e.g., selection and fitness) is an integral, yet largely unknown, part of the conservation planning process.

Typical quantitative approaches for linking sage-grouse habitat selection with fitness consequences use available habitat as a covariate in traditional survival (Aldridge and Boyce, 2007; Aldridge and Boyce, 2008; Casazza et al., 2011) or regression-type analyses (Baruch-Mordo et al., 2013). These approaches treat survival probability as a function of habitat in a time-independent or static manner and often link overall estimates of resource selection at the individual level to demographic performance at the population level without explicit consideration of the frequency or timing of encounters with landscape features (such as pinyon-juniper among different phases). Time-dependent analyses, in contrast, can be more computationally complex, yet they can more clearly establish a linkage between how individual sage-grouse respond behaviorally and demographically to specific encounters of pinyon-juniper among different phases. Studies that integrate both time-dependent and independent analyses can better identify mechanisms driving selection and survival.

Herein, we employ a novel two-stage Bayesian modeling approach to link estimated probability of avoidance of different pinyon-juniper cover classes with concomitant changes in annual probability of survival, while accounting for confounding factors and uncertainty in parameter estimation. Most importantly, this Bayesian approach incorporates the range of behavioral heterogeneity among individual sage-grouse, which allows for uncertainty in behavioral choices and their consequences for fitness. To provide target values for conifer removal to reach survivability, we then carried out a post hoc analysis that estimated survival directly as a function of time-dependent use of pinyon-juniper cover under conditions with varying levels of primary plant productivity. We focus on sage-grouse within the Bi-State Distinct Population Segment (Bi-State DPS) along the central border of California and Nevada, which has been recently ruled unwarranted for protection under the Endangered Species Act. Pinyon-juniper encroachment has been identified as the primary threat to the Bi-State DPS (CFR, 2015b), and the listing decision was informed in large part by planned implementation of large-scale treatments (thinning or removal) of encroaching pinyon-juniper across thousands of acres of habitat (Bi-State Action Plan, 2012).

Study Area

The Bi-State DPS comprises 18 325 km² split along the border of Nevada and California at the interface of the Sierra Nevada Mountains to the west and the Great Basin to the east (Fig. 1A; lat 119°11′1.94′′N, long 38°6′30.80′′W). We collected sage-grouse data from seven subpopulations: Bodie Hills (BH), Long Valley (LV), Parker Meadows (PM), Pine Nut Mountains (PN), Mount Grant (MG), Desert Creek/Fales (DF), and White Mountains (WM), which comprise three distinct subregions: northern (consisted of PN), central (consisted of BH, LV, PM, MG, and DF), and southern (consisted of WM) (see Fig. 1A). The subregions were delineated based on distinct differences in vegetative communities, topography, and sage-grouse population genetics (Oyler-McCance et al., 2014). We considered subregions in this analysis to be separate because they represent empirically supported subdivision of the populations within the Bi-State area.

Figure 1.

Greater sage-grouse telemetry data across study subregions separated by A, north, central, and southern region, and B, mapped pinyon-juniper classes (1, 2, and 3) at the 1-m resolution within the Bi-State Distinct Population Segment during 2003–2005 and 2011–2013.

The Bi-State area was topographically diverse, and several mountain ranges separated the northern and southern ends. Elevations ranged from 1 386 to 4 344 m with numerous rugged mountain ranges intermixed with broad valleys. General vegetation communities predominantly consisted of mountain big sagebrush (Artemisia tridentata vaseyana) interspersed with areas of low sagebrush (A. arbuscula), and Wyoming big sagebrush (A. t. wyomingensis). Silver sagebrush (A. cana) and basin big sagebrush (A. t. tridentata) occurred locally. Other common shrub species included snowberry (Symphoricarpos spp.), currant (Ribes spp.), bitterbrush (Purshia tridentata), green rabbitbrush (Chrysothamnus viscidiflorus), rubber rabbitbrush (Ericameria nauseosa), and Mormon tea (Ephedra viridis). Primary grass species included needle grass (Hesperostipa comata), squirrel tail (Elymus elymoides), and Indian rice grass (Achnatherum hymenoides). Cheatgrass (Bromus tectorum) was present but uncommon. Dominant forbs included phlox (Phlox spp.), lupine (Lupinus spp.), buckwheat (Eriogonum spp.), and hawksbeard (Crepis spp.) (Kolada et al., 2009). Singleleaf pinyon (Pinus monophylla) and Utah juniper (Juniperus osteosperma) woodlands occurred at elevations of 1 850–3 000m.

Methods

Linking Avoidance of Pinyon-Juniper Encroachment to Survival

We conducted a novel, two-stage Bayesian analysis to estimate the effects of pinyon-juniper within sagebrush ecosystems on sage-grouse habitat selection and survival. Our analysis quantified avoidance of pinyon-juniper across each cover class (CC1, CC2, and CC3) using mixed effects logistic regression models that contrasted patterns of used versus available resources (Boyce et al., 2002; Manly et al., 2002; Johnson et al., 2006), and then linked those patterns (particularly avoidance across each cover class) to sage-grouse survival. Specifically, we evaluated whether sage-grouse tendencies to avoid pinyon-juniper cover influenced a component of their fitness (i.e., survival), ultimately informing the effects of pinyon-juniper on habitat suitability. Specifically, the two-stage process consisted of 1) modeling an “avoidance index of pinyon-junipeR&Rdquo; using predicted negative log odds ratios while accounting for potentially confounding variables and deriving posterior distributions of estimated avoidance parameters, β (slope coefficients), for each individual sage-grouse; and 2) fitting the sampled parameters of individual-level avoidance as covariates for survival using a frailty analysis. This two-stage Bayesian model ran in parallel at each stage and used full posterior distributions to account for uncertainty in individual estimates of avoidance, allowing the empirical evaluation of whether or not sage-grouse that avoided pinyon-juniper were more likely to survive. Detailed explanation for each stage follows.

We employed a used-available design to evaluate avoidance of pinyon-juniper cover by sage-grouse similar to procedures described in Casazza et al. (2011). Specifically, we employed a Design II approach (Manly et al., 2002), meaning “use” of pinyon-juniper was measured at the individual level, whereas “availability” was measured at the population level (i.e., Erickson et al., 2001). Sage-grouse telemetry locations were imported into a geographic information system (GIS), and all predictor variables were extracted. To characterize the availability of pinyon-juniper cover, we extracted the same environmental attributes at random locations (Manly et al., 2002). We calculated a minimum convex polygon (MCP) of the combined sage-grouse relocations for each of the subpopulations to represent the boundaries of available habitat. We generated a random location for each used location, which resulted in a range of random locations within each monitored subpopulation as 751–3 279 (mean=1 501). We considered this an appropriate number of random locations to allow us to characterize availability and yet be manageable computationally using our Bayesian analytical approach. Before our modeling approach, we carried out a prerequisite analysis to identify the most parsimonious ecological scale to represent land cover variables (e.g., CC1, sagebrush) using single variable models. Specifically, we compared evidence among single variable models for each variable consisting of the three different spatial scales, and carried forward the scale that best represented the variable on the basis of the lowest value of Deviance Information Criterion (DIC; Spiegelhalter et al., 2002) and the posterior probability of nonzeroness derived from a stochastic search variable selection method (George and McCulloch 1996). DIC is particularly useful in this circumstance because of its simplicity, ease of calculation using MCMC samples, and similarities to AIC for generalized linearmodels constructed within a Bayesian framework (Hooten and Hobbs, 2015).

Figure 2.

Maps depicting A, sagebrush cover and B, resilience and resistance classes, defined by soil and temperature models (Maestas et al., 2016)within the Bi-State Distinct Population Segment during 2003–2005 and 2011–2013.

During the first stage, we modeled the pinyon-juniper cover classes (CC1, CC2, and CC3)while accounting for potential confounding variables. We developed separate models for each cover class so that effects specific to each cover class could be identified more clearly. All other variables remained the same in each model. We followed the procedures for modeling resource selection with random effects described in Gillies et al. (2006), where the mixed effects logit model, g(x), took the form:

denoting a population level intercept (β0),with coefficients β for covariates pinyon-juniper (PJ; continuous, percent cover), sagebrush cover (SB; continuous, percent cover), and productivity of area (PR; categorical, more or less productive). As a diagnostic test, we estimated correlation coefficients to reduce potential effects of multicollinearity between PJ and SB covariates. Because variables did not show evidence of covarying (r ≤ |0.65|),we retained both covariates in the model. We fit random effects to account for repeated measures within individual (bird[ςi]), spatial (subregion [κ_j]), and temporal (year [η_k]) correlation (Gillies et al., 2006). Subscripts h, i, j, and k reference telemetry location, bird, subregion, and year, respectively. x_hijk represents a data vector consisting of (PJ_hijk, SB_hijk, PR_hijk). Further definitions of parameters and their prior specifications are listed in Table 1. Negative log odds ratios indexed avoidance while positive log odds ratios indexed selection. We did not exponentiate predicted log-odds to estimate resource selection probability functions because availability was estimated and not censused (Boyce et al., 2002). All parameters were assigned objective priors. We followed procedures for the random-coefficients model described in Kéry (2010) with correlation between intercepts and slopes, which allowed estimation of different levels of pinyon-juniper avoidance among individuals described conceptually in Gillies et al. (2006). Thus, we derived posterior distributions for β_PJ for each bird (i), which were used as an index of avoidance behavior of pinyon-juniper (PJAI) and inputted into the second frailty model stage. In other words, this random effects structure allowed individual slope estimation for each bird to then relate to survival.

Table 1

Median estimates (with 2.5 and 97.5 percentiles in parentheses) of posterior distributions based on modeling the effects of pinyon-juniper cover classes (CC1, CC2, and CC3)with priors on three responses: 1) avoidance of each cover class (selection model); 2) influence of the magnitude of avoidance on survival (linked model); and 3) a time-varying effect of exposure to each cover class directly on survival (frailty model) of greater sage-grouse within the Bi-State Distinct Population Segment during 2003–2005 and 2011–2013.

In a simultaneous second stage, a frailty model (Halstead et al., 2012) was used to fit individual-level avoidance of each pinyon-juniper cover class (PJAI) on individual survival probability. Hence, avoidance was formally linked to survival by randomly sampling the posterior distribution of slope coefficients (β_PJ,i) for individual sage-grouse. The major advantage of this two-stage approach was to make inferences on survival as related to heterogeneity in the range of behavioral tendencies related to encountering pinyon-juniper among individual grouse, which aligns more with a biological evaluation of habitat suitability. Additionally, this analysis accounts for multiple levels of uncertainty inherent in parameter estimation such as 1) individual-level estimates of avoidance, 2) individual-level effects of avoidance on survival, and 3) population-level effects of avoidance on survival considering variation in avoidance among individuals. To prepare the data for this component of the analysis, we first developed encounter histories for sage-grouse using 10-day incremental relocation data. Ten days was a conservative measure to minimize the number of censored intervals (intervals missing relocation data) and align the telemetry sampling frequency with the scale of inference. We considered right censoring to be random and occurred under multiple circumstances, including when 1) the study period ended before death was classified; 2) the transmitters failed; or 3) sage-grouse were lost prior to follow-up monitoring. Thus, all individuals either died or were eventually right censored. We modeled survival as a continuous process but observed at the discrete 10-day interval. We assumed that mortality risk was independent among individuals considering that sage-grouse do not engage in social defense strategies. Importantly, numerous telemetry points were required to estimate slope coefficients for individual grouse and thereby precluded fitting avoidance as a time-varying covariate in this specific analysis (see Post Hoc Analysis). However, this analysis effectively evaluates whether or not the overall tendency of an individual to avoid pinyon and juniper influences survival. The model was expressed as:

denoting a change in unit hazard (UH; 10-day risk of mortality) of magnitude β for a unit change in the PJAI. Using UH, a cumulative hazard (CH) was calculated as:

where thirty-six 10-day intervals derived an annual survival parameter (S) as:

For both avoidance and survival models, we modeled each cover class separately and evaluated evidence of the cover class variable using 95% credible limits of the sampled posterior distribution. For illustrative purposes, we predicted the relative percent change in annual survival as a function of the relative probability of avoiding each cover class.

Post Hoc Survival Analysis

We carried out a post hoc survival analysis and fit sage-grouse associations with cover classes as time-varying covariates following similar procedures as Halstead et al. (2012). This analysis provided direct estimates of relationships between pinyon-juniper and survival that are often more useable for managerial decisions. For example, such estimates can be used to predict survival probabilities across various pinyon-juniper cover values and identify differences between cover classes. Additionally, continuous variables can be fit for each 10-day interval as time-varying individual covariates, allowing these variables to change across each interval as sage-grouse moved from one location to another. Such time-varying covariates also allow for a more thorough investigation of complex ecological relationships. Accordingly, we investigated how sage-grouse use of areas with varying levels of productivity (i.e., more versus less) altered the effect of pinyon-juniper on survival. Because we measured variables only at time intervals when sage-grouse were relocated, we implemented a first-order Markov process to impute data for covariates at intervals when sage-grouse were not relocated as described in Halstead et al. (2012). Because sage-grouse largely avoided areas with > 10% pinyon-juniper cover, we pooled CC2 and CC3 into one grouping (CCP) for this post hoc survival analysis.

We fit a frailty model in this analysis similar to the linked avoidance-survival analysis. However, because we modeled survival directly as a function of pinyon-juniper as an environmental variable, and without avoidance, we fit other environmental variables directly into the model that might otherwise confound the effects of pinyon-juniper. The model was expressed as:

where the unit hazard (UH) represented a constant hazard function, γ, with estimated coefficients for sagebrush (β_SB), transmitter type (β_TT, which represents β_VHF or β_GPS in accordance with type), month (β_M, which represents β_M1,..., β_M12 in accordance with month of 10-d time interval, h), and βPJ. To investigate the effect of PJ, we specified a shared frailty across the two levels of productivity, referenced as l. This specification allowed us to estimate differences in the PJ effects on survival between areas with low and high productivity. Sagebrush and transmitter type were included as covariates to account for potential confounding effects. Month was included because sage-grouse have been shown to have unequal survival probabilities across their annual lifecycle (Blomberg et al. 2013a, Moynahan et al. 2006). Subscripts h, i, j, and k reference 10-d time interval, bird, subregion, and year, respectively. Models were specified with random effects for subregion (κ_J) and year (η_K) to account for spatial and temporal intraclass correlation, respectively. Similar to the previous analyses, as a diagnostic testwe estimated correlation coefficients to reduce potential effects of multicollinearity between PJ and SB covariates. Variables did not show evidence of covarying (r ≤ |0.65|), so we retained them in the model. Definitions of parameters and their prior specifications are listed in Table 1. Priors for all parameters were selected to be uninformative. We evaluated evidence of each cover class variable across the two productivity classes using 95% credible limits of the sampled posterior distribution. For illustrative purposes,we predicted the cumulative annual survival for every unit increase in cover within each cover class by productivity class. Lastly, also for illustrative purposes, we mapped pinyon-juniper areas that are at highest risk to sage-grouse populations within the Bi-State DPS, based on our model results. All models were run on three chains of 10 000 iterations each following a burn-in period of 40 000 iterations and thinned by a factor of 5. Model convergence was assessed visually using history plots and the R-hat statistic (Gelman et al., 2004). We did not find a lack of convergence among any of the parameters monitored (maximum R-hat = 1.1). Posterior probability distributions for each model procedure were estimated using Program R version 3.1.1 (R-Core-Team, 2014) and package “rjags” (Plummer et al., 2015).

Results

We captured and marked 162 sage-grouse (129 females and 33 males) across the 6 study years (2003–2005 and 2011–2013), of which we attached 139 VHF and 23 GPS transmitters and obtained 6 942 and 3 565 relocations, respectively. We obtained 3 137, 6 619, and 751 locations from49, 95, and 18 sage-grouse within the northern, central and southern regions, respectively. Overall, the majority of transmitters (55%) were deployed within the central region of the Bi-State DPS, while the southern region only received VHF transmitters and had the least amount of representation (12%). Within the Bi-State DPS boundary, we estimated that CC1, CC2, and CC3 composed 33.2%, 11.6%, and 9.8% of all land cover, respectively (northern subregion = 36.7%, 13.3%, and 17.5%; central subregion= 35.7%, 12.6%, and 9.4%; and southern subregion = 27.1%, 9.2%, and 7.1% of CC1, CC2, and CC3, respectively). Within defined sage-grouse habitat within the Bi-State DPS boundary (BLM and USFS, 2015), we estimated CC1, CC2 and CC3 represented 46.2%, 8.9%, and 6.1% of land cover, respectively (northern subregion = 16.4%, 4.3%, and 6.2%; central subregion = 20.9%, 3.9%, and 2.0%; and southern subregion = 4.3%, 0.7%, and 0.2% of CC1, CC2, and CC3, respectively).

Modeling Pinyon-Juniper Avoidance and Effects on Survival

Based on the prerequisite analyses, the ecological scale of 1 451 m garnered the most support from the data (lowest DIC and highest SSVS) for the sagebrush covariate, while a 439-m scale garnered the most support for all pinyon-juniper cover classes (Table S1 in the online version at  http://dx.doi.org/10.1016/j.rama.2016.09.001). Thus, these scales were carried forward in subsequent modeling procedures. Each model in stage 1 of the analysis consisted of a cover class covariate (CC1, CC2, or CC3) and the sagebrush and productivity covariates as additive fixed effects. These variables were evidenced by the data and retained in the model on the basis of SSVS values (see Table S1 in the online version at  http://dx.doi.org/10.1016/j.rama.2016.09.001). Posterior distributions of parameters for potentially confounding variables were reported in Table 1. Sage-grouse exhibited evidence of avoidance of pinyon-juniper across all three conifer cover classes (Fig. 3A–C). However, we found clear differences in evidence among those classes. Specifically, although the CC1 model indicated evidence of avoidance from the data, the data also provided the least support for consistent strong avoidance of this cover class (median estimate = -4.18, 95% CI=-8.48 to -0.25; see Fig. 3A). Model parameter interpretation indicated that for each 1% increase in the amount of CCI, the relative avoidance increased by 4.2% (95% CI; 2.5–8.8%). Sage-grouse also used areas with proportionately less CC1 in comparison with the proportion available (Table 2). Although a large majority of sage-grouse demonstrated evidence of avoidance of CC1 (i.e., population-level trend), a portion of individuals demonstrated proportionate use to availability, while others demonstrated selection for CC1. Evidence of avoidance for CC2 and CC3 was substantial. For example, a 1% increase in the amount of CC2 and CC3 decreased the probability of selection by 34.9% (95% CI; 18.3–56.2%) and 34.0% (95% CI; 11.0–63.4%; see Fig. 3B and C), respectively. Similarly, sage-grouse used areas with proportionately less CC2 and CC3 in comparison with the respective amounts available (see Table 2). Furthermore, we found substantially more variation among individuals in their avoidance of CC1 (see Fig. 3A) than that of CC2 and CC3 (see Fig. 3B and C, respectively).

Figure 3.

A-F, Summary of Stage 1 (left panels) and Stage 2 (right panels) analyses by pinyon-juniper cover class (CC) for greater sage-grouse populations within the Bi-State Distinct Population Segment during 2003–2005 and 2011–2013. Stage 1 inferences the log odds response (avoidance) of sage-grouse to increasing pinyon-juniper cover across three classes (CC1, CC2, and CC3). Samples were derived from posterior distribution of individuals (thin lines; red depicts avoidance trend and blue depicts selection; log-odds of 0 equates to proportionate use to availability). Median avoidance trends (thick black lines) represent the population-level pattern. Stage 2 describes the effect of relative probability of pinyon-juniper avoidance on change in survival from a baseline hazard function. Samples were randomly chosen from posterior distributions of individuals (thin red lines) to illustrate variation among individuals. Median (thick black lines) survival probabilities are depicted across continuous cover class values. The dashed black line indicates the threshold of effect of avoidance on survival (zero slope).

Using the predictive response curve (median values) of the relative probability of selection across different proportions of the cover classes, we identified thresholds in the proportion of each cover class where sage-grouse exhibited avoidance. For example, at the population level, sage-grouse were tolerant of CC1 until the value reached between 30% and 40% (see Fig. 3A); hence, sage-grouse would be expected to exhibit avoidance once that threshold was surpassed. The 30–40% CC1 threshold value calculates to ∼1.5–2.0% of actual pinyon-juniper cover that occurs within the CC1 class. Thus, at the population level, evidence of avoidance (negative log odds, see Fig. 3A) was reached within this threshold range. Sage-grouse are likely to use areas with cover below 1.5% because they represent values with positive odds (better than 1:1 odds, see Fig. 3A). For CC2 and CC3, we found a precipitous increase in avoidance as the proportion increased from zero (see Fig. 3B and C) with no clear threshold.

Table 2

Mean proportion of area estimates (95% confidence intervals) of continuous variables grouped by used and available locations, as well as locations where sage-grouse were alive and dead within the Bi-State Distinct Population Segment during 2003–2005 and 2011–2013.

We recorded 87 mortalities from the 162 telemetered sage-grouse. In the frailty stage of the linked model consisting of CC1, we found 95% CIs (1.002–1.051) for the 10-day hazard function of avoidance of pinyon-juniper did not overlap 1.0, indicating that sage-grouse that chose to avoid CC1 were more likely to survive (see Fig. 3D). Nearly all posterior samples exhibited positive relationships between relative avoidance of CC1 and an increase in annual survival. Given the estimates from the posterior distribution, a 10% increase in relative avoidance of CC1 (based on odds)was associated with an increase in annual survival of 1.91% (95% CI=1.54–2.02%; see Fig. 3E). The proportion of CC1 at locations where sage-grouse were detected alive (mean=0.36; 95% CI= 0.32–0.40)was less than the proportion here they were detected dead (mean=0.47; 95% CI=0.40–0.54; see Table 2). These descriptive summaries of raw data corroborate model inferences and provide evidence of a threshold near 30–40% of CC1 (∼ 1.5–2% actual cover). A model that consisted of CC2 and CC3 showed similar trends, but 95% CIs for changes to annual survival probability included zero (CC2: 95% CI = -1.01 to 0.81%; CC3: 95% CI=-0.14 to 2.21%). Furthermore, 95% CIs for proportions of CC2 and CC3 at locations where sage-grouse were alive versus dead overlapped (see Table 2). Variation among individuals in avoidance of areas with isolated and scattered trees suggests that sage-grouse that avoided such areas were at less risk of dying (see Fig. 3D). Although overall avoidance was most evident for areas with denser tree canopy (CC2 and CC3), benefits of avoiding these areas in terms of consistent increase in survival were not supported strongly by the data (see Fig. 3E-F).

Post Hoc Survival Analysis

In the post hoc survival analysis, we found that the time-varying effect of CC1 of pinyon-juniper on sage-grouse survival was dependent on the primary plant productivity of areas sage-grouse used. The inclusion of this productivity term (see Equation 5) to represent this more complex effect was supported by SSVS (Table S1 in the online version at  http://dx.doi.org/10.1016/j.rama.2016.09.001) and the 95% credible interval for the parameter describing this term did not include zero, indicating substantial support from the data. Namely, areas that were more productive were associated with increased hazard within CC1 (Fig. 4). The 95% CIs for the 10-day hazard function for CC1 within more productive areas did not overlap 1.0 (95% CI = 1.1–2.2). On the basis of estimates from the posterior distribution, a 10% increase in CC1 within relatively more productive areas (i.e., equivalent to a 0.5%, on average, increase in actual cover conifer) was associated with a decrease in annual survival by 4.9% (95% CI = 2.9–5.2%; see Fig. 4). Thus, our model predicted that a 1% reduction of actual pinyon-juniper cover within areas of CC1 and higher primary productivity equated to an approximate 10% change in annual survival. However, in less productive areas, ourmodel did not supportmeasurable effects on survival as CC1 increases (see Fig. 4). On the basis of the intersection of predictive response curves (median values) for both levels of productivity, we identified a threshold of 38% cover of CC1 where the benefits to survival from inhabiting high productive areas were nullified by increases in CC1 cover. For example, at relatively low levels of CC1, survival was greater in more productive areas than those that were less productive. However, survival was greater in the areas with lower rather than higher productivity once the threshold was surpassed. Thus, areas with > 38% cover of CC1 (corresponding to ∼2% actual cover) contribute to lower than average survival probabilities, as a result of increased hazard within productive areas usually associated with higher elevations. In contrast, a post-hoc survival analysis that evaluated the time-varying effect of CCP was not supported by the data, given that 95% CIs for the hazard function for CCP and by habitat productivity included 1.0 (see Table 1).

Figure 4.

The relationship between pinyon-juniper cover class 1, primary plant productivity (based on soil moisture and temperature data; Maestas et al., 2016), and survival of greater sage-grouse populations within the Bi-State Distinct Population Segment during 2003–2005 and 2011–2013. Model inferenced greater sage-grouse annual survival probabilities (median values as thick lines and samples from posterior distribution as thin lines) across increasing area of cover class 1 within areas of relatively higher primary productivity (blue lines) compared with those of less productivity (red lines).

Discussion

This study identified clear and predictable relationships between pinyon-juniper cover and sage-grouse demography. The effects were especially apparent where pinyon-juniper cover was relatively low. Sage-grouse avoided environments that consisted of pinyon-juniper trees but exhibited distinct differences in the strength of avoidance based on the approximate phase of encroachment (Phase I, II, and III; Miller et al., 2005, Falkowski and Evans, 2012). Specifically, the estimated variation in avoidance behavior among individuals was the greatest for the lowest canopy cover (i.e., CC1), which indicated selection of these sparsely encroached areas by some sage-grouse, whereas others avoided them. However, areas with greater pinyon-juniper cover (i.e., CC2 and CC3) were consistently and strongly avoided by all individuals, indicating a population-level effect with much less variation among individuals. For CC1, trees were more likely to be scattered across the landscape with < 10% canopy cover and approximated Phase I encroachment, whereas for CC2 and CC3, trees were relatively clustered, representing > 10% canopy cover and approximated Phase II and III, respectively. Interestingly, the relative importance of cover class was opposite for survival compared with avoidance, whereby disproportionate use of CC1 imposed greater risks to survival compared with CC2 and CC3. This finding was facilitated by our approach that linked decisions by sage-grouse directly to a component of their fitness and demonstrated that individual sage-grouse that avoided CC1 were more likely to survive. In contrast, the effects on survival by avoiding CC2 and CC3 were not clear because confidence limits for the change in annual survival probability included zero. Upon further investigation, we identified differences in the effect of CC1 on hazard rates in relation to our proxy for plant productivity using time-dependent models. Survival-enhancing resources afforded normally by areas with high plant productivity were nullified (or reversed) when encroached by CC1 pinyon-juniper areas. Hence, high levels of CC1 encroachment in productive habitats can lead to substantial reductions in sage-grouse survival.

Heterogeneity in pinyon-juniper cover also creates a landscape mosaic, where areas with varying levels of tree cover and overall productivity likely have distinctly different influences on local distributional patterns and population viability of sage-grouse. For example, because CC2 and CC3 (clustered trees) were largely avoided by sage-grouse, these areas are not used widely and may confine space use by sage-grouse to sagebrush communities with few or no trees. Although areas of fewer scattered trees (i.e., CC1) were generally avoided at the population level, many individuals selected these areas or used them proportional to their availability, despite the relatively greater risk of mortality. In addition, we discovered that this adverse effect pertains to sage-grouse inhabiting areas with increased overall plant production that occur typically at higher elevations (Chambers et al., 2014a) and should otherwise contribute to increased survival in the absence of trees. Encroachment of CC1 into productive areas comprises large portions of the Bi-State DPS and is occurring within close proximity of breeding leks (Fig. 5), so this pattern may strongly influence population-level demographic rates.

Figure 5.

Map depicting encroachment of pinyon-juniper CC1 in areas of relatively high primary productivity based on soil moisture and temperature data (Maestas et al., 2016) and juxtaposition to breeding leks within the Bi-State Distinct Population Segment.

We demonstrated that in the absence of pinyon-juniper, sage-grouse using high-productivity habitat had greater annual survival compared with those using low-productivity habitat (y-intercepts, see Fig. 4). However, the pattern of decreasing annual survival with increasing proportions of CCI encountered in high, but not low, productivity habitats supports the notion that sparsely distributed trees in productive habitats (particularly at CC1 > 30%, actual cover > 1.5%) may facilitate development of an ecological trap for sage-grouse. Ecological traps can occur when a habitat conveys attractive traits that act to decouple behavioral cues that an individual uses to assess habitat quality and ultimately negatively impacts a component of their fitness. Unlike population sinks, traps involve individual preference for particular resources and can be highly ephemeral with environmental change (Battin, 2004; Robertson and Hutto, 2006). Some studies of sage-grouse indicate that habitat-mediated trade-offs in life history- specific vital rates can occur, such as decreases in annual survival that are offset by increased fecundity in productive habitats (Blomberg et al., 2013b; Caudill et al., 2014). Thus, productive areas despite the presence of isolated trees might ultimately contribute to population growth. However, evidence of this scenario in the context of conifer encroachment is lacking and a recent study demonstrated that tree removal can improve reproduction for sage-grouse (Sanford et al., 2017-this issue). Furthermore, a recent range-wide meta-analysis indicated that changes in female annual survival had the most influence on annual population growth rate (Taylor et al., 2012). While not discounting unmeasured impacts on fecundity in this study, it follows that decreases in annual survival facilitated by encroaching pinyon-juniper in habitats that sage-grouse will select can have nontrivial impacts on a component of their fitness and, perhaps, influence population growth rates.

Our telemetry data and analyses provided a demographically-based empirical explanation for previously reported large-scale adverse effects of pinyon-juniper on sage-grouse population dynamics using lek data (Baruch-Mordo et al., 2013). For example, the predictive curves reported in Baruch-Mordo et al. (2013) demonstrate that little to no lek activity was largely related to areas with scattered trees with relatively low canopy cover and are consistent with our model predictions of avoidance and survival. Specifically, our findings indicate an avoidance response of pinyon-juniper by sage-grouse at cover values that correspond to CC1 = 30–40%, which correspond to ∼1.5–2.0% actual tree canopy cover. In Baruch-Mordo et al. (2013), b50.0% probability of lek activity corresponded to ∼2.0% actual tree cover. In addition, median survival predictions in our study were identified at ∼40.0% (2.0% actual cover) of CC1, which also aligned closely with the point where survival in relatively high productive areas decreased below the average value as tree cover increased. Furthermore, a similar telemetry-based study in Oregon found that sage-grouse habitat selection was negatively associated with > 3.0% conifer cover (Severson et al., 2017-this issue). From a management perspective, the target value of 3.0–4.0% actual canopy cover suggested by Baruch-Mordo et al. (2013)was consistent with patterns of high avoidance and increased risk of mortality in our study. However, our results indicate that pinyon-juniper treatments that allow > 2.0% tree cover to remain standing may still lead to unintended reductions in sage-grouse distribution, survival, and perhaps population persistence. For example, a non-trivial proportion of sage-grouse in our study demonstrated use of CC1where cover exceeded 80% (corresponding to ∼4.0% actual tree canopy cover). Based on modeled parameter estimates, sage-grouse that use areas with 80% CC1 were predicted to have a median reduction of ∼ 10% in annual survival. Furthermore, grouse that use these lightly encroached habitats in areas of relatively high plant productivity incur a median reduction in annual survival by ∼ 30% in comparison to those that use the same high productivity areas with no trees (Fig. 4) (10% less than the median). Thus, given the pattern of consistent results from these studies and ours, land and wildlife managers might consider focusing efforts in areas of CC1 and targeting removal in areas relevant to sage-grouse as lowas 1.5–2.0% tree canopy cover, which approximates 30–40% cover of CC1, especially in areas of high plant productivity.

The consensus among these studies in estimated effects was reached despite differing methodologies that included measurements of explanatory variables, types of species response, and overarching statistical approaches. For example, Baruch-Mordo et al. (2013) calculated canopy cover using estimated crown areas from a spatial wavelet analysis (Falkowski et al., 2006) and modeled lek activity as a function of pinyon-juniper and other habitat covariates using a Random Forest nonlinear modeling approach (Evans et al., 2011). In contrast, our study estimated canopy cover using automated extraction algorithms within Feature Analyst that utilize spatial context and spectral information, and used a Bayesian approach to estimate avoidance of pinyon-juniper from use versus available data and survival from time dependent frailty models. The latter two estimates also required individual bird-level data with encounter histories. Despite these differences, consistent projections among studies indicate that estimated effects are not spurious artifacts of a specific modeling technique. Instead, our study complements Baruch-Mordo et al. (2013) by providing the explanation that lek inactivity was likely a deterministic result of low sage-grouse survival and describes functional differences between phases of pinyon-juniper encroachment. Our results point to the importance of reducing the scale of observation to identify processes driving broader patterns (see Prochazka et al., 2017) and, when coupled with results from Baruch-Mordo et al. (2013), indicate that pinyon-juniper encroachment leads to population-level extirpation by adversely affecting demographic rates. Furthermore, a similar study using extensive GPS data across the Great Basin demonstrated shifts in movement speed of sage-grouse when encountering pinyon-juniper that lead to increased hazard rates, particularly for juveniles (Prochazka et al., 2017-this issue).

These collective results align with findings in sage-grouse ecology that view predation as the leading source of sage-grouse mortality (Hagen, 2011; Blomberg et al. 2013), and the differences in mortality risk among cover classes and plant productivity are likely attributable to greater exposure to predators. For example, multiple species of large Buteo hawks, which are effective predators of full-grown sage-grouse (Schroeder et al., 1999), often prefer relatively flat, open shrublands that interface pinyon-juniper woodlands, especially outlier trees from the main woodlots (Bechard and Schmutz, 1995; Coates et al., 2014b). Ferruginous (Buteo regalis) and Swainson's (B. swainsoni) hawks that nest in sagebrush ecosystems prefer nesting in single trees or isolated groups of trees surrounded by open areas with sagebrush cover (Coates et al., 2014b) and tend to build nests on large flattopped junipers. Another effective Buteo predator of sage-grouse, the red-tailed hawk (B. jamaicensis), often uses areas with less tree canopy cover to accommodate their relatively large size and wingspan (Leyhe and Ritchison, 2004). Faster and riskier movements by sage-grouse in sagebrush areas with pinyon-juniper, as indicated in Prochazka et al. (In press), likely increase sage-grouse susceptibility to visually acute predators and may be a behavioral mechanism underlying part of the pattern observed in our study. Additionally, large- and medium-sized mammalian predators can be a strong source of predation (Blomberg et al., 2013a), and may also vary in abundance with tree canopy cover and primary productivity as a result of indirect links with numbers of their principle food source, small mammals. For example, within the range of precipitation observed in the Bi-State, total species richness of small mammals should increase with productivity (Reed et al., 2006). The abundance of small mammals is greater in areas where trees were thinned compared with those without treatment (Severson, 1984;Willis and Miller, 1999), largely as a result of an increase in herbaceous understory (Miller et al., 2000; Miller et al., 2005). Our study focused on sage-grouse survival, but future studies that address nest and brood success in relation to trees would be beneficial. For example, common ravens (Corvus corax) (Coates and Delehanty, 2010) also tend to avoid foraging and nesting in juniper woodlands (Coates et al., 2014a; Howe et al., 2014) but are highly likely to nest in single trees within shrub-dominated environments (Dunk et al., 1997; Coates et al., 2014a; Howe et al., 2014).

Implications

We recognize that these findings represent new trade-offs and questions (e.g., Tausch et al., 2009) that managers should consider when evaluating treatment options for pinyon-juniper to conserve sage-grouse populations. The response curves presented in this study provide managers with support for such decisions across different canopy cover classes and levels of productivity based on R&R classes. Most importantly, the curves demonstrate that relatively high-productivity areas that consist of > 30% CC1 (equivalent to > 1.5% actual trees isolated and scattered) offer resources attractive to sage-grouse but are risky environments with respect to predation. This should be of conservation concern within the Bi-State given that ∼29% of CC1 may be found in high-productivity habitats (i.e., R&R class 1) (Fig. 5). In addition, nearly half (48%) of CC1 within identified sage-grouse habitat in the Bi-State (BLM and USFS, 2015) also comprise high-productivity habitats. Importantly, these areas likely comprise upland wet meadow areas that are vital habitat during late summer periods, especially for brooding females (Casazza et al., 2011), and might prove to be important for restoration Success.

Additionally, prioritizing removal of all trees in CC1 over thinning of those in CC2, particularly in productive areas, would likely benefit sage-grouse population vital rates most because the overall hazard for the population was predicted to decrease with less pinyon-juniper encounters. From a perspective of restoring sagebrush ecosystems, management practices aimed at thinning pinyon-juniper trees in CC2 and CC3 (without removing all trees) are commonly conducted with intent of reducing woody fuel loads and wildfire risk (Bates et al., 2005; Roundy et al., 2014a, 2014b), increasing soil moisture retention (Roundy et al., 2014b), improving understory conditions (Bates et al., 2005; Miller et al., 2014; Roundy et al., 2014a), and improving health and vigor of individual trees (Page, 2008). Because CC2 areas are also codominated by sagebrush, complete removal of trees in these areas would likely make more sagebrush habitat available to sage-grouse, assuming space formerly occupied by trees is not invaded by annual grass (Davies et al., 2011). However, from a perspective of sage-grouse population ecology, our models indicate that thinning trees (e.g., allowing isolated and scattered trees to remain standing) in areas of CC2 with relatively high primary production (i.e., R&R class 1) could have adverse consequences on specific sage-grouse population vital rates by creating riskier conditions, effectively similar to CC1 (e.g., ecological traps). For example, we demonstrated that as low as 2.0% isolated and scattered tree cover (40% CC1) can lead to avoidance and reduced survival of sage-grouse.

Importantly, our results do not inform how sage-grouse might respond to additional variation imposed by the specific type of pinyon-juniper treatment employed. Sagebrush ecosystem processes respond variably to mechanical, chemical, and controlled burning treatments of pinyon-juniper (Tausch et al., 2009; Davies et al., 2011; Strand et al., 2013; Pyke et al., 2014), and sage-grouse would be expected to respond variably as well. However, certain environmental conditions may make some treatment options more tenable, thus benefiting both sage-grouse and sagebrush ecosystems. In particular, high-elevation mountain sagebrush habitats have become increasingly invaded by pinyon-juniper, in part due to suppression of wildfire (Brooks et al., 2015). Burning of sagebrush and associated pinyon-juniper in low-elevation areas with warm and dry conditions has been clearly demonstrated to increase risk of annual grass invasion and promote an annual grassfire cycle that consumes and destroys more sagebrush. However, wet and cool high-elevation areas have inherently greater resilience to disturbance (such as fire) and resistance to annual grass invasion (Chambers et al., 2014a, 2014b). In the short term, burning of sagebrush across elevation gradients associates with immediate reductions in population growth by decoupling positive benefits of precipitation (Blomberg et al., 2012; Coates et al., In press). Over longer time periods, however, burning high-elevation CC2 encroached areas that are highly avoided by sage-grouse but costly to remove completely by mechanical means may not have the same immediate adverse effects but instead help increase population growth. Management strategies are system dependent, and the expected effects of such strategies on sage-grouse populations and restoration success will likely vary in response to the underlying environmental constraints.

To the extent that the effects of pinyon-juniper in the Bi-State region are somewhat representative of those experienced by other sage-grouse populations, our results can be generalized to other portions of sage-grouse range. A simple comparison of the proportion of CC1 and CC2 within sage-grouse habitats between Bi-State DPS to proportions of CC1 and CC2 mapped across the rest of NV and northeastern CA reported by Coates et al. (2015) indicated that the Bi-State consisted of 15.7% and 1.6%, less CC1 and CC2, respectively. This comparison, Coupled with our findings of adverse effects of CC1 on survival, indicates the threat of pinyon-juniper encroachment to sage-grouse habitat is likely similar or more pronounced in other areas of the Great Basin. Lastly, our analytical approach for linking selection to survival can be applied within existing conservation planning tool (CPT) frameworks on the basis of resource selection functions (Farzan et al., 2015) to make spatially explicit predictions for where pinyon-juniper removal might best benefit sage-grouse, which would now take into account the relative mortality risk (hence, fitness consequences) caused by different phases and densities of pinyon-juniper.