Biofilm collection

For each stream site, we collected two biofilm samples (one from the downstream and one from the upstream boundaries of the pool), extracted DNA, and used Illumina sequencing to evaluate a 16S rDNA molecular marker that delimits representative biofilm communities. The biofilm samples (2/site; N = 16) averaged 153 mg/sample (wet weight), with significantly greater amounts from downstream sections of pools as compared to those upstream (average lower = 174.4 mg, average upper = 131.1 mg; F = 7.09, P < 0.011, one-way ANOVA [21]). DNA concentration per sample averaged 39.8 ng/μl and did not differ by site or pool location (results not shown).

Univariate analyses of site, hydrology, land use, and stream chemistry

Microbial community composition

We performed Illumina sequencing of a 16S rDNA marker as a means of identifying and quantifying microbial biofilm communities at each site. De-replication (i.e., merging of identical reads) condensed the data by 89% [i.e., from 761,914 reads into a unique set of 83,441 OTUs (operational taxonomic units)]. Elimination of singletons (i.e., OTUs that occurred but once) further reduced the total to 48,802 (a 41.5% reduction). The removal of chimeric sequences (i.e., hybrid sequences consisting of multiple OTUs) eliminated an additional 3753 (7.7%). A comparison of sequences against a reference database excluded an additional 50 (0.1%), and alignment with the core set database [22] removed an additional 345, yielding 6965 unique OTUs as a final total.

We generated rarefaction curves that estimated alpha-diversity for each site to determine whether depth of sampling and sequencing were sufficient to adequately capture microbial community diversity. These curves approached horizontal asymptotes when plotted against number of sequence reads, suggesting sufficient sequencing depth (Fig. 2). A total of 46 phyla were represented, with 24 of these found across all samples. Average per sample = 36 (range = 32–39), with several phyla dominating across all samples: Cyanobacteria (37.4%); Proteobacteria (31.7%); Bacteroidetes (7.6%); Planctomycetes (5.3%); and Actinobacteria (4%) [21].

A total of 141 microbial classes were also represented, with 54 found across all samples. Those with average abundance >2% (N = 20) are presented in Fig. 3. Of these, Alphaproteobacteria was the most dominant, averaging 18.9% across samples, with Betaproteobacteria averaging 8.4%. A total of 310 genera were subsequently identified, with 116 (37%) identified across all sites and 297 (95.8%) found at ≥ 4 sites [21].

Univariate analyses of biofilm communities

The class Synechoccophycideae was represented by six genera, listed in descending abundance as: Arthronema, Acaryochloris, Leptolyngbya, Pseudanabaena, Paulinella, and Synechococus. In turn, seven genera composed the class Oscillatoriophycideae: Microcystis, Chroococcus, Cyanobacterium, Chroococciddoipsis, Phoridium, and Planktothrix. Microcystis was particularly elevated at PICZ-site 2-D (at 5.54%) [21]. Of the five Spartobacteria genera, two were identified as Xiphinematobacter and Chthoniobacter (family Chthoniobacteraceae), while the remaining three were not identified to genus. The implications with regard to the abundances of these microbial classes and genera between Groups are discussed below.

Multivariate comparisons among group

In Fig. 4a, a bi-plot depicts relationships within and between Groups based upon the first two principle components (PCs) of the stream morphology, anthropogenic land use, and water chemistry variables. Sites are identified according to Group (number) and Site (letter) with MICZ-sites in blue (1-A through 1-D), and PICZ-sites in red (2-A through 2-D), respectively (per Table 1). PC-1 accommodated 60% of the variation in the data, and PC-2 absorbed an additional 17% (77% total). MICZ-sites clustered to the positive (right) side of PC-1 with congruent loadings for ‘slope,’ ‘elevation,’ and ‘%-forest.’ Separation on PC-2 was more prominent for MICZ-sites, largely due to the negative values that associated sites 1-B and 1-C with ‘%-forest’ and ‘slope.’ On the positive side of PC-2, MICZ-sites 1-A and 1-D were and allied with ‘elevation.’

PICZ-sites grouped instead to the far left of the PC-1 axis, quite distinct from MICZ-sites. They still separated into quite distinct pairs, with sites 2-A and 2-C on the negative side of this axis and consistent with vectors depicting ‘watershed’ size and ‘%-urban.’ PICZ-sites 2-B and 2-D fell more distant on the positive side of the PC-2 axis, and in alliance with vectors depicting ‘conductance,’ ‘total nitrogen’, total phosphorus,’ and ‘%-pasture.’ The acute angles of these four vectors reflected their close correlation. In this regard, PICZ-site 2-D was more strongly affected than 2-B. Scores on PC-1 differed significantly by Group (P < 0.003; results not shown), whereas those for PC-2 did not.

