Global biogeography of chemosynthetic symbionts reveals both localized and globally distributed symbiont groups

We sequenced, assembled, and binned metagenomes from 47 individual clams. This resulted in 63 MAGs (metagenome assembled genomes), 53 of which were considered high-quality with greater than 90% completeness and less than 10% contamination (18) (Dataset S1). All newly assembled MAGs from this study are available under GenBank accession numbers SAMN16825223 to SAMN16825285, including raw reads (SRA database project number PRJNA679177). All 63 MAGs were assigned to the genus Ca. Thiodiazotropha. Sixteen metagenomes yielded two distinct Ca. Thiodiazotropha MAGs in one gill. The majority of lucinids we investigated hosted a single symbiont species abundant enough for its genome to be assembled from gill metagenomes. Of these, 15 Clathrolucina costata and 1 Loripes orbiculatus contained two distinct MAGs (Dataset S1). Using DESMAN, no more than two strains were detected within any MAG, indicating that there is also limited symbiont strain diversity in each host individual (Dataset S1).

We reconstructed the phylogenetic relationships of 63 Ca. Thiodiazotropha MAGs we generated to previously described lucinid symbionts and their closest free-living relatives, Sedimenticola thiotaurini and Sedimenticola selenatireducens, using a concatenated amino acid alignment of 43 universally conserved marker genes, with Allochromatium vinosum as an outgroup [alignment on Figshare (19)]. The maximum likelihood phylogenetic reconstruction yielded 10 unique clades, each with an average nucleotide identity (ANI) above 95%, indicating that each clade likely corresponds to a distinct species, two of which have not been previously identified (Dataset S2). Apart from the symbionts of Phacoides pectinatus, Ca. Sedimenticola endophacoides, all other lucinid symbionts belong to a clade that most likely represents a single genus, named Ca. Thiodiazotropha. The first of the two species we identified was represented by 36 MAGs assembled from lucinid hosts from across the Pacific, Atlantic, and Indian Oceans, making it the most widely distributed lucinid symbiont known to date. We propose the name Candidatus Thiodiazotropha taylori for this species after Dr. John Taylor (SI Appendix, Supplementary Discussion), to honor his enormous contribution to understanding the biodiversity and evolution of the Lucinidae. Ca. T. taylori was associated with eight host species from three different subfamilies of Lucinidae (Lucininae, Leucosphaerinae, and Codakiinae), making this the most diverse group of hosts that any lucinid symbiont species is known to associate with (Fig. 1). The second species, Ca. T. sp. “RUGA,” was found in a single Rugalucina munda metagenome. The Ca. T. sp. “RUGA” MAG shared an ANI of 94.543% with all Ca. T. taylori MAGs, making it the closest identified relative to the Ca. T. taylori group.

Reanalysis of L. orbiculatus metagenomes containing a previously described symbiont species, Ca. T. endoloripes, revealed that Ca. T. endoloripes is actually two distinct lineages with an ANI of 91.938%, which indicates they are separate species (Dataset S2). We propose the names Ca. T. weberae and Ca. T. lotti for these two species after Dr. Miriam Weber and Christian Lott, who discovered the population of L. orbiculatus in Elba, Italy. Its symbionts Ca. T. weberae and Ca. T. lotti both share an ANI of roughly 89% with Ca. T. taylori and Ca. T. sp. “RUGA,” their next-closest relatives (Fig. 2). The lucinid symbionts from the genus Ca. Thiodiazotropha formed two major clades, one comprising the four above-mentioned species and another containing Ca. T. endolucinida and the Thiodiazotropha-like species associated with Ctena orbiculata from Florida, previously reported by Lim et al. (20) (Fig. 2). Like Ca. T. taylori, Ca. T. endolucinida was also found in multiple host species (five species from subfamilies Lucininae, Codakiinae, and Leucosphaerinae) but its geographic distribution appears to be restricted to the Caribbean based on samples that are so far available (February 2021).

