Introduction

In the Great Basin of western North America, contemporary encroachment of pinyon (primarily Pinus monophylla) and juniper (primarily Juniperus osteosperma) (hereafter, “pinyon-juniper”) woodland is contributing to the degradation of sagebrush ecosystems and habitat quality for sagebrush obligate species. Greater sage-grouse (Centrocercus urophasianus; hereafter, sage-grouse) are one such obligate whose populations have declined concomitantly with the loss and fragmentation of sagebrush ecosystems and now occupy slightly more than half their former range (Schroeder et al., 2004; Miller et al., 2011). This trend has contributed to multiple proposed listings of the species under the Endangered Species Act (CFR, 2015). Pinyon-juniper encroachment, which we define as both expansion and infill, contributes to sagebrush degradation by disproportionately reducing shrub and herbaceous habitat components that provide forage and cover to sage-grouse (Miller et al., 2005, 2011). Subsequent increases in woody biomass can also contribute to larger wildfires that kill sagebrush and propagate invasion of annual grasses such as cheatgrass (Bromus tectorum) that fuel more wildfires (Getz and Baker, 2008; Balch et al., 2013; Chambers et al., 2014). It follows that many restoration actions focus on thinning or removal of encroaching pinyon-juniper in sagebrush ecosystems, with the intent of enhancing sage-grouse populations by improving habitat quality (Tausch et al., 2009; Davies et al., 2011; Bi-State-Action-Plan, 2012; McIver and Brunson, 2014; Pyke et al., 2015).

Identifying how animals move in relation to habitat features and the demographic consequences of thosemovements can increase the effectiveness of conservation measures, but it is difficult to quantify these patterns given the inverse relationship between scale of observation and degree of stochasticity (Levin, 1992; Morales et al., 2010). As such, standard measures of animal-habitat interactions (e.g., habitat avoidance or selection) are quantified typically at the population level (Doherty et al., 2008; Freese, 2009; Casazza et al., 2011; Baruch- Mordo et al., 2013; Knick et al., 2013). However, higher-order patterns such as population persistence or distribution stem from individualbased behavioral decisions (e.g., move or stay) made in response to internal and external stimuli (Gurarie et al., 2009), and management strategies that rely on information gathered at larger scales must take into consideration how information is being transferred across scales (Levin, 1992). In the case of sage-grouse, Baruch-Mordo et al. (2013) documented a cessation of lek activity in response to tree cover that exceeded 4%, yet it was unknown whether lek cessation was a product of mortality, emigration, or both. If birds are emigrating to neighboring leks and stabilizing or improving population performance as a result, the negative effects of pinyon-juniper encroachment are likely to be ascribed to reductions in sage-grouse distribution alone. However, if pinyon-juniper encroachment also leads to a reduction in individual survival, the negative effects experienced by the metapopulation could be greater. In a similar analysis, Coates et al. (2017 this-issue) used a temporally higher resolution dataset, which incorporated mortality information (i.e., individual-level data) and found that tree canopy cover of 4% in productive habitat with dispersed trees was associated with a substantial reduction in the annual probability of survival. In other words, by decreasing the scale of observation, Coates et al. (2017 thisissue) identified a likely demographic mechanism for the result observed in the Baruch-Mordo et al. (2013) study. However, the behavioral mechanisms, such as speed and frequency of travel, which might drive lower survival in pinyon-juniper encroached habitats, were not identified. Decreasing the scale of observation to an even finer resolution may help identify these behavioral mechanisms. For example, avoidance patterns observed at coarser scales may reflect faster movements through less desirable habitats. General theories of optimality predict that organisms should spend a greater portion of their time in habitats where resources are plentiful but leave habitats and the resources they provide more quickly when risk of predation is high (Morris and Davidson, 2000). Movement is a temporal function, and therefore it logically extends that less desirable habitats should correlate with faster movements. However, in the case of sage-grouse, increases in movement rate may also increase susceptibility to visually acute predators and could serve as the mechanism driving the demographic consequences documented in Coates et al. (2017 this-issue). Identifying groups of individuals or life stages that demonstrate greater movement rates and, in turn, experience greater risk during encounters with pinyon-juniper encroached habitats could provide managers with a better informed approach to tree-removal strategies. For example, if adult females were to experience the greatest risk of mortality during trips between lek and nest sites, it would be beneficial for managers to focus tree-removal efforts on areas that intersect lek and nest site locations. Such efforts could reduce predator subsidies associated with pinyon-juniper encroachment and reduce movement rates in sagegrouse by decreasing the distance between desirable resources.

