Sediment Microbial Diversity in Urban Piedmont North Carolina Watersheds Receiving Wastewater Input

Urban streams are heavily influenced by human activity. One way that this occurs is through the reintroduction of treated effluent from wastewater treatment plants. We measured the microbial community composition of water, sediment, and soil at sites upstream and downstream from two Charlotte treatment facilities. We performed 16S rRNA gene sequencing to assay the microbial community composition at each site at four time points between the late winter and mid-summer of 2016. Despite the location of these streams in an urban area with many influences and disruptions, the streams maintain distinct water, sediment, and soil microbial profiles. While there is an overlap of microbial species in upstream and downstream sites, there are several taxa that differentiate these sites. Some taxa characteristics of human-associated microbial communities appear elevated in the downstream sediment communities. In the wastewater treatment plant and to a lesser extent in the downstream community, there are high abundance amplicon sequence variants (ASVs) which are less than 97% similar to any sequence in reference databases, suggesting that these environments contain an unexplored biological novelty. Taken together, these results suggest a need to more fully characterize the microbial communities associated with urban streams, and to integrate information about microbial community composition with mechanistic models.


Introduction
Quantifying microbial species diversity is critical for an understanding of complex ecological systems since microorganisms are important participants in nutrient cycling, energy flow, and contaminant fate. Microbial diversity has been described for a variety of freshwater ecosystems using multiple techniques [1][2][3][4][5], and there exists a similarity in taxa across streams and within similar habitats [6,7]. In rivers, microbial communities are dominated by Actinobacteria, Cyanobacteria, and Bacteriodetes regardless of watershed location, suggesting a common community of organisms [2].
Despite these similarities, there are multiple environmental discontinuities that can result in differences in species composition depending on the stressors associated with individual river systems [8][9][10][11].
Microbial communities in urban streams are subjected to multiple environmental stressors related to land use change, and urban infrastructure such as sewer lines and stormwater control measures. Microbial diversity varies with land use and similar to other taxa, intolerant species in streams are lost at approximately 12% impervious cover such as roads and other hard surfaces [5,12]. These human influences have resulted in the occurrence of human-dominated taxa, such as Bacteroidetes, Acinetobacter, Aeromonas, Trichococcus, and Proteobacteria [13] in urban streams, and a focus on using stream microbial communities to develop indicators of stream health or human impacts [5].
Treated water from wastewater treatment plants (WWTPs) are common sources of point-source pollution in urban streams and different members of the microbial community have different responses to perturbations from treatment plants. Effluent results in the decrease of some taxa while increasing other taxa immediately downstream of the input [5]. In some cases, previously distinct microbial populations upstream to wastewater discharge become indistinguishable downstream from other streams receiving effluent, with proportional reductions in diversity [14]. Communities closer to the outfall resemble WWTP communities and this similarity decreases further downstream [15]. Overall, effluent may act to homogenize communities near the treated wastewater input [14].
The impact of wastewater effluent on stream microbial communities has also become a research focus due to the emergence of antibiotic resistance as a public health threat. Several studies have quantified low concentration of antibiotics in effluent that can decrease microbial diversity, shift bacterial composition, and increase antibiotic resistance in bacterial communities [15][16][17][18][19][20][21][22].
To understand these impacts and feedbacks we have quantified stream microbial diversity in upstream, wastewater effluent, and downstream environments. In a previous paper by Lambirth et al. [18] we reported changes in microbial abundance and ambient antibiotic concentrations in urban stream water in response to the introduction of wastewater effluent. The current study expands on this previous work by including a characterization of the impact of wastewater effluent on the surrounding soil and sediment microbial communities. We hypothesized that downstream microbial communities will undergo changes as a result of continuous long-term exposure to treated wastewater, which contains significantly elevated concentrations of antibiotics relative to upstream water, as well as evidence of human gut and wastewater treatment plant-associated microbes [18]. Integrating stream sediment communities into our understanding of the impacts of wastewater effluent is important as stream sediment microorganisms represent most of the microbial function and are part of the base of the food web to higher trophic levels.

