1 INTRODUCTION

Nitrogen (N) pollution in coastal ecosystems due to excessive anthropogenic N inputs has become a serious environmental issue on regional and global scales, which leads to eutrophication and associated deleterious ecological changes. These changes included hypoxia and anoxia (Camargo and Alonso, 2006), increased harmful algal blooms (Anderson et al., 2002), alteration of community structure (Bürgi and Stadelmann, 2002), and loss in biotic diversity (Bürgi and Stadelmann, 2002). Therefore, increasing concerns regarding pathways for N loss in coastal environments have been raised for decades (Galloway et al., 2008).

Microbial mediated denitrification and anammox are two major pathways of N removal in marine environments (Ward, 2013). Denitrification reduces nitrate (NO₃^-) sequentially to dinitrogen gas (N₂) coupling to oxidizing organic matters, while anammox combines ammonium (NH₄⁺) and nitrite (NO₂^-) to yield N₂. These two pathways account for about 70% of fixed N loss in the marine N cycle (Codispoti, 2007; Ward, 2013), whereas their contributions vary widely over space and time (Hietanen and Kuparinen, 2008; Brin et al., 2014). Considering the importance of denitrification and anammox for nitrogen removal, it is critical to understand the dynamics and distributions of the relative functional microbes in coastal ecosystems.

Denitrification is performed by a diverse assemblage of microorganisms, during which different types of metabolic enzymes are produced (Zumft, 1997). The reduction of N₂O to N₂ is catalyzed by nitrous oxide reductases (Nos), and this is an important step in the denitrification process because greenhouse gas N₂O is converted into N₂ and complete denitrification is performed during this step (Zumft, 1997). Therefore, the nosZ gene is usually used as a biomarker to study the ecological behavior of denitrifying microorganisms in coastal environments (Scala and Kerkhof, 1999; Magalhães et al., 2008; Wyman et al., 2013; Wang et al., 2014; Yang et al., 2015).

The diversity of anammox bacteria has been explored using specific 16S rRNA, hzo (hydrazine oxidoreductase), Annirs (anammox nitrite reductase), and hzsA (hydrazine synthase) genes as molecular markers (Li et al., 2010, 2011b; Hou et al., 2013; Bale et al., 2014; Shehzad et al., 2016). The known anammox bacteria are affiliated to the order Candidatus Brocadiales within the phylum Planctomycetes and include five candidate genera: Ca. Brocadia, Ca. Kuenenia, Ca. Anammoxoglobus, Ca. Jettenia, and Ca. Scalindua (Schmid et al., 2003; Kartal et al., 2007, 2008; Humbert et al., 2010). Scalindua typically dominates in marine settings (Schmid et al., 2007; Woebken et al., 2008), while non-Scalindua mainly appeared in freshwater, reactors, estuarine and coastal environments (Dale et al., 2009; Dang et al., 2013; Hou et al., 2013). Despite extensive investigations of genetic diversity and community composition of denitrifiers and anammox bacteria in various coastal habitats, little is known about their spatial and seasonal patterns in margin basins, where the spatiotemporal heterogeneity in hydrology, sedimentary characteristics, and anthropogenic influences determined complex compositions of denitrifiers and anammox bacteria.

Various environmental factors have been suggested to affect distributions of denitrifying and anammox bacterial communities, including availability of nitrogen, temperature, oxygen, trace metals, salinity and organic matters (Dang et al., 2010; Hou et al., 2013; Babbin et al., 2014; Zhang et al., 2014). Recently, Lipsewers et al. (2016) found that seasonal hypoxia and elevated sulfide concentration in coastal bottom waters impacted distributions of denitrifiers and anammox bacteria. In the Mai Po nature reserve, a higher diversity of anammox bacteria was observed during summer due to numerous anthropogenic and terrestrial inputs bringing in Kuenenia (Li et al., 2011a). However, the seasonality and environment drivers of these functional groups have not been well illustrated in coastal ecosystems.

