Specimens of P. echinospica and Sclerolinum sp. tubeworms were sourced from the Haima cold seeps in the South China Sea (Since we could not determine the Sclerolinum species, we consistently apply “Sclerolinum sp.” in this study.). The collection of P. echinospica was carried out in April 2023 from a site of 1,385 meters depth (110°27′5,310″E, 16°43′5,017″N). Sclerolinum sp. samples were collected in July 2022 from a depth of 1,428 meters. The Remotely Operated Vehicle (ROV) is powered by the mother ship through a cable, allowing operators aboard the mother ship to remotely control the ROV and collect tubeworms. The tubeworms were placed into an insulated “bio-box” with a sealed lid to minimize temperature fluctuations of the water within the container. Upon boarding the vessel, they were immediately stored in a −80°C freezer.

A single specimen each of P. echinospica and Sclerolinum sp. was processed for metagenomic DNA extraction using (Supplementary Figure S1). The P. echinospica was dissected using sterile instruments, with the plume (a gill-like structure), vestimentum (mainly composed of muscle), and trophosome (an organ housing symbionts) (Supplementary Figure S1). Since the symbiotic bacteria filled the entire trophosome (Nussbaumer et al., 2006), we selected trophosome tissue of approximately 1 cm³ in size for the extraction experiment. Additionally, because the tubes of Sclerolinum sp. are narrow and soft, making direct dissection difficult, we minced the tubes, along with the trophosome tissue, directly for the extraction experiment (Supplementary Figure S1). The samples were thoroughly rinsed with DNA extraction buffer (100 mM Tris–HCl, 100 mM EDTA, 100 mM Na₂HPO₄, 1.5 M NaCl, 1% CTAB, pH 8.0) to eliminate surface impurities and preservatives. The trophosome tissue was immersed in 1 mL of the same buffer, minced with sterile scissors, and incubated at room temperature for 20 min. After the larger tissue fragments had settled, the supernatant containing microbial cells was separated. High-speed centrifugation (12,000 rpm, 13,201 g) was performed at room temperature for 3 min to pellet the microbes for further analysis.

A magnetic bead-based soil and fecal genomic DNA extraction kit (TIANGEN) was used for DNA extraction and purification as per the manufacturer’s guidelines. The purified DNA was then subjected to quantitative and qualitative analysis using a Qubit 4.0 fluorimeter and 1% agarose gel electrophoresis, respectively (Lee et al., 2012). The VAHTS Universal Plus DNA Library Prep Kit for Illumina was employed to prepare metagenomic libraries according to the manufacturer’s instructions. Finally, the metagenomic libraries were sequenced on the Illumina Novaseq 6000 platform using the PE150 strategy at Shanghai Personalbio Bioscience and Technology (Shanghai, China) with 10 Gb for each library.

The raw sequencing data were quality controlled using Fastp v0.23.2 (Chen, 2023) (parameters: -q 20 -u 20 -l 50 -w 40–5 -M 30 -g -D -dup calc accuracy 6) to remove adapters and low-quality reads. The clean reads obtained from each sample were individually assembled using MEGAHIT v1.2.9 (Li D. et al., 2015; Li et al., 2016) (parameters: -t 20 -m 0.96 --no-mercy --kmin-1pass) with k-mer sizes of 21, 29, 39, 49, 59, 69, 79, 89, 99, 99, 109, 119, 129, and 141 bp. The assembled metagenome contigs (≥2,000 bp) were binned using the binning module within MetaWRAP v1.3.2 (parameters: -m 500 -t 40, −c 50 -x 10 –quick), which also selected the best representative genomes through its Bin refinement module (Uritskiy et al., 2018). The quality of the metagenome-assembled genomes (MAGs) was assessed using Quast v5.2.0 (Mikheenko et al., 2018) and checkM2 (Parks et al., 2015) (parameters: --allmodels -x fasta --force -t 8). Prodigal v2.6.3 (Hyatt et al., 2012) (parameters: -o .fna -a .faa -p meta) was employed to predict the proteins and genes for all MAGs. The protein sequences were functionally annotated using the Kyoto Encyclopedia of Genes and Genomes (KEGG) database through KofamScan (v.1.1.0) (Aramaki et al., 2019). Additionally, functional classifications were assigned based on Clusters of Orthologous Groups using eggNOG-mapper v2 (Cantalapiedra et al., 2021).

