Population subdivision of hydrothermal vent polychaete Alvinella pompejana across equatorial and Easter Microplate boundaries

The Equator and Easter Microplate regions of the eastern Pacific Ocean exhibit geomorphological and hydrological features that create barriers to dispersal for a number of animals associated with deep-sea hydrothermal vent habitats. This study examined effects of these boundaries on geographical subdivision of the vent polychaete Alvinella pompejana. DNA sequences from one mitochondrial and eleven nuclear genes were examined in samples collected from ten vent localities that comprise the species’ known range from 23°N latitude on the East Pacific Rise to 38°S latitude on the Pacific Antarctic Ridge. Multi-locus genotypes inferred from these sequences clustered the individual worms into three metapopulation segments — the northern East Pacific Rise (NEPR), southern East Pacific Rise (SEPR), and northeastern Pacific Antarctic Ridge (PAR) — separated by the Equator and Easter Microplate boundaries. Genetic diversity estimators were negatively correlated with tectonic spreading rates. Application of the isolation-with-migration (IMa2) model provided information about divergence times and demographic parameters. The PAR and NEPR metapopulation segments were estimated to have split roughly 4.20 million years ago (Mya) (2.42–33.42 Mya, 95 % highest posterior density, (HPD)), followed by splitting of the SEPR and NEPR segments about 0.79 Mya (0.07–6.67 Mya, 95 % HPD). Estimates of gene flow between the neighboring regions were mostly low (2 Nm < 1). Estimates of effective population size decreased with southern latitudes: NEPR > SEPR > PAR. Highly effective dispersal capabilities allow A. pompejana to overcome the temporal instability and intermittent distribution of active hydrothermal vents in the eastern Pacific Ocean. Consequently, the species exhibits very high levels of genetic diversity compared with many co-distributed vent annelids and mollusks. Nonetheless, its levels of genetic diversity in partially isolated populations are inversely correlated with tectonic spreading rates. As for many other vent taxa, this pioneering colonizer is similarly affected by local rates of habitat turnover and by major dispersal filters associated with the Equator and the Easter Microplate region.


