Globally distributed marine Gemmatimonadota have unique genomic potentials

Abstract

Background

Gemmatimonadota bacteria are widely distributed in nature, but their metabolic potential and ecological roles in marine environments are poorly understood.

Results

Here, we obtained 495 metagenome-assembled genomes (MAGs), and associated viruses, from coastal to deep-sea sediments around the world. We used this expanded genomic catalog to compare the protein composition and update the phylogeny of these bacteria. The marine Gemmatimonadota are phylogenetically different from those previously reported from terrestrial environments. Functional analyses of these genomes revealed these marine genotypes are capable of degradation of complex organic carbon, denitrification, sulfate reduction, and oxidizing sulfide and sulfite. Interestingly, there is widespread genetic potential for secondary metabolite biosynthesis across Gemmatimonadota, which may represent an unexplored source of novel natural products. Furthermore, viruses associated with Gemmatimonadota have the potential to “hijack” and manipulate host metabolism, including the assembly of the lipopolysaccharide in their hosts.

Conclusions

This expanded genomic diversity advances our understanding of these globally distributed bacteria across a variety of ecosystems and reveals genetic distinctions between those in terrestrial and marine communities.

Video Abstract

Introduction

Microorganisms are crucial components of ecosystems on Earth, playing important roles in global food webs and driving biogeochemical cycling. Cultivation-based research advanced our understanding of fundamental roles of microbes in the environment. However, our understanding of microbial ecology has dramatically changed in recent years due to the application of high-throughput sequencing technologies [1]. Single gene-based surveys and metagenomic data from environmental samples have greatly expanded the tree of life and changed our understanding of biogeochemical cycling in the environment. For example, these approaches have led to the expansion of the tree of life as seen in candidate phyla radiation (CPR) and Asgard archaea [2,3,4].

The Gemmatimonadota phylum (formerly BD group; KS-B group; Gemmatimonadetes) was initially described using isolates [5] and recovery of 16S rRNA gene sequences from deep-sea sediment, soil, and bioreactor sludge [5,6,7,8,9]. Both rRNA sequences and metagenomic data have revealed Gemmatimonadota are ubiquitous residents of various environments (soils, freshwater, wastewater treatment plants, and oceans) with relative abundances at around 1% [10]. Prior work on Gemmatimonadota mostly focused on terrestrial environments and suggested that Gemmatimonadota prefer dry soil environments, with high proportions of Gemmatimonadota identified in semiarid/arid soil and deserts [11,12,13]. Their abundance increased during drought conditions [14], suggesting an active response to the change in soil moisture. In coastal sediments, Gemmatimonadota showed strong correlations with different key genes involved in sulfur, nitrogen, and oxygen metabolism [15]. At present, Gemmatimonadota contains six cultured species which are chemoorganoheterotrophs, two of which have the purple bacterial reaction centers for anoxygenic photosynthesis [16, 17]. Gemmatimonadota aurantiaca is capable of reducing N₂O, a potent greenhouse gas [18]. Gemmatimonadota isolates are naturally resistant to some antibiotics, such as ampicillin, penicillin, and chloramphenicol [17, 19, 20]. Moreover, metagenomic data from soil samples shows that they have many biosynthetic gene clusters (BGCs) [21]. Collectively, these studies suggest Gemmatimonadota likely play a significant role in the environment.

Despite the previous findings about the metabolic potential and ecology of Gemmatimonadota, their ecological patterns and metabolic diversity remain unclear in the ocean. To address this knowledge gap, we obtained 495 Gemmatimonadota metagenome-assembled genomes (MAGs) from five different marine environments. An updated phylogeny revealed four distinct phylogenetic groups and seven distinct clusters based on protein composition. These marine Gemmatimonadota possess genes involved in different processes of biogeochemical cycling, including carbon, nitrogen, and sulfur, as well as biosynthetic gene clusters (BGCs) for secondary metabolites to regulate the microbial community. This sheds new light on the genomic diversity and ecological roles of Gemmatimonadota in marine environments, which are distinct from terrestrial genotypes.

Results and discussion

Phylogeny of newly constructed genomes and distribution of Gemmatimonadota

Gemmatimonadota is monophyletic with the Fibrobacterota, Chlorobi, and Bacteroidota (FCB) superphylum [22, 23]. In this study, 495 Gemmatimonadota MAGs with completeness > 50%, and < 10% single gene duplications (based on CheckM [24]) were reconstructed from coastal sediments in the Bohai Sea, China (BS, 427 MAGs); coastal sediments in San Francisco Bay (SFB, 26 MAGs), USA; hydrothermal sediments in Guaymas Basin (GB, 31 MAGs), Gulf of California, Mexico; biomat and Fe oxyhydroxide precipitating from low-temperature hydrothermal fluid (Fe oxyhydroxide) samples of Longqi hydrothermal vents in the Indian Ocean (IO, 7 MAGs); and cold-seep sediments in the South China Sea (SCS, 4 MAGs), China (Supplementary Table 1). These bacteria represent < 5% relative abundance in the metagenomic assembled community in deep-sea environments (GB, SCS, and IO) (see “Methods,” Supplementary Fig. 1). However, they are more abundant in coastal environments, specifically up to ~ 11% and 16% in SFB and BS, respectively. This may correlate with their ability to catalyze denitrification, which is known to be a dominant process in SFB and BS sediments [25]. Interestingly, the relative abundance of Gemmatimonadota increased with depth at all three sampling stations in BS, while it decreased with depth in GB (Supplementary Fig. 1).