Identifying variation in sage-grouse movement behavior resulting from encounters with pinyon-juniper trees in sagebrush landscapes and their direct linkages to fitness has also been hampered by technological and analytical limitations. Previous studies have utilized relatively coarse spatiotemporal datasets such as traditional VHF telemetry data to quantify avoidance patterns of pinyon-juniper by sage-grouse (Doherty et al., 2008; Freese, 2009; Casazza et al., 2011, Coates et al. 2017 this-issue) or lek counts (Baruch-Mordo et al., 2013) and presence/absence of leks (Knick et al., 2013) to quantify effects of pinyon-juniper on integrated demographic processes. From the technological standpoint, recent advances in Global Positioning Systems (GPSs) small enough to place on small-game-size animals have increased the temporal accuracy of location data by reducing the length of time between sampling intervals and removing “the human-monitoring element” by transmitting data via satellite or cellular network. From the analytical standpoint, these large volumes of GPS-derived data led to the adoption of more advanced space-use estimators capable of dealing with serial autocorrelation and complex data error structures (Horne et al., 2007; Gurarie et al., 2009; Fischer et al., 2013). Yet these techniques often need to be adapted and modified in novel ways to address specific management questions and improve understanding of complex systems (Ovaskainen, 2004; Horne et al., 2008; Hebblewhite and Haydon, 2010; Morales et al., 2010).

Movement is a continuous process observed at discrete intervals. Therefore, relating movement processes to land-cover characteristics requires probabilistic space-use estimators. Brownian bridge movement models (BBMMs) have garnered considerable attention for their ability to estimate utilization distributions (UDs) along autocorrelated movement paths. However, before calculating BBMMs it is imperative that the data used to inform the model originate from a uniform movement process or that the UD estimated along a bridge reflect the heterogeneity across relocations (Kranstauber et al., 2012). This becomes problematic when working with temporally long datasets, as it is well known thatmobile organisms are capable of displaying a variety ofmotility patterns and thus introduce heterogeneity to relocation datasets. Constructing a BBMM from data depicting ≥ 2 movement processes without the use of a dynamic BBMM (Kranstauber et al., 2012) would result in an “averaged” value of the Brownian motion variance (BMV), ameasure of animalmobility (Horne et al., 2007), and in turn an underestimation of the UD during periods of rapid movement and an overestimation during periods of slower movement. If the end goal is to extract habitat data on the basis of an encounter probability, the previously described scenarios would result in the false exclusion and inclusion, respectively, of habitat features.

To better link changes in behavior to environmental stimuli, and subsequent demographic consequences of those changes in behavior,movement paths with underlying heterogeneous behavior should first be segregated statistically into homogenous sections that allow construction of more accurate and behaviorally specific BBMMs. State-space modeling is a popular approach to this issue, which allows for separation of a process model (e.g., movement behaviors) from an observation model (i.e., observation error) but also requires the number and type of behavioral states to be known a priori (Jonsen et al., 2005; Royer et al., 2005; Patterson et al., 2008). Gurarie et al. (2009) provided an alternative likelihood-based method for identifying changes in animal movement called behavioral change point analysis (BCPA) that does not require a priori knowledge of the number and types of behavioral states that an organism can transition into or out of, and is robust enough to deal with observation error,multidimensionality, and autocorrelation. Using the BCPA to identify shifts in the underlying movement process and constructing UDs from BBMMs that are informed by behaviorally homogenous location data provide a meaningful space usemetricwithin amovement context. Hierarchical models that predict movement patterns in relation to encounters with external stimuli (e.g., habitat characteristics extracted to the UDs within a geographic information system) can then be developed. Once established, consequences to fitness can be modeled directly by relating when an animal dies to how its behavior changed when encountering different levels of habitat-related external stimuli (Gurarie et al., 2009; Halstead et al., 2012).