Sample Collection
Sites chosen for this study were identical to those as previously described [18]. Briefly, soil and sediment were collected at two locations upstream and two locations downstream from treated wastewater release for both Sugar and Mallard Creek facilities ( Figure 1). Soil was collected from the stream bank at the midpoint from the waterline to the upper wetline, while sediment was collected at the stream's midpoint sediment bed. Sterile core collectors were fashioned from 15 mL spin tubes with their conical bases removed and were used to collect soil and sediment cores from the surface to approximately 12 cm down. Three replicates were collected from each site during each of the four time points beginning with late winter (time point 1), early and late spring (time points 2 and 3, respectively), and mid-summer (time point 4) of 2016. Tubes were stored in sterile collection bags before transportation back to the laboratory for DNA extraction. All samples referenced to the same time point were collected on the same day.

DNA Extraction
Whole-genome DNA extraction was carried out using DNeasy PowerMax Soil DNA extraction kits (MoBio) according to the manufacturer's instructions. Approximately 10 g of soil or sediment from each sample replicate was loaded into 50 mL extraction tubes and vortexed for 10 m at maximum speed. Soil and sediment were pelleted by centrifugation at 2500× g for 5 m, and the supernatant was passed through the on-column capture filter and washed. The DNA was eluted with 5 mL of sterile water, and quantification was completed using a Qubit fluorometer (Thermo-Fisher) and Nanodrop spectrophotometer (Thermo-Fisher) prior to sequencing library generation. For detailed methods regarding the handling of water samples, see Lambirth et al. [18]

16S rRNA Library Preparation and Illumina Sequencing
All sequencing was performed by the core laboratory at the David H. Murdock Research Institute (DHMRI). Ribosomal amplicon libraries were generated from collected DNA templates using V6 bacterial 16S universal primers. Each time point consisted of 50 samples for a total of 200, each of which was uniquely indexed, and all samples were pooled together in equimolar proportions and sequenced with 125 bp paired-end reads on an Illumina HiSeq 2500 flow cell. Sequences were demultiplexed and de-barcoded by the DHMRI team prior to delivery.

Sequence Trimming and Quality Filtering
Initial quality filtering and demultiplexing was performed using the integrated Illumina software package, with subsequent additional quality filtering using the Trimmomatic version 0.36 [23]. The resulting forward and reverse .fastq read files with matching mate pairs were merged into a single file per sample using PEAR [24]. The merged .fastq files were then converted to .fasta files using the FastX Toolkit [25] for downstream analyses.

16S Sequencing Analysis
A divisive Amplicon Denoising Algorithm (DADA2) [26] was used for 16S community analysis across all samples. Refinement of Unique Amplicon Sequence Variants (ASVs) relied on DADA2 chimera removal as well as a custom threshold of presence in at least 3 samples (as we had 3 technical replicates for each sample) and an abundance of 10 reads or greater for each ASV based on the rankabundance curve (supplemental Figure S1). Each unique ASV identified by DADA2 was aligned using BLASTN 2.3 to the best match in the SILVA 132 database which had a minimum length of 85%

DNA Extraction
Whole-genome DNA extraction was carried out using DNeasy PowerMax Soil DNA extraction kits (MoBio) according to the manufacturer's instructions. Approximately 10 g of soil or sediment from each sample replicate was loaded into 50 mL extraction tubes and vortexed for 10 m at maximum speed. Soil and sediment were pelleted by centrifugation at 2500× g for 5 m, and the supernatant was passed through the on-column capture filter and washed. The DNA was eluted with 5 mL of sterile water, and quantification was completed using a Qubit fluorometer (Thermo-Fisher) and Nanodrop spectrophotometer (Thermo-Fisher) prior to sequencing library generation. For detailed methods regarding the handling of water samples, see Lambirth et al. [18] 2.3. 16S rRNA Library Preparation and Illumina Sequencing All sequencing was performed by the core laboratory at the David H. Murdock Research Institute (DHMRI). Ribosomal amplicon libraries were generated from collected DNA templates using V6 bacterial 16S universal primers. Each time point consisted of 50 samples for a total of 200, each of which was uniquely indexed, and all samples were pooled together in equimolar proportions and sequenced with 125 bp paired-end reads on an Illumina HiSeq 2500 flow cell. Sequences were demultiplexed and de-barcoded by the DHMRI team prior to delivery.