We compared the predicted functional capabilities of Ca. T. taylori, the most widespread symbiont in our dataset, to 1) its closest relatives (Ca. T. weberae, Ca. T. lotti, and Ca. T. sp. “RUGA”) and 2) the sympatric Ca. T. endolucinida, which can cooccur with Ca. T. taylori in the same host gill, to investigate whether functional potential reflects symbiont geography or phylogeny. All five symbiont species shared most core metabolic functions (Table 1). A maximum of 7.2%, or 170 of the 2,361 protein families (Pfams) found in these genomes, were enriched in Ca. T. endolucinida, i.e., they were significantly more frequent within this species group (P_adj < 0.05), while the closely related species Ca. T. taylori, weberae, and lotti had far fewer enriched Pfam functions, suggesting that phylogenetic relationships are to some extent reflected in genomic functional potential (details in SI Appendix, SI Results, Table S3, and Fig. S2 A and B and Dataset S5). In other words, the more closely related the symbionts, the more likely they are to share genome content.

The Ca. Thiodiazotropha symbionts rely on multiple pathways to oxidize reduced sulfur compounds to sulfate through a polysulfur or elemental sulfur intermediate and there were no major variations in sulfur oxidation capabilities across all the MAGs (Table 1 and SI Appendix, Table S2 and Dataset S3). The energy generated from sulfur oxidation is used to power inorganic carbon fixation through the CBB cycle (Table 1). The MAGs of Ca. T. taylori, T. weberae, T. lotti, T. sp. “RUGA,” and T. endolucinida each contained a gene encoding Ribulose-1,5-bisphosphate carboxylase/oxygenase (RuBisCO) form I (Table 1). Ca. T. endolucinida was the only species that possessed two distinct types of RuBisCO, forms I and II. All five species also encoded the capacity for heterotrophic growth (Table 1). The genus name Thiodiazotropha was proposed following the recent discovery that lucinid symbionts are capable of fixing inorganic nitrogen from the atmosphere (14, 15). Indeed, all five symbionts encoded large gene clusters involved in nitrogen fixation, including the widely used functional marker for nitrogen fixation, dinitrogenase reductase subunit (nifH), the structural genes (nifD and nifK), ferredoxins, and regulatory factors (Dataset S3). All species had a complete denitrification pathway for reducing nitrate and nitrite to nitrogen gas, and thus all had the potential to use nitrate and nitrite as alternative electron acceptors (Table 1).

Despite sharing most metabolic functions that sustain the nutritional lucinid symbiosis, we observed notable variations in symbiont carbon and nitrogen metabolic pathways. Ca. T. taylori, Ca. T. endolucinida, and Ca. T. sp. “RUGA” possessed genes required for growth on reduced one-carbon compounds (20), including lanthanide-dependent methanol dehydrogenase (xoxF-like), genes for the synthesis of its cofactor pyrroloquinoline quinone (pqqABCDE), and the full tetrahydromethanopterin-dependent pathway for formaldehyde oxidation, which was found in close proximity to methanol dehydrogenase (21, 22). The presence of this cluster of genes suggests these symbionts are able to convert methanol to aldehydes, which can be subsequently converted to biomass or utilized as energy. All five species were capable of assimilating nitrogen gas and ammonia but there were notable differences in their ability to assimilate urea (Table 1), a potential waste product of lucinid clams. Only MAGs of Ca. T. taylori and Ca. T. weberae encoded genes for the uptake and conversion of urea to ammonia (ureACDEFGJ). MAGs of Ca. T. taylori, Ca. T. lotti, Ca. T. weberae, and Ca. T. “RUGA” also encoded genes for an alternative pathway for urea metabolism through urea amidolyase (SI Appendix, Table S2). Finally, Ca. T. endolucinida was the only species with MAGs encoding genes for the respiratory nitrate reductase complex (narGHJI), indicating that it can reduce nitrate for generating a proton motive force. These observations were based on the presence/absence of entire operons/gene clusters and were consistent across all the high-quality MAGs examined, of which we had multiple samples from each symbiont species (minimum of nine).

