Demographic characteristics of an avian predator, Louisiana Waterthrush (Parkesia motacilla), in response to its aquatic prey in a Central Appalachian USA watershed impacted by shale gas development

We related Louisiana Waterthrush (Parkesia motacilla) demographic response and nest survival to benthic macroinvertebrate aquatic prey and to shale gas development parameters using models that accounted for both spatial and non-spatial sources of variability in a Central Appalachian USA watershed. In 2013, aquatic prey density and pollution intolerant genera (i.e., pollution tolerance value <4) decreased statistically with increased waterthrush territory length but not in 2014 when territory densities were lower. In general, most demographic responses to aquatic prey were variable and negatively related to aquatic prey in 2013 but positively related in 2014. Competing aquatic prey covariate models to explain nest survival were not statistically significant but differed annually and in general reversed from negative to positive influence on daily survival rate. Potential hydraulic fracturing runoff decreased nest survival both years and was statistically significant in 2014. The EPA Rapid Bioassessment protocol (EPA) and Habitat Suitability Index (HSI) designed for assessing suitability requirements for waterthrush were positively linked to aquatic prey where higher scores increased aquatic prey metrics, but EPA was more strongly linked than HSI and varied annually. While potential hydraulic fracturing runoff in 2013 may have increased Ephemeroptera, Plecoptera, and Trichoptera (EPT) richness, in 2014 shale gas territory disturbance decreased EPT richness. In 2014, intolerant genera decreased at the territory and nest level with increased shale gas disturbance suggesting the potential for localized negative effects on waterthrush. Loss of food resources does not seem directly or solely responsible for demographic declines where waterthrush likely were able to meet their foraging needs. However collective evidence suggests there may be a shale gas disturbance threshold at which waterthrush respond negatively to aquatic prey community changes. Density-dependent regulation of their ability to adapt to environmental change through acquisition of additional resources may also alter demographic response.


Introduction
The rapid development of hydraulic fracturing techniques in the last decade has allowed the expansion of development for unconventional drilling activity, hereafter shale gas development [1]. The Marcellus-Utica shale basin is one of the largest natural gas plays underlying part of the northeastern United States with substantial growth in gas production [2−3]. As of activity was high [29]. However, Wood et al. [29] spanned only one breeding season with a limited sample size (n = 12 watershed samples) at the watershed-scale, meriting further evaluation with increased sampling efforts at other spatial scales (territory and nest).
As a follow-up to these previous studies, in 2013-2014 we: 1) evaluated the congruence between aquatic prey and riparian quality indices used to gauge waterthrush habitat (i.e., US EPA Rapid Bioassessment Protocol, EPA; Habitat Suitability Index, HSI), 2) evaluated if the amount of shale gas disturbance or potential hydraulic fracturing runoff in a territory or at a nest influence aquatic prey, and 3) quantified waterthrush demographic response to aquatic prey changes. We hypothesized that aquatic prey should be positively linked with riparian quality habitat scores. EPA and HSI scores were negatively affected by shale gas development [28]. As a consequence of habitat degradation, we expected a negative relationship between aquatic prey metrics and the amount of shale gas disturbance or potential hydraulic fracturing runoff in a territory or at a nest. We also hypothesized that clutch size, number of fledglings, and territory density would have a positive association with aquatic prey metrics. Annual territory length increased as territory densities decreased [28], so we expected smaller territories to be indicative of higher quality aquatic prey and stream quality (e.g., [30]). Nest survival was minimally affected by aquatic prey in 2011 [29] but we hypothesized that any stream impairment effects on the aquatic prey would affect nest survival.

