Active predation, phylogenetic diversity, and global prevalence of myxobacteria in wastewater treatment plants

Myxococcota and Bdellovibrionota were active constituents of activated sludge microbiota

To explore the predating activity and diversity of predatory bacteria in activated sludge, ¹³C-labeled Escherichia coli and Pseudomonas putida cells (determined as 97.09 and 97.30 atom% ¹³C, respectively) were added to the sludge microcosms for maximumly eight days of incubation, and ¹³C incorporation was examined using rRNA-SIP to identify prokaryotic and eukaryotic microorganisms involved in actively consuming the ¹³C-labeled prey cells. Bacterial 16S rRNA gene amplicon sequencing-based analysis indicated the relative contribution of 47.9% and 42.7% of the obtained sequences by the added biomass upon amendment in the ¹³C-E. coli (Fig. 1A) and ¹³C-P. putida (Fig. 1B) microcosms, which dropped below 1.0% after 16 h and eight days of incubation, respectively. The overall bacterial community structure at the steady state was highly comparable to that of the control microcosms (Fig. 1C), indicating that the prey cell amendments did not induce too strong fluctuation in the microbiota structure during the SIP experiment that prevented disentangling the indigenous community dynamics.

The overall prokaryotic communities were obtained by 16S rRNA gene amplicon sequencing of the total DNA from the activated sludge microcosms amended with ¹³C-E. coli (A) and ¹³C-P. putida (B) cells, and the control group (C) without amendment. The structure of the active prokaryotic communities was inferred based on amplicon sequencing of the light rRNA fractions from the microcosms amended with ¹³C-E. coli (D) and ¹³C-P. putida (E) cells. The temporal change in the proportion of produced ¹³CO₂ in total CO₂ indicated the mineralization of the ¹³C-labeled cells of E. coli and P. putida in duplicate microcosms (F). Relative sequence abundance of the ten most abundant prokaryotic phyla, together with the genera Escherichia-Shigella and Pseudomonas, was shown.

The metabolically active bacterial communities, as inferred by 16S rRNA gene transcripts of the light rRNA fractions from the microcosms, were rather consistent throughout the experiment (Fig. 1D, E), but they showed clear compositional differences compared to the overall prokaryotic communities inferred by 16S rRNA gene amplicon sequences (Fig. 1A, B). Myxococcota and Bdellovibrionota species showed an average relative abundance of 17.5 (±0.7) % and 2.7 (±0.2) % in the 16S rRNA gene transcripts, respectively, which were significantly higher than 5.4 (±0.6) % and 1.3 (±0.1) % in the 16S rRNA genes of bacterial communities (p < 0.001, Mann–Whitney U test), suggesting their relatively high metabolic activity in the microcosms. The same activity pattern was observed for Cyanobacteria (3.8% vs. 0.4%), Proteobacteria (37.4% vs. 28.4%), and archaeal Halobacterota (8.6% vs. 0.2%). In contrast, members of Acidobacteriota, Bacteroidota, Chloroflexi, Nitrospirota, and Patescibacteria were much less abundant (average rRNA /rRNA gene ratio: 0.10–0.49) in the active communities compared to in the overall microbiota. The relative proportion of ¹³CO₂ derived from the ¹³C-labeled biomass peaked already in 16 h or 1 d after cell amendment, indicating the rapid flow of the added biomass ¹³C into the microbial food web (Fig. 1F).

For the micro-eukaryotic communities, the added prey bacterial cells, as expected and afore-demonstrated for bacterial communities, did not exert large change in community structure inferred from 18S rRNA gene amplicon sequences, whether for ¹³C-E. coli or ¹³C-P. putida (Supplementary Fig. S2). Several micro-eukaryotic groups were metabolically active members in the whole community. For example, Ciliophora contributed to 55.2–78.4% of the 18S rRNA gene transcripts in all the microcosms (Supplementary Figs. S2D, E), while they were much less abundant (6.5–37.0% of the 18S rRNA gene sequences) in the whole micro-eukaryotic communities (Supplementary Figs. S2A, B). By contrast, fungal Cryptomycota accounted for 37.3–60.6% in the micro-eukaryotic communities (Supplementary Fig. S2A, B), while they constituted only 1.7–4.6% of the 18S rRNA gene transcripts (Supplementary Figs. S2D, E).

Myxococcota dominated ¹³C-biomass incorporators in activated sludge

