Hidden Historical Habitat-Linked Population Divergence and Contemporary Gene Flow of a Deep-Sea Patellogastropod Limpet

Abstract Hydrothermal vents and hydrocarbon seeps in the deep ocean are rare oases fueled by chemosynthesis. Biological communities inhabiting these ecosystems are often distributed in widely separated habitats, raising intriguing questions on how these organisms achieve connectivity and whether habitat types facilitate intraspecific divergence. The deep-sea patellogastropod limpet Bathyacmaea nipponica that colonizes both vents and seeps across ∼2,400 km in the Northwest Pacific is a feasible model to answer these questions. We analyzed 123 individuals from four vents and three seeps using a comprehensive method incorporating population genomics and physical ocean modeling. Genome survey sequencing and genotyping-by-sequencing resulted in 9,838 single-nucleotide polymorphisms for population genomic analyses. Genetic divergence and demographic analyses revealed four habitat-linked (i.e., three seep and one vent) genetic groups, with the vent genetic group established via the opportunistic invasion of a few limpet larvae from a nearby seep genetic group. TreeMix analysis uncovered three historical seep-to-vent migration events. ADMIXTURE and divMigrate analyses elucidated weak contemporary gene flow from a seep genetic group to the vent genetic group. Physical ocean modeling underlined the potential roles of seafloor topography and ocean currents in shaping the genetic connectivity, contemporary migration, and local hybridization of these deep-sea limpets. Our study highlighted the power of integrating genomic and oceanographic approaches in deciphering the demography and diversification of deep-sea organisms. Given the increasing anthropogenic activities (e.g., mining and gas hydrate extraction) affecting the deep ocean, our results have implications for the conservation of deep-sea biodiversity and establishment of marine protected areas.


