Bathymodiolus marisindicus mussels were collected from four hydrothermal vent fields on the CIR ( Table 1 , Figure 1 ). The mussels were obtained from the Onnuri and the Solitaire vent fields using a video-guided hydraulic grab (Oktopus GmbH, Germany, Dive number: GTV1702, GTV1807, and GTV1809) and from the Edmond and Kairei vent fields using a remotely operated vehicle, Jason I (Dive number: JL301, JL296). All samples were frozen at −80°C on board, transported to the laboratory on dry ice, and stored at −80°C until required for downstream analyses.

The genomic DNA of the mussels was extracted from each specimen’s foot or mantle tissue using the Qiagen DNeasy Tissue kit (Qiagen Inc., Hilden, Germany). One mitochondrial gene (mtCOI) and two nuclear genes (EF1α and Col-1) were amplified using gene-specific primer sets with primer-specific PCR conditions previously designed for bathymodioline mussels ( Supplementary Table S1 ). The PCR was performed in 20 μl volumes including 2 μl of 10× Taq polymerase buffer, 1.5 μl of 2.5 mM stock solution of dNTPs, 1.5 μl of each primer (10 μmol/l), 1 μl of 1 mg/ml bovine serum albumin, 1 μl of extracted DNA (30–150 ng), 0.625 units of IP Taq polymerase (Cosmo Genetech Inc., South Korea), and sterile H₂O to reach the desired volume. The PCR products were purified using the Dr. Prep Kit (Cat. No. MK02020, MGmed, South Korea). The DNA sequences were then determined using a Big Dye Terminator V3.1 Cycle Sequencing Kit on an Applied Biosystems 3730xl DNA Analyzer (Applied Biosystems Inc., South Korea). All sequences from the host mussels were deposited in the GenBank database under accession numbers OK509086–OK509159 (mtCOI), OK524315–OK524440 (EF1a), and OK524441–OK524584 (Col).

The genomic DNA of symbiotic bacteria was extracted from the gill tissues of individual mussels using the Qiagen DNeasy Tissue Kit (Qiagen Inc., Hilden, Germany). We targeted single-copy genes including 16S rRNA and six functional genes with different functional categories used in previous studies on sulfur-oxidizing symbionts of deep-sea mussels of the genus Bathymodiolus: sulfur metabolism (dsrB), protein folding (dnaK), glycolysis (pfk, pgi, and pykF), transcription regulation (rpoD) (Fontanez and Cavanaugh, 2014; Ho et al., 2017). Libraries for Illumina-MiSeq sequencing were prepared via two-step PCR. First, fragments of seven genes were amplified using specific primer sets with the following conditions ( Supplementary Table S2 ): 95°C for 3 min; 30 cycles at 95°C for 40 s, 50°C–55°C for 30 s, and 72°C for 90 s, and 72°C for 7 min. The 16S rRNA gene was amplified using a universal prokaryotic primer set targeting the V3 and V4 regions. Functional genes were amplified using a thiotroph-specific primer set designed based on the previously deposited sequences of B. marisindicus sulfur-oxidizing endosymbionts and newly generated sequences by PacBio RS II (unpublished data sequenced from gill tissue of specimen “S10” collected at the Solitaire field in 2017). In the second PCR, unique pairs of Illumina Nextera Indices (N7/S5) were attached to the end of the amplified fragments to separate individual hosts. After PCR purification using AMpure XP magnetic beads (Edwards, 2012), the amplicons were pooled at equal quantities (80 ng). Finally, the DNA of the symbiont communities was sequenced on an Illumina MiSeq sequencer (600PE) following a 2 × 300 bp paired-end protocol at the National Instrumentation Center for Environmental Management (NICEM) of Seoul National University (Seoul, Korea). The raw Illumina amplicon data were deposited in the GenBank database under accession number PRJNA769986.