The 495 MAGs were classified as Gemmatimonadota in Genome Taxonomy Database (GTDB) (Release 202, Supplementary Table 2), and this was confirmed by a maximum likelihood tree based on a concatenated alignment of 120 bacterial marker genes defined in GTDB-Tk (Supplementary Fig. 2). Glassbacteria, a phylum curated in NCBI, was classified as a class within Gemmatimonadota in GTDB. All Gemmatimonadetes genomes in NCBI are classified as Gemmatimonadetes class within the Gemmatimonadota phylum in GTDB. Based on the phylogenetic tree, these Gemmatimonadota MAGs were split into four groups: Group 1, Group 2, Group 3, and Group 4. These correspond to the JACCXV01, Gemmatimonadales, KS3-K002, and Longimicrobiales orders in GTDB. The genome sizes of the 495 MAGs range from 1.43 to 9.92 Mbp (average 3.68 ± 1.15 Mbp) (Supplementary Table 2). The wide range of genome sizes is likely associated with their evolution and ecology, rather than genome completeness. In support of this, different genome size ranges are associated with distinct habitats. For example, Group 1 has the smallest average genome size and was only recovered from deep layers (30–62 cm) in BS samples (Supplementary Fig. 1). Group3 has a wider range of genome sizes than Group 1 because of their prevalent distribution in different layers in BS samples, some GB samples, and Great Barrier Reef samples. Group 2 and Group 4 are distributed in diverse marine and terrestrial habitats and have a wider range of genome sizes than the other two groups, which were only recovered from marine environments (Supplementary Fig. 1).

We selected 245 MAGs which have completeness > 80% and contamination < 5%, ranging from 2.24 to 6.67 Mbp (average 3.91 ± 0.99 Mbp) for phylogenomic analyses (Supplementary Table 2). To further confirm the taxonomy of these MAGs, we constructed two phylogenetic trees (Fig. 1 and Supplementary Fig. 3) of the 245 MAGs and 211 reference genomes. The phylogenies were constructed based on the concatenated protein alignment of 120 single-copy markers in GTDB (Fig. 1) and 37 concatenated ribosomal protein encoding genes identified using PhyloSift (Supplementary Fig. 3; see “Methods”). Both trees supported the classification of the four groups, which contained 10, 85, 100, and 50 MAGs for Group 1, Group 2, Group 3, and Group 4, respectively. However, the 37 marker gene trees showed that Group 1 was phylogenetically closer to Glassbacteria than the rest of the three groups (Supplementary Fig. 3).

Fig. 1

Phylogeny of Gemmatimonadota and an overview of their metabolic potential. A maximum likelihood phylogenetic tree (IQ-TREE, based on concatenation of 120 single-copy proteins in GTDB) of 456 genomes including the 245 metagenome-assembled genomes (MAGs) described in this study. The four groups are marked in different background colors with black dots indicating the newly recovered MAGs, the inner ring indicating the environmental source of each genome, and the outer ring indicating seven protein clusters derived by using metagenomic entropy-based scores (MEBS) protein (Pfam) content of each genome. The metabolic potential of newly reconstructed genomes is shown in the heatmap for nitrogen (N), iron (Fe), oxygen (O), carbon (C), and sulfur (S), determined by using MEBS. The outer bars represent the number of biosynthetic gene clusters (BGCs) per genome. Bootstraps are shown in grey circles (≥ 75)

Average amino acid identity (AAI) analysis revealed that Gemmatimonadota MAGs are distinct from other phylogenetically related phyla (at most 45.9% identity to Fibrobacterota and 49.8% identity to Glassbacteria) (Supplementary Table 3 and Supplementary Fig. 4). AAI supported the classification of the four groups, which share a maximum 59.7% AAI between each other (Supplementary Table 4). The 16S rRNA gene phylogeny we constructed here is generally consistent with those previously reported [10] and the 37-marker ribosomal protein and 120-marker phylogenies from this study (Fig. 1). 16S rRNA genes in Group1 are classified as the class AKAU4049 (JACCXV01 order in GTDB) (Supplementary Fig. 5). Group 2 and Group 3 16S rRNA genes are classified as classes Gemmatimonadaceae and PAUC43f marine benthic group, which correspond to GTDB orders Gemmatimonadales and KS3-K002, respectively. Group 4 16S rRNA genes belong to the BD2-11 terrestrial group (Supplementary Fig. 5). However, based on our 37-marker ribosomal protein and 120-marker phylogenies, we suggest that Longimicrobiaceae, S0134 terrestrial group, and BD2-11 are within the Longimicrobiales order (Group 4, Fig. 1).

The recovery of Gemmatimonadota MAGs in this study further confirms their global distribution [10, 26,27,28,29,30,31,32,33] (Fig. 2). As of Nov. 11, 2021, 324 Gemmatimonadota genomes were available in NCBI. One-hundred fifty-one of these have been recovered from marine environments, including the Pacific Ocean, Atlantic Ocean, and Indian Ocean; however, most Gemmatimonadota to date have been recovered from terrestrial environments (173 of 324 genomes). Group 1 and Group 3 are primarily composed of MAGs recovered in this study.

Fig. 2

The global and environmental distribution of Gemmatimonadota bacteria. a The global distribution of 819 Gemmatimonadota MAGs (324 publicly available genomes from NCBI and 495 MAGs recovered from this study). Colors represent habitats, shapes represent general types for habitats (circle: marine, triangle: terrestrial, and star: study sites in this study), and shape size represents number of genomes obtained from each site. b Distribution of Gemmatimonadota in different habitats with numbers in stars representing the number of genomes recovered from this study and numbers in circles representing the number of Gemmatimonadota genomes obtained from publicly available databases

Protein-level comparison across the Gemmatimonadota

In order to resolve how these bacteria compare at the predicted protein level, we clustered all of the genomes based on their Pfam profiles (see “Methods”). This approach has proven to be an effective way to identify guilds of bacteria that share common ecological capabilities [34]. This revealed these bacteria fall into seven distinct protein clusters. Group 1 forms a unique Pfam cluster, while Group 2 is divided into three clusters (Fig. 1). Half of Group 2 MAGs were deposited in the database, which were mainly recovered from terrestrial environments. However, Group 2 MAGs recovered from this study were phylogenetically distinct with those recovered from terrestrial environments (Fig. 1). Moreover, the distribution of metabolic proteins based on the presence/absence of protein families in these newly recovered Group 2 MAGs was different from those curated in the database (Fig. 1). The interlaced Pfam clusters in Group4 (Fig. 1) together with their diverse habitats suggest more frequent horizontal gene transfer in Group 4 than the other three groups. The worldwide distribution of this phylum reveals they are overlooked and of ecological significance.

