Ecological distribution and potential roles of Woesearchaeota in anaerobic biogeochemical cycling unveiled by genomic analysis.

Woesearchaeota as a newly established member of the superphylum DPANN (Diapherotrites, Parvarchaeota, Aenigmarchaeota, Nanoarchaeota and Nanohaloarchaea) are surprisingly abundant and diverse in a wide variety of environments, including deep oil reservoir, sulfuric springs and anoxic aquifers, indicating a high diversity of their roles in global biogeochemical cycles. However, ecological functions of them remain elusive. To fill up this gap, we analyzed and compared the global distribution patterns of Woesearchaeota using the genomes available publicly. As a result, both ecological distribution patterns and metabolic predictions support a key role of woesearchaeotal lineages in cycling of carbon, nitrogen, and sulfur. Multivariate regression analysis reveals that Woesearchaeota might function in consortium with methanogens in the cycling of carbon in anaerobic environments, particularly in soils or sediments. Moreover, comparative genomic analysis and ecological distribution suggest the potential roles of Woesearchaeota in the processes of denitrification, nitrogen fixation, and dissimilatory nitrite reduction, especially in the wastewater treatment systems; and also uncovered the potential capability of sulfate reduction, sulfide oxidation and thiosulfate oxidation in sulfuric or sulfidic-rich environments. Our findings add more information into the ecological roles of archaea in the anoxic environment.

Introduction

Archaea are ubiquitous microbial members composed of more than 28 phyla and occupy one of the three domains in the tree of life [1], [2], [3]. Archaeal lineages have been detected in a wide variety of environments, ranging from extreme hydrothermal vents, hypersaline environments and polar permafrost soils to freshwater, wastewater treatment plants and coastal wetland [4], [5], [6], [7], [8]. Moreover, they are abundant but widely distributed in soils, marine or freshwater sediments and water columns, wetlands and estuaries [9], [10], [11], [12]. The rapid development of culture-independent sequencing technologies and bioinformatics tools have revolutionized our understanding of archaeal diversity and metabolic potentials in different ecosystems, which allows the discovery of the new archaeal lineages and their ecological roles [13]. Therefore, archaeal diversity, evolution and ecology have been of great interest in the field of microbial ecology.

The widely-distributed archaeal lineages play crucial roles in the global biogeochemical cycles [14]. For example, in global carbon cycling, many archaea have been implicated to participate in the methane metabolism that affects both the nutrient cycling and the global climate change [15]. In deep oceans, archaea with widespread occurrence are capable of scavenging of diverse organic substrates, such as carbohydrates, fatty acids and lipids as carbon and energy sources [16]. Asgard archaea can oxidize hydrocarbons, methane and butane in particular, driving anaerobic hydrocarbon cycling in the oceanic subsurface [17]. In nitrogen cycling, ammonia-oxidizing archaea (AOA) have been reported to play a crucial role in the oxidation of ammonia in a wide variety of environments [18], [19], [20], such as wastewater treatment systems [21], [22], sediments [23], stone cultural heritage [24], [25], [26], [27]. Moreover, genes involved in nitrogen metabolism have also been detected to be abundant in temperate marine sediments [28], hyperthermophilic environments [29], brine or salted pools [30], etc. In sulfur cycling, varieties of archaeal lineages are capable of using elemental sulfur and reduced inorganic sulfur species as energy source or electron donor in many sulfuric-rich environments [31], [32], [33], such as volcanic hot springs [34], peatlands [35] and hydrothermal vents [36]. Importantly, elemental sulfur can be used as not only the electron donor by aerobic archaea (e.g., Acidianus and Sulfolobus) but also the electron acceptor by anaerobic chemolithotrophic archaea (e.g., Acidianus and Pyrodictium) [37]. Therefore, archaea become one of essential participants in the global biogeochemical cycling of major elements. Exploration of novel metabolic potential of archaea will help understand their ecological roles in the global biogeochemical cycles.

Woesearchaeota as one of the members of the archaeal superphylum DPANN are widely distributed and abundant in soils, freshwater, sediments, wetlands, estuaries and hydrothermal vents [11], [38], [39], [40]. High diversity of woesearchaeotal lineages yields more than 26 subgroups at the class level based on the phylogenetic analysis [39]. High diversity and wide ecological distribution imply the important roles of woesearchaeotal lineages in the global biogeochemical cycles [41]. However, current studies have only revealed the roles in anaerobic carbon cycling [38], [39], and thus the major functions of Woesearchaeota remain still unknown. Identification of novel metabolic potential of woesearchaeotal lineages has great importance to unveil their ecological roles in the global biogeochemical cycles of major elements.

