INTRODUCTION

Fundamental to population ecology is a thorough understanding of population structure and dynamics at varying spatiotemporal scales (Allen and Hoekstra 1992, Levin 1992). Trends of population growth and decline, as well as underlying demographic processes that guide them, often vary across hierarchical levels of organization and their corresponding scales (Allen and Star 1982). A major goal in wildlife conservation is to align processes and patterns to their appropriate spatiotemporal scale and further identify biotic and abiotic factors governing such patterns. Macro-ecologists have recently stressed the importance of documenting and understanding cross-scale interactions, whereby environmental factors at regional spatiotemporal scales that strongly influence an ecological response can operate differently when interacting with factors occurring at local (and nested) scales (Peters et al. 2007, Soranno et al 2014). It follows that management actions at one ecological scale (e.g., representing a local subpopulation or population) might influence outcomes differently at other scales (e.g., regional; Hobbs 1998). However, incorporating spatiotemporal scales, particularly those occurring at regional extents (Soranno et al. 2014) to the management of wildlife populations can be challenging. Failure to account for scale-specific processes can result in significant misinterpretation of observed patterns and subsequent actions misaligned with drivers of population change (Bissonette 2016).

Integrated population models (IPMs) provide an empirically driven framework to investigate population dynamics by incorporating multiple sources of data within a single unified framework (Schaub and Abadi 2011, Kéry and Schaub 2012). For example, IPMs can include information on abundance through time and space using survey information (e.g., time-series count data), as well as information on population processes derived from more in-depth demographic information (e.g., life stage data). An IPM framework that combines data from multiple sources provides better parameter estimation and an understanding of processes and patterns, even in circumstances when data are missing or disparate (Schaub and Abadi 2011). IPMs also provide a platform to evaluate density-dependent effects on each population vital rate (Weegman et al. 2016) and recover parameters about specific demographic rates where data might be missing or sparse (Abadi et al. 2010, Kéry and Schaub 2012). For a detailed review of other benefits of IPMs see Schaub and Abadi (2011) and Kéry and Schaub (2012).

We focus on 2 properties of IPMs that lend themselves well to applied management. First, the inherent hierarchical properties of IPMs allow for estimation of vital rates and abundance across nested spatial and temporal scales, which would otherwise be inestimable and cause gaps in time-series data (Abadi et al. 2010, Kéry and Shaub 2012, Duarte et al. 2017). This property is particularly important because intrinsic and extrinsic factors influencing population dynamics at more localized extents representing subpopulations may differ from those at more regional extents representing larger metapopulations (Saether 1997). Demographic estimates spanning a sufficient temporal duration are also important for populations that exhibit cyclical dynamics that respond to resource availability (Bjornstad and Grenfell 2001, Ahrestani et al. 2016). Hence, contrasting IPM-derived estimates of vital rates, and rates of population change, among local and regional extents through time can disentangle factors likely affecting subpopulations at smaller spatial extents (e.g., anthropogenic disturbances) from those affecting metapopulations at larger extents (e.g., region-wide wildfire and drought).

This type of comparative framework would be extremely valuable for conservation and management of species' populations, allowing wildlife and land managers to make decisions that align with patterns and processes at the appropriate spatiotemporal scale (Bissonette 2016). For example, long-term population trends for single species can be evaluated at scales that encompass the species' distributional range to help inform federal regulatory agencies tasked with status assessments and protection policies (e.g., U.S. Fish and Wildlife Service) such as the Endangered Species Act (ESA; 16 U.S.C. § 1531 et seq.). Whereas, understanding processes driving population dynamics at relatively small spatial scales can provide critical information for conservationists and managers determining actions for local populations within a species' range.

Second, IPMs can be fit with different hierarchical, random effects (e.g., population and year) that allow “borrowing of strength from the ensemble” (Kéry and Schaub 2012) such that data can be shared across space and time to improve vital rate parameter estimates within multiple, subcomponent demographic models. The structuring of random effects for space and time (i.e. additive or nested) influences more finite estimation of grand means for parameters (Kéry and Schaub 2012, Nowak 2015). Different random effect structures, however, create tradeoffs between precision and accuracy of parameter estimates, which affect joint likelihood estimates of abundance. Among local subpopulations of birds (i.e. with spatial random effects included), data that inform specific demographic vital rate estimates often are limited. In these instances, information can be shared among subpopulations, which increases precision and reduces uncertainty. On the other hand, effects of sharing could be relaxed to increase subpopulation-specific estimates that are independent when sample sizes are robust among all subpopulations. Similarly, data sharing can occur across years, which allows for derivation of parameters annually and across longer time periods. Thus, the demographic component of IPMs can be improved by fitting the most parsimonious random effects that balance precision and accuracy.

Here, we describe an IPM for a Distinct Population Segment (DPS) of Greater Sage-Grouse (Centrocercus urophasianus, hereafter referred to as “sage-grouse”) found along the border of California and Nevada, known as the Bi-State DPS. Sage-grouse are sagebrush obligates (Braun et al. 1976) that have declined in distribution and abundance largely as a result of human disruption of sagebrush ecosystems (Schroeder et al. 2004, Knick and Connelly 2011). Sage-grouse have undergone multiple evaluations for listing under the ESA and, specifically, the Bi-State DPS has been evaluated separate from other sage-grouse populations throughout the species' range (U.S. Fish and Wildlife Service 2013, 2015a, 2015b). The Bi-State DPS may be at increased risk of population declines because of substantial loss of sagebrush habitat and lack of connectivity among subpopulations comprising the DPS (Oyler-McCance et al. 2005, 2014). The primary threats within the Bi-State are habitat loss as a result of conifer encroachment (Coates et al. 2017, Prochazka et al. 2017), wildfire (Coates et al. 2016), and annual grass invasion (U.S. Fish and Wildlife Service 2015a).