The Bohai Sea (BS) is the innermost basin of China, with an average depth of 18 m and a very long water exchange half-life of 17 to 21 months (Wei et al., 2002). The BS and its coast are known as a "golden necklace" in North China. Because of rapid developments, the BS environment is affected increasingly by human activities, especially inorganic N inputs. It receives roughly 2.5×10⁴ t per year of dissolved inorganic nitrogen (DIN) from more than 40 tributary rivers, mainly the Huanghe (Yellow) River (SOA, 2016), of which the fluxes in the flood season (July–September) account for 70%–80%. Up to 1/3 area of the BS was eutrophicated in 2015, and the average seawater N:P ratio reached 67:1 (SOA, 2016), was significantly higher than the Redfield Ratio (16:1) of clean seawaters (Redfield, 1958). Heavy eutrophication mostly occurred in the coastal regions and estuaries (Wang et al., 2009). From 2015 to 2017, more than 20 algal blooms occurred in the BS during May to September, with impact on an area of over 1 500 km² (SOA, 2015–2017). Through the narrow Bohai Strait, the BS connects with the outer basin North Yellow Sea (NYS). Relative to the BS, the NYS is opener and cleaner, with an average depth of 40 m. The most sediments of the two basins are composed of clayey silt and silt sediments and deeply influenced by sediment inputs from the Huanghe River (Qiao et al., 2017). The hydrographic conditions of the two basins are governed by coastal currents and the Yellow Sea Warm Current (YSWC). The YSWC transports warm and saline waters into the BS through the NYS, which prevails in winter and weakens in summer (Xu et al., 2009).

In the present study, we focus on (1) the heterogeneities in diversity and composition of denitrifiers and anammox bacteria in basins (BS and NYS), seasons (May and November), and regions (coastal and central), and (2) the key environmental factors affecting the distribution of denitrifiers and anammox bacteria. The results highlight the comparative ecological roles of the two functional microbes in complex and heterogeneous coastal environments.

2 MATERIAL AND METHOD 2.1 Sampling and physicochemical analysis

Four sampling sites (B66, B41, BF01, and B24) in the BS and NYS (Fig. 1) were selected during the R/V Dong-Fang-Hong 2 cruises on May 2–24 and November 1–20, 2012. Sediments were box-cored and three random replicated surface sediments (top 0–5 cm) were collected at each site during each cruise. At last, 24 sediment samples were obtained. The sediment samples were homogenized in sterile plastic bags, and aliquots were put into cryo-vials and stored immediately in liquid nitrogen for subsequent DNA extraction and physicochemical analysis.

Physicochemical parameters such as salinity (Sal), temperature (Temp), and water depth were estimated in situ by Seabird 911 Conductivity-Temperature-Depth (CTD). Sediment pore waters were obtained by centrifuging at 6 000 r/min, and NO₃^-, NO₂^-, NH₄⁺ and PO₄^3- in pore waters were determined with a nutrient AutoAnalyser (Seal, Germany). The geographic distance from the sampling site to the coast (offshore distance, L-dist), bottom water Chlorophyll a (Chl-a), dissolved oxygen (DO) and pH, as well as sediment grain size (GS), total organic carbon (TOC) and nitrogen (TN), and trace metals (Pb, Cr, Mn, Fe, Co, Ni, Cu, Zn, As, and Cd) was measured according to methods described previously by Zhang et al. (2018).

2.2 DNA extraction, amplification, and cloning

Approximate 0.5 g of sediment was used to extract total genomic DNA by the FastDNA spin kit for soil (MP Biomedical, USA) according to the manufacturer's protocol with slight modifications. The concentration of extracted DNA from sediment samples was measured by a NanoDrop 2000 Spectrophotometer (Thermo Scientific, Wilmington, DE, USA).

The nosZ gene was amplified using the primers nosZF/nosZR (Throbäck et al., 2004) to generate about 700-bp fragments. Anammox 16S rRNA gene was used to amplify the desired gene fragments (about 477 bp) by a nested PCR technique described previously (Hou et al., 2013) using the primers PLA46f/1390r-AMX368f/820r. The obtained products were separated by electrophoresis on a 1% agarose gel and purified using the Agarose Gel DNA Recovery Kit (Tiangen, Beijing, China). Three replicates of each site were mixed well, ligated into the pTZ57R/T vector (Thermo, USA) and transformed into Escherichia coli DH5α competent cells (Tiangen, China). Finally, 8 clone libraries were constructed for nosZ and anammox 16S rRNA gene, respectively. The positive clones were carefully selected using X-Gal-IPTG LB indicator plates supplemented with 100-μg/mL ampicillin, and re-amplified using M13F/R primers.

2.3 Sequencing and phylogenetic analysis