Study area
We studied waterthrush along 58.1 km of 1st-and 2nd-order headwater stream tributaries (n = 14) at Lewis Wetzel Wildlife Management Area (LWWMA) located in northwestern West Virginia (Fig 1). Our waterthrush aquatic prey study in 2013-2014 was part of a waterthrush demography study over a six year period (2009)(2010)(2011)(2013)(2014)(2015); [28]). The study area overlays the Marcellus-Utica shale region and occurs where waterthrush reach their highest densities within the central Appalachians [31]. The LWWMA is part of a regional core designated as a priority conservation planning area for both aquatic and terrestrial targets [32].
During our study, shale gas development activities included building of conventional (shallower formations) and Marcellus well pads, timbering for yet unbuilt well pads, the expansion of existing road and pipeline infrastructure, and the construction of new infrastructure. In 2008, the LWWMA was 95.3% forested and had 0.4% shale gas land cover; the first shale gas well development began in 2007 [33]. Between the 2010 and 2011 breeding seasons, shale gas development activities that occurred since 2007 accelerated across the study area and began to increase especially on ridgetops. In 2011, study area-wide shale gas land cover was 1.3% and increased to 2.7% in 2013-2014 ( [33]; Table 1). Starting in 2013, shale gas development abated study area-wide and instead became concentrated to specific streams and ridgetops. Clearing for additional new well pads occurred late (June-July) in the 2013 breeding season with well pad completion in 2014, in addition to re-drilling of an existing well pad. There was no new shale gas development or activity in the 2015 breeding season. In 2015, the LWWMA was 90.8% forested and 3.9% in shale gas development land cover, with 83.1% of shale gas development resulting in direct forest loss [33]. In summary, 2013 disturbances slowed and affected streams more noticeably late in the breeding season, while in 2014 shale gas activity peaked again, particularly at Slabcamp Run, but did not achieve 2011 levels of activity (S1 Table).

Mapping of streams and shale gas disturbance
Within 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 disturbance within the study area for each year of the long-term study, including years 2013-2014 of the aquatic prey study (see Frantz et al. [28] for full description). All forest canopy disturbances were classified as shale gas related (i.e., well pads and associated road and pipeline infrastructure, frequent truck traffic, and visual stream sedimentation) or as being unrelated or pre-existing (i.e., forest roads, recent even-aged timber harvests, and various types of existing clearings). We classified a few conventional impacts (i.e., stream-side vertical pump jacks) as related to shale gas development because their 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 [34]. Gas well records [35] were used to verify target shale formations, drilling status, and start dates for all well disturbances.
Lengths of each study stream (average length 4.1 ± 0.54 km, range 0.95-7.4 km) were calculated in GIS using a 3D functional surface length tool and a 3 m resolution digital elevation model to account for topography, and study streams were defined to have a drainage basin of 9.0 hectares (i.e., <100 ha; [36]) to delineate the uppermost headwater reaches (24 k scale or higher resolution; e.g., [37]). To describe and model waterthrush demography and riparian habitat quality as a function of shale gas disturbance, we created four continuous and one binary variable based on disturbance categories at the stream, territory, and nest scale. The first (termed StreamGas) described mostly localized streamside disturbance indicative of the presence of any shale gas infrastructure or activity. A section of stream was considered disturbed when well pads, infrastructure, or frequent vehicular activity were within 60 m of the stream centerline, which is the typical extent of waterthrush streamside use (i.e., 60 m; [38]). When a stream had visually observable sedimentation that resulted from shale gas development at any distance from the stream, we classified the entire stream network downstream of the sedimentation beginning point as disturbed. Streams were frequently and extensively ground-truthed each season, so there were no stream reaches where sedimentation events were likely to be missed.
We created a second shale gas disturbance category (termed StreamRunoff) that focused solely on potential run-off into streams from shale gas contaminants. A stream was considered disturbed from at and below a well pad or retaining pond (similar to Latta et al. [15]), resulting in the whole downstream network classified as at risk for surface pollution based on elevational maps and ground truthing. This 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 then calculated the proportion of each stream disturbed for each of these two disturbance categories.
We calculated the proportion of each waterthrush territory (a 60-m buffer around each territory vector) that was disturbed by StreamGas and called this metric TerrGas. The proportion Table 1. Louisiana Waterthrush annual demographic, riparian habitat quality, and shale gas disturbance metrics (mean ± SE) at Lewis Wetzel Wildlife Management Area, WV at peak (2011) and later stages (2013-2014) of shale gas development. Our study associated waterthrush response to aquatic prey community changes in relation to shale gas disturbance. All metrics are a subset of those originally reported in Frantz et al. [28] excepting % shale gas land cover which is cited from Farwell et al. [33]. Variable names are defined in S2