Introduction
The discoveries of deep-sea hydrothermal vents in the late 1970s and hydrocarbon seeps in the early 1980s have significantly changed our understanding of how life has evolved on Earth. Unlike shallow-water ecosystems that are mainly driven by photosynthesis, deep-sea vent and seep ecosystems distributed in tectonically active areas and along continental margins are primarily supported by chemosynthesis, with many species strictly endemic to either vent fields or seep areas (Van Dover et al. 2002). Nevertheless, a number of species are capable of thriving in both habitats (Wolff 2005;Watanabe et al. 2010), raising intriguing questions on how their populations achieve connectivity between widely separated habitats and whether differences in habitat types facilitate intraspecific divergence. Studying species that inhabit both vents and seeps may allow us to address these issues by shedding light on their historical colonization, migration dynamics, and adaptive evolution (Vrijenhoek 2010;Baco et al. 2016). Such information will in turn contribute to the evaluation of resilience of vent-and seep-associated macrofauna after anthropogenic disturbances, management of deep-sea resources, and designation of marine-protected areas (Baco et al. 2016;Breusing et al. 2016;Levin et al. 2020).
With over 90 hydrothermal vents and 70 hydrocarbon seeps discovered, and over 20% macrobenthic species shared across these two types of habitats, the Northwest Pacific has been recognized as an evolutionary hotspot for chemosynthesis-based biota Beaulieu and Szafra nski 2020). Previous population genetic and genomic studies of macrobenthos dwelling in both hydrothermal vents and hydrocarbon seeps in this region have mainly focused on the bathymodioline mussel Gigantidas platifrons and the galatheoid squat lobster Shinkaia crosnieri. Bathymodioline mussels, at least the species for which the dispersal mode has been studied, produce planktotrophic larvae that can disperse a long distance in the upper ocean currents (Arellano et al. 2014;Laming et al. 2018). A population genomic study based on genome-wide single-nucleotide polymorphisms (SNPs) revealed two semiisolated lineages of G. platifrons in the Northwest Pacific, with one mainly inhabiting the Jiaolong Ridge seep in the marginal South China Sea and the other distributed in vent fields in the Okinawa Trough and the Off Hatsushima seep in the Sagami Bay (Xu et al. 2018). Moreover, a minor genetic subdivision was unveiled between its local populations in the southern Okinawa Trough and those in the middle Okinawa Trough and the Sagami Bay. These results reinforced the Luzon Strait and the intra-trough grabens between the southern and the middle Okinawa Trough as potential geographic barriers that may have accelerated the local adaption and genetic differentiation of G. platifrons (Xu et al. 2018). Unlike bathymodioline mussels, the galatheoid squat lobster S. crosnieri produces lecithotrophic larvae that are expected to mainly drift in the middle to the deeper water layers (Miyake et al. 2010;Xiao et al. 2020). Population genetic and genomic studies based on gene markers, genome-wide SNPs, and transcriptome-wide SNPs consistently disclosed a clear genetic division between its seep population in the South China Sea and vent populations in the Okinawa Trough, also highlighting the barrier effect of the Luzon Strait and the environmental adaptabilities of S. crosnieri to the local conditions (Shen et al. 2016;Yang et al. 2016;Cheng et al. 2020;Xiao et al. 2020) . These studies have broadened our knowledge of the population connectivity and genetic divergence of two representative macrobenthos in vent and seep ecosystems. Nonetheless, how other vent and seep dominant macrobenthos achieve connectivity amongst geographically and geochemically distinct habitats, and whether they display a habitat-linked demographic history have not been critically examined.
The patellogastropod limpet Bathyacmaea nipponica is amongst the most widely distributed macrobenthos across chemosynthesis-based ecosystems in the Northwest Pacific, and can therefore serve as a model organism to fill such knowledge gaps. Bathyacmaea nipponica mainly grazes bacterial film on various hard substrates for nutrition, and this species has been considered to produce pelagic-dispersing lecithotrophic larvae as is typical in patellogastropod limpets (Chen et al. 2019;Ponder et al. 2020). To date, B. nipponica has been found in vent fields in the Okinawa Trough and hydrocarbon seeps in the South China Sea, the Ryukyu Arc, and the Sagami Bay ( fig. 1A), with depths ranging between 625 to 1,684 m (http://www.godac.jamstec.go.jp/bio-sample/ index_e.html, last accessed January 15, 2021) Zhao et al. 2020). Here, we applied a comprehensive approach incorporating population genomics and physical ocean modeling to examine the demographic history and population connectivity of B. nipponica. Through this study, we aimed to 1) understand the habitat colonization and evolutionary divergence of its vent and seep populations; 2) test whether these populations exhibit a habitat-linked genetic differentiation; and 3) investigate the migration dynamics and gene flow of B. nipponica under the intricate interactions amongst local habitats, seafloor topography, and ocean currents.

Assembly of B. nipponica Survey Genome and Detection of Genome-Wide SNPs
The Illumina whole-genome sequencing generated 111.2 Gb of raw reads. After removing low-quality reads, a total of 102.9 Gb reads were retained for genome assembly using Platanus, which produced 2,733,936 scaffolds with a total size of 950.  1F). In particular, the population divergence pattern displayed along the first two eigenvectors (i.e., G1 and G2; G3; G4) agreed with the ADMIXTURE results at K ¼ 3 ( fig. 1D). Additionally, the G1 and G2 genetic group was separated into two distinct genetic groups G1 and G2 along the third eigenvector, with G1 containing all individuals from the Jiaolong Ridge seep, and G2 including all individuals from the Off Hatsushima seep along with two individuals from the Hatoma Knoll vent and three individuals from the Izena Hole Field vent in the Okinawa Trough. These results were consistent with the reconstructed phylogenetic network ( fig. 1C) and the genetic structure revealed by ADMIXTURE (K ¼ 4; fig. 1D).

Population Genetic Statistics
The Okinawa Trough vent genetic group had the greatest number of polymorphic nucleotide sites (6,554; 66.6%) and private SNPs (2,886; 29.3%), as well as the highest values of the observed heterozygosity (H obs ), nucleotide diversity (p), and inbreeding coefficient (F IS ), compared with the other three seep genetic groups (table 1 and supplementary table S5, Supplementary Material online).

Migration Dynamics
TreeMix analysis uncovered four migration events (m ¼ 4; fig.  2A). Amongst them, three seep-to-vent migration events were statistically significant, including two from the phylogenetical ancestry of the Jiaolong Ridge and the Off Hatsushima seep populations to the Hatoma Knoll vent population in the southern Okinawa Trough (P < 0.00001) and to the Izena Hole Field vent population in the middle Okinawa Trough (P < 0.00001), as well as one from the Kuroshima Knoll seep population to the phylogenetical ancestry of the southern Okinawa Trough vent populations (P < 0.005). divMigrate analysis unveiled intense gene flow between the Jiaolong Ridge seep genetic group and the Off Hatsushima seep genetic group, weak gene flow from the Off Hatsushima seep genetic group to the Okinawa Trough vent genetic group, and very limited gene flow between the Kuroshima Knoll seep genetic group and the other two seep genetic groups in the Jiaolong Ridge and the Off Hatsushima ( fig. 2B). Extensive gene flow was also observed between the vent local population pairs within the Okinawa Trough, with an overall stronger gene flow from the shallower to the deeper vent fields, especially for the vent local populations within the same region (i.e., within either the southern or the middle Okinawa Trough; fig. 2C).