Sequence Trimming and Quality Filtering
Initial quality filtering and demultiplexing was performed using the integrated Illumina software package, with subsequent additional quality filtering using the Trimmomatic version 0.36 [23]. The resulting forward and reverse .fastq read files with matching mate pairs were merged into a single file per sample using PEAR [24]. The merged .fastq files were then converted to .fasta files using the FastX Toolkit [25] for downstream analyses.

16S Sequencing Analysis
A divisive Amplicon Denoising Algorithm (DADA2) [26] was used for 16S community analysis across all samples. Refinement of Unique Amplicon Sequence Variants (ASVs) relied on DADA2 chimera removal as well as a custom threshold of presence in at least 3 samples (as we had 3 technical replicates for each sample) and an abundance of 10 reads or greater for each ASV based on the rank-abundance curve (supplemental Figure S1). Each unique ASV identified by DADA2 was aligned using BLASTN 2.3 to the best match in the SILVA 132 database which had a minimum length of 85% of the query sequence. These matches were combined in an identity abundance curve to quantify novel ASVs. Taxonomic classification of the ASVs relied on DADA2 built in RDP with SILVA database, version 132 [27] as the training set.

Statistical Methods for Microbial Community Structure
In order to compare the distribution of taxa across pairs of sites, 2 × 2 tables were constructed in which the number of samples containing each taxa for each site pair was summed by a nonzero threshold of reads to signify presence and absence ( Figure S4) [28]. These tables were generated at a taxonomic level from phylum to genus where reads from lower classification levels were merged and normalized by a log frequency calculation [29]. Correction of false positive results from multiple hypothesis testing for all these Fisher tests were performed with the Benjamini-Hochberg method at a 5% False Discovery Rate (FDR) [30].
Quantification of beta diversity differences from Bray Curtis transformations [31] was employed using the Wilcoxon signed test [32] to identify significant coordinates. Beta Diversity from the R package vegan [33] was then visualized using a boxplot of pairwise dissimilarity distances, measuring every combination of inter-sample distances in the transformed matrix space. Additional figures were rendered using R packages ggplot2 [34] and VennDiagram [35].