In Fig. 4b, a second biplot depicted relationships within and between Groups, but in relation to the composition of their bacterial communities, with MICZ-sites in blue and PICZ-sites in red (as above). PC-1 accommodated 58% of the variation in the data, and PC-2 absorbed an additional 25% (83% total). Of the 20 bacterial classes evaluated, 16 clustered quite closely with one another and were represented by an ellipse in the plot. Four bacterial classes clearly separated from the ellipse, with arrows designating the magnitude and direction of their trajectories. PICZ-sites 2-B, 2-C, and 2-D aligned with vectors depicting classes Synechoccophycideae and Oscillatoriophycideae, whereas site 2-A grouped within the ellipse. MICZ-site 1-C was well separated and in conjunction with the class Spartobacteria, whereas class Planctomycetia separated but little from the ellipse. MICZ-site 1-D fell at the edge of the ellipse, but sites 1-B and 1-A were more distant, with 1-B particularly so.

Linear iscriminant analysis Effect Size (LEfSe) analyses corroborated much of the above by delineating 20 OTUs with an LDA score > 2.7. Subclasses Oscillatoriophycideae and Synechococcophycideae (Cyanobacteria), and Roseiflexales (a filamentous bacteria often found with Synechococcophycideae and deemed tolerant of eutrophication and/or poor water quality) were abundant at PICZ-sites. Bioindicators of healthy streams [i.e., families Rhodocyclaceae (Proteobacteria), Stigonematales (Cyanobacteria), and Rivulariaceae (Cyanobacteria)] were abundant at MICZ-sites, as were ‘negative bioindicators’ (N = 5; primarily Spartobacteria) whose abundances covary negatively with particular impacts such as elevated pH or salt concentrations.

Ecological variation among study sites

Instead, we found significant differences among several test variables, as evaluated by Group. For example, MICZ-sites reflected catchments with significantly greater ‘%-forested’ area, but significantly less ‘%-pasture’ (Table 2). Of interest is the fact that several other variables showed elevated but not significantly different values, as gauged by the Bonferroni-corrected probability value for multiple comparisons. We comment on this situation below.

Significant environmental differences between the two groups were also noted when multivariate analyses incorporated the ten variables across the three categories. Sites separated along PC-1, with strong positive (MICZ) and negative (PICZ) loadings manifested according to stream morphology, anthropogenic land use, and water chemistry (Fig. 4a). There was also considerably more variance among PICZ-sites, with paired catchments (i.e., 2-A/2-C and 2-D/ 2-B) well separated on PC-2. Somewhat surprisingly, these sites also showed a strong and concerted response to those variables not deemed significant in the univariate analyses. MICZ-sites displayed much less variability, yet were similarly separated on PC-2 according to a composite of variables that were significant (i.e., ‘%-forest’) and non-significant (i.e., ‘slope’ and ‘elevation’).

The univariate statistics provided differentiation by Group according to individual variables evaluated singly whereas the multivariate analyses provided broader patterns much more interpretable at the ecosystem level, yet not apparent from the separate univariate analyses. This was due largely to the reduced degrees of freedom in the univariate analyses, as constrained by small sample sizes parsed between groups and gauged with Bonferroni-adjusted probabilities. Although our multivariate analyses did not provide statistical probabilities within a hypothesis-testing framework, they more easily depicted the disparity within- and among-Groups, as promoted by watershed, land use, and water chemistry.

The variability in biofilm communities among study sites

Having established the environmental context for our sampling sites by Group, we could then contrast their biofilm communities (Fig. 4b). Groups again separated in multivariate space, albeit less distinctively and with a greater spread among PICZ-sites along PC-1. PICZ-sites also exhibited less variation than did MICZ-sites along PC-2. Clearly, microbial composition varies both among- and within-sites, but with different microbial taxa driving this result in each Group.

For example, Spartobacteria (non-significant in the 1-way ANOVA) clearly associated with MICZ-site 1-C, and differentiated it from all others, whereas site 1-B (associated with 1-C in Fig. 4a) was arrayed quite distantly from other MICZ-sites on PC-2. Furthermore, one site from each Group (i.e., 1-D/ 2-A) fell close to the origin of the PC-axes, suggesting a low overall diversity in their microbial composition (results not shown).