Demographic History
Posterior probabilities evaluated based on the emphdirect and logistic regression methods consistently indicated that scenario 1 provided the most likely demographic history of B. nipponica  3D). This result implied that demographic parameters of scenario 1 estimated herein (table 2) fitted the data well (Cornuet et al. 2014).
Under this demographic scenario, genetic divergence between the deeper seep (i.e., Jiaolong Ridge and Off Hatsushima) genetic group and the shallower seep (i.e., Kuroshima Knoll) genetic group occurred $24,400 generations ago; opportunistic invasion by a small number of limpet larvae ($650) from the Kuroshima Knoll seep genetic group to vent fields in the Okinawa Trough occurred $11,800 generations ago; establishment of the Okinawa Trough vent genetic group occurred $6,410 generations ago, and further genetic separation of the deeper seep genetic group into the Jiaolong Ridge and the Off Hatsushima seep genetic groups occurred $1,220 generations ago (table 2).
Amongst the four genetic groups, the Jiaolong Ridge seep genetic group had the smallest effective population size of $10,700, whereas the Okinawa Trough vent genetic group had the largest effective population size of $98,200 (table 2).

Physical Ocean Modeling
Seasonal-mean lateral ocean currents in the study region at selected depths during the five model years (2011)(2012)(2013)(2014)(2015) are illustrated in supplementary figure S1, Supplementary Material online. The time-mean lateral velocity along (or within) the Okinawa Trough decayed rapidly in the vertical direction, from $1 m/s at the sea surface to $0.1 m/s below 1,000 m depth, regardless of the seasons considered. The vertical decline of the lateral flow was also reflected in the nu-    5N and O and supplementary movie 12). Furthermore, no numerical particles released at 800 m depth from the Off Hatsushima seep coordinates on July 1, 2011 were found to drift into the Okinawa Trough (supplementary movie 11). In all cases, numerical particles released from the Off Hatsushima seep  Hatsushima seep genetic groups were inferred to be hybrid descendants between individuals from these two genetic groups. Therefore, these individuals were excluded from the demographic analysis to avoid potential bias.

