The CW system consists of four identical stages made from gray PVC plastic with each stage being 100 cm in height and 10 cm in diameter with a total treatment surface area of 0.0.032 m² (0.008 m² each) with 80 liters influent tank filled with 60 liters to feed the system and 80 liters effluent tank to collect the final effluent (Figure 1). The experiment was constructed outdoors at the School of Engineering, Cardiff University. The experiment was run throughout the varying seasonal conditions with the daily mean air temperatures ranging from 5–25°C between winter and summer periods respectively. Water samples temperature fluctuated less significantly ranging from 10–20°C between winter and summer periods, respectively. Gravel was used as the main substrate in all four stages. In each stage, a depth of 10 cm in the bottom layer was filled with coarse gravel (20–25 mm). This served as the supporting and drainage layer; the following layer was filled with gravel (4–9 mm) as the main substrate layer with a depth of 50 cm; a 10 cm top layer of gravel (10–19 mm) was added to facilitate the dispersion and the distribution of wastewater and the growth of plants. The porosity of the gravel was 40% combined for the bottom and main substrate layer. Each stage was planted with Phragmites Australis at the beginning of the experiment, and good growth with lush vegetation was observed after 2 months by feeding the system with synthetic wastewater. Phragmites Australis were chosen as these are the most commonly used reference plant for CWs in Europe and are able to survive in most conditions (Brix, 1994; Sun et al., 2005; Kadlec and Wallace, 2008; Vymazal, 2010). Moreover, these plants provide a comparatively high oxygen transfer into the rhizosphere, which facilitates the aerobic degradation of pollutants (Barbera et al., 2009; Wang et al., 2012). The system was fed with synthetic domestic wastewater prepared freshly each week via the peristaltic pumps from the influent tank. About 60 liters of synthetic wastewater were required to feed the system for 1 week. The system was operated with three batch cycles per day, with each cycle entailing 2 h saturation and 6 h unsaturation giving a total of 8 h per cycle. The synthetic wastewater simulates typical domestic wastewater with a high concentration of organic carbon source and nitrogen to obtain approximately 700 mg/l of COD and 60 mg/l of NH4+-N. Two liters of room temperature (circa 20°C) synthetic wastewater were pumped into the system in each cycle, totaling 6 liters per day being actively pumped into the system. The synthetic wastewater was batch loaded to the first stage and sequentially passed through the other stages, generating alternate wet/dry periods in individual stages. Before starting the experiment and loading the synthetic wastewater to the CW, the system was inoculated with activated flocs obtained from the aeration basin of a local domestic wastewater treatment plant for about 2 weeks to provide seed microorganisms for the system.

The experiment was divided into four phases: Phase 1 with only tidal flow, Phase 2 with tidal flow and step-feeding ratio (80:20), Phase 3 with tidal flow and step-feeding ratio (70:30), and Phase 4 with tidal flow and step-feeding ratio (60:40). Step-feeding ratios during experimental Phase 2 to 4 were distributed from the influent reservoir into the third stage of the CW as shown in Figure 1. In these phases, a portion of the flow will be allowed into the first and second stage as normal, and the remainder will be added directly from the reservoir into the system once the synthetic wastewater finishes passing through from the first and second stages to reach the third stage. Samples were collected once a week from the influent tank and the effluent of each stage and analyzed directly in situ for pH and temperature using a pH/EC/TDS meter (HANNA HI 991301). Chemical oxygen demand (COD), nitrite-nitrogen (NO₂⁻-N), nitrate-nitrogen (NO3--N), ammonium-nitrogen (NH4+-N) and total nitrogen (TN) were analyzed using a Hach DR/3900 spectrophotometer and digester in the laboratory.

