Local and population-level responses of Greater sage-grouse to oil and gas development and climatic variation in Wyoming
- Published
- Accepted
- Received
- Academic Editor
- Xavier Harrison
- Subject Areas
- Conservation Biology, Ecology, Natural Resource Management, Environmental Impacts
- Keywords
- Lek counts, Greater sage-grouse, Pacific Decadal Oscillation, Climate, Population dynamics, Oil and gas
- Copyright
- © 2018 Ramey et al.
- Licence
- This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
- Cite this article
- 2018. Local and population-level responses of Greater sage-grouse to oil and gas development and climatic variation in Wyoming. PeerJ 6:e5417 https://doi.org/10.7717/peerj.5417
Abstract
Background
Spatial scale is important when studying ecological processes. The Greater sage-grouse (Centrocercus urophasianus) is a large sexually dimorphic tetraonid that is endemic to the sagebrush biome of western North America. The impacts of oil and gas (OAG) development at individual leks has been well-documented. However, no previous studies have quantified the population-level response.
Methods
Hierarchical models were used to estimate the effects of the areal disturbance due to well pads as well as climatic variation on individual lek counts and Greater sage-grouse populations (management units) over 32 years. The lek counts were analyzed using generalized linear mixed models while the management units were analyzed using Gompertz population dynamic models. The models were fitted using frequentist and Bayesian methods. An information-theoretic approach was used to identify the most important spatial scale and time lags. The relative importance of OAG and climate at the local and population-level scales was assessed using information-theoretic (Akaike’s weights) and estimation (effect size) statistics.
Results
At the local scale, OAG was an important negative predictor of the lek count. At the population scale, there was only weak support for OAG as a predictor of density changes but the estimated impacts on the long-term carrying capacity were consistent with summation of the local impacts. Regional climatic variation, as indexed by the Pacific Decadal Oscillation, was an important positive predictor of density changes at both the local and population level (particularly in the most recent part of the time series).
Conclusions
Additional studies to reduce the uncertainty in the range of possible effects of OAG at the population scale are required. Wildlife agencies need to account for the effects of regional climatic variation when managing sage-grouse populations.
Introduction
If we study a system at an inappropriate scale, we may not detect its actual dynamics and patterns but may instead identify patterns that are artifacts of scale. Because we are clever at devising explanations of what we see, we may think we understand the system when we have not even observed it correctly (Wiens, 1989).
Effective conservation of a species requires an understanding of how human activities influence its distribution and abundance. Although much of science proceeds by experimental studies to understand the causal links between actions and responses, ethical, practical and statistical considerations typically prevent population-level experiments on species of concern. Consequently, many conservation-based ecological studies are forced to infer the population-level consequences of anthropogenic alterations from local gradients (Fukami & Wardle, 2005) in density (Gill, Norris & Sutherland, 2001), movement, habitat use, physiology, genetics, reproductive success or survival. However, local gradients may not accurately predict the population-level response (Gill, Norris & Sutherland, 2001; Fodrie et al., 2014).
The Greater sage-grouse (Centrocercus urophasianus, hereafter sage-grouse) is a large sexually dimorphic tetraonid that is endemic to the sagebrush (Artemisia spp.) biome of western North America (Knick & Connelly, 2011). Each spring, adult males aggregate in open areas called leks where they display for females. Fertilized females then nest on the ground among the sagebrush (Holloran & Anderson, 2005). Initially, the chicks feed on insects before switching to forbs. The adults predominantly feed on sagebrush, especially in the winter. Most males begin lekking 2 years after hatching. Mean peak counts of males on leks are commonly used as an abundance metric (Connelly & Braun, 1997; Doherty, Naugle & Evans, 2010; Fedy & Aldridge, 2011).
A multitude of studies have reported local negative effects of oil and gas (OAG) development on sage-grouse lek counts, movement, stress-levels and fitness components. The most frequently-reported phenomenon is the decline in lek counts with increasing densities of well pads (Walker, Naugle & Doherty, 2007; Doherty, Naugle & Evans, 2010; Harju et al., 2010; Green, Aldridge & O’Donnell, 2016). Reductions in fitness components such as lower nest initiation rates (Lyon & Anderson, 2003) and lower annual survival of yearlings reared in areas where OAG infrastructure is present (Holloran, Kaiser & Hubert, 2010) have been detected using radio-tracking. The development of global positioning system telemetry methods has facilitated the fitting of more sophisticated and realistic spatially-explicit habitat use models which suggest that nest and brood failure is influenced by proximity to anthropogenic features (Dzialak et al., 2011). More recently, experimental studies have suggested that noise alone can reduce lek attendance (Blickley et al., 2012b) and increase stress hormones (Blickley, Blackwood & Patricelli, 2012a).
However, to date no-one has examined whether sage-grouse population-level responses are consistent with the local studies. Although Green, Aldridge & O’Donnell (2016) state that they modeled sage-grouse populations, they use their population dynamic models to analyze the effects of OAG on changes in abundance at individual leks. Even authors such as Walker, Naugle & Doherty (2007) and Gamo & Beck (2017) who analyzed aggregations of leks, group their leks by level of OAG development as opposed to population boundaries.
Although it has received less attention than OAG, climatic variation has also been shown to influence sage-grouse lek counts, survival, clutch size and nesting success (Blomberg et al., 2012, 2014, 2017; Coates et al., 2016; Gibson et al., 2017). This is not surprising, as there is a long and ecologically important history of studies on the influence of climatic variation on the population dynamics of tetraonids (Moran, 1952, 1954; Ranta, Lindstrom & Linden, 1995; Lindström et al., 1996; Cattadori, Haydon & Hudson, 2005; Ludwig et al., 2006; Kvasnes et al., 2010; Selås et al., 2011, Viterbi et al., 2015; Ross et al., 2016). Consequently, the current study also includes annual variation in regional climate as a potential predictor of sage-grouse population dynamics.
Previous studies of the effect of climatic variation on sage-grouse have used local temperature and precipitation data with mixed results (Blomberg et al., 2012, 2014, 2017; Green, Aldridge & O’Donnell, 2016; Coates et al., 2016; Gibson et al., 2017). However, large-scale climate indices often outperform local data in predicting population dynamics and ecological process (Stenseth et al., 2002; Hallett et al., 2004). The Pacific Decadal Oscillation (PDO), which is derived from the large-scale spatial pattern of sea surface temperature in the North Pacific Ocean (Mantua et al., 1997), is potentially the most important climatic process influencing the sagebrush biome (Neilson et al., 2005). Consequently, the PDO index was chosen as the climate indicator.
Wyoming was selected for the current study because it contains approximately 37% of the recent range-wide population of sage-grouse (Copeland et al., 2009; Fedy et al., 2012), has experienced substantial levels of OAG development dating to the late 1800s (Braun, Oedekoven & Cameron, 2002) and because the lek location and count data were available for research.
Methods
Data preparation
Sage-grouse data
When multiple counts exist for the same lek in a single year, almost all authors take the maximum count (Holloran, 2005; Walker, Naugle & Doherty, 2007; Harju et al., 2010; Fedy & Aldridge, 2011; Fedy & Doherty, 2011; Garton et al., 2011, 2015; Blickley et al., 2012b; Blomberg et al., 2013; Davis et al., 2014; Coates et al., 2016; Fremgen et al., 2016; Monroe, Edmunds & Aldridge, 2016; Green, Aldridge & O’Donnell, 2016). The justification for using the maximum count is articulated by Garton et al. (2011) who state that,
…counts over the course of a single breeding season vary from a low at the beginning of the season, to peak in the middle, followed by a decline to the end, which necessitates using the maximum count from multiple counts across the entire season as the index.
However, as noted by Johnson & Rowland (2007), this results in a substantial upward bias at leks with multiple counts. To avoid this bias, several alternative approaches are available: exclude early and late counts and then either include the repeated counts in the model (Gregory & Beck, 2014) or take the mean of the repeated counts and/or explicitly model the change in attendance through time (Walsh et al., 2004) as is done for spawning salmon (Hilborn, Bue & Sharr, 1999). We excluded early and late counts and took the rounded mean of the repeated counts. However, as discussed below we also assessed the sensitivity of the results to the use of the rounded mean as opposed to maximum count.
The sage-grouse lek count and location data were provided by the State of Wyoming. After excluding male lek counts with unknown counts or dates or those before 1985 there were 88,771 records. To reduce potential biases, only the most reliable male lek counts were included in the analyses. In particular, only ground counts from leks that were checked for activity and were part of a survey or count were included (as per Wyoming Game and Fish guidelines). This reduced the number of records to 79,857. To ensure counts were close to the peak (see above), only data that were collected between April 1st and May 7th were included. This reduced the number of records to 65,439. Finally, lek counts for which the number of individuals of unknown sex were ≥5% of the number of males (suggesting unreliable identification) were excluded which left a total of 42,883 records. The leks with at least one remaining count are mapped in Fig. 1 and the associated mean male lek counts are plotted in Fig. 2.
The State of Wyoming utilizes eight regional sage-grouse working groups to facilitate local population management and data reporting, including lek counts and hunting harvest (Fig. 1, Christiansen & Belton, 2017). For the purposes of the current study, we also treat these working groups as if they are separate populations. The population densities (males per lek) were calculated by averaging the mean counts for individual leks for each working group in each year.
Oil and gas data
Wyoming Oil and Gas Conservation Commission’s conventional, coal-bed and injection well pad location and production data were downloaded from the Wyoming Geospatial Hub (http://pathfinder.geospatialhub.org/datasets/) at May 25, 2018 02:13 UTC. Well pads without a provided spud date were excluded as were well pads constructed before 1900 or after 2016. The included well pads are mapped in Fig. 1.
The intensity of OAG development was quantified in terms of the proportional areal disturbance due to well pads within a specific distance of the leks. The areal disturbance was calculated at lek distances of 0.8, 1.6, 3.2 and 6.4 km with the areal disturbance of each well pad considered to have a radius of 60 m (Green, Aldridge & O’Donnell, 2016). The annual areal disturbances for individual leks with lek counts at 3.2 km are plotted in Fig. 3.
Climatic data
The PDO index (Trenberth & Hurrell, 1994; Mantua et al., 1997) data were queried from the rpdo R package (Fig. 4).
Statistical analysis
Local models
The individual lek counts were analyzed using generalized linear mixed models (Bolker et al., 2009) with the standardized areal disturbance due to OAG and the PDO index as fixed effects and year and lek as random effects. The areal disturbance and PDO index were standardized (centered and divided by the standard deviation (SD)) to facilitate comparison.
More formally, the lek count model is described by the following equations (1) $${M}_{i,y}~\text{NegBin}({\text{\mu}}_{i,y},\varphi )$$ (2) $$\mathrm{log}({\mu}_{i,y})={\text{\beta}}_{0}+{\text{\beta}}_{A}\cdot {\text{AREA}}_{i,y}+{\text{\beta}}_{P}\cdot {\text{PDO}}_{y}+{\text{\alpha}}_{L}{}_{i}+{\text{\alpha}}_{Y}{}_{y}$$ (3) $${\text{\alpha}}_{L}{}_{i}~\text{Normal}(0,{\text{\sigma}}_{L})$$ (4) $${\text{\alpha}}_{Y}{}_{y}~\text{Normal}(0,{\text{\sigma}}_{Y})$$ where M_{i,y} is the rounded mean count of males for the ith lek in the yth year, β_{A} and β_{P} are the fixed effects of the standardized areal disturbance due to well pads (AREA_{i,y}) and PDO index (PDO_{y}) on the expected count (μ_{i,y}), σ_{L} and σ_{Y} are the SDs of the random effects of lek and year. In our parameterization of the negative binomial the parameter ϕ controls the overdispersion scaled by the square of μ, i.e., (5) $$\text{SD}[M]=\sqrt{\text{\mu}+\varphi \cdot {\text{\mu}}^{2}}$$ Key model parameters are also described in Table 1.
Parameter | Description |
---|---|
β_{0} | The intercept for the log lek count or log population density. |
β_{D} | The effect of population density on β_{0}. |
β_{P} | The effect of the standardized Pacific Decadal Oscillation index on β_{0}. |
β_{A} | The effect of the standardized areal disturbance due to well pads on β_{0}. |
ϕ | The overdispersion term. |
${\text{\sigma}}_{{\text{\eta}}_{g,y}}$ | The SD of the process error. |
β_{N} | The intercept for the SD of the process error. |
β_{L} | The effect of the log of the number of leks counted on β_{N}. |
σ_{G} | The SD of the random effect of working group on β_{D}. |
σ_{L} | The SD of the random effect of lek on β_{0}. |
σ_{Y} | The SD of the random effect of year on β_{0}. |
To identify the most important spatial scale (distance from each lek when calculating the areal disturbance) and temporal lags, a total of 64 models were fitted to the lek count data representing all combinations of the four lek distances (0.8, 1.6, 3.2 and 6.4 km) and independent lags of 1–4 years in the areal disturbance (Walker, Naugle & Doherty, 2007; Doherty, Naugle & Evans, 2010; Harju et al., 2010; Gregory & Beck, 2014) and PDO index. The relative importance of each spatial scale and temporal lag as a predictor of individual lek counts was assessed by calculating its Akaike’s weight (w_{i}) across all 64 models (Burnham & Anderson, 2002).
Once the model with the most important spatial scale and temporal lags was identified, the relative importance of β_{A} and β_{P} was quantified by calculating their Akaike’s weights across the selected full model and the three reduced variants representing all combinations of the two parameters (Burnham & Anderson, 2002) and by calculating their effect sizes with 95% confidence/credible intervals (Bradford, Korman & Higgins, 2005; Claridge-Chang & Assam, 2016). The effect sizes, which represent the expected percent change in the lek count with an increase in the predictor of one SD, were calculated for the final full model and by averaging across all four models (Burnham & Anderson, 2002; Turek, 2015). In addition, the conditional and marginal R^{2} values were calculated for the final full model on the observational (original) scale (Nakagawa & Schielzeth, 2013).
Population models
The calculated annual population densities (mean males per lek) in each working group were analyzed using Gompertz population dynamic models (GPDMs) (Garton et al., 2011) with the standardized areal disturbance and PDO index as fixed effects and year and group as random effects. GPDMs were used because they incorporate density-dependence (Dennis et al., 2006; Knape & De Valpine, 2012) and have performed well in explaining rates of change for sage-grouse in general and for Wyoming sage-grouse in particular (Garton et al., 2011).
The population model is described by the following equations (6) $$\mathrm{log}({M}_{g,y})~\text{Normal}(\mathrm{log}({\text{\mu}}_{g,y}),{\text{\sigma}}_{\eta}{}_{g,y})$$ (7) $$\mathrm{log}({\text{\mu}}_{g,y})={\text{\beta}}_{0}+({\text{\beta}}_{D}+1+{\text{\alpha}}_{G}{}_{g})\cdot \mathrm{log}({M}_{g,y-1})+{\text{\beta}}_{A}\cdot {\text{AREA}}_{g,y}+{\text{\beta}}_{P}\cdot {\text{PDO}}_{y}+{\text{\alpha}}_{Y}{}_{y}$$ (8) $${\text{\alpha}}_{G}{}_{g}~\text{Normal}(0,{\text{\sigma}}_{G})$$ (9) $${\text{\alpha}}_{Y}{}_{y}~\text{Normal}(0,{\text{\sigma}}_{Y})$$ (10) $$\mathrm{log}({\text{\sigma}}_{\eta}{}_{g,y})={\text{\beta}}_{N}+{\text{\beta}}_{L}\cdot \mathrm{log}({\text{LEKS}}_{g,y})$$ where M_{g,y} is the density at the gth group in the yth year, μ_{g,y} is the expected density, β_{D} is the typical density-dependence and α_{Gg} is the group-level random effect on the density-dependence, ${\text{\sigma}}_{{\text{\eta}}_{g,y}}$ is the expected process error (Dennis et al., 2006), β_{N} is the intercept for the log process error and β_{L} is the effect of the log of the number of leks surveyed (LEKS_{g,y}) on β_{N}. The other terms are approximately equivalent to those in the lek count model. The equivalence is only approximate as the terms in the population model act on the change in density (as opposed to density).
The carrying capacity, which represents the long-term expected density around which a population fluctuates (Dennis et al., 2006), is given by (11) $$\mathrm{log}({N}_{\infty})=\frac{-({\text{\beta}}_{0}+{\text{\beta}}_{A}\cdot {\text{AREA}}_{g,y}+{\text{\beta}}_{P}\cdot {\text{PDO}}_{y})}{{\text{\beta}}_{D}+{\text{\alpha}}_{G}{}_{g}}.$$
Preliminary analyses considered Gompertz state-space population models (Dennis et al., 2006; Knape & De Valpine, 2012) which estimate both process and observer error (Maunder, Deriso & Hanson, 2015). However, the models were unable to reliably estimate both error terms. As the observer error was estimated to be smaller than the process error and because ignoring process error can bias Akaike’s information criterion based tests toward incorrectly accepting covariates (Maunder, Deriso & Hanson, 2015), we followed Garton et al. (2011) in assuming no observer error. The preliminary analyses indicated that fixing the observer error at zero had little effect on the results. The process error was allowed to vary by the logarithm of the number of leks surveyed, which varied from one to 214, to account for the additional stochasticity associated with smaller populations and/or lower coverage.
The primary question this study attempts to answer is whether the sage-grouse population-level responses to OAG are consistent with the local studies. Consequently, the average areal disturbance in each working group was calculated at the spatial scale that was most important in the local analyses. However, as the timing of effects could differ between the local count models and the population dynamic models, the Akaike’s weight for each lag of 1–4 years in the areal disturbance and 1–4 years in the PDO index was calculated across the 16 models representing all lag combinations. Once the full model with the most important temporal lags had been identified, the relative importance of β_{A} and β_{P} was once again quantified from their effect sizes with 95% confidence intervals and their Akaike’s weights across the full model and the three reduced variants.
To assess the sensitivity of the population model outputs to the inputs, a local, qualitative sensitivity analysis (Pianosi et al., 2016) was conducted. More specifically, the effect sizes in the full model were estimated using (1) the maximum of the repeated counts (Garton et al., 2011) as opposed to the mean (Johnson & Rowland, 2007) and (2) the data from 1997 and 2005 onward. The former period was used during preliminary analyses due to the availability of hunter-harvested wing count data (Braun & Schroeder, 2015). From 2005 onward regulatory and technological developments were introduced to reduce the impacts of OAG on sage-grouse and to make operations more efficient (Applegate & Owens, 2014). The effect sizes were adjusted for any differences in the SDs of the data. Once again, the conditional and marginal R^{2} values were calculated for the final full model on the observational scale (Nakagawa & Schielzeth, 2013).
Predicted population impacts
To examine whether the population-level responses are consistent with the lek-level results, the expected effect of OAG on the long-term mean densities in each working group was calculated for the local and population models. In the case of the local model, the predicted population impacts represent the percent difference in the sum of the expected counts for the observed levels of OAG versus no OAG across all leks in the working group after accounting for annual and climatic effects. In the case of the population model, the predicted population impacts represent the percent difference in the expected carrying capacity with the observed levels of OAG versus no OAG accounting for annual and climatic effects.
Statistical methods
For reasons of computational efficiency, the initial 64 local and 16 population-level models were fit using the frequentist method of maximum likelihood (ML; Millar, 2011). The Akaike’s weights were calculated from the marginal Akaike’s information criterion values corrected for small sample size (Burnham & Anderson, 2002; Vaida & Blanchard, 2005; Greven & Kneib, 2010). Model adequacy was assessed by plotting and analysis of the standardized residuals from the final full ML model (Burnham & Anderson, 2002) with the most important spatial scale and lags. As both the local and population models used log-link functions, the effect sizes (percent change in the response for an increase in one SD) were calculated from exp(β)—1 where β is the fixed effect of interest or its upper or lower confidence limit. The ML effect sizes were calculated from the full model and averaged across the full model and three reduced variants (Lukacs, Burnham & Anderson, 2010).
To allow the predicted population impacts to be estimated with credible intervals, the final full models were also fitted using Bayesian methods (Gelman et al., 2014). The prior for all primary parameters was an uninformative (Gelman et al., 2014) normal distribution with a mean of 0 and a SD of 5. A total of 1,500 MCMC samples were drawn from the second halves of three chains. Convergence was confirmed by ensuring that Rhat was ≤1.01 (Gelman et al., 2014) and the effective sample size was ≥1,000 (Brooks & Gelman, 1998) for each structural parameter.
Software
The data preparation, analysis and plotting were performed using R version 3.5.0 (R Core Team, 2017) and the R packages TMB (Kristensen et al., 2016) and rstan (Stan Development Team, 2016). The clean and tidy analysis data and R scripts are archived at https://doi.org/10.5281/zenodo.837866. The raw sage-grouse data, which provide the lek locational information, are available from the Wyoming Department of Fish and Game. The raw data are not required to replicate the analyses.
Results
Local models
The Akaike weights for the spatial scales indicate that 3.2 km is unanimously supported (w_{i} = 1.00) as the most important lek distance for predicting individual lek counts from the areal disturbance due to well pads (Table S1). The Akaike weights for the lags in the areal disturbance provided close to unanimous support for a single candidate with the lag of 1 year receiving a weight of 0.99 (Table S2). The situation with the PDO index lags was less clear-cut (Table 2), although a lag of 2 years received the majority of the support (w_{i} = 0.73). Consequently, the local model with a lek distance of 3.2 km and lags of 1 and 2 years in the areal disturbance due to well pads and the PDO index, respectively, was selected as the final model. The standardized residuals, with the exception of a small number of high outliers, were approximately normally distributed and displayed homogeneity of variance. Most leks had an expected count of male sage-grouse in the absence of OAG of approximately 10 birds (Fig. S1).
PDO lag (year) | Models | Proportion | w_{i} |
---|---|---|---|
2 | 16 | 0.25 | 0.73 |
3 | 16 | 0.25 | 0.18 |
4 | 16 | 0.25 | 0.05 |
1 | 16 | 0.25 | 0.04 |
The Akaike weights for β_{A} (w_{i} = 1) and β_{P} (w_{i} = 0.98) across the final full model and the three reduced models indicate that both are very strongly supported as predictors of individual lek counts. The effect size estimates (Fig. 5) indicate that OAG and the PDO have large negative (Fig. 6) and positive (Fig. 7) impacts of similar magnitudes (just under 20%) on the lek counts and that the estimates are insensitive to the statistical framework (ML of Bayesian) or model-averaging (Table S3). Despite the inclusion of the PDO index as an important predictor, there was still substantial remaining annual cyclical variation in the lek counts (Fig. S2) which was modeled by the random effect of year. The conditional and marginal R^{2} values were 64% and 2%, respectively.
Population models
Based on the results of the local models, the level of OAG development in each working group was calculated in terms of the average areal disturbance due to well pads within 3.2 km of each lek (Fig. 8). The Akaike weights for the lag in the areal disturbance (Table 3) were largely indifferent (0.29–0.20) although a lag of 1 year had the most support. The Akaike weights for the PDO index (Table 4) provided the majority of the support for a lag of 1 year (w_{i} = 0.53). The model predictions provided a reasonable fit to the annual mean lek counts which exhibit large cyclical fluctuations (Fig. 9). The residuals were approximately normally distributed with homogeneity of variance. The carrying capacities in the absence of OAG varied between approximately 13 males per lek in the Upper Snake River to approximately 35 males per lek in the Upper Green River (Fig. S3).
Area lag (year) | Models | Proportion | w_{i} |
---|---|---|---|
1 | 4 | 0.25 | 0.29 |
2 | 4 | 0.25 | 0.27 |
3 | 4 | 0.25 | 0.24 |
4 | 4 | 0.25 | 0.20 |
PDO lag (year) | Models | Proportion | w_{i} |
---|---|---|---|
1 | 4 | 0.25 | 0.53 |
2 | 4 | 0.25 | 0.32 |
3 | 4 | 0.25 | 0.08 |
4 | 4 | 0.25 | 0.08 |
The Akaike weights for β_{A} and β_{P} (Table 1) across the final full model and the three reduced models indicate that while the PDO index received moderate support (w_{i} = 0.74) as a predictor of population changes, there was only weak support (w_{i} = 0.41) for the areal disturbance. In addition, the effect size estimates (Fig. 10), which are sensitive to model-averaging but not the statistical framework (Table S4), indicate that while the PDO has a moderate positive influence (effect size of approximately 8%) the effect of OAG on the subsequent year’s density is relatively small (effect size of −2.5%, Fig. 10). However, despite its relatively small effect size, OAG may have a substantial effect (−20% reduction) on the long-term carrying capacity in the most impacted working groups (Fig. 11). The reasons why are discussed below. The effect of the PDO on the carrying capacity (Fig. S4) is comparable to its effect on the counts at individual leks although there is much more uncertainty (Fig. 7). As for the local models, the random effect of year accounted for substantial unexplained annual cyclical variation (Fig. S5). The effect of density on the subsequent year’s density (i.e., density-dependence) is relatively minor: in a typical working group a reduction in the density to 10 males (half the carrying capacity) results in an average of just 13 males the following year (Fig. S6). The conditional and marginal R^{2} values were 83% and 57%, respectively.
The qualitative sensitivity analysis indicates that using the maximum as opposed to the rounded mean of the repeated lek counts has a negligible effect on the effect size estimates (Fig. 12). In contrast, the estimated effect of the PDO is strongest in the most recent data.
Predicted population impacts
The predicted population impacts from the Bayesian full models indicate that although OAG has a relatively minor influence on the change in the population density, the effect on the carrying capacity could be consistent with summation of the local lek-level impacts (Fig. 13).
Discussion
Climatic variation
A key conclusion of this paper is that regional climatic variation is responsible for the inter-decadal population fluctuations experienced by sage-grouse in Wyoming. More specifically, the PDO index is an important predictor of changes in sage-grouse numbers at both the lek and population level. This is perhaps unsurprising as the PDO has previously been used, in combination with the Atlantic multi-decadal oscillation and El Nino Southern Oscillation, to predict drought, drought-related fire frequency, and precipitation trends in the western USA and Rocky Mountains (McCabe, Palecki & Betancourt, 2004; Schoennagel et al., 2007; Kitchen, 2015, Heyerdahl, Morgan & Riser, 2008).
Although the current study does not identify the causal pathways through which sea surface temperatures in the North Pacific affect the sage-grouse population dynamics, we note that in Wyoming a positive PDO correlates with cooler, wetter weather, while a negative phase tends to produce warmer, drier conditions (McCabe, Palecki & Betancourt, 2004). We also note that given the relatively poor performance of local precipitation and temperature metrics (Blomberg et al., 2012, 2014, 2017; Green, Aldridge & O’Donnell, 2016; Coates et al., 2016; Gibson et al., 2017), the causal pathways may be complex and involve other organisms such as parasites (Cattadori, Haydon & Hudson, 2005; Taylor et al., 2013). In fact, the complexity of such pathways is one of the reasons that large-scale climate indices such as the PDO often outperform local weather data in predicting population dynamics and ecological process (Stenseth et al., 2002; Hallett et al., 2004). Additional studies to assess the explanatory value of the PDO index across the species range are needed (Doherty et al., 2016). It is noteworthy that the effect of the PDO index appears to be stronger in the most recent part of the time series at least at the population level. It is also worth noting that the strength of the population-level response appears to vary between working groups (i.e., the fluctuations in population density in Bates Hole are much larger than those in the Upper Snake River) and may have a latitudinal or altitudinal component.
The finding that the PDO index is an important driver of sage-grouse abundance in Wyoming has major implications for our understanding and conservation of the species. At the very least it is expected that any long-term population trends, including those due to OAG, will be better understood in the context of the PDO (Ballard et al., 2003; McClure et al., 2012). At best, it should allow regulators to account for and predict (Stenseth et al., 2003) the effects of climatic variation on sage-grouse population fluctuations, and therefore more effectively balance conservation efforts.
Oil and gas
Ours is not the first study to consider whether local impacts potentially extend to population declines. Based on local declines, Copeland et al. (2013) estimated that sage-grouse populations in Wyoming will decrease by 14–29%, but that a conservation strategy that includes the protection of core areas could reduce the loss to between 9% and 15%, while Copeland et al. (2009) estimated that future OAG development in the western USA will cause a long-term 7–19% decline in sage-grouse numbers relative to 2007. As argued by Doherty, Naugle & Evans (2010), estimation of population-level impacts is important because it provides a biologically-based currency for quantifying the cost of OAG as well as the benefits of mitigation or conservation. Our study is the first to examine whether the actual population-level response is consistent with the local impacts.
Interpretation of the local lek-level results is relatively straightforward. The areal disturbance from OAG is a well-supported strongly negative predictor of male attendance at individual leks. The effect is much stronger at a lag of 1 year and when considering wells within a radius of 3.2 km. However, this should not be taken to imply that there are no delayed effects or disturbances from more distant wells. There is little uncertainty in the magnitude of the local declines: an areal disturbance of 3% is associated with an average local decline of 50%, while a 6% areal disturbance is associated with a local decline of 75%. When scaled up, the lek-level results suggest that by 2016 the impact of OAG was equivalent to the loss of 20% of the birds in the Northeast and 30% in the Upper Green River.
Interpretation of the population-level results is more complicated. The mean areal disturbance within 3.2 km of all the leks was a weakly supported predictor of the annual density change with a small effect size. Yet, the predicted population-level impacts were consistent with summation of the local impacts. This apparent paradox it due to three statistical phenomena. The first is that when the statistical power is low, an important variable can have low predictive value and therefore receive a low Akaike’s weight. The second is that when density dependence is weak, a small effect on the expected density change can have a more substantial impact on the long-term carrying capacity. The third is that a standardized effect size can provide a misleading summary of the scale of the possible impact if the variation is highly skewed. The first two phenomena are related in that weak density dependence allows large population fluctuations to occur which obscure the inference of relationships and lower the power. All three factors are at play when dealing with sage-grouse in Wyoming. In summary, scaling up the local lek-level impacts produces estimates of population declines similar to those for the individual populations, however there is much uncertainty over the magnitude, if any, of the actual population-level response.
It may be possible to reduce the uncertainty in the population-level effects of OAG development through: the incorporation of additional variables (Ramey, Brown & Blackgoat, 2011), the incorporation of spatial variation in density-dependence (LaMontagne, Irvine & Crone, 2002), by modifying the population dynamic model so that it more closely matches the sage-grouse life history, and/or expanding the analysis to incorporate data from additional populations across a broader section of the species’ range. To enable this, sage-grouse lek count data should be made available to researchers by all states and provinces, which is currently not the case.
Supplemental Information
Fig. S1. Bayesian estimates of the frequency of leks by count of male sage-grouse in a typical year with no oil and gas.
Fig. S2. Bayesian estimates (with 95% credible intervals) of the effect of year on the expected count of male sage-grouse at a typical lek after accounting for the Pacific Decadal Oscillation index and oil and gas.
Fig. S3. Bayesian estimates (with 95% confidence intervals) of the carrying capacity in a typical year with no oil and gas by working group.
Fig. S4. Bayesian estimates (with 95% credible intervals) of the effect of the Pacific Decadal Oscillation index on the expected carrying capacity at a typical working group.
The effect is the percent change in the expected carrying capacity relative to a Pacific Decadal Oscillation index value of 0.
Fig. S5. Bayesian estimates (with 95% credible intervals) of the effect of year on the density the subsequent year after accounting for the Pacific Decadal Oscillation index and oil and gas.
Fig. S6. Bayesian estimates (with 95% credible intervals) of the effect of density on the density the subsequent year with no oil and gas.
Table S1. The relative importance (w_{i}) of spatial scale as a predictor of the count of males sage-grouse at individual leks.
The relative importance is across all models with the areal disturbance due to well pads and the Pacific Decadal Oscillation index both independently lagged one to four years.
Table S2. The relative importance (w_{i}) of the lag in areal disturbance due to well pads as a predictor of the count of males sage-grouse at individual leks.
The relative importance is across all models with a lek distance of 0.8, 1.6, 3.2 and 6.4 km and the Pacific Decadal Oscillation index lagged one to four years.
Table S3. The parameter estimates for the final lek count models with lower and upper 95% confidence/credible intervals.
The estimates are for a lek distance of 3.2 km, areal disturbance due to well pads of one year and Pacific Decadal Oscillation index lag of two years. Model parameters are described in Table 1.
Table S4. The parameter estimates for the final population models with lower and upper 95% confidence/credible intervals.
The estimates are for a lek distance of 3.2 km, areal disturbance due to well pads of one year and Pacific Decadal Oscillation index lag of two years. Model parameters are described in Table 1.