Kranstauber et al. (2012) developed a streamlined approach to computing dynamic Brownian bridges by combining a BCPA and BBMMinto a single R package (R Core Team 2014). However, that framework does not allow the user to adjust tuning parameters within the BCPA, extract movement statistics fromthe location data, or provide biological knowledge of the organismas it pertains to daily movement patterns.While it is possible to use the Brownian motion variance as a metric for habitat use (Byrne et al., 2014), we felt that themovement ratemetric calculated from the stand-alone BCPA was a more intuitive measure of an organism's interaction with its environment. As such, our approach to deconstructing homogenous movement paths mirrored Kranstauber et al. (2012), albeit a manual version, whereas the movement statistic chosen for posteriori analyses followed Buchin et al. (2015).

Given the economic costs of pinyon-juniper removal across thousands of acres of encroached sagebrush habitat (SageSTEP, 2011), successful implementation of conservation strategies would be enhanced by knowing the behavioralmechanisms underlying space-use dynamics and their demographic consequences to sage-grouse when traversing landscapes with variable levels of pinyon-juniper encroachment. For example, identifying individuals or groupswithin a population that utilizes resources within or adjacent to pinyon-juniper stands and quantifying their movement rates and survival risk associated with those encounters could provide managers with a finer-resolution treeremoval strategy and reduce the need for more costly, broad-scale applications. Hence, the objectives of our study were twofold. First, we estimated the movement response of sage-grouse across age classes to probabilistic pinyon-juniper encounters using movement statistics as a behavioral metric. Second, we estimated the survival risks associated with the behavioral response that probabilistic pinyon-juniper encounters elicit.We predicted that 1) an increase in the probability of encountering pinyon-juniper would result in faster movement rates by sagegrouse, 2) males and younger individuals would demonstrate the fastest movement rates in response to pinyon-juniper encounters based on an untethered (i.e., lacking brood and nest responsibilities) motility pattern in the former and a lack of familiaritywith the environment in the latter, and 3) risk of mortality incurred by each age class would vary as a function of the unique behavioral strategies observed in the movement rate analysis.

Methods

Statistical Analyses

Behavioral Change Point Analysis

An initial transformation of raw location data was required at the outset of the BCPA and involved the conversion of absolute positions and their associated time stamps into velocity (i.e., movement rate) and turning angle estimates. A second transformation was applied in order to produce a persistence velocity (PV) estimate, which Gurarie et al. (2009) describe as the tendency and magnitude of a movement to persist in a given direction in a given time interval. For every location in space-time the PV, which serves as a measure of the underlying movement process, can be described by three parameters: 1) a mean (µ), 2) a variance (σ), and 3) a continuous autocorrelation (ρ). Changes in the µ, σ, and ρ parameters correspond to different interpretations of the underlying movement process (Gurarie et al., 2009). With respect to PV, an increase in µ will correspond to a movement process that is both faster and more directed, whereas increases in σ will correspond to a movement process that is more variable in both speed and direction. Increases in ρ are indicative of more directedmovement, irrespective of speed.

In the subsequent step of the analysis a moving “window” was swept over all locations, per individual, shifting one location forward in time. Before shifting, eight models were constructed from the locations inside the window, and a most likely break point (MLBP) was selected on the basis of model fit, whichwas informed by the PV statistics calculated for each location during the previous step. The eight possible models corresponded to a significant change in 1) no parameters (i.e., null); 2) µ; 3) σ; 4) ρ; 5) µ and σ; 6) µ and ρ; 7) σ and ρ; or 8) all three parameters. Bayesian information criterion (BIC) was used in model comparison, and the lowest BIC value corresponded to the best model fit for that window, given a value of 2 for the sensitivity parameter (K). The sensitivity parameter is a constant by which the loglikelihood is multiplied and in that way serves as a penalty against model complexity-lower values chosen for the sensitivity parameter result in the selection of a simpler model. The sensitivity parameter, window size, and number of times an individual location was selected as the MLBP (threshold value) are referred to as “tuning knobs” and do impose a certain level of subjectivity to the analysis. However, using a simulated dataset, Gurarie et al. (2009) found that accurate estimates of significant change points could be produced using a window size of 50, a threshold value of 7, and sensitivity parameter equaling 2. For that reasonwe adopted the same values for our analysis.MLBPs chosen ≥ 7 times during the windows weep were stored for later use.Model assumptions (i.e., Gaussian error structure, exponential decay in autocorrelation) were verified by visual inspection of normal Q–Q plots and discrete autocorrelation functions. All BCPA analyseswere conducted in R version 3.1.1 (R Core Team, 2014) using the package bcpa (Gurarie, 2014).