Approximate 100 positive clones were randomly selected from each gene library for sequencing (Sangon Biotech, Shanghai, China). Analyses of nosZ sequences were carried out by translating into amino acid sequences using the BioEdit software (Hall, 2011). The nosZ amino acid sequences and anammox 16S rRNA gene sequences were aligned by BioEdit and grouped into operational taxonomic units (OTUs) at 80% and 95% identity, respectively, using the DOTUR program (Schloss and Handelsman, 2005). The closest matches of each OTU identified by BLASTn were retrieved from the GenBank. The phylogenetic trees were constructed by the maximum likelihood method using the program RAxML 8.0 (Stamatakis, 2014). The optimum substitution models for nosZ and anammox 16S rRNA were determined by the program ProtTest 2.4 (Abascal et al., 2005) and Modeltest 3.7 (Posada and Crandall, 1998), respectively, and WAG+I+G and GTR+I+G models were the best fits for nosZ and anammox 16S rRNA, separately. We determined the confidence in tree topology using a bootstrap analysis with 1 000 restarts. Sequences were deposited in GenBank under the accession numbers: MH152718–MH153503 (nosZ gene) and MH121704–MH122515 (anammox 16S rRNA gene).

2.4 Statistical Analysis

Coverage of each gene library was calculated as C=[1–(n/N)]×100, where n is the number of unique OTUs and N is the total number of clones in a library. By using DOTUR, the numbers of observed OTUs and alpha diversity indices (Shannon H, Simpson 1/D and Evenness J) were calculated for each gene library.

To compare the spatiotemporal heterogeneities in diversity and composition of denitrifiers and anammox bacteria, all samples were divided into BS and NYS groups based on the geographic separation between the two basins, the May and November groups based on seasonality, as well as coastal and central groups based on anthropogenic disturbance intensity. Student's t test was conducted to test the differences in environmental factors and alpha diversity estimators between groups. Pearson's correlations were conducted to investigate the associations between alpha diversity estimators or relative abundance of a specific cluster and environmental parameters. All of the analyses were conducted using SPSS v.19.0 (Chicago, USA). To visualize the differences of denitrifying and anammox bacterial community in all samples, nonmetric multidimensional scaling (NMDS) was conducted based on a Bray-Curtis similarity matrix using the PRIMER (v.6) software package (Primer-E, UK), and ANOSIM was used to test pairwise community structure differences among groupings of samples. To investigate the relationship between environmental parameters and denitrifying or anammox bacterial assemblages, CCA or RDA analysis was performed in R v.3.4.3 with the vegan package.

3 RESULT 3.1 Environmental factors of the BS and NYS

The environmental conditions of the studied area have been described in a recent study of Zhang et al. (2018). In brief, the water depth of our sampling sites was distinctly deeper in the NYS (42.4±5.10 m) than in the BS (19.4±4.05 m) (P=0.01; Table S1). Influenced by freshwater and nutrient discharges from the Huanghe River, the site B66 held the lowest salinity and highest Chl-a in bottom waters. Porewater nutrients showed significant seasonal and regional variations (P≤0.05; Table S1). Higher NO₃^- and NO₂^- concentrations were observed in May, especially at the coastal sites BF01 and B66, but on the contrary, the NH₄⁺ concentration was higher in November, particularly in the BS. The concentration of PO₄^3- was significantly higher in the central regions (5.54±1.70 μmol/L) than in the coastal regions (1.15±0.65 μmol/L) (P=0.05). However, N:P (the molar ratio of DIN to PO₄^3-) showed no any spatial or seasonal difference. No significant difference in sediment grain size (GS) was observed across seasons, basins, and regions, but the finest sediment was found at the site B66. Sediment TOC% and TN% were significantly higher in the NYS than those in the BS (P≤0.03), but C:N (the mass ratio of the total organic carbon to total nitrogen), in the range of 2.18 to 5.54, was similar in the two basins. Trace metals in sediments showed strong seasonal trends (P < 0.01) with higher values in November, except for metal Mn, which was much higher in the BS compared with the NYS (P < 0.01) (Table S1).

3.2 Spatiotemporal variations in the alpha diversity of denitrifiers and anammox bacteria

The obtained 786 nosZ sequences were 30.5%– 98.6% identical with each other and 77.4%–98.7% identical to the top-hit GenBank sequences at the amino acid scale. The obtained 812 anammox 16S rRNA sequences were 68.6%–99.7% identical with each other and 88.1%–100% identical to the top-hit GenBank sequences at the DNA scale. Finally, 112 nosZ and 51 anammox 16S rRNA OTUs were identified by DOTUR and the coverages of all libraries ranged from 76.6% to 100.0% (Table 1), indicating that the majority of the nosZ and anammox 16S rRNA sequence types was captured.