Discussion
In the present study, we examined the demographic history and population connectivity of the deep-sea patellogastropod limpet B. nipponica by using an integrated population genomic and physical ocean modeling approach. Our results illustrated how the vent and seep populations of B. nipponica have diversified and connected under the complex interactions amongst local habitats, seafloor topography, and ocean currents.
The seven studied populations of B. nipponica formed four genetic groups that were primarily differentiated by habitat types, including three seep (i.e., Jiaolong Ridge, Off Hatsushima, and Kuroshima Knoll) and one vent (i.e., Okinawa Trough) genetic groups. These results were contrasted with those of a previous study based only on the nuclear histone H3 and the mitochondrial cytochrome c oxidase subunit I (cox1) genes, in which B. nipponica in the   (Fujikura et al. 2000). Nevertheless, we cannot rule out the possibility that other factors (e.g., local currents and geofluid chemistry) may have contributed to this as well. Located in the fore-arc regions, the Off Hatsushima seep in the Sagami Bay and the Kuroshima Knoll seep in the Ryukyu Arc exhibit a high level of fluid seepage from the seafloor (Masuzawa et al. 1992;Takeuchi et al. 2007). In comparison, the Jiaolong Ridge seep in the South China Sea has displayed a weak supply of methane over the last two thousand years (ka) (Feng and Chen 2015). The relatively small pairwise F ST values ( fig. 1B and supplementary table S4, Supplementary Material online), high genetic similarities ( fig. 1C, D, and F), and extensive gene flow ( fig. 2B) between the Jiaolong Ridge and the Off Hatsushima seep genetic groups were in accordance with the results of DIYABC analysis, which revealed a very recent divergence ($1,220 generations ago) of these two deeper seep genetic groups (table 2). The accurate generation time of Bathyacmaea remains unknown. With the assumption that B. nipponica has a generation time of one to two years (previous live-rearing experiments in aquaria under atmospheric   (Feng and Chen 2015). Such diminished geofluid released from the Jiaolong Ridge seep resulted in the vanishing of chemosynthetic communities on the ridge flanks (Feng and Chen 2015), which may also be responsible for the small effective population size of B. nipponica inhabiting there (  fig. 2A) and the geographical proximity of these two regions further indicated that the Okinawa Trough vent genetic group may have been established following the historical invasion of limpet larvae from the Kuroshima Knoll seep. This inference was also supported by the results of DIYABC analysis, which unveiled that a small number of limpet larvae ($650) from the Kuroshima Knoll seep might have invaded the Okinawa Trough $11,800 generations ago, adapted to the local environments, and thereafter formed the Okinawa Trough vent genetic group $6,410 generations ago (table 2). Assuming a generation time of one to two years for B. nipponica is appropriate (Chen et al. 2019), these patterns of population divergence would fit in a time window (i.e., larval invasion from the Kuroshima Knoll seep to the Okinawa Trough: $23.6 to 11.8 ka; population establishment in the Okinawa Trough vents: $12.8 to 6.4 ka) corresponding to the possible historical variation in strength and main course of the Kuroshio Current. Previous studies suggested that potentially due to the decrease of sea level and/or emergence of a Ryukyu-Taiwan land bridge, the Kuroshio mainstream had shifted toward the east of the Ryukyu Arc during the last glacial maximum and then swung back to the Okinawa Trough $14 to 7 ka (Jian et al. 2000;Ujii e et al. 2003;Lim et al. 2017). Such a shifting of the Kuroshio Current may have thus facilitated the historical migration of B. nipponica larvae from the Kuroshima Knoll seep into the Okinawa Trough, and afterwards accounted for their successful population expansion and local adaptation in vent fields within the trough. Nevertheless, it should be noted that our estimation of the demographic events suffered from an uncertainty of the limpet's generation time. Determining its accurate generation time via laboratory larval culture experiments, although challenging for deep-sea organisms, is desired in the future to verify the inferences raised here.
Population genetic analyses of F ST estimation, phylogenetic network reconstruction, ADMIXUTRE (K ¼ 3 and K ¼ 4), and PCA concordantly revealed that all the vent local populations in the Okinawa Trough belonged to one metapopulation (fig. 1B-D and F and supplementary 2C). Physical ocean modeling uncovered a significant decline in the time-mean lateral velocity over depths within the Okinawa Trough (supplementary fig. S1, Supplementary Material online). Moreover, numerical particle release experiments demonstrated that water transport between the Okinawa Trough and the open Northwest Pacific was topographically constrained, especially at and below 800 m depth ( fig. 4J-R and supplementary movies 4-6). These different lines of evidence implied that the larvae of B. nipponica may mainly disperse along the intermediate to the deeper ocean currents, which is similar to the larvae of S. crosnieri (Miyake et al. 2010;Xiao et al. 2020). Therefore, the overall concave topography of the Okinawa Trough may have trapped the dispersal of B. nipponica larvae (Mitarai et al. 2016), and afterwards led to a high level of genetic homogeneity (e.g., small pairwise F ST values between the vent local population pairs), a high degree of inbreeding (i.e., the highest F IS value amongst the four genetic groups), and extensive contemporary gene flow with a direction from the shallower to the deeper vents for individuals inhabiting the Okinawa Trough. Such a larval dispersal potential of B. nipponica would also help to explain its different population segregation patterns in comparison with those of the bathymodioline mussel G. platifrons. This mussel has been considered to produce planktotrophic larvae that may mainly disperse in the upper currents, with the local populations in vent fields in the Okinawa Trough and the Off Hatsushima seep in the Sagami Bay sharing high genetic similarities (Xu et al. 2018). Moreover, the largest number of private SNPs identified in the Okinawa Trough vent genetic group implied the adaptive capabilities of B. nipponica to the local vent environments, which may have also contributed to the highest heterozygosity of this genetic group compared with all the other genetic groups.
In addition, occurrence of weak contemporary gene flow from the Off Hatsushima seep to the Okinawa Trough vent genetic groups of B. nipponica was accordantly unveiled by the phylogenetic network, ADMIXTURE results (K ¼ 4), and PCA plots; that is, five Okinawa Trough vent individuals inferred to be immigrants from the Off Hatsushima seep, and two Okinawa Trough vent individuals inferred to be hybrid descendants between individuals from the Okinawa Trough vent and the Off Hatsushima seep genetic groups ( fig. 1C, D, and F). Consistent with the results of these population genomic analyses, physical ocean modeling deciphered the intrusion of some water from the Off Hatsushima seep coordinates into the Okinawa Trough at 500 and 800 m depths mainly via the channel east of Taiwan and the Kerama Gap, although such invasion was largely restricted at 800 m depth due to the Xu et al. . doi:10.1093/molbev/msab278 MBE topographic constraints of the Okinawa Trough ( fig. 5 and  supplementary movies 7-12). This water body belongs to the North Pacific Intermediate Water (Bostock et al. 2010;Nakamura et al. 2013). Hence, transport of the North Pacific Intermediate Water might have facilitated the contemporary dispersal of a few B. nipponica larvae from the Off Hatsushima seep into the Okinawa Trough. These larvae then settled in different vent fields in the Okinawa Trough, and subsequently resulted in hybridization between individuals belonging to the Off Hatsushima seep and the Okinawa Trough vent genetic groups.
At present, deep-sea vent and seep ecosystems are facing an increasing threat from ongoing or upcoming anthropogenic disturbances, most noticeably natural resources exploitation (e.g., minerals and ores in vents and methane hydrates, oil, and gas in seeps) (Ramirez-Llodra et al. 2011;Levin et al. 2020). These disturbances have been predicted to cause adverse effects on the entire marine ecosystems, leading to degradation and loss of habitats, reduction in biodiversity, and alteration of community structure and ecological function (Ramirez-Llodra et al. 2011;Levin et al. 2020). Extraction of deep-sea ores has already been tested in hydrothermal vents in the Okinawa Trough (Okamoto et al. 2019), making this area potentially one of the world's first deep-sea mining sites in the near future. Therefore, the results generated herein would have practical implications for the conservation of vent-and seep-associated organisms, especially for those exhibiting similar dispersal potential as B. nipponica (Mills et al. 2009): Although the same species dwells in both vents and seeps, existence of habitat-linked genetic structure implies that quick recolonization from one habitat to another cannot be assumed. As a consequence, individual sites within the same habitat and/or geographic region, such as some vent fields in the Okinawa Trough, should be identified and assigned as protected areas to ensure a high resilience is maintained within the local population groups. Furthermore, recent drilling campaigns have deepened our understanding of the distribution and storage of gas hydrate in seep ecosystems in the South China Sea, at the same time, drawn our attention to monitor and protect the chemosynthesis-based communities in this region . Detectable genetic divergence between the South China Sea and the Northwest Pacific populations underlines the barrier effect of the Luzon Strait on larval dispersal of B. nipponica, as has been demonstrated in other deep-sea organisms inhabiting these two areas (e.g., Xu et al. 2018;Xiao et al. 2020). Due to the marginal geographic feature of the South China Sea and the genetic distinctness of its seep communities, we therefore suggest that seep ecosystems in the South China Sea should also be given a high priority of biodiversity conservation when designating deep-sea reserves and marine protected areas in the global ocean (Baco et al. 2016;Levin et al. 2020).

Conclusions
Using an integrated approach combining population genomics and physical ocean modeling, four distinct genetic groups with a habitat-linked pattern were disclosed for the deep-sea patellogastropod limpet B. nipponica in the Northwest Pacific, including three seep genetic groups inhabiting hydrocarbon seeps in the Jiaolong Ridge, the Kuroshima Knoll, and the Off Hatsushima, as well as one vent genetic group dwelling in four hydrothermal vents in the Okinawa Trough. The assumed lecithotrophic larvae of B. nipponica were considered to mainly disperse along the intermediate to the deeper ocean currents. The most likely demographic history of B. nipponica was deciphered to have occurred in the following sequences: 1) an initial split into two seep genetic groups, one deeper (i.e., Jiaolong Ridge and Off Hatsushima) and one shallower (i.e., Kuroshima Knoll); 2) a small number of limpet larvae from the shallower seep genetic group invaded vent fields in the Okinawa Trough, adapted to the local environments, and thereafter established a distinct vent genetic group; and 3) two deeper seep populations diverged very recently into two distinct genetic groups with detectable genetic differentiation. Assuming a generation time of one to two years for B. nipponica, formation of its vent genetic group may have been linked to the historical shifting of the Kuroshio Current. Additionally, diversification of the two deeper seep genetic groups might be associated with the recent reduction of seepage intensity in the Jiaolong Ridge seep along with the barrier effect of the Luzon Strait. Physical ocean modeling results elucidated the seafloor topography of the Okinawa Trough may have contributed to the genetic homogeneity of B. nipponica among different vent local populations in this region, and the North Pacific Intermediate Water might have facilitated the opportunistic migration of a few limpet larvae from the Off Hatsushima seep in the Sagami Bay into the Okinawa Trough, and thus led to their hybridization in the Okinawa Trough vents. Altogether, this study has enhanced our knowledge on the historical population divergence and contemporary gene flow of deep-sea organisms inhabiting hydrothermal vents and hydrocarbon seeps under the intricate interactions amongst local habitats, seafloor topography, and ocean currents. More importantly, these results could serve as the scientific basis for a better conservation of chemosynthesis-based ecosystems and an effective establishment of marine protected areas. MBE preserved in 99.5% ethanol. Genomic DNA was extracted from the foot of each individual using the cetyltrimethylammonium bromide method (Stewart and Via 1993). Purity, integrity, and quality of DNA were examined using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE), 1.0% agarose gel electrophoresis, and a Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA), respectively.

