The bottom water and surface sediments have been gathered from 16 stations in the Yangtze Estuary and adjacent sea from 1 to 6 March 2022 (Figure 1). At each station, surface sediments (0–5 cm) was extracted from the box cores, mixed homogeneously, and placed into air-tight sterile plastics bags. The sediment samples were then subdivided into two parts for different analyses. One part immediately stored at 4°C for denitrification and anammox rates incubation. Another was stored on board (-20°C) and then transferred to laboratory storage (−80°C) for DNA extraction and molecular analysis. The bottom water was filtered through a filter (0.45 μm) and stored in polyethylene bottles (-20°C). The salinity was determined onsite by HQ40d Water Analyzer (Hach, Loveland, USA). The sediment water content was calculated using wet-dry weight method (Klute, 1986). The NH₄⁺, NO₃⁻ and PO₄³⁻ in bottom water were analyzed using QuAAtro nutrient autoanalyzer (Seal Analytical, Germany). 2 M KCl was used to extract the exchangeable NH₄⁺, NO₂⁻ and NO₃⁻ from the sediments and then determined by QuAAtro nutrient autoanalyzer (Seal Analytical, Germany). Sediment TN and TOC were detected by the elemental analyzer (Vario EL cube, Germany). A laser particle sizer analyzer (LS 13320, Beckman, USA) was used to measure sediment grain size (MФ). The sediment Fe(II) and Fe(III) contents were measured by phenanthroline spectrophotometric method (Lovley and Phillips, 1987). Under the dark condition of 4°C for 24 hours, Sediment Chla can be extracted with 80% (v/v) acetone. A UV-5200 visible spectrophotometer was used to measure the absorbance of Chla after centrifuging these extracts at 2200 g for 10 min.

The potential denitrification and anammox rates were measured by ¹⁵N isotope-tracing technique (Lin et al., 2017). Briefly, the laboratory slurry incubation experiments were made with the in situ benthic sterilized seawater and the sediments at a seawater/sediment mass ratio of 7:1, and the homogeneous mixture was fed with helium gas for 30 min before transferring into 12-mL vials (Labco Exetainers, UK). Subsequently, these 12-mL vials were in dark preincubated for 24 h to eliminate residual dissolved oxygen, nitrite, and nitrate. After the preincubation, three treatments were given to these vials as follows: (1) ¹⁵ NH₄⁺ (100 μmol N L⁻¹), (2) ¹⁵ NO₃⁻ (100 μmol N L⁻¹), and (3) ¹⁵ NH₄⁺ + ¹⁴ NO₃⁻ (100 μmol N L⁻¹). After 0 h and 8 h incubation, these vials were spiked with ZnCl₂ solution (200-μL 50%) to stop microbial activity. Produced ²⁹N₂ and ³⁰N₂ concentrations during incubation were determined by the membrane inlet mass spectrometry (MIMS, Hiden Analytical Ltd, UK), and denitrification and anammox rates were calculated according to previous methods (Thamdrup and Dalsgaard, 2002; Song et al., 2016).

Total DNA extraction was performed from 0.5 g of homogenized sediment samples by PowerSoil^® DNA Isolation Kit (MoBio, USA). The obtained DNA was used as the template for real-time qPCR analysis with the primer pairs cad3Af/R3cd for the nirS gene, and FlaCu/R3Cu for nirK gene, and Amx-808-F/Amx-1040-R for Anammox 16S rRNA gene (Supplementary Table 1). Reactions were carried out using ABI PRISM^®7500 Sequence Detection System (Applied Biosystems, Foster City, CA, USA). The 20-μL qPCR mixture reaction system contained FastStart Universal SYBR Green Master (ROX) (10 μL, Roche, Germany), ddH₂O (6.6 μL), each primer (0.6 μL, 10 μM), bovine serum albumin (0.2 μL, 20 mg·mL⁻¹), and template (2.0 μL). The qPCR reactions of nirS, nirK, and Anammox 16S rRNA gene were performed as the following: 50°C for 2 minutes, 95°C for 10 minutes, 45 cycles of 95°C for 30 seconds, and 58°C for 40 seconds, 72°C for 60 seconds (Throbäck et al., 2004; Hamersley et al., 2007). The significant linear relationships (R² > 0.996) were produced with a 10-fold dilution series (10²–10⁹ copies) of standard plasmids containing the key functional genes, and the qPCR efficiency was 84%, 110%, and 92%, respectively. The abundance of the key functional genes can be calculated by the standard curves. The melting curves contained a single peak, and non-specific amplification was not found by agarose gel electrophoresis, demonstrating amplification specificity.