Phylum and Family Level Community Composition 16S Analysis Reveals Strong Microbial Differences Across the Watershed
An overview of the dominant taxa at the phylum level ( Figure 2) shows patterns of microbial phyla that are consistent with previous literature [36]. While many phyla are common throughout all samples, changes in abundance obviously distinguish soil and sediment from water samples, and to a lesser extent, downstream samples from upstream samples. All samples were dominated at the phylum level by Proteobacteria, as is common in environmental samples [4,36,37] where this large, diverse phylum has many niches. Verrucomicrobia, Nitrospirae, and Acidobacteria were abundant in soil and sediment fractions and were largely absent in water, while Actinobacteria were abundant in the water fraction. Human gut-associated Firmicutes were seen in highest abundance in directly sampled WWTP effluent and in downstream water samples.
At the family level, Principle Coordinate Analysis (PCoA) ordination showed strong separation by sample type (soil, sediment, or water) rather than creek or location in creek ( Figure 3). The unique characteristics of waterborne microbial communities appear to contribute to the most separation along the first principal coordinate. Principal coordinate axis 2 also captures smaller differences between riparian soil and sediment samples, and the expected differences between internal stages of the wastewater treatment process and the stream communities. Univariate permanova tests across all taxonomic levels (Table S1) showed no statistically significant differences between the two stream systems of Mallard Creek and Little Sugar Creek or between the four time points. Multivariate dissimilarity analyses using Bray-Curtis distances at the family level also illustrated the differences between water, sediment, and riparian soil communities ( Figure 4A). The distance between sediment and soil samples was only marginally greater than inter-sediment and inter-soil sample distances ( Figure 4B). At the family level, Principle Coordinate Analysis (PCoA) ordination showed strong separation by sample type (soil, sediment, or water) rather than creek or location in creek ( Figure 3). The unique characteristics of waterborne microbial communities appear to contribute to the most separation along the first principal coordinate. Principal coordinate axis 2 also captures smaller differences between riparian soil and sediment samples, and the expected differences between internal stages of the wastewater treatment process and the stream communities. Univariate permanova tests across all taxonomic levels (Table S1) showed no statistically significant differences between the two stream systems of Mallard Creek and Little Sugar Creek or between the four time points. Multivariate dissimilarity analyses using Bray-Curtis distances at the family level also illustrated the differences between water, sediment, and riparian soil communities ( Figure 4A). The distance between sediment and soil samples was only marginally greater than inter-sediment and inter-soil sample distances ( Figure 4B).     Figure 2 showing dissimilarity among the groups of all samples of the same type (A) and dissimilarity distances between pairs of sample collection sites separated by location (B). In these plots, higher distance scores approaching 1 imply more changes in the microbial community between samples, and distance scores close to 0 indicate more similar microbial communities. Downstream water (DSW), upstream water (UPW), downstream sediment (DSS), upstream sediment (UPS), and effluent which was exclusively water (EFF).