Nucleotide differences in closely related bacterial sequences can be caused by two primary mechanisms, either homologous recombination (HR), where a fragment of the recipient’s genome is replaced by that of a donor in a single generation (genetic exchange), or by mutations from one generation to the next (23). In the absence of HR, all differences in genomic sequences reflect a clonal genealogy. We calculated average HR to mutation rates (R/m) for each species group and for groups containing multiple species. The alignment lengths of shared genes used to calculate R/m ranged from 1,751,700 to 3,149,280 base pairs (bp), which for some groups covered more than half of the total length of the MAGs (Table 2). The alignment of shared genes was shortest in Ca. T. taylori due to the strict alignment parameters chosen (i.e., aligned genes must be present in all samples) and a greater genetic diversity within this species group. Just as Ca. T. taylori is unique among lucinid symbiont species for its extensive geographic range, its average R/m ratio was at least 10 times higher than the average R/m ratios in all the other species groups (0.814; 0.809 to 0.819; 95% confidence intervals; Table 2). The average R/m ratios for Ca. T. taylori populations from different locations were fairly consistent and ranged from 0.859 in Florida to 1.521 in Mauritania (SI Appendix, Table S4). The lowest R/m values were between Ca. T. endolucinida and Ca. T. taylori (0.053; 0.052 to 0.054), even though these can cooccur within host individuals. The lengths of all alignments, average number of recombination events, average length of events, and all R/m ratios and their 95% confidence intervals are shown in SI Appendix, Table S4. The average nucleotide diversity outside the recombining gene sequences was 3.6, 2.4, 1.9, and 0.6 single-nucleotide variants per kbp for Ca. T. taylori, Ca. T. lotti, Ca. T. weberae, and Ca. T. endolucinida, respectively (Dataset S6). Although MAGs consist of consensus sequences of the most abundant representatives of an entire bacterial population, the low strain estimates from DESMAN (two or fewer strains) across all MAGs in this study provide confidence that these results are not confounded by strain differences within the gill metagenomes (SI Appendix, Table S1).

We investigated the diversity and the predicted functional variability of the symbionts associated with lucinids on a global scale. Fresh samples were complemented with specimens from museum collections to expand the scope of the study beyond the most intensively studied lucinid habitats in the Caribbean and Mediterranean Seas (15, 20, 24). We identified two lucinid symbiont species, Ca. T. taylori and Ca. T. “RUGA.” Ca. T. taylori is found in a remarkable eight different lucinid species from three different subfamilies (Figs. 1 and 2), which makes Ca. T. taylori the most promiscuous lucinid symbiont described thus far and the first chemosynthetic endosymbiont species with a globally distributed population (25, 26). We similarly found Ca. T. endolucinida in five divergent host species from three distinct lucinid subfamilies (Fig. 1). These findings corroborate previous reports that the same symbiont 16S ribosomal RNA (rRNA) gene sequence, identical to that of Ca. T. endolucinida (SI Appendix, Fig. S3), was present in four distinct host species at one location in Guadeloupe (25). We predict that future surveys of lucinid symbiont diversity are likely to reveal similar instances of promiscuity. The remarkable flexibility in the association between lucinids and Ca. Thiodiazotropha species could thus be an important feature underlying the evolutionary success of this ancient and widespread symbiosis.

Genetic exchange plays a critical role in bacterial genome evolution and can either drive the divergence or homogenization of a population (27). HR drives homogenization by maintaining genomic cohesion of bacterial clades even across distant locations (28, 29). We reconstructed recombination events in four symbiont species groups—Ca. T. taylori, T. endolucinida, T. weberae, and T. lotti—to investigate the role of HR in maintaining genomic cohesion of Ca. T. taylori as a single species-level group. The rate of recombination to mutation (R/m) in Ca. T. taylori (0.814) is above the theoretical threshold for preventing population divergence (0.25) (29–31). In contrast, the R/m rates of Ca. T. endolucinida, T. weberae, and T. lotti are well below this 0.25 threshold (0.04 to 0.08), a pattern that correlates with the limited geographic range of these three species compared to Ca. T. taylori (Fig. 1). With such a large distribution range spanning many diverse habitats, selection pressures or neutral drift could cause local Ca. T. taylori populations to diverge, which might result in gene content and functional differences (Table 1). This is supported by our phylogenetic reconstruction of the Ca. T. taylori clonal frame (Fig. 2), which reveals that clams from each distinct geographic location are colonized by site-specific phylogenetic lineages that were poorly resolved in the phylogenomic tree (Fig. 1). Our recombination analyses indicate that HR has a cohesive effect that maintains the integrity of Ca. T. taylori as a single cosmopolitan species-level group. With an ANI of about 94.5%, the Ca. T. taylori and Ca. T. “RUGA” MAGs represent distinct, albeit closely related, genetic units (Fig. 1 and Dataset S1). Given that HR rates drop exponentially with decreasing sequence similarity, declining steeply between 90 to 95% ANI (32), we speculate that this seemingly small amount of nucleotide divergence likely prevents recombination between Ca. T. “RUGA” and Ca. T. taylori. Indeed, an ANI of 95% is widely used as a threshold for delineating bacterial species, as this is thought to correspond to a level of divergence that prevents genetic exchange by HR (33–35). One could speculate that Ca. T. RUGA and Ca. T. taylori previously belonged to a single species-level group but were relatively recently separated and that if the local geographic Ca. T. taylori populations encountered a barrier to HR they too might become new and distinct species-level groups. It is currently unclear which barriers could “remove” local populations from the cohesive forces of HR and allow them to diverge independently. Greater efforts to sample and sequence more Ca. T. “RUGA” and Ca. T. taylori symbionts will allow us to test these and other theories about the emergence of symbiont diversity in future.