11 sampling stations have been picked from 16 stations according to the class of muddy and sandy sediments (the value of median grain size (> 64 μm) was classified to sandy sediment). The PCR amplicons of key functional genes were from the same primer pairs with qPCR (cd3aF/R3cd and FlaCu/R3Cu for denitrifiers, and Amx-808-F/Amx-1040-R for anammox bacteria). The PCR conditions were as follows: initial denaturation was performed at 98°C for 2 min, followed by control at 98°C for 15 s, 55°C for 30 s and 72°C for 30 s (25-30 cycles), and finally extended at 72°C for 5 min. The Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) and Quant-iT PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA) were used for the purification and quantification of PCR amplicons, respectively. The Illumina novaSeq PE250 platform (Shanghai, China) was used for the paired-ended sequenced of high-throughput sequencing of nirS/nirK and AMX 16S rRNA genes. The QIIME2 (2019.4) was used to incorporate, quality-filter, and dereplicate the sequences. High quality sequences with 97% similarity level were clustered into OTUs based on Vsearch (v2.13.4_linux_x86_64), and then representative OTUs were output. The raw sequences obtained in this study for nirS, nirK and AMX 16S rRNA was kept in NCBI under the accession numbers of PRJNA891470, PRJNA891480, and PRJNA891491.

The map of the sampling stations and the distribution of environmental parameters were drawn using Ocean Data View (ODV 5.2.1). The differences of nirS, nirK and AMX 16S rRNA gene abundances in the same station were studied using one-way analysis of variance (ANOVA). Correlation between the three microbial communities and environmental factors were identified with CCA/RDA analyses with Monte Carlo tests by Canoco 5 software. The Pearson and Spearman correlation analysis of denitrifiers and anammox bacteria with environmental parameters was created by IBM SPSS statistic 26 and R (version 4.2.1, R package corrplot) (Wei and Simko. 2021). MEGA11, TBtools and Origin 2022b software were used for graphics. Finally, the phylogenetic tree was visualized and beautified by Evolview online tool (http://www.evolgenius.info/evolview-v2/).

The physicochemical parameters of bottom water and surface sediments were given in Figure 2. The salinity ranged from 27.3 to 37.1 (Figure 2A). Bottom water NH₄⁺ (0.1–3.15 μM), NO₃⁻ (3.55–23.93 μM), PO₄³⁻ (0.21–0.71 μM) values varied across different sampling sites (Figures 2B–D). The average concentrations of NH₄⁺ and NO₃⁻ in the sediments varied from 0.32–1.20 μg N g^–1, 0.04–0.86 μg N g^–1, respectively (Figure 2E, F). The sediment Fe(II)/Fe(III) ratios were in a range of 0.29–6.19 (Figure 2G). The contents of sediment TN and TOC ranged from 0.23 to 1.19 mg N g⁻¹ with an average of 0.51 ± 0.25 mg N g⁻¹ and from 1.05 to 6.97 mg C g⁻¹ with an average of 2.62 ± 1.65 mg C g⁻¹, respectively (Figures 2H, I). The sediment median grain size showed a greater change. The sampling stations of A1-1, A1-8, A3-2, and A4-3 were muddy sediments (8.83–44.27 μm), while other stations were sandy sediments (93.32–333.41 μm) (Figure 2L). Meanwhile, the C/N ratios (TOC/TN) were higher in muddy sediments (A1-1, A1-8, A3-2, and A4-3) than in sandy sediments (Figure 2J). The average content of Chla in sediment was in a range of 0.12–1.75 μg g⁻¹, with obvious spatial pattern (Figure 2K).

The nirS/nirK gene and AMX 16S rRNA gene were utilized to quantify denitrifiers and anammox bacteria abundances in surface sediments, respectively. The abundances of nirS gene were ranging between 9.87 × 10⁵ to 2.71 × 10⁷ copies g⁻¹ with the mean value of (1.28 ± 0.93) × 10⁷ copies g⁻¹ (Figure 3A, B). The nirS gene abundance had negative correlation with NO₃⁻ concentration in the bottom water (r = -0.576, p< 0.05) (Figure 4). The nirK gene abundances ranged from 4.64 × 10⁵ to 1.26 × 10⁷ copies g⁻¹ with the mean value of (4.52 ± 3.34) × 10⁶ copies g⁻¹ (Figure 3A, B). Compared to nirS/nirK gene, the abundances of AMX 16S rRNA gene were obtained between 1.03× 10⁵ to 2.97 × 10⁶ copies g⁻¹with the mean value of (1.40 ± 0.84) × 10⁶ copies g⁻¹ (Figure 3A, B). In all surface sediment samples, nirS abundances were the highest, followed by nirK abundances and then AMX 16S rRNA gene.

The potential rates of denitrification and anammox were analyzed in all sampling site of the surface sediments (n = 16), which varied considerably, ranging from 0–0.56 nmol N g⁻¹ h⁻¹ and 0–0.091 nmol N g⁻¹ h⁻¹, respectively (Figure 3C, D). The highest denitrification rate was detected at A1-1, whereas the lowest denitrification rate was observed at A3-8. The potential denitrification rates had positive association with sediment water content, NH₄⁺, TN, TOC, C/N, and Chla, while negatively regulated by sediment NO₃⁻ and sediment grain size (Figure 4, Supplementary Table 5). Compared to denitrification rates, the highest anammox rate was detected at A2-3, whereas the lowest anammox rate was observed at A3-4. Anammox was positively correlated with NH₄⁺ and Chla, and negatively correlated with NO₃⁻ (Figure 4, Supplementary Table 5). Moreover, Pearson’s correlation analyses showed that anammox was significantly related to denitrification (r = 0.710, p< 0.01) (Figure 4, Figure 5A; Supplementary Table 5). The contribution of denitrification to N removal (DEN%) was higher than anammox in most sampling sites except A3-4, A3-6, A3-8, and A4-1 (Figure 5B), with the contribution of 73.5% on average. In addition, DEN% was positively related to water content, NH₄⁺, Fe(II)/Fe(III), TN, TOC, C/N, and Chla, and negatively regulated by NO₂⁻.

Totally, 1,219,102 filtered, high-quality nirS sequences were acquired from 11 surface sediments samples, ranging from 73,783–172,218 for each sample (Supplementary Table 2). The effective sequences were clustered into 7,536 OTUs based on the 97% similarity level. Additionally, Good’s coverage (99.73–99.91%) indicated that most of the nirS-encoding denitrifiers had been obtained. The highest OTU numbers (1,186) and Chao1 richness estimators (1267.16) were observed in the samples A4-7. Higher Shannon diversity index and Simpson indices were found in samples A2-3, A4-5 and A4-7 (Supplementary Table 2). Moreover, the diversity of nirS gene decreased with the increasing the NH₄⁺ concentration in the bottom water (Shannon: r = -0.732, p< 0.05; Simpson: r = -0.897, p< 0.01) (Supplementary Table 6). Compared to the nirS gene, 763,468 high-quality reads of nirK gene sequences were acquired and 3,151 OTUs were identified (Supplementary Table 3). The diversity and richness of nirK-encoding were lower than nirS-encoding denitrifiers in all samples. In additionally, the richness of nirK gene declined with the increasing NO₃⁻ concentrations in the bottom water (Chao1: r = -0.606, p< 0.05) and increased with the sediment total nitrogen concentration (Chao1: r = 0.603, p< 0.05) (Supplementary Table 6). The rarefaction curves showed that the majority of the nirS- and nirK-encoding denitrifiers was captured in the current study (Supplementary Figures 1A, B).

A total of 50 dominant OTUs of the nirS gene were phylogenetically classified into 10 different clusters (Figures 6A, B; Supplementary Figures 2, 3). Within these 10 clusters, only C9 belonged to α-proteobacteria and C5, C6, and C7 affiliated to β-Proteobacteria, whereas the other 6 clusters were belonged to γ-Proteobacteria. C3 only contained OTU01, which was existed exclusively in the A3-2 station, closely affiliating to Pseudomonas (Figures 6A, B, Supplementary Figure 3). C8 (belonging to Marinobacter), C9 (belonging to Paracoccus) and C10 (belonging to Chromohalobacter) were observed almost in muddy sediments (Figures 6A, B). According to a redundancy analysis (RDA) with the dominant OTUs, the first two axes (RDA1 and RDA2) explained 79.97% of the cumulative variance in the nirS-encoding denitrifiers communities (Figure 6C). The bottom water NH₄⁺ (pseudo-F = 3.1 and p< 0.01), NO₃⁻ (pseudo-F = 2.7 and p< 0.01), and C/N (pseudo-F = 2.6 and p< 0.05) were the major factors influencing the nirS gene community. At genus levels, Pseudomonas was dominant genera only in A3-2, while Marinobacter was dominant genera and observed in all sampling stations (Supplementary Figure 4B). Spearman’s correlation analysis showed that Pseudomonas was positively related with in the bottom water (p< 0.01) (Supplementary Table 7), and Marinobacter was positively correlated with sediment NH₄⁺ (p< 0.05).

Dominant OTUs (top 50 OTUs) of nirK gene were classified into 9 clusters, all belonging to the Proteobacteria (Figures 7A, B, Supplementary Figures 5, 6). Among the 9 clusters, C2, C5, and C7 belonged to the β-Proteobacteria, whereas the other 6 clusters were intimately related to the α-Proteobacteria (Figure 7B, Supplementary Figure 6). Moreover, C7 (with the lowest abundances of top 50 OTUs) only contained OTU32 which was detected in the A3-2 and A3-4 sample, showing some similarities with the sediment denitrifying microorganisms in the Yangtze lakes (Jiang et al., 2017). There were no obvious differences in the relative abundance of OTUs between muddy sediments and sandy sediments (Figure 7A). From the CCA (Figure 7C), the first two axes explained 36.10% of the cumulative variance, and NH₄⁺ in the bottom water (pseudo-F = 2.0 and p< 0.01), PO₄³⁻ (pseudo-F = 2.1 and p< 0.05), sediment NO₃⁻ (pseudo-F = 2.2 and p< 0.05), and MФ (pseudo-F = 2.0 and p< 0.05) had significant influences on the nirK-encoding denitrifiers community. At genus level, Mesorhizobium was the dominant genera in A2-5 and A3-1 sample with no other genera observed, however, Bradyrhizobium was detected in most sampling stations except A2-5 and A3-1 (Supplementary Figure 7B). Spearman’s correlation analysis showed that Mesorhizobium was positively related with sediment Fe(Ⅱ)/Fe(III) (p< 0.05) (Supplementary Table 7).

For the anammox 16S rRNA gene (AMX 16S), a total of 479,390 high-quality reads were gained, which ranged from 40,134–47,445 for each sample (Supplementary Table 4). Based on the 0.03 genetic distance (97% as cutoff), 23,950 OTUs (2177 ± 440 OTUs per samples, n=11) were detected. The highest OTU numbers, Chao1 and Shannon index were found in station A4-7, showing the higher diversity and richness of the anammox 16S rRNA gene in surface sediment of station A4-7 than the other stations. In addition, Pielou’s evenness (from 0.76 to 0.84) and Good’s coverage values (from 98% to 99%) can be observed from each sample, demonstrating the dependability of the sequencing accuracy. The community diversity of the AMX 16S rRNA gene decreased with the increasing concentration in the bottom water (Chao1: r = -0.681, p< 0.05; Observed species: r = -0.710, p< 0.05), decreasing sediment NH₄⁺ (Chao1: r = 0.677, p< 0.05; Observed species: r = 0.678, p< 0.05), and decreasing Chla concentration (Chao1: r = 0.730, p< 0.05; Observed species: r =0.722, p< 0.05) (Supplementary Table 6). The rarefaction curves indicated that the majority of the anammox bacteria was captured in the current study (Supplementary Figure 1C).

The dominant OTUs (top 10 OTUs) of AMX 16S rRNA gene and their related sequences from the GenBank were phylogenetically analyzed to 4 clusters, including Ca. Anammoximicrobium moscowii and the other uncultured clusters (Figures 8A, B; Supplementary Figures 8, 9). The OTUs (OTU173, OTU159, and OTU65) were closely related to Ca. Anammoximicrobium moscowii, and the OTU abundances were higher in muddy sediment than those in sandy sediment. The other clusters were related to uncultured anammox bacteria, showing a very wide range and matched with marine habitats and activated sludge (Figure 8B). According to a redundancy analysis (RDA), the link between AMX 16S communities and environmental factors may be described by the first two axes (65.75% of the total variance) (Figure 8C). The sediment MФ (pseudo-F = 6.1 and p< 0.01), Chla (pseudo-F = 2.3 and p< 0.05), and bottom water NO₃⁻ (pseudo-F = 3.4 and p< 0.05), had huge impact on the anammox community. At genus level, Candidatus Anammoximicrobium was the dominant genera in all sampling stations except A3-2, while the dominant genera in A3-2 was Rhodoplanes (Supplementary Figure 10B). Spearman’s correlation analysis revealed that Candidatus Anammoximicrobium was negatively related with sediment water content (p< 0.05), TOC (p< 0.05), and positively correlated with sediment grain size (Supplementary Table 7).

In coastal aquatic ecosystems, denitrification in sediments can reduce bioavailable nitrogen from anthropogenic and natural inputs, thus, sediments are recognized as hotspot regions of N loss (Seitzinger et al., 2006). Lin et al. (2017) reported that the contribution of denitrification to N loss (2.2 × 10⁶ t N yr ^–1) was higher than anammox (4.6 × 10⁵ t N yr^–1) in the surface sediments from the East China Sea. The biogeochemical cycle of nitrogen is driven by microorganisms, and the microbial communities are constantly changing to adapt to the complex coastal environments (Dang and Lovell, 2016; Hutchins and Fu, 2017). Thus, a further understanding of the microbial processes of N loss is still highly significant. Here, potential rates of denitrification and anammox, and their related gene abundances and community composition, along with influencing environmental factors were explored in the surface sediments from the Yangtze Estuary and the adjacent sea.

In the present study, the average abundance of nirS gene (1.28 × 10⁷ copies g⁻¹) in the surface sediment was greater than that of nirK gene (4.52× 10⁶ copies g⁻¹), which were similar to the observations in the San Francisco Bay (1.50 × 10⁷ copies g⁻¹ for nirS and 3.20× 10⁶ copies g⁻¹ for nirK) (Mosier and Francis, 2010). Some previous studies showed that the mRNA abundance of nirS was also higher compared with nirK, suggesting a main role in N removal processes of nirS-encoding denitrifiers (Graham et al., 2010; Liu et al., 2021). The differential abundance between nirS and nirK may be due to the competition for the substrates required by different enzyme synthesis (Huang et al., 2011). Moreover, compared to nirK-encoding denitrifiers, the variety of nirS-encoding denitrifiers was also greater, which is consistent with previous studies (Shi et al., 2019b; Ming et al., 2021). In general, nirS- and nirK-encoding denitrifiers had different respond to environmental parameters, and nirS denitrifiers had better adaptation to the coastal environments (Jones and Hallin, 2010; Likhitrattanapisal et al., 2021).

In consistent to our study (Figure 3), the abundance of nirS and nirK have the same change trends across the sampling sites (Zhang et al., 2019), which was affected by both the microbial interactions between nirS and nirK- denitrifiers, as well as the different preferences in structure components of corresponding nitrite reductases (Cutruzzolà et al., 2001). Nitrite reduction is the key step in denitrification, and both the nirS and nirK- encoded nitrite reductase are functionally identical while being evolutionarily distinct (Shrewsbury et al., 2016). The previous studies reported that both enzymes are often absent in the same bacterium, and the estimated number of denitrifiers in the soil sites was determined using the summation of nirS and nirK copy counts (Philippot et al., 2011; Zhang et al., 2019). However, Pseudomonas was found in nirS- and nirK-encoding denitrifiers at genus level in our study (Supplementary Figures 4B, 7B). Furthermore, there other different genera were detected in nirS- and nirK-encoding denitrifiers, which have capabilities of denitrification (Kim et al., 2008; Borrero‐De Acuña et al., 2017; Gao et al., 2021). Thus, marker genes of nirS and nirK should be used to detect denitrifiers together. In addition, there are another four different genes involved in the denitrification pathway, nitric oxide reduction (norB and qnorB), and nitrous oxide reduction (nosZI and nosZII). A single gene does not capture the entire denitrifying microbial community, due to several maker genes involved in denitrification pathway (Shrewsbury et al., 2016). Furthermore, Zhang et al. (2019) found that the contribution of denitrification to N removal was positively related to the relative abundance of nosZ/nirS, implying that marker genes from different denitrification steps should be considered simultaneously.

According to phylogenetic analysis, the dominant OTUs of nirS showed some similarities to multiple estuarine and marine sediments (Dang et al., 2008; Hong et al., 2014; Zheng et al., 2015; Gao et al., 2016). In contrast to nirS, the nirK-encoding denitrifiers were found to originate from several habitats, such as marine sediments, fresh sediments, agricultural soil, and wastewater treatment systems (Figure 7B, Supplementary Figure 6). Numerous studies have revealed that nirS and nirK- denitrifiers occupy different niches or have distinct habitat preferences in the environment (Smith and Ogram, 2008; Hallin et al., 2009; Knapp et al., 2009). The communities of nirS and nirK gene were observed to clearly separate in saline environments, and nirS-encoding denitrifiers predominated in marine environment (Jones and Hallin, 2010), in consistent with our study. Furthermore, nirK denitrifiers might exhibit more habitat selection, showing positive correlations with organic matter content (Smith and Ogram, 2008). In the present study, RDA/CCA demonstrated that C/N ratio and MФ was the key factor affecting the community of nirS and nirK, respectively (Figure 6C, Figure 7C). High level of C/N ratio may promote the growth of denitrification bacteria, as previous studies showed that the denitrifiers OTUs and Shannon index decreased with decreasing organic matter and N contents (Chen et al., 2018; Wu et al., 2021). Meanwhile, the dominant OTUs of nirS- and nirK-encoding denitrifiers were affiliated with Proteobacteria (Figures 6B, 7B, Supplementary Figures 3, 6), which was consistent with previous reports (Kim et al., 2016; Xie et al., 2020). Furthermore, γ- Proteobacteria and α-Proteobacteria was the dominant class of nirS- and nirK-encoding denitrifiers, respectively. Some studies showed that denitrification can be performed by phylogenetically different microbes, and the most of denitrifiers early reported are belonging to subphyla of the Proteobacteria (Peralta et al., 2010; Wu et al., 2021). This phylum of Proteobacteria are widely distributed in aquatic environments because of high phylogenetic and phenotypic diversities (Liu et al., 2015; He and Zhang, 2016; Shi et al., 2019a). With the Proteobacteria community, Pseudomonas and Mesorhizobium were the dominant genera of nirS- and nirK-encoding denitrifiers, respectively. Pseudomonas stutzeri strain XL-2 and strain XL-2 were identified as the main denitrifiers which can remove 97.9% of nitrate (Dang and Lovell, 2016). In nirS bacterial community, Pseudomonas was mainly affected by NH₄⁺ in the bottom water; while in nirK bacterial community, Mesorhizobium was significantly influenced by with sediment Fe(II)/Fe(III). In addition, it has been reported that the higher abundance of Mesorhizobium was observed in sandy soil (Li et al., 2019), in consistent to our study that Mesorhizobium was the dominant genera in A2-5 and A3-1 with high median grain size.

In the present study, the abundance of AMX 16S rRNA gene was ranging from 1.03× 10⁵ to 2.97 × 10⁶ copies g⁻¹ (Figures 3A, B), which was lower than that of nirS and nirK gene in all sampling sites. The gene abundances of AMX 16S rRNA gene were like to the values in the Pearl River Estuary (Wu et al., 2020), but 1-2 orders magnitude lower than those in paddy soil (Zhu et al., 2011) and the coastal wetland sediments of China (Hou et al., 2015). Furthermore, anammox bacterial abundance was no obvious spatial heterogeneity at the sampling sites (Figure 3A).

The high-throughput sequencing of the 16S rRNA gene of anammox bacteria could provide adequate data, as compared to data from cloning libraries, to identify community changes and even recognize unusual bacteria with low relative abundance (Zhao et al., 2013; Shen et al., 2016). Based on the OTUs and the Shannon index, the AMX 16S genes in the surface sediment represented a very varied community of anammox bacteria (Supplementary Table 4). A wide variety of anammox bacteria was found in the surface sediment of Yangtze Estuary and the adjacent sea, in similarity to other estuarine sediments (Dale et al., 2009; Hong et al., 2014; Lisa et al., 2014). Salinity is a significant environmental factor that affects the diversity of anammox bacteria, according to earlier research (Schmid et al., 2007; Hamersley et al., 2009; Hu et al., 2012). However, in this present study, the Chao1 had significantly correlated with NH₄⁺ concentration in bottom water and sediment (Supplementary Table 6). These results implied that the substrates contents might be the main limiting factors for anammox bacteria.

In this study, the dominant OTUs were closely related to Ca. Anammoximicrobium moscowii and Candidatus Anammoximicrobium was the dominant genera at genus level (Supplementary Figures 9, 10B). In contrast to previous studies (Schmid et al., 2007; Fernandes et al., 2016; Fu et al., 2019; Wu et al., 2021), Ca. Scalindua are the dominant anammox bacteria in marine environments. Some studies reported that Ca. Scalindua has a relatively high tolerance for salinity and it can adapt to high hypersaline environment (Oshiki et al., 2017; Speth et al., 2017). Recently, the other genera have been found and Ca. Scalindua was not the only genus in the marine environments (Wu et al., 2019). Ca. Anammoximicrobium moscowii also can carry out anaerobic ammonium oxidation, which was firstly isolated in the Moscow region (Khramenkov et al., 2013), and the key role of this genus in N loss process cannot be ignored. Some studies reported that the different anammox bacteria lineages can become the dominant populations based on the slight change in the environment (Wang et al., 2021). In addition, the community composition study of anammox bacteria can be directly impacted by primer sets and sequencing depth. In the present study, considering the primer efficiency, Amx808F and Amx1040R was the appropriate choice for capturing anammox bacteria with good coverage, however, the specificity was not strong for Ca. Scalindua detection. From the RDA analysis, the MФ and bottom water NO₃⁻ significantly affected anammox bacteria community. In addition, according to Spearman’s correlation analysis, Candidatus Anammoximicrobium was negatively related by TOC.

For microbial mediated N transformation, the surface sediment are usually the most vibrant with greater potential N loss rates, particularly the denitrification rates; however, the potential rates of denitrification and anammox were lower than previous studies (Hou et al., 2015; Lin et al., 2017; Wu et al., 2021), denitrification and anammox can exhibit strongly spatiotemporal variabilities, and the rates might be affected by biotic and abiotic factors in sampling sites. The bulk of the samples belong to sandy sediments, with the low organic matter content and low N removal rates. In our study, denitrification was the main contributor of N loss, dominating the N loss process in surface sediment in most sampling sites (Figure 5B). These results were in accordance with those in the previous studies in estuarine and coastal habitats (Fernandes et al., 2012; Cao et al., 2017; Lin et al., 2017; Huang et al., 2021). Moreover, Huang et al. (2021) found that denitrification rates were substantially correlated with organics and median grain size (MФ), implying that both MФ and organic matter were the crucial factors regulating denitrification activities. High surface-to-volume ratios in fine-grained sediments typically result in greater organic matter absorption (Keil et al., 1994). Meanwhile, organic matter (as an electron donor) can directly or indirectly boost N removal processes by providing energy via the oxidation of organics (Devol, 2015). Additionally, we verified that denitrification rates were substantially correlated with MФ and the C/N ratio (Figures 5C, D). Higher specific surface area in sediment with little grains was more suited to microbial colonization, and the required energy for the growth of heterotrophic denitrifiers were provided by organic matter Moreover, the contribution of anammox was higher than denitrification in few sampling stations with lower C/N ratio and higher median grain size, further implying that sediment texture can control the N loss rates via affecting the residence time of organic matter (Huang et al., 2021).

Anammox rates and denitrification rates have a positive correlation in this study (Figure 5A), suggesting that denitrification was probably the main source of NO₂⁻ for anammox. Coherent to many previous researches in the aquatic environment (Brin et al., 2017; Zhao et al., 2020a), anammox metabolism requires the substrates of NO₂⁻, which are main supplied by denitrification rather than nitrification in these anaerobic environments. This coupling between anammox and denitrification may be mediated by an imbalance in NO₃⁻ and organic carbon needed for denitrification. Furthermore, in the present study, NO₂⁻ was negatively related to DEN%, but positively correlated to ANA%, implying the coupling of denitrification and anammox. However, there also studies showing that anammox coupled with nitrification rather than denitrification in estuarine sediments in winter (Hou et al., 2015; Lin et al., 2017). The uncertainties in the relationships between anammox and denitrification might be attributed to organic acids, and anammox bacteria can convert NO₃⁻ to NO₂⁻ in the presence of organic acids (Brin et al., 2014). Thus, coupled anammox and nitrification or denitrification play a vital role in N removal processes from the surface of marine sediments.

Previous studies have revealed that NH₄⁺ and NO₃⁻ play the crucial and central role in N removal (Lin et al., 2017; Huang et al., 2021), which were as substrates in denitrification and anammox processes. In consistent to our study, NH₄⁺ and NO₃⁻ were the key factors to influence denitrification and anammox rates. In the present study, Chla was the key environmental regulator for denitrification rates (p< 0.01) (Figure 4, Figure 5E, Supplementary Table 5). The sediment Chla is an indicator for the availability of labile organic matter which is obtained from phytoplankton deposited from upper waters, along as phytobenthos and benthic algal debris in the shallow marginal system (Franco et al., 2007; Kitidis et al., 2017). The sediment organic matter correlates positively with Chla in sediment and water column (Li et al., 2013), in similarity to our study (Figure 4, Supplementary Table 5). The growth of phytoplankton may be aided by light intensity and water temperature in March (Shi et al., 2019b), indicating that many phytoplankton cells might be destroyed to release labile DOC (Zhang et al., 2018), which could enhance denitrification by sinking down to anaerobic sediment layers. Meanwhile, anammox rates also positively correlated with Chla (p< 0.01) (Figure 4, Figure 5F, Supplementary Table 5). The main reason for that was the substrates (NH₄⁺ and NO₂⁻) of anammox originated from the mineralization of organic matters (Zhang et al., 2018). Thus, both the quality and amount of organic matter in sediment are essential to manage estuarine N-cycling because fresh organics has a higher lability and simpler to be utilized by NO_x⁻ removal processes (Westrich and Berner, 1984; Henrichs and Doyle, 1986).

There were no significant correlations between N removal rates and related gene abundances in this study. The gene abundance obtained from the DNA levels has limitations. On the one hand, since DNA extracts include a wide variety of cells (e.g., active cells, inactive cells and dead cells) and its degraded or lysed extracellular DNA, gene abundance does not always indicate activity (Gaidos et al., 2011). On the other hand, the number of microbials just indicates potential microbial metabolism and has no direct impact on the activities of N transformation (Dai et al., 2008). It has been reported that the related functional gene transcript abundance is reflective of N transformation process potential, further revealing N transformation potential had a close relationship with active microbials (Li et al., 2022). Thus, to gain an insight into changes in denitrification and anammox potential, the functional transcript abundance combining with N removal rates are needed.