Metabolic flexibility of Gemmatimonadota enables their wide distribution in marine environments

To understand the metabolic potential of the 245 MAGs (completeness > 80% and contamination < 5%, Supplementary Table 2), we compared their predicted proteins against six databases (see “Methods”). We determined that metabolic pathways for polysaccharide and detrital protein degradation, as well as sulfur, nitrogen, and iron utilization, are common in these bacteria.

Large molecule organic matter degradation

Gemmatimonadota MAGs encode over 21,000 potential carbohydrate-active enzymes (CAZymes) classified as glycoside hydrolases (GHs), carbohydrate esterases (CEs), polysaccharide lyases (PLs), glycoside transferases (GTs), and carbohydrate-binding modules (CBMs) (Supplementary Tables 5 and 6 and Supplementary Fig. 6). Among ~ 11,100 glycoside hydrolases (GHs), carbohydrate esterases (CEs), and polysaccharide lyases (PL), 115 are predicted to be extracellular, contributing to the degradation of polysaccharides outside the cell. For example, extracellular CAZyme genes belonging to families/subfamilies CE1, GH16_3, and GH18, contributing to the degradation of xylan, chitin, and laminarin, were commonly identified in Groups 2, 3, and 4. Several MAGs in Group 2 and Group 4 encode multiple types of extracellular CAZymes (CE1 and GH16_3 in Group 2; GH13_32 and GH16_3 in Group 4) for the degradation of different substrates, e.g., pectin and laminarin (Supplementary Fig. 7). One MAG in Group 3 (M8-44_Bin_110) encoded five types of extracellular CAZyme genes (GH16_3, GH30_1, GH136, PL31, and GH0), enabling them to degrade both complex and relatively simple carbohydrates, e.g., degrading laminarin by laminarinase (GH16) or releasing lacto-N-biose from oligosaccharide by lacto-N-biosidase (GH136). The released monosaccharides could benefit the community as a whole by supplying organic matter to other microorganisms.

Similar to CAZyme genes, Group 1 encodes the least diversity of peptidases. Of over 41,000 identified peptidase sequences, 1150 are predicted to be extracellular, suggesting that detrital proteins are degraded outside the cell and later taken up for consumption. Most of the MAGs recovered here have multiple extracellular peptidase genes (Supplementary Fig. 8). For example, they have genes predicted to produce extracellular peptidases belonging to family M28 (aminopeptidase and carboxypeptidase) [35] and S8 (serine endopeptidase subtilisin) [36]. These are nonspecific peptidases that release amino acids for assimilation or dissimilation (Supplementary Fig. 8). Family M4, which are primarily secreted peptidases [36], were identified across three groups (Groups 2, 3, and 4) for degrading extracellular proteins and peptidases. The wide distribution of extracellular peptidase genes in marine Gemmatimonadota suggests these bacteria are important players in the degradation of detrital proteins. Additionally, these marine Gemmatimonadota also encode key genes for the transport, activation, and cleavage of fatty acids through beta oxidation [37] (Fig. 3). The capability of degrading different types of large molecules, especially those extracellular degradations which release more readily degradable substrates, suggests that Gemmatimonadota may provide simple energy sources to support the entire microbial community.

Fig. 3

Overview of biogeochemical metabolic pathways in the four phylogenomic groups of Gemmatimonadota. Within each color wheel, colored segments, gray, and blank segments represent gene presence in over 50%, less than 50%, and gene absence, respectively, within a group. Red-dashed arrows indicate only partial of the known subunits are present in the MAGs. Gene annotations were based on KEGG assignments as summarized in Supplementary Table 9. The number and letter inside each circle represent the pathways in Supplementary Table 9

Nitrogen, sulfur, and hydrogen cycling

Metabolic inference using Multigenomic Entropy-Based Scores (MEBS) [38] (see “Methods”) indicates Gemmatimonadota have pathways for nitrogen and sulfur utilization (Fig. 1, Supplementary Table 9). Two-hundred three of 245 MAGs, belonging to all four groups, are capable of incomplete denitrification and encode genes for the reduction of nitrate, nitrite, and nitrous oxide (N₂O), but not nitric oxide, as well as the oxidation of hydroxylamine (NH₂OH) to nitric oxide (NO). Most MAGs (201/245) encode membrane-bound nitrate reductase (NarGHI, present in all four groups but only in over 50% MAGs in Group1) (Supplementary Fig. 9) and/or periplasmic nitrate reductase (NapAB, present in Group 2, Group 3, and Group 4) (Supplementary Fig. 10) suggesting Gemmatimonadota play a key role in nitrate reduction, the first step of denitrification. MAGs that encode NarG were commonly recovered from deep BS sediments (below 30 cm) and were rare at other depths. Some MAGs (40/201, belonging to Group 2, Group 3, and Group 4) encode both NarG and NapA, and these are predominantly from the BS (28–30 cm and 42–44 cm at M3 and 42–44 cm and 56–62 cm at M8). A phylogeny of NarG in Gemmatimonadota MAGs indicates that NarG is monophyletic and thus may have been present in the last common ancestor of Gemmatimonadota (Supplementary Fig. 9). In contrast to the widespread presence of genes for nitrate reduction in Gemmatimonadota, dissimilatory nitrite reduction via NirK/S (62/245 MAGs) for NO production and NrfAH (8/245 MAGs) for ammonia production is less common in these bacteria (Supplementary Table 9). The 64 NirK/S genes occurred in all four groups from all the sampling sites, while NrfAH was mainly distributed in Group 4 recovered from BS and SFB (Supplementary Table 9). For denitrification in BS and SFB, Gemmatimonadota likely relies on metabolic handoffs to complete denitrification, due to the lack of nitric oxide reductase, reducing nitric oxide to nitrous oxide. All Gemmatimonadota groups (125/245 MAGs) encode genes for periplasmic nitrous oxide reductase (NosZ), which reduces N₂O to N₂ (Supplementary Table 9). Phylogenetic analyses show that Gemmatimonadota NosZ are atypical type NosZ sequences (Supplementary Fig. 11), which are associated with microorganisms that are not complete denitrifiers [39, 40]. N₂O is a potent greenhouse gas and degrades ozone in the atmosphere [41]. The wide distribution of Gemmatimonadota suggests that they may have key roles in reducing N₂O fluxes in marine environments [42]. Moreover, Gemmatimonadota have different transporters for small molecules, including nitrate, nitrite, and ammonium (Fig. 3). Gemmatimonadota has also been reported to hydrolyze urea as an energy source in wastewater treatment sludge [43, 44]. However, we only identified urease (UreABC) from a single MAG in Group 4 (M3-22_Bin_219), suggesting this is not important in marine environments. Collectively, our findings suggest Gemmatimonadota may play an important role in nitrogen cycling in marine sediments, especially in the coastal zones.

