Metabolic handoffs between multiple symbionts may benefit the deep-sea bathymodioline mussels

Abstract

Bathymodioline mussels rely on thiotrophic and/or methanotrophic chemosynthetic symbionts for nutrition, yet, secondary heterotrophic symbionts are often present and play an unknown role in the fitness of the organism. The bathymodioline Idas mussels that thrive in gas seeps and on sunken wood in the Mediterranean Sea and the Atlantic Ocean, host at least six symbiont lineages that often co-occur. These lineages include the primary symbionts chemosynthetic methane- and sulfur-oxidizing gammaproteobacteria, and the secondary symbionts, Methylophagaceae, Nitrincolaceae and Flavobacteriaceae, whose physiology and metabolism are obscure. Little is known about if and how these symbionts interact or exchange metabolites. Here we curated metagenome-assembled genomes of Idas modiolaeformis symbionts and used genome-centered metatranscriptomics and metaproteomics to assess key symbiont functions. The Methylophagaceae symbiont is a methylotrophic autotroph, as it encoded and expressed the ribulose monophosphate and Calvin-Benson-Bassham cycle enzymes, particularly RuBisCO. The Nitrincolaceae ASP10-02a symbiont likely fuels its metabolism with nitrogen-rich macromolecules and may provide the holobiont with vitamin B12. The Urechidicola (Flavobacteriaceae) symbionts likely degrade glycans and may remove NO. Our findings indicate that these flexible associations allow for expanding the range of substrates and environmental niches, via new metabolic functions and handoffs.

Introduction

Chemosynthetic symbioses allow animals and some protists to colonize extreme environments including hydrothermal vents and cold seeps in the deep sea, as well as key productive habitats in coastal areas, such as the seagrass meadows [1, 2]. Chemosynthesis, that is, the assimilation of single-carbon molecules, such as carbon dioxide and methane, using the energy stored in reduced compounds, primarily reduced sulfur (sulfide, thiosulfate, elemental sulfur) and methane, fuels these symbioses. Chemosynthetic symbioses are usually characterized by low diversity and high fidelity of key nutritional symbionts at the species level [3]. They exhibit a broad range of transmission strategies, including vertical, horizontal and mixed modes [4], that may determine the fidelity of these associations, as well as the strength of environmental selection for the fittest symbionts and genetic diversity within the symbiont populations [3, 5–7]. As in free-living bacterial populations, genetic diversity broadens the functions of symbionts, allowing the holobiont to acquire new metabolic capabilities and environmental niches [8, 9].

Symbiont taxonomic and functional diversity varies markedly between lineages of bathymodioline mussels, which include the chemosynthetic genera Bathymodiolus, Gigantidas and Idas among others [10]. The key symbionts of bathymodioline mussels include the sulfur-oxidizing Thioglobaceae (gammaproteobacterial order PS1, also known as the SUP05 clade), and the Methyloprofundus (gammaproteobacterial order Methylococcales) methane-oxidizing bacteria. Some bathymodioline mussels, such as Bathymodiolus earlougheri, B. billschneideri and B. nancyschneideri host single-species populations of thiotrophic symbionts [11], whereas Gigantidas childressi and G. platifrons host mainly the methane-oxidizing bacteria [12–14]. Others host both [15]. There are also extreme cases of symbiont diversity loss and gain in bathymodioline symbioses. For example, a small deep-sea bathymodioline mussel Idas argenteus appears to lack symbionts, which may have been gained in the past and recently lost [16]. Others, such as B. heckerae, from the deep Gulf of Mexico, host several taxa of symbionts with distinct functions: apart from two species of methane-oxidizing symbionts and two species of sulfur-oxidizing symbionts, B. heckerae has additional bacterial symbionts, including the Methylophagaceae sp. methylotrophs, as well as Cycloclasticus that catabolize short-chain alkanes [17–19].

Here we focus on another extreme case of multi-species symbioses in bathymodioline mussels, namely Idas mussels from the Mediterranean Sea. Several previous studies noted the co-occurrence of at least six gill symbiont phylotypes in an Idas species from the east Atlantic and the Mediterranean Sea; this Idas species was first identified as Idas sp. MED, and is now called Idas modiolaeformis [20–23]. These small mussels thrive in seeps and wood falls and the diversity of their symbionts appears to be linked to the availability and composition of reduced fluids [20]. Early studies identified the following six key phylotypes: the Methyloprofundus sedimenti-related type I methanotrophs (M1), two phylotypes of Thioglobaceae sulfur-oxidizing symbionts (S1 and S2), Methylophagaceae methylotrophs (M2), as well as two phylotypes that lacked the potential for being chemosynthetic – Bacteroidetes (CFB) and a gammaproteobacterial lineage (G) [21]. Hereafter we refer to the methanotrophs (M1) and thiotrophs (S1 and S2) as the primary symbionts, and the others as secondary symbionts (M2, CFB, and G). These phylotypes were found to be associated with bacteriocytes in the symbiont-bearing gill tissue using fluorescence in situ hybridization (FISH), confirming the symbiotic nature of their association with I. modiolaeformis [21, 22]; however, these observations did not provide conclusive evidence for the intracellular localization of all the symbionts [22]. Other variants of the gammaproteobacterial secondary symbionts, were discovered later [20, 23]. While the functions of the primary symbionts and the methylotroph M2 were estimated based on the detection of marker genes, the functionality of the symbiotic methylotrophs has not been studied in detail using omics, and little is known about the role of other symbiont lineages.