Taxa with <97% Identity to the SILVA Database Indicate Biological Novelty in Our Samples That Are Distributed Differently at Upstream and Downstream Sites
To generate a more detailed view of the microbial community, for each ASV identified by DADA2, we used BLASTN 2.3 to align the ASV sequence to the best match in the SILVA 132 database for the upstream ( Figure 5A Figure 2 showing dissimilarity among the groups of all samples of the same type (A) and dissimilarity distances between pairs of sample collection sites separated by location (B). In these plots, higher distance scores approaching 1 imply more changes in the microbial community between samples, and distance scores close to 0 indicate more similar microbial communities. Downstream water (DSW), upstream water (UPW), downstream sediment (DSS), upstream sediment (UPS), and effluent which was exclusively water (EFF).

Taxa with <97% Identity to the SILVA Database Indicate Biological Novelty in Our Samples That Are Distributed Differently at Upstream and Downstream Sites
To generate a more detailed view of the microbial community, for each ASV identified by DADA2, we used BLASTN 2.3 to align the ASV sequence to the best match in the SILVA 132 database for the upstream ( Figure 5A), downstream ( Figure 5B), WWTP effluent ( Figure 5C), and downstream sediment ( Figure 5D) samples requiring the BLAST alignment to span at least 85% of the query ASV length. Sequences with less than 97% identity to a known sequence would not be placed in a cluster by a closed reference algorithm [26,37]. We therefore considered high abundance/low identity ASVs to be potentially novel if they met this 97% threshold as well as if they had an average read depth of ≥100 sequences (purple symbols, Figure 5).
Using these threshold parameters for novelty and minimal abundance, we found that downstream samples had the highest number of these ASVs, with 59 such sequences detected compared to 47 taxa in upstream samples. In contrast, 26 such taxa were specific to the WWTP effluent, with 2 of these being present at an average sequencing depth of over 1000 sequences (Tables S2-S5).
Compared to the WWTP samples, all stream samples (soil/sediment/water both upstream and downstream) had both more ASVs and more novel ASVs. Using the Fisher test, however, ( Figure S2), we found that there was significant enrichment of potentially novel ASVs in the WWTP compared to both upstream (p = 0.007) and downstream (p = 0.03). In contrast, there was not a significant enrichment of potentially novel ASVs when directly comparing upstream and downstream (p = 0.50) samples. Furthermore, the upstream and downstream samples shared more ASVs than the WWTP effluent did with either upstream or downstream samples ( Figure S3; Tables S6 and S7). Overall, this analysis demonstrates that the WWTP has a reservoir of potentially novel taxa that is distinct from both upstream and downstream samples.
Water 2020, 12, 1557 7 of 16 sediment ( Figure 5D) samples requiring the BLAST alignment to span at least 85% of the query ASV length. Sequences with less than 97% identity to a known sequence would not be placed in a cluster by a closed reference algorithm [26,37]. We therefore considered high abundance/low identity ASVs to be potentially novel if they met this 97% threshold as well as if they had an average read depth of ≥100 sequences (purple symbols, Figure 5). We only considered an ASV as possibly novel if the average read depth was >100 and the percent identity to the SILVA reference database was <97% (purple dots).
Using these threshold parameters for novelty and minimal abundance, we found that downstream samples had the highest number of these ASVs, with 59 such sequences detected compared to 47 taxa in upstream samples. In contrast, 26 such taxa were specific to the WWTP effluent, with 2 of these being present at an average sequencing depth of over 1000 sequences (Tables  S2-S5). We only considered an ASV as possibly novel if the average read depth was >100 and the percent identity to the SILVA reference database was <97% (purple dots).

ASV Analysis Reveals Exclusion of Distinct Taxa in Different Habitats
We subsequently considered all ASVs regardless of their novelty between the different sites and found a larger number of total ASVs in sediment and soil than in the water samples ( Figure 6). analysis demonstrates that the WWTP has a reservoir of potentially novel taxa that is distinct from both upstream and downstream samples.

ASV Analysis Reveals Exclusion of Distinct Taxa in Different Habitats
We subsequently considered all ASVs regardless of their novelty between the different sites and found a larger number of total ASVs in sediment and soil than in the water samples ( Figure 6). Despite large differences in the number of ASVs in each location, the differences in the number of overlapping taxa between water, sediment, and soil were not statistically significant according to multiple Fisher tests comparing the downstream and upstream overlap between sample types.
Individual sequence variants were classified to family using RDP and SILVA 132 as the reference database. To identify individual taxa which were present at one site but completely absent in another, we compared the presence and absence of taxa in all of the upstream samples against all of the downstream samples (Figures 7, S4, and S5). Out of the 180 taxa observed at the family level, three families were differentially present at a 5% FDR threshold after multiple hypothesis correction. Of these, Nocardiaceae and Neisseriaceae were of particular interest due to being significantly present in Despite large differences in the number of ASVs in each location, the differences in the number of overlapping taxa between water, sediment, and soil were not statistically significant according to multiple Fisher tests comparing the downstream and upstream overlap between sample types.
Individual sequence variants were classified to family using RDP and SILVA 132 as the reference database. To identify individual taxa which were present at one site but completely absent in another, we compared the presence and absence of taxa in all of the upstream samples against all of the downstream samples (Figure 7, Figures S4 and S5). Out of the 180 taxa observed at the family level, three families were differentially present at a 5% FDR threshold after multiple hypothesis correction. Of these, Nocardiaceae and Neisseriaceae were of particular interest due to being significantly present in downstream sediment and nonexistent in upstream sediment. This difference between the upstream and downstream sediment is maintained even though Nocardiaceae are abundant among waterborne microbes regardless of the location in the system. The last significant family member, Dongiaceae, has been found to be associated with the roots of wetland plants [38]. In our study, members of this taxa were more abundant in both soil and sediment along the entire waterway, while appearing elevated in water only downstream of the treatment plant.
downstream sediment and nonexistent in upstream sediment. This difference between the upstream and downstream sediment is maintained even though Nocardiaceae are abundant among waterborne microbes regardless of the location in the system. The last significant family member, Dongiaceae, has been found to be associated with the roots of wetland plants [38]. In our study, members of this taxa were more abundant in both soil and sediment along the entire waterway, while appearing elevated in water only downstream of the treatment plant.  Figure S4) and taxa were considered significantly different at a 5% FDR threshold with p values adjusted by Benjamini-Hochberg correction (C,D).