Unlike denitrification genes that are prevalent in all Gemmatimonadota groups, sulfur cycling genes are limited to specific Gemmatimonadota groups (Fig. 1). A clade within Group 3 (21/100 MAGs), which has a unique protein composition (cluster1 in Fig. 1), has gene clusters for sulfate reduction (including DsrAB, SAT, AprAB, and QmoABC). Some Group 2 MAGs (16/85) also appear to be capable of reducing sulfate. These DsrAB sequences were mostly associated with MAGs recovered from deep sediments (below 30 cm) in BS, where high concentrations of sulfate were detected (> 22 mmol/L in pore water) [15]. DsrAB genes (Supplementary Figs. 12 and 13) do not appear to have been horizontally transferred from different phyla [45], suggesting sulfite reduction may be an ancient function within Gemmatimonadota.

Three deep branching groups, consisting of five MAGs from SFB, four from BS, and one from BS in Group 2 (10/85 MAGs), contain genes encoding sulfhydrogenase I complex (HydADGB) [46] for coupling sulfur reduction with H₂ oxidation (Supplementary Table 9). However, the majority of Gemmatimonadota are capable of oxidizing different sulfur substrates, e.g., sulfide and sulfite. Specifically, 167/235 MAGs (excluding 10 Group 1 MAGs) have sulfide-quinone oxidoreductase (SQR) for sulfide oxidation, and 14 of those 167 MAGs encode both DsrAB and SQR (Supplementary Table 9). Phylogenetic analysis indicates that these SQRs belong to the membrane-bound type I, type II, and type III SQRs (Supplementary Fig. 14). Interestingly, 194/245 MAGs have homologs to eukaryotic thiosulfate/3-mercaptopyruvate sulfurtransferase (TST) [47], which could transfer thiosulfate and cyanide to sulfite and thiocyanate. Group 1 (8/10 MAGs) has homologs to eukaryotic sulfite oxidase (SUOX) [48, 49], a type of molybdopterin-dependent oxidoreductase, for sulfite oxidation with oxygen as the electron acceptor (Supplementary Table 9). However, Group 1 was recovered from the deep layer of BS sediments (below 30 cm), suggesting that this sulfite oxidase may also use cytochrome c as the final electron acceptor [49]. Interestingly, all Group 1 MAGs (10/10), as well as 2/85 Group 2 and 5/100 Group 3 MAGs recovered from deep layers of BS sediments (below 30 cm), encode methanethiol oxidase to aerobically oxidize methanethiol. Methanethiol is a key intermediate for global organosulfur compounds, e.g., dimethylsulfoniopropionate (DMSP) and dimethyl sulfide (DMS) cycling [50, 51]. Moreover, 34 (25/85 Group 1, 8/100 Group 2, and 1/50 Group 4) MAGs are predicted to produce DMS via methanethiol S-methyltransferase (MddA) from methylate L-methionine or methanethiol (MeSH) under oxic conditions [52]. In addition, 117 (62/85 Group 2, 42/100 Group 3, and 13/50 Group 4) MAGs encode genes for the large subunit of thiosulfate dehydrogenase (DoxD), which may convert thiosulfate to tetrathionate. Six (1/85 Group 2, 1/100 Group 3, and 4/40 Group 4) MAGs also have genes for the catalytic subunit of tetrathionate reductase (TtrA), which may reduce tetrathionate to thiosulfate. Thus, Gemmatimonadota likely play important roles in a variety of intermediate steps in marine sulfur cycling (Fig. 3).

Hydrogen metabolism is crucial in energy cycling in marine environments [53]. Gemmatimonadota, except for Group 1, have different types of [NiFe] hydrogenases (Supplementary Fig. 15) and few [FeFe] hydrogenases (mainly in Group 2) (Supplementary Fig. 16), suggesting hydrogen is coupled to metabolic pathways in these bacteria [54, 55]. Hya hydrogenase (HyaABCD, [NiFe] type) was widely distributed in Group 2, Group 3, and Group 4 (Supplementary Table 9). Hya hydrogenase is resistant to oxidative stress (e.g., superoxide and hydrogen peroxide), which may enable Gemmatimonadota to oxidize H₂ in the presence of oxygen [56, 57]. A subgroup of Group 3 and Group 2 also have the F420-nonreducing hydrogenase (MvhADG), belonging to Group 3c [NiFe] hydrogenase (Supplementary Fig. 17). This F420-nonreducing hydrogenase links with heterodisulfide reductase (HdrABC), by providing reducing equivalents without reacting with F420, i.e., transporting electrons using H₂ as an electron donor [58]. Additionally, 10/100 Group 3 MAGs have the HoxFHUY operon (Supplementary Fig. 17), a bidirectional [NiFe] hydrogenase mainly described in Cyanobacteria. The Hox operon serves as a regulator for maintaining a proper redox state in the cell [59], which could be important for the metabolic versatility of Gemmatimonadota.

Iron, mercury, and arsenic utilization

