Phylogenomics of expanding uncultured environmental Tenericutes provides insights into their pathogenicity and evolutionary relationship with Bacilli.

BACKGROUND: The metabolic capacity, stress response and evolution of uncultured environmental Tenericutes have remained elusive, since previous studies have been largely focused on pathogenic species. In this study, we expanded analyses on Tenericutes lineages that inhabit various environments using a collection of 840 genomes.
RESULTS: Several environmental lineages were discovered inhabiting the human gut, ground water, bioreactors and hypersaline lake and spanning the Haloplasmatales and Mycoplasmatales orders. A phylogenomics analysis of Bacilli and Tenericutes genomes revealed that some uncultured Tenericutes are affiliated with novel clades in Bacilli, such as RF39, RFN20 and ML615. Erysipelotrichales and two major gut lineages, RF39 and RFN20, were found to be neighboring clades of Mycoplasmatales. We detected habitat-specific functional patterns between the pathogenic, gut and the environmental Tenericutes, where genes involved in carbohydrate storage, carbon fixation, mutation repair, environmental response and amino acid cleavage are overrepresented in the genomes of environmental lineages, perhaps as a result of environmental adaptation. We hypothesize that the two major gut lineages, namely RF39 and RFN20, are probably acetate and hydrogen producers. Furthermore, deteriorating capacity of bactoprenol synthesis for cell wall peptidoglycan precursors secretion is a potential adaptive strategy employed by these lineages in response to the gut environment.
CONCLUSIONS: This study uncovers the characteristic functions of environmental Tenericutes and their relationships with Bacilli, which sheds new light onto the pathogenicity and evolutionary processes of Mycoplasmatales.

Background

The phylum Tenericutes is composed of bacteria lacking a peptidoglycan cell wall. The most well-studied clade belonging to this phylum is Mollicutes, which contains medically relevant genera, including Mycoplasma, Ureaplasma and Acholeplasma. Almost all reported mollicutes are commensals or obligate parasites of humans, domestic animals, plants and insects [1]. Most studies so far have focused on pathogenic strains in the Mycoplasmatales order (which encompasses the genera such as Mycoplasma, Ureaplasma, Entomoplasma and Spiroplasma), resulting in their overrepresentation in current genome databases. However, Tenericutes can also be found across a wide and diverse range of environments. Recently, free-living Izemoplasma (the new name proposed by the Genome Taxonomy Database (GTDB)) and Haloplasma were reported in a deep-sea cold seep and brine pool, respectively [2, 3]. Based on their genomic features, the cell wall-lacking Izemoplasma were predicted to be hydrogen producers and DNA degraders. The Haloplasma contractile genome encodes actin and tubulin homologues, which might be required for its specific motility in deep-sea hypersaline lake [4]. These marine environmental Tenericutes exhibit metabolic versatility and adaptive flexibility. This points out the unwanted limitation that we must take into account at present when working on isolates of marine Tenericutes representatives. The paucity of marine isolates currently available has limited further mechanistic insights. Using culture-independent high-throughput sequencing techniques, Tenericutes have been detected in the gut and gonad microbiomes of fish, sea star, oysters and mussel [5-7]. As seafood consumption rises [8], there are greater concerns about food safety and control. Aside from Salmonella and Vibrio pathogens transmitted from aquaculture products [9], there are also other unknown pathogenic Mycoplasma isolates from marine animals, such as those causing ‘seal finger’ [10]. These pathogens from the ocean may be natural or human pollutants. Millions of tons of untreated sewage and sludge are dumped into the ocean yearly. Within these wastes, highly abundant Tenericutes have recently been discovered [11]. But, the spread and diversity of the Tenericutes species in oceans remain unclear.

Environmental Tenericutes might be pathogens and/or mutualistic symbionts in the gut of their host species. For example, mycoplasmas and hepatoplasmas affiliated with Mycoplasmatales play a role in degrading recalcitrant carbon sources in the stomach and pancreas of isopods [12, 13]. Spiroplasma symbionts discovered in sea cucumber guts possibly protect the host intestine from invading viruses [14]. Tenericutes were also found in the intestinal tract of healthy fish and 305 insect specimens [15, 16]. Recently, over 100 uncultured Tenericutes displaying high phylogenetic diversity were discovered in human gut metagenomes [17], irrespective of age and health status. It remains to be determined whether these novel lineages found in the human gut are linked to the maintenance of gut homeostasis and microbiome function. As a consequence of the host cell-associated lifestyle, the Tenericutes bacteria show extreme reduction in their genomes as well as reduced metabolic capacities, eliminating genes related to regulatory elements, biosynthesis of amino acids and intermediate metabolic compounds that must be imported from the host cytoplasm or tissue [18]. Beyond genome reduction, evolution of pathogenic Mycoplasmatales species has also been accompanied by acquisition of new core metabolic and virulence factors through horizontal gene transfer [19-21]. A well-studied virulence factor is hydrogen peroxide produced during the metabolism of glycerol [22]. Other virulence factors include secreted toxins, surface polysaccharides and sialic acid catabolism [23], although the mechanisms of the infection pathogenesis are largely unclear. These factors are probably obtained in the process of adaption to the hosts of Tenericutes through genomic modification. Therefore, a comparison of the genetic profiles between environmental lineages and pathogens is needed to obtain insights into the adaptation of beneficial symbionts and the emergence of new diseases.

