Key genes for microbial nitrogen metabolism were searched for in the metagenomes of two types of sedimental samples (i.e., with submerged plants, SP; with no plants, NP), and differences between SP and NP were clear regarding gene inventories. In total, more than 70 functional genes involved in nitrification, denitrification, assimilatory nitrate reduction, dissimilatory nitrate reduction, nitrogen fixation, anammox, organic degradation and synthesis were annotated in this study (Supplementary file 1). More than two-thirds of the functional genes involved in the pathways of nitrogen fixation, assimilatory nitrate reduction and denitrification were more abundant in the SP sediments than in the NP sediments (Fig. 1a). Specifically, genes encoding the glutamate synthase (NADPH/NADH) small chain (gltD) accounted for more than 20% of the total abundance of functional genes involved in the nitrogen cycle and had significantly higher abundance in SP than in NP (corrected p < 0.05). Similarly, the abundance of genes encoding enzymes of the nitrate reduction process, including assimilatory nitrate reductase catalytic subunit (nasA), ferredoxin-nitrite reductase (nirA), and nitrite reductase (NADH) (nirB/D), was significantly higher in SP than NP. In contrast, the abundance of genes encoding nitric oxide reductase (norB) and hydroxylamine dehydrogenase (hao) had significantly higher abundances (corrected p < 0.05) in NP than SP. For the nitrogen fixation pathway, genes encoding the nitrogenase molybdenum-iron protein (nifD/K), nitrogenase iron protein (nifH), and nitrogenase-stabilizing/protective protein (nifW) were more abundant in SP than NP. In addition, the abundance of genes encoding urease subunit alpha/gamma (ureC/A), the key enzyme in the conversion of organic nitrogen to ammonium nitrogen, was higher in the SP than in the NP sediments (Fig. S1 in Supplementary file 2). The correlation variations based on redundancy analysis (RDA) indicated that total phosphorus (TP), total nitrogen (TN) and C: N ratio (C/N) were the major factors shaping microbial communities (Fig. S2a in Supplementary file 2), and explain 59.66% of the community variances by the first two axes. The relationships between environmental factors and microbial functional genes involved in nitrogen metabolic pathways were further investigated by partial Mantel tests. We found the genes involved in assimilatory nitrate reduction (r = 0.512), denitrification (r = 0.425) and dissimilatory nitrate reduction (r = 0.372) were significantly (p < 0.05) correlated with C/N (Table 1).
Figure 1. Abundance comparison of key genes involved in nitrogen cycling (a), sulfur cycling (b), and methanogenesis (c) between sediments with (SP) and without (NP) submerged plants. The p values are based on Welch's exact test and corrected by Benjamini-Hochberg false discovery rate (FDR) method
Gene category pH Conductivity (μs/cm) TP (g/kg) TN (g/kg) TC (g/kg) C/N Nitrogen-cycling pathway Anammox − 0.124 (0.792) − 0.311 (0.954) 0.104 (0.263) − 0.035 (0.578) − 0.320 (0.956) 0.110 (0.313) Assimilatory nitrate reduction 0.254 (0.055) 0.347 (0.017) − 0.220 (0.919) − 0.438 (0.990) 0.050 (0.412) 0.512 (0.006) Denitrification 0.184 (0.094) 0.380 (0.014) 0.009 (0.428) − 0.027 (0.560) 0.299 (0.019) 0.425 (0.008) Dissimilatory nitrate reduction 0.421 (0.020) 0.446 (0.007) − 0.198 (0.932) − 0.086 (0.723) 0.322 (0.016) 0.372 (0.019) Nitrification 0.040 (0.384) 0.206 (0.109) 0.266 (0.100) − 0.076 (0.642) 0.284 (0.077) 0.085 (0.274) Nitrogen fixation 0.453 (0.007) 0.476 (0.003) − 0.270 (0.968) 0.207 (0.151) 0.321 (0.024) 0.049 (0.304) Organic degradation and synthesis − 0.081 (0.681) 0.093 (0.294) 0.127 (0.248) − 0.256 (0.915) − 0.004 (0.505) 0.331 (0.049) Sulfur-cycling pathway Assimilatory sulfate reduction 0.122 (0.222) 0.331 (0.028) − 0.208 (0.922) − 0.361 (0.987) 0.007 (0.511) 0.603 (0.005) Dissimilatory sulfate reduction and oxidation 0.180 (0.120) 0.257 (0.069) − 0.191 (0.838) − 0.273 (0.961) − 0.050 (0.603) 0.346 (0.013) Link between inorganic and organic sulfur transformation − 0.142 (0.812) 0.066 (0.389) − 0.056 (0.536) − 0.285 (0.906) − 0.104 (0.703) 0.426 (0.050) Organic sulfur transformation 0.026 (0.475) 0.147 (0.238) − 0.212 (0.878) − 0.301 (0.946) − 0.091 (0.703) 0.463 (0.021) SOX systems − 0.090 (0.717) 0.059 (0.413) 0.038 (0.364) − 0.349 (0.968) − 0.012 (0.535) 0.456 (0.016) Sulfur disproportionation − 0.016 (0.559) 0.069 (0.342) 0.088 (0.268) − 0.074 (0.676) 0.064 (0.331) − 0.061 (0.594) Sulfur oxidation 0.132 (0.167) 0.286 (0.024) − 0.110 (0.751) − − 0.060 (0.650) 0.097 (0.306) 0.379 (0.026) Sulfur reduction 0.077 (0.354) 0.137 (0.227) 0.122 (0.209) -0.226 (0.849) 0.128 (0.292) 0.222 (0.139) Methanogenesis Core methanogenesis − 0.295 (0.977) − 0.084 (0.692) 0.164 (0.214) − 0.221 (0.871) − 0.122 (0.71) 0.313 (0.101) Values are correlation coefficients with p values in brackets (p < 0.05 in bold)
TP: Total phosphorus; TC: Total carbon; TN: Total nitrogen; C/N: Total carbon/Total nitrogen
Table 1. Summary statistics for partial Mantel tests of correlation between functional genes and environmental factors
More than 200 functional genes involved in assimilatory sulfate reduction, sulfate reduction, dissimilatory sulfate reduction and oxidation, sulfur reduction, sulfur oxidization (SOX) systems, sulfur oxidation, sulfur disproportionation, organic sulfur transformation, and the link between inorganic and organic sulfur transformation, were identified in both NP and SP samples (Supplementary file 3). Among them, 26 genes encoding the key enzyme of sulfur cycling showed significantly higher (corrected p < 0.05) abundances in the SP than in the NP sediment samples (Fig. 1b). For the assimilatory sulfate reduction pathway, most relevant genes involved in encoding for sulfate adenylyltransferase subunit (cysD/N), phosphoadenosine phosphosulfate reductase (cysH), sulfite reductase (NADPH) flavoprotein alpha-component (cysJ) and sulfite reductase (ferredoxin) (sir) were significantly more abundant in SP than NP sediments. Also, SP had significantly higher abundances of genes encoding thiosulfate sulfurtransferase (glpE), sulfite dehydrogenase (quinone) subunit (soeA), sulfide-quinone oxidoreductase (sqr), and the key enzyme of sulfide oxidation (Fig. 1b). For the SOX systems, similar gene abundances were found in NP and SP, except for soxZ. In contrast, the NP microbial communities showed significantly higher (corrected p < 0.05) abundances of genes encoding 3-(methylthio)propanoyl-CoA dehydrogenase (dmdC; catalyzing organic sulfur oxidation), anaerobic dimethyl sulfoxide reductase subunit B (dmsB) and sulfide dehydrogenase subunit alpha (sudA). According to the RDA plots, the environmental factors TN, TP and C/N had significant influences on the relative abundances of functional genes involved in sulfur cycling, and explained 58.07% of the total variation by the first two axes (Fig. S2b in Supplementary file 2). Furthermore, C/N showed significant positive correlations with the functional genes involved in assimilatory sulfate reduction (r = 0.603, p < 0.05), dissimilatory sulfate reduction and oxidation (r = 0.346, p < 0.05), organic sulfur transformation (r = 0.463, p < 0.05), SOX systems (r = 0.456, p < 0.05), and the sulfur oxidation pathway (r = 0.379, p < 0.05) (Table 1).
Genes involved in methanogenesis are summarized in Supplementary file 4. The detected methanogenic processes included methanol, methylamine, acetate, carbon dioxide, core carbon and core methanogenesis. Gene encoding for the methyl coenzyme M reductase system, component A2 (mcrA2/K00400) was more abundant in the SP than NP sediments. Similarly, abundances of genes encoding the trimethylamine-corrinoid protein Co-methyltransferase (mttB) and 5, 10-methylenetetra-hydro-methanopterin reductase (mer) were significantly higher (corrected p < 0.05) in SP sediments than in NP sediments. In contrast, the abundance of gene associated with heterodisulfide reductase subunit A2 (hdrA2) was significantly higher in the NP than the SP sediments (corrected p < 0.05) (Fig. 1c).
To reveal the diversity of metabolic potential for methanogenesis, we chose relatively high-quality genome bins (completeness > 80%) and performed bin annotation to identify the metabolic pathways of hydrogenotrophic, acetoclastic, and methylotrophic methanogenesis in urban wetland sediments (Fig. 2). Specifically, SP_bin28 had genes encoding conserved core enzymes of hydrogenotrophic methanogenesis, including Fwd, Ftr, Mch, Mtd, Mer, Mtr and Mcr. NP_bin6 and SP_bin28 contained genes encoding Acs, Cdh, Mtr, and Mcr enzymes for the utilization of acetate, which showed potential for acetoclastic methanogenesis. SP_bin4 contained genes encoding methyl-compound methyltransferase, such as Mts, Mta, Mtb, Mtm and Mtt. These genes showed a potential for methane production through the methylotrophic pathway. However, no complete methanogenesis pathways were identified in NP_bin8 and SP_bin24 due to the incompleteness of their genomes.
The 56 retrieved MAGs (25 from NP and 31 from SP) could be assigned to more than Ten phyla, including one archaeal phylum that lacks cultivated representatives (Table S1 in Supplementary file 2). For NP metagenomes, 11 bins were assigned to the class Deltaproteobacteria, which was identified as the dominant taxon with a total bin abundance of 27.8 (the number of genome copies per million sequences). NP_bin1 was affiliated to Gammaproteobacteria with a bin abundance of 14.3. Two archaeal bins (NP_bin6 and NP_8) were assigned to the class Methanomicrobia with a total bin abundance of 7.4. NP_bin13 was affiliated to the order Methanomethylicales, with bin abundance of 1.9 (Fig. 3). In the SP metagenomes, the total abundances of bins affiliated to Betaproteobacteria (six bins), Deltaproteobacteria (eight bins) and Gammaproteobacteria (five bins) were 15.5, 13.6 and 5.9, respectively. SP_bin26 assigned to the order Burkholderiales was the most abundant bin (7.3). For archaea, four bins were affilitaed to the phyla Euryarchaeota, Thaumarchaeota and Candidatus Verstraetearchaeota. Phylogenetic analysis indicated that the dominant archaeal bin (SP_bin24, bin abundance 1.8) was affiliated to the genus Methanoculleus (Fig. 4).
Figure 3. Phylogenetic characterization of 25 moderate/high-quality MAGs recovered from sediments without submerged plants. The phylogeny was generated from orthology proteins using Orthofinder software. 63 microbial genomes using IQ-TREE with a best fit LG + F + R9 model using the Bayesian Information Criterion (BIC). Bootstraps were based on 1000 replicated trees. The alignment of 63 genomes was generated using MAFFT
Figure 4. Phylogenetic characterization of 31 moderate/high-quality MAGs recovered from sediments with submerged plants. The phylogeny was generated from orthology proteins using Orthofinder software. 66 microbial genomes using IQ-TREE with a best fit LG + F + R8 model using the Bayesian Information Criterion (BIC). Bootstraps were based on 1000 replicated trees. The alignment of 66 genomes was generated using MAFFT
To improve understanding of the uncultured archaeal phylum Verstraetearchaeota, a phylogenetic analysis of a potentially divergent cluster of mcrA genes from these Verstraetearchaeota genomes was conducted (Fig. S3 in Supplementary file 2). Results showed that sequences obtained from this study clustered with natural environmental clones. For example, mcrA from MAG NP_bin13 was close to sequences from hot spring sediments, and mcrA from MAG SP_bin4 clustered with sequences from freshwater spring sediments.