To predict predator-prey interactions in activated sludge, the degree of ¹³C-labeling of each prokaryotic or micro-eukaryotic taxon was quantitively determined relative to the ¹²C-control group by calculating the taxon-specific enrichment factor (EF) (see Method for calculation). Most of the ¹³C-labeled prokaryotic ASVs (as a proxy for species) belonged to Myxococcota, followed by Bdellovibrionota (Supplementary Fig. S3A). Among the bacterial genera with relative abundance >1% in the ¹³C-heavy fractions, strong ¹³C-labeling was found for the as-yet-uncultivated myxobacterial mle1-27 clade (average EF 2.66 across time and treatments), which contributed to 10.3% to 38.9% of the 16S rRNA gene transcripts in the ¹³C-heavy fractions, indicating its high metabolic activity in consuming the ¹³C-labeled biomass of both E. coli and P. putida. Comparatively, Haliangium spp. and uncultured Polyangiaceae belonging to Myxococcota, as well as the as-yet-uncultivated OM27 clade belonging to Bdellovibrionota, also exhibited strong ¹³C-labeling (maximum EF across time: 2.4–39.5), but almost exclusively in the microcosms amended with ¹³C-E. coli cells (Fig. 2A). The as-yet-uncultivated myxobacterial VHS-B3-70 clade exhibited the strongest enrichment (average EF 16.67 across time and treatments) but made up only 0.2% to 2.3% of 16S rRNA gene transcripts of the ¹³C-heavy fraction (Fig. 2A). Overall, our microcosm experiment tracking added ¹³C-labeled prey bacterial cells with rRNA-SIP suggested prominent predatory activity of Myxococcota and Bdellovibrionota lineages including largely as-yet-uncultivated ones (e.g., the mle1-27, VHS-B3-70, and OM27 clades) in activated sludge.

¹³C-labeled prokaryotic (A) and micro-eukaryotic (B) genus-level taxa were identified by SIP in the microcosms added with E. coli and P. putida after one, two, and four days of incubation. Enrichment factor (EF) was calculated for microorganisms using heavy and light rRNA gradient fractions of the ¹³C- and ¹²C-microcosms to infer ¹³C-labeling. Taxa with an EF > 0.1 in at least one of the treatment groups at one sampling time point was considered labeled. The area of circles indicates the relative sequence abundance of the labeled taxa in heavy ¹³C-rRNA. The negative EFs and positive EFs < 10^−1.5 were arbitrarily set to 10^−1.5 for a neat graphical display. Heatmap of the ¹³C-incorporating prokaryotic (C) and micro-eukaryotic ASVs (D) with relative sequence abundance >1% in the heavy rRNA fractions of at least one of the ¹³C-E. coli and ¹³C-P. putida microcosms at a sampling point.

Myxococcota and Bdellovibrionota predated more selectively than protists

For the micro-eukaryotes, several taxa belonging to Ciliophora, especially Cyrtophoria spp. and Telotrochidium spp., and also Peritrichia spp., Vaginicola spp., Aspidisca spp., and Epistylis spp., were highly enriched (maximum EF across time and treatments: 0.9–6.7) in the ¹³C-heavy rRNA fractions (Fig. 3B), in agreement with the dominance of Ciliophora in the micro-eukaryotic rRNA gene transcripts (Fig. 2B). The Candida–Lodderomyces clade and Cyberlindnera–Candida clade within Ascomycota, Magnoliophyta spp. within Phragmoplastophyta, and Poteriospumella spp. and unclassified Chromulinales within Ochrophyta were also strongly labeled (maximum EF: 13.5–242.5, Fig. 2B). Moreover, the ¹³C-biomass incorporation by micro-eukaryotes was independent of whichever prey bacteria (Fig. 2B, D), revealing no detectable prey preference in the metabolically active micro-eukaryotic predators. On the contrary, differential labeling by ¹³C-E. coli and ¹³C-P. putida cells was frequently observed for the predatory bacteria (Fig. 2A, C). The most obvious example was the OM27 clade ASVs belonging to Bdellovibrionota, which were found to incorporate ¹³C-labeled biomass exclusively of E. coli (Fig. 2C). Comparatively, Haliangium-affiliated ASV27 and ASV63 were labeled only by ¹³C-E. coli, ASV57 labeled by both ¹³C-E. coli and ¹³C-P. putida, while ASV72 and ASV76 were also labeled by ¹³C-P. putida, but only at a later sampling point (Fig. 2C). These results on the divergent labeling patterns with the tested prey bacteria together strongly implied population-specific predating behaviors of predatory bacteria in activated sludge.