Background
The past 25 years of population genetic studies have revealed a number of physical and biological processes that shape the geographical structure, interpopulation connectivity and genetic diversity of deep-sea hydrothermal vent species (reviewed in [1]). Extrinsic factors, such as the geomorphology of oceanic ridges, deep oceanic currents and the temporal stability of vents, influence the genetic structure of vent species, and intrinsic factors, such as taxon-specific differences in larval development, larval duration, motility and behavior, affect connectivity [1][2][3][4][5][6][7].
The most intensively studied invertebrate animals, in these regards, inhabit the southeastern Pacific ridge systems (Fig. 1), composed of the northern and southern East Pacific Rise (NEPR and SEPR), the Galápagos Rift (GAR) and the northeastern Pacific Antarctic Ridge (PAR). Three metapopulation patterns have emerged from past studies (see Figure 4 in [1]); (1) the metapopulations typically exhibit geographical subdivision involving one or two partitions along the contiguous ridge axes, with relatively low genetic differentiation within metapopulation segments; (2) geographical boundaries of the partitions Fig. 1 Map of the eastern Pacific ridges and Alvinella pompejana sampling locations along the ridge axes. Colored dots indicate sampled vent field locations: blue, green and red dots, sampling sites within the NEPR, SEPR and PAR regions. Gray dots represent records for A. pompejana but were not included in this study. The numbers in parenthesis following the locality names indicate tectonic spreading rate (mm/year). Inset photo of A. pompejana was supplied courtesy of G. W. Rouse, Scripps Institution of Oceanography, UCSD. Photo was taken from sample collections of the Alarcón Rise can vary among the co-distributed species; and (3) the taxon-specific effects of these boundaries as dispersal filters can range from complete isolation and speciation (vicariance) to limited (or no) dispersal impedance. The most consistent boundaries identified to date coincide with the Equatorial region, separating the NEPR + GAR axes from the SEPR axis, and the Easter Microplate region, separating the SEPR and PAR axes.
This study examined geographical population structure and connectivity of the Pompeii worm, Alvinella pompejana Desbruyères and Laubier, (Fig. 1) [8]. With a known distribution spanning 8300 km, this annelid is among the pioneer-species that settle first on newly formed hydrothermal chimneys on the southeastern Pacific ridge systems [9]. Its extraordinary thermal tolerance has attracted the attention of vent researchers for the past two decades [10][11][12][13][14][15]. The worms are covered by a "dense fleece" of epsilon proteobacterial episymbionts that contribute to their nutrition and thermal protection [16,17]. Production of large (~200 μm) lecithotrophic larvae that arrest development in cold abyssal waters allows the worms to disperse great distances and be among the first to colonize nascent hot vents [9]. Its life history and behavioral traits appear to be optimized for exploiting the patchily distributed and highly ephemeral eastern Pacific hydrothermal vents [1]. Nonetheless, we have only rudimentary knowledge about the effects of these traits on the geographical structure and genetic connectivity of A. pompejana metapopulations.
Previous studies of A. pompejana identified distinct metapopulation segments separated by the Equatorial boundary. An examination of mitochondrial cytochromec-oxidase subunit-I (mtCOI) sequences (710 bp) in samples that ranged between 21°N to 32°S latitude identified distinct NEPR and SEPR metapopulation segments, but found no evidence for a distinct segment occupying the PAR axis [18]. A subsequent examination of mtCOI, multi-locus allozymes and four nuclear genes in A. pompejana samples ranging between 13°N and 21°S confirmed the distinct NEPR and SEPR segments [19,20]. To date, however, multi-locus genetic markers have not been examined in samples from the extended range of A. pompejana, 23°N on the EPR to 38°S on the PAR, reported for the first time in this study.
Comprehensive geographical sampling and the application of multi-locus genetic markers have often improved or even contradicted previous inferences about genetic structure and the demographic history of vent species. For example, a multi-locus investigation revealed that SEPR populations of the siboglinid polychaete Tevnia jerichonana exhibited a broad zone of intergradation between distinct metapopulation segments occupying the NEPR and PAR axes [21]. Inferences based on mtCOI evidence alone reached different conclusions [4,18]. Coykendall et al.'s [22] multi-locus study of the siboglinid polychaete Riftia pachyptila did not corroborate an earlier mtCOI study that concluded SEPR and PAR populations were partially isolated across the Easter Microplate boundary [18]. Multi-locus data also revealed a hybrid zone at the Easter Microplate boundary [23,24] that was not recognized in more limited samples of the vent mussels Bathymodiolus thermophilus and B. antarcticus [2].
Computer simulations revealed that inferences about population subdivision and isolation-by-distance are limited by the number of independent gene loci examined and the completeness of population sampling [4]. Sampling gaps can create false evidence for subdivision and examinations of mitochondrial DNA alone occasionally provide signals for population subdivision that are discordant with evidence provided by independent nuclear genes (e.g., [25,26]). To assess potential dispersal barriers associated with the Equatorial and Easter Microplate boundaries, we examined DNA sequences from mtCOI and 11 nuclear genes in geographical samples that extend the known range for A. pompejana northward to 23°N latitude in the Alarcón Basin and southward to 38°S on the Pacific Antarctic Ridge. The present study reports the most comprehensive geographical sampling and genetic analysis of this species, to date.

Samples
Samples were obtained during oceanic expeditions that spanned 21 years (Table 1, Fig. 1) with robotic manipulators or slurp guns on the human occupied vehicle (HOV) Alvin (Woods Hole Oceanographic Institution, WHOI) and the remotely operated vehicles (ROVs) Tiburon and Doc Ricketts (Monterey Bay Aquarium Research Institute, MBARI). Upon recovery at the surface, the samples were briefly stored in cold (2°C) filtered seawater prior to dissection and tissue removal. Tissue samples were frozen at −70°C on board the vessels and subsequently stored at −80°C in the land based laboratories. DNA sequencing was conducted with subsamples of individuals from each of the sample localities. Genomic DNA was extracted from muscle tissue with the Qiagen Blood and Tissue kit, following manufacturer's protocols (Qiagen, Hilden, Germany).