The raw reads of the functional genes were processed using Cutadapt v. 2.3 (Martin, 2011) and the ‘DADA2’ R package taking into account the DADA2 algorithm designed to detect a fine-scale variation from both 454 and Illumina data with more accuracy (Rosen et al., 2015; Callahan et al., 2016). The paired-end reads were demultiplexed according to gene region, and primer sequences were trimmed using Cutadapt v. 2.3. The DADA2 R package was used for filtering by quality, denoising, merging paired reads, and removing chimeric sequences. These procedures were conducted using the default settings except when removing chimeric sequences to prevent excessive deletion of sequences, in which case the ‘minfoldParentOverAbundance (default = 2)’ option was modified to 4. Subsequently, the sequences with <90% similarity with other reads and single reads of each locus were excluded in subsequent analyses.

Raw 454 pyrosequencing data from gill-associated symbionts of bathymodiolin mussels at seven vent fields in the eastern Pacific Ocean were obtained from a previously published study (Ho et al., 2017). Our study analyzed four common loci (dnaK, pgi, pykF, and rpoD) used for molecular analyses of symbionts of the CIR. The sequencing data from the CIR and the eastern Pacific symbionts were obtained by two different sequencing techniques (Illumina MiSeq and 454 pyrosequencing platforms, respectively). We used the same pipeline, DADA2, for the NGS data process to minimize the potential bias caused by the usage of different denoising algorithms. The DADA2 pipeline was conducted with the appropriate settings for the pyrosequencing dataset to denoise gap errors in homopolymer sites and restrict the number of insertions that the 454 platform frequently generates (HOMOPOLYMER_GAP_PENALTY = −1, BAND_SIZE = 32). Likewise with the processing of the Illumina dataset, the sequences with <90% similarity and single reads of each locus were excluded in subsequent analyses. The pipelines for the analysis of raw sequence data from Illumina MiSeq and 454 pyrosequencing platforms are schematized in Supplementary Material S1 .

The diversity and relative abundance of each bacterial strains in the gill tissue of B. marisindicus were characterized by ribotyping of amplicons of the V3–V4 regions of the 16S rRNA gene. The raw reads were first quality-checked, after which low-quality (<Q25) reads were filtered using Trimmomatic v. 0.32 (Bolger et al., 2014). Afterward, the paired-end sequence data were merged with VSEARCH v. 2.13.4 (Rognes et al., 2016) with the default parameters. Primers were then trimmed with the alignment algorithm of Myers and Miller (1988) at a 0.8 similarity threshold. Non-specific amplicons that do not encode 16S rRNA were detected using nhmmer 4 in the HMMER software package v. 3.2.1 with hmm profiles. Unique reads were extracted, and redundant reads were clustered with the unique reads using VSEARCH. The EzBioCloud 16S rRNA database (Yoon et al., 2017) was used for taxonomic assignment followed by more precise pairwise alignment using VSEARCH. The cutoff values of sequence similarity were taken from Yarza et al. (2014). Chimeric reads were filtered based on a <97% similarity threshold via reference-based chimeric detection using the UCHIME algorithm (Edgar et al., 2011) and the non-chimeric 16S rRNA database from EzBioCloud. Finally, operational taxonomic units (OTUs) with single reads were omitted from further analysis. The pipelines for the analysis of the 16S rRNA raw sequence data are schematized in Supplementary Material S2 .

A total of 5,776 to 26,558 gill-associated bacterial 16S rRNA sequence reads were obtained per host mussel from the CIR ( Supplementary Table S3 ). The bacterial reads were identified as Gammaproteobacteria, most of which were classified as sulfur-oxidizing bacteria within the genus Candidatus Thioglobus (>80%). Most reads of sulfur-oxidizing bacteria had high sequence similarity (>99%). Next, a small number of reads were identified as bacteria within the genus Kistimonas of the family Endozoicomonadaceae regardless of the sampling location ( Figure 3 , Supplementary Table S3 ). Likewise, the Candidatus Endonucleobacter within the same family Endozoicomonadacea was detected broadly in the Onnuri, Solitaire, and Kairei fields despite the considerably low abundance (<1.28% of reads). Most of the reads were identified as Candidatus Endonucleobacter bathymodiolin based on sequence similarity >97%.