A tidal flow strategy was generated in each stage using peristaltic pumps, which were controlled by specific programmable timers (Williamson Pumps Ltd & Williamson Manufacturing Company Ltd.; model: CM type variable speed cased pump). This process repeatedly allowed the mesocosm CWs to be filled with synthetic wastewater until the main media layer fully submerged (wet saturated conditions) and subsequently drained after a desired time (dry unsaturated conditions). Whilst, the additional oxygen availability resulted in improved nitrification it creates less favorable conditions for the denitrification step (Vymazal and Kröpfelová, 2011; Li et al., 2015). Therefore, a step-feeding strategy was adopted at later stages to improve the TN reduction (denitrification step) by introducing the influent synthetic wastewater to the nitrified liquid. Such approaches may provide a more efficient use of the available carbon source that could enhance the denitrification step (Miyaji et al., 1980; Fillos et al., 1996; Puig et al., 2004; Hu et al., 2012).

Gravel samples were taken from the lower (L) and upper (U) position of the main substrate from each of the four stages of the multistage CW system after 420 days (at the end of experimental Phase 4) and stored at −80°C. Genomic DNA was extracted from the gravel samples (1L, 1U, 2L, 2U, 3L, 3U, 4L, and 4U) using a Meta-G-Nome™ DNA Isolation Kit (Cambio Ltd). Gravel (10 g) was placed in a 20 ml sterile tube with 1 ml of 0.2% Tween 20 Wash Buffer (Cambio Ltd) and shaken for 10 min on a wrist action shaker at maximum speed (to remove the gravel microbial biofilm). The cell suspension was then transferred to a sterile 1.5 ml NoStick tube (Alpha Laboratories), centrifuged at 14,000 × g for 2 mins, and DNA extracted from the cell pellet according to the manufacturer's protocol. DNA was evaluated for quantity and size using a Qubit fluorometer (Qubit dsDNA BR Assay Kit; Invitrogen) and Agilent Tape Station (High Sensitivity D1000 ScreenTape and reagents; Agilent Technologies Inc) and stored at−20 °C until required for molecular analysis. Details for bacterial 16S rRNA gene PCR and Illumina sequencing are provided in Supplementary Text S1. Principal coordinates analysis (PCA) was used to visualize the difference in the microbial community calculated from the computed distance matrix and to visualize the performance characteristics of the system in terms of nutrient removal. R (v4.1.3) software was used to generate correlation and PCA analysis using packages corrplot (v0.92) and factoextra (v1.0.7) respectively (Alboukadel and Mundt, 2020; Wei and Simko, 2021).

The experiment was carried out over 420 days and evaluated the multistage CW system on nitrogen removal using both tidal flow (Phase 1) and tidal flow with step-feeding (Phases 2 to 4) as shown in Figure 2. COD removal rate was uniform across the experiments (ANOVA, p > 0.05) consistently reaching 90%, due to the system utilizing effective aerobic conditions induced by tidal flow operation (Figure 2A). However, a lower performance (70–75% COD removal) was observed in Phase 1 during a longer period with cold temperatures when water temperature approached 10°C. Significant COD removal obtained in all experiments (Phases 1 to 4) was predominantly due to enhanced oxygenation efficacy of the tidal flow system (Zhao et al., 2004; Hu et al., 2012; Chang et al., 2014) and the intensive microbial activities these conditions promote (Dušek et al., 2008). Previously, it has been reported that in most cases, insufficient oxygen supply is the main reason for poor biological COD removal in CW systems (Korkusuz et al., 2005; Ayaz et al., 2012; Wu et al., 2015). Nitrification is restricted even more than COD reduction because oxygen is utilized for carbon oxidization before nitrification due to the faster growth rate of heterotrophic organisms compared with that of nitrifiers (Wu et al., 2011; Bassin et al., 2015; Ge et al., 2015). However, well aerated systems provide unfavorable conditions for other process such as denitrification (Vymazal and Kröpfelová, 2011).