The cosmopolitan distribution of lucinid clams hosting Ca. T. taylori suggests this symbiont species is able to disperse over great geographic distances. How Ca. T. taylori achieves this feat of dispersal remains a mystery but hitchhiking on their bivalve hosts during the planktonic larval phase appears unlikely because the lucinid larvae studied thus far have been found to be aposymbiotic (11). Furthermore, the range of this symbiont species far exceeds the limited distribution range of any lucinid host species. Ca. T. taylori therefore likely migrates during its free-living phase rather than its host-associated phase. The dormant endospores of thermophilic Firmicutes achieve extensive distributions through long-distance passive dispersal in oceanic currents, but even these tough endospores with enhanced survival capacities face more substantial dispersal limitations than Ca. T. taylori (36). It seems unlikely that free-living Ca. T. taylori cells would be able to survive harsh open-ocean conditions long enough to traverse the globe. To our knowledge, there are no molecular data to date indicating Ca. Thiodazotropha spp. are present at a high abundance in coastal sediments or the water column, which suggests they are members of the rare biosphere present only at low relative abundance (37). Troussellier et al. put forward the intriguing hypothesis that macroorganisms could sustain the rare biosphere by serving as dissemination vectors for marine microbes (38). A recent meta-analysis of publicly available amplicon sequencing data revealed the presence of Ca. Thiodiazatropha 16S rRNA gene sequences associated with the roots of various seagrass species around the globe, a discovery the authors subsequently verified with microscopic imaging (37). Vegetative seagrass fragments (shoots and rhizomes) can reestablish in distant locations after detaching from their parents and have great potential for long-distance dispersal (39, 40). It is interesting to speculate that lucinid clams and seagrasses from across the globe could form a network of source habitats facilitating the dissemination and dispersal of Ca. T. taylori vectored by oceanic circulation, thereby overcoming potential barriers to dispersal between geographically distant bodies of water.

Despite its global distribution range, we may have identified one potential barrier to dispersal of Ca. T. taylori. Rigorous sampling of L. orbiculatus along the Atlantic and Mediterranean coasts of Europe did not reveal a single instance of Ca. T. taylori associated with L. orbiculatus in any of these locations (Fig. 1). Instead, L. orbiculatus in these locations all hosted the closely related sister species Ca. T. weberae and T. lotti (Fig. 1). The colonization of L. orbiculatus by Ca. T. lotti from the United Kingdom all the way south to Kotor, Montenegro suggests the strait of Gibraltar is unlikely to pose a barrier to the dispersal of Ca. T. taylori. Nor are these distribution patterns driven by host-symbiont specificity, because L. orbiculatus in Mauritania hosts Ca. T. taylori (Fig. 1). A possible explanation is that environmental factors associated with a temperate climate prevent Ca. T. taylori from establishing in clams along the coasts of Europe. This is consistent with the tropical distribution of all lucinid species so far found to host Ca. T. taylori (Fig. 1). An alternate nonmutually exclusive explanation could be that Ca. T. weberae and T. lotti are better adapted and able to outcompete Ca. T. taylori for establishment in lucinid hosts throughout temperate Europe. Further sampling of other lucinid species and seagrasses along European coasts will be necessary to address these questions.

Lim et al. (20) recently reported multiple Ca. Thiodiazotropha 16S rRNA gene amplicon sequence variants and MAGs of diverse Ca. Thiodiazotropha spp. associated with C. orbiculata individuals from Florida, providing the first indications that multiple symbiont species from the same genus may coexist in the same host gill. We were able to assemble and bin MAGs of two distinct symbiont species, Ca. T. taylori and T. endolucinida, from the metagenomes of 15 C. costata individuals sampled from across the Caribbean. Using the same methods, a reanalysis of new and previously published gill metagenome data (LVJZ00000000) from L. orbiculatus (Elba, Italy) similarly revealed that MAGs from two distinct symbiont species, Ca. T. lotti and T. weberae, can be assembled and binned from metagenomes of single host individuals. Our fluorescence in situ hybridization (FISH) results showed that in C. costata, Ca. T. taylori and T. endolucinida both inhabit gill epithelial cells. These findings add to a growing number of studies reporting the coexistence of closely related chemosynthetic symbiont species/strains within the same invertebrate host and indicate that this is more common in lucinids than previously assumed (41–43).