Microbially mediated iron cycling has been linked with many crucial marine processes, such as carbon storage, greenhouse gas emission, and primary production in the ocean [60]. We identified a variety of genes potentially involved in cryptic iron cycling in marine Gemmatimonadota, including iron acquisition, storage, oxidation, and reduction. Specifically, we identified two clusters of MAGs in Group 3 and Group 4 recovered from a wide range of depths (ranging from 0 to 62 cm) in BS sediments encoding genes for sulfocyanin [61, 62] (Supplementary Fig. 18), a putative iron oxidase. These bacteria may link iron oxidation via sulfocyanin with nitrate reduction via periplasmic nitrate reductase (NapAB) or nitrous oxide reduction via nitrous oxide reductase (NosZ) (Supplementary Fig. 18). We also identified three MtrABC operons [63] (M3-44_Bin_97, M3-38_Bin_128, and M3-30_Bin_133) in Group 2, suggesting they may be capable of reducing iron in anoxic sediments (Supplementary Fig. 18). Other widely annotated potential iron cycling gene homologs, such as cytochrome-c Cyc2 and DFE_461-465, in these MAGs suggest that Gemmatimonadota may actively participate in iron cycling; however, it is difficult to distinguish iron reduction and iron oxidation based on the current annotation.

Gemmatimonadota also encodes mercury and arsenic detoxification systems. They are capable of transforming the extremely toxic Hg(II) to metallic Hg(0) via mercuric reductase (MerA), potentially detoxifying mercury (Fig. 3). All four groups are capable of reducing arsenate to arsenite via arsenate reductase (ArsC) through thioredoxin [64] (Fig. 3). Resistance and detoxification of heavy metals may enable Gemmatimonadota to be widely distributed from coastal sediments to deep oceans [65], where Hg and As have accumulated from anthropogenic pollution [66, 67] or released via hydrothermal activity and volcanic eruptions [68, 69].

Extensive genetic potential for secondary metabolite biosynthesis in Gemmatimonadota

Microorganisms produce secondary metabolites to interact with other community members and their environment. The importance of biosynthetic gene clusters (BGCs) in Gemmatimonadota has been described in soil environments [21]. This has not been examined in marine Gemmatimonadota [70] due to the limited representatives in public databases. We identified a diverse genetic potential for secondary metabolite biosynthesis, including nonribosomal peptide synthetase (NRPSs) and polyketide synthases (PKSs) (Fig. 1). Combined gene clusters consisting of different NRPS, PKS, and hybrid NRPSs/PKS were identified in 69 MAGs in four groups (Fig. 1). NRPS and PKS are known to synthesize a diversity of antibiotics, antifungals, and immunosuppressants with pharmaceutical potential [71], while the majority of these NRPS and PKS have unknown end products [72].

The most common type of BGCs identified in Gemmatimonadota is involved in the biosynthesis of terpenes, including carotenoid, isorenieratene, and N-tetradecanoyl tyrosine, and was found in 174 MAGs in these bacteria. Terpenes can have antibacterial properties [73], participate in bacterial-fungal interactions [74], and provide colorful pigments [75]. However, the ecological functions of different terpenes remain poorly understood. BGCs encoding lasso peptides, a class of ribosomally synthesized and posttranslationally modified peptides (RiPPs) [76, 77], were identified in 29 MAGs mainly from Group 3 (Fig. 1). The antibacterial properties of lasso peptides produced by Gemmatimonadota suggests a potential role of affecting the abundance of the other community members. Bacteriocins (TIGR03798, Nif11-related peptide) experience intensive posttranslational modifications to generate antimicrobial peptides which are toxic to the strains of closely related species [78]. Genes encoding bacteriocin have been particularly prominent in Gemmatimonadota in soil environments [79]. We annotated genes for microcin, a type of bacteriocin [80], in 20 MAGs exclusively within Group 2 (Fig. 1). Specifically, six MAGs (recovered from below 30 cm at station M3 and M8, BS) within a monophyletic group in Group 2, have multiple copies of microcin genes that may mediate Gemmatimonadota population size [81].

A broad diversity of Gemmatimonadota have the potential to produce different secondary metabolites, which may play a critical role in the survival and adaptation of the microbial community and result in their prevalence across different habitats. Perhaps most strikingly, there are clades in Group 4 associated with corals that are enriched in bacteriocin, terpene, and type I polyketide synthase (T1PKS) genes (Fig. 1). These genotypes have unique protein composition comprising Pfam Cluster 7, suggesting these bacteria use secondary metabolites to interact with other organisms in reef communities. BGCs with low levels of similarity to known databases can be used to mine novel BGCs and point to new compounds [82]. The Gemmatimonadota phylum may thus represent a reservoir for the discovery of secondary metabolites, which could also be useful in medicine and biotechnology.

Potential Gemmatimonadota viruses

In total, 6,611 double-stranded DNA (dsDNA) viral metagenome-assembled genomes (vMAGs) of high- and medium-quality were identified from 15 BS samples (see “Methods”). We identified three CRISPR-Cas systems (Supplementary Fig. 19) and 639 CRISPR spacer sequences (Supplementary Table 10) in 156 of 245 Gemmatimonadota MAGs. However, only one Gemmatimonadota could be linked with vMAGs via the CRISPR spacer sequences. Using CRISPR spacers, tRNA matching, 6-mer oligonucleotide frequency, and whole genome matching, we identified 32 vMAGs ≥ 10 kilobases in length that potentially infected Gemmatimonadota (see “Methods”) (Fig. 4 and Supplementary Table 11). Among these are 15 viruses that could not be assigned taxonomy, while the other 17 of the 32 viruses were classified as Caudovirales belonging to Myoviridae (12), Podoviridae (2), and Siphoviridae (3) (Supplementary Table 12). However, none of these viruses was clustered with known viral genomes at the genus level based on shared-gene content (Fig. 4). To understand the viral roles in host metabolism, we assigned functions to the gene content of these 32 vMAGs, revealing a variety of putative auxiliary metabolic genes (AMGs) that may “hijack” and manipulate host metabolism (Fig. 4 and Supplementary Table 13). We identified D-beta-D-heptose 7-phosphate kinase in three Myoviridae viruses from different BS samples, mainly associated with Group 3 and Group 2, suggesting these viruses may contribute to the assembly of the lipopolysaccharide in their hosts [83]. In addition, one unknown taxonomy virus and four Myoviridae viruses encode heptosyltransferase, a class of glycosyltransferases (GTs) that may modify heptose residues on lipopolysaccharides to affect viral-host interactions [84, 85].