Figure 1.

Map of A, study sites and sage-grousemanagement zones (MZs) and B, pinyon-juniper distribution across northernNevada and northeastern California. Inlay is provided to represent scale and extent of the pinyon-juniper raster layer.Datawere collected fromGlobal Positioning System-marked greater sage-grouse at 12 study siteswithin the Great Basin during 2012–2015.

Table 1

The distribution of Global Positioning System-platform transmitter terminal units by site, state, sage-grouse management zone (SGMZ), age class, and sex. Unitswere deployed on greater sage-grouse at 12 study sites within the Great Basin during 2012–2015

Brownian Bridge Movement Models

Results of the BCPA were used to divide the behaviorally heterogeneous location history for each individual into behaviorally homogenous segments (Fig. 2), thus fulfilling the assumption of behavioral homogeneity within a single BBMM. Locations chosen as breakpoints were assigned to location histories preceding, but not following, them. This resulted in temporal but not necessarily spatial gaps between bridges (see Fig. 2) and circumvents the constant aggregation of behaviorally dissimilar locations,which occurswhen you incorporate a breakpoint into the location histories preceding and following it. In instanceswhere the BCPA failed to insert a behavioral change point between two separate diurnal periods, we did so manually. For example, if a behaviorally “homogenous” segment included locations between 12:00 hours 01 January and 12:00 hours 02 January, we split that segment at 00:00 hours 02 January. We also removed all nocturnal locations from the data set because sage-grouse are active during diurnal periods and inactive at night. Finally, we removed all locations greater than 1 hr apart from a preceding or subsequent location to lessen the chance of violating assumptions of conditional random movement, which increases with increased time elapsed between locations (Horne et al., 2007; Fischer et al., 2013). This final step resulted in the removal of approximately 1.08% of the data and likely had little impact on the final results of our analysis.

The BBMM requires an estimated error for all location data and a grid cell size for the output of the UD. We chose a uniform location error of 20 m based on GPS-PTT manufacturer specifications and a grid cell size of 1min order tomatch the resolution of our pinyon-juniper raster layer. The resulting raster layer represents a BBMM UD with an unrestricted confidence interval (Fig. 3A). The pinyon-juniper binary layer (values=0 or 1)was cropped tomatch the extent of the BBMMUD raster layer (values = 0–1) (Fig. 3B). We then multiplied the two raster layers together (Fig. 3C), resulting in probability values greater than zero only if the underlying pinyon-juniper raster contained a value equaling one (i.e., presence of pinyon-juniper). This new raster layer represented the probability of encountering pinyon-juniper for each individual BBMM. We cropped the probability of encountering pinyonjuniper raster layer to an “extent raster” created in a separate part of the analysis, which represented the 95th percentile of the original BBMMraster layer (Fig. 3D).We summed the probability of encountering pinyon-juniper raster values, which resulted in a single probability value of encountering pinyon-juniper for a given, behaviorally homogenous segment. All BBMMs were calculated in R version 3.1.1 (R Core Team 2014) using the package BBMM (Nielson et al., 2013). Raster layers for each BBMM were constructed using the raster package (Hijmans, 2015).

Figure 2.

Example of Brownian bridge movement models (BBMMs) constructed from two sets of behaviorally homogenous location data identified during the behavioral change point analysis in a pinyon-juniper-encroached sagebrush landscape. The “star” in the inset map indicates location of the BBBMs.

Figure 3.

Example visual representations. A, Brownian bridge movement model (BBMM). B, Pinyon-juniper raster layer clipped to the extent of the BBMM. C, Product of the BBMM and pinyon-juniper raster layer. D, Probability of encountering pinyon-juniper within a given BBMM, restricted to the 95th percentile of the original BBMM.

Modeling Movement Rate as a Function of Conifers

Fromthe first transformation in the BCPA analysis,we estimated velocity (i.e., movement rate), which provides an easy interpretation of one of the most fundamental behaviors in animal movement (Turchin, 1998). We then fit linear mixed-effects models (Kéry, 2010) using the single probability value of encountering pinyon-juniper for a given behaviorally homogenous segment as a predictor variable for movement rate, which took the form:

denoting a change in movement rate (y) of β for a unit change in pinyon-juniper (PJ; continuous, percent cover) and age class (AG; categorical; juvenile, yearling, and adult).We additionally fit randomeffects to account for repeated measures (i.e., BBMMs) within individual (bird[ς_i]), spatial (site [κ_j]), and temporal (year [η_k]) correlation (Gillies et al. 2006). Subscripts i, j, and k reference bird, site, and year, respectively. A second model included estimation of a sex parameter instead of age. Two additional models were developed, one that assumed that the effects of pinyon-juniper (slopes) differed among individuals shared for each age class and then a second model developed to be shared across sex classes. Essentially, these models extended the mixed linear model to include random coefficients without correlation between intercepts and slopes (Kéry, 2010), similar to concepts described in Gillies et al. (2006).

Each BBMM was constructed from a sequence of location data, and as such each BBMM had multiple associated movement rate estimates. In order to balance the data sets (i.e., movement rate, probability of encountering pinyon-juniper), we averaged the movement rate values across each BBMM. Movement rate estimates were normalized by applying a Box-Cox transformation (λ = 0.2). Models were constructed using Bayesian analysis with specified vague priors for all parameters in each model. Parameter definitions were specified in Table 2. Models were run on three chains of 10 000 iterations each following a burn-in period of 10 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). Models were compared against each other using deviance information criterion (DIC). DIC is a goodness-of-fit statistic used to rank competing models within a Bayesian framework and is similar to the maximumlikelihood- based Akaike information criterion (AIC) (Burnham and Anderson, 2002). Results are presented as DIC and ΔDIC, where ΔDIC represents the difference between the top model and the model with which it is being compared.

Shared Frailty Model

Unlike typical survival models, shared frailty models utilize random effects to account for intraclass correlation and allow for frailties (i.e., mortality risk) to be estimated across groups for a particular set of covariates (e.g., movement rate, probability of encountering pinyon-juniper) (Gutierrez, 2002).We utilized the same randomeffects structure as in themovement rate analysis and estimated frailties across age classes in the presence and absence of pinyon-juniper. We modeled survival as a continuous process observed at discrete intervals (1 day) and assumed an equal probability of survival across all intervals (i.e., constant hazard model) (Halstead et al., 2012). The model was expressed as:

where the unit hazard (UH) represented a constant hazard function, , with intercepts that varied across age classes (k), estimated coefficients formonth (which represents in accordancewith month as a categorical variable) movement rate ), probability of encountering PJ (MR), and an interaction between MR and P (. To investigate the effect of MR and P, we specified a shared frailty across the three levels of age class, referenced as k. This specification allowed us to estimate differences in the effects on survival between juveniles, yearlings, and adults. Because individuals can graduate into a higher age class with time, we fit age as a time-varying covariate within the model. On 1 March (approximate start to lek activity) of each year, individuals within the juvenile and yearling age classeswere graduated into the yearling and adult age classes, respectively. It was not possible to graduate beyond the adult age class. Month was included because sage-grouse have been shown to have unequal survival probabilities across their annual lifecycle (Blomberg et al. 2013a,b Moynahan et al. 2006). Subscripts h, i, and j, referenced bird, site, and year, respectively. Models were specified with random effects for year () and site () to account for temporal and spatial intraclass correlation, respectively. Covariate values were converted to daily averages in order to match the frailty interval length. Because we modeled movement rate and the probability of encountering pinyon-juniper as time-varying covariates, both variables could only be assessed on days in which the individual was observed. Therefore, we specified a first-order Markov process to impute missing values for days on which individuals were not observed as described in Halstead et al. (2012). A normal distribution with a standard deviation based on observed values was used to impute missing values. Movement rate and probability of encountering pinyon-juniper estimates were normalized by applying a Box-Cox transformation (λ = 0.2). Using the UH, a cumulative hazard (CH)was calculated thatwas expressed as:

where 365, 1-day intervals derived an annual survival parameter (S) as:

However, the duration of the study spanned 1 152 days (1 April 2012 to 27 May 2015). Because individuals entered into the study (i.e., captured and marked) at different times, we chose a study start date based on the date on which the first individual was captured (1 April 2012) and transformed all other start dates to days since the beginning of the study.

The model was run on three chains of 5 000 iterations each following a burn-in period of 10 000 iterations and thinned by a factor of 5. Convergence (i.e., a stationary posterior distribution) 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). Because the chance of stochastic variation increases with analytical scale resolution (Levin, 1992), posterior distributions for odds-ratios were evaluated at the 85% credible interval.

