The experiment was conducted at the Jinjing Catchment, Changsha, Hunan Province, China (28°30′ N, 113°18′ E). This site has a subtropical monsoon climate with an annual average rainfall of 1330 mm and an annual average temperature of 17.5°C. A surface-flow CW was established in March 2014, three cells were set up in parallel (three replicates). Each cell was divided into three identical segments (length, 5 m; width, 2 m; and water depth, 0.2 m), named CW1, CW2, and CW3 (Supplementary Figure 1). M. aquaticum was planted on the wetland surface at an initial density of 3 kg m^-2. The cell without planting the M. aquaticum was performed as the control. The influents from a storage tank and a fresh water tank were mixed in a settling pond installed in front of the CW, and the effluent of the CW was applied to the land. The characteristics of the swine wastewater are listed in Supplementary Table 1. Each cell was operated in an intermittent flow regime with a total of 0.18 m³ d^-1 wastewater, and the hydraulic retention time was 11 days in each segment.

Water samples were collected two or three times per month and were analyzed immediately at the laboratory. The dissolved oxygen (DO) and pH were measured in the field by a multi-parameters water quality monitoring instrument (HACH, HQ30D, United States). Chemical oxygen demand (COD) was measured by a standard potassium dichromate titration method (Federation and Association, 2005). NH4+-N, NO2−-N, NO3−-N, and total N (TN) concentrations were measured by a continuous flow analyzer (AA3; Seal Analytical, Norderstedt, Germany). Microbial samples were sampled at 24 ± 2°C during the stable operation periods. Water (0–5 cm depth) and sediment (0–5 cm depth) samples were taken in the second segment of the CW (CW2, 5–10 m from the inlets), which exhibited the highest N removal efficiency during the monitoring period. In September 2015, sediment cores were taken in 10 replicates from within monospecies of M. aquaticum and from unvegetated sediment for molecular analysis. Unvegetated sediment samples were obtained from the upper 5 cm, which included the oxic–anoxic interface. Vegetated sediment samples were obtained from all over the root zone by shaking off sediment that adhered to the roots. All the samples were stored in an ice cooler and brought back to our laboratory. The water samples were filtered onto a 0.22-μm cellulose membrane, and the sediments samples were freeze-dried. Samples for DNA extraction were stored at -80°C until use.

Total DNA were extracted using the FastDNA SPIN Kit for Soil (MP Biomedical, Santa Ana, CA, United States) according to the manufacturer’s protocol and then it was visualized using 1% agarose gel electrophoresis and stored at -80°C until use.

Quantification of V4 region of the rrs gene, which encodes the 16S rRNA, archaeal and bacterial amoA, nirK, nirS, and nosZ, were performed on a real-time PCR system (Bio-Rad, Hercules, CA, United States) using SYBR green as a fluorescent dye. The thermal cycling conditions and primers used for each reaction are described in Supplementary Table 2. Gene copy numbers in unknown samples were determined based on standard curves obtained from 10-fold serial dilutions of plasmids containing the target genes. The employed results with correlation coefficient and amplification efficiency were greater than 0.98 and 98%, respectively. All samples were run on an agarose gel after the reaction to confirm the size of amplicons, and the specificity of the amplification products was confirmed by a melting curve analysis.

The V4 region of the bacterial 16S rRNA gene was amplified using the primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) (Caporaso et al., 2012). Each pair of primers used to amplify a certain water or sediment sample was barcoded with a sample-identifying 12-base barcode on both the forward and reverse primers. The 50 μL reaction system contained 20–30 ng of DNA, 1.5 μL of each primer (at a final concentration of 0.3 μM), and 25 μL of Premix Ex Taq (TaKaRa, Dalian, China). The PCR conditions, performed in a PCR instrument (Eastwin, China), consisted of an initial denaturation step at 94°C for 1 min, followed by 30 cycles at 94°C for 20 s, 57°C for 25 s and 72°C for 30 s, and a final extension step at 72°C for 10 min. The PCR amplicons were purified with the Gel Extraction Kit (D2500-02, Omega Bio-Tek, Norcross, GA, United States), and DNA concentrations were quantified by a spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States). True positive amplicons were combined equally, and a DNA library was obtained according to the MiSeq Reagent Kit Preparation Guide (Illumina, San Diego, CA, United States), and then the DNA library was sequenced using the Illumina MiSeq platform according to the manufacturer’s instructions. All the raw sequences have been deposited into NCBI Short Read Archive (SRA) with accession no. SRP095281.