Discussion
Antibiotics at sublethal concentrations are capable of driving adaptation of bacteria and facilitating the development of antibiotic resistance [39]. In a previous study, we showed that treated water released in the Charlotte metro area is driving elevated concentrations of multiple antibiotics in downstream waters [18]. In this study, we expand on our previous work with the impact of that antibiotic-enriched environment on the relatively stationary microbial communities in downstream sediment. While the results from 16S microbial sequencing are interesting, they also give rise to many questions for further study.

Is Urbanization Homogenizing Aquatic Microbes?
Mecklenburg County, where the study sites are located, contains over 3000 miles of streams, most of which lie within the Charlotte metropolitan area and are to some degree impacted by urbanization. All of the assessed mainstem stream reaches in Charlotte are impaired, meaning they are not clean enough for their intended recreational use (Charlotte Stormwater Services, 2020, personal communication). The homogenization of urban landscapes as a result of urban land use change has been the subject of much study in recent years [40][41][42]. The focus of these studies has generally been on animal and plant species, but some recent studies using next generation sequencing approaches indicate that urbanization may affect the biodiversity of urban soils at the microbial community level [43] as well. However, it has been observed in other contexts that urban streams do maintain microbial community differences from location to location [44]. We previously observed  Figure S4) and taxa were considered significantly different at a 5% FDR threshold with p values adjusted by Benjamini-Hochberg correction (C,D).

Discussion
Antibiotics at sublethal concentrations are capable of driving adaptation of bacteria and facilitating the development of antibiotic resistance [39]. In a previous study, we showed that treated water released in the Charlotte metro area is driving elevated concentrations of multiple antibiotics in downstream waters [18]. In this study, we expand on our previous work with the impact of that antibiotic-enriched environment on the relatively stationary microbial communities in downstream sediment. While the results from 16S microbial sequencing are interesting, they also give rise to many questions for further study.

Is Urbanization Homogenizing Aquatic Microbes?
Mecklenburg County, where the study sites are located, contains over 3000 miles of streams, most of which lie within the Charlotte metropolitan area and are to some degree impacted by urbanization. All of the assessed mainstem stream reaches in Charlotte are impaired, meaning they are not clean enough for their intended recreational use (Charlotte Stormwater Services, 2020, personal communication). The homogenization of urban landscapes as a result of urban land use change has been the subject of much study in recent years [40][41][42]. The focus of these studies has generally been on animal and plant species, but some recent studies using next generation sequencing approaches indicate that urbanization may affect the biodiversity of urban soils at the microbial community level [43] as well. However, it has been observed in other contexts that urban streams do maintain microbial community differences from location to location [44]. We previously observed differences among sediment microbial communities across eight urban watersheds (single point, summer only) in Mecklenburg County (Clinton, unpublished data). Those watersheds did not have WWTP input at the point where they were assayed.
The stream systems examined in this study, Mallard Creek and Little Sugar Creek, are adjacent but separate watersheds, and have superficially different patterns of land use, with Mallard Creek being located in a more suburban and more recently urbanized area than Little Sugar Creek. Little Sugar Creek has 37% developed open space, and 11% high intensity development, compared with the 30% developed open space and 4% high intensity development for Mallard Creek [45]. Here, we observe that the main differences in microbial diversity that we see are along the environmental gradient from upstream to downstream of the WWTP, and that differences between the two streams' microbial communities are minimal. While we would require more observations of Piedmont streams to reach a more conclusive understanding of this phenomenon, the upstream-to-downstream shift in the WWTP-adjacent waterway is a much stronger signal than stream-to-stream differences in this study. This suggests that homogenization at the microbiome level may be occurring in Mecklenburg County's urban streams, and that WWTP input may be a factor in homogenization of urban stream sediment communities. Where loss of diversity occurs, it becomes a concern in that it may exacerbate the spread of antimicrobial resistance [46].