Since Tenericutes were recently reclassified by GTDB into a Bacilli clade of Firmicutes [24], the discovery of environmental Tenericutes renovates the question regarding the boundary between Tenericutes and other clades of Bacilli. RF39 and RFN20 are two novel Tenericutes lineages of Bacilli, reported in the gut of humans and domestic animals [25, 26]. Environmental lineages of Bacilli and Tenericutes are expected to represent close relatives but their genetic relationship has not been studied. This is important to address, as uncultured environmental Tenericutes and Bacilli may potentially emerge as pathogens. In this study, we compiled the genomes of 840 Tenericutes and determined their phylogenomic relationships with Bacilli. By analyzing the functional capacity encoded in these genomes, we deciphered the major differences in metabolic spectra and adaptive strategies between the major lineages of Tenericutes, including the two dominant gut lineages RF39 and RFN20.

Results

Phylogenetic tree of 16S rRNA genes and phylogenomics of Tenericutes

We retrieved all available Tenericutes genomes from the NCBI database (April, 2019). A total of 840 genomes with ≥50% completeness and ≤ 10% contamination by foreign DNA were selected (Additional file 1). From these, 685 16S rRNA genes were extracted and clustered together when displaying at > 99% identity, resulting in 227 representative sequences. Approximately 70% of the non-redundant sequences were derived from the order Mycoplasmatales (highly represented by the hominis group), which was largely composed of commensals and pathogens isolated from plants, humans and animals. Together with 33 reference sequences from marine samples, a total of 260 16S rRNA genes were used to build a maximum-likelihood (ML) tree. Using Bacillus subtilis as an outgroup, Tenericutes 16S rRNA sequences were divided into several clades (Fig. 1a). Acholeplasma and Phytoplasma were grouped into one clade, while Izemoplasma and Haloplasma were closer to the basal group. Tenericutes species were detected across a range of environments, including mud, bioreactors, hypersaline lake sediment, and ground water. The non-human hosts of Tenericutes included marine animals, domestic animals and fungi. Sequences isolated from fungi and Hemoplasma were associated with longer branches, indicating the occurrence of a niche-specific evolution. Hepatoplasma identified as a novel genus in Mycoplasmatales is also exclusively present in the gut microbiome of amphipods and isopods [12, 27]. Spiroplasma detected in a sea cucumber gut has been described as a mutualistic endosymbiont [14], rather than a pathogen. These isolates from environmental hosts were distantly related to others in the tree, indicating a high diversity of Mycoplasmatales across a wide range of hosts and their essential role in adaptation and health of marine invertebrates. Analyses of 135 16S rRNA amplicon datasets and 141 Tara Ocean metagenomes [28] from marine waters revealed the presence of mycoplasmas from the hominis group and other sequences from the basal groups of the tree in more than 21.7% of the samples. Four of the five representative 16S rRNA sequences from the hominis group were similar (95.9–99.3%) to that of halophilic Mycoplasma todarodis isolated from squids collected near an Atlantic island [29]. The finding of the Tenericutes isolated from humans and other animal hosts in the marine samples indicates that they may be spreading possibly through sewage. The relative abundance of the 12 representative 16S rRNA genes from the marine waters was low (< 0.1%) in the microbial communities of the oceans. However, considering the tremendous body of marine water, the oceans harbor a massive Tenericutes population composed of undetected novel lineages. We detected two major clades of human gut lineages (hereafter referred to as HG1 and HG2) that were placed between Mycoplasmatales and Acholeplasmatales (Fig. 1a). These two lineages have been revealed recently as encompassing many previously unknown species in the human gut [17]. However, their contribution to human health and the core gut microbiome stability remains unclear.

Fig. 1