The genomes of 24 tubeworm symbionts (Figure 1; Supplementary Table S1) downloaded from the National Center for Biotechnology Information (NCBI) database were analyzed and compared with the two symbiont genomes recovered in this study. We used GTDB-TK v2.2.6 with the default settings to perform taxonomic classification for the MAGs (Parks et al., 2018; Chaumeil et al., 2019). The ANIb model in PyAni v0.2.10 (parameters: -m ANIb -g) was used to calculate average nucleotide identity (ANI) values for all genome pairs (Goris et al., 2007; Richter and Rosselló-Móra, 2009). The pheatmap package in R was employed to generate heatmaps. After KEGG annotation, Venn diagrams were created using the Venn tool in the Origin and UpSetR packages to represent shared and unique genes between symbionts. Phylogenetic trees were constructed based on 112 conserved proteins extracted using GTDB-TK and automatically selected best model of IQ-TREE v2.2.6 and branch support estimated using 1,000 ultrafast bootstrap replicates (Minh et al., 2020). The iTOL editor was utilized to visualize the evolutionary trees (Letunic and Bork, 2021). Homologous gene cluster identification of symbionts was performed using Orthofinder (parameters: -t 24 -a 8 -M msa -S blast -A mafft -T iqtree) (Emms and Kelly, 2019).

Pangenome analysis was carried out using the PPanGGOLiN tool¹ to characterize the pangenome across 26 tubeworm symbiont genomes. PPanGGOLiN partitions the nodes using an Expectation–Maximization algorithm based on a multivariate Bernoulli mixture model coupled with a Markov Random Field (Bazin et al., 2020; Gautreau et al., 2020). This approach considers the graph’s topology and the presence or absence of genes in pangenomes to classify tubeworm symbiont gene families into persistent, cloud, and one or several shell partitions (Gautreau et al., 2020). Persistent genomes contain gene families found in most genomes; shell genomes have them at moderate frequencies, and cloud genomes have them at low frequencies. Specific definitions of the three partitions can be found in the Supplementary material.

Tubeworm symbionts, P. echinospica (Hai6C03) and Sclerolinum sp. (SCTW1H), were collected from the Haima cold seeps in the South China Sea to construct a metagenome. Genome quality assessment using CheckM2 indicated that the Hai6C03 symbiont genome assembly exhibited the highest completeness and least contamination compared to other P. echinospica symbiont genomes (Supplementary Table S2). SCTW1H also exhibited the lowest level of contamination compared to similar genomes. The assembled Hai6C03 genome (Figure 2a) had a size of 4,171,724 bp, a GC content of 54%, a completeness of 99.33%, and a contamination of 0.3%. The assembled SCTW1H genome (Figure 2b) had a size of 2,832,488 bp, a GC content of 50%, a completeness of 98.38%, and a contamination of 0.48%.

All tubeworm symbiont genomes were collected from public databases as of July 2023 to analyze the heterogeneity among deep-sea tubeworm symbiont genomes. The ANI results indicate that Vestimentifera tubeworms inhabiting hydrothermal vents in the Pacific Ocean have symbionts that are highly similar to one another. In contrast, Vestimentifera tubeworms from cold seeps exhibit greater variation in their symbiont similarities (see Supplementary Figure S2). Notably, the L. anaximandri tubeworms can live in both hydrothermal and cold seep environments and host different species of symbionts. The ANI results revealed that the symbionts of L. anaximandri found in hydrothermal vents (GCA_016756855 and GCA_016756895) were similar to those in cold seeps (GCA_016756955), with a similarity exceeding 0.98, suggesting they belong to the same species. Additionally, the GCA_016756865 and GCA_016756835 symbionts from the cold seep showed a similarity of 0.999, indicating they also belong to the same species. However, the similarity between the GCA_016756955 and GCA_016756865 symbionts from the cold seep was only 0.93, suggesting they are distinct types of symbionts. Considering the ANI diversity among L. anaximandri symbionts, we speculate that adaptive evolution occurs as these symbionts colonize the tubeworms, enabling them to better adapt to their host environments. More detailed ANI information is available in the Supplementary material.

To examine the differences in ANI similarity and the evolutionary relationships of deep-sea tubeworm symbionts, this study constructed a phylogenetic tree of siboglinid symbionts, which formed five branches (Figure 3). In line with previous research (Yang et al., 2020), siboglinid symbionts within the Vestimentifera lineage exhibited habitat-related patterns, forming two distinct clades: one associated with hydrothermal vents and one associated with cold seeps. Notably, regardless of whether the L. anaximandri symbionts are found in hydrothermal or cold seep environments, they belong to the cold seep clade (Figure 3). This result suggests that for this class of tubeworms, which can inhabit hydrothermal vents and cold seeps, the relatives remain similar and are part of the same evolutionary lineage despite the presence of two species of symbionts.