Is There Potential for Novel Biodiversity in Urban Stream Environments?
Over the past decade, microbial reference databases such as SILVA have grown rapidly as microbial environments are sequenced and characterized. However, in all of our samples, including upstream samples, we observe numerous ASVs which are less than 97% similar to any known sequence in the SILVA database, but are present at an average sequencing depth of over 100 copies. In the downstream and wastewater treatment plant effluent, we observe dozens of ASVs with a sequence identity below 97% and present in over 1000 copies. Due to their overall abundance and the use of chimera filtering in the analysis, these sequences are unlikely to represent sequencing errors.
It has previously been observed that WWTPs are sources of biological novelty [47], but these low identity/high abundance ASVs exist in downstream and upstream environments as well as inside the WWTP, and the low identity/high abundance components of each community have minimal overlap with each other or with the WWTP. This suggests that the urban streams themselves, and in particular, the downstream environments that are under selection pressure by multiple antibiotics and other discharged compounds, may also be sources of biological novelty. It is interesting to consider the potential role of these taxa. Are they highly adapted to stressed urban stream sediment environments, or to the WWTP downstream environment in particular? Do the downstream-adapted strains carry broad or specific antibiotic resistance elements? What are the implications for the roles of these abundant strains in the downstream sediment microbiome's environmental roles, such as nutrient cycling? We previously observed a suppression of nutrient scavenging and other metabolic functions in downstream waterborne microbial samples using pathway abundance analysis of shotgun metagenomic sequences [18], and while shotgun metagenomic data is not included in the current study, it points to an obvious next step in the investigation of these stressed urban sediment communities.

Are Upstream and Downstream Communities Diverging over Time?
Charlotte's Mallard Creek wastewater treatment facility was built in 1979, and while some form of wastewater treatment facility has existed on Little Sugar Creek since 1929, modern facilities were not put into place until the latter half of the 20th century. Charlotte's downstream urban stream microbes have been under environmental stress from treated wastewater release for 40 years or more, which is ample time for microbial communities to undergo adaptations as a result of exposure [48]. In addition to the potentially uncharacterized taxa discussed above, we also observed considerable non-overlap when presence or absence of taxa in the upstream and downstream environments were compared. Both upstream and downstream sediments had components of their microbial community which were not shared and were also not directly attributable to the WWTP effluent. Of course, there are many conceivable reasons for this, and many possible sources for bacteria in downstream sediments [22], including variations in land use patterns in upstream and downstream locations, resulting in different microbial sources in stormwater runoff.

Are There Indicator Taxa for Downstream Sediment Impact by WWTP Exposure?
One goal of analyzing the sediment microbial communities is to identify sentinel taxa, whose loss or gain in downstream sediments would be an indicator of the degree of impact of WWTP discharge on the stream. Water quality contamination associated with wastewater effluent is often quantified using fecal indicator bacteria such as Escherichia coli and Enterococcus spp. [49][50][51] We compared changes in presence and absence of taxa between upstream and downstream sites at the family and genus level using the Fisher exact test to investigate the impact of the WWTP effluent inputs. Our data provide examples of human-associated taxa that are only present downstream of the WWTP that could be added to other indicator species or used as specific targets to quantify human impacts on downstream environments.
At the family level, we observe that only Nocardiaceae, a family of organisms associated with human respiratory infection and necrosis, especially in immunosuppressed individuals, and Neisseriaceae, a family of gram-negative proteobacteria associated with the mucosal surfaces of animals, follow the pattern of being significantly present in downstream sediment samples while also being absent upstream. Nocardia genus is also involved in the formation of biofilms and foaming in WWTPs [52], and we did observe this genus as a member of the WWTP community in Lambirth et al. [18] Differential abundance of Neisseriaceae in downstream wastewater-exposed environments has been observed previously by Lekunberri et al. (2018) in a shotgun metagenomic study of upstream and downstream sediments from the Ter River in Catalonia, Spain [53], as well as by Chu et al. (2018) in the sediments of Lake Michigan [15] downstream of WWTPs in Manitowoc and Sheboygan WI. These taxa are not necessarily merely indicators of urbanization, as they do not appear in recent studies focused specifically on urbanization [5,54,55].