Genome Survey Sequencing and Assembly
To provide a reference genome for SNP identification, genome survey sequencing was performed for an individual of B. nipponica (i.e., JR-5) in Novogene (Beijing, China). A DNA sequencing library with an insert size of 350 bp was constructed using the NEBNext DNA Library Prep Kit and paired-end sequenced on an Illumina NovaSeq 6000 sequencer with a read length of 150 bp. Raw reads were filtered using Trimmomatic v.0.38 (Bolger et al. 2014) to remove adapters and low-quality reads with the following settings: ILLUMINACLIP: Adapters.fas: 2:30:10, LEADING ¼ 10, TRAILING ¼ 10, SLIDINGWINDOW ¼ 4:15, and MINLEN ¼ 40. Obtained high-quality reads were assembled using Platanus v.1.2.4 (Kajitani et al. 2014) under the default parameters. Assembly statistics were evaluated using the Perl script assemblathon_stats.pl (Bradnam et al. 2013), and assembly completeness was accessed using BUSCO v.3.0.2 (Simão et al. 2015) based on the metazoa_odb9 database.

Genotyping by Sequencing and Genome-Wide SNP Identification
Genotyping by sequencing was performed in Novogene (Beijing, China). In brief, genomic DNA from each individual was digested using the restriction enzymes MseI and MspI. Pooled DNA libraries with insert sizes ranging from 240 to 265 bp were paired-end sequenced on an Illumina NovaSeq 6000 sequencer with a read length of 144 bp. Raw reads of each individual were filtered to remove adapters and poorquality reads using fastp v.0.20.0 ) under the following settings: -detect_adapter_for_pe and -trim_po-ly_x, set -qualified_quality_phred to 15, -unqualified_per-cent_limit to 30, -n_base_limit to 0, and set -max_len1, -max_len2, and -length_required all to 130. Obtained highquality reads of each individual were mapped to the survey genome of B. nipponica using the MEM algorithm implemented in BWA v.0.7.17 ). The resultant .sam files were individually converted to .bam files and sorted using SAMtools v.1.9 ). Duplicate products of polymerase chain reaction and reads mapped to different positions of the survey genome were removed using Sambamba v.0.7.1 (Tarasov et al. 2015). Calling of SNPs was performed under a combination of the mpileup (-q 10 -Q 15 -skip-indels) and the call options (-multiallelic-caller -variants-only) implemented in BCFtools v.1.9 (Li 2011). Filtration was conducted using VCFtools v.0.1.16 (Danecek et al. 2011) to retain SNPs that meet the following criteria: 1) with two alleles (-min-alleles 2 -max-alleles 2); 2) with a sequencing coverage depth between 5 and 200 (-minDP 5 -maxDP 200); 3) with a genotype quality ! 20 (-minGQ 20); 4) presence in ! 80% of all individuals (-max-missing 0.8); 5) with a global minor allele frequency ! 0.02 (-maf 0.02); 6) SNP with an observed heterozygosity (H obs ) < 0.5 calculated using popStats (Skoglund et al. 2015) to avoid inclusion of potential paralogues (Hohenlohe et al. 2011). Based on the statistics of SNP per genomic locus retained after the sixth filtering criterion, all SNPs located in genomic locus with more than nine SNPs were further discarded to reduce potential sequencing or genotyping errors. Formatting of the resultant final SNP data set for downstream analyses was achieved using PGDSpider v.2.1.1.5 Excoffier 2012) andStacks v.2.5 (Rochette et al. 2019).