Phylogenetic trees of Tenericutes. The maximum-likelihood phylogenetic trees were constructed by concatenated conserved proteins (a) and 16S rRNA genes (b). The bootstrap values (> 50) are denoted by the dots on the branches. The colors of the inner layer indicate the positions of the different environmental lineages and groups of Tenericutes in the trees. Sources of the environmental lineages are shown as shapes in different colors in the outer layer

A phylogenomics analysis of Tenericutes was performed using concatenated conserved proteins from 840 Tenericutes genomes and three Firmicutes genomes. Interestingly, the topology of the phylogenomic tree coincides with that of the phylogenetic tree based on 16S rRNA genes. However, 67.6% of the genomes were derived from Mycoplasmatales, indicating a strong bias of Tenericutes genomes towards commensals, pathogens and disease-inducing isolates. The human gut lineages HG1 (n = 87) and HG2 (n = 21) were found to be neighboring clades of Mycoplasmatales as well (Fig. 1b). The genetic distance between the genomes of the gut lineages was much higher than that between the species in Mycoplasmatales, except for hemoplasmas found in infected blood and those hosted by fungi. Acholeplasma and Phytoplasma were within a clade composed of uncultured environmental Tenericutes lineages from ground waters, hypersaline sediments and mud, suggesting an environmental origin for the two genera.

By calculating the relative evolutionary divergence (RED) value of the genomes of several Tenericutes lineages [24], the average RED values for HG1 and HG2 were 0.94 ± 0.03 and 0.91 ± 0.07, respectively. Considering an expected RED value of 0.92 at the genus level, these two lineages can be considered new genera in Tenericutes. The RED value for the sequences from hypersaline lake sediments was 0.70, which supports the presence of a new order or family in Tenericutes.

Phylogenomic position of Tenericutes in bacilli

Tenericutes were recently integrated into the Bacilli clade within the Firmicutes phylum in GTDB [24]. To examine the phylogenetic positions of the new Tenericutes lineages and Bacilli, we used representative genomes of the orders within Bacilli collected by GTDB and those in Tenericutes available on NCBI. The topology of the phylogenomic relationships was supported by two ML methods. In the phylogenomic tree, four Bacilli orders, namely Staphylococcales, Exiguobacterales, Bacillales, and Lactobacillales, were clearly split from those of Tenericutes. Newly described orders RF39, RFN20 and ML615 in Bacilli, as defined by GTDB, clustered with HG1, HG2, and uncultured Tenericutes from bioreactors, respectively. This suggests that most of uncultured environmental Tenericutes submitted to the NCBI / INSDC database are probably also novel Bacilli orders, and that the genomic boundary between Tenericutes and Bacilli is thus uncertain. RF39, RFN20 and ML615 were also affiliated with Tenericutes if the boundary of Tenericutes on the tree was set at Haloplasmatales. Although RF39 and RFN20 are part of the HG1 and HG2 lineages, they have also been detected in domestic animals [30]. Interestingly, the Erysipelotrichales order was phylogenetically placed between the two human gut lineages (Fig. 2). Since all Erysipelotrichales species described in the literature so far possess a cell wall [31], their phylogenomic affinity to cell wall-lacking Tenericutes is unexpected.

Fig. 2

Phylogenetic positions of Tenericutes families in Bacilli. Representative genomes from orders of Bacilli were used to construct the phylogenomics tree using concatenated conserved proteins by IQ-TREE and RAxML. The bootstrap values were shown as triangles (50–90) and dots (> 90) with a red color for the results of RAxML and deep blue for those of IQ-TREE, respectively. The red clades represent the orders of Tenericutes. The Bacilli genomes for Erysipelotrichales and the other orders in purple were selected from GTDB. RFN20, RF39, ML615 were environmental clades named in GTDB and were phylogenetically placed within the NCBI clades consisting of human gut lineages 1, 2 and bioreactor group, respectively

We investigated the genome structure of Tenericutes and Erysipelotrichales species by calculating genome completeness, size and GC content (Additional file 3: Fig. S1). Most of the high-quality genomes (> 90% completeness and < 5% contamination) were assigned to Mycoplasmatales and Acholeplasmatales. In contrast to the rather stable genomes of the commensals and pathogenic species, the genome sizes of the uncultured Tenericutes species differed from each other and almost all were smaller than 2 Mb. Haloplasmatales genomes were the largest on average. Most of the Tenericutes genomes have a low GC content (< 30%), whereas the average GC content of those from a hypersaline lake was about 50%, consistent with a selection pressure exerted by ionic strength on the DNA double helix [32, 33]. Notably, GC content calculated on 1 kb intervals in Tenericutes genomes from ground water and HG1 (specifically RF39) varied from 20 to 70%, suggesting great plasticity and frequent gene transfers.