After assigning each sequence to its sample according to its barcode, allowing up to one mismatch, a total of 675,089 reads from both ends were obtained as a partitioned run for the 40 samples. Paired-end reads of sufficient length (those with at least a 30 bp overlap) were combined into full-length sequences with targeting 253 bp as average fragment length by FLASH (Magoč and Salzberg, 2011). Thereafter, quality trimming was performed using Btrim with a quality score of 20 as cutoff and 5 as window size (Kong, 2011). After removing the sequences with ambiguous bases, the sequences with lengths out of 245–260 bp were removed. Thereafter, UPARSE program (Edgar, 2013) was used to remove the chimeras and classify the sequences into operational taxonomy units (OTUs) at the 97% similarity level. The taxonomic annotation of individual OTUs was performed through Ribosomal Database Project (RDP) classifier (Wang et al., 2007) with a minimal 50% confidence score. Finally, random re-sampling was performed with 10,000 sequences per sample, and this resampled OTU summary table was used for further statistical analysis (Supplementary Data Sheet 2).

One-way analysis of variance (ANOVA) followed by Student–Newman–Keuls test was used to check for quantitative differences between treatments using Statistics for Windows version 20.0 (IBM Corp., Armonk, NY, United States). The dataset generated by 16S rRNA gene sequencing was further analyzed with all of the following statistical methods. (1) α-diversity comparison analysis was performed using R vegan package. (2) Comparison tests between two groups were based on unpaired Student’s t-test in Excel software. (3) The canonical correspondence analysis (CCA) was used to evaluate the linkages compositions and environmental factors (Legendre and Legendre, 2012). To select environmental factors in the CCA modeling, we used variation inflation factors to examine whether the variance of the canonical coefficients was inflated by the presence of correlations with other attributes. If a factor had a variation inflation factor value greater than 20, we deemed it to be affected by other attributes and consequently removed it from the CCA biplot (Yang et al., 2014). Multivariate testing, based on 999 Monte Carlo permutations, confirmed the significance of the two canonical axes (p ≤ 0.001). The importance of the environmental attributes was verified by Mantel tests with 999 permutations based on their significant correlation (p ≤ 0.001). (4) Phylogenetic molecular ecological networks (pMENs) were constructed by regularly a random matrix theory (RMT)-based interface approach to elucidate the microbial interactions in the sediments and water with and without M. aquaticum. pMENs were generated based on a resampled OTU table at Molecular Ecological Network Analysis Pipeline (Zhou et al., 2010, 2011; Deng et al., 2012)¹. The constructed pMENs were visualized using Cytoscape 3.30 software.

During the operation (600 days; Figure 1), the three-stage CW achieved COD removal efficiency of 88.2% under an average inflow concentration of 484 mg L^-1. More than 75% of removal occurred in the first-stage wetland, which had the highest removal efficiency of 66.5%. The average inflow NH4+-N concentration was 107 mg L^-1, accounting for more than 70% of the N inputs to the wetland. The average NH4+-N removal efficiency was 98.3%. The average TN concentration of 150 mg L^-1 in the wastewater inflow was reduced to an average effluent level of 5.9 mg L^-1, which achieved removal efficiency of 95.8%. The highest removal efficiencies of NH4+-N and TN (78.6 and 71.6%, respectively) were achieved in the second-stage CW. These performances were stable during the operation periods, even with the high average COD, NH4+-N, and TN removal efficiencies of 92.3, 99.5, and 97.9%, respectively, in winter.