Population Divergence Analyses and Genetic Statistics Estimation
Genetic differentiation between the local population pairs represented by pairwise F ST was estimated using Arlequin v.3.5.2.2 (Excoffier and Lischer 2010), with 10,000 permutations performed to test for significance. A phylogenetic network was constructed with the NeighborNet method implemented in SplitsTree v.4.15.1 (Huson 1998), based on the normalized distances estimated with the uncorrectedP algorithm and the option of Handle Ambiguous States set to MatchStates. Population genetic structure was investigated using two methods with only one SNP per locus retained to reduce potential bias derived from linkage disequilibrium. A maximum likelihood (ML) estimation method implemented in ADMIXTURE v.1.3.0 (Alexander et al. 2009) was applied to examine population structure and individual ancestry, with the number of genetic groups represented by K predefined from one to seven. The number of genetic groups that best fit the data (i.e., the optimal K) was assessed using the crossvalidation (CV) procedure implemented in ADMIXTURE with the -cv flag. The barplot function in R was used for result visualization with the selected K. The PCA implemented in the R package SNPRelate (Zheng et al. 2012) was applied to further investigate genetic divergence amongst all individuals. Several population genetic statistics, including expected heterozygosity (H exp ), observed heterozygosity (H obs ), nucleotide diversity (p), and inbreeding coefficient (F IS ), were calculated using the POPULATIONS module implemented in Stacks v.2.5 (Rochette et al. 2019).