Waterthrush riparian habitat quality
Riparian habitat quality was assessed using the Habitat Suitability Index specifically designed for waterthrush (hereafter HSI; [25]) and the US EPA Rapid Bioassessment Protocol for high gradient streams (hereafter EPA [39]) in the same manner as Wood et al. [29] and Frantz et al. [28]. The HSI is a broad-scale evaluation of waterthrush instream foraging habitat, nesting habitat, and upland habitat suitability [25]. The EPA assesses stream quality based primarily on instream characteristics that relate to the abundance and composition of aquatic organisms, and therefore may indicate relative quality of instream foraging habitat for waterthrush [29]. The HSI and EPA indices were quantified in a 50-m stream reach centered on each nest location monitored to make the indices sensitive to habitat immediately surrounding waterthrush nests.

Waterthrush demographic monitoring
We quantified annual waterthrush territory length (m), territory density (# territories/km), and nest survival for our 14 study streams as described in Frantz et al. [28]. Waterthrush territories were delineated as linear vectors along each stream during April 1-June 29 using standardized territory mapping (�6 stream visits [40][41]). Nest searching and monitoring occurred concurrently with territory mapping. Locations of waterthrush observations and nests were recorded with a WAAS-enabled Garmin 60CSX GPS unit with accuracy �5 m in 2013-2015.
To calculate daily survival rate (DSR) for nest survival, we monitored nests typically every three-four days initially and more frequently as fledging approached [42]. We assumed an undamaged empty nest had fledged if the nest was active the day before and had approached the predicted fledge date. Nest sites were revisited at least one more time to verify either no activity or renesting if the nest was not active prior to the expected fledge date. We counted number of eggs to determine clutch size of nests with complete clutches, and the number of fledglings for each successful nest was the count of nestlings in the visit prior to fledging. Nests were considered successful if they produced at least one waterthrush fledgling, including nests parasitized by Brown-headed Cowbirds (Molothrus ater).

Aquatic prey sampling
Aquatic prey occurring in riffle habitat adjacent to nest site locations were sampled once per nest using a Surber sampler. Nest site samples (n = 178) were collected shortly after the nest fledged, failed, or had been abandoned (May 22 -July 28, 2013; June 16 -July 6, 2014) to assess relative prey availability near the time a nest contained fledglings. During sample collection, we scrubbed rocks (>8 cm in diameter) and disturbed sediment 3-cm below the stream bed within the Surber frame for a total of 3 minutes [27]. We separated aquatic prey from detritus for each sample in the field and stored all organic matter in 70% or 95% ethanol.
Post-field season, aquatic prey in benthic samples were sorted, counted, and identified to genus level. Body lengths were also measured to estimate biomass (crayfish excluded). To summarize the aquatic prey composition for each sample, we calculated a family level multimetric IBI called the West Virginia Stream Condition Index (WVSCI [43]), and a genus level multimetric IBI called the Genus Level Index of Most Probable Stream Status (GLIMPSS), version CF, which does not require the genus-level identification of Chironomidae or Oligochaeta [44]. The values we calculated for both indices are based on sampling methods that are slightly modified [45] from the standard methods (i.e., Surber samples and all individuals used in calculations). Thus, they are not strictly interpretable as indicators of stream ecosystem health as originally intended. Nevertheless, they still quantify the composition and integrity of the aquatic prey resource available to waterthrush. We additionally calculated overall aquatic prey density and biomass using length-mass regressions [46]. In total, we selected six aquatic prey metrics to relate to waterthrush demography: WVSCI, GLIMPSS, biomass, density, EPT richness (component of WVSCI), and number of intolerant genera (component of GLIMPSS where pollution tolerance value is <4); S2 Table).