Table 2

Description of posterior distribution formedian estimated parameter (2.5 and 97.5 percentiles in parentheses) inmodeling the effects ofmovement and survival from Global Positioning System-marked greater sage-grouse at 12 study sites within the Great Basin during 2012–2015

Results

We captured and marked 233 sage-grouse and obtained 282 954 GPS telemetry locations across all study sites (see Table 1). The sex ratiowas ∼70% (n=162) female to ∼30%(n=71) male, and the structure of age-classes at time of capture was 148 adults (∼64%), 43 yearlings (∼18%), and 42 juveniles (∼18%). We were not able to deploy as many GPS-PTT units on juveniles and yearlings as on adults because of body mass limitations. The disparity among sexes was not based on logistic constraints or study objectives but is representative of the naturally occurring ratio that exists between male and female sage-grouse (Connelly et al., 2011).

The mean number of locations per bird was 1 214.4 ± 71.7. The mean number of days that an individual was recorded (≥ 1 location per day) during the study was 204.9 ± 11.0. Sage-grouse capturing and marking were conducted during spring and fall months to ensure that locationswere being collected throughout the annual cycle. Certain factors (e.g., feather coverage of solar array, shortened photoperiod, sage-grouse behavior) reduced the rate of battery recharge or inhibited GPS-PTT units from acquiring a fix and resulted in longer gaps between locations and a lower mean daily location rate (6.36 ± 0.15) as compared with the preprogrammed duty cycle (10–13 locations per day).

Shared Frailty Model

The daily risk of mortality across months showed considerable variation, with the lowest estimated hazard ratios at 0.23 (95% CI; 0.14–0.86) for October and 0.33 (95% CI; 0.14–1.33) for December (see Table 2). Both months also demonstrated the least amount of variation in their hazard ratio estimates. Conversely, the greatest estimated hazard ratios occurred duringmonths that are strongly associated with the reproductive life phases of sage-grouse (Coates et al., 2013). In particular, the daily risk ofmortality in Augustwas 5.59 (95% symmetric CI; 2.75–7.30) times greater than the baseline hazard function. April, May, August, and October were the only months for which the posterior distributions (95% CI) did not overlap 1.

Table 3

Model selection for movement rate. Data were collected from Global Positioning System-marked greater sage-grouse at 12 study sites within the Great Basin during 2012–2015

Figure 4.

Effect of the probability of encountering pinyon-juniper on sage-grouse movement rates, as a function of age class. Six thousand samples from the posterior distributions of model parameters were used to represent the range of uncertainty in these relationships. Colored lines represent the uncertainty in the estimates for juveniles (green), yearlings (orange), and adults (blue), respectively. Black lines represent the median effects. Data were collected from Global Positioning System-marked greater sagegrouse at 12 study sites within the Great Basin during 2012–2015.

For every one unit increase (6.04%, 6.21 km/hr) in the modeled parameters (probability of encountering pinyon-juniper, velocity), the daily risk of mortality was 1.56 (85% CI; 1.15–1.79), 1.42 (85% CI; 1.05–1.71), and 1.16 (85% CI; 1.03–1.31) times greater than the baseline hazard function for juveniles (Fig. 5A), yearlings (Fig. 5B), and adults (Fig. 5C), respectively. Posterior distributions (85% CI) for juvenile, yearling, and adult birds did not overlap 1. Predictions were based on the median value of the posterior distribution of the baseline hazard function for each age class. The average rate of travel for juveniles, yearlings, and adults was 74.8, 57.0, and 64.7 m/hr, respectively. The average probability of encountering pinyon-juniper across age classes was 0.34%, 0.56%, and 0.90% for juveniles, yearlings and adults, respectively. We found support for an interaction between probability of encountering pinyon-juniper and movement rate on survival, meaning the effect of pinyon-juniper encounters on survivalwas dependent on movement rate. However, this effect differed among age classes (Fig. 5A, B, and C). For example, under relatively fast movement rates (i.e., those that involve flight presumably), mortality risk increased as encountering pinyon-juniper increased across all age classes. In contrast, slower movement rates (i.e., average hourly movement rates) yielded similar adverse effects of pinyon-juniper on survival of juveniles and yearlings but no effects were supported by the data for adults.