We aimed to characterize the metabolism and physiology of the Idas symbionts, using genome-centric metagenomics, metatranscriptomics and metaproteomics. We collected 14 Idas individuals from hydrocarbon seeps and brine pools off the shore of Israel [24, 25]. We hypothesized that the secondary symbionts may contribute to holobiont fitness by providing important functions via metabolic handoffs. Alternatively, these potential heterotrophs, which are most abundant in mussels that colonize plant debris, such as sunken wood, may extract nutrients from organic substrates [20, 23]. We thus investigated the metabolic potential of the Idas symbionts, providing a snapshot of host-symbiont-symbiont interactions in this specific multi-member symbiosis.

Materials and methods

Sample collection, processing, and sequencing

Hydrocarbon seeps are found at the toe of Palmahim Disturbance, a large deformation feature at the Israeli margin [26, 27] (Supplementary Fig. S1). In 2011 we collected one Idas specimen from a carbonate sample that was retrieved, from the Palmahim pockmark 32°09.6′N and 34°10.0′E at a depth of 1036 m, using Hercules remotely operated vehicle (ROV), based on the Nautilus E/V. In 2021, multiple individuals (n = 14 studied herein) were collected from the edges of a brine pool at the toe of Palmahim disturbance, at 1,150 m water depth (32° 13.4’ N 34° 10.7’ E) [25], using the work-class SAAB Seaeye Leopard ROV, based on R/V Bat Galim (Supplementary Fig. S1).

Metagenomics

We preserved a single individual collected in 2011, as well as 5 individuals that were collected in 2021, at −20 °C. DNA was extracted from these individuals using the DNeasy Blood & Tissue Kit (Qiagen), following the manufacturer’s instructions. DNA libraries were constructed for these six individuals at HyLabs, Israel, and sequenced using 120 million 2 × 150 bp paired-end reads per sample with an Illumina NovaSeq at GENEWIZ, Germany.

Metatranscriptomics

An additional 4 individuals were preserved in RNAlater. DNA, RNA and proteins were extracted from these individuals using the AllPrep DNA/RNA/Protein Mini Kit (Cat. No. 80004, Qiagen). RNA libraries were constructed for these individuals at Novogene, Singapore, and sequenced using 100 million 2 × 150 bp paired-end reads per sample with an Illumina NovaSeq, following library construction with the NEBNext Ultra RNA Library Prep Kit for Illumina (Cat No. 7530) and ribosomal RNA depletion with the Ribo-Zero Plus rRNA Depletion Kit (Bacteria) (Cat No. 20037135) & Ribo-Zero Magnetic Kit (Plant Leaf).

Protein extraction and peptide preparation for metaproteomics

We resuspended the protein precipitates from the 4 individuals extracted with the AllPrep DNA/RNA/Protein Mini Kit kit in 60 µl of SDT lysis buffer [4% (w/v) SDS, 100 mM Tris-HCl pH 7.6, 0.1 M DTT] and heated to 95 °C for 10 min. The SDT protein mixture was cleaned up, reduced, alkylated and digested using the filter-aided sample preparation (FASP) protocol as described previously [28]. We performed all centrifugation steps mentioned below at 14,000 x g. We combined lysates (60 μl) with 400 μl of urea solution (8 M urea in 0.1 M Tris-HCl pH 8.5) and loaded it onto 10 kDa MWCO 500 μl centrifugal filters (VWR International) followed by centrifugation for 20 min. We washed filters once by applying 200 μl of urea solution followed by 20 min of centrifugation to remove any remaining SDS. We performed protein alkylation by adding 100 μl IAA solution (0.05 M iodoacetamide in urea solution) to each filter and incubating for 20 min at room temperature followed by centrifugation for 20 min. The filters were washed three times with 100 µL of urea solution with 15 min centrifugations, followed by a buffer exchange to ABC (50 mM ammonium bicarbonate). Buffer exchange was accomplished by three cycles of adding 100 μl of ABC buffer and centrifuging for 15 min. For tryptic digestion, we added 1 μg of MS grade trypsin (ThermoFisher Scientific) in 40 μl of ABC buffer to each filter and incubated for 16 h in a wet chamber at 37 °C. We eluted tryptic peptides by adding 50 μl 0.5 M NaCl and centrifuging for 20 min. Peptide concentrations were determined with the Pierce Micro BCA assay (ThermoFisher Scientific) following the manufacturer’s instructions.

LC-MS/MS

