Effect of freshwater mussels on the vertical distribution of anaerobic ammonia oxidizers and other nitrogen-transforming microorganisms in upper Mississippi river sediment

Anammox 16S rRNA gene quantification

Microbial culture from a sidestream deammonification process (Hampton Roads Sanitation District, Virginia Beach, VA, USA) served as a source of anammox genetic material for qPCR standard curve construction. PCR products (primers A483f (5′-GTCRGGAGTTADGAAATG-3′) and A684r (5′-ACCAGAAGTTCCACTCTC-3′) (Sonthiphand & Neufeld, 2013)) of the anammox 16S rRNA gene was purified with Qiaquick PCR purification Kit (Qiagen Inc.; Valencia, CA, USA), and cloned into the pCR 2.1-TOPO® vector using the TOPO® TA cloning Kit (Invitrogen Corp.; Carlsbad, CA, USA). Clones were Sanger sequenced at the University of Iowa Institute of Human Genetics with M13F (5′-TGTAAAACGACGGCCAGT-3′) and M13R (5′-CAGGAAACAGCTATGAC-3′) primers to ensure anammox 16S rRNA PCR products were inserted into the vector. Nucleotide sequences were aligned using the Standard Nucleotide Basic Local Alignment Search Tool (Altschul et al., 1997) (GenBank Accession: KU047953) and classified as Candidatus Brocadiales (of the Planctomycetes phylum) with a 95% confidence threshold using RDP Naïve Bayesian rRNA Classifier Version 2.10 (Wang et al., 2007). Plasmid DNA concentration was quantified with Qubit® Fluorometer 1.0 (Thermo Fisher Scientific, Inc.; Waltham, MA, USA), serially diluted, and used to construct qPCR calibration curves.

The anammox 16S rRNA gene from batches 1 and 2 was quantified (Wang et al., 2015) with qPCR using QuantStudio™ 7 Flex Real-Time PCR System (Thermo Fisher Scientific, Inc.; Waltham, MA, USA) with primers A483f and A684r (Sonthiphand & Neufeld, 2013) and analyzed with QuantStudio™ Real-Time PCR Software (Thermo Fisher Scientific, Inc.; Waltham, MA, USA). The threshold cycle (C_t) curves were satisfactory (slope = − 3.374, Y-int = 36.702, R² = 0.998, and amplification efficiency = 97.99%), and PCR product dissociation curves revealed singe peaks centered at a melting temperature of 83°C. The statistical significance of 16S rRNA gene copies was determined via a one-way, repeated measures analysis of variance (ANOVA) (SigmaPlot 13.0, Systat Software, Inc., Chicago, IL, USA) between the 4 treatment groups (n = 20) following a passed normality test (p = 0.826, Shapiro–Wilk) and an equal variance test (p = 0.073, Brown-Forsythe). Pairwise multiple comparison procedures were completed via the Holm-Sidak method with a significance level of 0.050 and a power of 0.990.

Non-targeted amplicon sequencing of the 16S rRNA gene