Genomic and functional divergence among environmental Tenericutes, commensals and pathogens

Erysipelotrichales and Tenericutes genomes were functionally annotated to characterize their metabolic pathways and stress responses that might determine the versatility and niche-specific evolution of different orders and lineages in Tenericutes. The annotation results against the Kyoto Encyclopedia of Genes and Genomes (KEGG) [34] and the clusters of orthologous groups (COGs) databases were used to calculate the percentages of the genes in the genomes (Additional file 2). Based on the frequency of all the COGs, Erysipelotrichales and Tenericutes were split into two major agglomerative hierarchical clustering (AHC) clusters. Mycoplasmatales and Phytoplasma formed AHC cluster 1, while the remaining formed cluster 2.

Using Mann-Whitney test, 203 KEGG genes and 420 COGs showed a significant difference (p < 0.01) in frequency between the two AHC clusters (Additional file 2). We selected 62 of the genes to represent those for 16 functional categories that were distinct in environmental adaptation and carbon metabolism between the two clusters (Additional file 3: Table S1 and Fig. 3). Sugars such as xylose, galactose and fructose might be fermented to L-lactate, formate and acetate by Tenericutes. The sugar sources and fermentation products differed between the groups (Fig. 3). Phosphotransferase (PTS) systems responsible for sugar cross-membrane transport were encoded by most of the genomes of Spiroplasma, Entomoplasma (including Mesoplasma) [35], Haloplasmatales, Erysipelotrichales, mycoides, and pneumoniae groups. Although most of the environmental Tenericutes genomes did not maintain PTS systems, sugar uptake might be carried out by ABC transporters. Almost all of the Tenericutes groups in the AHC cluster 2 (containing all the environmental lineages) were found to encode genes involved in starch synthesis (glgABP) and carbon storage, except for HG1. These Tenericutes groups also encoded the pullulanase gene PulA involved in starch degradation. Autotrophic pathways were present almost exclusively in environmental Tenericutes genomes. CO2 is fixed by two autotrophic steps mediated by the citrate lyase genes that function in reductive citric acid cycle (rTCA) and the 2-oxoglutarate/2-oxoacid ferredoxin oxidoreductase genes (korABCD) that encode enzymes for reductive acetyl-CoA pathway. The resulting pyruvate might be further stored as glucose and glycan via reversible Embden–Meyerhof–Parnas (EMP) pathway. Pyruvate orthophosphate dikinase (PPDK) is the key enzyme that controls the interconversion of phosphoenolpyruvate and pyruvate in prokaryotes [36]. Among all the environmental lineages and Erysipelotrichales, ppdK gene was frequently identified (73.8–100%) except for Haloplasmatales and Acholeplasmatales.

Fig. 3

Distribution of genes and pathways in the Tenericutes lineages. Tenericutes lineages were grouped using an agglomerative hierarchical clustering on the basis of the distribution of COGs within each group. The color and size of each dot represent the percentage of genomes within each lineage that carries the gene. The functions of these genes are shown in Additional file 3: Table S1

Aromatic biosynthesis pathway was lost in Mycoplasmatales, indicating their complete dependence on hosts for aromatic amino acids. Acquisition of amino acids by some environmental Tenericutes was likely conducted by peptidases (pepD2) and cross-membrane oligopeptide transporters. Glycine was also probably an important carbon and nitrogen source for the environmental Tenericutes, as a high percentage of their genomes (76.3–100%) contained the glycine cleavage genes gcvT and gcvH.

Glycerol is a key intermediate between sugar and lipid metabolisms and is imported by a facilitation factor GlpF. Phosphorylation of glycerol by a glycerol kinase (GK) is followed by oxidation to dihydroxyacetone phosphate (DHAP) by glycerol-3-phosphate (G3P) dehydrogenase (GlpD), which is further metabolized in the glycolysis pathway [37]. More than 95% of the genomes of Mesoplasma, pneumoniae, mycoides and wastewater groups contained the glpD gene; in contrast, Phytoplasma and Ureaplasma genomes lacked a glpD gene. 62% of RFN20 genomes harbored the glpD gene, while it was only found in 2% of RF39. RF39 genomes also lacked the GK-encoding gene, which suggests that RF39 cannot utilize glycerol from diet or the gut membrane. Hydrogen peroxide (H2O2) is a by-product of G3P oxidation, and has deleterious effects on epithelial surfaces in humans and animals [22]. On the other hand, these H2O2 catabolism genes were more frequently identified in uncultured environmental Tenericutes (Fig. 3).