Analysis
We used spatial generalized linear mixed models (hereafter SGLMMs) to assess relationships between waterthrush demography and aquatic prey as well as between riparian habitat quality and aquatic prey for each year (i.e., 2013, 2014) separately. SGLMMs accounted for possible effects of spatial autocorrelation and were modeled using corrHLfit within the spaMM package [47][48] in R [49]. Model residuals were evaluated graphically, extreme or influential data outliers identified graphically and with packages car [50] and stats [49] for potential removal, and other data exploration diagnostic tools were used [51] to ensure model assumptions were met. We used x-y coordinates as a spatial random effect in a Matern correlation model and included a stream random effect. For all SGLMMs, we determined statistical significance of fixed effects using a likelihood ratio test and set significance at α = 0.10 to be cautiously moderate in our assessment of biological significance [52].
We evaluated the degree to which the six aquatic prey metrics were related to riparian habitat quality (i.e., EPA and HSI scores) in individual SGLMM models. We also tested the relationship between the six aquatic prey metrics and the shale gas disturbance metrics for the nest and territory scales (TerrGas, TerrRunoff, and NestGas) as G. Merovich (pers. comm.) found differences in benthic communities up and downstream of shale gas development. We modeled WVSCI, GLIMPSS, biomass, and density using a Gaussian distribution with biomass and density receiving a log10 transformation to approximate normality. We removed an outlier from our benthic density dataset because it was identified as a strong influential outlier not representative of other samples (6422.2 m 2 vs. 354.7 ± 31.3 per m 2 average density) by using the outlierTest and influence.measures functions with packages car and stat. EPT richness and number of intolerant genera were modeled using the Conway-Maxwell-Poisson (COM-Poisson [53][54]) distribution that generalizes the Poisson distribution to handle a wide range of under and over-dispersion typically found in ecological count data [55]. If a COM-Poisson model could not converge, we substituted with a Poisson distributed model (n = 12 models). Because Wood et al. [29] did not assess EPA and HSI in relation to aquatic prey metrics for their 2011 data, we completed a retrospective analysis of their 2011 data. We used Pearson (i.e., EPA) and Spearman (i.e., HSI) correlation coefficients to relate 2011 EPA and HSI to aquatic prey metrics with package psych [56] for comparison to our 2013-2014 results.
We additionally assessed the relationship between waterthrush demography (number of fledglings, clutch size, territory length, and territory density) and the six aquatic prey variables in individual SGLMM models as above. Number of fledglings, clutch size, and territory density were modeled using the COM-Poisson distribution. Territory length was modeled using a gamma distribution. We removed an outlier sample from the fledging dataset where only one fledgling was represented due to Brown-headed Cowbird (Molothrus ater) parasitism. Because the Wood et al. [29] aquatic prey study did not assess territory length, we used a Pearson correlation on data collected in 2011 in the same manner as riparian habitat quality above to relate the six aquatic prey variables to territory length for comparison to our 2013-2014 results.
We used program MARK 7.1 (Colorado State University, Ft. Collins, Colorado, USA, [57]) to estimate daily survival rate (DSR) of monitored waterthrush nests in 2013 and 2014. We removed 63 nests that did not meet the assumption requirements of program MARK and 11 nests that had no benthic information, leaving 107 nests for analysis. We assumed a 29-day nesting period (egg-laying 5, incubation 14, nestling 10 days) based on the chronology of nests monitored on our study area [28].
We developed a set of 7 a priori candidate models [58] that we hypothesized might influence DSR of waterthrush nests based on the results of Wood et al. [29] and Frantz et al. [28]. All covariates are defined in S2 Table. All a priori models included 3 temporal covariates and a shale gas covariate that influenced nest survival in our study area [28]; they included nest age (NestAge), quadratic effect of time of nesting within the breeding season (TT), average daily rainfall (Rain), and percent potential hydraulic fracturing runoff within a territory (TerrRunoff). Instead of an intercept model with no covariates, these 4 covariates formed our base null model given their known importance [28], allowing us to assess whether aquatic prey also influenced nest survival by accounting for them. Nest age indicates vulnerability as the nest ages [59] and within-season trends in DSR reflect dynamic activity patterns of nest predators (e.g., [60]). Mean daily rainfall (in mm) by influencing prey availability affects waterthrush nest survival [28,38] as headwater riparian systems are subject to seasonality and annual changes in rainfall [61]. For each nest, we averaged daily rainfall estimates across the period in which an active nest was under observation [38]. Precipitation estimates were pooled from four Weather Underground, Inc. network stations closest (avg. 36 km) to the study area. Six additive models included the null model plus each of our aquatic prey covariates of interest.
We used Akaike's Information Criterion for small sample sizes (AICc) to evaluate support for candidate models [62] in program MARK. We modeled the binomially distributed data with the user-defined, logit-link function while simultaneously considering associations with the covariates of interest. We considered the model with the lowest AICc value to be the bestsupported model given the data, and any models with ΔAICc <2 were considered plausible. We used regression coefficients and 85% confidence intervals (hereafter CIs) to infer biological importance of covariates in plausible AICc models [63]. We model-averaged NestAge, TT, Rain, and TerrRunoff across all models [62].