Molecular methods
Primer pairs were previously described for one mitochondrial and three nuclear protein-coding genes: mtCOI, cytochrome-c-oxidase sub-unit I [18]; SAHH, S-adenosylhomocysteine hydrolase; GlobX, globin X; and PGM, Phosphoglucomutase [20]. The SAHH marker includes part of an exon, whereas GlobX and PGM markers are composed entirely of introns.
We developed primer pairs for sequencing eight non-coding nuclear regions. High concentration DNA (100μg/μl) extracted from one individual from 14°S was sent to the National Instrumentation Center for Environmental Management (NICEM) at Seoul National University, Seoul Korea, for pyrosequencing on a Roche/ 454 Life Sciences Genome sequencer GS FLX machine. We obtained 28,222 reads with a total length of 9,804,791 base pairs. The pyrosequencing sequence reads were partitioned into protein-coding genes vs. non-coding DNA fragments identified through BLAST. Primer pairs were then designed for a subset of non-coding fragments containing > 450 base pairs (bp) and no poly (n) microsatellite repeats. Of these, we chose eight polymorphic loci, based on their PCR-efficacy, for further study: AP_NC1, AP_NC3, AP_NC8, AP_NC20, AP_NC22, AP_NC28, AP_NC32 and AP_NC43. The primer pairs used for nested PCRs are provided (Additional file 1: Table S1).
PCRs for population screening of 12 markers were performed under the following conditions: the first PCR initial denaturation temperature at 94°C/1 min followed by 35 cycles at 92°C/30 sec, annealing at different temperatures depending on the markers from 45 to 58°C/1 min, 72°C/1 min with a single final extension at 72°C/7 min, and second PCR repeated same conditions with first PCR, except for annealing temperature for 35 more cycles (Additional file 1: Table S1). All PCR products were sequenced with an ABI 3730xl automatic sequencer (Applied Biosystems, Foster City, CA).

Statistical methods
We used PHASE v. 2.1.1 [27,28] to resolve the phase of sequences with two or more heterozygous sites, and set the thresholds to 60 % with recombination model and stepwise mutation model. A number of sequences that could not be resolved in this manner needed to be cloned with three to five clones per individual. We used the pGEM-T easy vector (pGEM-T Easy Vector System, Promega, Madison, WI, USA) and DH5α (DH5α chemically competent E. coli, Enzynomics, Korea) for cloning and prepared the products for sequencing with the Hybrid-Q Plasmid rapidprep kit (GeneAll, Korea). All DNA sequences obtained in this study were deposited in GenBank (accession numbers: KX187433-KX189058, KX233878-KX233915 and KR868948-KR868986).
We used SITES [29] program to detect putatively nonrecombinant segments within the DNA sequences from all eleven genetic markers for Isolation with Migration Analysis (IMa2), Tajima's D, and Fu's F S test. Putatively recombinant fragments were excluded and the longest remaining fragment was used for subsequent analyses (as previously recommended, [30]). We used ARLEQUIN v. 3.5 [31] to estimate haplotype (H d ) and nuclear (π) diversity indices, Tajima's D [32] and Fu's F S [33]. Significance levels for Tajima's D and Fu's Fs were corrected by Bonferroni method. Allelic richness and private allelic richness was estimated by rarefaction methods implemented in HP-RARE v. 1.0 [34]. Tests for linkage disequilibrium (LD) between markers were conducted with GENEPOP v. 4.2 [35]. A haplotype network for each marker was constructed with HAPSTAR v. 0.7 [36] based on minimum spanning network resulting from ARLEQUIN v. 3.5 [31].