The DNA mismatch repair machinery components MutS and MutL were almost entirely absent from Mycoplasmatales and Phytoplasma genomes. RFN20 genomes also had a low percentage of the DNA repairing genes (33.3% for mutS and 57.1% for mutL). This lack of DNA repairing genes might have generated more mutants in small asexual microbial populations capable of adapting to new environments due to Muller’s ratchet effect [38].

In Mycoplasma species as in mitochondria, tRNA anticodon base U34 can pair with any of the four bases in codon family boxes [39]. To make this ability more efficient U34 is modified in some organisms by enzymes using a carboxylated S-adenosylmethionine. The SmtA enzyme, also known as CmoM, is a methyltransferase that adds a further methyl group to U34 modified tRNA for precise decoding of mRNA and rapid growth [40, 41]. The high frequency of smtA gene in the environmental Tenericutes genomes indicates a capacity to regulate their growth under various conditions. OmpR is a two-component regulator tightly associated with a histidine kinase/phosphatase EnvZ for regulatory response to environmental osmolarity changes [42]. Its presence in most of the environmental Tenericutes genomes (> 70.4%) suggests its involvement in regulating stress responses in these organisms. The genomes of two gut lineages RFN20 and RF39 also contained a high percentage of the ompR gene. In contrast, almost all Mycoplasmatales and Phytoplasma genomes lacked the ompR gene.

The cell division/cell wall cluster transcriptional repressor MraZ can negatively regulate cell division of Tenericutes [43]. The mraZ gene that is thus responsible for dormancy of bacteria is conserved in Erysipelotrichales and Mycoplasmatales. Further studies are needed to examine whether this gene can be targeted to control pathogenicity of the bacteria in the two orders.

The Rnf proton pump system evolved in anoxic condition and is employed by anaerobes to generate proton gradients for energy conservation [44]. In single-membrane Tenericutes, proton gradients can hardly be established by the Rnf system due to the leakage of protons directly to the environment. However, this system was well preserved in genomes from Izemoplasmatales and the wastewater group. The Rnf system in these species was likely used for pumping protons out of the cell to balance cytoplasmic pH.

Metabolic model of gut lineages RFN20 and RF39

A recent study reported the genome features of RFN20 and RF39, the two main clades comprising uncultured Tenericutes [25]. The major findings on these two lineages were their small genomes and the lack of several amino acid biosynthesis pathways. After correction for genome completeness in this study, we found that the RF39 genomes were indeed significantly smaller than those of RFN20 genomes (t-test; p = 0.0012). We selected four nearly complete genomes of RFN20 and RF39 for annotation and elaborated their metabolic potentials (Table 1). The genome sizes were between 1.5 Mb–1.9 Mb, smaller than those from Sharpea azabuensis belonging to the order Erysipelotrichales. TGA is a stop codon for RFN20 genes, unlike Mycoplasmatales genes that use TGA as a tryptophan codon [23]. Coding regions of RFN20, represented by genomes HG2.1 and HG2.2 (Table 1), could be correctly predicted by using TGA as a stop codon. This was evidenced by a 20-aa unnecessary extension of the predicted translation initiation factor IF-1 in HG2.1 and HG2.2, compared with the orthologs when TGA was used as a tryptophan codon. Similar cases were observed for the other RF39 and RFN20 genes.

Table 1

Representative genomes of RFN20 and RF39. RF39 (HG1) was represented by HG1.1 and HG1.2 from the Tenericutes downloaded from NCBI; RFN20 (HG2) was represented by HG2.1 and HG2.2. S. azabuensis was a species in Erysipetrichales

ID	HG1.1	HG1.2	HG2.1	HG2.2	Sharpea azabuensis
Accession	UQAI01000000	UQAG01000000	UPZX01000000	UQBB01000000	JNKU00000000
Genome size (bp)	1,690,546	1,911,898	1,525,481	1,699,832	2,411,783
%GC	30	29.5	30.1	30.4	37.1
No.contigs	109	71	31	16	94
%Complete	98.7	98.7	98.9	98.5	99.1
%Contaminant	0	0	0	0	0.9
No. tRNA	38	35	34	45	57
No. rRNA	0	2	1	0	10
%Coding density	92	90.8	92.5	91.6	89
No. CDSs	1548	1834	1488	1570	2424