Ecological models predict that cooccurring symbionts with the same resource requirements will compete for limited resources and that this competition is detrimental to the symbiosis (e.g., ref. 44). The cooccurring symbiont pairs we identified, Ca. T. taylori/endolucinida and Ca. T. weberae/lotti, shared core metabolic functions fundamental to their symbiosis with lucinid clams (15), including identical pathways for oxidizing sulfur, fixing inorganic carbon, and fixing nitrogen. Ca. T. taylori and Ca. T. endolucinida both had the additional potential to utilize methanol as a source of energy and carbon with an xox-type methanol dehydrogenase and the serine pathway for C1-carbon incorporation into biomass. Compartmentalizing coexisting symbiont species into separate bacteriocytes could prevent direct competition by allowing the host to partition resources and discriminate cooperative symbionts from potential cheaters that might destabilize the symbiosis (45). Consistent with this, our FISH analyses showed that although symbiont species cooccurred in host individuals they never cooccurred in single host bacteriocytes (Fig. 2). Despite their prolific productivity, habitats abundant with lucinids, such as seagrass meadows and coral reef lagoons, tend to occur in oligotrophic waters that are nitrogen-limited (14). It is noteworthy that the ability to hydrolyze urea to ammonia, which can subsequently be assimilated, is a conserved function of Ca. T. taylori and Ca. T. weberae that is absent in Ca. T. endolucinida and Ca. T. lotti (Table 1). The precise composition of lucinid nitrogenous waste remains unknown, but some bivalves do excrete urea as a waste product and the urea transporter DUR3 is highly expressed in the gills of L. orbiculatus (46, 47). Future studies are required to investigate whether the ability to utilize urea as a nitrogen source confers any advantage to Ca. T. taylori and Ca. T. weberae in the host or the external environment.

We observed some additional predicted metabolic differences between Ca. T. taylori and Ca. T. endolucinida. This symbiont pair has one major difference: Ca. T. endolucinida MAGs encode both RuBisCo forms I and II, while Ca. T taylori only encoded form I (Table 1). Ca. T. endolucinida also has the genes encoding a respiratory nitrate reductase protein complex, which indicates the ability to use nitrate instead of oxygen as an electron acceptor for respiration (Table 1). This combination of functional traits suggests Ca. T. endolucinida may be better adapted to survival in a lower-oxygen environment than Ca. T. taylori. This could be highly beneficial for the symbiosis as the endosymbionts of both the lucinid clam Lucinoma aequizonata and the vent tubeworm Riftia pachyptila are able to respire nitrate as an adaptation to living in deep-sea environments with fluctuating oxygen availability (48–50). Whether shallow-water lucinids in the Caribbean are able to exploit this unique metabolic capability of Ca. T. endolucinida remains unknown but there are several observations suggesting it does not. First, the vanishingly small nitrate reductase activity in the gills of Codakia orbicularis, an established host of Ca. T. endolucinida, and the absence of nitrates in the tissues of this host species together suggest nitrate is not used as an electron acceptor in this symbiosis (51). Second, nitrate is only sporadically present at low concentrations in the pore water of Thalassia testudinum sediments and is undetectable in the overlying waters (14, 51). Third, oxygen is abundantly available in Caribbean lucinid habitats, especially during periods of photosynthesis, and lucinid clams construct burrows leading to the surface to obtain oxygenated water from above the sediment layer (52). These observations are inconsistent with Ca. T. endolucinida’s requiring these metabolic adaptations while they are housed in the buffered environment of the lucinid bacteriocytes. Moreover, our analysis of the C. orbicularis metatranscriptome indicates that the genes encoding RuBisCO form II and respiratory nitrate reductase are expressed at a much lower level compared to RuBisCo form I and the assimilatory nitrate reductases, respectively (Dataset S7). These data are further supported by the much lower abundance of RuBisCo form II proteins compared to form I in the C. orbicularis proteome (14). Based on these observations, we hypothesize that Ca. T. endolucinida relies on these metabolic functions primarily during its free-living phase and we speculate that Ca. T. endolucinida could occupy an external environmental niche—one characterized by low oxygen availability— distinct from that of Ca. T. taylori.