These differences indicate a disparity between biofilm community composition and environmental variability among sites, a result that parallels those from other studies. For instance, microbial assemblages will often group according to historical (i.e., phylogenetic) events [28], but also in response to more contemporary biochemical conditions associated with streams [29]. In this study, numerous factors (per Tables 1, 2 and 3) obviously impact the relationship between biofilm communities and their environment. This result confounds any ‘cause-and-effect’ scenarios for the observed patterns. However, we draw two strong conclusions from the multivariate analyses: PIZC-sites display greater ecological variability within the environmental matrix (Fig. 4a), yet are much less variable when embedded within the biofilm community matrix (Fig. 4b).

An alternate approach would be to contrast our results with those from studies in other shale plays that employed biofilm communities as bioindicators of fracking. Yet many of the latter are tangential to the present study, in that they examined biofilm communities in either flowback [30] or produced waters [30, 31]. Few evaluated stream catchments into which groundwater from fracked sites would eventually percolate, as herein. One study that did so evaluated headwater streams in the Marcellus shale (PA) (Fig. 1a), and found significantly lower species richness and evenness values at sites impacted by fracking [24]. Several of these sites also contained an abundance of bacterial OTUs that correlated positively with decreasing pH, suggesting more acidic stream environments. The diversity of carbon sources available in a stream promotes the functioning of its biofilm communities, and as diversity decreases, so does the community [32]. In this sense, a reduction in carbon sources would be an ecological explanation for the observed reduction in species richness at these sites, although this was not stated as such.

The richness and evenness of species within biofilm communities

In our study, the species richness, evenness, and number of OTUs in biofilm communities were not significantly different when compared between MICZ- and PICZ-sites. Yet such comparisons often mask the interactions among OTUs within these communities. For example, a decrease in abundance of some taxa can also stimulate growth in others normally more rare, a situation that would promote rather than depress evenness [29, 32]. Those streams with lowest values for evenness in each of our Groups [i.e. 1-B and 2-D; 21] may indeed reflect this consideration. For example, 1-B is a headwater stream (stream order = 1) with the greatest ‘%-forest’ in the study (=96%) both of these environmental aspects would promote deposition of leaf litter into the stream that, in turn, must be decomposed. This similarly constrains the biofilm community.

The two least diverse streams in an ecological sense (i.e., 1-A and 1-B) also had low numbers of OTUs [21], again suggesting the potential for a reduction in available nutrients [19]. In a similar vein, PICZ-site 2-D had the highest value for ‘%-pasture’ in the study, and was also associated with elevated levels of available phosphates and nitrates (Fig. 4a), both of which can promote a few dominant species. This was represented at PICZ-sites by the elevated abundances of two Cyanobacterial classes (i.e., Synechococcophycideae and Oscillatoriophycideae). Cyanobacteria are primary producers that seemingly track the significantly elevated levels of nitrogen and phosphorus found in these streams.

A second but related limitation with regard to species richness and evenness is the strong competition among bacteria and hyphomycetes (stream fungi), as promoted by the reduction in dissolved organic matter (DOM) [27]. Dissolved nitrogen primarily exists as nitrates within ground and surface waters, and must be transformed by microbes before entering into and moving through the ecosystem. This, in turn, could promote bacterial OTUs more strongly competitive, at the expense of those less competitive, a situation that would also constrain biofilm community diversity. In addition, and as a second consideration, elevated nitrogen levels are often associated with UNG well sites [19].

Additionally, the removal of pollutants can, paradoxically, also reduce microbial diversity and evenness [33], suggesting (as above) that external sources of carbon can promote the development of OTUs normally more rare. These caveats, in turn, provide numerous potential corollaries to explain the low values for evenness at sites, particularly when carbon sources have become more limited due to fracking [24, 32].

Biofilm communities as bioindicators

The function of many bacterial lineages is not well understood at the ecosystem level, despite their abundances in soil and aquatic systems, and this in turn makes it more difficult to ascertain their status as potential bioindicators. Despite this, general functions are indeed assignable to some clades. Many Synechococcophycideae, for example, employ unique metabolic pathways that allow them to persist in highly acidic environments such as volcanic seeps. The Oscillatoriophycideae is an equally diverse clade that can also serve a bioindicator for organic pollutants. For example, Microcystis (a genus of Oscillatoriophycideae) is abundant at PICZ-sites, and its presence may point to the presence of elevated polycyclic aromatic hydrocarbons (PAHs) that in turn promote its growth [34].