Shannon, Simpson and Evenness indices showed that the greatest nosZ gene biodiversity was observed at the site B66 near the Huanghe River estuary, while the lowest value occurred at the central NYS site B24, especially in November (Table 1). Generally, the diversity of the nosZ gene was significantly higher at coastal sites than at central sites (P≤0.01). However, no obvious seasonal and basinal difference in the nosZ gene diversity was observed (P > 0.05) (Table 1).

Overall, the diversity of anammox bacteria was much lower than that of nosZ-denitrifiers. The greatest anammox 16S rRNA gene biodiversity was also found at the site B66 during November, while central sites B41 and B24 held the lowest anammox 16S rRNA gene diversities, where only 1 or 2 OTUs were found (Table 1).

3.3 Community compositions of nosZ- denitrifiers and anammox bacteria

A wide range of nosZ sequence divergence was observed in the BS and NYS sediments. The maximum likelihood phylogenetic tree showed that all nosZ sequences fell into 10 (Ⅰ–Ⅹ) clusters (Fig. 2). These sequences were closely affiliated with other environmental nosZ clones retrieved from the South China Sea sediments (Yu et al., 2018), Arable Land coastal marine sediments, New Jersey marine sediments (Scala and Kerkhof, 1998), Laizhou Bay sediments, wastewater treatments, and soils. However, the sequences from cluster Ⅴ, the most numerically dominant genotype identified (44.1% of all nosZ clones), exhibited low overlap with known denitrifiers in the database, suggesting that they might be unique to the BS and NYS sediments (Fig. 2).

All nosZ sequences putatively derived from Alpha-, Beta-, and Gamaproteobacteria, and Alphaproteobacteria-related sequences were highly represented in our samples, accounting for 78.5% of all nosZ sequences detected. The most abundant OTU, E119 (19.7% of all clones), with a sequence identity of 89% to the nosZ amino acid sequence of Nitratireductor indicus, was frequently detected at the central sites B24 and B41 (Fig. 2). The second abundant OTU, H92 (10.3% of all clones), with a sequence identity of 92% to the nosZ amino acid sequence of Ruegeria pomeroyi, represented higher relative abundance in May than in November (Fig. 2). However, Betaproteobacteria and Gamaproteobacteria-related sequences were exclusively retrieved from coastal sites BF01 and B66 (Fig. 2), and among Betaproteobacteria, the dominant OTU, C37, with 90% similarity to the nosZ amino acid sequence of Thiobacillus denitrificans, occurred only at the site B66 (Fig. 2).

Three known anammox bacterial genera including Ca. Scalindua, Ca. Brocadia and Ca. Kuenenia were detected in the studied area (Fig. 3). Ca. Scalindua was predominant (78.3% of all clones) in the anammox bacterial libraries, and five (SI–SV) distinctive Scalindua clusters were identified. Clusters SI, SIII, and SV are affiliated with Scalindua marina with 96.4%–99.8% sequence identity, and clusters SII and SIV show 98.5%–99.4% identity to the sequences of Scalindua wagneri. The most abundant OTU, H100, has an identity of 99.8% to the sequence of Scalindua marina from the sediments of the Changjiang (Yangtze) River estuary (Hou et al., 2013). Most Brocadia clusters belonged to Brocadia fulgida (98.3%–98.7% identity), including BI, BII, BIV, and BV, while cluster BIII shows 95.8% identity to the sequence of Brocadia anammoxidans. The dominant Brocadia OTU, C1, and D74 has an identity of 100% to the sequences of anammox clones in the suspended sediments of Huanghe (Yellow) River, and 99.9% to the sequences of Brocadia fulgida in the sediments of the Changjiang River estuary (Hou et al., 2013). Altogether 34 sequences are affiliated to Kuenenia, with 96.2%–98.9% identity to the sequences of Kuenenia stuttgartiensis. In addition, the phylogenetic analysis showed a potential anammox cluster, sharing less than 93% similarity with sequences from all the other clades, which only occurred at the site B66.

3.4 Distributions of nosZ- denitrifiers and anammox bacteria