Discussion

Results from our modeling approach provide useful insight into the behavioral mechanisms underlying sage-grouse movements and their subsequent consequences for fitness. In terms of movement, our findings of sage-grouse increasing rates of travel with increased probabilities of encountering pinyon-juniper encroached habitats corroborate previous studies that have identified pinyon-juniper as an aversive plant community (Commons et al., 1999; Atamian et al., 2010; Casazza et al., 2011; Baruch-Mordo et al., 2013; Knick et al., 2013). Our results extend further to indicate that the response of individuals encountering pinyon-juniper encroached habitats varies as a function of age class. For example, adults experienced probabilities of encountering pinyonjuniper that were nearly twice that of yearlings and three times greater than juveniles. The lower probabilities calculated for the juvenile age class may in part be explained by their movement rates in response to pinyon-juniper encroachment, which were nearly 1.5 times greater than yearlings and adults. In essence, as age increases, the response (i.e., movement rate) to pinyon-juniper encroachment decreases while the encounter rates increase. This should ultimately produce greater amounts of time spent in higher-density stands of pinyonjuniper for adults as compared with the other two age classes.

Figure 5.

Annual survival probabilities (%) for A, juveniles, B, yearlings, and C, adults, as a function of the interaction between movement rate and the probability encountering pinyon-juniper (%). Green and red areas of the plane represent scenarios inwhich survival probabilities are better andworse than the average survival probability outside of pinyonjuniper communities, respectively. Black areas of the plane indicate an annual survival probability of 0. Data were collected from Global Positioning System-marked greater sage-grouse at 12 study sites within the Great Basin during 2012–2015.

We also identified significant differences in movement parameters among age classes from the BCPA. Specifically, adult birds produced more ballistic (i.e., fast and determined) movements compared with the other two age classes, irrespective of plant community type, which indicates retention of a cognitive map of their surroundings (Kamil and Jones, 1997) used to make predetermined movement decisions (Stankowich and Blumstein, 2005), and is corroborated by interannual nest site fidelity (Schroeder and Robb, 2003). Conversely, the movement processes among yearlings were highly variable in both direction and speed, which is characteristic of frequent switching among behavioral states observed among many organisms (Bell, 1990; Morales and Ellner, 2002; Fauchald and Tveraa, 2003). Organisms encountering novel environments generally demonstrate combinations of movements as they alternate between searching for resources and taking advantage of existing resources (Fauchald and Tveraa, 2003). Accordingly, BCPA patterns displayed by yearling sage-grouse suggestmore unfamiliar situations in comparison with adult birds. In contrast, juvenile birds appeared to adopt an intermediate strategy during a brief period (September–March) when they congregate communally (Beck, 1997); hence decisions are influenced more by group behavior.

On the basis of survival analysis, our findings indicate that juveniles experience the greatest risk when encountering pinyon-juniper encroached habitats, regardless of movement speeds. Low survival probabilities among the youngest individuals within a population is a common trend across a wide range of taxa, and inherently lower fitness characteristics may explain why individuals belonging to the juvenile age class are more susceptible to the adverse conditions associated with pinyonjuniper encroached plant communities. As part of a learning experience, juveniles may be testing pinyon-juniper encroached habitats as potential resource opportunities and the faster movements that they display may reflect the paucity of resources that those plant communities provide. In turn, their behavior (i.e., greater movement rates in response to pinyonjuniper encroached communities) may be soliciting their whereabouts to visually acute predators that favor those areas. Adult and yearling birds, which have obtained prior experience with their surroundings, should have a cognitive map of these habitats, and the slower movement rates that they display in the presence of pinyon-juniper could reflect resource utilizations in small highly productive pockets, which happen to be adjacent to pinyon-juniper woodlands (Coates et al., 2016). Adult and yearling age classes are also adversely affected by pinyon-juniper encounters, evidenced by posterior distributions for their hazard ratios that were greater than one. However, the slowermovement rates exhibited by yearling and adult age classes,may be insulating themfromnegative predator encounters by maintaining a sufficient level of crypticism in an otherwise risky habitat.