The results of GTDB-TK classification aligned with the phylogenetic tree, dividing siboglinid symbionts into five groups at the genus level (Supplementary Table S3; Figure 3). The classification of siboglinid symbionts revealed significant diversity based on host species and habitat characteristics. The 26 siboglinid symbionts are distributed across three families: Sedimenticolaceae, Sulfurovaceae, and SZUA-299, with Sulfurovaceae being the second symbiont type found in L. satsuma. These symbionts are further classified into five genera: QGON01 (16 symbionts), Candidatus Endoriftia (five symbionts), Candidatus Vondammii (three symbionts), Sulfurovum (one symbiont), and SZUA-229 (one symbiont). The names QGON01 and SZUA-229 are provisional genus-level clusters derived from the GTDB genomic classifications (Parks et al., 2022). QGON01, Candidatus Endoriftia, and Sulfurovum are members of the Vestimentifera lineage, while Candidatus Vondammii is classified within the Sclerolinum clade, and SZUA-299 falls under the Frenulata group. The symbionts associated with the tubeworm L. anaximandri from both hydrothermal vent and cold seep environments belong to the QGON01 genus (Figure 3). In conclusion, the analysis suggests that the similarity among siboglinid symbionts is primarily driven by their habitat and host species. The findings suggest that host species exert a greater influence on symbiont evolution than habitat, highlighting the importance of host selection for symbiont survival.

To investigate the influence of varying environments and host selections on symbiont genomes, we assessed whether the four key metrics (genome size, GC content, gene density, and coding density) differed significantly among the five genera (Figure 4). We used the Kruskal-Wallis test to calculate the significance levels for genome size, GC content, gene density, and coding density of each genus (Theodorsson-Norheim, 1986). Notably, QGON01 symbionts from cold seeps had a significantly larger average genome size than Candidatus Endoriftia symbionts collected from hydrothermal vents (Kruskal-Wallis test, p < 0.01) (Figure 4A). Additionally, the GC content of Vestimentifera tubeworm symbiont genomes displayed an environment-dependent pattern. Candidatus Endoriftia symbionts collected from hydrothermal vents exhibited a significantly higher GC content compared to other genera (Kruskal-Wallis test, p < 0.01) (Figure 4B). The relationship between genome size and gene density is inverse in siboglinid symbionts, with larger genomes exhibiting smaller gene densities (Kruskal-Wallis test, p < 0.05) (Figure 4C). Widely distributed bacteria have large genomes that help them adapt to various environments (Zhou et al., 2022). The Vestimentifera tubeworms are found in cold seeps globally (Figure 1), which has led to the QGON01 symbiont exhibiting a large genome. Additionally, genome size is correlated with nutrient limitation (Giovannoni et al., 2014), which suggests that symbionts residing in nutrient-limited environments gain survival advantages by minimizing metabolic costs. Generally, the smallest and most compact genomes exhibit higher gene density due to a reduced allocation of genomic space to non-coding DNA (Giovannoni et al., 2005; Bobay and Ochman, 2017). Adaptation to extreme environments is a significant factor influencing variations in GC content. High GC content in Candidatus Endoriftia symbionts indicates their adaptation to varying environmental conditions caused by thermal gradients and oxidative stress (Teng et al., 2023), leading to higher energy consumption and enhanced thermal stability. In contrast, symbionts from other genera exhibit lower GC content, a consequence of their adaptation to cold seep environments (Teng et al., 2023). The growth rate of tubeworms in unstable hydrothermal vent environments is much higher than that of those in cold seep environments (Bright et al., 2013). This difference suggests that symbionts in hydrothermal vent environments must have a faster metabolic rate to provide sufficient nutrients to their hosts. The Candidatus Vondammii symbiont genome was found to have the highest gene density but the lowest coding density (Kruskal-Wallis, p < 0.05) (Figure 4D). This is due to the loss of many genes in the Candidatus Vondammii symbionts, leading to a simplified metabolic capacity (Gao et al., 2023).