The abundance of the abundant genera belonging to Myxococcota and Bdellovibrionota in aerobic and anaerobic sludge were compared according to amplicon sequencing-based analysis of bacterial 16S rRNA gene V3-V4 region. The top 10 abundant genus-level taxa across samples collected from eight samplings are shown, with the putative predators identified by SIP in the microcosm experiment highlighted. The asterisk denotes significant difference in relative abundance between aerobic and anaerobic sludges (p < 0.05, paired samples Wilcoxon test, n = 8).

Diverse Myxococcota and Bdellovibrionota inhabited wastewater treatment plants

The prominent activity of Myxococcota and Bdellovibrionota species in consuming the added ¹³C-labeled biomass in the microcosm experiment led us to examine their in situ abundance and diversity in activated sludge of local and global WWTPs. First, we profiled activated sludge and anaerobic sludge microbiota in the aeration and anaerobic tanks at a local WWTP (WWTP01) over two years and surprisingly found that Myxococcota (6.5 ± 1.3%) accounted for 5.7 times higher relative abundance than Bdellovibrionota (1.0 ± 0.2%) in the aerobic activated sludge samples (Fig. 3A). In total, ten myxobacterial genera were found with an average relative sequence abundance >0.1% in the activated sludge of WWTP01, including the putative predators identified in the microcosm experiment, i.e., Haliangium spp. (2.8 ± 0.7%) which represented the most abundant myxobacterial lineage in the activated sludge, uncultured Polyangiaceae (0.4 ± 0.1%), and the mle1-27 clade (0.2 ± 0.0%; Fig. 3). Moreover, Pajaroellobacter (1.2 ± 0.2%), Nannocystis (0.4 ± 0.1%), Phaselicystis (0.3 ± 0.1%), and several other myxobacterial clades, although not identified as putative predators in the microcosm experiment, were among the abundant myxobacteria in situ in the activated sludge. Although the myxobacterial genera showed comparable relative abundance in the anaerobic tanks, fed by returned activated sludge, to their counterparts in the aerobic tanks, the obligately aerobic myxobacteria were presumably metabolically inactive in the anerobic sludge. Unlike Myxococcota, members of Bdellovibrionota altogether showed significantly higher relative abundance in the aerobic sludge (1.0 ± 0.2%) than in the anaerobic sludge (0.6 ± 0.1%, paired samples Wilcoxon test p < 0.05, n = 8) (Fig. 3). The OM27 clade of Bdellovibrionota identified as transcriptionally active was abundant predators (0.5 ± 0.2%) in the aerobic activated sludge. Bdellovibrio spp. (0.5 ± 0.1%) represented another abundant genus belonging to Bdellovibrionota, although it was not enriched by ¹³C-labeled biomass tested in this study.

To reveal the diversity of myxobacteria residing in global wastewater treatment systems, we next constructed an integrative myxobacterial 16S rRNA gene dataset by combining the full-length 16S rRNA gene sequences assigned to Myxococcota that were generated from the local WWTP01 (507 representative OTU sequences) according to the SILVA taxonomy with sequences from the MiDAS 4 (2 521 sequences clustered to 1 010 OTUs), a well-established 16S rRNA gene reference database for bacteria in WWTPs globally. Phylogenetic analysis showed that Haliangiaceae and Polyangiaceae, the most diverse [49] and abundant [50] families in soil and other habitats, were well represented in wastewater treatment systems (Fig. 4A). The as-yet-uncultivated mle1-27 clade as well as members of four families, i.e., Myxococcaceae, Sandaracinaceae, Nannocystaceae, and Phaselicystidaceae, also contributed significant proportions of the myxobacterial diversity in global WWTPs (Fig. 4A).