Analysis of population genetic structure
We used the PGDSPIDER [37] data format to prepare a table of multi-locus genotypes of nuclear genes for each individual. The individual genotypes were then examined with STRUCTURE v. 2.3.4 [38] to estimate the most probable number of discrete clusters (K) [38,39]. We let K range from 1 to 10 and repeated the simulations at least five times for each value of K using admixture model with correlated allele frequencies among populations [40]. All simulations included 5 × 10 6 MCMC (Markov chain Monte Carlo) generations after excluding the first 5 × 10 5 as 'burn-in'. The most probable value of K was  [41]. We used BAYESASS v. 3 [42] to assess recent immigration events. Each run included 1 × 10 7 iterations (−i) with random number seed (−s) and burn-in of 1 × 10 6 (−b), a sample interval of 100 (−n), allele frequencies of 0.3 (−a), and inbreeding coefficients of 0.4 for (−f).
We used GENODIVE v. 2.0b23 [44] to examine correlations (Mantel's r) between pairwise genetic differentiation (F ST ) based on nuclear genes and mtCOI gene separately, and geographic distances (from GEODATASOURCE <http://www.geodatasource.com/distance-calculator>). Because hierarchical subdivision can generate an apparent Isolation-by-Distance (IBD) pattern, we partitioned the samples into three regions, NEPR, SEPR, and PAR and applied a stratified Mantel test. The "Stratified" option in GENODIVE randomly permutes the data within partitions to test for the residual correlations between the genetic and geographic matrices.
Finally, we examined correlations between degrees north latitude, tectonic spreading rates (mm/year) for the sampled localities, and average haplotype diversities (H d ) for the twelve loci. Seafloor spreading rate was estimated under the NUVEL-1A model as implemented by the website < http://ofgs.aori.u-tokyo.ac.jp/~okino/ platecalc_new.html>. Correlation analyses were implemented with SPSS v.21 (IBM Inc.).

Isolation with migration analyses
We used the Isolation with Migration (IMa2) method [45,46] to infer six demographic quantities of model parameters of population divergence from DNA sequence data. First, we applied the IMa2 model to analyze neighboring pairs of population clusters (NEPR vs. SEPR and SEPR vs. PAR) previously identified with the STRUCTURE analysis. Because IMa2 assumes no recombination, we used non-recombinant segments of each of the 11 nuclear loci resulting from the SITES [29] program with mtCOI sequence data. We applied the Hasegawa-Kishino-Yano (HKY) [47] mutation model for mtCOI (inheritance scalar, H = 0.25) and the Infinite Sites (IS) model [48] for the nuclear markers (I = 1.0). Each analysis included at least 1.0 × 10 7 MCMC steps, and the first 1.0 × 10 5 steps were discarded as burn-in, with 40 attempts of chain swapping per step, 40 chains with geometric heating, and h1 and h2 values of 0.975 and 0.75, respectively. The mutation scaled model parameters of IMa2 were transformed into corresponding demographic quantities (t as years, m as migration rate per generation, and N as effective population size) using a mutation rate of mtCOI. Due to the absence of confirmed mutation rates in Alvinella pompejana, we borrowed a substitution rate of 1.0-2.0 % per million years for mtCOI, as estimated for marine taxa partitioned across the Isthmus of Panama [49,50]. Significance of migration rates between pairs of groups were analyzed with simple loglikelihood ratio tests [51].
The pairwise results were then used to analyze a three-population model (NEPR, SEPR and PAR) as described in [30,46]. We used splitting times obtained from the previous analyses as prior information for the three-population model. Running conditions assumed the same mutation models for mtCOI and nuclear markers as in the 2-population analyses. After 5.0 × 10 5 burn-in steps, at least 4.5 × 10 6 MCMC steps proceeded with 100 chains with geometric heating, and h1 and h2 values of 0.99 and 0.75, respectively, for 200 attempts of chain swapping per step. We conducted this process multiple times with the same options using random seeds, and genealogies were saved every 100 Markov chain steps (default). Finally, we combined 2.0 × 10 5 genealogies which were produced by multiple runs in L-mode, and tested significance with log-likelihood ratio test [52]. Demographic quantities of model parameters were visualized with the IMFig program [46].

Genetic diversity
The present suite of twelve genetic markers was polymorphic in nearly all of the vent samples (Table 2). Pairwise test for linkage disequilibrium among the eleven nuclear markers were non-significant for all the samples (P-value: 0.06-1.00). Noticeable regional differences existed in frequencies of the nuclear alleles and mitochondrial haplotypes. For example, the mtCOI haplotype network (Fig. 2a) exhibited haplotype clusters (shades of blue, green and red) that segregated between the NEPR (Fig. 2b, blue shades) and southern (SEPR + PAR) regions (Fig. 2b, red and green shades). Five mutational steps separated the NEPR and southern clusters. The predominant mitochondrial haplotypes in each region were complemented with numerous singletons (Fig. 2b, gray shades) and rare variants (i.e. q i ≤ 0.05), but the NEPR variants tended to differ from predominant haplotypes by more mutations, leading to substantially greater nucleotide diversity in the NEPR samples (π for mtCOI in Table 2).
Haplotype networks for eleven nuclear loci did not exhibit comparable fixed differences between the NEPR and the southern (SEPR + PAR) groups (Fig. 2a). Instead, allelic frequencies varied among populations in a geographically structured manner (Fig. 2b), for which the blue-shaded (NEPR) alleles yielded to green-shaded (SEPR) alleles, which yielded, in turn, to red-shaded (PAR) alleles. Although the frequencies varied among populations, they did not show congruent gradient among different loci. SAHH and AP_NC43 had the simplest haplotype networks, with a single dominate allele (q i ≥ 0.7) and variants that differed by only one or two mutations. AP_NC1 also had a single dominant allele, but some of NEPR variants (blue) differed by as many as 11 mutational steps. The remaining eight loci exhibited more diverse networks. In general, the NEPR samples had more singletons and rare alleles than SEPR and PAR samples. These regional differences were reflected in estimates of genetic diversity (Table 2). Aside from the regional effects, differences in diversity also existed among loci. For example, AP_NC32 exhibited the greatest overall nucleotide diversity (π = 0.0130, Table 2), whereas AP_NC43 had the smallest (π = 0.0014). Haplotype diversity (H d ) ranged from 0.44 (SAHH) to 0.87 (mtCOI) for the proteincoding markers, and from 0.41 (AC_NC1) to 0.87 (AC_NC22) for the non-coding markers. Estimators of genetic diversity declined with southern latitudes ( Table 2). For example, haplotype diversity (H d , r = −0.893, P < 0.001), rarefaction estimates of allelic richness (r = −0.912, P < 0.001), and private alleles (r = −0.867, P = 0.001) all decreased with southern latitudes. Estimates of allelic richness and of private alleles are correlated (r = 0.909, P < 0.001). Although H d is expected to be more sensitive to the evenness of allelic frequencies, in this case H d is almost perfectly correlated with richness (r = 0.997, P < 0.001). Tectonic spreading rates for the sampled localities ( Fig. 1) also increased with southern latitudes (r = 0.721, P = 0.019); consequently, spreading rates also were inversely correlated with genetic diversity H d (r = −0.724, P = 0.018) and richness (r = −0.740, P = 0.014).

Geographical structure
The STRUCTURE analysis identified three geographical clusters (K = 3, ΔK = 1431.57): NEPR (23N, 21N, 13N, and 9N); SEPR (11S, 14S, 17S, and 18S); and PAR (32S and 38S), north to south in order (Fig. 3a, Additional file 2: Table S2). Samples from the recently discovered vents at 23N and 38S clustered with those from the neighboring NEPR and PAR vents, respectively. The STRUCTURE analysis identified very limited evidence for mixed ancestries within the three clusters, without any noticeable gradient around the geographical boundaries. BAYESASS is capable of identifying recent immigrants from the other clusters ( Fig. 3b) Figure S1: for nuclear genes, r = 0.463, P = 0.002; and for mtCOI, r = 0.705, P = 0.003). However, assignment of each sample to its respective geographical region (i.e. strata), and application of the stratified Mantel test revealed no significant residual correlation (See figure on previous page.) Fig. 2 Haplotype networks of the twelve genetic loci. a The circles represent each unique haplotype which is scaled proportionally by relative frequency. The smallest white circles represent missed haplotypes in DNA samples and the line between the circles represents a single nucleotide difference between them. The three colors of the circles correspond to three geographic regions of the eastern Pacific ridges as illustrated in Fig. 1: blue = NEPR; green = SEPR; red = PAR. b Pie diagrams represent frequencies of twelve different loci in each sampling locality. For the sake of simplicity, some colors are omitted in the above haplotype network, but used here to represent different haplotypes. Gray color represents private alleles seen only at that single locality, and different private alleles are separated by a white line Fig. 3 Multi-locus genotypic assignments for nuclear markers. a STRUCTURE analysis with K = 3. Vertical bars indicate probability that an individual could be assigned to the NEPR (blue), SEPR (green) and PAR (red) clusters. b BAYESASS assignments of migrant ancestries. Color codes indicate posterior probabilities that individuals had native versus immigrant ancestries with geographical distance (P = 0.120 for nuclear genes; and P = 0.626 for mtCOI gene).