Fig. 4

Genomic diversity and composition of viruses that have infected Gemmatimonadota. a, b vMAG statistics, including sampling location, quality, length, lifestyle, and number of genes. c A gene-sharing network of viral sequences between vMAGs from this study and prokaryotic virus genomes in the GenBank database. Nodes represent individual genomes, and edges indicate similarity among genomes within a viral cluster. d Number of vMAGs identified using different tools. e Annotation of AMGs in 27 vMAGs. f Potential association between vMAGs and Gemmatimonadota genomes

We also identified viral genes involved in genome replication, nucleotide metabolism, and posttranscriptional modifications, including ribonuclease H (RNaseH-like domain), ATP-dependent DNA ligase (ligD) [86], and peptidases (Supplementary Table 14). Four vMAGs contain genes predicted to encode ribonucleotide reductase, which is important for nucleotide metabolism in nucleocytoplasmic large DNA viruses (NCLDVs) [87]. Putative Gemmatimonadota viruses also contain genes for DNA methylation and glycosylation that may be important for host interactions. We identified genes for methyltransferase and endonuclease (Supplementary Table 14), suggesting the viruses may be involved in epigenetic modification via autonomous DNA methylation [88]. Finally, we identified genes encoding pyruvate-ferredoxin oxidoreductase in one vMAG (Supplementary Table 14), suggesting that this virus may contribute to host anaerobic metabolism by generating acetyl‐coenzyme A, carbon dioxide, and reduced ferredoxin (Fd^2-) [89].

Potential phototrophic and autotrophic capabilities in Gemmatimonadota

Gemmatimonadota in order Gemmatimonadales have recently been shown to possess photosynthetic gene clusters (PGCs) [16]. However, none of the newly reconstructed MAGs recovered here codes for PGCs, suggesting horizontal gene transfer of PGCs is not common among Gemmatimonadota [16] (Supplementary Fig. 20). We did not identify any key genes for the type II photosynthetic reaction center (puf, bch, and acsF genes) in our MAGs, as found in terrestrial environments, e.g., the isolate from freshwater Swan Lake in the Gobi Desert in China, the Cock Soda Lake, and Lake Baikal in Siberia [16, 90, 91]. Therefore, marine Gemmatimonadota appear to lack phototrophic metabolism. However, as stated above, Gemmatimonadota encode bacteriocins (TIGR03798, Nif11-related peptide) and carotenoids, which are associated with photosynthetic Cyanobacteria [80, 92], and the latter is thought to contribute to adaptation to low light conditions [93] or UV exposure [94]. Thus, the phototrophic metabolism may be occurring in shallow marine environments.

Carbon fixation genes via Calvin–Benson–Bassham (rbcS, rbcL, and prk) have been reported in Gemmatimonadota from soda lakes [91, 95]. However, the soda lake Gemmatimonadota MAGs are phylogenetically distinct from our marine groups (Supplementary Fig. 20). Moreover, only the large subunit of ribulose 1,5-bisphosphate carboxylase/oxygenase-like protein (RLP, form IV RuBisCO), potentially important for sulfur metabolism rather than CO₂ fixation (Supplementary Fig. 21), was annotated in 43 MAGs in this study. Additionally, we did not find any complete autotrophic pathways (Wood-Ljungdahl pathway, Calvin–Benson–Bassham, reductive tricarboxylic acid, 3-hydroxypropionate bicycle, 3-hydroxypropionate-4-hydroxybutyrate, and dicarboxylate-4-hydroxybutyrate cycles) in our marine Gemmatimonadota MAGs (Supplementary Table 9). There has been no physiological confirmation of autotrophic metabolism in Gemmatimonadota, and thus, they are likely heterotrophs.

Ecology of Gemmatimonadota

Gemmatimonadota is estimated to be the eighth most abundant bacterial phylum in soils, with relative abundance of ~ 1% of soil bacteria worldwide [26]. They are globally distributed with low abundance (< 2%) in marine environments [96, 97] and are estimated to be over 10% relative abundance in deep-sea sediments [98]. Marine clades are phylogenetically distinct from terrestrial clades, where Group 1 and Group 2 members described in this study are distinct from their terrestrial sister groups, and Group 1 was only recovered from deep sediments (38–62 cm) in two sampling sites. This suggests a potential unique ecological role of marine Gemmatimonadota. The marine genotypes described here are metabolically diverse, and many are capable of partial denitrification and organic carbon degradation. A diversity of nitrous oxide reductases suggests marine Gemmatimonadota may mediate the reduction of nitrous oxide to nitrogen gas for the removal of the most potent greenhouse gas, which is increasing due to increased anthropogenic activities [41, 99, 100], and a vital process in ocean biogeochemistry. These organisms encode proteins for the degradation of different complex carbon compounds, including pectin, laminarin, and fatty acids. Marine pectin and laminarin are produced by photosynthetic marine microalgae [101], diatoms, macrophytes [102], and terrestrial plants [103]. Thus, Gemmatimonadota are likely players in organic matter degradation in the oceans. Also, the protein repertoire of these MAGs suggests they participate in arsenic and mercury cycling/detoxification. Interestingly, there are clades associated with coral reefs that are enriched in BGC genes; these genotypes also have unique protein profiles (Pfam Cluster 7, Fig. 1). This suggests that they produce metabolites for interactions in reef ecosystems.