Results
Stream disturbance due to shale gas (i.e., StreamGas) was 32.3% in 2011, dropped to 17.3% in 2013, and elevated to 21.5% in 2014, reflecting different levels of shale gas activity despite study area-wide shale gas land cover not changing between 2013 and 2014 (Table 1; S1 Table; Fig 2). The potential for hydraulic fracturing runoff within streams (i.e., StreamRunoff) increased from 19.3% in 2011 to 24.2% in 2013 and 2014 ( Table 1). The percent of each waterthrush territory disturbed by shale gas (i.e., TerrGas) had the same patterns as StreamGas while TerrRunoff increased each year (Table 1). Only 20.3% of territories (39 of 192 total from 2011, 2013-2014) had their full territory length (100%) disturbed by TerrGas or TerrRunoff.
In 2013, aquatic prey biomass and density increased with increasing EPA score, while in 2014 intolerant genera increased with increasing EPA score (Table 2; Fig 2). No relationships were statistically significant in 2013 between HSI and aquatic prey, but in 2014 intolerant genera and WVSCI (approaching significance) increased with increasing HSI score ( Table 2; Fig  2). For 2011 data, aquatic prey biomass had a statistically significant, positive correlation with EPA (R 2 = 0.67, P = 0.02) and HSI (Rho = 0.51, P = 0.09). In 2013, EPT richness increased with increasing TerrRunoff, but in 2014 EPT richness decreased with increasing TerrGas (Table 3; Fig 2). In 2014, intolerant genera decreased with increasing TerrRunoff, TerrGas, and NestGas (Table 3; Fig 2). Raw data used to analyze the relationship between riparian habitat quality and aquatic prey can be found in S1 Dataset (years 2013-2014) and S2 Dataset (year 2011), and S3 Dataset for the relationship between shale gas disturbance and aquatic prey.
All tests for the relationships between clutch size, number of fledglings, and territory density with aquatic prey metrics were statistically non-significant (Table 4). Territory length decreased with increasing aquatic prey density and number of intolerant genera in 2013 (Table 4; Fig 2). For 2011 data, territory length had a statistically significant, negative correlation with GLIMPSS, EPT richness, and number of intolerant genera (R 2 = -0.65, -0.68, -0.67, P = 0.02; Fig 2), respectively. Raw data used to analyze the relationships between waterthrush demography and aquatic prey can be found in S4 Dataset.
Of 7 a priori nest survival models (Table 5) Table. https://doi.org/10.1371/journal.pone.0206077.g002 (wi = 0.17) and GLIMPSS the most in 2014 (wi = 0.18). Regression coefficient 85% CIs overlapped zero for all aquatic prey covariates indicating little, no, or highly variable influence on DSR, but the direction of the relationship between nest survival and aquatic prey switched from negative to positive for 5 of the 6 aquatic prey covariates from 2013 to 2014 (Table 6). In the null base model Rain had positive influence on DSR in 2013 and 2014, while TerrRunoff had negative influence on nest survival in 2014 (Table 6). MARK-formatted files (.inp file extension) used to analyze the relationship between waterthrush nest survival and aquatic prey are the S5 Dataset and S6 Dataset.