Phylogenetic tree inferred from full-length 16S rRNA gene sequences within the phylum Myxococcota originated from global WWTPs (A) and the myxobacterial genera with presumed predatory behavior therein (B). The tree was constructed with 507 representative OTU sequences from a local WWTP (WWTP01) and 1 010 representative OTU sequences from the MiDAS 4 database [30]. The inner-ring color and labels indicate family classification with the SILVA database [46]. The external bars showed the average relative abundance of the OTUs at the WWTP01 across eight aerobic activated sludge samples collected over two years (maximum 0.07%). The scale bar corresponds to 0.1 substitutions per nucleotide position. Red-colored branches indicate sequences ≥94.5% identical to 16S rRNA gene sequences of myxobacteria with reported predatory behavior (suggesting potential predators at the genus level), the number of which was summarized for each myxobacterial genus in B. The numbers above bars indicate the proportion of the potential predators in each genus. Phylogenetic trees were inferred for Haliangium (C), the mle1-27 clade (D), and Pajaroellobacter (E) across different habitats. The trees were constructed with representative OTU sequences from WWTP01 (174, 134, and 69 sequences, for Haliangium, the mle1-27 clade, and Pajaroellobacter, respectively) and sequences obtained from the SILVA database [46] with known isolation source (507, 49, and 102 sequences, respectively).

The MiDAS and WWTP01 OTUs were classified as potential predators in this study if their full-length 16S rRNA gene sequence showed at least genus-level identity (i.e., ≥94.5%) to those of ly verified predatory isolates of Myxococcota (none of them from wastewater treatment systems, see a full list in Supplementary Table S2). The identified potential predators in global WWTPs mostly belonged to Polyangiaceae (130 out of 209 OTUs), followed by Haliangiaceae (37), Nannocystaceae (33), Myxococcaceae (6), and Phaselicystidacea (2) (Fig. 4A, B). Among them, 37 MiDAS OTUs classified as potential predators were taxonomically assigned to Haliangium, which was metabolically identified as putative predators in our microcosm experiment, substantiating their important roles in the microbial food webs in global WWTPs. Moreover, all the 173 Haliangium-assigned OTUs of WWTP01 were only remotely related to (maximumly 94.3% sequence identity of full-length 16S rRNA gene) myxobacterial predators reported in literature (Fig. 4B), suggesting the existence of possibly novel predatory Haliangium species and ecotypes in this local Chinese WWTP. Although members of Nannocystis, Polyangium, and Pajaroellobacter were not enriched by ¹³C-E. coli and ¹³C-P. putida cell amendment, these three myxobacterial genera were predicted as yet-to-be-verified predators in WWTPs based on full-length 16S rRNA gene sequence identity (94.8–98.1%) to known myxobacterial predators (Fig. 4B). Further phylogenetic analysis applied to myxobacterial 16S rRNA gene sequences from the local WWTP01 combined with sequences from the SILVA database originated from different environments (Fig. 4C–E and Supplementary Fig. S4) revealed that myxobacteria residing in WWTPs tended to form unique clades, phylogenetically distinct from those isolated from soil or other environments, especially for Haliangium, mle1-27, and Nannocystis, suggesting habitat filtering on myxobacteria in activated sludge.

For Bdellovibrionota, Bdellovibrionaceae and Bacteriovoracaceae species constituted important fractions of their diversity in global WWTPs (310 and 167 out of 1 338 OTUs, respectively; Supplementary Fig. S5), while over half of the collected sequences were assigned to the as-yet-uncultivated 0319-6G20 clade (Supplementary Fig. S5) which warrants future in-depth exploration.

Myxobacteria were ubiquitous in activated sludge across global wastewater treatment plants

To depict the biogeographic distribution patterns of predatory bacteria (e.g., members of Myxococcota and Bdellovibrionota) and explore their ecological roles and environmental drivers in activated sludge on a global scale, we exploited a 16S rRNA gene amplicon dataset of 1 186 activated sludge samples from 269 WWTPs worldwide published by the Global Water Microbiome Consortium [2]. To facilitate effective cross-study comparison, bacterial taxonomy was re-assigned to the representative sequences of OTUs using the SILVA SSU database (version 138), following the same bioinformatics procedures and cutoffs as applied in our study. Overall, we found that myxobacteria were ubiquitously detected (in 1 183 samples) in activated sludge of global WWTPs, constituting a non-negligible proportion (5.4 ± 0.1%) of activated sludge microbiota, and being consistently more abundant than Bdellovibrionota (1.1 ± 0.0%; Fig. 5A). While both Myxococcota and Bdellovibrionota lineages showed the lowest average relative abundance in the South America (2.7 ± 0.3% and 0.7 ± 0.1%, respectively; n = 80), they were the most abundant in North America (6.5 ± 0.2%, n = 616) and Africa (1.5 ± 0.1%, n = 36), respectively, suggesting a geographic distribution (p < 10⁻¹⁰ for both, Kruskal–Wallis test; Fig. 5A). Similar to the observation at the local WWTP, Haliangium, members of which were identified among the main putative predators in microcosm experiment, was the most abundant genus of Myxococcota, accounting for an average relative abundance of 2.4 (±0.1) % in the global activated sludge samples, followed by the mle1-27 clade (0.5%), Nannocystis (0.4%), and Pajaroellobacter (0.3%) (Fig. 5B). Comparatively, Bdellovibrio was the most abundant genus of Bdellovibrionota, with an average relative abundance of 0.4% (Fig. 5B).