We used the PPanGGOLiN software to conduct a pangenomic analysis and identify trends in the gene families of tubeworm symbionts related to persistent, shell structure, and cloud formation (Bazin et al., 2020; Gautreau et al., 2020). Heatmaps illustrating the distribution of gene families among tubeworm symbionts in the persistent, cloud, and shell categories reveal significant differences across various genera (see Supplementary Figure S4). Notably, the gene family distribution of L. anaximandri symbionts found in hydrothermal vents is consistent with that of L. anaximandri symbionts in cold seeps. This finding suggests that the tubeworm host has a more pronounced effect on the distribution of symbiont gene families than the habitat itself.

We predict evolutionary trends in siboglinid symbiont gene family partitioning by rarefaction curves (Supplementary Figures S5, S6). The “persistent” curve of the siboglinid symbionts plateaued, indicating that genes essential for key metabolic pathways and biosynthetic capabilities remained largely constant in the symbiont genome. The “shell” curve changed minimally, suggesting that genes acquired through horizontal gene transfer, such as those encoding functions related to environmental adaptation, pathogenicity, virulence, or secondary metabolites, exhibited minor variations. The “cloud” curve showed an upward trend, implying that genes acquired through horizontal gene transfer, such as phage-related genes, antibiotic resistance genes, and plasmid genes, increased continuously. We speculate that the trend of core metabolic genes is stable in siboglinid symbionts, indicating that the symbionts can steadily produce energy to meet the nutritional requirements of both themselves and their hosts. At the same time, the related secondary metabolites showed small fluctuations. The increasing trend of phage and other related genes suggests that phages have complex interactions with the tubeworm symbionts. Specifically, phages play an important role in the evolution of tubeworm symbionts by influencing their metabolism, environmental adaptation, and facilitating nutrient release (Han et al., 2024; Huang et al., 2024). For example, genes involved in carbon metabolism (korABCD, sucDC, sdhB, mdh, por, rbcL, PGK, GAPDH, TKT, and PRK), nitrogen metabolism (napAB, norBC), and sulfur metabolism (aprAB, soxYZB, and dsrAB) are primarily located in the persistent partition. Additionally, genes related to the biosynthesis of valine, leucine, and isoleucine are also present in the persistent partition. In the shell partition, genes involved in pantothenate and CoA biosynthesis, biotin metabolism, and terpenoid backbone biosynthesis are predominantly distributed. In the cloud partition, genes related to prokaryotic defense systems, the biosynthesis of ansamycins, the biosynthesis of enediyne antibiotics, carbapenem biosynthesis, and streptomycin biosynthesis are found. A more detailed distribution of these genes is provided in the Supplementary Table S4.

Siboglinid symbionts possess typical metabolic pathways for nutrient generation, including amino acid metabolism, cofactor/vitamin metabolism, and carbohydrate metabolism. These symbionts can produce 11 amino acids and 12 cofactors/vitamins (Supplementary Table S6).

Amino acids are essential for symbionts to synthesize nutrients and energy. Glycine not only supports growth but also plays a crucial role in various metabolic processes. However, high levels of glycine can lead to the production of toxic metabolites, such as aminoacetone and methylglyoxal (Kim et al., 2015). The glycine cleavage system regulates glycine catabolism, resulting in the production of CO₂, NH₃, NADH, and 5,10-methylene-THF (Kume et al., 1991; Nakai et al., 2005; Liu et al., 2021). In this study, key genes involved in the complete glycine cleavage system (gcv (PA, PB, and T) and DLD) were identified in 23 symbiont genomes, primarily in the SZUA-229, Candidatus Vondammii, Candidatus Endoriftia, and QGON01 genera (Supplementary Table S6). These four types of symbionts can utilize glycine as a carbon source for energy and biomass synthesis, as well as for detoxification.