Discussion
Shale gas disturbances on our headwater stream ecosystem varied with the intensity of shale gas development that year ( [28]; Table 1; S1 Table). Our follow-up study was able to establish how shale gas alterations to riparian habitat quality and the food web can lead to potential effects at a higher trophic level in an apex predator. By also documenting waterthrush demography decline (Table 1, [28]) and shifts in aquatic prey community structure ( [29]; G. Merovich pers. comm.) due to shale gas development, our study establishes the extent of dependency of waterthrush demographic response and adaptation due to the integrity of ecosystem conditions at the aquatic-terrestrial interface.

Waterthrush foraging resources
Our study builds a connection for decreasing riparian habitat quality due to shale gas altering, at least in part, waterthrush foraging resources. The EPA riparian habitat assessment has been Table 2. Association between waterthrush riparian habitat quality indices (i.e., EPA and HSI) and aquatic prey metrics in spatial generalized linear mixed models. In 2013, aquatic prey biomass and density increased with increasing EPA score, while in 2014 intolerant genera increased with increasing EPA score. No relationships were statistically significant in 2013 between HSI and aquatic prey, but in 2014 intolerant genera and WVSCI (approaching significance) increased with increasing HSI score. Results with P are from a Poisson model. P values of variables that are statistically significant are bolded. Variable names are defined in S2 successfully used in other studies in conjunction with waterthrush occupancy to explain biotic integrity [27]. Higher EPA index and HSI scores were indicative of a larger and healthier aquatic prey community in our system although not with all metrics and statistical significance was dependent on year (Table 2). Additionally, EPT richness and intolerant genera were negatively related to shale gas disturbance, mainly in 2014 (Table 3); this was important to establish since waterthrush riparian habitat quality was negatively affected by shale gas [28]. Overall, HSI was less reliable than EPA for describing aquatic prey, which may be due to HSI consisting of not just food (i.e., stream habitat and quality) scores, but also scores tabulated for waterthrush habitat cover, nesting, and a surrounding landscape classifier [25].

Waterthrush demographic responses
Most demographic responses to aquatic prey were variable or statistically non-significant. Even so, general demographic responses were negatively related to aquatic prey in 2013 then shifted to a positive response in 2014 when shale gas disturbance had a stronger negative influence on aquatic prey and nest survival ( Table 4, Table 5). On streams acidified by mine drainage, waterthrush establish larger territories and forage on peripheral and novel prey items (e.g., terrestrial salamanders) to acquire sufficient prey resources [31]. We saw a similar effect where territory length increased with measures indicating poorer aquatic prey base (e.g., low EPT genus richness). However, territory length in 2014 tended to increase with increasing aquatic prey metrics, opposite of previous years (Table 4). Waterthrush likely had the ability to compensate for loss of food resources by foraging in undisturbed parts of their territory, in addition to increasing territory length, as only 20.3% of territories had their full territory length disturbed by TerrGas or TerrRunoff. The waterthrush's compensation ability in combination with the decline in annual territory density likely contributed to the disassociation between territory length and aquatic prey in 2014. In contrast, given the stronger response and higher territory densities in 2011, under normal territory density conditions (�1.5 territories/ Table 4. Association between waterthrush demographic response (i.e., clutch size, number of fledglings, territory length and territory density) and its aquatic prey in spatial generalized linear mixed models. All tests for the relationships between clutch size, number of fledglings, and territory density with aquatic prey metrics were statistically non-significant. Territory length decreased with increasing aquatic prey density and number of intolerant genera in 2013. Results with P are from a Poisson model. P values of variables that are statistically significant are bolded. Variable names are defined in S2 km) the hypothesis of smaller territories indicate higher quality habitat and foraging resources [30] likely still holds true.