The average NH4+-N removal was 91 ± 8%, 89 ± 1%, 74 ± 5%, and 82 ± 5% for experimental Phases 1 to 4, respectively as shown in Figure 2B. TN removal was limited by denitrification (Figures 2C, D) and showed significant variability across the experiments (ANOVA, p < 0.05). The average removal of TN in Phase 1, where only tidal flow was employed, was 71%. After applying the step-feeding strategy alongside tidal flow in Phases 2, 3 and 4 the average TN removal was 66–81% as shown in Figure 2D. The nitrogen removal rate for Phases 2 and 4 were 77 gN/m3.d and 70 gN/m3.d respectively, higher than observed in Phase 1 (67 gN/m3.d). This is due to the hydraulic loading rate (HLR) and application of step-feeding which re-introduced external carbon source to enhance the denitrification process. Interestingly, only 63 gN/m3.d nitrogen was removed in Phase 3. Figure 2 shows that each time the new feeding ratio was introduced (Phase 2–4) the nitrogen removal was dropping and then gradually increasing with the duration of the applied step-feeding phase. This could indicate microbial community response and adaptation to new conditions when the nutrient availability was altered. The step-feeding ratio of 70:30 occurs as breakthrough point for nutrient distribution and availability in the system as it is clearly visible that the following Phase 4 although also experiencing initial drop, shows much faster rebound and consequent increase in the N removal over the studied period.

The operational condition applied in Phase 1 showed that, stages 1 and 2 were responsible for a majority (>85%) of COD removal (Supplementary Figure 1), hence depleting the remaining part of the system (stages 3–4) from the required source of carbon to deliver effective denitrification as indicated by the high nitrate concentration in the final effluent. This shows a need to improve carbon pool distribution in the system. Therefore, a step-feeding strategy was introduced in Phases 2–4 with varying feeding ratios split between stages 1 and 3. It can be seen however that the internal addition of a carbon source did not result in locally (stage 3) elevated COD concentration (Figure 3). This might suggest that the added carbon load was efficiently processed and quickly entered various microbial metabolic pathways including nitrogen cycling.

In Phase 1 (Figure 3A), the influent NH4+-N was reduced by 91% to 5.3 mg/l and the system showed significant (p < 0.05) increase of nitrification between stages 1–3 while plateauing and reaching similar (p > 0.05) levels at stages 3 and 4. The first two stages operated with 28 and 52% of nitrification rate respectively while NH4+ reduction at stage 3 and 4 reached 59 and 63%, respectively. Nitrate accumulation varied between the stages reaching 49% at the stage 4 and only 4% at the stage 2. This shows great disproportion of denitrification potential across the treatment stages of the system and a stepwise exhaustion of carbon sources from inflow to effluent of the system.

Upon introduction of step-feeding, nitrogen dynamics were significantly (ANOVA, p < 0.05) affected by the tested step-feeding ratios displaying marked differences both between the phases but also within internal, stage-wise dynamics. The influent NH4+-N was reduced by 89, 74, and 82% with the bulk of reduction, a third of reduction and slightly more than half of the reduction occurring in the first three stages in Phase 2 (Figure 3B), Phase 3 (Figure 3C), and Phase 4 (Figure 3D), respectively. For all three phases, nitrate accumulation started from stage 1 indicating that denitrification became the limiting process for TN reduction due to the carbon deficiency. In the Phase 2 of the experiment, step-feeding introduction of carbon source activated denitrification potential at the carbon dosing point (stage 3) and showed lowest NO3--N accumulation rate (9%) while at the same time delivering simultaneous very effective NH4+-N removal (77%) (Lai et al., 2020; Gupta et al., 2022). For step feeding ratios of 70:30 and 60:40 (Phases 3 and 4 respectively) internal introduction of carbon source showed to have initially inhibit nitrification rates (circa 15% NH4+-N reduction) and lead to increased nitrate accumulation (circa 60%) due to insufficient denitrification. This could be a sign of growing competition for carbon source between nitrogen cycling microbial consortia and the rest of bacterial community in the system as mentioned previously via TN and COD dynamics analysis. Overall, the differences in NO3--N accumulation in the system displayed denitrification capabilities of the individual stages coupled with step-feeding. High nitrate accumulation rates (45–49%) observed for final effluent at the Phases 1 and 2 were due to low inflow NH4+-N levels (>10mg/L) and relatively high residual NO3--N carried over from the previous stage (3 mg/L) in comparison to conditions observed at Phases 3 and 4 (25–40 mg NH4+-N /L). Nevertheless, despite different internal dynamics observed across all experiments (Phase 1–4) the differences in final effluent NO3--N concentration were insignificant (p > 0.05).