All proteomic samples were analyzed by 1D-LC-MS/MS as previously described [28]. We loaded 1.2 μg of peptide from each sample onto a 5 mm, 300 µm ID C18 Acclaim® PepMap100 pre-column (Thermo Fisher Scientific) with an UltiMateTM 3000 RSLCnano Liquid Chromatograph (Thermo Fisher Scientific) in loading solvent A (2% acetonitrile, 0.05% trifluoroacetic acid). Elution and separation of peptides on the analytical column (75 cm × 75 µm EASY-Spray column packed with PepMap RSLC C18, 2 µm material, Thermo Fisher Scientific; heated to 60 °C) was achieved at a flow rate of 300 ml min⁻¹ using a 140 min gradient going from 95% buffer A (0.1% formic acid) and 5% buffer B (0.1% formic acid, 80% acetonitrile) to 31% buffer B in 102 min, then to 50% B in 18 min, and finally to 99% B in 1 min and ending with 99% B. The analytical column was connected to a Q Exactive HF hybrid quadrupole-Orbitrap mass spectrometer (Thermo Fisher Scientific) via an Easy-Spray source. Eluting peptides were ionized via electrospray ionization (ESI). Carryover was reduced by a wash run (injection of 20 µl acetonitrile, 99% eluent B) between samples. MS1 spectra were acquired by performing a full MS scan at a resolution of 60,000 on a 380 to 1600 m/z window. MS2 spectra were acquired using data-dependent acquisition, selecting for fragmentation the 15 most abundant peptide ions (Top15) from the precursor MS1 spectra. A normalized collision energy of 25 was applied in the HCD cell to generate the peptide fragments for MS2 spectra. Other settings of the data-dependent acquisition included: a maximum injection time of 100 ms, a dynamic exclusion of 25 s, and the exclusion of ions of +1 charge state from fragmentation. About 50,000 MS/MS spectra were acquired per sample.

Protein identification and quantification

