Causal mechanisms for negative impacts of energy development inform management triggers for sagebrush birds

Abstract

Estimated population trends can identify declining species to focus biological conservation, but monitoring may fail to illuminate causes of population change and strategies for reversing declines. Monitoring programs can relate trends with environmental attributes to test causal hypotheses, but typical analytical approaches do not explicitly support causal inference, diluting available data for informing conservation. The U.S. Bureau of Land Management (BLM) extended Integrated Monitoring in Bird Conservation Regions with a quasi-experimental sampling design over a 10-year period (2010–2019) to evaluate the impacts of oil and gas development on sagebrush birds within the Atlantic Rim Natural Gas Field Development Project in southern Wyoming. We analyzed resulting data using a multiscale community occupancy model to estimate trends in species occupancy and richness relevant to management triggers. Additionally, we employed path analysis to evaluate mechanisms underlying observed trends to inform potential management responses. Fine-scale occupancy for sage thrasher (Oreoscoptes montanus) declined within the high-development stratum at a rate sufficient to meet an a priori management trigger established by the BLM. Two additional sagebrush-associated species, Brewer's (Spizella breweri) and sagebrush sparrow (Artemisiospiza nevadensis), exhibited negative development relationships with trend, as did overall species richness, and richness of grassland, sagebrush, and generalist guilds. We identified well pad density and invasive plants associated with energy development as causal factors contributing to these negative development impacts. We demonstrate an analytical approach for both estimating occupancy trends and identifying underlying causes to inform conservation action. Reducing the development footprint, including well pad density and associated invasive plants, could help reduce or limit impacts on birds within this landscape.

INTRODUCTION

Habitat loss, fragmentation, and degradation resulting from anthropogenic development are primary causes of biodiversity declines worldwide (Sala et al., 2000; Wilcove et al., 1998). Ecological monitoring critically informs conservation efforts (Lindenmayer & Likens, 2010) by identifying population trends so conservation can focus on declining species. Nevertheless, monitoring trends often fails to illuminate mechanisms underlying population change, and thus poorly informs conservation action (Nichols & Williams, 2006; Yoccoz et al., 2001). Having identified declining species, conservation practitioners may lack time for studying underlying processes prior to extirpation (Nichols & Williams, 2006).

Monitoring programs can avoid this pitfall by collecting additional data for evaluating mechanisms underlying population change (Barrows et al., 2005; Lindenmayer & Likens, 2010). Mechanisms can help distinguish anthropogenic effects from potentially confounding factors inherent to complex systems. Additionally, knowledge of mechanisms provides a stronger foundation for formulating management intervention (Barrows et al., 2005). Monitoring populations and determining drivers of population change are important for making state-dependent decisions and evaluating conservation success within an adaptive management framework (Lyons et al., 2008).

Along with relevant data, evaluating population trends and their underlying mechanisms requires appropriate analysis. Hierarchical models have vastly improved knowledge production by allowing population estimation while explicitly accounting for the observation process (Kéry & Royle, 2016, 2020). With appropriate models, analysts can relate population parameters with covariates that serve different roles in representing underlying process. When evaluating response to environmental perturbation (e.g., anthropogenic development, management intervention, or natural disturbance), population models could include treatment effects representing levels of perturbation intensity, covariates representing potential causal factors underlying ecological responses, and potentially confounding environmental attributes (Latif, Ivan, et al., 2020; Latif, Truex, et al., 2020). Quasi-experimental designs permit stronger inferences than retrospective analysis of observational data (MacKenzie et al., 2018), but resolving cause and effect relationships may additionally require causal modeling (e.g., path analysis) to adjust design-based estimates for what would be expected under a true randomized and controlled experiment (Arif & MacNeil, 2022).

With growing energy demand and societal interest in curbing reliance on energy imports, oil and gas development (hereafter “energy development”) increasingly impacts sagebrush (Artemisia spp.) ecosystems in North America (Knick et al., 2003). United States federal agencies, particularly the Bureau of Land Management (Department of Interior; BLM), manage most remaining sagebrush-dominated landscapes. Institutional mandates require the BLM to manage public lands in a manner that accommodates both energy development and biological conservation.

Fulfilling this multiuse mandate has become increasingly challenging with growing evidence for negative development impacts on wildlife (Gilbert & Chalfoun, 2011; Ingelfinger & Anderson, 2004; Knick et al., 2003; Mutter et al., 2015). Development of energy infrastructure, including well pads, roads, and pipelines, reduces and degrades habitat (Walker et al., 2020). Roads for transporting resources can fragment habitat and facilitate spread of invasive vegetation (Ingelfinger & Anderson, 2004). Invasive plants alter habitat structure both directly and by altering disturbance regimes (Knick et al., 2003). Additionally, invasive plants can result in changes to food availability and nutritional content (Nelson et al., 2017). Structures associated with development can provide perches for avian predators (Gilbert & Chalfoun, 2011; Ingelfinger & Anderson, 2004). Noise associated with increased traffic volume and extraction infrastructure may interfere with aspects of avian communication (Ingelfinger & Anderson, 2004; McClure et al., 2013).