The V4 region of the rrs gene, which encodes for the 16S rRNA, archaeal and bacterial amoA, and nirK/nirS/nosZ genes, was quantified by qPCR to demonstrate the abundances of total bacteria, ammonia-oxidizing archaea (AOA), ammonia-oxidizing bacteria (AOB), and denitrifiers. In the sediment, the results showed that the abundance of the bacterial 16S rRNA gene was increased significantly from 1.61 × 10⁷ to 1.44 × 10¹⁰ copies per gram dry soil (862-fold) by M. aquaticum (Figure 2). The abundance of archaeal amoA gene increased significantly from 8.85 × 10⁴ to 1.93 × 10⁷ copies per gram dry soil (218-fold), and the abundance of the bacterial amoA gene increased from 6.65 × 10⁴ to 5.75 × 10⁶ copies per gram dry soil (86-fold). The archaeal amoA gene abundance was more increased than the bacterial amoA abundance, and it outcompeted the bacterial amoA gene in all sediments containing M. aquaticum. However, there were no significant differences between the abundances of archaeal and bacterial amoA genes in the unvegetated sediment. For denitrifying genes, M. aquaticum significantly increased the abundance of nirK gene from 2.03 × 10⁵ to 2.16 × 10⁸ copies per gram dry soil (1,063-fold), while the abundance of nirS gene increased significantly from 7.73 × 10⁴ to 3.58 × 10⁸ copies per gram dry soil (4,620-fold). In the unvegetated sediment, the nirK outcompeted nirS gene. However, a greater abundance of the nirS gene, compared with that of the nirK gene, was detected in sediments containing M. aquaticum. M. aquaticum significantly increased the abundance of nosZ gene from 6.84 × 10⁵ to 8.31 × 10⁸ copies per gram dry soil (1,210-fold), while in the water (Supplementary Figure 2), M. aquaticum significantly decreased the abundance of the bacterial 16S rRNA gene from 1.66 × 10¹² to 3.79 × 10¹⁰ copies per liter water (44-fold), and the bacterial amoA decreased by one order of magnitude. The abundances of archaeal amoA and denitrifying genes in the water showed no significant response to M. aquaticum (p > 0.05).

MiSeq high-throughput sequencing of the bacterial 16S rRNA gene was performed to explore the effect of M. aquaticum on the structure of functional communities. In the unvegetated zone, α-diversity indices (Supplementary Table 3) did not differ between sediments and water (p > 0.05). The sediments containing M. aquaticum showed significantly higher bacterial species richness compared with the unvegetated sediment. The composition of bacterial communities in the sediment with M. aquaticum differed significantly from those in the unvegetated sediment based on a principle coordinates analysis (Supplementary Figure 3).

The relative bacterial community abundances were analyzed at the phylum level (Supplementary Figure 4). Firmicutes was the most abundant phylum in the unvegetated sediment and water, accounting for 43.1 and 50.6%, respectively, of the total effective bacterial sequences. In the sediment and water containing M. aquaticum, Proteobacteria was significantly enriched (p < 0.001), and it was the most abundant phylum, accounting for 23.73 and 46.65%, respectively, of the total effective bacterial sequences. This result is similar to the previous findings of bacterial communities in constructed and natural wetlands, in which Proteobacteria was also the most dominant community (Wu et al., 2016). Within the Proteobacteria, Alpha-, Beta-, and Gamma-Proteobacteria were all enriched substantially in sediment and water containing M. aquaticum (Supplementary Figure 5). The other dominant phyla in the sediment containing M. aquaticum were Firmicutes (13.1%), Acidobacteria (11.9%), Planctomycetes (11.0%), Chloroflexi (11.5%), Actinobacteria (9.1%), Euryarchaeota (3.1%), Bacteroidetes (1.9%), and Nitrospirae (0.6%).