We constructed a protein sequence database for protein identification using the protein sequences predicted from the metagenome-assembled genomes obtained in this study. To identify peptides from the host, we used the annotated protein sequences of the bathymodioline mussel Bathymodiolus childressi obtained in a previous study (PRIDE PXD008089-1) [29], since no annotated Idas protein sequences were available. We added sequences of common laboratory contaminants by appending the cRAP protein sequence database (http://www.thegpm.org/crap/). The final database contained 49,401 protein sequences and is included in the PRIDE submission (see data access statement) in fasta format. Searches of the MS/MS spectra against this database were performed with the Sequest HT node in Proteome Discoverer 2.3.0.523 as previously described [30]. The peptide false discovery rate (FDR) was calculated using the Percolator node in Proteome Discoverer and only peptides identified at a 5% FDR were retained for protein identification. Proteins were inferred from peptide identifications using the Protein-FDR Validator node in Proteome Discoverer with a target FDR of 5%. To estimate species abundances based on proteinaceous biomass using the metaproteomic data we followed the previously-described approach [31] with the added filter criterion of requiring two protein-unique peptides for a protein to be retained for biomass calculations.

Bioinformatics

Metagenomes were assembled using SPAdes V3.15 with –meta k = 21,33,66,99,127 parameters [32], following adapter trimming and error correction with tadpole.sh, using the BBtools suite following read preparation with the BBtools suite (Bushnell, B, sourceforge.net/projects/bbmap/). Downstream mapping and binning of metagenome-assembled genomes (MAGs) were performed using DAStool, Vamb, Maxbin 2.0 and Metabat2 [33–36] within the Atlas V2.9.1 framework [37], using the genome dereplication nucleotide identity threshold of 0.975. MAG quality was verified using Checkm2 [38] and QUAST [39]. Taxonomy was assigned using GTDB-tk V2.1 [40]. Functional annotation was carried out using the Rapid Annotations using the Subsystems Technology RAST server [41], and key annotations were verified by BLASTing against the NCBI database. RNA reads were quality-trimmed and mapped to the MAGs using BBmap (Bushnell, B, sourceforge.net/projects/bbmap/), and read counts were assigned to coding sequences using FeatureCounts [42]. The counts were normalized as transcripts per million (TPM), separately for each MAG [43].

Phylogenetic analyses were performed with IQ-TREE 2 [44] or MEGA11 [45], based on the best model determined by ModelFinder, following sequence alignment with MAFFT [46]. For phylogenomics, we used GTOtree [47] and Fasttree [48], as implemented in GTOtree. Mitochondrial genomes were assembled with MitoZ [49], rotated with MARS [50] and aligned with MAFFT [46].

Fluorescent in-situ hybridization

Four additional individuals were fixed in 2% paraformaldehyde in 1x phosphate-buffered saline (PBS) for at most 12 h at 4 °C, rinsed three times in 1x PBS, and stored at 4 °C in 0.5x PBS/50% ethanol. Whole tissue from Idas individuals was dehydrated, embedded in Steedman’s wax, and cut into 8 μm sections as previously described [19]. To identify the secondary symbionts, we used the Cy3-labeled BhecM2-822, ImedaG-193 and CF319 probes, targeting the Methylophagaceae, Nitrincolaceae and Flavobacteriaceae symbionts, respectively [22]. Cy5-labeled EUB338 targeting most bacteria was used as a positive control [51] and the Cy3-labeled NON338 probe was used as a control for background autofluorescence [52]. Photomicrographs were acquired with Nikon A1R Confocal Laser Scanning Microscope, using the Plan Apochromat VC 61.5x objective with oil immersion. The brightness and contrast of the images were adjusted with Adobe Photoshop (Adobe Systems, Inc., USA).

Results and discussion

Idas modiolaeformis from the eastern Mediterranean hosts at least six symbiont genotypes

The mussels were identified as Idas modiolaeformis, based on the analysis of mitochondrial cytochrome c oxidase I (MT-CO1) sequences from the six metagenomes, using the Barcode of Life Data System identifier. These sequences were 98.4–98.6% similar to those of the I. modiolaeformis in the database (for example, the eastern Atlantic individuals, KT216487 in NCBI). The best hits, 98.7–98.9% similarity, were to the sequences of Idas individuals recovered from the Nile Deep Sea Fan mud volcano at the depth of 1693 meters (FM212787). Hereafter we refer to our samples as I. modiolaeformis, as among other known Idas species, the next closest hit was Idas macdonaldi with a lower MT-CO1 similarity of 95%. Yet, given that the individuals from the Plamahim Disturbance and the Nile Deep Sea Fan are not genetically identical, it is plausible that genetic exchange between these populations is low. We observed some variation of MT-CO1 and full mitochondrial sequences among the Palmahim individuals corresponding to the occurrence of different haplotypes. Genetic distance may increase with the geographical distance or time, as the individual that we collected in 2011 from a distinct seep located several hundreds of meters from the brine pool site where other individuals were obtained, was the most diverged (Supplementary Fig. S2). Yet, as we discuss below, the host phylotypes were not linked to the symbiont diversity, in agreement with previous findings [20].

We generated metagenome-assembled genomes (MAGs) for the six key I. modiolaeformis symbiont phylotypes. We note that besides these six lineages, metagenomic binning resulted in the discovery of MAGs for two additional taxa – a Rhizobiales alphaproteobacterium, as well as a Bdellovibrionaceae species. Given that the nature of the association between these species and I. modiolaeformis is unclear, as they have not been previously described by amplicon sequence, we hereafter focus only on the six previously detected lineages. For these, good or high-quality genomes were assembled and binned (Table 1).

In the previous amplicon-based sequencing studies, the symbionts, particularly the secondary ones, were only classified to class and phylum taxonomic levels. We used our MAGs to obtain improved taxonomic assignments through the Genome Taxonomy Database (GTDB) inference tool and validated their taxonomic placement using phylogenies based on the reference collections of closely related genomes from GenBank. We confirmed that the methane-oxidizing symbionts belonged to the Methyloprofundus clade, also known as marine methylotrophic group 1 or MMG1, and verified their placement using a PmoA protein phylogeny, as genomes of most Methyloprofundus symbionts are missing from the databases (Supplementary Fig. S3). The sulfur oxidizers belonged to two distinct clades, Candidatus Thioglobus A and Candidatus Thiodubilierella (GTDB inference, hereafter Thioglobus and Thiodubilierella). The Thiodubilierella genus name has been suggested for the symbionts of B. septemdierum [53], and was integrated into the GTDB taxonomy. Yet, Thiomultimodus has been proposed as an alternative name for both the Thioglobus and Thiodubilierella clades, with two SUP05 A and B subclades [54]. Whereas Thioglobus (SUP05 clade A) is often represented by free-living species and clam symbionts, Thiodubilierella (SUP05 clade B) has been described only in mussel symbioses [54]. Hereafter we refer only to Thioglobus and Thiodubilierella nomenclature for the two thiotrophic symbiont clades. The co-occurrence of these clades in an individual has only been documented in B. heckerae [17–19] and I. modiolaeformis [22].

The new MAG-based taxonomic assignments greatly improved the classification of the secondary symbionts and allowed us to make predictions about their physiology based on previously characterized relatives. The methylotrophic symbionts were assigned to the Methylophagaceae clade GCA-002733105. Phylogenetic placement based on a tree using single-copy genes indicated that this clade lacks a cultivated representative, and includes lineages commonly found at cold seeps (e.g., the genome of the most closely related lineage, Methylophaga GLR1851, was curated from a sample at the Hikurangi Margin gas hydrate deposits), as well as the only other Methylophagaceae symbionts that occur in B. heckerae (Supplementary Fig. S4).

Phylotype G belonged to the Nitrincolaceae family ASP10-02a (Table 1, Supplementary Fig. S5). Nitrincolaceae are rarely found in symbioses; however, the Rs1 and Rs2 symbionts of the bone-eating worm Osedax belong to the Nitrincolaceae family (unnamed genus Rs1) [55–57]. Nitrincolaceae are prominent degraders of nitrogen-rich compounds such as amino acids [55, 58, 59], and the ASP10-02a clade is prominent in the water column, often exhibiting high abundance and activity during algal blooms [60, 61]. The symbionts Idas and Osedax belong to distinct clades within Nitrincolaceae, their 16 S rRNA gene sequences are only ~90% similar and their average nucleotide identity (ANI) is <70%. Given that the hosts thrive on different substrates, and that the symbionts occupy different tissues, we hypothesize that the functionalities of these symbionts may differ.

The Bacteroidota (CFB) phylotype belonged to the genus Urechidicola (Flavobacteriaceae) (Table 1, Supplementary Fig. S6). The first cultivated Urechidicola genus representative, U. process, was isolated from the intestine of a marine spoonworm, Urechis unicinctus [62]. Similar to other Flavobacteriaceae, such as the closely related Lutibacter, Urechidicola can degrade multiple organic compounds, including DNA and starch [62, 63]. The closest relative of the Urechidicola symbiont in Idas was found in a bone-degrading biofilm, and it is represented by GenBank assembly accession GCA_016744415.1. Taxonomic affiliation suggests that both of these secondary symbionts are likely copiotrophs that can catabolize both difficult-to-degrade and protein-rich organic matter, such as polysaccharides, nucleic acids and peptides, in the marine environment.

The metagenomic read abundances of the primary symbionts varied among the six Idas individuals that were analyzed with metagenomics (Fig. 1). A single individual (Idas 16) lacked the methane-oxidizing and methylotrophic symbionts. In this individual, the Nitrincolaceae symbiont was the most abundant, and Urechidicola was found (metatranscriptomics and metaproteomics were, unfortunately, not performed for this individual). We performed metatranscriptomics and metaproteomics on four individuals, which were distinct from the ones used for metagenomics. In these four individuals, the relative abundances of symbionts were consistent between metatranscriptomics and metaproteomics measurements, showing that the sulfur oxidizers were most abundant in sample 3 (5% and 2% biomass for Thioglobus and Thiodubilierella, respectively, Fig. 1B, C). It is important to note that the relative abundances obtained from metagenomics cannot be compared directly with the relative abundances obtained with metaproteomics; metagenomics-derived abundances are a good estimator of cell numbers, whereas metaproteomics-derived abundances are a good estimator of species biomass in a microbial population [31]. The methane-oxidizing bacteria were the most abundant in the four samples analyzed with metatranscriptomics (~70–80%) and metaproteomics (~89–98% estimated biomass). The biomasses of the other symbionts were much lower, in the range of 1.6 to 7.2% for both thiotrophs, with only very few peptides detected for Urechidicola. In agreement with the previous study [20], our data suggest that the symbiont proportions differ between individuals, likely based on local environmental conditions experienced by individuals. Methylophagaceae symbionts may require methane-oxidizers to be present to supply their substrate methanol. Nitrincolaceae are consistently found among the samples with meta-omics, and are metabolically active, as their genes and proteins were expressed. Methylophagaceae, Nitrincolaceae and Urechidicola symbionts appear to co-occur with the primary symbionts in the gill tissue where the primary symbionts were present (Fig. 2), in agreement with previous data [22]. Yet, it is still not clear if the secondary symbionts are extra- or intra-cellular.

The relative abundances and biomass distributions of the six key symbiont genotypes in individual Idas modiolaeformis specimens.

Fluorescent in situ hybridization (FISH) shows the presence of the secondary symbionts in the Idas gill tissue.

Oxidation of methane and its products fuel the Idas symbiosis

Methane oxidizers are the key nutritional symbionts of I. modiolaeformis, as suggested by their abundance in the metatranscriptomes and metaproteomes. The particulate methane monooxygenase PmoCAB complex was most well-expressed at both mRNA and protein levels (Supplementary Fig. S7, here and hereafter see Supplementary Table S1 for the full list of genes and proteins expressed). The lanthanide-containing XoxF type methanol dehydrogenase was the predominant one; the calcium-dependent enzyme MxaFI was also found but expressed at lower levels. Carbon assimilation and energy conservation from formaldehyde oxidation function via the ribulose monophosphate (RuMP) pathway, using only the Entner–Doudoroff (ED) bypass, but not the Embden–Meyerhof–Parnas (EMP) variant [64], as the pfk gene encoding the pyrophosphate-dependent phosphofructokinase was not found in the genomes of the Methyloprofundus symbiont of I. modiolaeformis (Supplementary Fig. S7). This is similar to Bathymodiolus japonicus symbionts and in contrast to those of Bathymodiolus platifrons, which have both the ED and EMP variants [65]. The serine cycle is incomplete, as hydroxypyruvate reductase was missing; however, some central proteins of this pathway, such as the serine-glyoxylate aminotransferase and malyl-CoA lyase, were substantially expressed at both the RNA and amino acid levels, indicating that the partial serine pathway is functional (Supplementary Fig. S7) [15]. All the genes in the tricarboxylic acid (TCA, Krebs) cycle were found and expressed, therefore energy can be conserved through the oxidation of compounds derived from methane assimilation, such as pyruvate (Supplementary Fig. S7). Formaldehyde detoxification and oxidation to CO₂ via the dissimilatory tetrahydromethanopterin route is likely prominent, due to the high expression of this pathway, in particular, that of the 5,6,7,8-tetrahydromethanopterin hydro-lyase (fae gene, Supplementary Fig. S7). Both oxygen and nitrate respiration is feasible, given the presence of the respiratory NarGHIJ and NirK, yet the respiratory complex IV was expressed at much higher levels, indicating that oxygen is the key electron acceptor for energy conservation (Supplementary Fig. S7).

The metabolism of the Methylophagaceae symbiont is a rare case of methylotrophic autotrophy in chemosynthetic symbioses. Their co-occurrence with the methane oxidizers suggests that they might use methanol, formaldehyde and other metabolites produced by the methane oxidizers. Methane oxidizers excrete metabolites such as methanol, formaldehyde, formate, acetate, and succinate [66–68], which can be used by co-occurring organisms, in particular, methylotrophs [69]. Partnerships between methane oxidizers and methylotrophs are widespread in free-living communities [69, 70], and our results suggest that this association plays a role in Idas symbiosis.

The Methylophagaceae symbionts lack methane monooxygenases, yet highly express the pyrroloquinoline quinone (PQQ)-dependent methanol dehydrogenase, as well as the key enzymes of the RuMP pathway for formaldehyde assimilation (Fig. 3). As opposed to the primary methane-oxidizing symbionts, the methylotrophs encoded and expressed both the ATP-dependent EMP and ED variants of the RuMP pathway, which have distinct energetic demands [64]. Formate oxidation is also more flexible in the methylotrophs, as we identified the occurrence and expression of not only the two-subunit tungsten-containing formate dehydrogenase FdhAB, but also the respiratory molybdoenzyme dehydrogenase-O (FdoGHI) that catalyzes the oxidation of formate to carbon dioxide, donating electrons to the membrane soluble quinone pool [71–74]. This hints that the oxidation of formate alone can fuel the metabolism of the Methylophagaceae symbiont. In addition to originating from methane oxidation, the used formate could also come from mitochondrial metabolism [75]. In terms of terminal electron acceptors (TEA) for methanol and formate oxidation, we only found genes for the use of oxygen as TEA and no genes for other TEAs such as nitrate.

Central carbon metabolism in the Methylophagaceae symbiont of Idas modiolaeformis.

Alongside the RuMP pathway, Methylophagaceae symbionts have the genes for all steps of the Calvin–Benson–Bassham (CBB) cycle and highly express the key enzyme, form II ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO, Fig. 3). We note that although form I RuBisCo appeared to be widespread in Methylophagaceae, form II RuBisCo was rare in this clade (Supplementary Fig. S8). These closely related RuBisCo forms differ in their affinity to oxygen and CO₂, and form II enzymes have a low specificity factor, that is, function better under high CO₂ and low O₂ conditions [76–78]. Such conditions may exist in the host’s bacteriocytes, giving an advantage to bacteria with form II RuBisCO [79].