Migration Dynamics Analyses
Migration events amongst local populations were evaluated using TreeMix v.1.13 (Pickrell and Pritchard 2012). Specifically, an ML tree was constructed based on the SNP data set (only one SNP per locus retained) with migration events added to this tree one at a time to infer migration event models along with the -global plus -se options. Besides, migration patterns between pairs of genetic groups and pairs of the selected local populations were examined using divMigrate-online (Sundqvist et al. 2016) based on the G ST statistic as a measure of genetic differentiation.

Demographic History Reconstruction
Demographic history was reconstructed using the approximate Bayesian computation (ABC) analysis implemented in Xu et al. . doi:10.1093/molbev/msab278 MBE DIYABC v.2.1.0 (Cornuet et al. 2014). According to the results of pairwise F ST estimation, phylogenetic network, population structure, and migration dynamics (details refer to the Results section), three demographic scenarios were hypothesized to infer how the four identified genetic groups of B. nipponica had formed and diversified in the Northwest Pacific ( fig. 3A). Each scenario was predefined by a number of demographic parameters expressed as the number of generations back in time (i.e., T 1 < T 2 < T 3 < T 4 < T A ) with details summarized in supplementary table S8, Supplementary Material online.
In scenario 1, B. nipponica initially diversified into two seep genetic groups at T A . The first seep genetic group inhabited the shallower seep area in the Kuroshima Knoll and formed the Kuroshima Knoll seep genetic group. A few limpet larvae from the Kuroshima Knoll seep genetic group (i.e., N KKP ) invaded and settled in the Okinawa Trough vents by chance at T 3 , which later diverged to form the Okinawa Trough vent genetic group at T 2 . The second seep genetic group inhabited the deeper seep areas, which diverged into two distinct seep genetic groups in the Jiaolong Ridge and the Off Hatsushima seeps at T 1 . In scenario 2, B. nipponica inhabited the deeper seep areas was hypothesized to be ancestral. A few limpet larvae from this ancestral genetic group (i.e., N AP ) incidentally invaded and settled in the shallower seep area in the Kuroshima Knoll at T A , which gradually formed the Kuroshima Knoll seep genetic group at T 4 . Afterward, a few limpet larvae from the Kuroshima Knoll seep genetic group (i.e., N KKP ) occasionally dispersed and settled in the Okinawa Trough vents at T 3 , which thereafter formed the Okinawa Trough vent genetic group at T 2 . More recently, the ancestral deeper seep genetic group further diversified into two distinct seep genetic groups in the Jiaolong Ridge and the Off Hatsushima seeps at T 1 . In scenario 3, B. nipponica inhabited the shallower seep area in the Kuroshima Knoll was speculated to be ancestral. A few limpet larvae from the Kuroshima Knoll seep genetic group (i.e., N KKP1 ) fortuitously invaded the deeper seep areas at T A , which then established the deeper seep genetic group at T 4 . Furthermore, a few other limpet larvae from the Kuroshima Knoll seep genetic group (i.e., N KKP2 ) incidentally entered and settled in the Okinawa Trough vents at T 3 , which afterward formed the Okinawa Trough vent genetic group at T 2 . Lately, the deeper seep genetic group further differentiated to establish two distinct seep genetic groups in the Jiaolong Ridge and the Off Hatsushima seeps at T 1 .
To obtain robust results, a total of 3,000,000 simulations were performed as recommended by DIYABC for all three demographic scenarios. Summary statistics, including proportion of zero values, mean and variance of nonzero values, and mean of complete distribution, were calculated for each genetic group. In addition, summary statistics, including mean and variance of nonzero values of F ST , proportion of zero values and mean of complete distribution Nei's distances, were calculated between pairs of genetic groups. Posterior probabilities and 95% confidence intervals for each scenario were estimated based on the 500 simulated data closest to the observed data using the emphdirect method and 1% of the simulated data closest to the observed data using the logistic regression method. Posterior distribution of each demographic parameter of the best scenario was calculated according to 1% simulated data closest to the observed data with the logistic regression method. Fitness of the parameter priors to the best scenario was illustrated using PCA with the function of Perform Model Checking.