Live clams were collected from sites throughout the Caribbean and Mediterranean (SI Appendix, SI Methods and Dataset S1). Gills were dissected in the field and preserved in RNAlater (AM7020; Life Technologies) or DNA/RNA Shield (R1100-250; ZymoBiomics) according to manufacturer’s instructions and stored at −20 °C. Samples from locations in the Pacific and Indian Oceans came from the collections of the Natural History Museum (NHM) in London and the Florida Natural History Museum (FLMNH), Gainesville, FL (SI Appendix, SI Methods and Dataset S1). Access was permitted and organized by Dr. John Taylor (NHM) and Dr. Gustav Paulay and Dr. Amanda Bemis (FLMNH). SI Appendix, SI Methods and Dataset S1 list all the species used, sampling locations, dates, and sample sizes.

DNA was extracted from gill tissues using the Qiagen DNeasy Blood and Tissue kit (69506; Qiagen) or the animal tissue protocol from Analytikjena Innuprep DNA Mini Kit (845-KS-1041250) (SI Appendix, SI Methods). Samples were either treated with RNase or directly quantified before DNA libraries were prepared using Illumina-compatible library prep kits (Dataset S1 and SI Appendix, SI Methods and Table S5). All libraries were sequenced with Illumina technology to generate paired-end reads of 150 bp or 250 bp length (SI Appendix, SI Methods and Dataset S1).

Read libraries were trimmed, PhiX contamination-filtered, and quality-checked using BBMAP v37.61’s BBDUK feature (53); parameters used are in SI Appendix, SI Methods and the Jupyter notebook. Individual read libraries were assembled using SPAdes v3.13.1 (54, 55); parameters used are in SI Appendix, SI Methods. The resulting metagenomic assembly scaffolds were binned using a combination of Anvi'o v6.1 (56, 57) using CONCOCT v1.1.0 (58) and MetaBAT v2.15 (59) (details and parameters in SI Appendix, SI Methods and Dataset S1). The bins were then compared using dRep v2.4.2’s dereplicate workflow (60) (SI Appendix, SI Methods). The bins were checked for completion using Checkm’s taxonomy specific workflow and manually refined using “anvi-refine” (SI Appendix, SI Methods). MAGs that were determined to be 90% or more complete and less than 10% contaminated post refinement, referred to as high-quality MAGs, were used for further analyses. Potential strain numbers in individual clams/metagenomes were calculated using DESMAN v2.1.1 (61) and the snakemake workflow (62) on the program (SI Appendix, SI Methods).

Only the MAGs taxonomically assigned to Sedimenticolaceae or Chromatiaceae, using GTDB v0.3.3 (63) were used in this study. We also downloaded publically available MAGs of other lucinid symbionts (Ca.Thiodiazotropha spp. and Sedimenticola spp.), alongside A. vinosum as an outgroup, for this analysis (complete list of accession numbers in Dataset S1). The publicly available MAGs were quality-checked and filtered using the methods described above (and in SI Appendix, SI Methods). The CheckM v1.1.3 (64) lineage workflow was used to identify, align, and concatenate a default set of 43 universal marker genes from all the available MAGs (Dataset S1); concatenated marker gene alignment can be found on Figshare (19). This concatenated amino acid alignment was then submitted to the W-IQ-TREE server (65) using default settings (SI Appendix, SI Methods) and the resulting maximum likelihood tree was visualized using the Interactive Tree Of Life (iTOL) v5 (66). All MAGs were placed into species groups based on ANI values above 95% using FastANI v1.3 (34). The accuracy of these species boundaries was tested with ANIb and ANIm through the jspecies web server (67).