Similar to obligate autotrophs [80, 81], the TCA cycle appears to be incomplete in most cultivated Methylophagaceae, which lack the 2-oxoglutarate dehydrogenase activity [82–84]. We did not find the 2-oxoglutarate dehydrogenase-encoding genes in the high-quality MAG of the symbiotic Methylophagaceae, indicating that the TCA cycle was incomplete. The incomplete TCA allows for the production of intermediates from one-carbon (C1) compounds and prevents futile cycling, that is, the destruction of larger organic compounds via catabolism. This is typical of obligate methylotrophs, such as Methylobacillus flagellates [85]. In these bacteria, energy and reducing equivalents may largely result from the oxidation of formaldehyde to CO₂ [85] (Fig. 3).

Variability in terminal electron acceptors, metabolite uptake and vitamin usage may lead to niche differentiation among the thiotrophic symbionts

The two sulfur-oxidizing symbionts, Thioglobus and Thiodubilierella species, appear to be obligate autotrophs as the key enzymes of the TCA cycle, 2-oxoglutarate dehydrogenase and malate dehydrogenase, were not found in either MAG (Supplementary Fig. S9). This is in agreement with previous observations of symbiotic and free-living organisms from this clade [54, 86, 87]. Both sulfur-oxidizing symbionts have the genes to gain energy from sulfur compounds using the sulfide: quinone oxidoreductase (Sqr), the Sox sulfur oxidation system (SoxYZXAB) and the reverse dissimilatory sulfite reduction (rDSR) system, which catalyzes the oxidation of sulfide to sulfate via the adenosine-5’-phosphosulfate, using the adenylylsulfate reductase and sulfate adenylyltransferase for the two final steps of the pathway. The genes encoding these functions were among the most abundantly transcribed ones in the metatranscriptomes and were often detected at the protein level (Supplementary Fig. S9). Both lineages highly expressed the Calvin-Benson-Bassham cycle enzymes, in particular the two subunits of the type I RuBisCo (in particular, the rbclLS genes and respective proteins). These obligate autotrophs likely can import some organic compounds similar to the symbionts of Bathymodiolus azoricus [15], as TRAP-type dicarboxylate transport systems were encoded and expressed (C4-dicarboxylate and unknown substrate 6 in Thioglobus, only the latter in Thiodubilierella).