We built a schematic metabolic map for the representative RFN20 and RF39 species on the basis of the KEGG and COG annotation results. The two lineages were predicted to be acetogens since the four genomes encoded genes for acetate production (Fig. 4). We hypothesize that sugars are imported from the environment by ABC sugar transporters, while autotrophic CO2 fixation might occur via carboxylation of acetyl-CoA to pyruvate by the pyruvate:ferredoxin oxidoreductase (PFOR). Glycerol is imported and enters glycerophospholipid metabolism, which results in cardiolipin biosynthesis instead of fermentation through the EMP pathway. In some pathogenic mycoplasmas, glycerol can be taken into central carbon metabolism [37], as mentioned above.

Fig. 4

Schematic metabolism of RFN20 and RF39. We depicted the metabolic models based on gene annotation results of four representative genomes of RFN20 and RF39 (see Table 1). Solid squares indicate presence of the genes responsible for a step or a pathway. The products depicted in the MEP/DOXP pathway are 1-deoxy-xylulose 5-P, 2-C-methyl-D-erythritol 4-P, 4-(Cytidine 5′-PP)-2-C-methyl-erythritol, 2-P-4-(cytidine 5′-PP)-2-C-methyl-erythritol, 2-C-methyl-erythritol 2,4-PP, 1-hydroxy-2-methyl-2-butenyl 4-PP, dimethylallyl-PP, isopentenyl-PP, and farnesyl-PP

RFN20 and RF39 are probably mixotrophic since CO2 can be fixed to pyruvate and stored as starch, while central carbon metabolism is also connected with amino acid metabolism. After uptake of oligopeptides by the App ABC transporter system, an endo-oligopeptidase encoded by pepF yields amino acids for protein synthesis. Glycine and serine might feed into pyruvate metabolism. The peptidoglycan biosynthesis pathway was found to be complete in all four RFN20 and RF39 genomes here considered, but two genomes, namely HG1.1 and HG2.1 (Table 1), lacked the genes encoding the enzymes for UDP-N-acetylglucosamine (UDP-NAG) synthesis. Instead, these genomes harbored all the genes required for the subsequent synthesis steps to generate extracellular peptidoglycan. murG and mraY genes, which are involved in integration of UDP-NAG and UDP-N-acetylmuramate (UDP-NAM) into the peptidoglycan unit, respectively, were identified in the four genomes. With the addition of an oligopeptide, the peptidoglycan unit is secreted into the cell surface with the assistance of bactoprenol (C55 isoprenoid alcohol) [45, 46], which is formed by condensation of eight isopentenyl-diphosphate (IPP) units and one farnesyl-diphosphate (FPP). The uppS gene responsible for the bactoprenol formation was identified in the four RFN20 and RF39 genomes [47]. In bacteria, IPP can be synthesized by several metabolic steps. All the genomes contained the genes encoding the respective enzymes involved in the intermediate steps of IPP and dimethylallyl diphosphate (DPP) synthesis through MEP/DOXP pathway, except for ispD gene in one genome (Fig. 4). However, the polyprenyl synthetase gene (ispA), which is essential in the formation of FPP, was missing in three of the genomes. Given the loss of the ispA gene, the source of FPP for bactoprenol synthesis is unclear. Overall, 86.9 and 14.3% of the RF39 and RFN20 genomes contained the mraY gene, respectively, while 68.7 and 5.2% of the RF39 and RFN20 genomes had the murG gene, respectively. Therefore, most of the RFN20 genomes collected in this study lacked the complete pathway for peptidoglycan synthesis. The two essential genes for peptidoglycan synthesis were only frequently detected in Tenericutes genomes from the bioreactor group (75.0% for both genes) and Erysipelotrichales genomes (80.0 and 60.0% for mraY and murG, respectively). Therefore, the capacity of peptidoglycan synthesis is possibly deteriorating in the gut lineages, as a potential adaptive strategy to the gut environment. Similarly, the H. contractile was reported to possess the peptidoglycan synthesis genes in its genome [4], although it also lacks a cell wall. Our further examination of the genome found that the murEF genes involved in extending the oligopeptide attached on UDP-NAM were absent. This result suggests that losing the ability to synthesize a cell wall was an event that occurs independently in Mycoplasmatales and in Haloplasma. The synthesis of aminosugars NAG and NAM probably served as a mechanism of carbon and nitrogen storage for H. contractile.