Physical Ocean Modeling
Population genomic analyses conducted in the present study uncovered extensive contemporary gene flow of B. nipponica between its vent local populations in the Okinawa Trough, and weak contemporary gene flow from its Off Hatsushima seep to the Okinawa Trough vent genetic groups (details refer to the Results section). To establish a deeper understanding of these gene flow patterns from an oceanographic aspect, physical ocean modeling data were further analysed to uncover the 1) seasonal-mean lateral ocean current patterns and 2) topographic constraints and fluid dynamics in the study region via two sets of numerical particle release experiments.
Velocity fields of lateral ocean currents were extracted from the HYCOM þ NCODA Global 1/12 Reanalyses (experiment sequence: 53.X) modeling output (https://www. hycom.org, last accessed April 10, 2021) in the study region. This set of ocean modeling data has integrated multiple sources of ocean observational records, including satellite measurements of the sea surface height/temperature and in-situ hydrographic information documented by the XBTs/ Argo floats/moored buoys. Seasonal-mean lateral ocean current patterns at 0, 500, 800, 1,000, and 1,250 m depths were elucidated by averaging the daily velocity fields within each season (spring: March-May; summer: June-August; fall: September-November; winter: December-February) across an ensemble of five model years (2011)(2012)(2013)(2014)(2015). This calculation effectively filtered out high-frequency processes, such as transient mesoscale eddies, which allowed investigation of the present-day large-scale ocean circulations.
To understand the intensive contemporary gene flow between the vent local populations of B. nipponica within the Okinawa Trough (details refer to the Results section), the first set of numerical particle release experiments were performed by uniformly releasing numerical particles within a 25-km radius from the Futagoyama Field vent coordinates (i.e., the westernmost vent location in the Okinawa Trough included herein), which then flew passively with ocean currents for around three years since 2011. Based on the seasonal-mean lateral ocean current patterns in the study region and the migration dynamics of B. nipponica (details refer to the Results and Discussion sections), B. nipponica larvae were considered to mainly disperse in the intermediate to the deeper water similar to those of S. crosnieri (Miyake et al. 2010;Xiao et al. 2020). Therefore, a series of depths ranging from the intermediate to the deeper water were selected for numerical particle release experiments in the Okinawa Trough region to verify such an inference (supplementary table S9, Supplementary Material online).
Previous studies suggested that the North Pacific Intermediate Water can intrude into the Okinawa Trough Population Genomics of a Deep-sea Limpet . doi:10.1093/molbev/msab278 MBE mainly via the channel east of Taiwan and the Kerama Gap (Bostock et al. 2010;Nakamura et al. 2013), which may explain the weak contemporary gene flow of B. nipponica from the Off Hatsushima seep in the Sagami Bay to vent fields in the Okinawa Trough (details refer to the Results section). To examine this hypothesis, the second set of numerical particle release experiments were performed by uniformly releasing numerical particles within a 25-km radius from the Off Hatsushima seep coordinates at 500 and 800 m depths. These two depths were selected based on the seasonalmean lateral ocean current patterns in the study region (details refer to the Results section), topographic constraints of the Okinawa Trough inferred from the results of the first set of numerical particle release experiments (details refer to the Results section), and depth range of the North Pacific Intermediate Water (Bostock et al. 2010;Nakamura et al. 2013). Numerical particle release experiments unveiled the locations for numerical particles released from the Off Hatsushima seep coordinates to invade the Okinawa Trough varied with the release dates. Therefore, experiments performed on three representative release dates spanning 2011 at two selected water depths till the end of 2015 are presented herein to illustrate the intrusion of the North Pacific Intermediate Water into the Okinawa Trough (supplementary table S9, Supplementary Material online).
The two sets of numerical particle release experiments were analyzed based on the horizontal and time-varying ocean current velocity. To attain statistically robust results, a broader domain (15 N-40 N Â 115 E-180 E) than the study region was selected and over 10,000 numerical particles were released for each set of experiments. This approach accommodated long-range trajectories of particles following the Pacific Ocean currents (e.g., the Kuroshio Extension). Trajectories of these particles were calculated using a stepadapting fourth/fifth-order Runge-Kutta method with interpolations done using a cubit method (Wang et al. 2016). Seasonal-mean lateral ocean currents and distribution patterns of numerical particles were visualized using MATLAB v. R2018a (https://www.mathworks.com/products/matlab. html), aiming to reflect the topographic constraints and fluid dynamics of the study region. No data are available on the swimming capability of B. nipponica larvae. Therefore, the present study is not applicable for determining the exact travel time of B. nipponica larvae along ocean currents, such as from the Off Hatsushima seep in the Sagami Bay to vent fields in the Okinawa Trough.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.