The number of reads between Illumina and 454 pyrosequencing platforms was considerably different (mean number of reads per gen: Illumina MiSeq data = 94,946–144,842, 454 pyrosequencing = 4,939–7,260; Table 2 , Supplementary Table S4 ). Before the further examination of the genetic divergence of symbiotic bacteria, the nucleotide diversity (π) of each vent locality was measured and compared to each other to check if the different platforms have led to a possible systematic underestimate or overestimate of the population genetic index of D _xy. The nucleotide diversities were considerably low in the vents of the eastern Pacific ridges, of which symbiont sequences were generated by 454 pyrosequencing, compared to those from the vents of the CIR of which sequences were generated by Illumina MiSeq (mean nucleotide diversity: eastern Pacific vents = 0.0015–0.0051, central Indian vents = 0.0031–0.0153; Supplementary Figure S1 ). However, two vent fields, Edmond and Kairei, from the southernmost ridge segment of the CIR were similarly low in nucleotide diversity and equivalent to those of the eastern Pacific vents. In contrast, the symbiont populations in Onnuri and Solitaire vent fields exhibited higher nucleotide diversity than others. The low nucleotide diversity in the eastern Pacific by the 454 pyrosequencing was contrary to the expectation of higher sequencing error rates of this platform. Therefore, the possible introduction of systematic overestimation or underestimation into the present study due to the difference in sequencing platforms was unlikely.

The symbiont populations exhibited a distinct geographical differentiation pattern between the two different ocean basins ( Figure 4 ). The comparison between the CIR and the eastern Pacific Ocean showed the highest range of genetic differentiation D _xy (mean = 0.271–0.278; Figure 4 ; Supplementary Table S5 ) followed by the estimates between the East Pacific Rise (EPR) including the Galápagos Rift (GAR) and the Pacific Antarctic Ridge (PAR) symbiont populations (mean = 0.019–0.021). The D _xy values between CIR symbiont populations ranged from 0.007 to 0.017. Finally, the lowest D _xy values between symbiont populations within the same ridges (EPR including GAR, and PAR) in the eastern Pacific Ocean ranged from 0.002 to 0.004. The Mantel test indicated that the mean of genetic differentiation of symbionts in the CIR increased sharply and proportionally with the geographical distance, but this effect was not statistically significant (Mantel’s r = 0.912, P = 0.08). However, the Pearson test showed a statistically significant positive correlation between the two variables (Pearson’s r = 0.836, P = 0.011). Meanwhile, genetically homogeneous populations were observed in the same ridges in the eastern Pacific Ocean regardless of the geographic distance (partial Mantel’s r = 0.13, P = 0.268; Pearson’s r = 0.001, P = 0.871).

To the best of our knowledge, our study is the first to evaluate the quantitative composition of gill-associated symbiotic bacteria of B. marisindicus. Our findings revealed that the mussels harbor symbiotic communities in which sulfur-oxidizing bacteria dominate (minimum >80%) regardless of the geographical distribution along the CIR, which is consistent with previous findings (Won et al., 2008; Fontanez and Cavanaugh, 2014). On the other hand, from the bacterial community of the gill tissue, we unexpectedly identified several strains belonging to the genus Kistimonas of the family Endozoicomonadaceae, typically found in shallow marine invertebrates such as starfish, dead clams, and sandworms (Choi et al., 2010; Lee et al., 2012; Ellis et al., 2019). The genus Kistimonas was recently observed in the gills of the deep-sea mussel Gigantidas haimaensis (Xu et al., 2019). Recently, other bacterial taxa of the family Endozoicomonadaceae have been identified from various marine mollusks globally (Cano et al., 2020). It suggests that Kistimonas sp. is also widely distributed in association with marine invertebrate animals at different ocean depths. Meanwhile, the Candidatus Endonucleobacter bathymodioli of the Endozoicomonadaceae is known as parasitic bacteria in nuclei of Bathymodioline mussels (Zielinski et al., 2009). Before our research, the parasites were found in bathymodioline mussels in the Gulf of Mexico, Atlantic, and eastern Pacific oceans (Zielinski et al., 2009). Interestingly, the parasites in that study were not observed in the symbiont-infected cell. Therefore, it is considered that the sulfur-oxidizing bacteria of gill tissue seem to protect the host mussel. This interpretation is supported by the antibacterial activity of mussel gill homogenates and the presence of related genes from gill tissue (Bettencourt et al., 2007). The very low abundance of the Candidatus Endonucleobacter bathymodioli compared to the sulfur-oxidizing bacteria is consistent with the expectation according to the hypothesis. Still, the protection mechanism is not known yet.