RFN20 and RF39 are probably hydrogen producers, as the four genomes of HG1 and HG2 had [FeFe]-hydrogenase encoding genes. All the genomes carried the feo and fhu genes for ferrous iron uptake. Ferrous irons are taken by ABC transporters Feo into the cells when ferrous iron concentration is high in the environment. The Fhu receptor for ferrichrome absorption is required in iron-limiting condition such as the human gut [48]. The oxygen-sensitive [FeFe]-hydrogenases contain 4Fe-4S cluster and an H-cluster consisting of several conserved catalytic motifs involved in hydrogen production. Three distinct binding motifs of H-cluster in [FeFe]-hydrogenases, TSCXP, PCX2KX2E and EXMXCXGGCX2C [49], were present in the five hydrogenases encoded by all the four genomes (Additional file 3: Fig. S2). However, three of the hydrogenases from HG1 and HG2 harbor specific sites that differ from the others in some of the active sites. We have identified several orthologs with these distinct amino acids in the conserved motifs. These [FeFe]-hydrogenases formed a novel cluster in the phylogenetic tree. HG2.1 genome harbored two copies of the [FeFe]-hydrogenase genes, which were diversified as shown by their positions in the phylogenetic tree and the differences in conserved catalytic sites (Additional file 3: Fig. S2). In the human gut, three groups of [FeFe]-hydrogenases have been detected, and were proposed to be involved in methanogenesis, acetogenesis and sulfate reduction [50]. Lignocellulose-feeding termites also produce a high concentration of hydrogen in their guts, probably for degradation of wood [51]. Therefore, the HG1 and HG2 gut lineages are probably important for maintenance of a healthy gut microbial ecosystem and degradation of recalcitrant carbon.

Discussion

In the present study, we revealed phylogenomics relationships between Tenericutes lineages from different sources, and position of Tenericutes in Bacilli. Interestingly, Tenericutes lineages did not form a monophyletic clade in the phylogenomic tree of Bacilli because of the presence of the Erysipelotrichales species among them. As more environmental lineages of Bacilli are explored in the future, the taxonomic placing and monophyly of Tenericutes will probably be further challenged. Moreover, these results were dependent on the number of genomes considered from different sources and may be influenced by the quality of genome binning.

In this study, the genomic features of RFN20 and RF39 were shown to be highly dynamic among genomes from different sources. RF39 genomes lacked most of the genes for carbohydrate storage but maintained mutSL genes involved in DNA repair (Fig. 3). Except for this, there were no major differences between the two lineages, although a previous study claimed that RF39 were prone to be autotrophic [25]. Nonetheless, the predicted lifestyle of RFN20 and RF39 may vary among human populations. For example, we found that 68.7 and 76.2% of RF39 and RFN20 genomes, respectively, harbored the uppS gene for bactoprenol synthesis. However, the lack of high-quality, isolate genomes representing these lineages hinders the evaluation of their dynamics and evolutionary processes in the human gut. In deep-sea isopod gut, we also identified two types of Tenericutes bacteria, Mycoplasma sp. Bg1 and Bg2 [13]. Mycoplasma sp. Bg1 was able to degrade sialic acids probably by attachment to the host gut surface. The co-existence of two Tenericutes lineages in human and animal intestinal tracts is still enigmatic and warrants further investigations using microscopy and transcriptomics methods.

Conclusion

Our study revealed phylogenetic diversity of the Tenericutes groups and their phylogenomic relationships with Bacilli. In the environmental groups of Tenericutes, we uncovered novel lineages in human guts and marine environments, indicating the lack of environmental representatives for studies on their adaptive strategies, symbiosis and pathogenicity. Our finding of the gut lineages and their metabolic characteristics casts lights into unknown diversified mutualistic Tenericutes in gut microbiome.

Methods

Genome collection and quality check

A total of 857 Tenericutes genomes were downloaded from the NCBI database. These genomes included complete genomes of pure culture and metagenome-assembled genomes (MAGs). Three MAGs of deep-sea symbiotic Tenericutes were collected from the previous studies [13, 14]. Completeness and contamination of the genomes were evaluated by CheckM (v1.0.5) [52]. Those with > 10% contaminants and < 50% complete were removed. To explore variations of GC content in these genomes, GC content within 1-kb genome intervals were calculated. 16S rRNA genes were identified from these genomes using rRNA_HMM with default settings [53], and only those longer than 300 bp were extracted. If there was more than one 16S rRNA gene in a genome, the longest one was selected. The sequences were grouped with an identity cutoff of 99% using CD-HIT [54] and only the longest was retained as the representative. From each order of Bacilli, five genomes (see Additional file 1) were obtained from the GTDB [24]. They were selected from different families if possible.

Genome annotation and comparison