Details for statistical analysis results for the comparison of the mean removal efficiencies between different step-feeding ratios are provided in Supplementary Text S2.

The invested multi-stage TFCW operated in step-feeding strategy. The wastewater inflow was distributed between stage 1 and stage 3 of the system according to the 3 step-feeding ratios of 80:20, 70:30, and 60:40. Nitrogen removal performance was significantly different between the applied step-feeding ratios (ANOVA, p < 0.01). The system operated at a minimum of 60% TN removal, reaching maximum effectiveness of over 80% when applying different step-feeding ratio thus, step-feeding ratio has been proven to be a differentiating factor in achieving effective TN removal in the study.

However, the subsequently observed effectiveness variability between studied ratios was not expected and did not show a linear relationship (ratio vs. TN removal effectiveness). The minimum TN removal performance was obtained at 70:30 ration while both 80:20 and 60:40 ratios performed comparably well, with 60:40 ratio gradually reaching 80:20 ratio performance levels by the end of the monitoring period. A similar removal performance of 80% TN removal was achieved in a four-stage TFCW study optimizing step-feeding ratio and reported 80:20 step feeding ratio to outperform 90:10 and 85:10:5 dosing ratios (Hu et al., 2012). However, test in similar TFCWs showed only above 30% nitrification, including even lower TN removal at 80:20 step-feeding ratio (Yang et al., 2011). This indicates that in addition to step-feeding ratio optimization knowledge on carbon distribution and utilization in the system is required. Carbon to nitrogen (C/N) ratio is a key parameter in nitrogen cycle (Her and Huang, 1995; Ji et al., 2015). In wastewater treatment systems, carbon to nitrogen ratio is often reported in reference to carbon source using organic carbon fraction [i.e., biochemical oxygen demand (BOD), COD, total organic carbon (TOC) or dissolved organic carbon (DOC)] or total carbon fraction (TC). This study uses a COD parameter to represent the available carbon in the studied system. A complete denitrification requires a stoichiometric COD/N ratio equal 2.86 (Fu et al., 2009). However, practical considerations show that the effective nitrogen removal via denitrification occurs when the COD/N ratio reaches above 3. Nevertheless, operational conditions, quality of wastewater and carbon source composition requires denitrification process to be carried out at COD/N ratios often reaching 10 and above (Han et al., 2015; Pelaz et al., 2018; Deng and Shi, 2020). Meanwhile, anammox process becomes dominant and capable of delivering over 90% N removal when the COD/N ratio equals to 1 and lower (Wang et al., 2019; Sarvajith et al., 2020). In conventional CW designs (Vertical flow (VF)-, Horizontal flow (HF)-, Free water surface (FWS)-CW) operating with municipal wastewater, COD/N ratio ranging between 5–10 can deliver up to 50% nitrogen removal (Li et al., 2020; Rampuria et al., 2020; Zhu et al., 2021). High performing CWs specifically designed for N removal operate on higher COD/N ratios due to optimization of carbon dosage. Introduction of tidal-flow or spray aeration in CW can help to deliver up to 80% of nitrogen removal however a COD/N >6 is required (Zhi and Ji, 2014; Wang et al., 2020).