Several differences between Thioglobus and Thiodubilierella were observed in their terminal electron acceptor use. Only Thiodubilierella carried the periplasmic NapABGH for the use of nitrate as a terminal electron acceptor, confirming the modularity of nitrogen metabolism in the Thioglobaceae clade [54]. We also found some discrepancies in oxygen respiration were found: Thioglobus encoded and highly expressed only the ccoNOP genes encoding the cbb₃-type cytochrome c oxidase, whereas Thiodubilierella encoded three variants (cbb₃, bo₃, ba₃). These terminal oxidases differ in their affinity to O₂ and H₂S [88], likely highlighting adaptation to distinct redox conditions.

Secondary symbionts are heterotrophs that can use nitrogen-rich and difficult-to-degrade metabolites

Data from all three meta-omics approaches suggested that the Nitrincolaceae symbiont is a heterotroph that is capable of using numerous substrates, given the presence and substantial expression of the complete TCA cycle, as well as that of multiple transport systems for organic compounds, such as peptides, amino acids and nucleosides, in line with the metabolic reconstruction of Osedax Nitrincolaceae symbionts [55, 59]. Relevant examples include: (1) most enzymes for the degradation of branched-chain amino acids co-occurred in the genomes and were transcribed (Fig. 4); (2) a cluster 11 RfuABCD riboflavin/purine/nucleoside transporter, which clustered with nucleoside degradation enzymes, such as deaminases and phosphorylases (Cda, Ada, DeoABC); (3) the putrescine transport (PotABCD) and degradation system (PuuABCDE) (although the puuC gene was not found, likely due to an assembly issue, as these genes occurred at the edge of a contig). An additional gene cluster encoded the complete PaaA-K and HpaA-I machinery for the catabolism of phenylacetate and hydroxyphenylacetate[89]. This pathway allows for the degradation of recalcitrant, plant-derived organic aromatics, especially in marine bacteria that thrive under fluctuating oxygen conditions in the Mediterranean Sea [90, 91]. Some catabolic reactions may be carried out anaerobically, particularly given the fermentative potential of this bacterium, suggested by the presence of lactate dehydrogenases (some expressed at the mRNA level, Fig. 4). We also identified the genes for nitrate respiration to nitric oxide including the genes for the respiratory nitrate reductase NarGHIJ and copper-dependant nitrite reductase NirK. The nirK gene was among the most highly expressed genes at the mRNA level in this symbiont (Fig. 4). However, oxygen is likely the key electron acceptor for these symbionts as the cbb₃- and caa₃ complex IV-encoding genes were highly expressed at the mRNA level (Fig. 4).