The protein coding sequences in the genomes were predicted by Prodigal (v2.6.2) [55] (proteins from Tenericutes in particular were predicted with parameter –g 4). The proteins were then searched against the eggNOG database by eggNOG-mapper (v2) [56] (with parameters --seed_orthorlog_evalue 1e-10), KEGG [34] and COG databases by Blastp with E-value cutoff of 1e-05 and similarity threshold of 40%. The functions of essential COGs belonging to Tenericutes were referred to those for a synthetic bacterium JCVI-Syn3.0 with a minimal genome [57]. RED was calculated by PhyloRank (v0.0.27) using bac120 tree in GTDB as a reference [24].

The collected Tenericutes genomes were grouped by taxonomy and source (Additional file 1). The percentage of the KEGG genes and COGs in the genomes of each group was calculated. This was also accomplished for Erysipelotrichales genomes. To filter low-frequency genes, at least one of the groups had a target gene in > 30% of the genomes. The percentages of the genes used for Bray-Curtis dissimilarity estimates were calculated using the COG frequency table. AHC [58] analysis was conducted using the pairwise dissimilarities between groups. A Mann-Whitney test was performed using the percentages of COGs and KEGG genes between the AHC clusters. The KEGG genes with p value < 0.01 were clustered into functional modules on the KEGG website (www.kegg.jp).

Phylogenetic and phylogenomic analyses

The analyses on the datasets of 16S rRNA gene amplicons from marine samples were described in our previous study [59]. The representative reads of Tenericutes OTUs were recruited for this study. Raw metagenomic data from Tara Ocean project were checked by FastQC (version 0.11.4). Reads with low quality bases (PHRED quality score < 20 over 70% of the reads) were removed using the NGS QC Toolkit [60]. The quality-filtered reads were merged using PEAR (v0.9.5) [61] and those 16S rRNA fragments > 140 bp were identified and extracted with rRNA_HMM [53]. After taxonomic classification of the fragments using the Ribosomal Database Project (RDP) classifier version 2.2 against the SILVA 128 database [62, 63], those belonging to Tenericutes were collected for the following phylogenetic study.

The 16S rRNA genes from the genomes, the amplicons and the Tara project were first clustered by MUSCLE (v3.8) [64] and then trimmed by trimAl v1.4 (automated1) [65]. An ML phylogenetic tree of 16S rRNA genes was built by IQ-TREE (v1.6.10) [66, 67] (with parameters -m GTR + F + R10 -alrt 1000 -bb 1000). Conserved proteins of the Tenericutes genomes were identified by AMPHORA2 [68]. A total of 31 conserved proteins were used to construct the phylogenomic tree for Tenericutes. The conserved proteins were aligned with MUSCLE (v3.8) [64], concatenated and then trimmed with trimAl (v1.4) (automated1) [65]. The conserved proteins from Syntrophomonas wolfei (NC_008346), Thermacetogenium phaeum (NC_018870) and Desulfallas geothermicus (NZ_FOYM01000001) were combined with the dataset of Tenericutes as an outgroup. The phylogenomics tree for Tenericutes was built by IQ-TREE (v1.6.10) [66, 67] (with parameters -m LG + F + R10 -alrt 1000 -bb 1000). The phylogenomic tree for Bacilli and Tenericutes was constructed first with IQ-TREE (v1.6.10) using the same settings as that for the phylogenomics tree of Tenericutes and then with RAxML 8.1.21 using PROTGAMMA+BLOSUM62 model with 100 bootstrap replicates.

Prediction of metabolic models of RFN20 and RF39

Four Tenericutes genomes were selected as representatives of the RFN20 and RF39 lineages on the basis of their high genome completeness. The protein-coding sequences were predicted by Prodigal (v2.6.2) [55] with codon Tables 1 and 4, respectively. The proteins were then searched against COG database [58] by Blastp [69] with an E-value cutoff of 1e-05. KEGG annotation was conducted using the online BlastKOALA tool [34]. The two sets of predicted proteins using different codon tables were compared to search for terminal extension of over 20 amino acids in the proteins due to usage of TGA codon for tryptophan. Conserved functions were then searched in the COG annotation of these proteins. Length comparison with known conserved proteins was applied to determine the usage of TGA codon for a stop codon or a tryptophan codon.

Statistical analysis

T-test and Mann-Whitney tests were performed with SPSS19 and inhouse java codes.

Additional file 1.

Additional file 2.

Additional file 3: Figure S1. GC content, completeness and size of Tenericutes. Figure S2. Phylogenetic tree of [FeFe]-hydrogenases and their conserved mitifs. Table S1. Genes and their functions in Fig. 3.