We carried out FISH to visualize Ca. T. taylori in gill sections of L. orbiculatus (Mauritania, 2018) and C. costata (Panama, 2019). Probes were designed to target the 16S rRNA gene sequence of Ca. T. taylori using DECIPHER’s design probes web tool (SI Appendix, SI Methods) (68). A formamide series from 0 to 60%, in 10% steps, was carried out to optimize the probe hybridization conditions (SI Appendix, SI Methods; probe sequences used are in SI Appendix, Table S6). Nonsense sequences of the target-specific probes were also tested as negative controls (69). Dissected gills were preserved in 4% paraformaldehyde, dehydrated into 70% ethanol, and stored at 4 °C (SI Appendix, SI Methods). The gills were embedded in paraffin wax by the Histopathology Facility at Vienna BioCenter Core Facilities, Austria (SI Appendix, SI Methods). We cut the embedded gills into 5-µm sections with a Leica microtome and mounted the sections on SuperfrostPlus adhesion slides (Thermo Scientific). The sections were dewaxed in Roti-Histol (Carl Roth) following the manufacturer’s instructions. Probes were hybridized to the gill sections in a 50% formamide buffer (details on hybridization conditions and washing protocols are in SI Appendix, SI Methods and Tables S6 and S7). Following the hybridization, the samples were DAPI-stained (1 µg/mL) and mounted in ProLong Glass antifade mounting media (Thermo Fisher Scientific) (SI Appendix, SI Methods). Images were captured on a Leica TCS SP8 X confocal laser scanning microscope using a 63× objective (SI Appendix, SI Methods).

Anvi’o’s pangenomic workflow was used for orthologous group clustering and functional comparisons (SI Appendix, SI Methods). We annotated all features containing open reading frames (ORFs) using eggNOG-mapper v2 (70) with eggNOG database v5.0 (71). All ORFs were also annotated with Pfam domains (72) using HMMER v3.3 (73) (parameters in SI Appendix, SI Methods). All high-quality MAGs were also used to create pangenome in Anvi’o with an mcl inflation value of 8; pangenomic genomes database and profiles are on Figshare (74, 75). The MAGs were also annotated on the Rapid Annotation using Subsystem Technology (RAST) web server (https://rast.nmpdr.org/) using the RASTtk pipeline (76). Where necessary, genes of interest were manually screened using NCBI blast+ v2.8.1 (77) and BBMAP to search for genes potentially missing from the assemblies (SI Appendix, SI Methods). To determine if any annotated functions were present in a given species group or population at a higher frequency than expected under a uniform distribution, EggNOG and Pfam terms were statistically tested for enrichment across different species groups and populations through the “anvi-get-enriched-functions-per-pan-group” function in Anvi’o with an “adjusted q value” cutoff of 0.05 (SI Appendix, SI Methods).

Three C. orbicularis specimens were collected in Guadeloupe (2016), preserved in RNAlater, and stored at −80 °C. Total RNA was extracted using TRIzol (Thermo Fisher Scientific) according to the manufacturer's instructions. rRNA depletion and library construction were carried out at the Vienna Biocenter Core Facilities GmbH as described in Yuen et al. (47). The RNA-sequencing reads were trimmed and processed as previously described by Yuen et al. (47). We used BBMAP (slow = t, ambiguous = best, minid = 0.99) to align the processed reads to the MAG of Ca. T. endolucinida (GCA_001715975.1), an endosymbiont of C. orbicularis from Guadeloupe (53). FeatureCounts was used to quantify gene-level transcript abundances, which were subsequently converted to transcripts per million (78). Metatranscriptomic results can be found in Dataset S7.

To infer recombination events in bacterial genomes, we used the maximum likelihood implementation of ClonalFrame, ClonalFrameML (23). High-quality MAGs were aligned using progressiveMauve, in Mauve v2.0 (79). Nucleotide sequences (i.e., core genes) shared within the groups of MAGs analyzed were extracted with stripSubsetLCBs, a script previously described in ref. 80. These core genes were realigned with MUSCLE (81) and cleaned with trimAl (82) (parameters in SI Appendix, SI Methods). RAxML v8.2.10 (83) was used to build a phylogenetic tree from this new alignment as described in ref. 84. All alignments are available on Figshare (85–87). The resulting tree and alignment were fed into ClonalFrameML to calculate the ratio of recombination vs. mutation events with default parameters (See Dataset S1 for MAGs used in this analysis). Recombination events were computed and visualized in R v6.3.2 with the script cfml_results.R (SI Appendix, SI Methods). GUBBINS was used to estimate average nucleotide diversities inside and outside the clonal frame (88). We also used the clonal frame alignments to reconstruct and visualize the phylogenetic relationships of the Ca. T. taylori MAGs as previously described (SI Appendix, SI Methods). All scripts used can be found in the associated Jupyter Notebook (89).