Sage-grouse populations in the Great Basin inhabit sagebrush ecosystems where productivity is tied strongly to variation in precipitation and temperature owing largely to the cold-desert climate (Noy-Meir 1973). In the Bi-State DPS, population dynamics at larger spatial scales likely follow cyclical patterns that alternate between short bursts of primary productivity resulting from periods of above-average precipitation, and then often followed by periods of drought. Increased primary productivity likely increases resources needed for reproduction and survival, whereas drought imposes the opposite effects and drives populations downward (Fedy and Doherty 2011, Blomberg et al. 2013, Coates et al. 2016). It is imperative, therefore, to identify when and where perturbations decouple trends at smaller spatial scales from trends at larger spatial scales governed by environmental stochasticity. Different types of IPMs have been used to describe sage-grouse population dynamics in the more mesic northeastern portion of the species' range (McCaffery and Lukacs 2016) and for Gunnison Sage-Grouse (C. minimus) in Colorado (Davis et al. 2014). Other than preliminary reports (Coates et al. 2014), however, population dynamics of sage-grouse in the Bi-State DPS within the more xeric Great Basin have not been described with IPM methods used in this study.

The primary objective of this study was to develop a stochastic stage-based IPM to provide insight into the demographic processes responsible for the state of subpopulations and any decoupling from broader population trends. Specifically, we (1) evaluated different random effects (e.g., subpopulation and year) and density-dependence for each life stage to refine the demographic process models; (2) integrated demographic and lek count data using a joint likelihood process; (3) estimated demographic rates across each life stage for different age and subpopulation classes; and (4) derived λ at the subpopulation and metapopulation extent through time. We then extended the model to fit the climatic characteristics of precipitation as covariates on λ and compared these effects among spatial extents (subpopulation and metapopulation). We also compared annual estimates of λ derived from the individual state-space and demographic sub-models to those derived from the joint likelihood IPM to examine relative differences in variance and alignment of trends among the 3 model outputs.

METHODS

Data Analysis

Step 1: Evaluate random effect structures of vital rate models. Before developing the process model (i.e. demographic matrix model) of the IPM, we conducted an initial step that evaluated separate models composed of different random effect and density-dependent structures for each population vital rate. One purpose of this step was to evaluate the appropriate degree of data sharing (i.e. borrowing of strength) across space and time to estimate parameters for each vital rate. We collected demographic data to estimate parameters for clutch size, nest survival, second nest propensity, hatchability, chick survival, yearling survival, and adult survival. We relied on published data to help inform first nest propensity rates and juvenile survival (explained below). We modeled each vital rate separately and created 5 models with different random effect structures representing different spatiotemporal effects: null (no random effects), year only, subpopulation only, year and subpopulation additive, and year and subpopulation nested. A model with random effect (year) only allowed parameters to vary across time, but not space, and a model with random effect (subpopulation) only allowed variation across space, but not time. Both structures allow for increased degree of data sharing (e.g., estimation pulled toward a grand mean). The additive space and time structure allow estimates to vary across subpopulation and years, which should be more robust in cases where data are sparse (e.g., missing subpopulation or year combinations) and some degree of data sharing is needed. Lastly, a nested random effects structure offers the least amount of data sharing and should be most parsimonious with clear differences among subpopulations and years with robust sample sizes. Greater degrees of data sharing can help inform disparate data sets and increase precision in estimates (i.e. narrower credible limits). On the contrary, less data sharing might increase accuracy for individual subpopulation and year combinations where sample sizes are robust. However, where data are scarce, uncertainty in parameter estimation will increase. Thus, this method allows the model fit to indicate the balance between precision and accuracy given the data for each life stage.

Because λ is likely regulated by density-dependent effects similar to other species (Turchin 1995, 1999), we also fit density-dependent terms within each vital rate model to account for variation explained by population size. Several simple density-dependent models exist (Turchin 1995, 1999). Here, we considered the Ricker model that assumed constant linear decrease in the demographic rate as population size increases N (Ricker 1954), and the Gompertz model that assumed constant linear decrease as a function of logarithmic transformation of N. Thus, the second form favors stronger density dependence with relatively smaller N. To account for time-delayed dynamics (Jacobson et al. 2004), we also fit 1-yr lag effects (Garton et al. 2011, 2015). Each life stage model consisted of 8 candidate models to identify the most parsimonious random effect structure and density-dependent model form. We calculated Watanabe–Akaike information criterion (WAIC; Watanabe 2010) to compare models fitted with the suite of random effect structures. The WAIC is a fully Bayesian technique that is valid in hierarchical models (Watanabe 2013) and is generally recommended for Bayesian ecological analyses (Hooten and Hobbs 2015). We chose the model with the lowest WAIC score to represent each vital rate model. Model specifications that follow represent the most parsimonious model structure for a given life stage (i.e. lowest WAIC). Importantly, we derived estimates based on the model with the most support from the data. For example, if WAIC favored a nested random effect, then parameters for a particular vital rate were derived for each combination of subpopulation and year. If WAIC supported a model with additive or single random effect (no nested effect in the model), then the supported model was used to derive parameters across year or subpopulation accordingly. Lastly, if WAIC supported the null model, then parameters were derived across all subpopulations and years. The data preparation and most parsimonious model expressions (Appendix Table 4) for each life stage follow.

Yearling and adult survival models. Using telemetry data, we located individual monitored sage-grouse and modeled annual survival as a continuous process observed at discrete monthly intervals. We monitored survival at daily to weekly intervals during the breeding season (e.g., March–September) using ground-based telemetry, and at monthly intervals during the winter months using aerial telemetry. Therefore, encounter histories used to calculate survival of individual yearling and adult grouse were scored as alive or dead for each month of the year, which corresponded to our lowest temporal resolution for monitoring survival. We “right-censored” individuals from the data, meaning we did not always observe the date of death and only had knowledge that individuals survived beyond a certain point. We considered censoring to be a random process, because sage-grouse were censored under 3 circumstances: (1) the study period ended before death was determined, (2) the transmitter failed, or (3) sage-grouse were lost to follow-up monitoring. Thus, all individuals either died or were eventually right-censored. We also allowed for graduation of yearling sage-grouse to adults between years. Inference was based on a constant hazard model, meaning the probability of mortality did not change across months. We recognize that sage-grouse survival can vary markedly among seasons (Blomberg et al. 2013) and more complex error structures can be fit to account for uneven survival probabilities across the annual life cycle of sage-grouse. Because our objective was to estimate the probability that grouse lived or died for a given year, rather than when they lived or died within that year, our assumption of constant monthly hazard was appropriate.