The fact that the two NGS DNA data sets of eastern Pacific and CIR were sequenced by the 454 and Illumina platforms, respectively, might cause a concern of systematic errors in the estimation of D _xy between vent localities due to the difference in platforms. Indeed, the two sequencing platforms used in the present study have different sequencing error rates. The Illumina sequencing generates fewer sequencing errors than the 454 platform (Loman et al., 2012). Despite the substantial differences in read length and sequencing process of the two platforms, however, the 454 and Illumina sequencing platforms showed reliable quantitative results about the abundances of genes and genotypes (Luo et al., 2012). The algorithm, DADA2, used for sequencing error correction in this study, was proven for its high performance in identifying more real variants and fewer false-positive sequence variants than other methods (Callahan et al., 2016). The result of nucleotide diversity (π) of symbionts in each vent locality indirectly dispels the concern ( Supplementary Figure S1 ). The nucleotide diversities of the eastern Pacific symbiont populations assessed by the 454 platform were overall lower than those of the CIR that were by the Illumina platform, indicating that the suspected systematic errors cannot have affected our conclusion. In addition, the varying degrees of diversity estimated in the CIR mean that the Illumina platform did not seem to cause a systematic decrease or increase in the variation identification. Therefore, we conclude that the nucleotide diversity of populations and D _xy values between them reflect the actual attributes of the mussel symbionts in different mid-ocean ridges.

The comparison of genetic variation of vent mussel symbionts from the eastern Pacific and CIR allowed us to evaluate various factors influencing the population divergence of symbiotic bacteria. The level of genetic differentiation could be partitioned according to the following geographically hierarchical categories in decreasing order: 1) between oceans, the eastern Pacific, and CIR; 2) inter-ridge site between the EPR and PAR across the well-known physical barrier of Ester Microplate in the eastern Pacific; 3) intra-ridge vent fields in the CIR; and 4) intra-ridge vent fields in the EPR and PAR. These categories involve various geographical, geological, and hydrological features of mid-ocean ridges, such as distances between habitats, depths, ridge-axis offsets, transform faults, ridge spreading rates, microplates, and deep-oceanic currents, and the interaction of these factors.

First, the mussel symbiotic bacteria from the eastern Pacific and the CIR evolved independently for a long time. A phylogenetic tree based on 16S rRNA gene sequences showed that the two symbionts markedly diverged (Won et al., 2008). Accordingly, the D _xy estimates between the oceans are at least 10 times larger than the other subordinate comparisons. Secondly, at the inter-ridge scale of the eastern Pacific, the genetic differentiation was high. The Easter Microplate located at the boundary of the EPR and PAR was observed to impede the dispersal of diverse vent invertebrates and gave rise to a genetic isolation of several vent species (Guinot and Hurtado, 2003; Mateos et al., 2012; Johnson et al., 2013; Jang et al., 2016). The congruent vicariance among different taxonomic groups can be interpreted as a result of the orogeny of the Easter Microplate, 2.5–5.3 million years ago around the boundary of EPR and the PAR (Naar and Hey, 1991). Ho et al. (2017) reported that the EPR and PAR mussel symbionts exhibited complete geographical isolation independent of the hybridization of host species, B. thermophilus and B. antarcticus. For the first time, these previous studies (Johnson et al., 2013; Ho et al., 2017) indicated that the geomorphological process had driven the population divergence of both host animals and symbiotic bacteria in the eastern Pacific ridges.