The investigated system showed a strong correlation (p < 0.05) between COD/N ratio and TN removal observed the final effluent at the stage 4 (Supplementary Figure S2). The highest performing Phase 2 (80:20 ratio) delivered >80% TN removal in the effluent, maintaining an average 4.9 COD/N ratio while Phase 3 and Phase 4 both showed C/N ratios varying between 2–3 indicating a potential shortage of carbon source in the system. Interestingly, the highest performing Phase 2 had least favorable thermal conditions, with water temperatures on average 6°C lower than Phase 4 (11.5°C and 18.1°C respectively). PCA analysis (Figure 4) reveals overlapping clusters between Phase 3 and 4 while Phase 2 performance data was grouped separately further indicating supreme performance of Phase 2. High performance effectiveness of 80:20 ratio is even more visible when considering variability of thermal conditions the system was exposed to Supplementary Figure S2 and which would affects kinetics of nitrogen removal (Pang et al., 2015; Myszograj and Bydałek, 2016). Nevertheless, the strong inverse correlation between water temperatures in the system (temp_out, Supplementary Figure S2) and COD/N ratio indicated that the inhibitive effect of low temperature can be overcome with higher carbon availability and conversely at higher temperatures, excessive carbon surplus is not required to sustain effective N removal in the investigated system.

A combined total of 1,079,014 16S rRNA gene sequence reads were obtained from the extracted DNA from the upper and lower positions of the gravel substrate from the four stages of the CW system after 420 days (Phase 4; 1L, 1U, 2L, 2U, 3L, 3U, 4L, and 4U). Read numbers ranged from 104,118 (3U) to 191,610 (1U) sequences per sample with an average read count of 134,877 reads. Rarefaction curves (Supplementary Figure S3) and Good's coverage statistics (Supplementary Table S5) indicate that 16S rRNA gene libraries for each gravel sample were sampled sufficiently to capture the majority of the bacterial diversity. Interestingly, both diversity (Shannon and Simpson diversity) and species richness (S_Chao1) indices (Hugerth and Andersson, 2017) increased in the CW system from stage 1 to stage 4 (Supplementary Table S5, Supplementary Figure S3), which demonstrates that bacterial diversity increased at each stage throughout the system. This agrees with other multi-stage CW systems that also showed an increase in bacterial diversity through the system from inlet to outlet (Babatunde et al., 2016; Rajan et al., 2019) and is corroborated by PCA UniFrac analysis (Supplementary Figure S5).

Bacterial 16S rRNA gene sequences were assigned taxonomy to identify the different bacterial communities at each stage and were classified from phylum to genus level. The relative abundance of a given bacterial group was set as the number of sequences affiliated with that group divided by the total number of reads per sample (Figures 5, 6; Supplementary Table S5, Supplementary Figure S4). In the CW system, a total of 14 bacterial phyla representing 62 assigned genera were identified at greater than 0.1% of the community. However, in all samples a large fraction of the sequences could not be assigned at the genus level (50 to 61.2%; Supplementary Figure S4a), and this presumably represents a large proportion of novel and unknown genera.

Overall, the CW system was dominated by members of the phyla: Proteobacteria (44.5%), Firmicutes (10.9%), Planctomycetes (13.6), Bacteroidetes (5.7%), Chloroflexi (3.8%), Verrucomicrobia (3.3%), and Acidobacteria (3.3%), and lesser proportions of Synergistetes (1.2%), Actinobacteria (0.53%), Ignavibacterae (0.36%), Nitrospirae (0.35%), Hydrogenedentes (0.12%), Armatimonadetes (0.11%), and Spirochaetes (0.10%), although there were some clear differences between phyla at different stages (Figure 5). The five most commonly found phyla in our study have been reported previously in CW systems and are thought to be key phyla in the successful operation of CW systems (Verduzo Garibay et al., 2021). The most abundant phylum Proteobacteria showed clear differences at the class level, with higher numbers of Gammaproteobacteria being found in stages 1 and 2 and Betaproteobacteria becoming more prevalent at stages 3 and 4 (Supplementary Figure S4b). Alphaproteobacteria and Deltaproteobacteria were found consistently throughout, while unclassified Proteobacteria were abundant in stage 1 but rapidly declined after stage 2.