Only the Sulfurovum symbiont, which contains the key genes involved in threonine biosynthesis (lysC, asd, hom, thrB, and thrC), can produce threonine, the crucial precursor required for isoleucine biosynthesis. Although the key genes involved in isoleucine biosynthesis (ilvABCDEH) are present in all symbionts, we hypothesize that these genes primarily facilitate isoleucine production through the valine/isoleucine biosynthesis (ilvBCDEH) pathway, rather than the isoleucine biosynthesis pathway. Only the Sulfurovum symbiont can use both pathways to synthesize isoleucine. Notably, Candidatus Vondammii symbionts can also convert pyruvate to 2-oxobutanoate to regenerate isoleucine via the isoleucine biosynthesis pathway (Supplementary Table S6). Additionally, the results reveal that siboglinid symbionts can synthesize lysine through two distinct pathways. First, Candidatus Vondammii and QGON01 symbionts synthesize lysine through the DAP aminotransferase pathway. Key genes involved in this pathway, including lysAC, asd, dapABF, and LL-diaminopimelate aminotransferase, were identified in 19 genomes. Alternatively, SZUA-229 and Candidatus Endoriftia symbionts synthesize lysine via the succinyl-DAP pathway, with key genes (lysCA, asd, dapABDCEF, and argD) identified in six genomes (Supplementary Table S6). This diversity aligns with the general strategy in which siboglinid symbionts utilize different pathways to synthesize the same biomolecule depending on environmental and physiological needs (details on the metabolism of other amino acids and cofactor/vitamin metabolism are provided in the Supplementary material). However, the symbionts encode few substrate-specific transporters, which results in an inability to effectively transport nutrients to the host (Yang et al., 2020). On the other hand, digestive enzymes are present in the host’s trophosome (Wippler et al., 2016), enabling it to obtain nutrients by digesting the symbionts. The ability of siboglinid symbionts to metabolize multiple amino acids and cofactors/vitamins reflects their capacity to provide nutrients to their hosts and sustain themselves autotrophically in extreme environments.

Despite the high oxygen buffering capacity of the tubeworm’s hemoglobin (Arp and Childress, 1983), elevated oxygen levels may still be transmitted to the symbiotic bacteria in the trophosome, potentially subjecting them to oxidative stress. We identified the alkyl hydroperoxide reductase ahpC, involved in oxidative stress response, in the symbionts of all 26 siboglinid species. Besides, siboglinid symbionts can supply reducing power through the oxidative branch of the pentose phosphate pathway, which helps maintain intracellular homeostasis via suppressing oxidative stress (TeSlaa et al., 2023). In this work, key genes (G6PD, PGLS, and PGD) involved in the oxidation stage of this pathway were detected in 21 symbiont genomes across the Candidatus Vondammii, Candidatus Endoriftia, and QGON01 genera. The Sulfurovum symbiont also possesses key genes (rpe, rpiB, tktA, ALDO, and GPI) involved in the non-oxidative stage of the pathway, suggesting that it can provide phosphorylated carbohydrates as precursors for ribose-5P synthesis (Miclet et al., 2001). Under nutrient-deficient conditions, siboglinid symbionts sustain growth through conserving energy. This study identified genes related to glycogen biosynthesis (UGP2, glgC, glgA, and GBE1) in 25 symbiont genomes in the Sulfurovum, Candidatus Vondammii, Candidatus Endoriftia, and QGON01 genera, reflecting their capacity to convert glucose 1P into glycogen for intracellular glucose homeostasis and a consistent energy supply (Moslemi et al., 2010).

Microbial terpenoids, primarily known as infochemicals, play a crucial role in intra- and inter-specific communication and interactions (Avalos et al., 2022). Infochemicals also aid microorganisms in adapting to oxidative, nitrosative, and nutrient stress (Avalos et al., 2022). Isoprene serves as the structural foundation for terpenoids (Bouvier et al., 2005). Among the 14 symbionts of QGON1, crucial genes related to C5 isoprenoid biosynthesis (dxsr, ispDEFH, gcpE, and idi) and C10–C20 isoprenoid biosynthesis (idi, GGPS, and ispA) were identified (Supplementary Tables S5, S6). Exposure to reactive oxygen species (ROS) induces oxidative stress in bacteria. Studies have shown that terpenoid compounds enhance the antioxidant potential of bacteria (Ostrovsky et al., 2003). Low temperatures reduce enzyme activity, and terpenes support bacterial growth under cold stress by modulating the fluidity of the cell membrane (Seel et al., 2018). Under nutrient limitation, some bacteria can utilize terpenoid compounds as a carbon source (Trutko et al., 2008). Therefore, we hypothesize that terpenoid compounds may enable the symbiont QGON01 to better adapt to the host’s internal environment and sustain its growth and reproduction (Bouvier et al., 2005).

This study identified genes encoding UDP-glucose (glk, pgm, pmm-pgm, and UGP2) and UDP-GlcNAc (glk, GPI, and glmSM, and U) in 21 genomes, primarily in the Candidatus Endoriftia and QGON01 genera. These symbionts inhabit the Vestimentifera tubeworm, whose trophosome is a highly vascularized organ (Hand, 1987; Uchida et al., 2023). UDP-glucose can differentiate into multiple blood cell lines (Kook et al., 2013), suggesting that the synthesized UDP-glucose may facilitate trophosome formation in Vestimentifera tubeworms, enabling rapid symbiont colonization within the tubeworms.