Next, the symbionts of intra-ridge on the CIR exhibited a high degree of differentiation with geographical distance in contrast with the symbiont of intra-ridge sites on the eastern Pacific. In fragmented habitats such as deep-sea hydrothermal vents, the geographical distribution of habitats and life history features of species related to dispersal capability could be the main factors that affect genetic connectivity. This notion suggests that the different geographical distributions of vent habitats must have led to the contrasting genetic differentiation patterns between the eastern Pacific and CIR if we accept that the two symbionts may have similar mobility because they belong to the same taxa, Candidatus Thioglobus.

The seafloor spreading rates of ridges are relevant to the geomorphological features, including the distribution pattern of vent habitats along the ridge axes. Contrary to the relatively linear topology of the fast-spreading ridge, disconnection of adjacent ridges by lateral offset frequently occurs, and an extensive depth range often occurs along the slow-spreading ridge. Therefore, the geomorphological feature in the slow-spreading ridge tends to decrease the genetic connectivity of vent organisms and leads to population subdivision and even speciation (Van Dover et al., 2002; Vrijenhoek, 2010). For example, the Blanco Transform Fault (TF) in the Juan de Fuca Ridge, with approximately 450 km offset in the northeast Pacific, contributed to the genetic separation of the vent tubeworm, Ridgeia piscesae (Young et al., 2008), and the speciation of vent limpets of the genus Lepetodrilus between the northern and southern regions (Johnson et al., 2006). The bathymodioline mussels in the MAR also exhibited a substantial genetic subdivision reaching interspecies levels between the northern and southern MAR across numerous transform faults around the equator, including the largest Romanche TF with an offset of approximately 935-km width and 7-km depth (Van Der Heijden et al., 2012). Likewise, the disconnection of ridge axes by transform faults was observed in the CIR; the Argo and Marie Celeste transform faults between the Onnuri and Solitaire fields offset the ridge by about 100 and 200 km, respectively, and up to 5 km in depth (Pak et al., 2017). The ridge segmentation of the ridge also appears to interrupt the geographical connectivity of symbiotic bacteria along the CIR, especially between the Onnuri and the other fields. The details of population structure of the CIR organisms will be discussed with the genetic variation of host mussels below.

At the intra-ridge scale, the relatively distant vent habitats along the slow-spreading ridge seem to decrease the connectivity of vent organisms. Geological modeling based on the global survey dataset for the past 35 years estimated the frequency of vent field in the fast-spreading EPR up to three- or fourfold higher than the slow-spreading CIR (vent frequency per 100 km = 0.76–1.33 for CIR and 1.79–3.91 for EPR) (Beaulieu et al., 2015). The geological model and the present results predict that isolation by distance (IBD) becomes apparent in a slow-spreading ridge along with ridge-axis offsets by transform faults. For instance, the genetic variation of symbiotic bacteria observed in the slow-spreading MAR supports this hypothesis. The composition of symbiont assemblages of bathymodioline mussel, I. modiolaeformis, exhibited a statistically significant correlation with geographic distance (Laming et al., 2015).

A recent study has proved that about 13% of total genetic variation on symbionts of bathymodioline mussels distributed along the northern MAR is determined by geographical distance alone (Ücker et al., 2021). If the factor of geographic distance in that study was combined with other factors of environment and host species, the accountability for the total genetic variation increased by about 45%. Therefore, all findings above indicate that we should consider various factors, including the interaction with the host and natural selection related to the chemical and physical environment, to comprehensively understand the population divergence process of symbiotic bacteria.