The NMDS plots (Fig. 4) showed that the whole nosZ-denitrifier community was divided into two distinctive groups and the assemblages of B66 were clearly separated from those of other sites. Similarly, the anammox bacterial assemblages of B66S, B66W, and BF01S were distinctive from those of other samples. The ANOSIM results (Table 2) verified that the seasonal and basin-wise differences of both nosZ-denitrifiers and anammox bacteria were not significant (P > 0.05), but the regional difference was distinct (for nosZ, P=0.057 and for anammox 16S rRNA, P=0.029).

The distributions of specific clusters of nosZ-denitrifiers and anammox bacteria identified in the phylogenetic analysis are demonstrated in the heatmap plots (Fig. 5). The nosZ sequences from cluster Ⅴ are shared in all samples and occur in the central sites more frequently than at coastal sites. In addition, Cluster Ⅰ also presents in all samples and less in B66. Cluster Ⅶ occurred most at B24 in November. On the contrary, sequences from Clusters Ⅱ, Ⅵ, Ⅸ, and Ⅹ present higher relative abundance at coastal sites than at central sites (Fig. 5a). As for anammox bacteria (Fig. 5b), different Scalindua clusters have their respective niches. Sequences from SIII were primarily retrieved from B24 and B41, representing 75.5% and 95.2% of the sequences in the two sites, respectively. Cluster SII was found prevalent in the sample BF01S, while Cluster SI is prevalent in the sample BF01W. However, Cluster SIV presents frequently in samples B66S and B24S. The Brocadia cluster BV is highly represented in the sample B66W, while Cluster BIV mainly presents in the sample B66S. The Kuenenia clusters are restricted to the site B66.

3.5 Factors driving variations in nosZ- denitrifiers and anammox bacterial community diversity and structure

All alpha diversity estimators of nosZ-denitrifiers negatively correlated with the offshore distance (R≤ -0.88, P≤0.004), depth (R≤-0.74, P < 0.04), and PO₄^3- (R≤-0.70, P < 0.05), and positively with the ratio of N:P (R≥0.72, P≤0.04; Table 3). However, no any environmental factor showed significant correlations with the alpha diversity of anammox bacteria (P > 0.05; Table 3).

The whole nosZ-denitrifier community structure was significantly co-varied with bottom water salinity (P=0.005), DO (P=0.015), and the concentration of Pb in sediment (P=0.005) (Fig. 6a). These factors explained 84.1% of the total variance of the nosZ-denitrifier community-environment relationship. Nevertheless, only salinity (P=0.002) was identified as the significant environmental factor in correlation with the variation of the whole anammox bacterial community structure and spatial distribution (Fig. 6b), providing 35.1% of the total variance of the anammox bacterial community-environment relationship.

For cluster-specific correlation with environmental factors, the nosZ cluster Ⅴ responded positively to elevated offshore distance, depth, salinity, and PO₄^3- (R≥0.77, P≥0.02); Cluster Ⅰ positively correlated with DO (R=0.79, P=0.02) and negatively correlated with temperature (R=-0.80, P=0.02). Nevertheless, Cluster Ⅱ showed a significant correlation with NO₃^- and NO₂^- (R≥0.81, P < 0.02). Higher sedimentary TOC% and TN% seemed to favor the relative abundance of Cluster Ⅶ (R≥0.80, P < 0.02), while Cluster Ⅹ decreased with increasing of salinity, depth and sediment grain size (R≤-0.71, P < 0.05). Within anammox clusters, Cluster SIII was highly correlated with the offshore distance (R=0.91, P=0.002), Cluster SII was most strongly correlated with NO₃^- and NO₂^- (R≥0.88, P≤0.004), and Cluster SIV showed high correlation with the concentration of metal Pb (R=0.82, P=0.01). Salinity seemed to be the most important factor affecting the relative abundance of Brocadia and Kuenenia clusters (R≤-0.74, P≤0.03). Apart from that, Clusters BIV and KIII had a high correlation with the concentration of Chl-a in bottom waters (R≥0.71, P < 0.05) and high levels of tolerance to metal As in sediments (R≥0.85, P≤0.007) (Table 4).

4 DISCUSSION 4.1 Distribution and environmental drivers of benthic nosZ-denitrifiers

Sequencing of nosZ clones revealed 112 OTUs at 80% amino acid identity in the BS and NYS sediments. This inherent diversity of nosZ gene was much higher than those in other coastal and marine environments, such as Douro River estuary (Magalhães et al., 2008), Goa mangrove forest (Fernandes et al., 2012), and Pacific Ocean (Scala and Kerkhof, 1999, 2000). The result suggested that the anthropogenic perturbation-dominated setting of the BS and NYS might promote the biodiversity of the nosZ gene to adapt to the complex environments. Another possible reason was the PCR bias of the primer set nosZF/nosZR, which would likely recover more diverse environmental nosZ phylotypes than the primer sets (e.g. nosZ661F/ nosZ1773R and nosZ1211F/nosZ1897R) used in other coastal environments, although these primers captured a similar range of cultrue-based phylotypes (Throbäck et al., 2004). This possibility should be verified in further studies.

