Putative functions of overall microbial communities
The putative functions of the microbial communities in the ecosystems were analyzed by using FAPROTAX. In each community, 8.9% ~ 60.0% sequences were assigned to a certain or several functional processes. A total of 64 functional groups were detected, with 29 to 57 groups present in each community. The average RA distribution of the known functional taxa in each dataset is shown in Fig. 1. The functional populations with an average RA lower than 1% were classified into the group of ‘Others’. RAs of several functional categories of the communities in every two datasets showed significant differences, revealed by Dunn's Kruskal-Wallis multiple comparisons (Table 1).
In the carbon cycling process, chemoheterotrophy related taxa were the most abundant in each dataset except FW. Their RAs in the AD systems were the highest (43.7%±4.2%), followed by those of the soil systems (40.2%±1.9% in SA and 36.4%±2.0% in SB). On the one hand, aerobic chemoheterotrophy related taxa were with the lowest RAs in the AD systems (5.0%±2.8%) and with the highest RAs in the soil systems (34.2%±3.7% in SA and 34.8%±2.0% in SA). On the other hand, fermentation associated populations were with the highest RAs in AD (36.6%±6.6%), and took up 4.6%±2.9% and 1.2%±0.4% in SA and SB, respectively. These results indicate that chemoheterotrophic bacteria were relatively abundant in the AD systems and soil environments, with aerobic chemoheterotrophs abundant in the soil ecosystems and anerobic chemoheterotrophs abundant in AD, respectively. This can be possibly explained by the characteristics of the habitats. AD systems are designed to convert organic wastes into a methane-rich gas by microorganisms in the absence of oxygen. They thus harbor higher abundances of chemoheterotrophs and fermenters. SA and SB contained certain organic matters [3, 10] and can support a variety of chemoheterotrophs, especially aerobic chemoheterotrophs.
In the nitrogen cycling process, nitrate reduction populations were the most abundant in the WWTPs, accounting for 17.2%, 10.9%±1.5% and 7.1%±1.2% in FW, WA and WB, respectively. Aeration in WWTPs supports the process of nitrification, with nitrate as the final form of oxidized nitrogen. This can give an explanation to the high RAs of nitrate respiration populations in WWTP systems in this study. Also, the nitrogen/nitrate respiration populations were the most abundant in WWTPs, accounting for 17.0%, 9.9%±1.6% and 6.5%±1.2% in the FW, WA and WB, respectively. This may result from the high loading rate of nitrogen in WWTPs compared with other ecosystems. The ureolysis populations were the most abundant in soil systems, accounting for 2.9%±1.5% and 2.9%±1.3% in SA and SB, respectively, higher than those of other ecosystems (p < 0.05). Urea is usually added to soils as a fertilizer, which can be a possible reason for the high proportions of ureolysis populations in SA and SB.
In the sulfur cycling, the taxa participating in sulfur compound respiration/sulfate respiration were with the highest RAs (16.9%±13.0%) in CSS. Sulfate is a typical electron acceptor in marine sediment as it is with a high concentration in seawater (~28 mM) [35]. The biogeochemical sulfur cycle of marine sediments is intensively discussed and reviewed in previous studies [35, 36].
Identification of indicative OTUs of the habitats
A total of 683, 228, 208, 715, 776, 505 and 2551 OTUs were detected to be indicative to the AD, WA, WB, FW, SA, SB and CSS datasets, respectively. The indicative OTUs accounted for 18.0%, 10.8%, 4.2%, 60.5%, 9.6%, 7.2% and 41.8% of the total OTUs detected in each of the above dataset, respectively. The cumulative abundances of the indicative populations for the above ecosystems were 41.5%, 20.0%, 14.1%, 50.5%, 31.1%, 42.6% and 78.4% on average, respectively. This indicates the indicative OTUs were relatively abundant compared with other OTUs in most of the ecosystems. The high proportion of indicative community in CSS can be partly explained by its unique geographic location in the investigated sampling sites. The indicative communities were significantly correlated to pH (r = 0.1476, p < 0.05), temperature (r = 0.3471, p = 0.001), latitude (r = 0.0975, p < 0.05) and longitude (r = 0.4295, p = 0.001), respectively, as revealed by mantel test. In our previous study, however, overall communities did not show significant correlations to pH values or longitude. These together indicate that indicative communities can reflect the abiotic state environment well, consistent with the results of previous studies [14].
Comparisons of phylum composition of indicative and overall communities
Comparisons of the phylogenetic composition of the indicative community and the overall community of each dataset were made at phylum level. In the SA dataset, Crenarchaeota, Gemmatimonadetes and Verrucomicrobia were more abundant in the indicative communities (accounting for 4.50%±3.27%, 11.24%±4.49% and 4.34%±1.13%, respectively) than in the overall communities (accounting for 1.56%±1.19%, 4.77%±2.12% and 2.65%±0.59%, respectively) (p < 0.05, p < 0.001 and p < 0.001, respectively). In the SB dataset, Verrucomicrobia was greatly more abundant in the indicative communities (22.51%±8.18%) than in the overall communities (10.62%±3.05%) (p < 0.001). As discussed in our previous study [11], the phyla Gemmatimonadetes [37] and Verrucomicrobia [37, 38] are known typically as K-strategists capable of efficient utilization of scarce resources in the environment. Crenarchaeota has been reported to be present in heavy metal contaminated soil [27, 39]. This phylum was reported to show significant positive correlations with Cd, As and Zn concentrations in SA [27]. The higher RAs of Crenarchaeota in the indicative community reflect SA’s environmental characteristics that the investigated farmlands were with heavy metal contamination.
In the CSS dataset, Chloroflexi was also slightly more abundant in the indicative community than in the overall community, accounting for 8.64%±3.68% and 7.66%±3.27%, respectively (p < 0.001). In the AD dataset, this phylum was also with higher RAs in the indicative communities (14.53%±15.37%) than in the overall communities (7.58%±7.80%) (p < 0.01). The higher RAs of Chloroflexi in the indicative communities of ADs and sea sediments suggest the importance of this phylum in these ecosystems. Chloroflexi is reported to be present in a wide range of habitats including microbial mats, the deep oceans, fresh water aquifers and lakes, as well as wastewater treatment systems [40]. Chloroflexi in ADs was suggested to have the functions of carbohydrate degradation and cellular matter degradation [41]. Other isolates of this phylum from subseafloor sediment also show metabolic pathways of carbohydrate degradation [42].
In WB dataset, Proteobacteria was more abundant in the indicative communities than in the overall communities, accounting for 69.95%±11.55% and 56.12%±9.38%, respectively (p < 0.001). Proteobacteria (with an average rrn number of 4.9) [43, 44] are generally classified as r-strategists. Taxa from this phylum have been reported to degrade diverse carbohydrates [45] and participate in denitrification [46, 47]. The higher RAs of Proteobacteria in the WB indicative community may be attributed to the high pollutant concentrations of WWTPs.
At genus level, two genera Mycobacterium and Herbaspirillum were associated with six datasets (CSS, FW, SA, SB, WA and WB), indicating their strong ability of habitat adaption. In FAPROTAX database, Mycobacterium is shown to play a role in the processes of aerobic chemoheterotrophy and chemoheterotrophy. Herbaspirillum is suggested to be linked with the processes of ureolysis and nitrogen fixation.
Comparisons of functional composition of indicative and overall communities
Functional compositions of the indicative community and the overall community of each dataset were compared, with significant differentiations detected (Table 2). For each dataset, only the functional categories with average RAs over 5% in both overall and indicative communities and increased significantly by over 5% in the indicative community were listed. In the SA and SB datasets, aerobic chemoheterotrophy related microbes were more abundant in the indicative communities. The CSS indicative community showed higher RAs in sulfate respiration, respiration of sulfur compounds and fermentation. The AD indicative community showed higher RAs in fermentation and chemoheterotrophy. This can be explained by the environmental characteristics of each kind of habitat discussed in Part 3.1. The WB indicative community showed higher RAs in aerobic chemoheterotrophy, chemoheterotrophy and predatory/exoparasitic related bacteria. Noticeably, predatory/exoparasitic populations were greatly more abundant (40.06%±8.92%) than that of the overall community (5.60%±2.19%) (p < 0.001). Predatory/exoparasitic related population in the investigated communities were from the families of Burkholderiaceae, Bacteriovoracaceae, Bdellovibrionaceae, Haliangiaceae and Phaselicystidaceae. Bacteriovoracaceae and Bdellovibrionaceae are classified as Bdellovibrio-and-like organisms (BALOs) [48]. BALOs are obligate predators of a variety of gram-negative bacteria including human pathogens [49]. Bdellovibrio bacteriovorus is a typical BALO and serves as a model organism for bacterial predation [49]. B. bacteriovorus was reported to be present in activated sludge [50]. The RAs of the genus Bdellovibrio were 1.28±2.41% in the WB indicative communities and can be up to 11.19% in some sample, while its RAs were 0.47±0.52% in the WB overall community. Predatory bacteria are suggested to be a key factor in limiting waste activated sludge production [51], which can be an explanation of the abundant predatory/exoparasitic population in the WB overall community.
Microbial co-occurrence patterns of the ecosystems
Potential microbial interactions were shown by analyzing the co-occurrence patterns of the microbial communities from WB, SB and AD datasets through network construction based on a RMT algorithm. A network of 1,065 nodes and 2,086 links was obtained for the activated sludge (AS) microbial communities from the WB dataset, with 29.6% correlations positive and 70.4% negative. A network of 219 nodes and 349 links was obtained for the AD microbial communities from the AD dataset, with 53.3% correlations positive and 46.7% negative. A network of 849 nodes and 1758 links was obtained for the microbial communities in the broadleaved forest from the SB dataset, with 63.2% correlations positive and 36.8% negative. The networks of AD, SB and WB are displayed in Figure 2. The topological characteristics of the three networks are shown in Table 3. The constructed networks exhibited general topological features including scale free, small world and modularity (Supplementary Text S1 online).
The AD network showed the highest proportion of positive correlations. This can be explained by the microbiology of anaerobic transformation of organic wastes. This process involves several different bacterial species, such as hydrolytic, acid forming, acetogenic, and methanogenic bacteria with CO2 and CH4 as the main products [52, 53]. Acetogenic bacteria and methanogenic bacteria cannot directly metabolize polymers, without the help of hydrolytic bacteria who digest polymer into simpler soluble monomers. From this perspective, microbial cooperation widely exists in ADs. Thus, a high proportion of positive correlations can be detected between the taxa. The proportion of positive correlations of AD microbial communities detected in this investigation was comparable to the values of other networks built for AD communities [29].
The WB network showed the highest proportion of negative correlations. Negative correlations may possibly indicate competition between the taxa [54]. WWTPs are designed to remove pollutants from anthropogenic activities. The microbial inhabitants are confronted with environmental pressure such influent fluctuations and exposure to toxic substances [55], which may aggravate the competition between the taxa. Similar results were reported in previous studies which indicated that salinity and lower food to microorganism ratio (F/M) may possibly lead to increased negative links of microbial components in lakes [14] and wastewater treatment systems [28], respectively.
The nodes of AD, SB and WB networks showed an average PD of 0.51, 0.50 and 0.59, respectively, lower than the average PD of all the detected OTUs. Noticeably, the positively linked OTUs showed a lower PD value than the negatively linked OTUs in each network. In the AD network, the average PD of negatively correlated OTUs was 0.59, and that of positively correlated OTUs was 0.44. In the WB network, the average PD of negatively correlated OTUs was 0.60, and that of positively correlated OTUs was 0.56. In the SB network, the average PD of negatively correlated OTUs was 0.54, and that of positively correlated OTUs was 0.48. This is consistent with the results of other studies. These studies together suggest that positively correlated taxa are more phylogenetically close in general. This is usually attributed to niche overlap of phylogenetically similar taxa [56].
Two OTUs were present in all the networks constructed, both assigned to the order Actinomycetales. This indicates that Actinomycetales may play important roles in a wide range of habitats. One of them was from the genus Janibacter, the other with unknown affiliation at genus level. Janibacter is reported to show both antibiotic resistance/pathogenicity and the ability to degrade pollutants [57].
The nodes were classified into four subcategories by within-module connectivity (zi) and among module connectivity (Pi). In each network, most OTUs were peripheral nodes, accounting for 97.3%, 97.1% and 99.0% in the AD, WB and SB networks, respectively. This reveals a general topological feature of the microbial association networks that they are scale-free. This feature generally implies the presence of many taxa with only a few links and a few highly connected (hub) taxa.
In the AD network, a connecter and five module hubs were detected. The connecter was affiliated to the genus Aminomonas. Aminomonas was reported to be a dominant genus and characterized by anaerobic fermentation of amino acids in an anaerobic membrane bioreactor (AnMBR) treating urban wastewater and food waste [58]. Three module hubs (60%) were from Firmicutes, suggesting the importance of this phylum in the network. Two of them were from the genus Clostridium sensu stricto. Clostridium sensu stricto was shown to be abundant in anaerobic reactors in previous studies [59, 60]. Also, this genus was shown to be abundant exoelectrogens in anaerobic digestion with direct interspecies electron transfer stimulation through nickel supplementation [61], which may indicate intensive connections between Clostridium sensu stricto and other microorganisms.
In the SB network, seven connecters and 18 module hubs were detected. Among the module hubs, seven were from Proteobacteria (38.9%) and five were from Acidobacteria (27.8%). Among the connectors, three were from Proteobacteria (42.9%) and three were from Acidobacteria (42.9%). These results together suggest the importance of Proteobacteria and Acidobacteria in the SB network. These Acidobacteria OTUs were assigned to the genera Gp2, Gp5 and Gp6. Their high abundance in the soil samples was also reported in our previous investigation [11]. Similarly, Acidobacteria Gp2 was identified as hubs in most of subnetworks of the microbiomes in various environments in a previous publication [30].
In the WB network, four connecters and seven module hubs were detected. Three module hubs (42.9%) and two connectors (50%) were from Proteobacteria, indicating the importance of this phylum in the WB network. Two connectors in WB were assigned to the genera of Thauera and Acidobacteria Gp4, respectively. Thauera has been reported to be frequently detected in WWTPs [62, 63] and play important roles denitrification [46, 47] and carbon degradation [63, 64]. The wide presence of Acidobacteria Gp4 in WWTPs has been shown by several investigations [62, 65]. Gp4 was reported to be one of the characteristic genera for protein hydrolysis and short chain fatty acids (SCFAs) accumulation in waste activated sludge fermentation [66]. The high abundance and versatile functions of Thauera and Acidobacteria Gp4 may explain their important topologic roles in the WB network. One module hub was from Byssovorax, which was reported to be able to degrade cellulose [67].