Isolation with Migration (IMa2) analyses
Analysis of a three-population IMa2 model (Fig. 4, Additional file 4: Table S3) detected evidence that southward migration (from NEPR into SEPR, 2N s m s = 0.61) exceeded northward migration (2N n m n = 0.15). A series of nested model tests on this observation resulted in statistically supports for the presence of gene flow (Additional file 4: Table S4 and S5). However, these tests could not reject other alternative hypotheses: an equal migration rate for  the southward and northward migrations, and unidirectional migration. Analyses also detected greater southward migration from the SEPR/NEPR ancestral population into the PAR population (2 Nm = 1.13 vs. 0.12), although unidirectional migration was not rejected (Fig. 4, Additional file 4: Table S4 and S5). However, gene flow from PAR into SEPR was greater after the split between NEPR and SEPR about 0.79 Mya (95 % HPD: 0.07-6.67 Mya) (2 Nm = 0.30 vs. 0.83). This northward direction was statistically supported by the log-likelihood ratio test. The estimated time of population splitting between the SEPR-PAR pair was much older,~4.20 Mya (95 % HPD: 2.42-33.42 Mya) (Fig. 4). Estimates of effective population sizes (N) can be ranked in the following order: NEPR > SEPR > PAR.

Demographic stability
Tajima's D and Fu's F S statistics can be sensitive indicators of demographic processes [32,33]. Altogether, 120 tests of these metrics for each locus in population resulted in only two significantly negative values, following Bonferroni corrections for multiple tests (Additional file 5: Table S6). However, samples sizes per gene per locality were limited and statistical power was low. Thus, we pooled individual samples within the NEPR, SEPR and PAR regions. Pooling samples might bias D and F S metrics in a positive direction if samples are genetically heterogeneous within each region, but the vast majority of estimates were still negative, and only a few were statistically significant, following Bonferroni corrections (Additional file 5: Table S7). None of the positive estimates were significant. Consequently, no substantive evidence for recent demographic bottlenecks or expansions was found.