Although certain nosZ sequences appeared to be ubiquitous in coastal and marine sediments, for example, some of our sequences were closely affiliated with those in the South China Sea sediments (Yu et al., 2018), Arable Land coastal marine sediments, New Jersey continental shelf sediments (Scala and Kerkhof, 1998), and Laizhou Bay sediments, a big nosZ group (including more than 40% of all nosZ clones) was restricted in the BS and NYS setting (Fig. 2), suggesting that nosZ-denitrifiers could evolve inhabit-specifically not geographically. In addition, the alpha diversity and composition of our nosZ-denitrifiers presented a clear heterogeneity between coastal and central sites but no significant difference between seasons and basins (Table 1 and Table 2). Several studies have monitored the distribution of denitrifiers in marine sediments. Magalhães et al. (2008) found a similar result to us, that the composition of nosZ assemblages in the Douro River estuary sediments showed a site-specific difference but was stable over time. On the contrary, Scala and Kerkhof (2000) found that geographic distance (centimeters to kilometers) had a major influence on the structure of nosZ-denitrifiers in continental shelf sediments. In the Pacific coast of Mexico, Liu et al. (2003) observed that the distribution of nirS and nirK-denitrifiers was controlled by geographic location and biogeochemical conditions. Collectively, local environmental properties and geographic separation could have a combined effect on structuring marine denitrifier communities. The coastal sites, especially near the Huanghe River estuary, held higher diversity and distinct community structure of nosZ gene from central sites, probably because freshwater and sediment inputs from the river brought in terrigenous nosZ phylotypes, and on the other hand, some specific nosZ phylotypes could be evolved in the anthropogenic perturbation-dominated coastal setting.

Similar to previous studies of coastal and marine nosZ-denitrifiers (Magalhães et al., 2008; Mills et al., 2008; Fernandes et al., 2012; Yang et al., 2015; Yu et al., 2018), all nosZ sequences detected are related to proteobacteria, and Alphaproteobacteria formed the most dominant and ubiquitous putative nosZ group in the BS and NYS sediments, suggesting that this group was well adapted to coastal marine sediments and likely contributed importantly to nitrogen removal in the environments. Within Alphaproteobacteria, N. indicus-related (89% similarity) nosZ sequences primarily occurred in the central sites B24 and B41. In addition to denitrification (Labbé et al., 2004), some members of Nitratireductor spp. have been reported to be capacity of degrading crude oil, hydrocarbon or complex organic matters in marine environments (Lai et al., 2011a, b; El Hanafy et al., 2016). Site B41 located nearby the drilling platform Penglai 19-3, where a large oil spill incident occurred in November 2011 (Pan et al., 2015), resulting in serious crude-oil pollution in the sea. Site B24 located in the center of the NYS, in which large amounts of old and recalcitrant organic matters with high C:N ratio (0.73%–0.76%) existed, owing to the hydrodynamic forces constrained by cyclonic circulations (Hu et al., 2016). Responding to these complex organic matter compositions in sediments, N. indicus may be selected and enriched in these sites. R. pomeroyi-related (92% similarity) nosZ sequences were mainly detected during May. R. pomeroyi was a model marine Roseobacter bacterium, with the capacity not only to denitrify (Wyman et al., 2013; El Hanafy et al., 2016) but also to demethylate dimethylsulfoniopropionate (DMSP) (Gonzalez-Silva et al., 2017), an important compatible solute of marine algae. This organism often dominates in the phycosphere microenvironment and utilizes DMSP as a carbon source (Rink et al., 2007; Goecke et al., 2013). In the BS and NYS, algal blooms generally peak during May to July (Wei et al., 2004), and thus the organism is prevalent in that season. In contrast to other works (Magalhães et al., 2008; Mills et al., 2008; Fernandes et al., 2012; Yang et al., 2015; Yu et al., 2018), we observed that the Betaproteobacteria-related nosZ group comprised a substantial fraction of the overall nosZ community detected in the Huanghe River estuary site. Betaproteobacteria is commonly detected in the freshwater and estuarine systems, such as Columbia River estuary (Crump et al., 1999), Changjiang River estuary (Feng et al., 2009), and Zhujiang (Pearl) River estuary (Liu et al., 2015). In this study, most Betaproteobacteria-related nosZ sequences were similar to the terrigenous nosZ clones retrieved from wastewater treatments, lake waters and soils (Fig. 2), suggesting that Huanghe (Yellow) River diluted water contributed profoundly to the composition of nosZ-denitrifier in the region. The most abundant Betaproteobacteria nosZ phylotypes were closely related (about 90% similarity) to a known autotrophic denitrifier Thiobacillus denitrificans. Thiobacillus denitrificans is also a well-known sulfur-oxidizing bacterium (Shao et al., 2010) and its predominance may suggest the frequent occurrence of sulfur oxidization process in the coastal region.