On a regional scale, the present study examined the genetic diversity and connectivity of the host mussels, B. marisindicus, and their symbiotic bacteria along the CIR, including the recently discovered Onnuri field. Hosts’ genetic diversity was similar among vent fields and significantly higher at the Onnuri. Similarly, the genetic diversity of symbionts was higher at Onnuri and Solitaire than at Edmond and Kairei. To our interests, a similar pattern of population genetic variation was detected in the symbiotic bacteria of scaly-foot snail, Chrysomallon squamiferum, between the Solitaire and Kairei vent fields (Lan et al., 2022). The Edmond mussel symbionts were genetically homogeneous much more than other populations. The genetic diversity of symbionts with the free-living stage is associated with the environmental pool (Ansorge et al., 2019). Frequent habitat turnover by tectonic activity reduces genetic diversity by recurrent bottleneck and founder effect (Wade and Mccauley, 1988; Mccauley, 1991). The loss of genetic diversity of the Edmond mussel symbiont could be attributable to recurrent extinction/recolonization events of bacterial communities or small population size in the first place due to potentially low habitability. In this regard, it is notable that there is circumstantial evidence that the mussel patches in Edmond were very scarce (personal observation during 2001 exploration, which is reported in Van Dover et al., 2001). Additionally, the difference in depths of vent habitats is also a potential factor for the connectivity between vents’ low genetic diversity of the Edmond symbionts. Table 1 shows that the depth difference between Edmond and the two nearby adjacent vents, Solitaire and Kairei, ranges from 650 to 800 m. The oceanic currents flowing at different depths may reduce the gene flow between vents with different depths, increasing geographical isolation by disruption of dispersal (Fujio and Imasato, 1991; Young et al., 2008). The neutrality tests with Tajima’s D and Fu’s Fs statistics refute any significant demographic changes in the four symbiont populations in the CIR, which indicates that the symbiont populations from every vent field including Edmond are in selectively neutral and/or in equilibrium of population size ( Supplementary Table S6 ).

The AMOVA results indicated that symbionts were genetically heterogeneous at the intra-host level and homogeneous at the inter-host level within the same vent field regardless of the examined genes ( Supplementary Table S9 ). Wentrup et al. (2014) suggested two hypotheses regarding the source of symbionts during the development of host mussels: 1) continuous acquisition from the environment and 2) self-infection from own older gill tissues. However, these two alternatives may not be mutually exclusive and could be mixed. Several studies supported the self-infection scenario based on the high genetic variation of symbionts at the inter-host level even within the same fields (Ho et al., 2017; Picazo et al., 2019). On the other hand, the homogeneity of symbionts at the inter-host level observed in B. marisindicus is similar to the pattern of genetic variation observed from the genomic data of the sulfur-oxidizing symbionts of Bathymodiolus mussels in the MAR (Ansorge et al., 2019). These findings suggest that symbiont populations may be intermixed among adjacent host individuals within the same vents by repeated release to and uptake from the environment.

The multi-locus genetic data of symbionts revealed a relatively high level of genetic differentiation among three populations in the Onnuri, Solitaire, and Kairei vent fields ( Supplementary Table S7 ). In contrast, the geographically adjacent Edmond and Kairei symbionts showed relatively low differentiation. The present IM analyses revealed the population divergence order and times from a common ancestral population to the extant three populations representing the three ridge segments of the CIR. To our interest, the population tree topology of symbiotic bacteria was identical to that of the host mussels, although the time depths are different. The northernmost Onnuri mussel populations were phylogenetically distant from the Solitaire and Kairei populations. Meanwhile, the population divergence topologies of host mussel and symbiotic bacteria are consistent with the biogeographic province of vent fauna of the Indian Ocean (Zhou et al., 2022). To our interest, Zhou et al. (2022) reported a summary of faunal similarity in a dendrogram in which a subset of CIR populations reveals an identical tree topology with the present mussel hosts and symbiont populations. Remarkably, the faunal similarity was largely explainable by solely along-ridge distances between vents, about 47.5%. Another relevant finding came out from a phylogeographic study of the scaly-foot snails of the Indian Ocean (Lan et al., 2022), which showed a relatively close relationship between the Solitaire and Kairei symbiont populations compared to the distant populations in the Southwest Indian Ridge (SWIR) and the northernmost Carlsberg ridge. These two studies again emphasizes the strong effect of geographical distance on the divergence of vent invertebrate hosts and their symbiotic bacteria in the Indian Ocean ridges with ultra- to slow-spreading rates. Finally, the consistent pattern of the present study and the previous biogeographic studies altogether indicate the influence of common historical and physical constraints on the population divergence processes of vent invertebrate hosts and their symbionts in the CIR.