Batch 2 genomic DNA (20 µL, 1–50 ng/µL) was analyzed by the Argonne National Laboratory, Environmental Sample Preparation and Sequencing Facility (ESPSF) utilizing the Earth Microbiome Project protocol (http://www.earthmicrobiome.org/emp-standard-protocols/16s/). All samples were analyzed together in one batch. The v4 region of prokaryotic 16S rRNA gene (515F-806R) was amplified using the following conditions: 3 min at 94°C, 35 cycles of 94°C for 45 s, 50°C for 60 s, and 72°C for 90 s, followed by 10 min at 72°C (Caporaso et al., 2012). The PCR mixture consisted of 13.0 µL PCR grade water, 10.0 µL 5 PRIME HotMasterMix (Quanta Biosciences, Beverly, MA, USA), 1.0 µL genomic DNA, and 0.5 µL forward and reverse primers (10 µM). 16S rRNA gene amplicon libraries were sequenced by ESPSF using Illumina MiSeq paired end reads (2 × 151 bp) (Caporaso et al., 2012) and uploaded to MG-RAST (ID’s: 4705672.3–4705709.3) and NCBI (BioProject ID PRJNA374585).

Determining the operational taxonomic units (“OTUs”) in each sample from the raw 16S rRNA gene amplicon reads was accomplished using the default Quantitative Insights into Microbial Ecology (QIIME) open-reference pipeline (Navas-Molina et al., 2013). Briefly, the QIIME open-reference pipeline takes paired-end reads as input, which are then joined, demultiplexed, filtered, and clustered into OTUs with uclust (Edgar, 2010). Representative sequences from each cluster were aligned (Caporaso et al., 2010) to GreenGenes 13.5 reference database (DeSantis et al., 2006) with a 97% similarity threshold. RDP classifier (Wang et al., 2007) was used for taxonomy assignment, PyNAST (Caporaso et al., 2010) was used for multiple sequence alignment. Phylogenetic trees were constructed using FastTree2.1.3 with default settings (Price, Dehal & Arkin, 2010). The OTU table from QIIME open reference picking (‘otu_table_mc2_w_tax_no_pynast_failures_json.biom’ in the standard QIIME workflow) was imported into R using the phyloseq package (McMurdie & Holmes, 2013) for downstream analysis, along with the corresponding phylogenetic tree (‘rep_set.tre’) and a metadata mapping file. These datasets were merged to create a single ‘physeq’ object representing the experiment. Alpha-diversity was calculated on the unfiltered OTU abundance data using the Observed species, Chao1 (Chao & Chiu, 2001), and Shannon (Li et al., 2011) metrics. Beta-diversity was calculated using a matrix of Bray-Curtis (Bray & Curtis, 1957) intersample distances and ordination plots calculated with non-metric multidimensional scaling (NMDS). Differential abundance analysis was carried out using the DESeq2 (Love, Huber & Anders, 2014) R package with default settings (test type was “Wald,” fit type was “parametric”). Translating physeq objects into a compatible DESeq2 object was performed with the “phyloseq_to_deseq2” function. The complete data analysis R script can be downloaded from the public GitHub repository: https://github.com/mchimenti/black_chimenti_just_phyloseq/blob/master/phyloseq.r.

Analysis at the OTU level provided a fine scale resolution for significant differences in microbial ecology between mussel and no mussel treatments. To put these results into a biological context, the genus-level OTU file was used to compare relative abundances for N-cycle phylotypes. These groups include AOA genus Candidatus Nitrososphaera, nitrate-damo family “ANME-2d”, NOB genus Nitrospira, anammox genus Candidatus Brocadia, AOB family Nitrosomonadaceae, and nitrite-damo phylum “NC10”. Relative abundance counts for each N-cycle group was tested for statistical significance between treatments, using metadata groups “3 cm with-mussels” (n = 10), “5 cm with-mussels” (n = 10), “3 cm no-mussels” (n = 10), and “5 cm no-mussels” (n = 10). 1-way ANOVA’s of each N-cycle group was performed using the Kruskal–Wallis test (p < 0.05) with Dunn’s multiple correction test (Padj < 0.05) (GraphPad Prism 7.0; La Jolla, CA, USA). Similarly, multiple comparisons were made between all N-cycle phylotype groups and their respective treatments (n = 10); significant differences between relative abundances were tested using the Kruskal–Wallis test (P < 0.0001) and Dunn’s multiple comparison test (Padj < 0.05).

Anammox-targeted 16S rRNA gene quantification

The targeted 16S rRNA gene data from Batch 1 (n = 3, with-mussels) indicated an anammox bacterial gene copy maximum (∼3 ×10⁵ copies g⁻¹ sediment) between 3 cm and 7 cm sediment depth in the presence of mussels (Fig. 2A). The Batch 1 data was normally distributed between 1 cm and 15 cm (Shapiro–Wilk normality test, W-statistic = 0.954, p = 0.773). Only one sediment core went beyond 7 cm leaving anammox bacterial gene copy data at 11 cm and 15 cm without replicates. The Batch 2 data (n = 20 for 3 cm with-mussels, n = 20 for 5 cm with-mussels, n = 20 for 3 cm no-mussels, n = 20 for 5 cm no-mussels) showed that anammox bacteria experienced a 2.2-fold increase (p < 0.001) at 3 cm with-mussels compared to the no-mussels control (Fig. 2B). The anammox gene copies measured at 5 cm were statistically indistinguishable between the with-mussels and no-mussels treatments.

Non-targeted sequencing of the 16S rRNA gene

Summing across all samples, a total of 2,103,661 amplicon sequences were analyzed and about 76,000 unique OTUs were reported by QIIME. Of the unique OTUs, 18,777 had 10 or more reads and 3,916 OTUs had counts exceeding 100 reads. Mussel bed samples had read counts of 45,290 (±15,271) at 3 cm sediment depth and 52,451 (±7,044) at 5 cm sediment depth, while no-mussel samples had 48,920 (±7,517) read counts at 3 cm depth and 63,706 (±25,379) at 5 cm sediment depth (read depths depicted in Fig. S1). The top phyla in mussel bed sediments were Proteobacteria (40.7%), Nitrospirae (35.2%), Chloroflexi (5.9 %), Euryarchaeota (5.0%), Chlorobi (4.2%), and Bacteroidetes (2.3%). Proteobacteria decreased by about 6% with mussels while Nitrospirae increased by 10% with mussels. The most abundant taxonomic families in the Nitrospirae phylum were Thermodesulfovibrionaceae (55%), “FW” (33%), and Nitrospiraceae (13%), and were 5% less, 3% and 2% greater than in no-mussel samples, respectively. With mussels, Proteobacteria taxonomic classes consisted of the following proportions: 68% Deltaproteobacteria (8% less than without-mussels), 16% Gammaproteobacteria, and 15% Betaproteobacteria. A majority of these Deltaproteobacteria OTUs were from “BPC076”, Desulfarculales, and Syntrophobacterales taxanomic orders, while orders Burkholderiales and Xanthomonodales made up a majority of Betaproteobacteria and Gammaproteobacteria taxons.

Species richness was analyzed using three common measures: Observed species, Chao1 and Shannon indices (n = 20 with-mussels and n = 20 without mussels). Together, the three measures indicated a decrease in microbial community richness and evenness in the presence of mussels as compared to sediments without mussels (Fig. 3A). The observed decrease in alpha-diversity reached significance for each of the three measures tested (p = 0.0054 or lower). A similar result was obtained when calculating alpha-diversity measures in samples exclusively from 3 cm (n = 10) or exclusively from 5 cm (n = 10) depths in the presence and absence of mussels. However, the decrease in richness was more pronounced at 5 cm than at 3 cm depth (Figs. S2 and S3).

To compare intersample diversity in species abundances and community composition (“beta diversity”), we employed NMDS scaling to accurately visualize, in 2D space, the higher-order community structure between with-mussels and no-mussels samples (Fig. 3B). The NMDS model produced an excellent representation of the bray–curtis distances for all samples (convergence in 20 iterations, stress ∼ 0.06; shepard plot shown in Fig. S4). The beta diversity clearly differentiated as a function of mussel presence, but not sediment depth (Fig. 3B). Taken together, these data show that mussel presence had a pronounced influence on the microbial community evenness, richness, and composition within the sediment.

Differential abundances in OTUs did not reach significance for metadata values of sediment depth or comparisons between sediment cores. On the other hand, there were numerous differences in OTU abundances when comparing sediment with mussels and without mussels. We performed a differential abundance estimation with the DESeq2 R package using mussel presence status (n = 20 with-mussels, n = 20 no-mussels) as our covariate. 734 OTUs (or 0.94% of the 77,288 OTUs tested) reached significance with a false discovery rate of 0.01. The vast majority of OTUs belonging to the phyla Gemmatimonadetes, Actinobacteria, Acidobacteria, Plantomycetes, Chloroflexi, Firmicutes, Crenarcheota, and Verrucomicrobia decreased by at least 4-fold in the presence of mussels. In contrast, Proteobacteria showed a marked decrease in order Alphaproteobacteria, while showing mixed increasing and decreasing OTUs among Beta-, Delta-, and Gammaproteobacteria. Phylum Nitrospirae also had 38 OTUs which were differentially abundant with p-adj < 0.001. OTUs assigned to the GreenGenes taxonomic family of “0319-6A21” were the most abundant among those OTUs increasing without mussels, while families Thermodesulfovibrionaceae and “FW” were most abundant among those OTUs increasing with mussels.

Many of the Nitrospirae taxons that increased without mussels did so from a smaller average abundance (17 average counts for Nitrospira and up to 126 average counts for Thermodesulfovibrionaceae) relative to those that were increased with mussels (209 average counts for Nitrospira and up to 581 average counts for Thermodesulfovibrionaceae). This explains the 10% increase in Nitrospirae abundance when summing across all samples with mussels. Figure 4 shows the Log2FC categorized by phyla for OTUs with p-adj < 0.0001 (to enhance visual clarity). Significant differences within the Nitrospirae phylum were represented by increases of genus “HB118” in family Thermodesulfovibrionaceae (2.0Log2FC from a mean count of 52, p < 0.001) and unclassified Nitrospira species (0.8Log2FC from an average count of 209, p < 0.001) with mussels. No-mussel treatments showed increases in genus “LCP-6” from family Thermodesulfovibrionaceae (3.6Log2FC from an average count of 126, p < 0.001) and unclassified Nitrospira species (2.1Log2FC from an average count of 17, p < 0.001).

Despite seemingly even representation of phylum Thaumarchaeota between treatments, unclassified species from Candidatus Nitrososphaera were enhanced from an average abundance of 126 (1.73Log2FC, p < 0.001) without mussels, and AOA species, Candidatus Nitrososphaera gargensis increased from an average count of 16 (2.85Log2FC, p < 0.001) without mussels. One OTU classified in the anammox genus, Candidatus Brocadia, increased from an average count of 17 (3.72Log2FC, p < 0.001) without mussels, while another OTU classified as an unknown Candidatus Brocadia species increased from a mean count of 16 (1.2Log2FC, p = 0.001) with-mussels. Furthermore, OTUs belonging to the AOB family Nitrosomonadaceae increased from an average abundance of 6 (1.9Log2FC, p < 0.001) with mussels. Without mussels, taxonomic groups capable of nitrite-damo, phylum “NC10”, increased from average abundances up to 130 (4.4 Log2FC, p < 0.001), and nitrate-damo family “ANME-2d” increased from average abundances up to 59 (3.4 Log2FC, p < 0.001). A summary of Log2FC values for OTUs relevant to N-transformations are listed in Table S1.

N-cycle phylotypes were examined for statistically significant relative abundances between treatments of mussel presence and sediment depth (Table 1). Candidatus Nitrososphaera experienced a 2.6-fold decrease (p = 0.047) with mussels at 5 cm sediment depth. ANME-2d was three times greater (p = 0.049) at 5 cm sediment depth without mussels, compared to 3 cm sediment depth without mussels. Within the mussel bed, Nitrospira were 1.7 times greater (p = 0.0497) at 3 cm depth, and experienced a 1.9-fold increase (p = 0.025) with mussels at 3 cm sediment depth versus control. Candidatus Brocadia was three times greater (p = 0.013) at 5 cm depth without mussels versus 3 cm without mussels, and the 3 cm sediment showed a 2-fold increase (p = 0.002) with mussels versus control. Nitrosomonadaceae was 2.7 times greater (p = 0.015) at 3 cm with mussels versus 5 cm depth with mussels.

Relative abundances of N-cycle phylotypes were compared within each treatment (Figs. 5A, 5B, 5D, 5E) and between treatments (Figs. 5C, 5F, 5G–5I). Within 3 cm sediment samples with mussels (Fig. 5A), Nitrospira was statistically greater in abundance than Candidatus Brocadia, and ANME-2d was less abundant than Nitrospira. Sediment without mussels at 3 cm depth (Fig. 5B) contained statistically greater abundances of Candidatus Nitrososphaera than Candidatus Brocadia, and greater Nitrospira abundances compared to Candidatus Brocadia, Nitrosomonadaceae, and ANME-2d.

Relative abundance comparisons between mussel and no-mussel treatments at 3 cm depth (Fig. 5C) showed that Candidatus Nitrososphaera was reduced in the mussel treatment, while Nitrospira and Candidatus Brocadia were enhanced with mussels. Within mussel sediment samples at 5 cm depth, Nitrospira was more abundant than Candidatus Brocadia, ANME-2d, and Nitrosomonadaceae. (Fig. 5D). On the other hand, Candidatus Nitrososphaera and Nitrospira were both more abundant than Nitrosomonadaceae without mussels at 5 cm sediment depth (Fig. 5E). Comparing microbial communities at 5 cm depth between mussel and no-mussel treatments (Fig. 5F) revealed that Candidatus Nitrososphaera was less abundant with mussels versus the no-mussel population. Nitrospira and Nitrosomonadaceae phylotypes were more prominent with mussels in shallow sediment depths (Fig. 5G).

Overall, Nitrospira made up larger proportions of microbial communities with and without mussels compared to many N-cycle organisms, especially Candidatus Brocadia and Nitrosomonadaceae (Figs. 5C, 5F, and 5G–5I). Without mussels at 3 cm sediment depth, Candidatus Brocadia made up a smaller proportion of the N-cycling microbial community, especially when compared to Candidatus Nitrososphaera, ANME-2d, NC10, and Nitrospira in deeper sediments (Fig. 5H).

N-cycle microbial community

Our research revealed an increase in anammox bacteria abundance 3 cm below the water-sediment interface when mussels were present, shown for the anammox community using anammox-targeted qPCR (2.2-fold increase) and for Candidatus Brocadia using non-targeted 16S rRNA gene amplicon sequencing (2-fold increase). The significance of agreement between these techniques is finding that increases in the genus Candidatus Brocadia are representative for the anammox phylotype as a whole. Candidatus Brocadia may also make up a majority of the anammox community in UMR sediment, as amplicon sequencing did not detect anammox bacteria belonging to other genera. We are confident in these conclusions, as Candidatus Brocadia is often the dominant anammox genus in freshwater sediments (Humbert et al., 2009; Shen et al., 2016; Sonthiphand, Hall & Neufeld, 2014). One study showed that feeding of NH₄⁺, NO₂⁻, NO₃⁻, and acetate led to an 80% enrichment of Candidatus ‘Brocadia fulgida’, signifying that B. fulgida could outcompete anammox species in genera Candidatus Anammoxoglobus and Candidatus Kuenenia, species Candidatus ‘Brocadia anammoxidans’, and even denitrifiers when acetate is present (Kartal et al., 2008). This indicates that Candidatus Brocadia has a distinct ecological niche and can utilize intermediates from anaerobic degradation of organic C to reduce NO₃⁻ (Kartal et al., 2008). Therefore, it is possible that a portion of our observed increases in Candidatus Brocadia with mussels was attributable to C biodeposition in the UMR.

Our research also revealed a vertical distribution of anammox bacteria with higher abundances near the sediment surface, which reflects the vertical distribution found in an agricultural field (Shen et al., 2017), oxygen minimum zone (Galán et al., 2009), flooded paddy fields (Shen et al., 2017; Zhu et al., 2011), and an urban wetland (Shen et al., 2015). A vertical anammox distribution has been shown to coincide with NH₄⁺ presence and NO₂⁻ production (Oshiki, Satoh & Okabe, 2016; Shen, Xu & He, 2014; Shen et al., 2015; Sun et al., 2014) and anammox “hotspots” occur in zones of low, but not entirely absent, O₂ availability (Zhu et al., 2013). Anammox abundance in freshwater sediment can range between 7 × 10⁴ and 8 × 10⁶ gene copies g⁻¹ sediment (Shen et al., 2016), or between 10⁶ and 10⁷ gene copies g⁻¹ sediment in peak NO₂⁻ microniches at the oxic–anoxic interface (Nie et al., 2015; Shen et al., 2015; Zheng et al., 2016). Studies have shown anammox bacteria increase 1.5 to 2-fold within their niche (Nie et al., 2015; Zheng et al., 2016), similar to our findings of a 2.2-fold increase in anammox bacteria 3 cm below the water-sediment interfacewith mussels.

Co-occurrence of aerobic NH₄⁺ oxidation and anammox niches are likely due to linked NO₂⁻ oxidation and reduction, respectively (Shen, Xu & He, 2014). Interestingly, we found that mussels also enhanced taxa from the AOB family Nitrosomonadaceae and the OTUs made up a greater proportion of mussel bed sediment populations near the water-sediment interface. To this point, the pacific oyster (C. gigas) was found to increase porewater NH₄⁺ and elevate the concentration of NH₄⁺ oxidizing microorganisms (Green, Boots & Crowe, 2012). Furthermore, our previous research (Bril et al., 2014; Bril et al., 2017) showed elevated NH₄⁺ and NO₂⁻ in porewater of a similar depth below mussels. It makes sense that these groups of N-transforming bacteria co-occur where their substrate microniches overlap, and is likely enhanced by mussels periodically aerating the sediment (Chen & Gu, 2017). Intermittent aeration has shown to enrich microbial cultures in AOB and anammox bacteria in engineered partial nitritation-anammox processes (Shannon et al., 2015; Yang et al., 2015), and similar to our findings, enriches the anammox genus Candidatus Brocadia (Shannon et al., 2015).

On the other hand, we saw a decrease in Candidatus Nitrososphaera (AOA) with mussels at 3 cm (1.7-fold) and 5 cm (2.6-fold) sediment depths. It makes sense that mussels suppress abundance of AOA since these organisms typically dominate sediment niches with low NH₄⁺ concentrations (Hatzenpichler, 2012; Martens-Habbena et al., 2009). Furthermore, a group of OTUs suppressed by mussels were classified at the species level as Candidatus ‘Nitrososphaera gargensis’, which are partially inhibited by NH₄⁺ concentrations (3.08 mM) much lower than AOB (Hatzenpichler, 2012; Hatzenpichler et al., 2008; Nakagawa & Takahashi, 2015; Pester, Schleper & Wagner, 2011). Furthermore, nitrifier niche partitioning studies using agricultural soil showed that AOB increased in abundance and activity following the addition of urine-derived N, while AOA remained unchanged (Di et al., 2009; Hatzenpichler, 2012; Jia & Conrad, 2009). Therefore, it is possible that mussel biodeposits and an increased flux of agriculturally-fed water into sediment by mussel burrowing enhanced porewater NH₄⁺ composition such that Nitrosomonadaceae out competed Candidatus Nitrososphaera. Our results agree with Chen & Gu (2017), who found bioturbated sediment corresponded with a greater diversity of AOB and lower diversity of AOA microbial communities (Chen & Gu, 2017). On the other hand, our results of decreased abundance of Candidatus Nitrososphaera co-occurring with an increase in Nitrospira is in contrast to previous findings that these organisms may exhibit similar niche partitioning (Pester, Schleper & Wagner, 2011). For example, some species in Candidatus Nitrososphaera can adjust their metabolism for low oxygen availability (Zhalnina et al., 2014) and Nitrospira species are adapted to low oxygen concentrations (Maixner et al., 2006; Nowka, Daims & Spieck, 2015; Zhalnina et al., 2014) and microoxic environments (Schramm et al., 2000). Alternatively, our detected increased abundance of Nitrospira may include species with a variety of environmental niches.

Some Nitrospira species have shown to occupy a niche at oxic–anoxic interfaces, in opposition to NOB with higher O₂ tolerances such as those in genus Nitrobacter (Schramm et al., 2000). This supports our mussel-attributed increases in relative Nitrospira abundances (1.9-fold) at 3 cm sediment depths. Although we saw two different Nitrospira OTUs suppressed and enhanced by mussels, the mussel-enhanced OTUs had a larger mean abundance by about 12%. Different NOB OTUs enhanced with and without mussels further suggests that mussel bed sediments harbor specific NOB strains sensitive to microoxic niches. Despite OTU variability, we can conclude that mussels enhance the Nitrospira phylotype, especially near the water-sediment interface where Nitrospira were 1.7-times greater than deeper mussel bed depths.

On the other hand, we did not expect to see an increase in both NOB and anammox phylotypes due to competition of NO₂⁻ as a substrate. The co-occurrence of Nitrospira and anammox bacteria may be explained by the metabolic versatility of Nitrospira species, especially if mussel-derived urea provided an additional source of NH₃ and NO₂⁻ via reciprocal feeding between ammonia oxidizers and Nitrospira. Furthermore, these phylotypes have been shown to coexist in an oxygen minimum zone, where anammox bacteria obtained a majority of NO₂⁻ from NO₃⁻ reducers (Lam et al., 2009). The similar effect size of mussels on Nitrospira (1.9-fold) and Candidatus Brocadia (2-fold) at 3 cm depth suggests that mussels may exert similar influences on the niches of these phylotypes. It is possible that these anammox and Nitrospira phylotypes were functionally linked in shallow mussel bed sediment, which has been shown for microoxic niches (Van Kessel et al., 2015). Furthermore, it is possible that the Nitrospira co-occurring with Candidatus Brocadia were Nitrospira species with the genetic potential for comammox, as a fluorescence in-situ hybridization study confirmed the extensive aggregation of the 2 phylotypes in hypoxic conditions (<3.1 µM O₂) (Van Kessel et al., 2015). Despite Nitrospira comammox being identified in numerous aquatic environments (Chao et al., 2016; Daims et al., 2015; Pinto et al., 2016), we cannot conclusively identify comammox without sequencing the ammonia monooxygenase gene (Pinto et al., 2016; Van Kessel et al., 2015).

In contrast to studies which found significant N-reduction on both a marine mussel (Mytilus californianus) (Pfister, Meyer & Antonopoulos, 2010) and a freshwater mussel (Limnoperna fortunei) (Zhang, Cui & Huang, 2014), our results showed that mussels suppressed n-damo OTUs in phylum “NC10” (2.1-fold) and family “ANME-2d” (1.8-fold). One study determined NO₃⁻-damo was responsible for NO₃⁻ reduction and anammox for NO₂⁻ reductions in a bioreactor supplied with NH₄⁺, NO₂⁻, NO₃⁻, CH₄, and anoxic conditions, thus concluding anammox outcompeted NO₂⁻-damo (Hu et al., 2015). These findings make sense, because anammox bacteria have a higher affinity for NO₂⁻ (Luesken et al., 2011), anammox outperform n-damo in bioturbated sediments with higher NH₄⁺ and lower NO₂⁻ and NO₃⁻ (Chen & Gu, 2017), and anammox and n-damo communities have a competitive relationship in burrowed mangrove sediment (Chen & Gu, 2017). Furthermore, NC10 bacteria in a peatland were most prevalent at depths with porewater CH₄ concentrations near 300 µM, where NO₃⁻ consumption exceeds production, and in completely anoxic conditions (Zhu et al., 2012). According to the literature, it makes sense that we found UMR mussels enhanced Candidatus Brocadia and suppressed NO₂⁻ reducing-NC10. Perhaps n-damo organisms did not have a favorable niche in mussel bed sediment because biodeposition products created an excess of NH₄⁺ in sediment porewater (Winkler et al., 2015), or burrowing activity increased oxygen concentrations and made methane oxidation unfavorable (Van Bodegom et al., 2001). Our finding that no-mussel sediment contained three times more ANME-2d in deeper, and presumably anoxic sediment, further suggests that mussels broaden the oxic–anoxic interface niche (Chen & Gu, 2017; Luesken et al., 2011). However, we cannot extrapolate these findings to all denitrifying organisms, since denitrifying species are sporadically distributed among various taxonomic lineages, and are difficult to identify solely with16S rRNA amplicon sequencing (Ishii et al., 2011).

Although we observed greater relative abundances of Nitrospira than Candidatus Brocadia in a majority of treatments, both phylotypes increased by a factor of 2 with mussels at 3 cm depth. No-mussel samples contained a significantly smaller proportion of Candidatus Brocadia in shallow sediments compared to almost all N-transformers found in the deeper control sediments. Our phylotype-level analyses revealed similarities with the OTU-level differential abundance comparisons. For example, phylotype comparisons showed ANME-2d was less abundant than Nitrospira in 3 cm sediments with mussels, and Candidatus Nitrososphaera was more abundant than Candidatus Brocadia in 3 cm sediment samples without mussels. These results relate to DESeq2 OTU comparisons which found Candidatus Brocadia and Nitrospira enhanced with mussels while ANME-2d and Candidatus Nitrososphaera were suppressed with mussels.

Extending our focus beyond N-cycling organisms, we demonstrated that mussels promoted a large effect size for OTUs classified as Thermodesulfovibrionaceae (Nitrospirales order). In contrast to Nitrospira, the Nitrospirales genus Thermodesulfovibrio contains multiple sulfate reducing species (Kirchman, 2012; Sekiguchi et al., 2008) and can outcompete other anaerobic organisms when sulfate is present (He et al., 2015). These findings are corroborated by discoveries of significantly greater C and sulfate concentrations from mussel biodeposits and 63% greater sulfate reduction in sediments with mussels (McKindsey et al., 2011). Biodeposition products often lead to increasingly anoxic sediment and greater activity of anoxic microorganisms (Kellogg et al., 2013; McKindsey et al., 2011), presumably due to consumption of excretion products by oxygen-consuming microorganisms (Pollet et al., 2015). Interestingly, Fdz-Polanco et al. (2001) observed simultaneous N and sulfate removal in an anaerobic fluidized-bed reactor and proposed simultaneous anammox and sulfate reduction. Coupled biological sulfate reduction and anammox reactions are metabolically feasible (Schrum et al., 2009; Strous et al., 2002) and have been of interest in the recent history (Ali et al., 2013; Cai, Jiang & Zheng, 2010; Rikmann et al., 2016; Rios-Del Toro & Cervantes, 2016), therefore warranting further research. Therefore, we showed that Thermodesulfovibrionaceae are significantly increased in the presence of mussels which may affect sulfate reduction (Mahmoudi et al., 2015) in tandem with anammox reactions in UMR sediments.

As a whole, mussels do have an impact on microbial niches and lower the overall community diversity. Mussel-influenced changes in microbiological diversity may have larger ecosystem implications, such as macrobiota richness and diversity (Arribas et al., 2014; Borthagaray & Carranza, 2007). Native freshwater mussels are capable of increasing macrobiota diversity as a result of being keystone species (Hartmann et al., 2016) and ecosystem engineers (Chowdhury, Zieritz & Aldridge, 2016; Lopes-Lima et al., 2014). Mussel biogeochemical hotspots can lead to a bottom-up trophic cascade by enhancing N substrates normally limiting primary productivity, ultimately leading to increased richness (Atkinson et al., 2013) and biodiversity (Allen et al., 2012) of higher trophic levels.