To further validate the functional populations involved in N removal that are associated with M. aquaticum, a heat map of major genera was generated, as illustrated in Figure 3. For nitrifiers, the AOB Nitrosospira was enriched significantly by M. aquaticum in the sediment (p < 0.001). Two AOB Nitrosomonas and Nitrosococcus were also enriched. AOB in the vegetated sediment was dominated by Nitrosospira. The abundance of the AOA Nitrososphaera was significantly higher in the M. aquaticum-containing sediments (p < 0.001). Nitrososphaera was dominant among ammonia oxidizers. Two NOB Nitrospira and Nitrobacter were not detected in unvegetated sediment, but were significantly increased in the vegetated sediment (p < 0.001). Moreover, M. aquaticum significantly increased the abundances of the anammox bacteria Candidatus Brocadia and Candidatus Kuenenia (p < 0.001), while these genus were not detected in the unvegetated sediment. For the denitrifying groups, abundant denitrifiers, such as Bradyrhizobium, Rhodoplanes, Bacillus, Zoogloea, Rhizobium, Acidovorax, and Hyphomicrobium, were enriched significantly in the vegetated sediment (p < 0.001). Bradyrhizobium and Rhodoplanes were dominant among the denitrifiers. In the water containing M. aquaticum, among nitrifiers, the abundances of Nitrosococcus and Nitrosospira increased, whereas the abundances of Nitrosomonas and Nitrososphaera decreased. In addition, sequences were not obtained for NOB and anammox bacteria in the water with and without M. aquaticum. For the denitrifying group, a few denitrifiers such as Rhizobium, Rhodoplanes, and Rhodobacter were enriched significantly by M. aquaticum (p < 0.001).

Compared with the unvegetated sediments, NH4+-N, pH, and total organic carbon (TOC) decreased significantly (p < 0.001), whereas NO2−-N increased, in the sediments containing M. aquaticum (Table 1). CCA was used to examine the relationship between the environmental variables and the microbial community composition. The results of the CCA, a significant model at the confidence level of p ≤ 0.001, indicated that sediment NH4+-N and TOC were important environmental factors controlling the microbial community structure (Figure 4). The importance of these environmental attributes was verified by Mantel tests based on their significant correlation (p ≤ 0.001) as well (Table 2).

Interactions among different microbial species in the CWs were calculated using a relative abundance matrix of genus-level OTUs with 97% similarity, and they were constructed in a network to describe the effect of M. aquaticum on the microbial interactions (Figure 5 and Supplementary Figure 6). The structure of the identified networks with and without M. aquaticum differed substantially in terms of their network composition, interaction patterns (positive, zero, or negative), and node overlap, which demonstrated that M. aquaticum greatly altered the interactions among different microbial species. A small and simple connected module of nodes corresponded to N-cycling members in the unvegetated sediments, primarily including ammonia-oxidizing Nitrososphaera and denitrifiers such as Acinetobacter and Thauera, and these species did not show any interaction. However, in the sediment containing M. aquaticum, a large and densely connected module of nodes corresponded to members involved in N removal, primarily including ammonia-oxidizing Nitrososphaera, nitrite-oxidizing Nitrospira, anaerobic ammonia-oxidizing Candidatus Kuenenia and denitrifiers such as Rhodoplanes, Bradyrhizobium, and Hyphomicrobium. These groups were distributed widely in the vegetated sediments, and they showed positive ecological associations with different species. The positive interactions between the AOA Nitrososphaera and denitrifiers were stimulated by M. aquaticum. Nitrososphaera showed positive interactions with the anaerobic denitrifier Rhodoplanes. Heterotrophic denitrifiers, such as Rhodoplanes, Bradyrhizobium, and Hyphomicrobium, were enriched and distributed widely in the sediments containing M. aquaticum. Moreover, microorganisms carry out fermentation including Clostridium and Smithella were also enriched in these sediments (Ezeji et al., 2007; Wawrik et al., 2016). The enriched denitrifiers showed positive relationships with fermentative bacteria. For example, Rhodoplanes was frequently found to have positive interactions with the anaerobic fermentative bacteria Clostridium. Moreover, positive relationships among the denitrifiers were observed. However, in the vegetated water (Supplementary Figure 6), the most abundant genus, Rhizobium, displayed co-occurrence through positive correlations with some denitrifiers or potential denitrifiers, including Rhodoplanes, Bradyrhizobium, and Rhodobacter.