The Nitrincolaceae symbiont can import and catabolize multiple metabolites, including peptides, amino acids (in particular, branched-chain amino acids), nucleosides and some aromatics, such as phenylacetate.

A key feature of the metabolism encoded in the Urechidicola symbiont genome (Flavobacteriaceae, Bacteroidota) is its potential to degrade glycans (polysaccharides). Bacteroidota, Flavobacteriaceae in particular, are ubiquitous degraders of complex glycans and are found in many environments from human gut mucus to algal blooms in the ocean [92–95]. The diversity of glycans is very large and glycan degradation is carried out by an array of carbohydrate-active enzymes (CAZymes) that are often organized in polysaccharide utilization loci (PULs), which are typical of Bacteroidota [92]. The key features of PULs are the presence of SusCD transporters: the cell surface glycan-binding lipoprotein SusD and the outer membrane TonB-dependent transporter SusC, which can be found in multiple copies in Bacteroidota genomes [96, 97]. We identified four SusCD pairs in the Urechidicola symbiont, two of which were among the top 50 most abundantly transcribed genes in the metatranscriptomes, and one SusC protein was detected in the metaproteomes despite the low overall protein identification rate for this symbiont. We note that genomic coding sequence coverage was very low for Urechidicola in the metatranscriptome (17%) and metaproteome (5%, mostly with low confidence), and thus only the most active features were detected by these analyses. Given the considerable fragmentation of this genome (300 scaffolds), we couldn’t identify large PULs, yet multiple CAZymes were found and often clustered (Fig. 5). For example, a single 30,996 bp contig contained two SusCD pairs, as well as multiple CAZymes, such as DD-carboxypeptidase EC 3.4.16.4; L-Ala-D/L-Glu epimerase EC 5.1.1.20; glucosamine-6-phosphate deaminase EC 3.5.99.6; N-acetylmuramic acid 6-phosphate etherase EC 4.2.1.126; as well as an AmpG family muropeptide MFS transporter; and glycosyl hydrolases of family 18, family 10, and family 3 (Fig. 5). Thus, genomic and transcriptomic data provided evidence for the glycan degrading activity of Urechidicola, yet whether the substrates are derived from the gill (e.g., mucus), or the environment remains unknown. To date, FISH analyses do not show clearly if these bacteria are endo- or ectosymbionts of Idas gills. One piece of evidence suggesting that they might be ectosymbionts is that cell adhesion appears to play an important role for Urechidicola, as a sequence encoding the FAS1 (fascilicin) domain protein was among the top 20 most abundantly transcribed genes [98, 99]. One hypothesis that remains to be explored is whether there are metabolic handoffs between Urechidicola and Nitrincolaceae symbionts, for example, via the exchange of aromatic derivatives of polymer degradation, such as phenylacetate and hydroxyphenylacetate [100].

Polysaccharide utilization loci (PULs) in Urechidicola symbionts of Idas modiolaeformis.

Nitrincolaceae symbionts may contribute to holobiont fitness by producing vitamin B12

The Nitrincolaceae symbiont can produce adenosylcobalamin (vitamin B12), given the presence of almost all the genes needed for its synthesis (Fig. 4), similar to the Nintrincolaceae Osedax symbionts [55, 59]. The cobG gene was the only gene, out of 25 genes for vitamin B12 biosynthesis, that was undetected, likely due to the incompleteness of the MAG. Many of the vitamin B12 biosynthesis genes were transcribed and detected in the metatranscriptomes. We also found the genes for vitamin B12 biosynthesis in the methane-oxidizing symbiont, which is a likely key producer of adenosylcobalamin in I. modiolaeformis, given its dominant biomass in most individuals (Fig. 1). In hosts that lack the methane-oxidizing symbionts, B12 production by the alternative Nitrincolaceae symbionts may be crucial for the holobiont. Nitrincolaceae, specifically the ASP10-02a lineage to which the symbionts belong, are (i) dominant species in Earth’s cold oceans, (ii) are often associated with algal bloom degradation, (iii) appear to be the most prominent B12 producers in these habitats and (iv) supply B12 specifically to Methylophagaceae, suggesting that B12-based ASP10-02a-methylotroph/autotroph associations can be advantageous under some conditions [61, 101]. The B12-based dynamics play a key role in the water column [102], but also human gut [103] and insect symbioses [104], and thus also may be crucial for chemosynthetic symbioses. The outer membrane receptor BtuB of the vitamin B12 transporter was found in the genomes Methyloprofundus, Thioglobus, Urechidicola and Methylophagaceae symbionts. Whereas methane-oxidizing symbionts are likely prototrophs for cobalamin, but still can take it up from the environment, the other symbionts are auxotrophs and may depend on its uptake.