The analysis was performed based on microbiota datasets of 1 186 activated sludge samples from 269 WWTPs globally published by the Global Water Microbiome Consortium [2]. Stacked bars (A) showed the average relative sequence abundance of Myxococcota and Bdellovibrionota across the samples collected from each country. The potential impact of predatory bacteria on sludge performance was assessed by testing the correlation (B) between their abundance and the removal of BOD (biochemical oxygen demand), COD (chemical oxygen demand), ammonia nitrogen (NH₄-N), total nitrogen (TN), and total phosphorus (TP). The color and size of the circles indicate Spearman’s rank correlation coefficients, and circles were displayed only for significant correlation (p < 0.05, n = 529, 281, 423, 268, and 364, for BOD, COD, NH₄-N, TN, and TP removal, respectively, corrected with the Benjamini–Hochberg method for multiple testing). Only the genus-level taxa with average relative sequence abundance >0.1% across all the 1 186 samples are shown, with the putative predators identified by SIP in the microcosm experiment marked with a yellow background. Genera of Myxococcota and Bdellovibrionota are indicated with purple and yellow bars, respectively. Examples showed how the sludge properties potentially influence Myxococcota, including that conductivity was negatively correlated with the relative sequence abundance of Haliangium (C) and mle1-27 (D), whereas aeration tank hydraulic retention time showed positive correlation with Haliangium (E).

To discern the potential impacts of predatory bacteria on activated sludge process performance, Spearman’s rank correlation analysis was performed to identify significant associations (FDR-adjusted p < 0.05) between removal rates of organic carbon (BOD and COD) and nutrient (NH₄-N, TN and TP) pollutants and relative abundance of Myxococcota and Bdellovibrionota genera in global WWTPs by re-analyzing performance data provided by the Global Water Microbiome Consortium [2]. Haliangium spp. were found to significantly correlate positively with COD removal (r_s = 0.16), but negatively correlate with total phosphorus (TP) removal (r_s = −0.12; Fig. 5B). The second most abundant myxobacterial mle1-27 clade displayed significantly positive correlation with all the investigated removal parameters (r_s: 0.13 to 0.31). Consistently, the mle1-27 clade correlated positively with Tetrasphaera (r_s = 0.18; Supplementary Fig. S6), a polyphosphate-accumulating organism, and a variety of denitrifying bacteria, including Dechloromonas, Zoogloea, Thauera, Comamonas (r_s: 0.12–0.17), and especially Candidatus Accumulibacter (r_s = 0.39; Supplementary Fig. S6). On the contrary, the removal of BOD, COD, and TP was negatively related with the OM27 clade (r_s: −0.13 to −0.12) and Anaeromyxobacter (r_s: −0.17 to −0.33) belonging to Bdellovibrionota and Myxococcota, respectively (Fig. 5B). In addition, no significant correlation was observed between Bdellovibrio spp. and the investigated pollutant removal rates.

Further correlation analysis revealed the differential effects imposed by activated sludge process parameters on members of Myxococcota and Bdellovibrionota. The Myxococcota and Bdellovibrionota lineages could mostly benefit from longer hydraulic retention time (Fig. 5E and Supplementary Fig. S7), whereas higher pH restrained most of these bacteria (Supplementary Fig. S7). Conductivity was negatively correlated with Haliangium (Fig. 5C), the mle1-27 clade (Fig. 5D), Pajaroellobacter, and the 0319-6G20 clade (r_s < −0.1, all adjusted p < 0.001; Supplementary Fig. S4), but was positively correlated with Nannocystis (r_s = 0.17, adjusted p < 0.01; Supplementary Fig. S7).