At present, a great number of woesearchaeotal 16S rRNA gene sequences (ca. 5000) in public databases are originated from a wide range of environments. This dataset provides large information available currently to explore the global ecological patterns and to imply more on the relationship between the distribution patterns and the ecological functions/roles of Woesearchaeota in the global biogeochemical cycles. Here, we retrieved the public databases of woesearchaeotal lineages and meta-genomes to (i) characterize the global distribution of woesearchaeotal communities and unveil the driving environmental factors that affect the woesearchaeotal distribution among diverse habitats, (ii) identify the metabolic potential on the basis of comparative genomic analysis of the woesearchaeotal genomes, and (iii) link the ecological distribution patterns to their potential ecological roles in the cycling of carbon, nitrogen, and sulfur. We expect to expand our understanding of the ecological roles of archaea in the global biogeochemical cycles.

Materials and methods

Construction of cultivation-independent libraries

We recovered both the published literature and the GenBank database for the construction of woesearchaeotal 16S rRNA clone libraries. The detailed steps involved are available elsewhere [39]. Specifically, the Esearch Utility was employed to capture archaeal 16S rRNA gene sequences from the database as per the following string ‘16S AND 600:2000 [Sequence Length] AND archaea [Organism] AND rrna [Feature key] AND isolation_source [All fields] NOT genome OR chromosome OR plasmid’. To make sure the availability of further analysis, sequences with isolation source tags missing were removed. To guarantee only studies or libraries that include woesearchaeotal sequences were retained, the remaining 122,559 archaeal sequences from 217 libraries were BLASTed against the reference database SILVA SSU 128 [42], [43]. Libraries or studies with the number of 16S rRNA gene sequences <10 were discarded. Finally, we accomplished the screening with 133 studies or libraries, which contain 15,012 archaeal sequences, including 3584 woesearchaeotal sequences. Such studies or libraries were grouped into seven distinct habitats as follows: 32 libraries from freshwater sediment (Fsed), 30 from freshwater (Fwc), 22 from soil (S), 19 from marine sediment (Msed), 12 from hydrothermal vent (Hdv), 11 from hypersaline environment (Hsal), and 7 from marine water column (Mwc) (Table S1). Furthermore, a semi-quantitative environmental matrix was constructed according to the range of environmental gradients among these habitats: temperature (cold, temperate, and hot), life style (water, sediment, and soil), salinity (non-saline, saline, and hypersaline), oxygen concentrations (oxic and anoxic), and trophic state (oligotrophic, mesotrophic, and hypertrophic) (Table S1).

Taxon-based analysis

The total 15,012 archaeal sequences formed 11,323 operational taxonomic units (OTUs) at an identity of 97%, including 2180 OTUs that belonged to Woesearchaeota (Table S1). Next, these OTU sequences were compared with the reference database SILVA SSU 128 and assigned to the corresponding taxon once one OTU matches with a reference or identified sequence with an identity higher than 97%. If an archaeal OTU does not match with any reference sequence, this OTU will be named ‘unknown archaea’. After the taxonomy assignment, each named OTU was dated back to the specific libraries that were associated with a given habitat or environmental matrix for further analyses, such as relative abundance and occurrence frequency.

Multivariate regression analysis

A multivariate regression tree (MRT) was determined and generated by the R package mvpart to illustrate the relationship between the environmental matrix (Table S1) and the table of lineage relative abundances [44]. Meanwhile, we introduced the indicator value (IndVal) index, which combines relative occurrence frequency and relative abundance to identify archaeal lineages, like the concept of ‘indicator species’, among the habitats. The indicator species represents the predominance of a species in both relative abundance and relative occurrence frequency in a given habitat, indicating the important ecological niche of a species in such habitat.

Comparative genomic analysis

Genomes of Woesearchaeota were retrieved from the NCBI Genome Database and further confirmed with the authors/owners. The quality of the retrieved 19 woesearchaeotal genomes was estimated by CheckM and the parameters, such as completeness (%), sequencing depth (x), relative GC content (%), and number of protein coding genes were listed in Table S2. To infer the metabolic functions of specific subgroups of Woesearchaeota, only genomes that have a 16S rRNA gene sequence were used for analysis. Genomic BLAST analysis for all the genomes was performed to provide the reciprocal best BLAST hits (rBBHs). The MCL76 algorithm was then introduced to cluster rBBHs into protein clusters. Finally, the proteins were ready for further annotation and identification.

Annotation and identification of genes