Discussion
Deep-sea expeditions conducted during 2005 and 2015 obtained samples that extend the known range of Alvinella pompejana northward by 288 km to the Alarcón Basin in the Gulf of California (23°N latitude) and southward by 666 km to a northeastern segment of the Pacific Antarctic Ridge (38°S latitude). Multi-locus genotypes of individuals sampled from ten localities distributed across this range clustered into three metapopulation segments (Fig. 3). The Equator separates the NEPR and SEPR segments, and the Easter Microplate region separates the SEPR and PAR segments (Fig. 1). Large portions of the DNA sequence diversity in A. pompejana resided in the variation among individuals within each sample location: 46.6 % for mtCOI and 74.7 % for nuclear genes. Most of the remaining diversity resided in the differences among samples from the three regions: 52.8 % for mtCOI and 22.7 % for nuclear genes. Very small portions of the diversity resided in the differences among sample localities within regions: only 0.6 % for mtCOI and 2.6 % for nuclear genes. This lack of within-region differentiation reflected very high rates of gene exchange along each ridge axis (discussed below). Dispersal along contiguous segments of the three ridge axes appeared to be relatively unimpeded, despite high bottom-currents and distances of several hundred kilometers between active vents [20,53]. Developmental arrest of its embryos and delayed metamorphosis undoubtedly plays a significant role in A. pompejana's capacity for long-distance dispersal [54]. Disconnection and reconnection of active vent habitats over time due to frequent shifts of the magma supply might also contribute to the genetic homogeneity of A. pompejana populations distributed along a ridge axis [55].
Estimates of genetic differentiation (F ST values) increased significantly with geographical distances among the sample locations (Table 3; Additional file 3: Figure S1). This apparent Isolation-by-Distance (IBD) pattern was previously inferred to result from stepping-stone dispersal [4]. We observed the same pattern for nuclear and mitochondrial genes, but application of stratified Mantel tests, as developed by Meirmans [56], revealed that the pattern resulted from hierarchical subdivision. IBD-like patterns can result from a variety of processes including hierarchical structure, geographical selection gradients, secondary intergradation and range expansions (e.g., [57][58][59]).