In addition, the genomes of aquatic Spartobacteria encode for a diversity of glycoside hydrolases that are employed in the degradation of complex carbohydrates [35]. This physiological aspect also explains its common co-occurrence with Cyanobacteria, in that the former metabolizes the complex carbohydrates produced by the latter [12]. Spartobacteria should thus positively correlate with Cyanobacteria at PICZ-sites, but was instead found to be significantly reduced. Brine contamination is a well-known fracking by-product, and it continues to be pulled upwards from deeper strata long after drilling has subsided [36]. In addition, PICZ-sites also reflected greater conductance in their water chemistry. Spartobacteria has a pronounced intolerance for salt, and these environmental conditions at PICZ-sites would impede its expected proliferation.

In an attempt to gain a more comprehensive perspective, we can also contrast results from this study with those from earlier studies at the same sites. For example, UNG development had a definite impact on stream macroinvertebrate communities, with short-lived generalists being more abundant at those sites [20]. Yet these effects were difficult to parse across specific taxa, or to specifically associate with the benthic habitat found at PICZ-site.

A second study [19] found increased primary production and eutrophication at sites impacted by UNG activities, and this was interpreted as a potential response to the enhanced levels of nitrogen these sites displayed. Our data support these conclusions in that two classes of Cyanobacteria were clearly more abundant at impacted sites, suggesting the presence of an environment that is beneficial for primary production. Our statistical analyses also verified significant levels of ‘total nitrogen’ and ‘total phosphorus’ at these sites, as well as heightened conductance.

Overall, the observed differences between MICZ-sites and PICZ-sites may reflect the accessibility of sites chosen for UNG well construction, and as such, may add an additional consideration for the design controlled studies to gauge the effects of fracking (see Discussion).

Sampling sites and environmental data for catchments

Eight sites from an ongoing stream ecology project [19, 20] (Fig. 1b, c, d) were assigned to ‘Group’ using two parameters that relate to UNG-well activity: ‘well density’ and ‘inverse flow length’ (IFL). ‘Well density’ is defined as the number of UNG well sites within a 1-km² radius (=catchment area), whereas ‘IFL’ represents the length of flow from each well site to the stream channel, corrected for slope, and calculated for wells upstream of each sampling location via the flow length tool in the ‘Spatial Analyst Toolkit’ of ArcGIS [19]. The inverse of each flow length was summed across all well sites for each catchment area, such that wells more proximal had a higher value that corresponded to a greater potential effect. Sites with an IFL < 0.25 and Well Density (no./km2) < 0.5 were scored as ‘1,’ and designated as a ‘minimally-impacted catchment zone’ (MICZ), whereas those with an IFL ≥ 0.25 and a Well Density ≥ 0.5 were scored as ‘2’ and grouped as a ‘potentially-impacted catchment zone’ (PICZ).

We characterized ten variables distributed across three categories at each site so as to ascertain if designated Groups differed with regard to environmental or anthropogenic factors that could, in turn, affect microbial communities. The first category related to stream morphology and employed four variables (i.e., ‘elevation,’ ‘stream order,’ ‘%-slope,’ and ‘watershed area’). The second utilized three variables that summarized anthropogenic land use (i.e., ‘%-forest,’ ‘%-pasture,’ and ‘%-urban’). The third recorded three water-chemistry parameters (i.e., ‘total nitrogen,’ ‘total phosphorus,’ and ‘conductivity’) deemed important in gauging relationships between stream metabolism and bacterial communities [19]. Abundance of nutrients was measured as μg/L, whereas dissolved salt/ions was in microSiemens (uS)/cm, with higher values signaling an elevated presence of ions.

Biofilm collection, DNA extraction, and Illumina sequencing

At each site, a pool was identified peripheral to the greatest stream flow and a biofilm-covered rock was then selected at downstream (lower) and upstream (upper) boundaries and scrubbed with a sterile Nasco Whirl-Pak Speci-Sponge™. Sponges were immediately re-sealed in the sterile Whirl-Pak and placed onto dry ice for transport to the lab where they were stored at −80°C until processed. For DNA extraction, 20 ml of phosphate buffered saline solution (PBS; 137 mM NaCl, 2.7 mM KCl, 4.3 mM Na₂HPO₄, 1.47 mM KH₂PO₄, pH 7.4) was added to each sample, and the sponge squeezed manually for 5 min to suspend biofilm. Suspensions were transferred to individual centrifuge tubes and pelleted by centrifugation (8000 g for 20 min), with biofilm quantified via wet weight (mg). Standard laboratory protocols were used for all procedures to prevent sample contamination.