Modeling procedures for annual survival followed those described in Halstead et al. (2012). The unit hazard (UH) was modeled where each encounter interval represented a Bernoulli process per month. The additive effect model was expressed as:

where UH was a function of random effects for subpopulation α_i, year γ_j, and subpopulation by year ζ_ij (nested structure). Random variables were assumed to arise from normal distributions with mean zero, and variances , and , respectively. Expected change in the ln(hazard ratio) was estimated by age (β_age), where x_age represents the vector of age class data (yearling and adult) and the indicator was equal to one for adults. The hazard ratio represented the ratio of hazard rates (in this case, monthly risk of morality) between the 2 age classes. Subscripts h, k, i, and j reference grouse, month, subpopulation, and year. Following the modeling process, we derived the annual (an) survival parameter (s) as:

where CH represented the cumulative hazard (T=12 months represented annual survival). We also attempted to fit a first-order autoregressive hazard model; adding this complexity resulted in no convergence and stationary distribution.

First nest propensity. We did not develop a model for nest propensity (np1) for first nests based on data collected in the Bi-State DPS. Rather, an extensive study conducted by Taylor et al. (2012) provided information across numerous, range-wide sage-grouse studies and reported 0.96 (95% confidence interval [CI] = 0.94–0.97) and 0.89 (95% CI = 0.87–0.91) for adults and yearlings. These values were considered reliable because of the large number of studies used in the analysis. However, as a more conservative approach, we provided priors that were slightly wider than these published confidence intervals. We set informative priors for adults = Beta (97,5) and yearlings = Beta (90,12). In the absence of additional information, we assumed this proportion to be constant among subpopulations and years.

Clutch size. Data for clutch size (c1) of first nests were modeled using the Poisson distribution as follows:

Thus, the log expected count of clutch size μ_c1 of first clutch for grouse h and year j is a linear function of random year effects γ_j, which were assumed to arise from a normal distribution with mean zero and variances . We estimated a change in the expected count by age (β_age), where x_age represented the vector of age class data (yearling and adult) and the indicator was equal to one for adults. Parameter estimates for clutch sizes of second nests (c2) were derived in a separate modeling process that followed the same procedures.

Nest survival. Survival parameters of first (ns1) and second (ns2) nests were derived separately. Nest survival was modeled as proportional hazards as expressed in equation 1. However, for nest survival estimation, encounter histories were developed for each nest and modeled as a continuous process at discrete, daily intervals. In this case, T=38 was used to estimate survival during 10 days of laying and 28 days of incubation. Female age was used as a fixed effect in this model. This model did not consist of random effects based on lack of evidence during step 1.

Egg hatchability. Data on egg hatchability (h) in the Bi-State DPS were derived for most nests that were successful. Hatchability was modeled as arising from a Binomial distribution (logit-link function) and took the form: where the initial number of eggs laid in a clutch e represented the number of trials (n_h,eij) at subpopulation i and year j, and the number of hatched eggs represented successes (y_h,eij). The logit(p_h,eij) is a linear function of random subpopulation effects α_i, random year effects γ_j, as well as subpopulation and year effects combined ζ_ij, each of which were assumed to arise from normal distributions with mean zero, and variances , , and , respectively, and a change in the expected count of magnitude βage, where the indicator of age was equal to one for adults.

Second nest propensity. Parameters were derived for the second nest attempt propensity (np2) directly from data collected in the Bi-State DPS. Data to estimate second nest propensity were more reliable than those for the first nests because females were monitored intensively (relocated every 1–3 days) following first nest failures. Field technicians were able to determine whether or not a female attempted a second nest by approaching sage-grouse and confirming nesting when sage-grouse were located at the same coordinate as the previous location. Second nest propensity data were modeled as arising from a Binomial distribution as follows: where the number of unsuccessful nests at each subpopulation in each year were denoted by n_np2_,ij. In this model, y_np2_,ij represents the number of renests and logit(p_np2_,ij) is a linear function of random subpopulation effects α_i, random year effects γ_j, as well as random subpopulation and year effects combined ζ_ij, each of which were assumed to arise from normal distributions with mean zero, and variances , , , respectively. The influence of age and density dependence on np2 were measured as fixed effects with magnitude βage and βdd, where the indicator of age was equal to one for adults, and the density-dependent variable was the natural log of abundance with a 1-yr lag.

Chick survival. As described previously, chicks were not directly marked and followed using radio-telemetry to reduce researcher disturbance. Thus, chick survival (cs) probabilities were derived from the 2 brood counts with time-interval lengths that varied throughout the 10-yr study period. However, the number of days elapsed from nest hatch to brood count (survival periods) varied by study year (2003–2005, 50 days; 2007–2009, 35 days; 2010–2011, 28 days; 2012–2015, 50 days). Rather than choosing the minimum survival period (28 days), we accounted for variable survival periods to allow full use of data collected across all years of the study as follows. We modeled chick survival based on brood-count data as arising from a Binomial distribution where the initial brood size was scored as the number of trials and chicks that survived to days d were scored as successes and took the form: where d on the binomial probability p is d=d(j) and represents one of 3 survival periods depending on the year j of data collection of time period (t = 28, 35, or 50). For a 35-day interval, the probability of survival is modeled by this logistic relationship:

In this model, y_cs,bi represents the number of chicks that survived for each brood b at subpopulation i. The logit(p_cs,bi,35) is a linear function of random subpopulation effects α_i. The influence of age and density dependence on chick survival were measured as fixed effects with magnitude β_age and β_dd, where the indicator of age was equal to one for adults, and the density-dependent variable was the natural log of abundance with a 1-yr lag. We assumed a constant hazard function, and consistent with this assumption the probabilities of survival for the other intervals are related as follows:

Juvenile survival. Juvenile sage-grouse (js; post-fledging, >35 days and <1 yr old) were not radio-marked and tracked in the Bi-State DPS. However, we derived a posterior distribution of juvenile survival rates (js) during this period by using an informative prior of 0.75 (95% CI = 0.67–0.82) reported in Taylor et al. (2012) in the form of .

Step 2: Developing the IPM. The framework of the IPM was a state-space model, where the state process was informed by the demographic data and observation process was informed by the lek-count survey data. The state process of the IPM consisted of an age-structured demographic matrix based on annual and subpopulation parameter estimates of survival (s) and fecundity (f). Using the demographic model structures described previously, f was a derived parameter that took the form:

where i, j, and a represent subpopulation, year, and age class, respectively. Fecundity (f) was divided by 2 as this represents a female-based demographic model and we assumed equal sex ratios at hatch (Atamian and Sedinger 2010). Sub-models of population vital rates consisted of stochastic processes, where posterior distributions of population vital rate parameters were estimated.

Lek counts were appropriate for Markovian processes of time-series data where observations are imperfect (Baumgart 2011). Each lek was counted multiple times within a season, and the maximum male count was assumed to reflect peak male attendance at leks. The average of these maximum counts (y_obs) were linked with breeding male apparent abundance N at subpopulation i and year j using the equation

For sage-grouse lek counts, this Poisson distribution most appropriately represents error in observations because sage-grouse can be undercounted and overcounted during observations, and as the true population becomes larger the degree of error becomes larger, too. Undercounting (false negatives) can take place when grouse are difficult to observe and overcounting (false positives) is possible when grouse move around during the counting process. Initial N at each subpopulation was drawn from a discrete uniform distribution. We also compiled the number of leks, maximum number of males per lek, number of active leks, and percentage of active leks by subpopulation and year (Table 1).

TABLE 1.

Summary of lek survey data averaged across years for Greater Sage-Grouse in the Bi-State Distinct Population Segment study area, California and Nevada, 2003–2015. Number pairs in parentheses are lower and upper limits of the 95th percentile.

In a final step of the IPM, we used joint likelihoods, of which lek counts informed the observation process (equation 11) and demographic data informed the state process (equations 2 and 10). By fitting random effects, N was estimated through time and across subpopulations. We estimated unobserved N over the next 5 yr into the future by extending the loop of the state process and adapted the prior definitions of the demographic rates (Kéry and Shaub 2012). From the estimated posterior distribution of N we derive the finite rate of change (λ; Caswell 2001), which took the form:

where i represents subpopulation, j represents year, and j−1 represents a previous year. The instantaneous per capita rate of change (r) then was calculated by natural logarithmic transformation (Gotelli and Ellison 2006), and is expressed as:

An annual estimate of abundance was calculated for the metapopulation extent (Bi-State region) by summing the annual abundance estimates across the 6 subpopulations. Estimates of and r were calculated at the metapopulation extent using equations 12 and 13, respectively. For each subpopulation, we calculated the probability that the subpopulation was increasing, stable, and decreasing based on the posterior distributions of derived parameters. We then calculated the odds of increase from the probability values, where the odds of increase represented the probability of increase divided by the sum of the probability of decrease and stability. Likewise, we calculated the odds of decrease. We then created a ratio of the 2 odds (OR; increase:decrease) and applied natural logarithmic transformation to that ratio. The purpose of this procedure was to quantify evidence for population trends (when ratio equals 0, then the odds of increase are the same as decrease). We categorized values of ln(OR) between 0.9 and 1.1 as stable, values <0.9 as decreasing, and values >1.1 as increasing. We report and plot observed lek counts, N, λ, and ln(OR) for each subpopulation and the Bi-State DPS. Posterior distributions of parameters were summarized as median and 95% credible intervals (CRI), expressed as 0.025–0.975 quantile. Lek count observation data were summarized by mean, as well as 2.5^th and 97.5^th percentile of the data distribution. Vital rates were averaged across time and space when random effects were evidenced, and reported by age (Table 2). Detailed information of estimated parameters for each vital rate are available ( https://www.sciencebase.gov/catalog).

TABLE 2.

Summary of posterior distributions of derived population vital rate parameters for Greater Sage-Grouse in the Bi-State Distinct Population Segment study area, California and Nevada, 2003–2015.

Lastly, we derived λ using 3 sources of modeling for comparison: (1) only a demographic matrix model (DMM), (2) only a state-space model (SSM) using lek count data, and (3) IPM (as described above). The DMM was an age-structured stochastic process that derives λ using only likelihoods for S and f as described above (equations 2 and 10) and did not include count observations. The SSM model relied only on count observations and did not incorporate demographic likelihoods. In this case, the SSM separates process variance (that is, environmental) from observation error (Kéry and Schaub, 2012) by partitioning each variance component, which took the form:

Here, the state process (equations 14 and 15) was modeled while accounting for observation error (equation 16). Equation 16 maps the true state of the process onto the observed data (y_t), which in this case is an annual sum of averaged maximum counts across each of the 6 subpopulations. Similar to the observation process in the IPM, the errors in the counts were modeled using a Poisson distribution with a mean equal to the variance. We assigned vague priors to the initial (t = 1) population size (equation 17).

We used JAGS 3.4.0 (R version 3.1.1) using the package rjags to obtain posterior samples of parameters. We used Markov chain Monte Carlo methods and ran 5 independent chains of 100,000 iterations, following a burn-in of 99,000 iterations and an adaptive phase of 1,000 iterations. Chains were thinned by a factor of 500 because of storage limitations given the large number of parameters. We found no evidence for lack of convergence, observed by examining history plots. We could not compute the R-hat statistic (Gelman et al. 2004), typically used to measure convergence, because of limitations in memory owing to the large number of parameters in our models. Instead, we examined history plots which indicated no evidence for lack of convergence.

Post hoc precipitation analysis. We conducted a post hoc analysis using the IPM to investigate links between λ and climatic conditions. Because precipitation has been reported as an influential driver of recruitment (Blomberg et al. 2002) and population growth rates using lek count only data within the Great Basin (Coates et al. 2016), we chose to focus our analysis on precipitation. Spatially explicit data for local measurements of precipitation at a spatial resolution of 800 m from 2002 to 2015 were obtained from the PRISM Climate Group (http://www.prism.oregonstate.edu/). We restricted our analyses to include variables related to precipitation to limit the number of covariates in our climatic model set, and because effects due to precipitation are more readily explained from a perspective of sage-grouse life history (e.g., relations between rainfall, primary productivity, and available resources for grouse) than from extremes in temperature, particularly for a cold-adapted gallinaceous species.

We aggregated precipitation into seasonal and annual intervals based on a priori hypotheses that aligned with phenology of sage-grouse. For example, precipitation measured during March–May represented our hypothesis that resources from precipitation would benefit sage-grouse during nesting life stage (i.e. growth of perennial grasses and forbs). During June–August, precipitation would influence those resources required for successful brood-rearing (e.g., wet meadow productivity, forb growth, delay of plant senescence). Precipitation during September–November was thought to influence spurts of new growth and water requirements that influence late brood-rearing and juvenile survival during late summer for all sage-grouse. During December–February, precipitation (primarily as snow) might contribute snowpack for increased runoff in following seasons that was thought to drive more productive and possibly longer growing seasons. We formed 6 models that represented additive effects for our precipitation-based hypotheses, which were (1) snowpack runoff + nesting; (2) nesting + brood-rearing; (3) brood-rearing + extended brood-rearing/juvenile; (4) snowpack runoff + nesting + brood-rearing; (5) nesting + brood-rearing + extended brood-rearing/juvenile; and (6) all life stages (annual precipitation). Although heavy snowfall during winter can have adverse impacts on overwinter survival (Moynahan et al. 2006, Anthony and Willis 2009), generally winter survival rates are greatest among the seasonal survival rates (Blomberg et al. 2013). Therefore, we aligned winter at the beginning rather than the end of the precipitation year to represent our hypothesized positive carryover effects (e.g., snowpack melt and runoff) on successful reproduction and increased recruitment. We also investigated the influence of logarithmic transformation of the annual change in precipitation [ln(P_t/P_t−1)] with 1-yr lag on r (t + 1), extending on findings from previous analyses (Blomberg et al. 2013, Coates et al. 2016). Our final model set consisted of 20 candidate models.

The model was formulated as a 2-stage process: (1) estimating N from the IPM and deriving r for each subpopulation (sp) and year (equations 12 and 13); and (2) in a parallel stage we fit the derived parameter r sampled from the full posterior distribution as a function of different precipitation (pr) covariates, although total variance (parameter uncertainty + process uncertainty) may be slightly underestimated possibly owing to inflated sample size and requires further study. The model took the form:

where r represents samples from the posterior probability distribution of per capita rate of change for subpopulation i and year t, and β represents model coefficients. The model Posterior distributions were derived using Program JAGS within the rjags package (Plummer et al. 2015) in R version 3.1.1 (R Core Team 2014). Specifically, posterior distributions of parameter estimates were generated from 5 chains of 5,000 iterations each, after a burn-in of 5,000 using Markov-chain Monte Carlo (MCMC) methods. The chains were thinned by a factor of 5 and the number of adaptation iterations was 1,000. Convergence of MCMC output was assessed visually with history plots and the R-hat statistic, where values <1.1 indicated convergence (Gelman et al. 2004). Priors were selected to be uninformative. We compared models using WAIC (Hooten and Hobbs 2015) and chose the model with the lowest to represent seasonal precipitation effects. We provide plots displaying the relationship between λ and change in precipitation (Δprecipitation) across the Bi-State and varying estimated slopes of this relationship for each subpopulation. We also provide plots of N varying with Δprecipitation across each subpopulation to illustrate the decoupling of specific subpopulations through time.

RESULTS

Population Estimates and Trends

Survey data consisted of >2,000 independent lek counts resulting in 607 lek surveys that reflected annual maximum male attendance for each lek. We monitored 354 sage-grouse using 2 forms of telemetry (VHF, n = 330; GPS, n = 24) during the 13-yr study period. Sample sizes used for demographic parameter estimates varied across life stages (s, n = 354; c, n = 207; ns, n = 280; h, n = 126; cs, n = 113; np2, n = 49).

The derived posterior estimates of λ across the Bi-State region, as a whole, indicated that the sage-grouse population trend over the 13-yr study period did not exhibit evidence of a general pattern of decrease or increase (median = 0.98; 95% CRI = 0.69–1.25), but instead showed evidence of cycling within an otherwise stable trajectory during the 13 yr of study (Figure 1A). Based on the Bi-State–wide posterior distribution of λ for this period, the following probabilities were calculated: (1) increased population, 42.6%; (2) stable population, 3.1%; and (3) decreased population, 55.4%. Thus, the Bi-State as a whole showed similar evidence for the odds of increase to decrease (odds ratio; Figure 1B). We found that the estimates of SSM (without demographic data) and DMM (without count data) were generally consistent in the pattern of cycling, displaying similar peaks and valleys throughout the 13-yr period (Figure 2). The SSM exhibited the most amount of annual stochasticity, whereas the DMM revealed the least. Posterior distributions for the IPM (demographic and count data) exhibited the least amount of uncertainty (differences between 25^th and 75^th quartiles), whereas the DMM revealed the most. Within years of disagreement between estimates from the SSM and DMM, estimates from the IPM were generally evenly split between the model types.

FIGURE 1.

Average number of males per lek (A: dashed line represents observed counts and the solid line represents apparent abundance N) and log odds ratio (B: OR[increase:decrease]) of Greater Sage-Grouse population trend (dashed horizontal line represents a stable population and solid line is overall average for the 13-yr period) within the Bi-State Distinct Population Segment during 2003–2015. Gray shading represents 95% credible limits.

FIGURE 2.

Comparison in posterior distributions of rate of change in apparent abundance between estimates from Integrated Population Model (IPM), Demographic Matrix Model (DMM; telemetry data only), and State-Space Model (SSM; observed count data only) for data collected on Greater Sage-Grouse within the Bi-State Distinct Population Segment, California and Nevada, 2003–2015. Plot displays full range of posterior (dashed lines; minimum to maximum), first and third quartile (lower and upper edge of rectangle, respectively), and median value (horizontal line within rectangle).

Although the Bi-State, as a whole, exhibited a cyclic pattern with evidence of stability across the 13-yr period, we detected variation among years and among subpopulations. The average number of observed males per lek ranged from 4.5 (Parker Meadows) to 29.3 (Long Valley; Table 1). Bodie Hills had the greatest average λ (1.03, 95% CRI = 0.76–1.34) among all subpopulations and Parker Meadows had the least average λ (0.90, 95% CRI = 0.55–1.41). Most populations showed equal odds of decrease to increase (Figure 3A–F), however the odds of decrease for the Parker Meadows subpopulation were 2.5 times greater than the odds for an increase (Figure 3E). Nevertheless, our results indicate that the Bi-State DPS in its entirety was relatively stable during the 13-yr time frame, despite relatively short-term cycles and strong evidence for a declining trend in the Parker Meadows subpopulation. In other words, the Parker Meadows subpopulation, which consisted of 3 leks (one active) with <5 males on average (Table 1), had relatively little influence on the overall average population trend for the entire DPS.

FIGURE 3.

Log odds ratio (increase:decrease) of population trend evidence for Greater Sage-Grouse subpopulations Pine Nuts (A), Desert Creek (B), Fales (C), Bodie Hills (D), Parker Meadows (E), and Long Valley (F) within the Bi-State Distinct Population Segment, California and Nevada, 2003–2015. The dashed horizontal line represents a stable population. The solid horizontal line is overall average for the 10-yr period. Vertical bar represents final observations for all data types (right of bar is 5-yr model projections).

Precipitation Effects on λ

We investigated 20 candidate models that reflected seasonal precipitation effects (amount and change) on λ. The Δprecip during spring, summer, and fall months garnered the most support from the data in explaining more variation in λ than any other seasonal effect model (model 1, Table 3). This model was 7.2 WAIC units less than a null model (intercept only) and more parsimonious than a model consisting of precipitation across the entire year (model 2, ΔWAIC = 4.45). Change in λ cycled with Δprecip during the 13-yr period (Figure 4). The greatest negative change in precipitation occurred from 2011 to 2012, which corresponded proportionally with the greatest negative change in λ from 2012 to 2013. Similarly, the substantial decline in precipitation from 2007 to 2008 corresponded with a decline in λ from 2008 to 2009. Based on the estimated median of the posterior distribution, a 50% increase in precipitation corresponded to a 15.5% (95% CI = 5.4–26.9) increase in λ the following year. A model with the effect of Δprecip outperformed one with amount of precipitation while comparing among the same seasonal effect (model 5, ΔWAIC = 5.31). The effect of Δprecip varied across subpopulations (Figure 5). For example, the subpopulation with the greatest model support was Bodie Hills, where 50% increase in precipitation corresponded to 17.3% (95% CI = 11.0–24.5) increase in λ the following year. Parker Meadows had the weakest association, where a 50% increase in precipitation was associated with a 10.6% increase (95% CI = 3.9–18.1) in λ. All subpopulations cycled with Δprecip throughout the 13-yr period except Parker Meadows, where λ was decoupled from precipitation patterns during 2008 and failed to recover as precipitation increased (Figure 6). Other populations showed similar signs of correlation with Δprecip following the 2013 drought (Figure 6). A secondary model, the amount of precipitation across all months of the year, evidenced less support (2.60 WAIC units less than a null model). Based on a median estimate, an increase of 28.5 mm from one year to the next was associated with an increase in λ by 6.2% (95% CI = −0.2 to 12.8).

TABLE 3.

Comparison of seasonal precipitation effects on population growth rate of Greater Sage-Grouse within the Bi-State Distinct Population Segment, California and Nevada, 2003–2015.

FIGURE 4.

Relationship between changes in population growth (λ; solid line) and precipitation (▵; dashed line) between successive years for Greater Sage-Grouse within the Bi-State Distinct Population Segment, California and Nevada, 2003–2015. Population growth (per capita) and change in precipitation are on a logarithmic scale. Primary x axis reflects year of population growth change and secondary x axis represents year of precipitation change to illustrate the 1-yr lag effect. Horizontal dashed line reflects no change for either variable between years.

FIGURE 5.

Finite population rate of change (λ) as a function of change in precipitation (▵) between years within the Bi-State Distinct Population Segment (BS), as a whole, and among different subpopulations (BH = Bodie Hills, FA = Fales, DC = Desert Creek, LV = Long Valley, PN = Pine Nuts, PA = Parker Meadows) for Greater Sage-Grouse within the Bi-State Distinct Population Segment, California and Nevada, 2003–2015.

FIGURE 6.

Apparent abundance N (solid line) related to precipitation (dashed line) for Bodie Hills (A), Fales (B), Desert Creek (C), Long Valley (D), Parker Meadows (E), and Pine Nuts (F) for Greater Sage-Grouse within the Bi-State Distinct Population Segment, California and Nevada, 2003–2015. Primary x axis reflects year of population growth change and secondary x axis represents year of precipitation to illustrate the 1-yr lag effect. Gray shade represents 95% credible limits.

DISCUSSION

Our results indicate that population abundance of sage-grouse across the Bi-State metapopulation did not appear to increase or decrease throughout the 13-yr period and could be characterized as stable. However, variability among subpopulation trends indicated differing intrinsic and extrinsic factors influencing the subpopulations. Sage-grouse populations generally exhibit short-term oscillations (approximately 6–10 yr) in population abundance (Fedy and Doherty 2011, Garton et al. 2015). Thus, our 13-yr time series seemed adequate to encompass potential oscillations in abundance throughout a complete sage-grouse population cycle. In particular, the Bodie Hills subpopulation exhibited the strongest population growth, whereas Pine Nuts and Parker Meadows declined to very low numbers. Fales, Desert Creek, and Long Valley appeared stable throughout the 13-yr period.

The growth rate of the overall metapopulation cycled during relatively short periods and was correlated with climatic patterns. Specifically, annual changes in λ and the amount of precipitation were well correlated (Figure 4). In comparison to all other seasonal combinations evaluated, annual changes in precipitation during spring, summer, and fall were most strongly associated with variation in population growth. This effect of non-winter precipitation on λ across the Bi-State is consistent with results observed at larger spatial extents (e.g., Great Basin; Coates et al. 2016). The period of severe drought toward the end of the time series (peaking in 2013) allowed better estimation and insight to how climatic variation can influence population dynamics of sage-grouse.

Productivity of cold-desert sagebrush ecosystems is tied strongly to variation in precipitation (Noy-Meir 1973), and other studies have reported links between precipitation and sage-grouse population vital rates (Moynahan 2004, Blomberg et al. 2012, Guttery et al. 2013, Blomberg et al. 2014b, Gibson et al. 2017), and between precipitation and count-based changes in abundance across large spatiotemporal scales (Coates et al. 2016). Donnelly et al. (2016) demonstrated close spatial associations between sage-grouse and the distributions of mesic resources. Collectively, research indicates that cumulative precipitation during the sage-grouse reproductive period (spanning nesting to brood break-up) is important for delaying plant senescence and desiccation of mesic resources (i.e. upland springs and seeps) that provide critical resources during the brood-rearing stage as the summer progresses. However, the only vital rates influencing fecundity that were best explained by random effects (and would therefore indicate tracking of precipitation) were second nest propensity, clutch size, and hatchability. Chick survival, which would be most indicative of precipitation effects during late summer when drought conditions can be most exasperated, was best explained by subpopulation random effects. We certainly do not discount potential effects of precipitation during the breeding season on fecundity (Blomberg et al. 2014a, Blomberg et al. 2017); our failure to detect a year, or subpopulation by year, effect on chick survival likely was influenced by limited samples across subpopulations and years. However, although annual survival varied by subpopulation, the overall lowest survival estimates were during 2007 and 2012, which were the 2 driest years in the time series. Drought can negatively impact populations of sage-grouse through many mechanisms (Blomberg et al. 2012, Donnelly et al. 2016), yet a likely explanation for how drought precipitated sage-grouse declines in our study relates to the hypothesis that lack of water during the breeding season causes more movement by sage-grouse as they search the landscape for increasingly scarce mesic resources (Gibson et al. 2017). Furthermore, increased movements among individuals decreased their probability of survival (Prochazka et al. 2017), likely as a function of increased vulnerability to predation.

Although results from our study at the metapopulation extent implicate extrinsic (and density-independent) variation in annual precipitation as a strong driver of cyclic abundance, Fedy and Doherty (2011) implied that intrinsic feedbacks from delayed density dependence can drive sage-grouse population cycles. This is likely a functional response of sage-grouse due to predation or declining resource availability to sage-grouse as populations approach carrying capacity. Consistent with Garton et al. (2015) and Coates et al. (2016), who described Gompertz-type effects on λ, we found that probability of renesting and chick survival were influenced by a 1-yr lag effect of density. Thus, it is this lag effect that causes these demographics to be greatest following years with relatively low sage-grouse population abundance. These studies attempted to account for both density-dependent and density-independent effects, yet patterns illustrated in Figures 4 and 5 indicate a more recent paradigm—that the strength of density-dependence can be modified by interactions with environmental covariates, which then influences population stability (Saether 1997, Ahrestani et al. 2016). In our study, differences in the magnitude of sage-grouse population response to annual change in precipitation indicated that populations at lower densities (owing to declining λ) can respond more rapidly than those at higher densities to equivalent amounts of seasonal precipitation. Whereas per capita resources during periods of higher density and higher precipitation are essentially capped near carrying capacity, more resources are available to fewer individuals during periods of high rainfall following drought. Moreover, our results indicate that precipitation needs to be approximately 20% greater than average for population recovery following drought, consistent with results from the Great Basin during a longer time series and in the absence of wildfire (Coates et al. 2016). Hence, interactions between density-dependent and density-independent factors help fuel “boom-and-bust” cycles that are ubiquitous among sage-grouse populations, and particularly those inhabiting more arid environments such as the Great Basin.

Hierarchical models, such as IPMs, are well suited for disentangling local and regional drivers in cross-scale interactions (Soranno et al. 2014), and results from our analysis provided strong evidence of these interactions across spatial scales in the Bi-State DPS. We also observed large variation in the response to precipitation at the subpopulation extent, and that widespread effects of precipitation at the metapopulation extent were decoupled from the subpopulation extent for some subpopulations (Figures 5 and 6). For example, the relationship between changes in precipitation and λ was weakest for Parker Meadows and strongest for Bodie Hills, relative to the effect at the metapopulation extent (Figure 5). In fact, consistently low abundance at Parker Meadows was strongly invariant to precipitation throughout the time series (Figure 6). In addition, subpopulations whose median effects differed greatest from the metapopulation effect (in ascending order: Desert Creek, Long Valley, Pine Nuts, and Parker Meadows; Figure 5) failed to respond to more typical strong effects of high precipitation. This demonstrated increasingly greater odds of population decline compared with odds of population increase. In contrast, the subpopulation at Bodie Hills appeared buffered from the effects of drought, possibly because it responded more strongly to periods of high precipitation. Alternatively, this subpopulation was at higher elevation where ground moisture was likely more stable throughout the season and precipitation is greater.

Demographic rates among subpopulations can yield more detailed insight into the mechanisms that result in decoupling. Random effects analysis within the demographic models identified that hatchability influenced fecundity and ultimately λ within Parker Meadows, a subpopulation that is isolated and small. Estimated parameters for hatchability in this subpopulation were approximately one-half of those estimated among the others. Hatchability at Parker Meadows was also markedly lower than range-wide average estimates derived for Greater Sage-Grouse average (> 0.92; Taylor et al. 2012), and other species of prairie grouse such as Sharp-tailed Grouse (Tympanuchus phasianellus; > 0.90; McNew et al. 2017) and Greater Prairie-Chicken (Tympanuchus cupido; > 0.80; Peterson and Silvy 1996, McNew et al. 2012). Field observations concluded that decreased hatchability was caused by low egg fertility rates. Parker Meadows is likely experiencing an Allee effect, perhaps explained by limited opportunities for females to copulate with the few males that remain in that population (Stephens et al. 1999). Alternatively, an Allee effect could be explained by lower mitochondrial and nuclear genetic variation compared with other sage-grouse (Oyler-McCance et al. 2014) as a result of genetic drift (Luque et al. 2016). In addition, second nest propensity estimates help explain observed λ and low abundance at Parker Meadows and Pine Nuts.

Responses of different subpopulations to changes in precipitation elucidate potential reasons for decoupling from trends at the metapopulation scale. For example, changes in abundance at Long Valley were correlated with relative changes in precipitation through peak drought conditions in 2013, but remained depressed thereafter, despite a subsequent increase in precipitation in 2015. Subpopulations Pine Nuts and Desert Creek displayed similar patterns. Failure to respond positively to the cessation of drought conditions likely contributed to greater odds of decrease versus increase for these subpopulations. Along with Parker Meadows, leks comprising these subpopulations generally were situated in narrow bands of mesic and high-elevation habitat surrounded by more expansive arid areas at low elevation, using lek locations and soil moisture and temperature maps from Maestas et al. (2016). In contrast, leks comprising Bodie Hills and Fales were surrounded by more contiguous expanses of mesic and high-elevation habitat, and have either recovered (in the case of Fales) or were strongly buffered (in the case of Bodie Hills) from the effects of drought.

Preliminary results of subpopulation and metapopulation growth rates for Bi-State used a shorter data time series that ended prior to the height of the drought in 2013 (Coates et al. 2014). During the shorter time period, subpopulations at Pine Nuts and Long Valley had higher odds of increase than decrease, and subpopulations at Fales had slightly higher odds of decrease than increase. When the IPM described in this study was expanded to include longer time series that encompassed the full extent and partial recovery from the drought, odds of increase versus decrease reversed for these subpopulations. Interestingly, the lower annual survival described for Pine Nuts and Long Valley in relation to all other subpopulations during the major drought year of 2013 may point to increased risk of predation as adult sage-grouse search for scarce mesic resources. Regardless of the proximate mechanism, these patterns underscore the importance of multiyear data sets for evaluating vital rates and abundances through time series that encompass extremes in climatic variability in semi-arid ecosystems (i.e. the Bi-State DPS). Accordingly, the winter of 2015 to 2016 was one of the wettest on record, and ongoing intensive monitoring of Bi-State sage-grouse planned through 2020 will allow better understanding of demographic responses to periods of extreme precipitation at the subpopulation and metapopulation extents using the methods we described herein.

We compared estimates of λ from DMM (demographic data only), SSM (count data only), and IPM (combined data) because our unique dataset is “vital-rate rich”; the DMM is informed by data that estimate 9 different vital rates describing annual survival and fecundity, and encompasses the life-cycle of Greater Sage-Grouse (Dahlgren et al. 2016). In contrast to lek count data that informed the SSM, demographic data were not available for all years and across all vital rates. Demographic model estimates of λ reflect more within-year variation and generally have wider posterior distributions for 2 main reasons: (1) estimates of uncertainty for each component vital rate propagate at each life stage step, and (2) estimates shrink toward a grand mean during years when demographic information is relatively sparse. In contrast, the state-space model estimates of λ might be more precise, but could be biased because they are only based on counts subject to observation error. In cases where posterior distributions of λ from the demographic model are wide, the IPM appears to “favor” information from the SSM, and median estimates of λ from the SSM and IPM aligned closely. When demographic data were available, the IPM estimate relied on information from both sub-component models. In both cases, however, the demographic model provided information to estimate λ. Hence, sage-grouse populations can be monitored economically with lek counts, but estimates of λ can be refined even with sparse demographic data that can also help identify specific vital rates most responsible for variation in λ at the subpopulation extent. Dahlgren et al. (2016) conducted a similar study with life-cycle models instead of IPMs, and also found that demographic models helped refine estimates of λ, and matched changes in lek counts through time. Perhaps most importantly, patterns illustrated by our comparison of demographic, state-space, and IPM outputs of λ further the notion proposed by Nowak (2015)—that the fear of incomplete data should not hamstring application of IPMs.

One caveat is that IPMs assume independence among datasets for robust parameter estimation (Kéry and Schaub 2012). In practice, this assumption is frequently violated for demographic sub-model components, as it was for our study. This is because the number of monitored individuals required to estimate each vital rate independently is often logistically daunting and financially infeasible. However, any violations of nonindependence of data used in our study likely had little influence on parameter bias (Besbeas et al. 2009, Abadi et al. 2010). Additionally, females were used for vital rate estimation and lek counts were based on unmarked male sage-grouse. Thus, the datasets we used to inform demographic (as a whole) and abundance sub-models can be assumed to be independent.