Open reading frames (ORFs) were predicted on their genomic scaffolds using the metagenome mode of Prodigal and were annotated by the NCBI Prokaryotic Genome Annotation Pipeline (Table S3) and by the comparison with KEGG PATHWAY Database (https://www.genome.jp/kegg/pathway.html). Genome-based functional gene annotation was determined by searching all predicted ORFs in a genome with the presence of the corresponding enzymes. If a metabolic process is dependent upon more than one enzyme, it was confirmed based on the presence of all the necessary enzymes. To find out the metabolic potential of the specific woesearchaeotal subgroups, we assigned the genomes to the corresponding subgroups on the basis of their 16S rRNA gene sequences if available [39] (Table S2).

Results

Ecological distribution patterns of Woesearchaeota

Woesearchaeota over the retrieved 133 archaeal libraries had a relative abundance of 23% and dominated differentially in the anoxic habitats, particularly in estuary sediments, lake and wetland sediments, and cold sulfuric springs (Fig. 1 and Table S1). For example, the abundance of Woesearchaeota throughout the 5 libraries of estuary sediments reached up to 70%, whereas the abundance over the 5 libraries of cold sulfuric springs was more than 50%. They tended to occur more frequently in paddy soils, lake and wetland sediments, and hydrothermal vents, where oxygen availability is largely limited or none existence.

Fig. 1

Relative abundance of woesearchaeotal lineages (based on sequence abundance) within anoxic habitats. The number of libraries used for each kind of habitat is given above the bar. Error bars represent the standard deviation.

To link the abundance of the archaeal lineages to environmental factors, MRT was constructed (Fig. 2). The MRT analysis generated a twelve-leaf tree ordination based on oxic status for the first two branches, totally explaining 51.1% of the phylogenetic variance. Anoxia tended to combine S (soil with a depth of more than 10 cm) and anoxic water–sediment, Fsed and Fwc (sediments and water column from anoxic freshwaters or springs), as well as Hdv and Msed (hot- and cold-temperate anoxic marine sites), which agrees with previous study [4]. Woesearchaeota are the indicator lineage over the 7 leaves of the anoxic branch of the tree, covering 103 libraries. However, Woesearchaeota are detected only in 30 libraries on the oxic branch of 5 leaves, which represents only a small portion of the ecological distribution. Therefore, both distribution and abundance patterns indicate the predominance of Woesearchaeota in anoxic environments.

Fig. 2

MRT of the interaction between archaeal lineage abundance and environmental factors. The model explained 51.1% of the variance of the whole data set. Libraries clustered in the tree leaves according to their sources (Table S1). Pie charts under each leaf indicate how the relative abundance of archaeal lineages contributed to the structure of the leaves. The IndVal index for each leaf (data not shown) revealed that most core lineages are still indicator lineages (p < 0.01), which are labeled on each pie. Abbreviation of taxonomy: Woese, Woesearchaeota; Thermpl, Thermoplasmata; SCG, Soil Crenarchaeotic Group; Tergp, Terrestrial group; MG-I, Marine Group I; UncThaum, Uncultured Thaumarchaeota; Halo, Halobacteria; Bathy-8/-6, Bathyarchaeota subgroup-8/-6; Metcoc, Methanococci; Metmic, Methanomicrobia; Metbac, Methanobacteria; SAGMCG, South African Gold Mine Group 1.

Notably, the MRT analysis showed that Woesearchaeota are the only lineage to be indicator for all anoxic leaves and also the co-indicator lineage of Methanomicrobia in soil, anoxic water and sediments (Fig. 2), such as paddy soil, permafrost soil, oil and gas fields, and freshwater sediments (Fig. 1). Moreover, Woesearchaeota are a co-indicator lineage of Methanococci in temperate fresh sediments and Methanobacteria in anoxic water. Such findings might further provide ecological evidence for Woesearchaeota as a potential partner of methanogens in methanogenesis, which takes place in the anaerobic environments lacking sulfate.

Potential functions in nitrogen cycling

We identified some genes involved in the processes of denitrification (e.g., nirK, norD, and nosZ), nitrogen fixation (e.g., nifH), and dissimilatory nitrate reduction to ammonium (e.g., nrfA and nirB) as well as the genes relevant to the conversion of organic nitrogen (e.g., ure ABC) (Fig. 3). However, genes involved in the processes of anammox, aerobic ammonia oxidation, nitrate reduction, and nitrite oxidation are missing in Woesearchaeota. Such functional properties further agree with the distribution patterns that Woesearchaeota dominate in anoxic environments, particularly in the wastewater treatment systems (Fig. 1), where nitrogen is normally very rich and anaerobic nitrogen cycling processes are very active.

Fig. 3

Metabolic potential of Woesearchaeota in nitrogen cycling. Colorful nodes indicate the corresponding woesearchaeotal subgroups that might hold the metabolic function of the corresponding pathways. Only genes detected are listed for the corresponding pathways. Dash lines indicate the pathways that are not detected in any woesearchaeotal subgroups based on the current genomic information available, whereas the solid lines indicate the pathways that might occur. The two arrows outside the chart indicate the start of the pathways that take place under anoxic or oxic condition.

Potential roles in sulfur cycling

Comparative genomic analysis also revealed the potential roles of Woesearchaeota in sulfur cycling, particular for the reduction of sulfate (e.g., the genes sat, aprA, and dsrA) and thiosulfate oxidation (e.g., soxB) as well as sulfide oxidation (e.g., sqr) (Fig. 4). All of the genomes contained the genes sat and/or aprA that encode the enzymes that catalyze the reduction of sulfate to sulfite. Members of Woese-3, Woese-18 or the other ungrouped might complete the sulfate reduction to sulfide. Furthermore, members of Woese-3, Woese-5a or the other ungrouped might oxidize sulfide to sulfur. Only the members of Woese-18 show the potential to perform the thiosulfate oxidation. However, genes involved in the processes of thiosulfate reduction, sulfur oxidation, sulfur reduction or sulfite oxidation are not detected in the genomes of Woesearchaeota in this study. These functions in combination with the distribution patterns support the potential ecological role of Woesearchaeota in the sulfuric or sulfidic environments, particularly for hydrothermal vent and cold or hot sulfuric or sulfidic springs (Fig. 1), where sulfur is typically abundant and the element cycling is actively carried out.

Fig. 4

Metabolic potential of Woesearchaeota in sulfur cycling. Colorful nodes indicate the corresponding woesearchaeotal subgroups that might hold the metabolic function of the corresponding pathways. Only genes detected are listed for the corresponding pathways. Dash lines indicate the pathways that are not detected in any woesearchaeotal subgroups based on the current genomic information, whereas the solid lines indicate the pathways that might occur.

Discussion

Dominance of Woesearchaeota in anoxic environments

The ecological distribution patterns indicated that Woesearchaeota dominate in anoxic environments, such as soil, sediments and other anoxic habitats. Meanwhile, high biodiversity of Woesearchaeota in such habitats has also been observed in our previous work [39]. These results reflect the important ecological niche specificity of Woesearchaeota to such environments, because niche differences over space allow the formation of biogeographic patterns of species [45]. The properties of the environments would reveal the ecological functions of woesearchaeotal communities in anaerobic biogeochemical cycles of elements [38]. For example, paddy soils [12], oil or gas fields [46], [47], [48], [49], wetlands or sediments [10], [50], [51], [52] are the typical habitats of methane-producing organisms, which contribute greatly to the global biological methane emission. In wastewater treatment systems, nitrogen cycling is actively driven by nitrifiers or denitrifiers [22], [53], [54], [55]. In sulfur-rich hydrothermal vents and cold or hot springs, sulfur-cycling archaea and bacteria manipulate the reduction and oxidation of elemental sulphur [36], [56], [57], [58]. These ecological implications reveal the potential roles of woesearchaeotal lineages in the specific biotopes. However, the frequency of the habitats may not represent the real ecological distribution of Woesearcheota, given the small sample sizes or the choice of habitats depent upon scientists’ interest in the studies.

Potential partner of methanogens in anaerobic carbon cycling

According to the distribution patterns and metabolic modeling, it seems that Woesearchaeota tends to co-occur with methanogens by a mean of consortium or syntrophy [39]. The MRT analysis showed that Woesearchaeota are the co-indicator lineage of methanogens in many anoxic environments, which further supports the hypothesis that Woesearchaeota might function as a partner with methanogens biochemically in anaerobic carbon cycling. Although the direct evidence has not yet been obtained or confirmed, many indirect findings support this hypothesis. Woesearchaeota have been found to perform the fermentation lifestyle and convert organic substrates (e.g., starch) into acetic acid or hydrogen [38], which might favor the growth of acetotrophic or hydrogenotrophic methanogens [39]. The retroelement-guided protein diversification mechanisms widely-existed in DPANN super-phylum, including Woesearchaeota, have revealed the presence of a versatile tool that could help them adapt to a dynamic, host-dependent existence [59]. As methanogenesis has a great impact on the global carbon cycle and climate change, further confirmation of the participation by Woesearchaeota is extremely important. Further work should be focused on the two directions to determine the effects. On the one hand, given the fact that the identified number of woesearchaeotal genomes is still very small, it is necessary to expand the genomic information of woesearchaeotal communities to support the direct evidence by exploring the methyl-coenzyme M reductase (Mcr) complex, which is the key enzyme that catalyses the final step in methanogenesis and the first step in methanotrophy [15]. Meanwhile, the analysis using more undiscovered genomes from both archaea groups for metabolic models’ reconstruction may provide further evidence to reveal their partnership in C cycling in anaerobic environments. On the other hand, co-culture of Woesearchaeota and methanogens or other highly co-occurring dominant microbial members based on metabolic complementation could also help prove the hypothesis of the partnership or syntrophic lifestyle.

Participant in nitrogen cycling

Woesearchaeota are also found to be abundant in the libraries of wastewater treatment systems, such as anaerobic wastewater bioreactor, activated sludge, and anaerobic digester plan. Such environments are characterized by the active nitrogen cycling driven by nitrifiers or denitrifiers, which include archaea, bacteria and fungi [21], [22], [60], [61], [62]. A few genes, nirK and nosZ in particular, involved in nitrogen cycling are detected in Woesearchaeota. The genes are mainly responsible for biochemical processes of nitrogen cycle under anoxic conditions, such as denitrification from nitrite to dinitrogen (N2), N2 fixation, dissimilatory nitrite reduction to ammonium as well as the ammonification of organic nitrogen. Thus, both ecological distributions and metabolic predictions revealed the potential roles of Woesearchaeota as a participant in nitrogen cycling in anoxic environments. However, it is also possible that the denitrification process from nitrite to N2 is a typical mechanism of detoxification in Woesearchaeota [38]. If this hypothesis is true, it is difficult to explain the dominance of Woesearchaeota in such nitrogen-rich environments because they may need more energy to cope with the damage from toxic nitrite or other products generated in the nitrogen cycle [63], [64]. Therefore, further experimental tests are still required to elucidate the functions of Woesearchaeota in the nitrogen cycle. Isotopic labeling test with specific substrates would be an effective way to confirm the actual expression of the genes in woesearchaeotal genomes.

In addition, the comparative genomic analysis indicated some potential roles of Woesearchaeota in sulfur cycling, such as sulfate reduction to sulfide, sulfide oxidation to sulfur, and thiosulfate oxidation. Importantly, such functions again agree with the high abundance of woesearchaeotal lineages in sulfuric or sulfidic-rich environments, such as hydrothermal vents and cold or hot springs [58], [65], [66], [67]. Therefore, these findings would well support the potential roles of Woesearchaeota in sulfur cycling. If it is true, this information will greatly advance our understanding of archaea in the global biogeochemical cycling of sulfur. Investigation into sulfur metabolism in archaea could help explore many novel enzymes and pathways of the archaea, particularly from extreme environments. For example, two anaerobic chemolithotrophic archaeal genera, Acidianus and Pyrodictium, have been found to be able to use elemental sulfur as the electron acceptor in a short electron transport chain based on a sulfur reductase and a membrane-bound hydrogenase in hyperthermal habitats [37]. In contrast, the relevant mechanisms remain unclear in heterotrophic anaerobic archaea, such as Pyrococcus and Thermococcus, though they are also found to use element sulfur as the electron acceptor. In a low-salt, sulfide- and sulfur-rich spring, Halobacteriales may adopt a novel sulfur-metabolic mechanism different from that in their typical high-salt habitats [68]. Such metabolic mechanisms and pathways would evolve diversely in the specific environments and, thus, could be one reason to explain archaeal distribution and diversity with ecological function on Earth.

Further studies for functional confirmation are still required. On the one hand, cultivation of the samples in these environments or in situ experiments might help us to confirm the roles of Woesearchaeota. On the other hand, molecular test based on synthetic biology is an effective approach to verify the function of the specific gene. Moreover, a comparative analysis of metabolic potential of Woesearchaeota from different environments will provide more insightful details on their ecological roles and functions in the biogeochemical cycles, given the fact that more genomes are discovered from different environments.

Conclusions

In this study, we investigated the ecological distribution patterns of woesearchaeotal lineages and linked the habitat properties to the predicted functions. The agreement between ecological distribution patterns and metabolic functions reveals the potential roles of Woesearchaeota in anaerobic biogeochemical cycling of carbon, nitrogen and sulfur. These findings will advance our understanding the ecological functions of the new archaea in the global biogeochemical cycles. Since thre is no pure culture of Woesearchaeota available now, information derived for this group is based entirely on 16S rRNA genes and Metagenome Assembled Genomes of this group, and further experimental demonstration is still required to confirm their ecophysiology in the ecosystems.

CRediT authorship contribution statement

Xiaobo Liu: Conceptualization, Methodology, Software, Data curation, Writing - original draft. Ji-Dong Gu: Supervision, Writing - review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.