DNA from pelleted biofilm was extracted for all 16 samples (2 per site) using a MOBIO commercial kit (PowerBiofilm® DNA Isolation Kit) following manufacturer’s instructions. DNA was quantified (ng/ul) using a Qubit 2.0 Fluorometer (Invitrogen®). Extractions were subjected to PCR using primers that amplified the hyper-variable V4 region of the 16S structural subunit rRNA gene [37]. Multiplexed 16S metagenomic libraries were constructed using standard Illumina protocols, and were sequenced on an Illumina MiSeq platform. Raw Illumina reads were de-multiplexed (MiSeq Reporter software™) and downloaded from the Illumina BaseSpace® cloud.

Bioinformatics

Sequences were trimmed to 251 bp and quality filtered at an expected error of <1% using USEARCH v8.0 [38]. A pipeline developed by the Brazilian Microbiome Project [39] was employed to correct any Illumina formatting issues for subsequent analyses in QIIME v1.7 [40]. OTUs were selected with the UCLUST method (as implemented in QIIME) and taxonomy assigned using the Greengenes 16S rRNA gene database [22], with subsequent conversion into an OTU table (QIIME).

Univariate analyses

Prior to analyses, nine variables were transformed: Percentages (N = 4) were arcsin transformed to radians; areas (N = 1) reduced to square root; and quantitative variables (N = 5) transformed to log₁₀. ‘Stream order’ was evaluated as recorded. Each category was test by Group using a 1-way analysis of variance in R [41], with statistical significance assigned according to Bonferroni-corrected probabilities.

Shannon entropy was computed in QIIME to gauge the number of unique bacterial taxa in each community (i.e., richness) and the evenness of their distributions, with results compared by Group using a 1-way ANOVA in R. Species richness (with repeated subsampling) was then plotted by site as rarefaction curves, so as to estimate whether sampling at each site was of sufficient depth to accurately characterize biofilm communities. Analyses were carried out with the default number of Monte-Carlo permutations (N = 999) at a p-value of 0.05. UniFrac analyses (in QIIME) were used to derive beta (or between sample) diversity estimates using both unweighted data (i.e., OTU presence/absence) and weighted (by relative abundance) [42]. To identify potential bioindicators, a heat map was generated in QIIME using the 20-most abundant taxonomic classes of bacteria. Potential bioindicators were then identified and compared by Group using a 1-Way ANOVA in R with Bonferroni-corrected probabilities.

Multivariate analyses

A principal components analysis (PCA) was performed using a matrix of correlations among sites based on the ten variables across the three categories (i.e., stream morphology, anthropogenic land use, and stream chemistry) using library “prcomp” in R [41]. The first two principal components depicted relationships among the eight sites (i.e., PC-scores) and were contrasted against principal component loadings for the variables. Both scores and loading were visualized in a single plot (hence the term, ‘biplot’), so as to promote the interpretation of the component axes in relation to the variables. Those in the biplot were represented as vectors, and the angle at their origin(s) reflects pairwise correlations (i.e., the more acute the angle, the greater the correlation). We then compared the first six principal components by Group in R, using a 1-way ANOVA with Bonferroni-corrected probabilities.

A principal component analysis was also used to contrast densities of the 20-most abundant bacterial classes among study sites (using library “prcomp” in R [41]), with densities represented as arcsin-transformed percentages. The first two principal components depicted relationships among the eight sites (i.e., PC-scores) and were contrasted against principal component loadings for the 20-most abundant classes.

The biomarker discovery algorithm LEfSe (Linear discriminant analysis Effect Size) was used to designate potential bioindicators among biofilm communities [43]. The program employs a linear discriminant analysis (LDA) with effect size estimated by linking output to the level-6 (Kingdom to Genus) taxonomic summary in QIIME. Parameters employed were: an alpha value of 0.05 for the Kruskal-Wallis (KW) test, an LDA score threshold of >2.7, and a pairwise Group-comparison. Initially, LEfSe conducts the KW rank-sum test as a means of detecting OTUs that differed significantly in abundances between Groups. Biological significance was then investigated with the (unpaired) Wilcoxon rank-sum test. Finally, LDA was then employed to evaluate each OTU with an effect size > 2.7, and with biological indicator gauged via habitat and metabolism.