In bacteria, the key reaction that may depend on cobalamin is methionine synthesis, as cobalamin is required by methionine synthase; thus, the lack of this vitamin may hinder DNA synthesis, methionine regeneration and lead to homocysteine accumulation [105]. Thiodubilierella appears to lack the btuB gene needed for B12 import; however, unlike the Thioglobus symbiont, which encodes only the cobalamin-dependent methionine synthase (MetH; EC 2.1.1.13), Thiodubilierella has both the MetH, and the cobalamin-independent methionine synthase (MetE; 5-methyltetrahydropteroyltriglutamate–homocysteine methyltransferase; EC 2.1.1.14), and therefore may not depend on cobalamin for production of the essential amino acid methionine and thus growth [104]. We thus hypothesize, that similar to insect symbioses [104], the interplay between cobalamin requirement, synthesis and uptake may contribute to determining the complexity of Idas symbiosis. In our metagenomic dataset, the relative read abundances of Nitrincolaceae and Thiodubilierella, but not Nitrincolaceae and Thioglobus symbionts, appear to positively correlate (R² = 0.79 and R² = 0.03, respectively, for linear regression of centered-log transformed data). Although compositional data correlation should be treated with caution, our results hint at the interdependence of these two symbiont populations.

Urechidicola symbionts may denitrify to N₂, removing NO

Urechidicola appears to be the only bacterium capable of complete denitrification to N₂. Nitrate and nitrite reduction can be catalyzed by the periplasmic enzymes heterodimeric nitrate reductase (NapAB, no expression found) and the NrfAB cytochrome c552 nitrite reductase (both subunits were expressed at the mRNA level). The norBC genes that encode the respiratory nitric oxide (NO) reductase were highly transcribed. Most importantly, the nosZ gene encoding the periplasmic sec-dependant nitrous oxide reductase was the 11th most abundantly transcribed gene in this bacterium. Both the methane-oxidizing and methylotrophic symbionts can oxidize nitrite to nitric oxide but lack the genes needed to complete the remaining steps of denitrification. Symbiont nitrate respiration can contribute to holobiont fitness, as it fuels energy conservation during hypoxia induced by natural perturbations or when the shell valves are closed and reduce competition between host and symbiont for oxygen [15]. Removal of the toxic nitrite respiration product NO could benefit the symbionts [106]. Whereas the Urechidicola symbiont is only present in low abundance, complete denitrification is likely beneficial mainly for their population, which can take advantage of excess NO produced by the host or the other symbionts. We speculate that in some cases the Urechidicola symbiont may also contribute to NO removal at the holobiont level.

Conclusions

The I. modialiformis symbiosis has a higher number of symbiont species than the typical chemosynthetic associations in most other hosts. The metabolic flexibility in these symbioses expands the range of catabolized substrates and likely allows for the colonization of not only chemosynthetic environments but also organic substrates, such as wood. This nutritional plasticity may have played a role in the adaptation and evolutionary transition of these mussels from organic substrates to the deep sea [107]. Yet, most large bathymodioline mussels, such as Gigantidas and most Bathymodiolus, host only a few species of chemosynthetic bacteria, highlighting the advantage of limited symbiont diversity. Expanding the symbiont diversity and substrate range may lead to energetic costs, resulting in lower growth rates or reduced reproductive output [20]. These costs may be reduced as only the symbiont that can use the locally abundant substrate grows up to high abundances. For example, at seeps and brine pools, where the key substrate is methane, the methanotrophic symbiont dominates in terms of biomass. We hypothesize that even in the absence of relevant amounts of external organic substrates the secondary symbionts might provide additional benefits, which include the recycling and efficient use of resources such as NO, methanol and formate, as well as sharing goods such as vitamin B12 (Fig. 6). Some of these positive interactions have equivalents in free-living communities, suggesting that symbiotic interactions may have evolved based on the existing interactions.

Schematic representation of hypothetical key metabolic handoffs in the Idas symbioses.

Acknowledgements

We thank all individuals who helped during the expeditions, including onboard technical and scientific personnel, the captains and crew of the R/V Bat Galim and E/V Nautilus, and teams operating ROVs “Hercules” and “Yona”. This research used samples and data provided by the E/V Nautilus Exploration Program—expedition NA019. We thank Sharon Tal at the Histology Service Unit at the University of Haifa and Boris Shklyar at the Bioimaging unit, Faculty of the Natural Sciences, University of Haifa. We are grateful to Dr. Heather Maughan for editing and feedback on the manuscript. This study was funded by the Israeli Science Foundation (ISF, grant 913/19 to MR-B), the U.S.-Israel Binational Science Foundation (BSF, grant no. 2019055 to MR-B and MK), the Israel Ministry of Energy (grant no. 221-17-002 to MR-B) and the US National Science Foundation (grant IOS 2003107 to MK).

Data Availability

DNA and RNA raw reads, as well as symbiont metagenome-assembled genomes, were deposited to NCBI with project accession number PRJNA930646. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE [108] partner repository with the dataset identifier PXD040273 (Reviewer access at https://www.ebi.ac.uk/pride/login with Username: [email protected] and Password: 1Td7u0kv).

Author contributions

TZ-K, MR-B and DT conceived this study. TZ-K conducted microscopy. MR-B performed bioinformatics. SV and MK performed proteomics. MR-B and TZ-K compiled the data. MR-B, DT and MK acquired funding. TZ-K, MR-B and MK wrote the paper with the contributions of all co-authors.

Competing interests

The authors declare no competing interests.

References

Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s43705-023-00254-4.