The BLM administers the Atlantic Rim Natural Gas Field Development Project in southern Wyoming (hereafter “Atlantic Rim”), where they established an avian monitoring program aimed at evaluating energy development impacts on population and diversity trends to inform potential management responses (e.g., Garman, 2018; Kiesecker et al., 2009). The Integrated Monitoring in Bird Conservation Regions (IMBCR) program (Pavlacky et al., 2021, 2017) provided the foundation for strong inference via a quasi-experimental sampling design. Regional IMBCR monitoring provided reference data for distinguishing energy development impacts within the Atlantic Rim from broader population trends. Within this framework, the BLM established management triggers defined by occupancy trends for four sagebrush-associated species of conservation concern (sage thrasher [Oreoscoptes montanus], Brewer's sparrow [Spizella breweri], sagebrush sparrow [Artemisiospiza nevadensis], and green-tailed towhee [Pipilo chlorurus]; hereafter sagebrush-associated species) at two spatial scales, coarse-scale occupancy of 1-km² grid cells and fine-scale occupancy of 125-m radius plots nested within grid cells (BLM, 2007, 2013; Pavlacky et al., 2012). The BLM also sought to understand development impacts on trends for the broader bird community occurring within the sagebrush landscape and guilds representing various vegetative landscape components (i.e., sagebrush, shrubland [includes sagebrush], grassland, habitat generalist, montane, riparian, wetland, and woodland guilds).

We analyzed 10 years of Atlantic Rim monitoring data (2010–2019) using a path analysis coupled with a hierarchical community occupancy model to address these goals. Our objectives were to (1) evaluate whether trends in coarse- and fine-scale occupancy exceeded a priori management triggers, (2) evaluate energy development impacts on community trends, and (3) identify mechanisms underlying apparent development impacts. We hypothesized occupancy trends would be most negative in areas representing the greatest level of development intensity, and we considered support for this hypothesis as evidence for negative development impacts on birds. Moreover, we hypothesized that well pad density, invasive plants, and road density would represent causal factors contributing to negative development impacts.

METHODS Study area

Our study area included BLM-managed public lands interspersed with privately-owned land within the Wyoming portion of the Northern Rockies Bird Conservation Region (BCR 10; NABCI, 2007), reflecting land distribution history (BLM, 2022). The 1093-km² Atlantic Rim portion of the study area was located south of Rawlins, Wyoming, between Highways 789 and 71 and bordered to the south by Highway 70. Energy development within the Atlantic Rim accelerated in 1999–2008, and was largely completed by 2015. Beyond Atlantic Rim, the study area represented remaining BLM-managed lands within BCR 10, encompassing 64,582 km² characterized by more modest intensities of energy development. The study area was primarily shrub–steppe and semidesert shrubland basins dominated by mountain big sagebrush (Artemisia tridentata) and Wyoming big sagebrush (A. tridentata wyomingensis), and Utah juniper (Juniperus osteosperma) woodland in the uplands (1300–2700 m above sea level; Appendix S1).

Sampling design

Our sampling design followed a quasi-experimental design rooted in IMBCR. Primary sampling units were 1-km² grid cells, each containing a 4 × 4 array of evenly-spaced (250 m) point-count plots (hereafter points; Pavlacky et al., 2021, 2017). Our design was an impact level-by-time interaction design (Morrison et al., 2008, p. 250) with two distinct strata reflecting levels of energy development anticipated by the BLM in the Atlantic Rim: a high-development (HD) and a low-development (LD) stratum. Reference data represented IMBCR background monitoring of all BLM-managed lands in BCR 10 and Wyoming representing average conditions in the region (hereafter BG stratum). Within each stratum, we selected primary sampling units following a spatially balanced sampling algorithm (generalized random tessellation stratification; Stevens & Olsen, 2004) implemented in R (Dumelle et al., 2008). Our sample consisted of 1465 points within 99 grid cells surveyed across a 10-year (2010–2019) study period (Table 1). Although not all sampling units were surveyed every year, we followed established methods for selecting spatially balanced subsets of units to ensure the sample was representative of each stratum in each year (Pavlacky et al., 2021, 2017; Stevens & Olsen, 2004).

TABLE 1 Sample sizes for monitoring sagebrush bird communities in relation to varying intensities of energy development in Wyoming, USA, 2010–2019.

Stratum	Grid cells	Cell × year	Survey points	Point × year
High development	22	148	339	1980
Low development	24	136	320	1640
Background	53	384	806	5296
Total	99	668	1465	8916

Note: Survey points were evenly spaced (250 m apart) within grid cells. Unit × year values represent the total number of occasions units of the specified type were surveyed across all years. Background units represent surveys within lands managed by the Bureau of Land Management (U.S. Department of Interior) in Wyoming. Data collection

Bird surveys consisted of surveying points within grid cells for 6 min during the breeding season between 0.5 h before and 5 h after official sunrise. The seasonal timing of bird surveys varied by elevation and largely adhered to the following bounds identified by the IMBCR protocol (Hanni et al., 2011): May 20–July 1 at <2000 m, May 25–July 7 at 2000–2300 m, May 28–July 16 at 2300–2600 m, and June 17–July 19 at >2600 m. Timing of surveys within these elevational bounds was designed to balance logistical considerations (e.g., snow melt, lambing on private lands, opening of mountain roadways) with bird arrival on breeding grounds. Surveyors recorded all birds detected by species, radial distances (in meters) to detected individuals (measured with laser range finders), and time when individuals were first detected within the survey (0–6 min; for details, see RMADC, 2022). For our dataset, we only included detections within 125 m of the point, so sampling plots were effectively 125-m radius nonoverlapping circles centered on points. During 50% of sampling occasions, a subset of the 16 points were surveyed in a grid cell due to failure to obtain landowner permission, inclement weather, or inability to adequately hear bird vocalizations because of ambient noise (mean = 10.7 and min = 3 points per grid cell for those with <16 points surveyed).

Nine environmental metrics derived from remotely sensed data sources served as covariates of species occupancy (Table 2). Development was a categorical covariate distinguishing HD, LD, and BG strata. We derived well pad densities from coordinates provided by the Wyoming Oil and Gas Conservation Commission, and road densities from 2009 spatial data updated with 2019 aerial imagery. The study area largely supports Greater sage-grouse (Centrocercus urophasianus) populations and therefore some surfaces were subject to development restrictions. Nevertheless, strata rankings based on expected development intensity of high, low, and background corresponded with actual mean well pad densities within 900 ha of our sampling units (HD = 10.1, LD = 3.5, BG = 1.2). We interpreted remotely sensed annual herbaceous cover as extent of invasive cheatgrass (Bromus tectorum) cover (Maestas et al., 2020). Additional metrics of annually varying vegetation structure, productivity, topography, and survey timing provided covariates for occupancy and detectability (Table 2; for data sources, see Appendix S2).

TABLE 2 Covariates for community occupancy models estimating trends for sagebrush birds in relation to energy development in Wyoming, USA, 2010–2019.

Variable	Scale (resolution)a	Description	Model parameters
Development	Grid cell	3 strata: (1) high development, (2) low development, and (3) background	Ψ, Δ, ψ, δ, p
Well pad density	Grid cell (900 ha)	No. active well pads and well pads decommissioned after 1990 within a 9 km2 neighborhood	Ψ, Δ
	Grid cell (100 ha)	No. active well pads and well pads decommissioned after 1990 within the 1-km2 grid cell	ψ, δ
Road density	Grid cell (100 ha)	Length of roads (km) per km2	Ψ
	Point (5 ha)		ψ, δb
Cheatgrass	Point (5 ha)	Percent cover of annual herbaceous vegetation (primarily cheatgrass)	ψ, δb
Sagebrush	Point (5 ha)	Percent cover of sagebrush	ψb
Herbaceous cover	Point (5 ha)	Percent cover of perennial herbaceous vegetation	ψb
Pinyon–juniper forest	Grid cell (100 ha)	Percent area classified as pinyon and/or juniper dominated forest	Ψ
Normalized difference vegetation index (NDVI)	Grid cell (100 ha)	Maximum daily breeding season (May 1–July 31) NDVI (indexes vegetation potential each year).	Ψ
Topographic position index	Point (5 ha)	Min topographic position index (TPI) within neighborhood, where TPI is the pixel elevation minus the neighborhood mean	ψb
Day of yearc	Observation	No. days elapsed since January 1 during survey	p
Time since sunrise	Observation	No. minutes elapsed since sunrise when initiating the survey	p
Shrub cover	Point (1 ha)	Ocular estimate of percent shrub cover within 50 m of survey point, with shrubs defined as woody vegetation 0.25–3 m tall	p

Note: Model parameters are those with which covariates were related: occupancy probabilities (Ψ, ψ) and trends (Δ, δ) for 1-km² grid cells and 125-m radius point-count plots (respectively), and detection probability (p). a 900-, 100-, 5-, and 1-ha resolution metrics represented 3-km² grid-cell neighborhoods, 1-km² grid-cell neighborhoods, 125-m radius point-centered circular neighborhoods, and 50-m radius point-centered circular neighborhoods, respectively. b All continuous covariates initially compiled at the point level were averaged (mean) across points in each grid cell to accommodate grid-cell-level analysis. c We allowed a nonlinear relationship of detectability with day of year by specifying a quadratic effect (y ~ x + x²). Data analysis

We analyzed avian species occupancy and richness using a multispecies extension (Dorazio et al., 2011) of a multiscale occupancy model (Mordecai et al., 2011). We describe in detail the analysis model structure and fitting in Appendix S3. In short, the model estimated coarse-scale (grid cell) occupancy probability for each species (Ψ), fine-scale (point) occupancy (ψ; conditional on coarse-scale occupancy), occupancy trends at coarse and fine scales (Δ and δ, respectively), and detectability of the species (p). We modeled these primary parameters as logit-linear functions of covariates (Table 2). Occupancy at both scales and detectability also varied annually over and above variation explained by estimated trends and covariate relationships. Considering the short (6-min) timeframe for estimation, point occupancy likely reflected variability in local abundance for species with territories ⪆5 ha, whereas grid-cell occupancy quantified coarser scale species distributions (Latif et al., 2016; Steenweg et al., 2018).

Our model specified species-specific parameters (Ψ, ψ, Δ, δ, p, and covariate relationships) as random variables governed by community-level hyper-parameters. This structure allowed information sharing across species, improved species-specific parameter precision, and allowed us to derive species richness by summing occupancy probabilities across species (Dorazio et al., 2011; Kéry & Royle, 2016). We excluded raptors, grouse, cranes, and water birds not readily detectable with our survey methods, and species that did not regularly breed in our study area. We augmented data to represent all potential breeders to fully correct species richness estimates for imperfect detection (Dorazio et al., 2011). Our species list included 114 species detected during our surveys and 25 additional species detected in broader regional surveys (Latif et al., 2023). Thus, species richness estimates represented the expected number of relatively small-bodied, territorial, breeding landbird species present within a grid cell (100 ha) or point (5 ha). In addition to overall richness, we derived richness for species guilds defined by ecological association: habitat generalist guild (18 species), grassland guild (18 species), montane guild (34 species), riparian guild (18 species), sagebrush guild (4 sagebrush-associated species), shrubland guild (15 species, including the 4 sagebrush-associated species), wetland guild (5 species), and woodland guild (31 species; Latif et al., 2023). We considered the sagebrush-associated guild as a subset of the shrubland guild. Ninety percent of shrubland species detections were of sagebrush-associated species, however, and covariate relationships and trends were extremely similar for sagebrush versus overall shrubland guilds. We therefore only report richness results for the sagebrush guild, the guild of primary interest to management.

To assess management triggers, we compared sagebrush-associated species occupancy trends across development strata. Management triggers identified by the BLM concerned four sagebrush-associated species: Brewer's sparrow (S. breweri), sagebrush sparrow (A. nevadensis), sage thrasher (O. montanus), and green-tailed towhee (P. chlorurus). Trigger criteria were a 25% decline in coarse-scale occupancy or a 10% decline in fine-scale occupancy within the HD stratum over the 10-year monitoring period (2010–2019) that was also ≥2% more negative than the LD or BG trend at the same scale (BLM, 2007, 2013). Per BLM specification, trends needed to meet or exceed these criteria with 90% probability to trigger management. Accordingly, we considered a criterion met if the 80% credible interval for the relevant parameter (trend or difference in trends) exceeded the trigger value. Hereafter, reference strata refer to either LD or BG strata, and we describe cases where we found ≥90% probability of a HD trend being different from trends in either reference stratum as a “HD trend relationship.” To additionally inform avian species and community responses to development, we report HD trends for all other species and for overall community and guild richness.

Well pad density, road density, and cheatgrass represented mechanistic covariates capable of explaining HD trend relationships (hereafter mechanistic covariates). Additional covariates controlled for potentially confounding environmental variation (hereafter control covariates). We related control and mechanistic covariates with occupancy (Ψ and ψ), but only mechanistic covariates with trend (Δ and δ). Road density relationships with coarse-scale trend (Δ) were invariant in the initial analysis, indicating insufficient information in the data for differentiating species-specific relationships with this covariate. We therefore treated coarse-scale road density (100 ha) as a control covariate (i.e., related with Ψ but not Δ).

By estimating trend relationships with mechanistic covariates, we were able to consider potential mechanisms underlying energy development effects using path analysis (Figure 1; Clough, 2012). Specifically, we calculated percent contributions of mechanistic pathways to negative HD trend relationships (Appendix S3). We only evaluated mechanisms where strata-specific differences in mechanistic covariate values were consistent with the hypothesized environmental impacts of development (i.e., negative due to greater well pad density, road density, and cheatgrass in the HD compared with reference strata), and where the mechanistic covariate relationship with trend was negative.

FIGURE 1. Model relationships to assess energy development effects on sagebrush bird species occupancy and richness trends, and underlying mechanisms. Solid lines represent trend relationships, whereas dotted lines represent relationships with occupancy. For covariate descriptions, see Table 2.

RESULTS

The five most frequently detected species were Brewer's sparrow, horned lark (Eremophilia alpestris), vesper sparrow (Pooecetes gramineus), green-tailed towhee, and sagebrush sparrow (Latif et al., 2023). As expected, well pad densities were greatest in the HD stratum, intermediate in the LD stratum, and lowest in the BG stratum (Appendix S4). Road densities and cheatgrass cover were similar between HD and LD strata, but both were greater than in the BG stratum (Appendix S4). Control covariates differed among strata and years (Appendix S4), and accounted for variability in occupancy for most species (Appendix S5). Posterior median detectability estimates ranged 0.36–0.99, and 42 species exhibited statistically supported covariate relationships with detectability (Latif et al., 2023).

Species occupancy trends

Twenty-seven species exhibited HD trend relationships (Figure 2); only 2 species (vesper sparrow and green-tailed towhee) exhibited positive HD trend relationships, whereas 24 species exhibited negative HD trend relationships consistent with hypothesized development effects. These 24 species represented all guilds except wetland species, with grassland and generalist guilds represented most frequently (6 and 7 species, respectively). HD trend relationships were most apparent for fine-scale occupancy. Only seven species exhibited coarse-scale HD trend relationships, all of which were negative (Figure 2).

FIGURE 2. Occupancy trends (posterior median estimates and 80% credible intervals) over 2010–2019 for 27 species with supported differences between high-development and reference (low-development or background) strata trends. Log-scale trends are displayed, but with x-axis labels indicating percentage of change in occupancy. Single asterisks indicate reference strata trends that differed from high-development trends with 90% confidence. Double asterisks indicate reference strata trends that exceeded the high-development trend by 2% with 90% confidence (required for management triggers). Fine-scale occupancy trends describe changes within occupied grid cells. Species codes are defined in Latif et al. (2023) and codes in parentheses indicate guild membership (habitat generalist [Gn], grassland [Gr], montane [M], riparian [R], sagebrush-associated [Sg], shrubland but not necessarily sagebrush-associated [Sh], or woodland [Wd]). Species are listed in taxonomic order.

One sagebrush-associated species, sage thrasher, exhibited a negative HD trend relationship that met performance criteria defining a management trigger (Figures 2 and 3). We estimated a 48% decline (posterior median) in fine-scale occupancy for sage thrasher from 2010 to 2019 in the HD stratum, whereas we estimated a 183% increase and little change (5% increase) in the LD and BG strata, respectively, over the same period. Two other sagebrush-associated species, Brewer's and sagebrush sparrow, also exhibited negative fine-scale HD trend relationships, although not sufficient to meet management triggers (Figures 2 and 3). Fine-scale occupancy for the fourth sagebrush-associated species, green-tailed towhee, increased in both HD (44%) and LD (63%) strata while remaining relatively constant (6% increase) in the BG stratum. The data provided less support for coarse-scale compared with fine-scale trends for sagebrush-associated species. Moreover, a negative coarse-scale HD trend relationship for Brewer's sparrow likely reflected bounds on occupancy—occupancy probability was near one in the HD stratum with no room to increase, unlike in the BG stratum (Figure 3).

FIGURE 3. Predicted occupancy (posterior median and 80% credible intervals) for sagebrush-associated bird species over the study period. Lines and ribbons show predicted mean occupancy across years at model-estimated covariate values for each stratum. Points and error bars also include annual deviations around the mean. Fine-scale values describe occupancy of points within occupied grid cells.

Seven grassland species exhibited HD trend relationships, six of which were negative, with only the vesper sparrow relationship being positive (Figure 2). The positive vesper sparrow HD trend relationship was of a much smaller magnitude than the six negative HD trend relationships observed for the other species. In particular, horned lark and western meadowlark, the first and third most frequently detected grassland species, respectively (Latif et al., 2023), exhibited 55% and 59% declines in fine-scale occupancy within the HD stratum over the study period. Western meadowlark (Sturna neglecta) also exhibited a 28% HD decline in coarse-scale occupancy compared with 19% and 22% LD and BG increases, respectively.

Seven generalist species exhibited negative HD trend relationships (Figure 2), of which American robin (Turdus migratorius) and mourning dove (Zenaida macroura) were the most frequently detected (Latif et al., 2023). All of these generalists except barn swallow (Hirundo rustica) exhibited differences in trend among strata that implied nonlinear relationships with development intensity, that is, LD trend was more positive than both HD and BG trends (Figure 2). Coarse-scale occupancy trends did not relate nonlinearly with development, however, which is why we did not see a similar pattern for fine-scale generalist guild richness (reviewed below).

Species richness trends

We found clear evidence for a negative HD trend relationship for overall species richness at coarse and fine spatial scales (Figure 4; Appendix S6). Fine-scale community-wide richness decreased by 21% in the HD stratum while increasing by 38% and 20% in LD and BG strata, respectively. Reflecting species occupancy trends (reviewed above), sagebrush, grassland, and generalist guild richness exhibited clear negative HD trend relationships. Grassland guild richness exhibited the strongest negative HD trend relationships of any guild. Grassland guild richness declined by 46% in the HD stratum while increasing by 31% and 10% in LD and BG strata, respectively (Figure 4; Appendix S6). Reflecting species-level relationships, fine-scale montane and riparian guild richness exhibited negative HD trend relationships (Figure 4; Appendix S6), but these guilds represented relatively small portions of the overall community, and montane guild richness also declined in the BG stratum.

FIGURE 4. Predicted species richness (posterior median and 80% credible intervals) for the overall community and community components over the study period. Lines and ribbons show predicted mean richness across years at model-estimated covariate values for each stratum. Points and error bars also include annual deviations around the mean. Fine-scale richness represents the sum of unconditional species occupancy probabilities (Ψ × ψ), and thus is not conditional on coarse-scale richness.

Mechanistic contributions to trends

We found statistically supported fine-scale trend relationships with at least one of the three mechanistic covariates for 30 species (Figure 5). Trends for 23 species related negatively with well pad density, whereas we only found a supported positive trend relationship for one species (vesper sparrow). Eleven species exhibited fine-scale trend relationships with road density, eight of which were positive. Five species exhibited fine-scale trend relationships with cheatgrass, two of which were negative and three positive. We found no supported coarse-scale trend relationships with well pad density.

FIGURE 5. Mechanistic covariate effects (posterior median and 80% credible intervals) on grid-cell (Δ) and fine-scale (δ) occupancy trends for 30 species with at least one statistically supported positive (red) or negative (blue) effect. Species codes are defined by Latif et al. (2023), and codes in parentheses indicate guild membership (habitat generalist [Gn], grassland [Gr], montane [M], riparian [R], sagebrush-associated [Sg], shrubland but not necessarily sagebrush-associated [Sh], or woodland [Wd]). Species are listed in taxonomic order.

We found evidence that mechanistic covariates represented causal factors contributing to negative HD trend relationships for 10 species (Table 3). Well pad density contributed to negative HD trend relationships relative to LD trends for eight species, and relative to BG trends for four species. Cheatgrass contributed to a negative HD trend relationship for horned lark. No species exhibited a negative HD trend relationship and a negative trend relationship with road density, so we were unable to confirm road density as a causal factor for negative HD trend relationships. In many cases, estimated percent contributions of mechanisms were small (posterior median <10% for 6 of 13 contributions) or uncertain even if credible intervals excluded zero. Because we found no statistically supported coarse-scale trend relationships with well pad density (Figure 5), we could not corroborate well pad density as contributing to negative coarse-scale HD trend relationships.

TABLE 3 Statistically supported mechanisms for development effects on fine-scale species occupancy and richness trends.

Species/guild	Mechanistic covariate	Contributions
		Low-development reference	Background reference
Common nighthawk	Well pad density (100 ha)	4.51 (0.32, 25.46)	…
Common raven	Well pad density (100 ha)	4.91 (0.02, 28.25)	…
Horned lark	Cheatgrass	…	6.74 (3.26, 11.71)
House wren	Well pad density (100 ha)	3.08 (0.37, 13.95)	…
Sage thrasher	Well pad density (100 ha)	1.18 (0.11, 2.96)	10.07 (0.84, 25.83)
American robin	Well pad density (100 ha)	14.52 (3.70, 48.94)	…
Brewer's sparrow	Well pad density (100 ha)	…	13.73 (4.90, 23.76)
Sagebrush sparrow	Well pad density (100 ha)	15.48 (5.23, 50.07)	…
Western meadowlark	Well pad density (100 ha)	2.10 (0.15, 6.37)	25.71 (2.71, 66.82)
Brewer's blackbird	Well pad density (100 ha)	19.41 (4.69, 68.97)	51.85 (9.24, 231.11)
All species	Well pad density (100 ha)	8.44 (4.09, 14.7)	24.30 (14.08, 39.46)
Sagebrush	Well pad density (100 ha)	13.02 (6.39, 24.44)	20.54 (11.57, 34.54)
Generalist	Well pad density (100 ha)	11.26 (3.93, 28.96)	…
Grassland	Cheatgrass	…	5.71 (2.04, 10.81)
Montane	Well pad density (100 ha)	7.78 (0.24, 39.35)	…
Riparian	Well pad density (100 ha)	8.15 (0.24, 37.84)	…

Note: Displayed values are percent contributions (posterior median and 80% credible intervals) of mechanistic pathways to the difference in trend between high-development and reference strata where 80% credible interval for high-development trend minus reference trend fell entirely below zero.

We corroborated the mechanisms for negative fine-scale HD trend relationships for the overall community and five of seven guilds (i.e., all except wetland and woodland species; Table 3, Figure 6). We corroborated well pad density as contributing to negative HD trend relationships for overall richness and four guilds. Additionally, although not clearly corroborated as a mechanism for the grassland guild, well pad density did contribute to negative HD trend relationships for three grassland species (common nighthawk [Chordeiles minor], Brewer's blackbird [Euphagus cyanocephalus], and western meadowlark). Percent contributions of well pad density to negative HD richness trend relationships ranged 8%–24% (Table 3). Cheatgrass contributed somewhat (6%) to a negative HD trend relationship relative to the BG stratum for grassland guild richness, likely reflecting the mechanism for horned lark (see above). Mirroring species-level results, we were unable to corroborate well pad density as a mechanism for negative coarse-scale HD trend relationships.

FIGURE 6. Predicted species richness (posterior median and 80% credible intervals) for the overall community and community components over the study period at alternate levels of mechanistic covariates. Development and mechanistic covariates other than those represented in each panel were set to levels representing the high-development stratum, and control covariates were set to their study-wide mean values. Fine-scale richness represents the sum of unconditional species occupancy probabilities (Ψ × ψ), and thus is not conditional on coarse-scale richness.

DISCUSSION

Negative impacts of Atlantic Rim energy development on fine-scale sage thrasher occupancy surpassed mitigation triggers developed by the BLM (BLM, 2007, 2013). We also found substantial evidence corroborating hypothesized negative development impacts for other sagebrush-associated species (Brewer's and sagebrush sparrow) and the broader bird community. Moreover, we corroborated well pad density as a causal factor contributing to these apparent impacts. Evidence for mechanisms further corroborates development impacts, while also suggesting foci for potential management.

Considering our results, well pad densities represent the most obvious potential focus for management aimed at limiting or offsetting energy development impacts. Supported well pad density relationships with fine-scale occupancy trends were almost universally negative, driving negative relationships with species richness trends. Managers could offset negative impacts by limiting or excluding well pads where doing so would ensure increases or no net loss of critical habitat (Kiesecker et al., 2009). Techniques such as horizontal drilling technologies (Garman, 2018) could additionally limit the negative impacts of energy development on birds while facilitating continued resource extraction.

Predominantly negative relationships of species distributions with well pad density are broadly consistent across studies, although reported relationships for individual species vary (compared with Gilbert & Chalfoun, 2011; Mutter et al., 2015). Unlike here, Gilbert and Chalfoun (2011) reported negative relationships with well pads for vesper sparrow, and no supported effects for sage thrasher. They suggest strong site fidelity may obscure negative development impacts on sage thrasher populations (see also Mutter et al., 2015), and accordingly predicted population declines in developed areas over longer temporal periods than represented in their study. Our quasi-experimental sampling design implemented over a longer study period allowed us to corroborate their prediction.

Steep population declines suffered by grassland birds worldwide make this guild a high conservation priority (Peterjohn & Sauer, 1999; Rosenberg et al., 2019). Unfortunately, mechanisms underlying development effects for grassland species were less clear than other guilds in our study. Trend relationships with well pad density for the three most common members (horned lark, vesper sparrow, and western meadowlark) were relatively weak. Cheatgrass contributed somewhat to diminished trends for horned lark. Horned lark typically occur in areas with sparse vegetation and extensive bare ground (Beason, 2020), whereas cheatgrass reduces bare ground where it has invaded extensively (Klemmedson & Smith, 1964).

Positive species occupancy and trend relationships with roads estimated here (see also Appendix S7) appear to contradict patterns reported elsewhere in sagebrush and other systems (Fahrig & Rytwinski, 2009; Gilbert & Chalfoun, 2011; Ingelfinger & Anderson, 2004; Mutter et al., 2015). We suggest well pads represent a localized process of habitat loss, whereas roads represent a linear process of habitat fragmentation. By evaluating the causal mechanisms of landscape change (Didham et al., 2012), we discovered that habitat loss from increasing well pad densities was more detrimental to avian diversity than habitat fragmentation from increasing road densities. Nevertheless, habitat fragmentation may interact with habitat loss such that negative fragmentation effects may not arise without sufficient habitat loss (Didham et al., 2012). Moreover, redistribution following habitat loss may concentrate populations in remnant areas (Ewers & Didham, 2006), some of which will inevitably include roads. Short-term crowding effects in recently fragmented landscapes may give way to an extinction debt, with trends eventually turning negative in regions with high road densities (Ewers & Didham, 2006).

Negative development effects on habitat generalists here contrast with the positive effects on resource generalists reported in other systems (e.g., Clavel et al., 2011). Negative impacts for habitat generalists along with specialists suggests negative energy development impacts extend beyond specialists to species typically more tolerant of other forms of human development (e.g., agricultural or urban), which provide ornamental fruit-bearing trees and shrubs and other supplemental food resources (Evans et al., 2009). Non-native grasses planted for reclamation of well pad sites can promote deer mice (Peromyscus maniculatus), which may elevate nest predation pressure for birds in general (Sanders & Chalfoun, 2018). The BLM requires native seed mixes for reclamation, which occurs both after well activation in disturbed areas not required for production and after well retirement, but non-native grasses can become established, especially on private lands.

Despite clear evidence for mechanisms, mechanistic pathways considered here often accounted for only a small percentage of negative HD trend relationships. Nonlinear and interactive effects of habitat loss, fragmentation, and population redistribution may account for a greater percentage of apparent relationships (see above). Additionally, areal extent of disturbance (i.e., complement of loss of native vegetation) may be more relevant for understanding response to development by sagebrush-associated species, but may not scale linearly with well pad or road densities (Hethcoat & Chalfoun, 2015). We assumed that the overall development footprint would largely track well pad and road densities, but other landscape features could also contribute and may be worth considering further (Walker et al., 2020). Moreover, although the mechanistic factors considered here represented a substantive range of possible development effect drivers, we could only consider those for which data were readily available. Additional factors not considered include levels of traffic, reclamation, and sound, light, and air pollution. Causal models trained on quasi-experimental data in this study could provide the basis for bird conservation objectives for adaptive management (Lyons et al., 2008) or mitigation of oil and gas development in semiarid rangelands (Kiesecker et al., 2009).

Study strengths

Our quasi-experimental sampling design and representation of bird distributions at multiple scales across a broad spatial and temporal extent provided a nuanced view of development effects not necessarily available to others. Given the logic of our design, we expected larger differences in trend between the HD and BG strata, and less difference between HD and LD. In contrast, we found many cases where species exhibited stronger positive LD compared with BG trends (e.g., sage thrasher), pointing to possible nonlinear effects of development (e.g., see discussion of habitat fragmentation and loss above). Evaluation of occupancy trends, rather than just bird distributions in space, also strengthened our inferences. Where high-quality habitat and desirable sites for development coincide, habitat quality may confound inference of development effects based solely on spatial analysis. By comparing population change over time at varying development intensity, trend analyses can provide inference analogous to before–after, control–impact studies (Popescu et al., 2012).

Analysis within a Bayesian framework allowed us to derive parameter estimates directly applicable to management triggers, which specified percent changes in occupancy probability for sagebrush-associated species. Because trends in occupancy probability are constrained when initial occupancy is near zero or one (e.g., the coarse-scale HD Brewer's sparrow trend), MacKenzie et al. (2018, p. 361) recommend evaluating trends in odds occupancy instead. Although we modeled trends on the log-odds scale, we derived posterior trend estimates on the probability scale, which are directly relevant to BLM-established management triggers but also subject to computational pitfalls described by MacKenzie et al. (2018, p. 361). We verified that our inferences would have been similar had we evaluated trends in odds occupancy (Appendix S7). Nevertheless, we recommend managers consider developing trend-based triggers on odds occupancy rather than probability of occupancy.

Broader implications

Ecologists advocate explicitly designing monitoring to evaluate how we think systems operate (Barrows et al., 2005; Lindenmayer & Likens, 2010; Nichols & Williams, 2006). Along with population surveys, remote sensing products increasingly track environmental factors expected to drive population change, providing data critical for evaluating hypothesized causal relationships. We demonstrate the use of path analysis combined with a quasi-experimental sampling design (Arif & MacNeil, 2022) to leverage such data toward evaluating our conception of an ecological system relevant to land management.

AUTHOR CONTRIBUTIONS

Quresh S. Latif developed and implemented data analysis, drafted the manuscript, and incorporated feedback from co-authors. Nicholas J. Van Lanen worked with BLM partners to develop the sampling design and funding agreement for this work. Nicholas J. Van Lanen supervised field crews and managed data curation. Eric J. Chabot processed and curated spatial data, generated the study area map (Appendix S1), and provided GIS support. David C. Pavlacky Jr. helped develop a prototype of the community occupancy portion of our analysis model, and helped develop research questions and inference for the study. Nicholas J. Van Lanen, David C. Pavlacky Jr., and Eric J. Chabot provided editorial support and feedback during manuscript development.

ACKNOWLEDGMENTS

This work was supported by the Bureau of Land Management (United States Department of Interior) and partners of the Integrated Monitoring in Bird Conservation Regions program. We thank F. Blomquist and H. Cline for their commitment and guidance toward implementing this study. Bird Conservancy of the Rockies administers the IMBCR program and implemented Atlantic Rim bird surveys. We also thank Bird Conservancy of the Rockies and Colorado State University for their support during manuscript development. We thank C. White and M. McLaren for assistance with project logistics and supervision. This manuscript benefitted from suggestions of two anonymous reviewers.

CONFLICT OF INTEREST STATEMENT

The authors declare no conflicts of interest.

DATA AVAILABILITY STATEMENT

All scripts and data (Latif et al., 2023) are available from Dryad: https://doi.org/10.5061/dryad.tb2rbp04r. Remotely sensed data used in our analyses are accessible through the links found in Appendix S2.