According to geochronological analysis, the age of the earliest measurable hydrothermal activity in the Kairei vent field was estimated to be approximately 95 Kya (Wang et al., 2012). From the IM analysis, the host populations of Kairei and Solitaire have undergone the divergence of approximately 75 Kya (95% highest posterior density (HPD) interval: 13–904 Kya). Despite the broad interval of HPD, the highest value is considerably close to the age of the first vent activity measured by isotope, which indicates that the tectonic activity in Kairei may be connected to the divergence of an ancestral mussel population of the Kairei and Solitaire populations. In the present IMa analysis, we could not resolve narrowly how earlier the Onnuri population diverged from the most recent common ancestral population than the subsequent splitting of the Solitaire and Kairei populations. Due to the broad intervals of HPD for the splitting times, the comparison of the divergence between the symbionts and host mussels is difficult. Nevertheless, the coincident pattern of geographical connectivity from the spatially and taxonomically more broad-scaled biogeographic studies suggests that the Onnuri mussel and symbiont populations have diverged substantially from the other populations for a long time.

Overall, both the symbiont and host populations are common in that genetic differentiation increases with geographical distance between habitats along the CIR. The IM analysis provided another insight into migration patterns among populations after divergence of an ancestral population to the descending two populations. We could not detect any significant gene flow between the branches of ancestral populations in both host and symbiont populations. We could only detect significant gene flow in some terminal branches, which implies a sort of secondary contact after divergence in both host and symbiont for a while. The secondary contact after isolation is often observed in deep-sea vent mussels, such as the secondary integration between B. azoricus and B. puteoserpentis in MAR (Breusing et al., 2017) and between B. thermophilus and B. antarcticus in the eastern Pacific (Johnson et al., 2013). In the IM analysis of symbionts, the northward gene flow from the southern two symbiont populations is consistent with the direction of the deep oceanic current at an approximately 2,500-m depth, which was estimated based on geostrophic circulation models (Reid, 2003). The northward migration of symbionts from Kairei to Solitaire was the highest (Nm = 1), indicating that the Kairei population provides genetic variation to Solitaire. The degree of nucleotide diversity at each CIR population agrees well with this convey of genetic variation ( Supplement Figure S1 ), showing a higher level in the receiving populations, Onnuri and Solitaire. On the other hand, the southward gene flow from the Onnuri to Solitaire vents was also detected in the host mussel. It is difficult to interpret this gene flow from Onnuri to Solitaire due to a lack of related environmental information in this area. Considering the daily and seasonal fluctuations of deep oceanic currents, long-distance dispersal of hosts and symbionts may be possible (Murty et al., 2001). Despite the similarities in the phylogenetic tree topologies, the different migration patterns between the hosts and symbionts might be related to different life-history traits and dispersal capabilities of them. The host populations exhibited a higher level of southward gene flow from the Onnuri to Solitaire vents (2Nm >1), which was much more significant than the estimates of the symbiont populations. Bathymodiolin mussels have long-distance dispersal capability due to the long duration of their planktonic larval stage (Lutz et al., 1980; Arellano and Young, 2009; Breusing et al., 2016) and vertical dispersal (Arellano et al., 2014). Therefore, the different patterns and extent of gene flow between the symbionts and hosts indicate that they disperse among vent localities independently.