The prevalence of Gemmatimonadota across various terrestrial environments has been shown in several studies [9, 11,12,13,14]. However, the metabolic potential and ecological roles of Gemmatimonadota in the ocean are poorly understood due to a lack of genomic sampling. Gemmatimonadota have versatile metabolisms and high abundance in coastal areas where they appear to be involved in the degradation of complex organic carbon, denitrification, sulfate reduction, and sulfide/sulfite oxidation. Interestingly, marine genotypes are distinct in their numbers of BGCs, as well as sulfur and iron metabolic genes. The expanded genomic biodiversity provided in this study is a framework to understand the roles of Gemmatimonadota on a global scale.

Materials and methods

Sampling procedures

The 495 Gemmatimonadota MAGs in this study were reconstructed from five marine environments including the following: coastal sediments in the Bohai Sea (BS), China (427 MAGs); coastal sediments in the San Francisco Bay (SFB), USA (26 MAGs); hydrothermal vents in the Guaymas Basin (GB), Gulf of California, Mexico (31 MAGs); biomat and Fe oxyhydroxide precipitating from low-temperature hydrothermal fluid (Fe oxyhydroxide) samples of Longqi hydrothermal vents in the Indian Ocean (IO) (seven MAGs); and cold-seep sediments in the South China Sea (SCS), China (four MAGs) (Supplementary Table 2). BS samples were collected from sediments at three stations: BHB10, M3, and M8 in the BS, China. Details of sampling sites and procedures were described previously [15]. SFB samples were collected from sediments in the San Francisco Bay, USA, across four seasons (July, October, January, and May) between July 2011 and May 2012 [104]. GB samples were collected from sediments in the Gulf of California, Mexico, during cruises in 2008 and 2009. Backgrounds of GB and the sampling details were described previously [34]. Details of the procedure of obtaining hydrothermal vent mat and Fe oxyhydroxide samples from the Indian Ocean and cold seep sediment samples from the SCS were described by Gong et al. [105]. Information about the samples from each environment is summarized in Supplementary Table 1.

DNA extraction, metagenomic sequencing, assembly, and binning