Based on multiple analyses, salinity had the most significant impact on nosZ-denitrifier diversity, community structure, and distribution. Separation of nosZ genes according to salinity was also observed in estuarine sediments along a salinity gradient (Magalhães et al., 2008; Yang et al., 2015). Additionally, at the global scale, salinity was the major driver of nirS and nirK-denitrifier communities in aquatic environments (Jones and Hallin, 2010). The ecological mechanism of how salinity governs denitrifier community remains unclear. A study in the Douro estuarine water column, in which NO₃^- concentrations co-varied with salinity (Magalhäes et al., 2005), and previous NO₃^- and salt addition experiments (Magalhäes et al., 2005) revealed that denitrification rates could be a function of NO₃^- availability and that salinity does not have a direct effect on the denitrification process. In this study, the nosZ cluster Ⅴ (Alphaproteobacteria) preferred high salinity, whereas Cluster Ⅹ (Betaproteobacteria) was exclusively detected in the low-salinity site, which is in line with the patterns of the whole bacteria along a salinity gradient with a shift in the dominance of Betaproteobacteria in ecosystems influenced by freshwater inputs into a predominance of Alphaproteobacteria in the high-salinity conditions (Bouvier and del Giorgio, 2002; Piao et al., 2012).

4.2 Distribution and environmental drivers of benthic anammox bacteria

It has been reported that functional gene marker hzo can describe anammox bacterial ecology more comprehensively than 16S rRNA because of its linking with anammox activity (Li et al., 2011a; Dang et al., 2013), however, hzo products are usually difficult to be obtained due to lower ratios of anammox bacteria compared with the whole bacteria. In the present study, no target hzo product was obtained, so only anammox bacterial 16S rDNA libraries were shown. Similar to nosZ-denitrifiers, the anammox bacterial community also exhibited significant regional patterns. At the central sites, the anammox bacterial community was dominated by Ca. Scalindua, while at the coastal sites, besides Ca. Scalindua, Ca. Brocadia and Ca. Kuenenia were also detected. In the BS sediments, Dang et al. (2013) have found Ca. Jettenia, but no Ca. Brocadia and Ca. Kuenenia. Other previous studies have also found the coexistence of non-Scalindua anammox bacteria along with Scalindua in coastal environments that received strong terrestrial influences (Amano et al., 2007; Dang et al., 2010; Li et al., 2010; Shehzad et al., 2016). Even so, the non-Scalindua clades generally accounted for a small fraction of anammox bacteria in marine environments (Dang et al., 2010, 2013; Li et al., 2010; Shehzad et al., 2016), and were considered to exist allochthonously without activities (Amano et al., 2007). In this study, most Brocadia sequences were similar to those in the suspended sediments of the Huanghe River, indicated that these anammox bacteria might be introduced from river runoffs. However, inconsistent with other results (Amano et al., 2007; Dang et al., 2010, 2013; Li et al., 2010; Shehzad et al., 2016), Brocadia fulgida accounted for more than 60% of total clones of the Huanghe River estuary site in this study, suggested that this phylotype was likely gradually adapted to the dynamic habitat, where river-sea interaction was intensive with lower salinity and richer nutrient and became a prominent group. Brocadia fulgida is generally hypothesized to be freshwater-adapted anammox bacteria and has extremely low tolerance to salinity (Gonzalez-Silva et al., 2017). Recently, nevertheless, Malovanyy et al. (2015) found that Brocadia fulgida was the major anammox phylotype and active in removing nitrogen when slowly adapt to the salinity of 15 in a wastewater bioreactor. Therefore, Brocadia fulgida would be highly possible to perform active nitrogen removing in the marine habitat not only survive and hibernate (Amano et al., 2007). This assumption could be supported by the recent investigation of Zhang et al. (2018), in which the high anammox activity was detected using ¹⁵N tracing technology at the same site.