Changes in bacterial diversity and community structure with stage and position may be linked to the different processes and rates of N removal occurring at each stage, coupled with the introduction of fresh carbon sources at stage 3. This is evident from the dominant assigned bacterial genera (Figure 6) and the top 20 16S rRNA gene operational taxonomic units (OTUs) (Supplementary Figure S6) identified as proxy for bacterial species. The dominant bacterial genera observed (Figure 6) in the whole CW system (representing >1% abundance) were Thiothrix, Planctomyces, Azonexus, Pseudoxanthomonas, Hydrogenophaga, and Gemmobacter, and all varied in abundance depending on CW stage. For example, the top three genera, Thiothrix ranged from 3.9–20.3% in stages 1 to 3 and but were <0.2% at stage 4. Conversely, Planctomyces were <0.7% in stages 1 to 3 but ranged from 10.4–17.1% at stage 4, while Azonexus steadily increased from 1.3–2.4% at stage 1 to 1.6–6.5% at stage 4. All dominant genera were represented in the top 20 OTUs, although the most common OTU (or bacterial species) was OTU0009, which was unassigned at the genus level but could be assigned at the family level and belonged to members of the Veillonellaceae.

Operational taxonomic unit OTU0009 was particularly prevalent at stages 1 to 3 and nearly absent at stage 4. The increased abundance of OTU008 in stage 3 and the upper portion of stage 4 correlates with the addition of the carbon source supplement as denitrification in CW systems by these organisms is known to be associated with high TOC (Wu et al., 2016). Pseudoxanthomonas are also known heterotrophic denitrifiers and OTU0045 was found throughout the system but particularly within stage 2. Other abundant denitrifiers found throughout included Gemmobacter (OTU0014) which can denitrify mixotrophically (Du et al., 2020) and Hydrogenophaga (OTU0056) which can utilize hydrogen (Willems et al., 1989).

Stage 4 was dominated by members of the Planctomycetes (Figure 5), including the genera Planctomyces and Pirellula (Figure 6), and other uncultured genera (Supplementary Figure S6). Many Planctomycetes conduct anammox metabolism, a process in which ammonia is oxidized by nitrite to nitrogen gas. Planctomyces OTU0059, and Planctomycetaceae OTUs (OTU0062, OTU0065, and OTU0080; Supplementary Figure S6) were well represented in the top 20 OTUs in high numbers at stage 4. Presumably very low oxygen (or anoxic) conditions and low concentrations of organic matter in the fourth stage allowed for these bacteria to proliferate and outcompete denitrifying Proteobacteria (Jin et al., 2012; Mosley et al., 2022).

However, further investigations are necessary to confirm the taxonomic identity of the putative anammox bacteria through further sequencing (16S rRNA and hydrazine synthase genes; Mosley et al., 2022) or anammox-specific probes (Tal et al., 2006).

Abundance of denitrifying Proteobacteria (e.g., Thiothrix, Thauera, Pseudoxanthomonas and Hydrogenophaga) in the CW system in stages 1 to 3 suggests that this is the major N removal process with both heterotrophic and autotrophic denitrification occurring. Evidence of autotrophic ammonia-oxidizing bacteria (Nitrosomonas and Nitrosospira) and autotrophic nitrite-oxidizing bacteria (Nitrospira and Nitrolancea) occurred throughout the system, although at very low abundance (<1.0%; Supplementary Figure S7) demonstrates that ammonia was continually oxidized to nitrate. The high abundance of Planctomycetes and the potential for anammox bacteria suggests that the anammox process was an alternatively occurring pathway for nitrogen removal, especially at stage 4 of Phase 4 when conditions were suitable. Clearly, nitrogen removal in the CW system was accomplished via diverse pathways, including autotrophic nitrification, heterotrophic denitrification, autotrophic denitrification, and possibly anammox. Similar collaborative microbial pathways for N removal have been found routinely in constructed wetlands (Wei et al., 2021; Zhang et al., 2021) and it is the abundance, consortia and distribution of these organisms that is key to the performance of a CW system (Zhang et al., 2021).