Conclusions
It is widely accepted that antibiotic overuse and release into the environment via a variety of pathways is highly concerning, as common antibiotics rapidly lose their effectiveness when human pathogens are able to acquire widespread antibiotic resistance. In urban stream sediments downstream of wastewater treatment plants where stationary communities of microbes live under the influence of sublethal concentrations of multiple antibiotics, we observe several phenomena which are of concern for the development of antibiotic resistance: evidence of homogenization of urban stream microbial communities, emergence of previously unclassified taxa specific to downstream sediments, and the presence of specific human-associated taxa only in sediments downstream of the wastewater treatment plant.
There are multiple hypotheses as to why microbial communities vary along environmental gradients, but there is no consensus in the literature as to the relative importance of community turnover, spatio-temporal variability in the dominance of core taxa, and the addition (or loss) of novel taxa [37,[56][57][58][59][60][61][62][63][64][65][66][67][68]. Further study is needed to identify strains emerging under selection in antibiotic-contaminated urban stream environments, to analyze antibiotic resistance-associated gene signatures which may belong to these organisms, and to understand the impact of shifts in the microbial community on the ability of downstream sediments to maintain critical ecological function.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4441/12/6/1557/s1, Figure S1: Average Sequencing Depth vs log10(Rank) of all sequence variants found through the DADA 2 pipeline. The red line indicates an average sequence depth of 10 sequences. This threshold was used to remove all ASVs below this threshold from all other analyses in this manuscript, Figure S2: Fisher 2 × 2 tables that were used to generate p-values comparing the number of Novel and Non-novel ASVs for Downstream vs Upstream, Downstream vs Wastewater Treatment Plant, and Upstream vs Wastewater Treatment Plant, Figure S3: The number of novel ASVs (defined at <97% identity to the SILVA database with an average read depth of ≥100 reads) between downstream and upstream (top panel) and downstream and WWTP, Figure S4: An example of how the Fisher test was used to generate p-values comparing ASVs between two sites. Two by two contingency tables were generated with the columns of each table representing two sites and the rows of each table representing whether a taxa was present or absent at that site. For the example presented, the taxa Nocardiaceae was present in 0 out of 16 samples upstream and 12 out of 16 samples downstream leading to the displayed contingency table. Using the function fisher.test in R, a p-value of 1.61 × 10 −5 was generated for the null hypothesis that samples from both locations had equal prevalence of the taxa. This corresponds to an FDR adjusted p-value of 0.00145 after correcting for all taxa tested at the family level, Figure S5: Using the Fisher test described in Figure S4, the taxa at the order level that were significant for the Downstream vs Upstream Comparisons at an FDR adjusted threshold of 0.05, Table S1: Permanova table of Univariate Models, Table S2: The ASVs that were considered novel for Upstream Samples (purple dots, Figure 2A), Table S3: The ASVs that were considered novel for Downstream Samples (purple dots, Figure 2B), Table S4: The ASVs that were considered novel for the WWTP samples (purple dots, Figure 2C), Table S5: The ASVs that were considered novel for Downstream Sediment Samples (purple dots, Figure 2D), Table S6: List of Overlapping ASVs between Upstream and Downstream from Figure S2 with closest Taxonomic Identification, Table S7: List of Overlapping ASVs between WWTP and Downstream from Figure S2 with closest taxonomic identification. NGS datasets created in this study are deposited in the SRA, accession number PRJNA415974.