Salinity was evidently a significant factor governing anammox bacterial distributions in this study. The results are not surprising given that previous studies have demonstrated that salinity influenced the geographical distribution of anammox bacteria in estuarine sediments (Dale et al., 2009; Hou et al., 2013). Recently, Sonthiphand et al. (2014) concluded that salinity drove the distribution of anammox bacteria at the global scale. The important role of salinity may be related to the different salinity tolerance of anammox bacteria. Nevertheless, the direct influence of salinity on anammox bacterial community is difficult to be confirmed in natural environments due to the covariation of salinity and other factors (e.g. NO₃^-, NO₂^- and NH₄⁺) (Sonthiphand et al., 2014). Scalindua adapts better to high-salinity habitats and can also attribute to its higher affinity of NO₂^- and NH₄⁺ that is commonly limited in highsalinity habitats (Sonthiphand et al., 2014).

4.3 Responses of nosZ-denitrifiers and anammox bacteria to coastal anthropogenic perturbations

In the eutrophic coastal ecosystems, in addition to natural factors, such as freshwater dilution, currents, tides, waves, upwelling, lateral transport, and water mixing, anthropogenic perturbations contributed importantly to the dynamics of marine denitrifying and anammox bacterial community diversity and composition (Dang et al., 2013; Babbin et al., 2016).

Highest primary productivity (represented as Chl-a) and highest contents of sediment heavy metal Pb, As, and Cd were found at the coastal site near the Huanghe River estuary, due to excessive anthropogenic nitrogen loading and wastewater discharge (Table S1). These factors were found significant in constraining the distributions of nosZ-denitrifiers and anammox bacteria (Table 4; Fig. 6). The primary productivity can supply carbon source for heterotrophic denitrifiers, and a few denitrifier phylotypes respond quickly to labile organic matter pulses and suit to highly productive conditions (Babbin et al., 2016). Anammox bacteria gain NH₄⁺ and NO₂^- from oxidizing organic materials, and some special species can directly utilize small organic compounds for alternative sources of NH₄⁺ (Van De Vossenberg et al., 2013). The biogeochemical cycling of heavy metals is coupled to nitrogen cycling in environments, and the significance of heavy metals to denitrifiers and anammox bacteria was also observed in other coastal ecosystems (Dang et al., 2010, 2013; Yang et al., 2015). Metal Pb, As, and Cd, no biological role, are potentially toxic to microorganisms. Interestingly, the non-Scalindua clusters were strongly-positively correlated with metal As, Cd, possibly because that these anammox bacteria were originally found in wastewater treatments with high tolerance to various heavy metals.

Additionally, a few big aquaculture farms located around the coastal site BF01, in which extreme high concentrations of NO₃^- and NO₂^- were detected during May owing to large discharges of aquaculture effluents. The high NO_x^- induced the shift of community structure of nosZ-denitrifier and anammox bacteria (Table 4), as reported in other coastal sediments (Dang et al., 2010). Especially, Scalindua wagneri took over Scalindua marina as the dominant anammox bacteria in that site and was strongly-positively correlated with the concentration of NO_x^-. Scalindua wagneri was first discovered in the wastewater treatment plant, which is typically associated with high nitrate and nitrite loads. This could explain why Scalindua wagneri adapted better to the high-NO_x^- habitat in this study.

5 CONCLUSION

In summary, the Bohai Sea and North Yellow Sea sediments harbored diverse nosZ-denitrifying and anammox bacteria assemblages with clear heterogeneity between coastal and central regions. Salinity was the most important environmental factor driving the distributions of nosZ-denitrifying and anammox bacteria community. In addition, anthropogenic perturbations (e.g. nitrogen overloading and consequent high primary productivity, and contaminant heavy metal discharges) contributed significantly to the communities of the two functional microbes and consequently changed the ecological functions of coastal systems.

6 DATA AVAILABILITY STATEMENT

The sequence data generated during the current study are available in the GenBank nucleic acid sequence database. The datasets analyzed during the current study are available from the corresponding author on reasonable request.

7 ACKNOWLEDGEMENT

We thank all the scientists and cruise members of the R/V Dong-Fang-Hong 2 for their assistance in sample collection during the expedition.