Fine-scale spatiotemporal studies can better identify relationships between individual survival and habitat fragmentation so thatmechanistic interpretations can be made at the population level. Recent advances in Global Positioning System(GPS) technology and statistical models capable of deconstructing high-frequency location data have facilitated interpretation of animal movement within a behaviorally mechanistic framework. Habitat fragmentation due to singleleaf pinyon (Pinus monophylla; hereafter pinyon) and Utah juniper (Juniperus osteosperma; hereafter juniper) encroachment into sagebrush (Artemisia spp.) communities is a commonly implicated perturbation that can adversely influence greater sage-grouse (Centrocercus urophasianus; hereafter sage-grouse) demographic rates. Using an extensive GPS data set (233 birds and 282,954 locations) across 12 study sites within the Great Basin, we conducted a behavioral change point analysis and constructed Brownian bridgemovementmodels fromeach behaviorally homogenous section.Wefound the probability of encountering pinyon-juniper among adults was two and three times greater than that of yearlings and juveniles, respectively. However, the movement rate in response to the probability of encountering pinyon-juniper trees was 1.5 times greater for juveniles. Parameter estimates indicated a 6.1% increase in the probability of encountering pinyonjuniper coupled with a 6.2 km/hour increase in movement speed resulted in a 56%, 42% and 16% increase in risk of daily mortality, for juveniles, yearlings, and adults, respectively. The effect of pinyon-juniper encounters on survival was dependent on movement rate and differed among age class. Under fast speed movements (i.e., flight), mortality risk increased as encountering pinyon-juniper increased across all age classes. In contrast, slower speeds (i.e., average) yielded similar adverse effects for juveniles and yearlings but not for adults. This analytical framework supports a behavioral mechanism that explains reduced survival related to pinyon-juniper within sagebrush environments,whereby encountering pinyon-juniper stimulates riskier movements that likely increase vulnerability to visually acute predators.
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.
Sage-grouse were captured and monitored at 12 study sites across two states (northern Nevada, n = 10; northeastern California, n = 2) that encompassed three sage-grouse management zones (MZ III- Southern Great Basin, MZ IV-Snake River Plain, MZ V-Northern Great Basin) with different climatic, elevation, topographic, soil, and geologic characteristics (Miller and Eddleman, 2001) (Fig. 1A). Sagebrush communities comprise a mixture of Wyoming big sagebrush (Artemisia tridentata wyomingensis), basin big sagebrush (A. t. tridentata), mountain big sagebrush (A. t. vaseyana), and low sagebrush (A. arbuscula). Sites within the drier and more southern portion of the study area (i.e., southern Great Basin) also possess a strong black sagebrush (A. nova) component (Miller and Eddleman, 2001). Other prevalent shrub species that occur across the entire study area include rubber rabbitbrush (Ericameria nauseosa), green rabbitbrush (Chrysothamnus viscidiflorus), antelope bitterbrush (Purshia tridentate), snowberry (Symphoricarpos spp.), andwestern serviceberry (Amelanchier alnifolia). Pinyon-juniper woodlands are more prevalent at higher elevations across the Northern and Southern Great Basin MZs and nearly absent within the portion of the Snake River Plain MZ that encompassed our study area. Both species are typically restricted to mountain slopes and demonstrate dominance patterns above (pinyon) and below (juniper) 2 100 m (Tausch et al., 1981; Miller et al., 2000). Tree density is greatest among the mountain ranges of east-central and northeastern Nevada with west-facing slopes containing slightly higher densities (Tausch et al., 1981).
We captured sage-grouse at each field site using spotlighting techniques (Wakkinen et al., 1992) at night near known active leks during spring and fallmonths from1 April 2012 to 27May 2015 (Table 1). Captured individuals were fitted with a rump-mounted GPS-platform transmitter terminal (PTT) unit (Northstar Science and Technology LLC, King George, VA). GPS-PTT unit mass did not exceed 3% of sagegrouse body mass. The GPS-PTT unit power source consisted of a solar arraywith rechargeable battery (operational life=2–3 yr).Duty cycles were programmed to record between 10 and 13 locations per day (06:00, 07:00, 08:00, 09:00, 10:00, 11:00, 12:00, 13:00, 14:00, 15:00, 16:00, 19:00, 21:00 PST) with the intent of capturing one nocturnal (i.e., inactive period) location and multiple diurnal (i.e., active period) locations. GPS-PTT units recorded location coordinates, date and time stamps, and location accuracy estimates. Location data were downloaded using Argos PTT decoding software (PTT Tracker) and manually screened for erroneous locations. Screening was conducted by a single individual in order to remove the potential for variation in screening behavior. Sage-grousewith fewer than 100 locationswere removed from the dataset in order to meet statistical assumptions in the BCPA, specifically only incorporating movement paths long enough to capture the exponential decay in the autocorrelation (Gurarie et al., 2009). Values reported are means ± SE.
Conifer cover wasmapped at a 1-mresolution using object recognition algorithms in Feature Analyst (Overwatch Systems, Sterling VA) software from 2013 National Agriculture Imagery Program (NAIP) imagery across the entire study area (US Geological Survey, unpublished data) (Fig. 1B).We chose to use this relatively high-resolution map because previously available datasets derived from Landsat imagery (US Geological Survey) were too coarse (30-m resolution) for our analysis. We could not distinguish all potential species of conifers (e.g., firs, subalpine pines) from pinyon and juniper using the Feature Analyst algorithm. However, pinyon and juniper dominate conifer species composition throughout most of the Great Basin (Tausch et al., 1981; Miller et al., 2005), so we assumed that conifers mapped and used in our analyses are primarily pinyon-juniper.
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).
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).
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: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: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:
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.
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
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).
Behavioral Change Point Analysis
In the initial step of the BCPA we transformed raw location data into velocity (i.e., movement rate) and turning angle estimates. Population level means for velocity and turning angle estimates were 1 795.65 m/day ± 48.34 and 0.0002° ± 0.0049, respectively. The second data transformation, which converted track statistics into the modeled parameters of PV, had mean values of -1,268.56 ± 181.79 (µ), 7301.27 ± 692.99 (σ), and 0.02250 ± 0.00007 (ρ).
Significant variation in modeled parameters existed among age and sex classes. Adult individuals exhibited the highest overall value for the µ parameter (-242.58 ± 30.15) and lowest overall value for the σ parameter (3 190.24 ± 130.73), which is indicative of a movement process that is faster and more direct. Yearling individuals demonstrated the opposite in PV parameters and exhibited the lowest overall value for the µ parameter (-5 010.69 ± 945.69) and highest overall value for the σ parameter (22 408.60 ± 3 630.21). As such, the movement process of yearlings, on average, was more variable in both speed and direction in comparison with the other two age classes. Juvenile birds were intermediate for both the µ and σ parameters, which suggest an averagemovement process that is more variable in speed and direction when compared with adults and less variable when compared with yearlings. The µ (-2,004.82 ± 357.16) and σ (10,201.06 ± 1,371.90) parameters for femaleswere nearly twice as lowand high, respectively, as those observed for males (-859.32 ± 115.36) (6 147.01 ± 464.83). Therefore, females are demonstratingmovement processes that are, on average, more variable in both speed and direction as compared with males. Considerable variation also existed in PV parameters across hours. Interestingly, there was significant within-hour correlation across the different PV parameters. Peaks and concurrent troughs were documented for the µ and σ parameters during the morning and evening hours, respectively. Between 07:00 and 15:00 hours there was a gradual decline in the µ parameter and concurrent increase in the σ parameter. In other words, the movement processes across all individuals decreased in speed and directionality from early morning to late afternoon. From late afternoon to early evening the directionality and speed of the movement process increased to values more similar to early morning hours. Thesemirrored movement processes following and preceding crepuscular hours and most likely represents movements to and from roost locations.
Brownian Bridge Movement Models
A total of 25 049 BBMMs were constructed from 153 236 locations and across 214 birds. The mean length of time (minutes) and number of locations per BBMMwere 305.87±0.99 and 6.12±0.02, respectively. Approximately 10.4% of the BBMMs constructed contained a probability value of encountering pinyon-juniper ≥ 1%. However, we found considerable variation in the probability of encountering pinyonjuniper among sex and age classes. Females demonstrated greater probabilities (0.0087 ± 0.0001) than males (0.0052 ± 0.0001), and an increase in the probability of encountering pinyon-juniper correlated with age. Adults demonstrated probability values (0.0091 ± 0.0001) that were nearly twice that of yearlings (0.0056 ± 0.0002) and three times greater than juveniles (0.0035 ± 0.0002). Probability values for encountering pinyon-juniper are low, as would be expected for an organism that is considered a sagebrush obligate and the inverse relationship between percent sagebrush and pinyon-juniper cover (Tausch et al., 1981; Miller et al., 2000).
Modeling Movement Rate as a Function of Conifers
All models demonstrated over whelming support over the nullmodel, as evidenced by the difference in DIC values (ΔDIC) (Table 3). Model comparisons using DIC revealed that an additive effect of age class and the probability of encountering pinyon-juniper explained the greatest amount of variance for these data. All predictions were based on the effects derived from this best-fit model. Input variables (i.e., probability of encountering pinyon-juniper) used for model predictions were truncated to the highest value observed from the BBMM analysis in order to restrict predictions to a realistic range. The posterior distributions of the intercept parameters, which signify the baseline movement rates (i.e., absence of pinyon-juniper) in meters per day, for juveniles, yearlings, and adults were 1 109.8 (95% CI; 949.7–1 286.9), 891.8 (95% CI; 776.4–1 015.8), and 862.6 (95% CI; 757.0–966.9), respectively.
Model predictions showed an increase in movement rate across all age classes in response to an increased probability of encountering pinyon-juniper, and this relationship appeared exponential (Fig. 4). However, responses to increased probabilities of encountering pinyonjuniper were not linear across age classes. For every 10% increase in the probability of encountering pinyon-juniper, there were concurrent increases of 892.8 m (95% CI; 653.0–1 158.2), 640.9 m (95% CI; 477.7–808.6), and 611.1 m (95% CI; 459.8–759.3) in the daily movement rates of juveniles, yearlings, and adults, respectively. This equates to a dailymovement rate that is 1.5 times greater in juveniles than yearlings and adults. Under a hypothetical scenario in which juvenile, yearling, and adult birds move from an area of no pinyon-juniper to the greatest observed probability value of encountering pinyon-juniper, the daily displacement estimated for each age class would be 6 428.1 m (95% CI; 4 071.5–8 338.8), 4 614.4 m (95% CI; 3 439.4–5 821.8), and 4 400.1 m (95% CI; 3 310.8–5 466.6), respectively.
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.
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
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.
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.
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.
We thank K. Andrle for coordinating and managing field work and numerous field technicians and volunteers for capturing and monitoring hundreds of sage-grouse.We thank B. Halstead, J. Yee for statistical support, C. Overton for expertise with spatial ecology, and W. Perry, E. Sanchez-Chopitea, K. Mauch, and T. Kroger for mapping of pinyonjuniper and GIS expertise. Manuscript drafts were reviewed by K. Klasing and 2 anonymous referees.
 ☆ Data collection for this project was supported logistically and financially by multiple federal (Bureau of LandManagement,US Fish andWildlife Service); state (Nevada Department of Wildlife, Nevada Sagebrush Ecosystem Program, California Department of Wildlife); and private (Ormat Technologies, LS Power, Midway Gold) cooperators. Funding for data analysis was provided by the USGS Western Ecological Research Center and the Sage-Grouse Initiative.