Spreading rates, disturbance and genetic diversity
Tectonic spreading rates have been interpreted as surrogates for frequencies of habitat turnover due to local extirpations from tectonic and volcanic events and colonizations of new or "reborn" habitats [1,60]. Increased frequencies of habitat turnover were expected to reduce genetic diversity within localities and increase the homogeneity among the localities [61,62]. As predicted, haplotype diversity in the A. pompejana samples decreased significantly as tectonic spreading rates increase in southern latitudes (r = −0.724, P = 0.018). The PAR and SEPR axes exhibit "superfast" spreading rates of 141-151 mm/yr [63]. Reduced genetic diversity of these southern populations of A. pompejana and the co-distributed siboglinid tubeworm, Riftia pachyptila, correspond with the rapid cycles of habitat extinction and rebirth in this region [22]. Such demographic instability is also expected to leave other genetic footprints. Previous population genetic studies of A. pompejana reported evidence for recent demographic expansion within individual populations [18,20]. In contrast, the present estimates of Tajima's D and Fu's F S from a larger sample of genes did not corroborate these conclusions (Additional file 5: Tables S6 and S7). The present gene networks (Fig. 2a) do not exhibit the star-like clusters of haplotypes typically associated with such events [18].
The IMa2 analyses provide additional information about demographic processes. The NEPR cluster appears to have increased in size from the hypothetical ancestral population, whereas the southern SEPR and PAR clusters may have become smaller (Fig. 4). The trend of decreasing population size with southern latitudes (NEPR > SEPR > PAR) corresponded with increasing tectonic spreading rates (Fig. 1). Nonetheless, relatively large effective population sizes of each regional group (Additional file 4: Table S3; Additional file 5: Table S8) suggest that A. pompejana has maintained high site occupancy within each region. Indeed, A. pompejana is one of the genetically most diverse vent invertebrates studied to date [1]. This pioneer species is among the first animals to colonize nascent hydrothermal vents [64]. It can persist in hydrothermal flows approaching 50°C [11] that would appear to exclude potential competitors like the siboglinid tubeworms Riftia pachyptila and Tevnia jerichonana, and the bivalves Bathymodiolus thermophilus and Calyptogena magnifica. Coupled with its exceptional colonization abilities, much lower rates of local extirpation probably explain the ability of this species to retain such higher levels of genetic diversity.