Shale gas disturbance influences on nest survival and aquatic prey
Models used to explain nest survival were also dependent on year (Table 5) with TerrRunoff significantly decreasing daily survival rate in 2014 (Table 6). Our study aligns with Wood et al. [29] in that aquatic prey likely is less influential on nest survival than temporal effects like rain or shale gas disturbance [28]. While our waterthrush-related shale gas disturbance metrics (i.e., TerrGas, TerrRunoff, and NestGas) suggest conflicting yet overall weak negative effects on aquatic prey (Table 3), aquatic prey community structure at our study area quantified upstream and downstream of shale gas at a subwatershed scale also mirrored shale gas activity: community changes differed the most in 2011 [29], were weaker in 2013, and then slightly stronger in 2014 but not as much as 2011 when shale gas activity was highest (G. Merovich, pers. comm.; S1 Table).

Implications
The year-to-year waterthrush demographic responses to aquatic prey in our study were not strongly proportional but instead followed relatively weaker patterning. Timing of benthic sampling in 2013 in relation to shale gas activity levels likely in part explain the lack of a clear signal between waterthrush demography and its aquatic prey. New shale gas activity in 2013 was not as evident until near or after sampling late in the breeding season (S1 Table), and shale gas well pad construction and drilling typical of our study site and elsewhere occur in "pulses" [5]. With our sampling design, we were able to detect benthic community responses as stronger in 2014 than 2013 (similar to G. Merovich, pers. comm.), but increased sampling efforts during appropriate sampling periods may be even more critical for higher food web organisms in overcoming variability of demographic response to shale gas. Our study, through collective evidence, suggests the potential for localized negative effects to aquatic prey from shale gas development, in particular EPT and intolerant taxa that are believed to be the waterthrush's preferred prey [24]. Additionally, the shift in demographic response in 2014 when shale gas disturbance had stronger negative effects on aquatic prey and nest survival may suggest a shale gas disturbance threshold ( [64] ; Fig 2) at which waterthrush respond and adapt to aquatic prey in the same manner aquatic prey community structure concurrently reflected levels of annual disturbance ( [29]; G. Merovich, pers. comm.). Waterthrush are most likely to forage in locations that have higher EPT and intolerant genera [29], making it important to maintain or improve riparian habitat quality linked directly with their aquatic prey. In consideration of population regulatory mechanisms (e.g., [65]) that may influence annual demographic response and adaptability, continued long-term monitoring will be required to discern if a threshold of shale gas disturbance exists that alters aquatic prey communities and, in turn, affects demography of higher-level trophic linkages [66−67]. To some degree, waterthrush in our system appear to have the ability to adapt and meet their foraging needs. However, based on the response in 2011 and 2014 when aquatic prey was the most affected by shale gas, mechanisms used to compensate when stream disturbance is conservatively � 25% (Fig 2) may be less dependable, altering demography. The fact that benthic communities even within pristine streams may be at risk when isolated within heavily Table 6. Annual waterthrush nest survival covariates found in the top supported (ΔAICc <2, n = 6) and unsupported (n = 1) AICc models based on regression coefficients, standard error (SE), and 85% confidence intervals. In the null base model Rain had positive influence on daily survival rate (DSR) in 2013 and 2014, while Terr-Runoff had negative influence on nest survival in 2014. Significant covariates with non-overlapping confidence intervals are bolded. Covariates included in every model to account for their influence (i.e., Rain, NestAge, TT, and TerrRunoff; [28]) have model-averaged regression coefficients and unconditional SEs. Variable names are defined in S2 impacted regions [68], and the tendency for upper reaches of Appalachian headwater streams to have resource extraction activities [69], warrants more attention to multi-dimensional wildlife community responses within aquatic-terrestrial linkages associated to shale gas development.
Supporting information S1 Table.