Total DNA was extracted from each sample and sequenced after quality control as described in the corresponding references [15, 34, 104, 105]. The sequences were assembled and binned with different protocols based on the batch of the samples [34, 105]. In total, 495 MAGs with completeness greater than 50% and contamination less than 10% with CheckM lineage_wf v1.0.5 [24] were identified as Gemmatimonadota based on the taxonomy assigned using GTDB-Tk v1.5.1 [106] with release 202. The relative abundance of each MAG was calculated using MetaGaia with default settings (https://github.com/valdeanda/MetaGaia). All reconstructed MAGs were grouped into five genome datasets based on five study sites. Raw reads from each sample were mapped to the matched genome dataset. The relative abundance of Gemmatimonadota in each sample is the ratio of reads mapped to Gemmatimonadota to the reads mapped to the genome dataset with normalized genome size.

Phylogenetic analyses

A concatenated protein sequence alignment of 120 marker genes was generated with GTDB-Tk v1.5.1 [24] from a set of genomes consisting of 495 Gemmatimonadota MAGs in this study and 356 reference genomes downloaded from NCBI (Supplementary Table 2). All selected reference genomes for phylogenetic analysis were with completeness greater than 50% and contamination less than 10% confirmed with CheckM lineage_wf v1.0.5 [24]. The alignment was refined using MAFFT v7.471 [107] with option –auto and trimmed using TrimAL v1.4 with the option − gappyout [108]. The alignment was manually checked to ensure the short sequences with large gaps were removed before the construction of the phylogenetic tree. A maximum likelihood tree was built using IQ-TREE v1.6.12 [109] with the parameters: -m LG + C60 + F + R -bb 1000 -nt AUTO -bnni -alrt 1000.

Only MAGs and reference genomes (defined here as MAGs and reference genomes with completeness > 80% and contamination < 5%) were used for further phylogenetic and metabolic analyses. In the end, 245 MAGs in total, i.e., 228, 9, 5, 2, and 1 MAGs from the BS sediments, the SFB sediments, GB sediments, biomat and Fe oxyhydroxide from the IO, and cold-seep sediments from the SCS, respectively, and 211 reference genomes were used for downstream analysis. To further resolve the phylogeny of Gemmatimonadota, the concatenated alignment of 120 marker genes was generated with GTDB-Tk v1.1.1 from the 245 Gemmatimonadota MAGs and reference genomes, except for a genome of Longimicrobium terrae strain: CB-286315 (GCA_013000925.1) with a completeness of 98.90%, and contamination of 7.69% was included due to the lack of representatives for this group. The alignment was refined and trimmed, and the phylogenetic tree was built the same way as described above. Another set of alignment with 37 universal single-copy marker genes encoding proteins was extracted using PhyloSift v1.0.1 [110] to further confirm the phylogenetic placement of the MAGs. The alignment was refined and trimmed, and the phylogenetic tree was built the same way as described above.

16S rRNA gene sequences were extracted from MAGs using Barrnap v0.9 (https://github.com/tseemann/barrnap) with the following parameters: –reject 0.03. Sequences with length over 700 bp were uploaded to ARB v6.0.6 [111] and aligned with the reference database (release 138) [112]. The alignment was manually checked, exported from ARB v6.0.6, and trimmed using trimAl v1.4. The alignment with 184 sequences in total: 35 16S rRNA sequences from MAGs and 149 from the reference database was used to build a maximum likelihood tree using IQ-TREE v1.6.12 with the following parameters: -m MFP -bb 1000 -bnni -nt AUTO. Trees were uploaded to Interactive Tree Of Life (iTOL) v5 [113] for visualization.

The average amino acid identity (AAI) values of 245 MAGs and 211 reference genomes were calculated using CompareM v0.1.2 (https://github.com/dparks1134/CompareM) to create an AAI matrix.

The distribution and habitat of Gemmatimonadota genomes were based on metadata of genomes downloaded from NCBI (Supplementary Tables 1 and 2).

Metabolic annotations

Gene predictions were performed using Prodigal v2.6.3 [114] with the default settings. Predicted genes were annotated using standalone version KofamScan v1.3.0 with the settings: –e-value 1e-5 [115] and further using the KAAS (KEGG Automatic Annotation Server) web server [116] with the “complete or draft genome” setting with parameters: GHOSTX, custom genome dataset, and BBH assignment method. Multigenomic Entropy-Based Score (MEBS v1.0) [38] was used with the mebs.pl script using the -comp option to scan against the Pfam v3.0 database to determine the presence/absence of a protein. The hierarchical cluster was based on the mebs_clust.py script with the parameters: –method ward –cutoff 0.5. The score of metabolism of sulfur, carbon, oxygen, iron, and nitrogen was calculated using the mebs_vis.py scripts in MEBS and visualized by importing the output to the iToL. To visualize the sharing of genomic contents between the genomic cluster and taxonomic classification, a Sankey chart was generated online (https://powerbi.microsoft.com/) based on the cluster produced by MEBS.

Iron-related genes were further identified using FeGenie v1.0 with default settings [117]. Biosynthetic gene clusters (BGCs) were identified using antiSMASH v5.1.2 [118] with the following parameters: –cb-general –cb-knownclusters –cb-subclusters –asf –pfam2go –smcog-trees –genefinding-tool prodigal. Carbohydrate-active enzymes (CAZymes) were annotated using dbCAN2 [119] against database v9 [120] with default settings. Peptidases were annotated using DIAMOND v2.0.11 [121] with the parameters: -e 1e-10 –subject-cover 80 –id 50 against the MEROPS database (release 12.2) [122]. Prediction of the localization of CAZymes and peptidases was performed with PSORTb v.3.0 [123] with the parameters: -n -o terse.

Characterization of different functional genes

Nitrate reductase (NarG and NapA), nitrous oxide reductase (NosZ), dissimilatory sulfite reductase (DsrAB), sulfide-quinone oxidoreductase (SQR), flavocytochrome c sulfide dehydrogenase (FCCD), and hydrogenase were identified using DIAMOND v2.0.11 against corresponding curated databases [34, 124, 125]. Identified sequences were manually checked with the annotation in KofamScan and KAAS. Hydrogenases were further confirmed with the web-based hydrogenase classifier (HydDB) [125]. The large subunit of ribulose-1,5-bisphosphate carboxylase-oxygenase (RuBisCO) (RbcL) sequences were downloaded from NCBI based on previously published data [126]. RbcL sequences in newly reconstructed MAGs were annotated using KofamScan and KAAS confirmed by manually BLASTp in NCBI. Sequences for each gene were aligned with the reference sequences using MAFFT v7.471 or Clustalw v2.1 [127] and trimmed using trimAl v1.4 or BMGE v1.12 [128] (Supplementary Table 15). Maximum likelihood trees were built using IQ-TREE v2.1.2 for each gene. A neighbor-joining phylogenetic tree of hydrogenases was constructed using MEGA 11 [129] as described previously [130]. All phylogenetic trees were visualized using the iTOL. Detailed software and parameters for generating the phylogenetic trees were listed in Supplementary Table 15.

Virus identification, classification, annotation, and virus-host match

We identified viral contigs with length ≥ 10 kb in the assemblies from BS samples using VIBRANT v1.2.1 [131] and VirSorter2 v2.2.3 [132] with settings (-min-score 0.5, -min-length 10,000). Quality of virus genomes was estimated using CheckV v0.8.1 [133]. In total, 6611 high- and medium-quality (> 50% completeness) viral contigs were identified as viral metagenome-assembled genomes (vMAGs) from BS samples. CRISPR arrays and Cas cassette genes in the MAGs were detected using CRISPRCasFinder v4.2.20 [134] with default settings. We further linked the vMAGs with Gemmatimonadota MAGs using four different methods. The first method to predict the virus-host interactions is based on CRISPR spacer match. Briefly, CRISPR arrays in MGAs were identified using CRISPRDetect v2.4 [135], CRISPR Recognition Tool (CRT) v1.2 [136], and CRISPRidentify v1.1.0 [137] with the default parameters. The identified spacer sequences in MAGs were subjected to search using BLASTn v2.11.0 + against the identified viral sequences with the blastn-short algorithm and parameters: identity ≥ 95%, coverage = 100%, word size = 8, and e-value ≤ 10⁻⁵). The second way to link the vMAG with MAG is based on the 6-mer oligonucleotide frequency match between MAGs and vMAG using VirHostMatcher v1.0 [138] with the default settings. The virus-host pairs were selected with the threshold of d₂* < 0.15. The third way is based on the tRNA matching. Briefly, tRNAs in MAGs and vMAGs were predicted using tRNAScan-SE v2.0.11 with the bacterial mode. The recovered tRNAs were compared using BLASTn v2.11.0 + with the threshold (100% identify and 100% coverage). The fourth method to link MAG and vMAGs is based on nucleotide sequence homology between MAGs and vMAGs using BLASTn with the thresholds: coverage of viral contig length ≥ 75%, similarity ≥ 70%, bit score ≥ 50, and maximum e-value ≤ 0.001 [139, 140]. Open reading frames (ORFs) in identified vMAGs were predicted using Prodigal v2.6.3 [114]. All predicted viral protein sequences were annotated using eggNOG-mapper (eggNOG database v5.0.0 [141]), KofamKOALA [115], DRAM-v v.1.2.2 [142], and VIBRANT v1.2.1 [131]. vMAGs were classified using vConTACT2 v0.11.3 [143] with standard parameters: -rel-mode Diamond and -db ProkaryoticViralRefSeq211-Merged. Networks based on shared protein sequences between identified vMAGs and reference prokaryotic viruses were visualized using Cytoscape v3.9.1 [144]. The taxonomy of each vMAG was also assigned using VPF-Class v1.0 [145] with the thresholds: membership ratio ≥ 0.5 and confidence score ≥ 0.2. The viral lifestyle was predicted using PhaTYP [146]. Viral AMGs were identified using DRAM-v v.1.2.2 [142] with the criteria as follows: “M” or “F” flags and auxiliary scores < 4. The heatmap shown in Fig. 4 was generated using TBtools v2.001 [147].