The equatorial boundary
A deep strong eastward current crossing the East Pacific Rise at the Equator generates northern and southern gyres [65] and might impede along-axis dispersal of vent species that produce pelagic larvae. However, surface currents are unlikely to affect A. pompejana, a species that produces negatively buoyant larvae with benthic dispersal [18,64]. Gaps in the spatial or temporal frequency of hydrothermal habitats are expected to disrupt the dispersal of vent species [4]. Active vents supporting A. pompejana have not been reported within 900 km of the Galápagos Triple Junction region (Fig. 1), creating a large contemporary gap in the distribution of this species. Nonetheless, we cannot exclude the possibility that the Triple Junction region hosts active vents, because this region is not well explored. Furthermore, the age of this putative gap is unknown. It might be coincide with formation of the Hess Deep formed about 1 Mya [66]. Auspiciously, our estimated time of separation between NEPR and SEPR population segments of A. pompejana was about 0.79 million years ago (95 % HPD: 0.07-6.67 Mya; Fig. 4, Additional file 4: Table S3). Using similar methods, Plouviez et al. [20] estimated a slightly older time of separation, 1.2-1.3 Mya for A. pompejana. This slight difference must be due to different data sets and substitution rates used in both studies.
The Equatorial region creates a semipermeable barrier to dispersal by A. pompejana and several other vent species. The IMa2 analysis (Fig. 4) provided evidence for weak but statistically significant gene flow across the filter (Fig. 4, Additional file 4: Table S4 and S5). This pattern was also evident in the analysis of BAYESASS which exhibited a few candidates of recent immigration (Fig. 3). Nevertheless, the extent of gene flow does not seem to be sufficient to prevent A. pompejana from diverging across the Equatorial filter. The siboglinid tubeworms Tevnia jerichonana and Riftia pachyptila also exhibit evidence for partial isolation across this boundary, but their subdivision might be due to historical range expansions or recolonizations of the SEPR axis from NEPR sources [21,22]. Several gastropod limpets also clustered into groups separated by this boundary [19,67]. In contrast, the polychaete annelids Branchipolynoe symmytilida and Hessiolyra bergi, and bivalve mussel Bathymodiolus thermophilus exhibit no substantive genetic differentiation across this boundary [18,19,23].

The Easter Microplate boundary
Formation of the Easter Microplate resulted in severe topographic changes that originated 2.5-5.3 Mya [68,69]. Our estimates of the time of splitting between the PAR and SEPR population segments were not well resolved, but the process of divergence may be at least 3.3 million year (My) old (Fig. 4, Additional file 5: Figure S2 and Table S8). Moreover, Won et al. [2] hypothesized that strong currents create a contemporary barrier to dispersal of deep-sea species across this boundary. Geostrophic models and empirical evidence indicate such a presence of strong cross-axis currents in the 22-25°S region of the East Pacific Rise [70,71]. Recently, McGillicuddy et al. study [72] provides more information on the characteristics of larval dispersal under more realistic hydrographic influences along the ridge axis of EPR. Using simulation of particle transporting along a virtual mid-ocean ridge of the EPR, they found that the dispersal distance of numerical larvae decreases as the height where they stay above the ocean bottom increased. The result arises from the fact that as larvae closer to the bottom tend to be more influenced by strong currents along the flanks of ridge, they are transported farther than other ones staying higher from the bottom. The larval transporting simulation clearly showed that deep-sea currents on the ridge axis influence on the direction and distance of larval dispersal. Marsh et al. [7] also reported a long-distance (~100 km) dispersal potential of the larvae of tubeworm, R. pachyptila, when it met a favorable deep-sea flow along the EPR. Otherwise, most proportion of them might be retained near to the source population. It is manifest that the hydrodynamic effects on the larval dispersal are important and can contribute to shape population genetic structures of vent animals. However, the structures are also products of other extrinsic and intrinsic factors and their interactions, including ridge geomorphology, temporal stability of vents, larval developments and behaviors [1].
The formation of Easter Microplate and its associated deep-sea currents appear to be the primary cause for the parallel divergence observed in several vent species codistributing around the Easter Microplate boundary. Like A. pompejana, the siboglinid tubeworm Tevnia jerichonana also exhibits evidence for partial isolation across this boundary [21]. The vent mussels Bathymodiolus thermophilus and B. antarcticus meet and hybridize at this boundary [23]. Sister species of vent crabs Bythograea laubieri and B. vrijenhoeki also separate across this boundary [73,74].