Shale gas development continues to outpace the implementation of best management practices for wildlife affected by development. We examined demographic responses of the Louisiana Waterthrush (Parkesia motacilla) to shale gas development during 2009–2011 and 2013–2015 in a predominantly forested landscape in West Virginia, USA. Forest cover across the study area decreased from 95% in 2008 to 91% in 2015, while the area affected by shale gas development increased from 0.4% to 3.9%. We quantified nest survival and productivity, a source–sink threshold, riparian habitat quality, territory density, and territory length by monitoring 58.1 km of forested headwater streams (n = 14 streams). Across years, we saw annual variability in nest survival, with a general declining trend over time. Of 11 a priori models tested to explain nest survival (n = 280 nests), 4 models that included temporal, habitat, and shale gas covariates were supported, and 2 of these models accounted for most of the variation in daily nest survival rate. After accounting for temporal effects (rainfall, nest age, and time within season), shale gas development had negative effects on nest survival. Population-level nest productivity declined and individual productivity was lower in areas disturbed by shale gas development than in undisturbed areas, and a source–sink threshold suggested that disturbed areas were more at risk of being sink habitat. Riparian habitat quality scores, as measured by a U.S. Environmental Protection Agency index and a waterthrush-specific habitat suitability index, differed by year and were negatively related to the amount of each territory disturbed by shale gas development. Territory density was not related to the amount of shale gas disturbance, but decreased over time as territory lengths increased. Overall, our results suggest a decline in waterthrush site quality as shale gas development increases, despite relatively small site-wide forest loss.
INTRODUCTION
The United States has 20 shale formations that contain unconventional natural gas deposits which are being developed using new horizontal drilling and hydraulic fracturing techniques (USDOE 2009, USEIA 2011). The central Appalachian region is experiencing the most rapid growth in unconventional drilling activity and development (hereafter, ‘shale gas development'; MCOR 2016) as the underlying Marcellus-Utica shale is the most expansive shale basin with the most potentially recoverable gas (USDOE 2009). In 2015, the state of West Virginia alone had 2,966 (20%) of 14,619 producing Marcellus wells (WVGES 2015, MCOR 2016), with upward of 40,000 wells projected to be operational in West Virginia by 2030 (USDOE 2010). Within the Appalachian region, nearly 75% of close to 3.1 million forested hectares are at high potential risk from energy development, primarily shale gas (Dunscomb et al. 2014). West Virginia is 1 of 2 states with the highest probability for development (21%; Dunscomb et al. 2014), and species with specialized habitat needs that overlap these forested areas will be the most vulnerable to multiple ecological risks (Brittingham et al. 2014).
Recent studies have examined how shale gas development may positively and negatively influence Appalachian songbird communities (Barton et al. 2016, Farwell et al. 2016), how predators and their aquatic prey respond (Wood et al. 2016), and how shale gas development may be associated with bioaccumulated contaminants (Latta et al. 2015). However, there has not been a mechanistic assessment of the interplay between shale gas disturbance and songbird demographic responses (Northrup and Wittemyer 2013). Baseline demographic data is needed in both terrestrial and aquatic habitats of the Appalachian shale gas basin to detect and understand changes as they begin to occur (Brittingham et al. 2014).
The proximity of shale gas development to water resources is of particular concern due to the potential for sedimentation runoff, reduced streamflow, and contamination of surface waters (Entrekin et al. 2011). Therefore, biological communities and organisms that use water resources downstream of shale gas development are at increased risk from shale gas activities near surface waters (Latta et al. 2015). The Louisiana Waterthrush (Parkesia motacilla; hereafter, ‘waterthrush') is a habitat specialist and species of conservation concern (USFWS 2008) that breeds along forested headwater streams and feeds primarily on benthic macroinvertebrates (Mattsson et al. 2009). Waterthrushes are well-established biological indicators of aquatic stream integrity (O'Connell et al. 2000, Mulvihill et al. 2008) and reach some of their highest abundances in the Marcellus shale region (Sauer et al. 2017). As such, where shale gas development is creating increasingly intense core forest disturbance in areas where headwater streams and this species co-occur, the waterthrush is an ideal organism to assess the potential demographic consequences of shale gas development (Drohan et al. 2012, Farwell et al. 2016). Because the waterthrush is a bioindicator species, we hypothesized that we would detect an inverse relationship between the amount of shale gas disturbance and the species' demographic metrics, despite the species' ability to compensate for resource loss (Mulvihill et al. 2008, Wood et al. 2016). Identifying how variability in habitat quality contributes to population surpluses (source habitat; Pulliam 1988) or deficits (sink habitat) is key to long-term conservation planning in landscapes undergoing development (Kirol et al. 2015). If source–sink dynamics were to exist in our local population, we hypothesized that productivity would differ between areas disturbed and undisturbed by shale gas development.
We examined how increased shale gas development influenced waterthrush demography during 2009–2011 and 2013–2015. We quantified waterthrush nest survival and productivity, an overall source–sink threshold, riparian habitat quality, and territory density and length. In addition to identifying demographic responses, our results should inform gas well siting guidelines for shale gas development to minimize risks to ecological resources.
METHODS
Study Area
We studied waterthrush demography along 58.1 km of first- and second-order forested headwater stream tributaries (n = 14) in the Lewis Wetzel Wildlife Management Area (LWWMA) in northwestern West Virginia, USA (39.490216°N, 80.650713°W; Figure 1), an area that supports the highest waterthrush densities within the central Appalachians (Sauer et al. 2017). The study area lies within the Permian Hills subdivision of the Western Allegheny Plateau Ecoregion, an area of deeply dissected topography and relatively continuous Appalachian oak and mixed-mesophytic forest (Woods et al. 1999) at elevations of 221–480 m that overlays the Marcellus-Utica shale region. In 2008, 95% of the LWWMA was forested and 0.4% of the land cover was under shale gas development; the first shale gas well development began in 2007 (Farwell et al. 2016). Shale gas development in our study area and within the surrounding region has since increased rapidly (WVGES 2015). In 2015, 91% of the LWWMA was forested and 3.9% was under shale gas development, with 83% of the shale gas development having resulted in direct forest loss (Farwell et al. 2016).
FIGURE 1.
Our study area included 14 forested first- and second-order headwater streams in the Lewis Wetzel Wildlife Management Area (WMA). It lies within the Marcellus-Utica shale basin, is within the core Louisiana Waterthrush (LOWA) range based on Breeding Bird Survey (BBS) data, and falls in the Appalachian Mountains Bird Conservation Region (BCR) in the eastern United States. LOWA relative abundance represents average Louisiana Waterthrush counts from nearby BBS routes (detections per 50 stops).
During our study, gas well development activities included the building of conventional (shallower formations) and Marcellus well pads, forest clearance for yet unbuilt well pads, the expansion of existing road and pipeline infrastructure, and the construction of new infrastructure. Early in the study (2009–2010), the majority of Marcellus wells and their water holding ponds were located along the main stem of Buffalo Run, into which the majority of our headwater study streams empty. Thus, although a few Marcellus well pads were located along our study streams, they tended to primarily affect the lower portions of the study streams. Between the 2010 and 2011 breeding seasons, shale gas development activities accelerated across the study area and began to increase, especially on ridgetops. As a result of ridgetop activity, the whole downstream network of some streams became disturbed by sedimentation and surface runoff which continued for the remainder of the study. The LWWMA experienced a 1.5% increase in the area affected by timber harvests in 2010–2011 (Farwell et al. 2016), but these and preexisting harvests did not result in complete forest canopy loss (Sheehan et al. 2014) and typically did not intersect or influence streams where we monitored waterthrushes. Shale gas development peaked in 2011, but started to abate site-wide in 2013 and in general became concentrated around specific streams and ridgetops. Clearing for additional new well pads occurred late (June–July) in the 2013 breeding season, with well pad construction and the redrilling of an existing well pad completed in 2014. There was no new shale gas development or activity in the 2015 breeding season. There were no ‘control' streams in our study, given that the majority of forest loss and fragmentation in the surrounding landscape resulted from recent shale gas activities (Farwell et al. 2016), but shale gas disturbance was concentrated on some streams more than others, as indicated by the large annual range in the percentage of stream disturbed (range: 0–67%; Table 1).
TABLE 1.
Annual percentage of each study stream's length that was disturbed by shale gas development or activity (StreamGas, SG; see Table 2) or potentially affected by runoff (StreamRunoff, SR) from shale gas well pads (not including pipeline or road disturbance) in 2009–2011 and 2013–2015 in an area undergoing shale gas development in the eastern U.S. The last new well pad construction occurred near study streams in 2014 (see Figure 1 for study stream locations).
TABLE 2.
Variables used in analyses evaluating the demographic response of Louisiana Waterthrushes to shale gas development. Nest survival is daily survival rate (DSR) over a 29-day nesting period (DSR29).
Mapping and Quantifying Disturbances
Employing a Geographic Information System (GIS), we used a sequence of leaf-on and leaf-off aerial photographs from the National Agriculture Imagery Program (NAIP) for 2011 and 2014, satellite Quickbird imagery for 2009, and extensive annual ground-truthing to manually digitize areas of forest canopy disturbance within the study area for each year of our study. All forest canopy disturbances were classified as related to shale gas development (i.e. well pads and associated roads and pipeline infrastructure) or as unrelated or preexisting development (i.e. forest roads, recent even-aged timber harvests, and various types of existing clearings) using FRAGSTATS 4 (McGarigal et al. 2012). We determined the percentage of forest canopy disturbance from shale gas development (GasFCD) and unrelated development (OtherFCD) within a 100-m radius of each waterthrush nest for use as habitat covariates in nest survival models (Table 2). We used a 100-m radius because forest edges may negatively affect the reproductive success of ground-nesting species at this scale (Flaspohler et al. 2001). We classified a few conventional impacts (i.e. stream-side vertical pump jacks) as related to shale gas development because their well pads were managed in conjunction with nearby shale gas infrastructure and because their targeted formation, even though they remained shallow after development, was listed as Marcellus (WVGES 2015). Gas well records (WVDEP 2015) were used to verify target shale formations, drilling status, and start dates for all gas well disturbances.
The surface (i.e. 3D) length of each study stream (average length: 4.10 ± 0.54 SE km; range: 0.95–7.40 km) was calculated in GIS using a 3 m resolution digital elevation model and defined to have a drainage basin of 9.0 ha (24-km scale or higher resolution) to delineate the uppermost headwater reaches. To describe and model waterthrush demography and riparian habitat quality as a function of disturbance from shale gas development, we created 4 continuous variables and 1 binary variable based on disturbance categories at the stream, territory, and nest scale. The first (StreamGas) described mostly localized streamside disturbance indicative of the presence of any shale gas infrastructure or activity (Figure 2A, Table 2). A section of stream was considered disturbed when well pads, infrastructure, or frequent vehicular activity (Figure 2E) were within 60 m of the stream centerline, the typical extent of waterthrush streamside use (Mattsson and Cooper 2009). When a stream had visually observable sedimentation from shale gas development (Figure 2F), we classified the entire stream network downstream of the sedimentation beginning point as disturbed. The condition of streams was frequently and extensively ground-truthed each season, so there were no stream reaches where sedimentation events were likely to have been missed.
FIGURE 2.
Examples of quantifying shale gas disturbances at the scale of the (A, B) study stream and (C, D) Louisiana Waterthrush territory, and (E, F) examples of types of disturbance in an area undergoing shale gas development in the eastern U.S. (E) A section of stream considered disturbed from frequent vehicular activity; and (F) Sedimentation that we traced upstream to its point of origin at an active shale gas well pad. See Table 2 for variable notation.
We created a second shale gas disturbance category (StreamRunoff) that focused solely on potential runoff into streams from shale gas contaminants (Figure 2B, Table 2). A stream was considered disturbed at and below a well pad or retaining pond (similarly to Latta et al. 2015), resulting in the whole downstream network being classified as at risk of surface pollution based on elevational maps and ground-truthing. The StreamRunoff category did not include pipeline or road disturbance and was a broader, distance-independent disturbance category describing potential water pollution. For each year of the study, we calculated the proportion of each stream disturbed for both the StreamGas and the StreamRunoff disturbance categories.
We calculated the proportion of each waterthrush territory (a 60-m buffer around each territory vector; Mattsson and Cooper 2009) that was disturbed by StreamGas and called this metric TerrGas (Figure 2D, Table 2). The proportion of each territory disturbed by StreamRunoff was termed TerrRunoff (Figure 2C, Table 2). We classified each waterthrush nest location as undisturbed or disturbed by StreamGas within 60 m of the nest and called this variable NestGas.
Territory Density and Length
We mapped waterthrush territories in 2009–2011 and 2013–2015 along 14 streams with varying amounts of disturbance from shale gas development and contaminant runoff (StreamGas and StreamRunoff; Table 1) to determine annual territory density and length. Territories were mapped for 11 streams in all years. Hiles and Huss Pen were mapped in 2014–2015, and Carpenters was mapped in 2009 and 2014–2015 due to restricted access. Territorial waterthrushes were target-netted and banded with an aluminum U.S. Geological Survey leg band and a unique combination of plastic color bands to allow identification of individuals. Waterthrush territories are typically described by length (m) rather than area, given that birds form linear territories along a stream reach (Mulvihill et al. 2008); therefore, all analyses were based on territory length. Waterthrush territory density is defined as the number of individual territories per km of stream monitored (Hallworth et al. 2011). Waterthrush territories were typically delineated from April 1 to June 29 each year. Standardized territory mapping (Robbins 1970, Bibby et al. 1992) was conducted in 2009 (8 visits per stream) and 2010 (5 visits per stream). In 2011, we switched to a more opportunistic approach that allowed us to map waterthrush locations and behaviors during nest searching of each stream (≥5 visits per stream). During 2013–2015, standardized territory mapping included ≥6 (average: 11.5 ± 0.6 SE) visits to each stream reach, with visits preceding peak incubation, and visits within 4 hr after sunrise to ensure high rates of detection (Mattsson and Cooper 2006). Observations were recorded using in-field spot-mapping on topographic maps during 2009–2011 and with a WAAS-enabled Garmin 60CSX GPS unit (Garmin, Olathe, Kansas, USA) with accuracy ≤5 m in 2013–2015. Given frequent and similar site visits each year and some of the same observers monitoring waterthrushes in multiple years of the study, we can reasonably assume that any changes in annual territory length were not due to using in-field spot-mapping vs. GPS territory mapping methods.
Nest Survival and Productivity
Nest searching and monitoring occurred concurrently with territory mapping to determine waterthrush nest survival and productivity. We typically monitored nests every 3–4 days initially and more frequently as fledging approached (Martin and Geupel 1993). We used nestling morphology to determine hatching date (Mattsson and Cooper 2009). We assumed that an undamaged empty nest had fledged if the nest had been active the day before and had approached the predicted fledging date. We attempted to verify fledging by looking for fledglings or adults carrying food if a nest was believed to have fledged. We counted the number of eggs to determine the clutch size of nests with complete clutches. We used the count of nestlings in the visit prior to fledging as the number of fledglings for each successful nest. Nests were considered successful if they produced at least one waterthrush fledgling, including nests that were parasitized by Brown-headed Cowbirds (Molothrus ater).
Riparian Habitat Quality Assessment
Riparian habitat quality was assessed using a habitat suitability index (HSI) specifically designed for waterthrushes (Prosser and Brooks 1998) and the U.S. Environmental Protection Agency's (EPA) rapid bioassessment for high gradient streams (Barbour et al. 1999). The HSI is a broad-scale evaluation of waterthrush instream foraging and upland habitat suitability on a scale of 0–1 (Prosser and Brooks 1998). The EPA index (range of 0–200) assesses stream quality based primarily on instream characteristics that relate to the abundance and composition of waterthrush aquatic macroinvertebrate prey, and therefore may indicate the relative quality of instream foraging habitat (Wood et al. 2016). Both the HSI and EPA index were quantified in a 100-m stream segment centered on each nest in 2009. In later years we used a 50-m segment centered on each nest to make the indices more sensitive to habitat immediately surrounding waterthrush nests. Mattsson and Cooper (2006) conducted EPA index assessments on stream reaches that were 20 × channel width in length. Our 50 m segment was approximately two-thirds of our average channel width (3.7 m) × 20, in line with our goal of making the indices representative of nesting habitat. We did not collect EPA index or HSI data in 2015 due to time constraints.
Analysis
Territory density and length. To model the effect of gas well development on waterthrush territory density, we used a generalized linear mixed model (GLMM) with study stream as a random effect and year and StreamGas as fixed effects. For these mixed models and all hereafter, we did not test more than one gas disturbance variable per analysis to avoid multicollinearity. The response variable was the number of territories along each stream in each year sampled, with the length of each stream included as an offset. We specified a Poisson distribution based on the absence of overdispersion in the fixed-effects version of this model (Zuur et al. 2009). Modeling was performed using the glmer function in the lme4 package (Bates et al. 2015) in R (R Core Team 2014). For this model, as well as the other mixed models mentioned below, model residuals were evaluated graphically, and we used various data exploration diagnostic tools detailed in Zuur et al. (2010) to ensure that model assumptions were met. Statistical significance (α = 0.05) was assessed via a likelihood ratio test (Zuur et al. 2009). If year was significant, a post hoc contrast Kruskal-Wallis sum rank test was completed with Bonferroni correction using the dunn.test package (Dinno 2016) in R to determine the years between which territory density differed.
We used a gamma family GLMM with stream as a random effect to test whether territory length differed among years using the glmmADMB package (Bolker et al. 2012) in R. Overall statistical significance and post hoc testing for year was done in the same manner as for territory density. To test the hypothesis that territory length would increase with a decrease in territory density (Lack 1954), we related territory length to territory density with an asymptotic Spearman rho correlation test using the R packages coin (Hothorn et al. 2008a, 2015a) and psych (Revelle 2017).
Nest survival, productivity, and source–sink threshold. We used program MARK 7.1 (White and Burnham 1999) to estimate the daily survival rate (DSR) of waterthrush nests in each year of the study. Of 364 total nests across all years, we removed 84 nests that did not reach the egg-laying stage, that were discovered after they had fledged or failed, or that had unknown fates and thus did not meet the assumptions of MARK. We assumed a 29-day nesting period (egg laying = 5 days, incubation = 14 days, nestling stage = 10 days) based on the chronology of nests monitored in our study area to calculate annual nest survival using DSR. We plotted annual nest survival ± standard error (SE) to graphically evaluate trends over time.
We developed a set of 11 a priori candidate models (Buckland et al. 1997) containing temporal, shale gas disturbance, and habitat covariates that we hypothesized might influence the DSR of waterthrush nests. We did not include random effects (i.e. stream) in any model due to the difficulty of modeling such effects in nest survival analyses, but recognize that the random effect of stream could have accounted for variability among study streams if present. All a priori models included temporal covariates to account for their influence on nest survival based on results reported in the literature: nest age, quadratic effect of time within season (TT), and average daily rainfall (Rain). We included nest age because nests may be more vulnerable as they age (Dinsmore et al. 2002, Grant et al. 2005, Burhans et al. 2010) and because the most-supported nest survival model in a 2011 waterthrush benthic aquatic prey study included a similar covariate called nest stage (Wood et al. 2016). We included TT because it was most parsimonious in a post hoc waterthrush nest survival model (Mattsson and Cooper 2009). We included mean daily rainfall (mm) because headwater riparian systems are subject to seasonality and annual changes in rainfall (Richardson and Danehy 2006) that can affect waterthrush nest survival rates (Mattsson and Cooper 2009). For each nest, we averaged daily rainfall estimates across the period during which an active nest was under observation (Mattsson and Cooper 2009). Precipitation estimates were pooled from the 4 Weather Underground (The Weather Company, Atlanta, Georgia, USA) stations closest (average: 36 km) to the study area (3 weather stations in 2009). We included an additional fixed effect of year as a variable of interest in some models because shale gas development increased over the study period (Farwell et al. 2016) and to account for annual variation in DSR associated with biotic and abiotic factors not included in our models. We did not include a model with year only because we a priori evaluated nest survival graphically to review trends and found some overlap in annual estimates.
The primary variables of interest in our models included 3 shale gas disturbance covariates (TerrGas, TerrRunoff, NestGas) and 2 habitat covariates (GasFCD, OtherFCD; Table 2). Gas disturbance covariates were not combined in an additive fashion in a single model because they were related metrics, and the habitat covariates were not combined in a single model as we wanted to distinguish whether the source of forest canopy disturbance was important. We chose GasFCD and OtherFCD as habitat covariates as we hypothesized that shale gas development through the removal and fragmentation of riparian forest cover could negatively influence waterthrush reproduction through modified predator assemblages and activity as well as altered stream hydrology and water quality (Petit and Petit 1996, Mulvihill et al. 2008, Mattsson and Cooper 2009) and because waterthrushes are known to be sensitive to removal of forest canopy cover (O'Connell et al. 2003).
We used Akaike's information criterion corrected for small sample sizes (AICc) to evaluate support for candidate models (Burnham and Anderson 2002) in program MARK. We modeled the binomially distributed data with the user-defined logit-link function while simultaneously considering associations with temporal, shale gas disturbance, and habitat covariates. We assessed the relative plausibility of each model in each model set by comparing Akaike weights (wi). We considered the model with the lowest AICc value to be the best-supported model given the data, and considered any models with ΔAICc < 2 plausible (Burnham and Anderson 2002). We used model-averaged regression coefficients (Burnham and Anderson 2002), and used 85% confidence intervals (CIs) to infer the biological importance of covariates in plausible models as 95% CIs with the information-theoretic approach can lead to variable selection uncertainty (Arnold 2010).
We quantified average overall individual productivity and average annual population-level productivity using an approach similar to that of Boves et al. (2015). Mean number of fledglings per successful nest per male (the per capita value) was multiplied by annual nest survival (DSR29) separately for areas undisturbed and disturbed by shale gas development. Areas undisturbed by shale gas development (n = 78) were categorized as territories with TerrGas = 0%, and areas disturbed by shale gas development (n = 55) were classified as territories with any amount of TerrGas (range: ∼3–100%). For population productivity, the per capita productivity value was calculated per year, multiplied by annual nest survival, and then multiplied by average annual territory density to determine whether average annual population-level productivity changed over time. The significance of individual productivity between areas undisturbed and disturbed by shale gas development and of population productivity across years was evaluated graphically by examining the overlap of 95% CIs for simple biological inference (Payton et al. 2003, MacGregor-Fors and Payton 2013). Productivity SEs used to construct the 95% CIs were SE values for the unadjusted mean number of fledglings to reflect the full range of variability for each metric (T. Boves personal communication).
Additionally, we assessed whether productivity could compensate for adult mortality (e.g., Robinson and Morse 2000) by calculating a source–sink threshold (Pulliam 1988). Since the threshold is the minimum number of fledglings needed to compensate for adult mortality, productivity above or below the threshold allowed us to evaluate whether habitat quality was sufficient for local populations to be maintained. In the manner of Robinson and Morse (2000), the source–sink threshold was the annual productivity per pair necessary to compensate for adult mortality, modeled as 2(1 − φ) / φ0, where φ is adult survival and φ0 is juvenile survival. We calculated overall adult mortality (2(1 − φ)) using the average of our separate estimates of male and female apparent survival (φ). Male survival (φ) was 0.56 ± 0.04 SE and female φ was 0.44 ± 0.08 SE (Frantz 2018). We assumed juvenile survival to be half of the adult value (Nolan 1978) as low estimates (≤0.30) are likely to be more accurate than previously thought (McKim-Louder et al. 2013). The source–sink threshold value was multiplied by annual nest survival (DSR29) to convert it to the same scale as individual productivity values. We then graphically evaluated whether average individual productivity in areas undisturbed and disturbed by shale gas development fell above or below the threshold value and based significance on overlapping 95% CIs for simple biological inference, where a proportion of overlap of ≤0.5 was considered significant (Cumming and Finch 2005).
Riparian habitat quality assessment. We used a beta family GLMM to compare HSI scores for nest-centered stream segments located in territories disturbed and undisturbed by shale gas development using the glmmADMB package (Bolker et al. 2012) in R. Models included TerrGas and year as a fixed effect and stream as a random effect. Because the 2009 HSI score was based on 100-m stream length segments and scores in following years were based on 50-m segments, we used Wilcoxon-Mann-Whitney tests to compare HSI scores from 2009 and 2010, 2 years in which the percentage of stream disturbed was the same. We found no differences (z78 = 1.45, P = 0.15), so we did not account for differences in stream segment lengths in our models. Our global model had marginal overdispersion (χ2240 = 276.81, ratio = 1.15, P = 0.05), so we added an observer-level random effect (OLRE) in which each observation received a unique random effect level to absorb extraparametric variability (Harrison 2015). Statistical significance (α = 0.05) was assessed via a likelihood ratio test (Zuur et al. 2009). If year was significant, we completed a post hoc contrast Kruskal-Wallis rank sum test to indicate the years between which HSI scores differed.
We used linear mixed effects modeling (LMM) in the lme4 package (Bates et al. 2015) in R to assess nest-centered EPA index scores in disturbed and undisturbed territories. Models included the percentage of shale gas territory disturbance and year as fixed effects and stream as a random effect. We did not test other gas disturbance variables within the same model to avoid multicollinearity. A significant difference between EPA index scores from 2009 (100-m stream segment assessment) and 2010 (50-m stream segment assessment) was indicated by t-tests (t78 = 6.12, P < 0.001). Therefore, we dropped 2009 data from our models to avoid variability due to the length of stream segment assessed, and assumed that 2010 was representative of initial gas disturbance to streams because the percentage of disturbance along streams was the same in 2009 and 2010. We used the afex package (Singmann et al. 2015) in R to retrieve P-values for the F-test assessment of fixed effects. We used a post hoc Tukey HSD test in the multcomp package (Hothorn et al. 2008b, 2015b) in R to determine the years between which EPA index scores differed. We set significance at α = 0.05 for all tests. Results are presented as mean ± SE except where noted.
RESULTS
Territory Density and Length
In 6 yr we monitored 400 waterthrush territories. Waterthrush territory density was not related to StreamGas (χ21 = 0.002, P = 0.97), but differed significantly by year and generally declined over time (χ25 = 13.424, P = 0.02; Table 3). A post hoc contrast Kruskal-Wallis sum rank test for year indicated that territory density was significantly higher in 2010 than in 2015 (χ25 = 3.05, P = 0.02). Across years, study streams had a mean of 23% ± 0% of their length disturbed by shale gas development (range: 0–67%; Table 1).
TABLE 3.
Louisiana Waterthrush demography (± SE) across 6 yr in response to shale gas development in the Lewis Wetzel Wildlife Management Area, West Virginia, USA. Population productivity is the mean number of fledglings per successful nest per year multiplied by annual nest survival and average annual territory density per km of stream. Also shown are average annual EPA index and habitat suitability index (HSI) scores ± SE. A larger EPA index or HSI score indicates better riparian habitat quality.
Territory length was significantly different by year and generally increased over time (χ25 = 59.44, P < 0.001; Table 3). A post hoc test for year indicated that territory length was greater in 2009 than in 2010 (z150 = 3.10, P = 0.01), but was less in 2009 than in 2014 (z134 = −3.82, P = 0.001) and 2015 (z125 = −3.84, P < 0.001). Territory length was less in 2010 than in 2011 (z145 = −2.79, P = 0.02), 2013 (z139 = −4.90, P < 0.001), 2014 (z143 = −6.95, P < 0.001), and 2015 (z134 = −6.83, P < 0.001). Territory length also was less in 2011 than in 2014 (z129 = −4.00, P < 0.001) and 2015 (z120 = −4.01, P < 0.001). Territory length significantly increased as territory density decreased (rho = −0.49, z394 = −9.66, P < 0.001).
Nest Survival, Productivity, and Source–sink Threshold
The daily nest survival rate (DSR) and annual nest survival peaked in 2010 and generally declined over time (Table 3). Overall DSR was 96.4% ± 0.3% and yielded average annual nest survival of 34% ± 3%. Across all years, 8 nests were parasitized by Brown-headed Cowbirds, primarily in later years of the study (Table 3).
Of 11 a priori models, 4 models that included habitat (GasFCD) or gas (TerrRunoff, TerrGas, NestGas) covariates were supported (ΔAICc < 2; Table 4). The 2 models that included TerrRunoff or TerrGas accounted for most of the variation in DSR (wi = 0.28 and 0.27). The 85% CIs of model-averaged regression coefficients did not overlap zero for Rain and GasFCD, which had positive influences on DSR, and also did not overlap for TerrGas, TerrRunoff, and NestGas, which had negative influences on DSR (Table 5, Figure 3). About one-third (∼30%) of monitored nests were disturbed by shale gas development (i.e. NestGas), and nest survival dropped from 37% ± 4% in undisturbed areas to 31% ± 5% in areas disturbed by shale gas development. Territories that contained nests had on average 25% ± 2% of the territory disturbed by shale gas development (TerrGas) and 22% ± 2% of the territory affected by contaminant runoff (TerrRunoff). Nests across years were predominantly in forested habitat (94% ± 1%), with 2.3% ± 0.3% of the forest canopy disturbed by shale gas development (GasFCD). Nest age and TT (quadratic effect of time within season) had regression coefficient 85% CIs that overlapped zero, indicating little or highly variable influence on DSR.
TABLE 4.
Model selection results for 11 a priori models examining the nest survival of Louisiana Waterthrushes in an area undergoing shale gas development in the eastern U.S. See Table 2 for model notation. K = the number of parameters in each model, ΔAICc = the difference from the top model in Akaike's information criterion corrected for small sample size, wi = Akaike weight, and Dev is the model deviance.
TABLE 5.
Estimates for nest survival covariates (n = 7) from the top 4 models (see Table 4) based on model-averaged regression coefficients, with unconditional standard errors (SE) and 85% confidence limits (85% CL). Significant covariates with 85% CL that do not overlap 0 are highlighted in bold font.
FIGURE 3.
Four nest survival covariates (with 95% CI) that significantly influenced Louisiana Waterthrush daily survival rate (DSR) in an area undergoing shale gas development in the eastern U.S. The binary response covariate NestGas also significantly influenced DSR (not depicted). See Table 5 for model results and Table 2 for covariate notation.
Overall population productivity was 2.3 ± 0.5 fledglings per km and generally declined from the earlier years of the study (2009–2011) to the later years (2013–2015; Table 3) based on 95% CIs. Overall individual productivity was 1.5 ± 0.1 fledglings per adult male. Individual productivity was higher in areas undisturbed (1.6 fledglings, 95% CI = 1.4–1.8) than disturbed (1.4 fledglings, 95% CI = 1.2–1.6) by shale gas development based on 95% CIs (Figure 4). The completed source–sink equation was 2(1 − 0.50) / 0.25, with 0.50 being the average of male and female adult survival and juvenile survival assumed to be half that value at 0.25, resulting in a threshold of 1.4 fledglings per pair after accounting for overall nest survival. The source–sink threshold of 1.4 fledglings per pair was below the individual annual productivity in areas undisturbed by shale gas (Figure 4), suggesting that these are source habitats. However, the 95% CI for individual productivity in areas disturbed by shale gas (on average, 57% ± 5% of waterthrush territory was disturbed by shale gas development in disturbed areas) overlapped the threshold, suggesting that these are borderline sink habitats.
Riparian Habitat Quality
HSI scores were negatively related to TerrGas (χ21 = 65.34, P < 0.001; Appendix Figure 5) and differed by year (χ24 = 34.84, P < 0.001; Table 3). A post hoc contrast Kruskal-Wallis sum rank test for year indicated that HSI scores were significantly higher in 2009 than in 2013 (χ24 = 4.03, P < 0.001) and 2014 (χ24 = 3.14, P = 0.01).
EPA index scores were also negatively related to TerrGas (F1,158 = 14.54, P < 0.001; Appendix Figure 5) and differed by year (F3,196 = 14.07, P < 0.001; Table 3). A post hoc Tukey HSD test for year indicated that EPA index scores were significantly higher in 2010 than in 2014 (z394 = 3.29, P = 0.005), higher in 2011 than in 2013 (z394 = −4.26, P < 0.001), and lower in 2013 than in 2014 (z394 = 6.18, P < 0.001).
DISCUSSION
Across our 6-yr study in the LWWMA, we saw general declines in waterthrush territory density, nest survival, nest productivity, and riparian habitat quality concurrent with a site-wide increase in disturbance related to shale gas development (Farwell et al. 2016). Our source–sink threshold suggested that individuals breeding in areas disturbed by shale gas development were potentially inhabiting sink habitat and that these populations were more at risk of decline than those in areas undisturbed by shale gas development. Declines in waterthrush demography occurred despite <5% forest cover loss in our predominantly forested study site (Farwell et al. 2016), which suggests that factors other than loss of forest cover also influenced demography (Wood et al. 2016). In general, all demographic parameters for waterthrushes appeared to be negatively affected by shale gas disturbances occurring in headwater stream ecosystems (Table 3). To our knowledge, ours is the first study to have established the potential for Marcellus-Utica shale gas development to affect the reproductive success and productivity of forest birds.
Waterthrush territory density declined across years, but was not explained by the percentage of the stream disturbed by shale gas development (StreamGas). Streams on average had less than one-quarter (23% ± 3%) of their length disturbed by shale gas development (Table 1), and no stream was ever completely disturbed (maximum: 67% disturbed). Consequently, undisturbed areas occurred along every stream, so that waterthrushes could shift their territories to forage and nest along undisturbed sections of streams. Waterthrushes on acidified streams in Pennsylvania, USA, used a similar strategy (Mulvihill et al. 2008). Despite waterthrushes exhibiting high site fidelity (O'Connell et al. 2003), we noted that, by the end of our study, our initial high territory densities of 1.5 km−1 had dropped to 1.0 territory km−1, lower than typical densities across the breeding range of the species (Mattsson et al. 2009). A headwater stream assessment based on waterthrush densities in Pennsylvania found that 0–1 territories km−1 indicated degradation and 1–2 territories km−1 indicated possible degradation (O'Connell et al. 2003). This suggests increased degradation of our study streams across our study period, as was also suggested by our declining HSI scores. The decline in riparian stream quality over time likely influenced the decrease in territory density and increase in territory lengths in our study. Increasing territory length in disturbed areas may be a mechanism that allows waterthrushes to compensate for poor habitat quality (Mulvhill et al. 2008). Waterthrushes that increase their territory lengths may need additional foraging resources to meet minimum breeding requirements, as suggested by greater territory densities in 2011 in areas where macroinvertebrate density, biomass, and stream quality were higher (Wood et al. 2016).
Nest survival was positively influenced by average daily rainfall, similarly to Mattsson and Cooper's (2009) finding of maximum daily survival rate at intermediate (3–10 mm) rainfall levels. This intermediate rainfall range is similar to what waterthrushes encountered during the active nesting period in our study site (range: 0–11 mm; average: 3.6 ± 0.1 mm). Rainfall in this range likely leads to increased prey availability and foraging efficiency, and therefore greater nest vigilance (Mattsson and Cooper 2009). Lack of sufficient water flow was likely a more influential factor than flooding in our headwater system as only 4 nests were confirmed to have failed from high water events. At the beginning of the breeding period our streams were typically flowing, but by late summer, when young were fledging, streams had intermittent or little flow. While not a documented threat to waterthrushes along our headwater study streams, shale gas operations withdraw large amounts of surface water and groundwater from small streams (Entrekin et al. 2011). As such, waterthrushes breeding downstream of water withdrawal operations have the potential to be negatively affected by altered hydrology in the same manner in which water withdrawals affect other species, such as brook trout (Salvelinus fontinalis; Weltman-Fahs and Taylor 2013).
After accounting for the positive influence of rainfall, waterthrush DSR had a significant negative relationship with 3 shale gas development covariates (Table 5). TerrRunoff was a measure of potential surface water contamination, while TerrGas and NestGas assessed the physical presence of shale gas infrastructure and included human activity and sedimentation (Table 2). All previous bird community studies of Marcellus shale drilling in the Appalachians have focused primarily on the presence of gas infrastructure, with less attention paid to noise and light levels (Davis 2014, Barton et al. 2016, Farwell et al. 2016). Waterthrushes in our study could have been directly affected by the presence of infrastructure, given the similar findings of negative effects of oil and gas development on bird species in other regions (Van Wilgenburg et al. 2013, Thompson et al. 2015, Hethcoat and Chalfoun 2015), but we must also consider indirect effects on stream and terrestrial food webs from possible contamination (Entrekin et al. 2011). For example, waterthrushes in areas disturbed by shale gas development had higher levels of barium and strontium in their feathers than waterthrushes in areas undisturbed by shale gas development (Latta et al. 2015). Because barium and strontium are 2 heavy metals associated with the drilling process (Chapman et al. 2012), and the LWWMA was a region sampled by Latta et al. (2015), this finding could be related to our modeling result that TerrRunoff negatively influenced DSR. Heavy metals can interfere with DNA methyl transfer (Hala et al. 2014), so one potential mechanism through which DSR can be affected is by differential methylation via epigenetics.
Previous studies have shown potential waterthrush vulnerability to forest habitat fragmentation (Robbins 1979, McIntyre 1995, Adams 2007) and declines in abundance after loss of ∼16% of forest from the landscape (Becker et al. 2015). In contrast, waterthrush DSR in our study showed a slight positive relationship with forest canopy disturbance caused by shale gas development (GasFCD). Clearance of land for shale gas development in some instances may increase net primary production in streams (Johnson et al. 2015) and may increase the abundance of certain types of aquatic prey, such as shredders (Barton 2016); thus, GasFCD may have potentially increased some aquatic prey taxa for waterthrushes. In addition, Davis (2014) found that, while nest survival was lower for Field Sparrows (Spizella pusilla) in the presence of gas wells, gas pipelines and access roads had a lower index of predation risk, possibly due to increased noise (Francis et al. 2012) or light (de Molenaar et al. 2006) levels, even though predation typically increases near forest edges (Paton 1994). However, these potential benefits from GasFCD could be offset by higher nest abandonment rates (Davis 2014) or cowbird nest parasitism in areas disturbed by shale gas development. Although few nests were parasitized (∼3% of nests) in our study area compared with other waterthrush studies (range: 0–81%; Mattson et al. 2009), we observed an apparent increase in parasitism rates of waterthrush nests across years as forest cover declined, and the majority of parasitized nests (75%) failed. Concurrently, cowbird detections increased, with detections at 2% of sample points in 2008 and at ∼28% of points in 2015 (Farwell et al. 2016). Parasitized nests had double the amount of forest canopy disturbance due to shale gas development (6% ± 2%) than unparasitized nests (2.7% ± 0.3%). Given that the average forest cover at nests was 94% ± 1% and that waterthrushes will not occupy areas with <40% forest cover (O'Connell et al. 2003), GasFCD may have played only a minor role, at least initially, among several factors influencing nest survival (i.e. rainfall and other disturbances associated with shale gas development).
Population productivity generally declined over time, and areas disturbed by shale gas development had lower individual productivity, levels of which broadly overlapped with the source–sink threshold (Figure 3). Corresponding with lower individual productivity, lower nest survival, and decreased riparian habitat quality with increasing disturbance, areas disturbed by shale gas development may have a greater risk of being sink habitat. Our individual productivity values matched those of other studies, which have shown productivity in optimal (or presumably undisturbed) source habitat to be barely above estimated source–sink thresholds (Morse 1996, Holmes et al. 1996). Headwater streams may need to be buffered from potential disturbances if they are only marginally source habitats even under ideal conditions (Robinson and Morse 2000).
Conclusion
Our study is one of the first to demonstrate that shale gas development can affect the reproductive success and productivity of a wildlife population, likely through the presence of shale gas infrastructure and by indirect negative effects on stream health and aquatic prey (Wood et al. 2016). Increasing overall aquatic ecosystem health will necessitate measures to protect water quality from upstream sediment load and pollutant sources (Cook et al. 2015), which will require watershed-scale (Merovich et al. 2013) habitat conservation efforts. Spills and erosion are the most commonly reported environmental violations (Rahm et al. 2015), which could be avoided with setbacks from streams and avoidance of building on steep grades (Evans and Kiesecker 2014). Development is outpacing the implementation of best management practices (Brittingham et al. 2014), so placing well pads farther away from water than currently permissible may be the most effective way to avoid multiple disturbances from shale gas development (Milt et al. 2016).
For effective mitigation strategies at a regional level, additional species- and area-specific studies are needed (Northrup and Wittemyer 2013), as well as clarification of the specific mechanisms involved in species' responses (Hethcoat and Chalfoun 2015) to shale gas disturbance. Our study results combined with postfledging survival data (Streby and Anderson 2011) and a cross-ecosystem evaluation of food web interactions (Soininen et al. 2015) with potential contaminants could fill important knowledge gaps. Lack of information regarding the full range and interdependence of waterthrush demographic responses to shale gas development should not negate immediate risk management activities (Loss 2016), especially if multiple lines of evidence suggest declines and negative demographic responses of this known important bioindicator of headwater stream ecosystems.
ACKNOWLEDGMENTS
West Virginia Division of Natural Resources (WVDNR) provided access to the study area and Wheeling Jesuit University provided access to field housing. We thank many field assistants and graduate students who collected data during our 6-yr study. In particular we thank J. Mizel, D. Becker, and K. Aldinger for assistance with data collection and analysis approaches. We are grateful to L. Farwell, S. Latta, A. Welsh, M. Strager, and S. Welsh for helpful comments that improved the manuscript. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.
Funding statement: Our research was funded by the West Virginia Division of Natural Resources, U.S. Department of Energy National Energy Technology Laboratory, West Virginia University, and National Aviary. None of our funders had any influence on the content of the submitted or published manuscript and only the U.S. Geological Survey (USGS) required approval of the final manuscript prior to publication as required by their Fundamental Sciences Practices protocols ( https://pubs.usgs.gov/circ/1367/).
Ethics statement: Banding was conducted under USGS banding permit #23412 and #23059. This study was completed under the auspices of West Virginia University IACUC protocols #04-0302 and #07-0303.
Author contributions: P.B.W. and G.G. conceived the original idea for the project, with P.B.W. supervising the research and securing most of the research funding. M.W.F., J.S., and G.G. conducted the research. M.W.F. and P.B.W. wrote the paper with input from J.S. Methods were developed by all authors. M.W.F. analyzed the data. P.B.W. and M.W.F. contributed substantial resources.
LITERATURE CITED
Appendices
APPENDIX FIGURE 5.
(A) EPA index and (B) habitat suitability index (HSI) scores in relation to the percentage of Louisiana Waterthrush territory disturbed by shale gas development (TerrGas). The EPA index (range of 0–200) indicates the relative quality of instream foraging habitat for waterthrushes (Wood et al. 2016). The HSI is a broad-scale evaluation of waterthrush instream foraging and upland habitat suitability on a range of 0–1 (Prosser and Brooks 1998).