The strong links we identified between rapid movement and reductions in survivalmay also indicate changes to underlying predator-prey functional responses, whereby decisions made by sage-grouse to travel faster through pinyon-juniper encroached habitat likely increased their susceptibility to predators.When sampling at discrete intervals, individuals that move along a ballistic trajectory (perhaps engaged in flight) will appear to move faster over a given period of time as compared with those that adopt a more diffusive pattern (i.e., slow and stealthy) (Visser and Kiørboe, 2006), and the more ballistic the trajectory, the greater the likelihood that a predatory encounter will occur (Lima and Dill, 1990; Visser and Kiørboe, 2006). All age classes of sage-grouse are at risk of avian predation from large raptor species (Schroeder et al., 1999) that have strong visual acuity (Jones et al., 2007) and will utilize pinyon-juniper encroached sagebrush environments as nesting and perching substrate (Dwight andMurphy, 1973; Coates et al., 2014). Predation pressure from raptors may be especially high during fall raptor migration (Smith et al., 2008; Blomberg et al., 2013a) and exacerbated further by perch subsidies provided by trees. Sage-grouse of all age- classes are also vulnerable to mammalian predators that use both visual and olfactory-based methods of prey detection (Hartzler, 1974; Blomberg et al., 2013a). Hence, sage-grouse in pinyon-juniper-encroached plant communities are likely exposed to higher predation risk frommultiple species of predators. Visual verification of mortalities did not occur within an acceptable temporal window (Blomberg et al., 2013a) during this study to definitively differentiate between avian and mammalian predation/scavenging events. Therefore,we cannot provide an explicit estimate for either.

Pinyon-juniper woodlands are predicted to continue expansion in sagebrush communities throughout the Intermountain West (Tausch et al., 1981; Miller et al., 2005; Romme et al., 2009). Given the results of this study, we should expect changes in sage-grouse movement behavior and subsequent consequences in survival. The large geographic distribution of sage-grouse (Schroeder et al., 2004) precludes the use of spatiotemporally coarse tree-removal strategies, which become logistically and financially infeasible at larger scales. Results from this analysis may help provide land and wildlife managers with cost-effective strategies for conservation planning. In particular, we discovered that sage-grouse are at a substantially greater risk of mortality during encounters with pinyonjuniper encroached habitats, especially younger birds based on their movement behavior. Thus, focusing tree-removal efforts in areas utilized by sage-grouse during brood-rearing and juvenile life phases will likely provide the greatest benefits to sustaining sage-grouse populations. However, high hazard ratios were observed for March, April, and May. This temporal period is marked by extremely overt movement behaviors, wherein females make large movements between courtship arenas (i.e., leks) and nest locations across multiple days (Schroeder et al., 1999). Where pinyon-juniper woodlands intersect or parallel lek and nest locations, decreased probabilities of survival should be expected for adult and yearling individuals. Given the reproductive contributions and, therefore, “worth” of adult and yearling female sage-grouse to population growth (Johnson and Braun, 1999), tree-removal scenarios that highlight areas utilized by those individuals should not only increase individual survival probabilities but are also likely to have a disproportionately greater contribution to population growth as well.

In summary, the objectives of this analysis were to describe the underlying process ofmovement in an avian species as a function of probabilistic encounters with an aversive plant community and to subsequently differentiate behavioral risk in the presence and absence of that community. The results indicate that 1) sage-grouse are moving faster through environments with the presence of pinyon-juniper, 2) rate of movement in response to pinyon-juniper encroached habitats decreases with age, 3) the probability of encountering pinyon-juniper habitats increases with age, and 4) increases in movement rate during pinyon-juniper-habitat encounters resulted in significantly higher risks to survival for all age classes. Variation in survival probabilities among age classes may be explained by differences in movement behavior as a function of encounterswith pinyon-juniper encroached habitats. Finally, our results linking movement with survival in relation to sage-grouse encounters with trees of different densities dovetail well with a companion article in this volume (Coates et al., 2017 this-issue) that established linkages between habitat selection and survival in relation to sage-grouse encounterswith trees at different densities and configurations (e.g., phases of encroachment). Collectively, these novel studies provide keen insight into how and why sage-grouse move through and select habitat within a pinyon-juniper-encroached landscape, as well as the resulting fitness consequences of those decisions.