A New Barrier to Dispersal Trapped Old Genetic Clines That Escaped the Easter Microplate Tension Zone of the Pacific Vent Mussels

Comparative phylogeography of deep-sea hydrothermal vent species has uncovered several genetic breaks between populations inhabiting northern and southern latitudes of the East Pacific Rise. However, the geographic width and position of genetic clines are variable among species. In this report, we further characterize the position and strength of barriers to gene flow between populations of the deep-sea vent mussel Bathymodiolus thermophilus. Eight allozyme loci and DNA sequences of four nuclear genes were added to previously published sequences of the cytochrome c oxidase subunit I gene. Our data confirm the presence of two barriers to gene flow, one located at the Easter Microplate (between 21°33′S and 31°S) recently described as a hybrid zone, and the second positioned between 7°25′S and 14°S with each affecting different loci. Coalescence analysis indicates a single vicariant event at the origin of divergence between clades for all nuclear loci, although the clines are now spatially discordant. We thus hypothesize that the Easter Microplate barrier has recently been relaxed after a long period of isolation and that some genetic clines have escaped the barrier and moved northward where they have subsequently been trapped by a reinforcing barrier to gene flow between 7°25′S and 14°S.


Introduction
Genetic structure is easier to detect and understand in a onedimensional system than in a two-dimension space [1]. However, even in a one-dimension space, detecting a genetic cline with a correlation between genetic differentiation and geographical distance does not always mean that populations are following an isolation-by-distance (IBD) model. Such a correlation can also be due to the presence of barriers to dispersal (e.g., [2]) or secondary contacts between previously isolated populations under expansion (e.g., [3]).
When detecting a genetic cline, one should consider the possibility that the location of this cline may be due to the presence of a natural barrier to dispersal because clines are expected to be trapped by such a barrier [4], [5]. Genetic clines depend on the relative impact of dispersal, selection (e.g., [6], [7]) and the recent demographic history of populations (e.g., [8]), and a natural barrier to dispersal impacts the balance among these parameters. As a result, clines often typify adaptive gradients [9], [10] or hybrid zones (i.e., regions containing recombinant individuals between genetically differentiated populations).
The geographic region in which a balance between dispersal and selection is maintained is defined as a tension zone. Tension zones tend to stabilize over natural barriers to dispersal [5], [7], [11], maintaining a cline around that barrier. They can also couple with a local adaptation cline and be stabilized at an environmental boundary [12]. If the origin of stabilized clines is difficult to establish from a single gene (e.g., isolation-by-distance without natural barriers to dispersal, hybrid zone), comparing allele frequencies between genes and gene divergences can help explain its emergence and, if a barrier is present, to identify its position more precisely.
Deep-sea hydrothermal vents are patchily distributed along mid-ocean ridges and back-arc basins. Along the East Pacific Rise (EPR), they follow a one-dimension pattern, ideal for testing an isolation-by-distance model of populations [13]. However, vent displacements along the ridge and eruptive phases leading to local faunal extinctions together with transform faults which likely impede gene flow are able to seriously alter expectations of such population models [14].
Comparative phylogeographic analyses of deep-sea hydrothermal vent species have previously shown the presence of a genetic break between the northern and southern regions of the EPR [15], [16], [17]. These studies established a shared vicariant event among species and suggested the emergence of a barrier to dispersal near the equator about 1.5 to 2 Mya [18]. However, the width and position of the barrier was not matching among species: some such as the tube-dwelling polychaete Alvinella pompejana displayed an abrupt separation of the populations across the equator [3] whereas others showed genetic patterns closer to an isolation-by-distance model [15]. Interestingly, the deep-sea mussel Bathymodiolus thermophilus exhibited a smooth clinal distribution of mitochondrial lineages along the EPR (13uN to 21u339S) [15], [17] and the presence of a cryptic species, B. aff. thermophilus further south (31uS-32uS) [15].
Bathymodiolus species have been supposed to be long-distance dispersers because of their planktotrophic mode of larval development [19] and the possibility of larvae reaching the upper (photic) layers of the water column [20]). The wide-dispersal capabilities of Bathymodiolus have been confirmed in several population genetic studies showing the absence of genetic differentiation across the Atlantic [21] and across the Gulf of Mexico (Mississippi Canyon and Alaminos Canyon, 550 km apart) [22]. Such dispersal characteristics could strongly favour population connectivity among geographically isolated sites. However, hybrid zones resulting in an abrupt change of allele frequencies over relatively short distances have also been observed in both the Atlantic [23] and the Pacific [24].
The likely cause of the clines observed at the EPR has been a source of debate. Although Plouviez et al. [17] suggested the observed genetic shifts might be the consequence of a natural barrier to dispersal that separated two interacting Bathymodiolus units, Audzijonyte & Vrijenhoek [25] proposed that the differentiation observed could simply be obtained under IBD. However, Johnson et al. [24] recently refuted IBD by describing a tension zone localised at the Easter Microplate between B. thermophilus and B. aff. thermophilus (renamed B. antarcticus by the authors). One of the loci, S-Adenosyl Homocysteine Hydrolase (SAHH), displayed a discordant cline position with fixed substitutions around the geographic zone previously identified by Plouviez et al. [17] to be a barrier to gene flow. The small likelihood of discovering fixed alleles over such a small spatial scale under IBD, and the fact that mitochondrial DNA exhibited an abrupt shift in allele frequency at the exact same position [17] prompted a reinvestigation of the hypothesis of a second barrier to gene flow.
In the present study, we obtained allozyme and DNA sequence datasets for four nuclear genes including the SAHH marker and analysed this new dataset together with mitochondrial haplotypes previously obtained from the EPR mussels to further investigate how genetic diversity is structured across the EPR and more specifically to test for the homogeneity of gene divergences across identified barriers. We obtained evidence for the existence of two barriers, including the one positioned at 7u259S-14uS. We thus propose that a pair of (semi-)permeable barriers along the EPR is likely responsible for the geographically discordant clinal distribution of alleles in this region.

Ethics statement
No specific permits were required to perform field studies described in this article. No specific permissions were required to access geographic localities and sample specimens (sampling sites belong to international waters). The locations are not privately owned or protected in any way. The field studies did not involve endangered or protected species.

Collection
Bathymodiolus thermophilus specimens were sampled from seven deep-sea hydrothermal vent fields along the East Pacific Rise (EPR) from 9u509N to 21u339S (Table 1) using the telemanipulated arm of the manned submersible Nautile operated from the oceanographic vessels Le Nadir and L'Atalante during three oceanographic cruises: at 9u509N during HOT 1996 and Mescal 2010 and from latitudes 7u259S to 21u339S during BIOSPEEDO 2004. During the three cruises, all fresh specimens were measured and dissected on board and tissues (mantle and muscle) were preserved in 80% alcohol. In addition, the anterior muscle of each individual was also frozen in liquid nitrogen for allozyme analyses during BIOSPEEDO 2004 and Mescal 2010.
The program Genetix 4.05.2 [27] was used to perform population genetic analyses on allozyme data. For each locus, allele frequencies, heterozygosities and Weir & Cockerham [28] f statistic (departure from Hardy-Weinberg equilibrium tested by a 1000-permutations test) were estimated for each population along the EPR. The overall genetic differentiation across populations was estimated at all loci using Weir & Cockerham's h estimator [28] and 1000 permutations used to determine significance. The isolation-by-distance model was tested with a Mantel Spearman test with 5000 permutations using Genepop 4.0.10 [29].

DNA sequencing
Genomic DNA was extracted using a CTAB-PVP extraction procedure following Jolly et al. [30]. Mitochondrial lineages of B. thermophilus were identified from cytochrome oxidase I gene (mtCOI) sequences previously obtained by Plouviez et al. ([17], Table 1). Sequences from three nuclear genes (GenBank accession numbers KC858658-KC858846, see Table S1) were obtained from the same individuals using the mark-recapture (MR) cloning technique developed by Bierne et al. [31] with two times the capture effort. MR-cloning was chosen over direct sequencing of PCR products because of the presence of insertions/deletions in intronic regions and to have access to linkage disequilibrium among polymorphic sites without the use of computer algorithms to determine allelic phase. Primers developed by Faure et al. [23] were used for amplification of introns in the S-Adenosyl Homocysteine Hydrolase (SAHH) and Lysozyme (Lyso) genes, previously obtained from EST sequences from a Bathymodiolus azoricus cDNA library. The SAHH gene contains a poly-A tract varying in length among individuals, the length polymorphism in this region was not included in the analyses due to the high potential error rate resulting from PCR and cloning. Specific primers for the Sulfotransferase (Sulfo) gene were also designed from the cDNA sequence data (BtSulfo-F: 59-TCTTTAAAGT-CAGGATCACATTGG-39, BtSulfo-R: 59-TAAGGCAAAGTG-GAACAACGAGACCGC-39). Sequences from Sulfotransferase were sorted in two paralogous genes called Sulfo1 and Sulfo2 from individual allele recaptures, respectively. Because of the low number of sequences recaptured from Sulfo2, only Sulfo1 sequences were used in this study.
Faure et al. [23] used a similar MR-cloning approach on two individuals from 37u709S for the three nuclear genes studied in the present paper. Faure et al. [32] also used a MR-cloning approach for a fourth nuclear gene (Elongation Factor 1a, EF1a), for which four of our populations were sampled and already sequenced prior to this population analysis. This gene was thus included in our analyses. In Johnson et al. [24], a study done in parallel to ours, the SAHH and EF1a genes were also used but alleles were obtained by direct sequencing. Because of the presence of insertion/deletion in the intronic region, Johnson et al. [24] were able to analyze only a quarter of the sequence length we obtained by MR-cloning. Consequently, this did not allow us to include these new datasets into our SAHH and EF1a analyses.
DNA sequences obtained from the MR-cloning method were visualized and edited using CodonCode Aligner 2.0.6 (http:// www.codoncode.com/aligner/). Sequence alignments were initially performed with ClustalW [33] and improved manually. The number of individuals for which at least one of the two alleles was recaptured (number of recaptured individuals) varied from one population to another is indicated in Table 1. Because of the random nature of the recapture, it was not possible to distinguish 'true' homozygotes from heterozygotes with the recapture effort. Table 1. Location and sample size of populations sampled for allozymes (n allo ), nuclear genes (n SAHH , n Lyso , n Sulfo1 , n EF1a ) and mitochondrial (n mtCOI ) gene.

Locality
Site name Latitude Depth n allo n Sulfo1 n Lyso n SAHH n EF1a n mtCOI * 110u879 W GR, Galapagos Rift; EPR, East Pacific Rise; PAR, Pacific Antarctic Rise. Depth is given meters. Sulfo 1, Sulfotransferase paralogue 1; Lyso, Lysozyme; SAHH, S-Adenosyl Homocysteine Hydrolase; EF1a, Elongation Factor 1a; mtCOI, cytochrome oxidase I. n SAHH , n Lyso , n Sulfo1 , n EF1a , total number of recaptured individuals for each nuclear gene. *, sequences from [15], [17].`, populations used for the Monmonier analysis. doi:10.1371/journal.pone.0081555.t001 Consequently, and because our main results is based on the coalescence theory, only the most recaptured allele was retained from each individual to avoid sample bias when performing demographic analyses and genetic diversity estimations (see :  Table 1). Multiple recaptures allowed us to discard intra-individual in vitro recombinants and putative artefactual/somatic mutations. Recombinants between different individuals (1-2% of the dataset for each population) from the same PCR set were detected (and removed) based on abnormal combinations of the 59-tails.

Polymorphism and divergence from DNA sequences
For each locus and locality, nucleotide diversity (p n ) and Watterson's theta (h w ) were estimated using DnaSP 4.10.3 [34]. Phylogenetic relationships among alleles were estimated using the median joining algorithm of the Network software (version 4.5.0.0; www.fluxus-engineering.com) [35] to detect potentially divergent clades. The geographic distribution of divergent clades was then examined to locate potential barriers to gene flow by plotting synthetic clade-specific allele frequency distributions for each locality. To test for gene divergence homogeneity across barriers, divergence time between clades was estimated using the formula T = D/2r) under the assumption of a local molecular clock (tested using the BEAUti/BEAST 1.4.8 package [36] with parameters previously described in [17]), where D is the average net divergence between geographic clades and r the mutation rate per site per million years [37]. The initial separation between B. thermophilus and B. azoricus (across the Isthmus of Panama, set to about 8-12 Mya for deep-sea fauna; [38]) was used as calibration point using a published dataset [23]. This calibration of deep-sea fauna separation across the Isthmus of Panama has been used successfully in other species [39].
Neutrality of loci was also tested using a multi-loci Hudson-Kreitman-Aguadé test (HKA, [40]) performed separately for each divergent clade, with one sequence of B. azoricus as an outgroup, via the software HKA ( J. Hey's web page: http://lifesci.rutgers. edu/,heylab/HeylabSoftware.htm#HKA). This test compares polymorphism and divergence at several loci to detect if at least one of these loci displays a departure from neutral evolution. This test was preferred to the McDonald-Kreitman [41] test because of the absence of fixed non-synonymous mutations between the two deep-sea mussel species.

Statistical evidence for a barrier to gene flow at 7u259S-14uS
A Monmonier algorithm was implemented using Barrier 2.2 [42] that compares matrices of multigene genetic distances (Fstatistics, w st , estimated with DnaSP 4.10.3 [34]) and geographic distances under the assumption of gene flow breaks. This Bayesian program was used to determine the geographic position of a potential barrier along the East Pacific Rise (insertions/deletions were coded as presence/absence). As the software Barrier 2.2 does not hold missing w st values in the matrix, the EF1a gene and some of the populations for which, at least one gene was missing were discarded from the analysis (see Table 1). The 38uS population was not included because of its small sample size. Localities from each side of the 7u259S-14uS barrier were grouped together to test for genetic differentiation across this barrier using w st [43] computed using DnaSP 4.10.3 (1000-permutations [34]).

Migration rates and demographic history of populations
Migration rates across the 7u259S-14uS barrier and the effective size of the southern and northern populations were estimated by fitting an isolation with migration model (IMa2 program, [44]) using previously described parameters [3], with the following exceptions: upper bounds of uniform priors set at h = 50 (population size), m = 5 (migration rate) and t = 30 (divergence time). Demographic and migration parameters were calibrated using divergence across Panama to inform mutation rates of loci (geometric mean among loci following a strict molecular clock) and a generation time of 2 years as previously estimated by Faure et al. [23].

Hybrid/introgressed individuals detected using SAHH RFLP analysis
A Hinf I restriction site polymorphism fixed between the two main sets of SAHH alleles was identified and used to check for the occurrence of putative introgressed/hybrid individuals between mitochondrial lineages by looking for individuals that possess one SAHH allele from each of the two distinct clades (called N/S individuals). An RFLP analysis was then performed using the Hinf I site to detect N/S individuals by screening all individuals from the Bathymodiolus collection. Incubation of PCR-products was done at 37uC for 1.5 hour in a 20 ml total volume containing 17 ml of PCR product, 1X buffer (supplied by the manufacturer) and 10 U of Hinf I (Ozyme TM ).

Geographic distribution of allozymes
A series of differentiation tests were performed on B. thermophilus populations located along the EPR. Among the eight allozyme loci, four (Mpi, Odh, Mdh 2 and Hk) were nearly monomorphic with the most frequent allele occurring at a frequency greater than 95% in all populations. The remaining four loci (Pgm, Lap, Gpi and Mdh 1) exhibited enough polymorphism to investigate their allelic distribution over the range of the EPR. Genetic differentiation was very low and not statistically significant overall (h = 0.016, P value . 0.05), showing the absence of any strong genetic structure at allozyme loci along this portion of the EPR. The Mantel Spearman test showed a slight but significant correlation (P value = 0.05) between the genetic distance (h/(1-h)) and the geographic distance between vent fields (Fig. 1).
Nuclear gene networks exhibited two clades separated by a pronounced net divergence (Fig. 2). The SAHH gene revealed the presence of two clades (2% divergence) well established from each part of the previously suggested 7u259S-14uS break (Fig. 2, Fig. 3, [17]). SAHH clade 1 contains the sequences from most individuals sampled at 9u509N or 7u259S whereas clade 2 corresponds to individuals sampled only from the southern EPR sites. Two 0.3%divergent clades were also found using the Lyso gene but are more difficult to attribute to a particular geographic area (Fig. 2, 3). However, a 1-bp deletion was only found in the intron of sequences sampled in the southern EPR sites below the latitude of 7u259S (Fig. 3), indicating that the 7u259S-14uS barrier is playing a role in impeding the spread of new mutations. The Sulfo1 gene also displayed a clear geographic structure between north and south of the barrier with a 0.5% divergence between the two clades (Fig. 2).

Statistical evidence for a 7u259S-14uS barrier
The Monmonier analysis identified a barrier between 7u259S and 14uS for all tested loci (Fig. 3) and w st values were significantly different from zero for all sampled genes across this barrier ( Table 2). An Isolation with Migration (IMa2) analysis was performed between the southern and northern populations across this barrier. The marginal posterior probability distribution of migration rate across the barrier overlapped with zero, indicating that the absence of migration cannot be ruled out. If present, migration across this barrier could have occurred in both directions, possibly slightly orientated from north to south ( Table 3). The estimated effective population sizes of the present populations were both greater than the ancestral (N A ) population size ( Table 3), indicating that present-day populations of the vent mussel may be expanding, possibly at a higher rate in the south (N N being slightly higher than N S ). However, expansion cannot be confirmed because of the overlap range of Highest Posterior Density among north, south and ancestral populations.

Gene divergence homogeneity across the barrier
Under the verified assumption of a molecular clock (BEAST analysis) and using a 8-12 Mya calibration time obtained from the splitting of B. thermophilus and B. azoricus across the Isthmus of Panama, divergence times between clades 1 and 2 were estimated at 0.6-0.9 Mya for mtCOI, conforming to geological estimates of the ages of the transform faults in the 7u259S-14uS area (Fig. 4). In contrast, divergence times between mtCOI clades 2 and 3 (3.0-4.5 Mya) as well as between Sulfo1 (4.4-6.6 Mya) and between SAHH (4.0-6.1 Mya) clades conformed to geological estimates for Easter Microplate formation (Fig. 4). Divergence times between clades for the Lyso (0.4-0.6 Mya) and EF1a (3.1-4.7 Mya) genes have to be interpreted with caution because of mutation rate heterogeneity among clades (departure from a strict molecular clock).
The multi-locus HKA test showed a significant departure from neutral evolution (P , 0.02) among genes within clade 1 (no departure among genes within clade 2), indicating that at least one gene could be under selection. Sulfo1 displayed the highest divergence to polymorphism deviation when compared to other genes. The removal of Sulfo1 (possible outlier) from the dataset led the HKA test to refit the model of neutral evolution (P . 0.31). This gene displayed a polymorphism/divergence ratio of 8.66: a value that is more than eight times greater than expected (1.05) with a very low divergence between B. thermophilus and B. azoricus, indicative of balancing or strong purifying selection.

Detection of introgressed individuals using the SAHH gene
The RFLP analysis of the SAHH marker revealed that only 13 individuals south of the 7u259S-14uS barrier displayed at least one nuclear allele corresponding to the northern clade. Except for one individual showing two northern-type alleles in the southern populations, all of these individuals possessed one allele of each clade. These N/S individuals were mainly located at 14uS and decreased abruptly in frequency further south (Fig. 5). No N/S individual was detected at 7u259S or further north. N/S individuals of SAHH had Sulfo1 alleles from clade 2 only but they were not associated with a specific mtCOI clade: the SAHH northern-type allele being found in association with haplotypes from both the mtCOI clades 1 and 2, indicating that these individuals are most probably introgressed rather than first generation hybrids which should result in cyto-nuclear disequilibrium [45], [46].

Discussion
The unpredictable nature of fluid circulation and the high level of fragmentation found at deep-sea hydrothermal vents should result in species with high dispersal capabilities in order to promote longdistance (re)colonization of new vent sites [47]. In accordance with this expectation, most species from the genus Bathymodiolus have high dispersal capabilities [20], [21], [22]. Data from B. thermophilus supports this hypothesis in showing a complete lack of differentiation at allozyme markers between the northern EPR (13uN, 11uN,  9uN) and the Galapagos sites [48] (but see [49]). Having high dispersal capabilities, however, does not necessarily imply a lack of geographic structure because gene flow is likely impacted by physical barriers to dispersal, great distances without suitable habitat and/or hybridization fronts [23].
Previous and present works strongly support the hypothesis that the population structure of Bathymodiolus thermophilus is impacted by two barriers to gene flow along the EPR: the so-called Easter Microplate barrier ( [15], [24], present manuscript) and the 7u259S-14uS barrier. This 7u259S-14uS barrier was first suggested by Plouviez et al. [17], but then challenged by Audzijonyte & Vrijenhoek [25] who proposed that frequency changes in the two most divergent mtCOI clades found in the mussel populations along the EPR were most likely due to isolation-by-distance (Mantel Spearman test performed on the Won et al.'s dataset found significant [15]) and not an actual barrier. Johnson et al.'s study [24], together with our present work, refutes the hypothesis that isolation-by-distance alone might be responsible for the genetic patterns observed along the EPR, instead these studies suggest that a tension zone across the Easter Microplate might explain the clinal distribution of alleles. Johnson et al. [24] proposed the occurrence of a hybrid zone at 23uS and confirmed that alleles from the Pacific-Antarctic mussel B. antarcticus were likely to introgress in the South EPR populations of B. thermophilus but did not discuss the possible presence of a barrier to gene flow further north (7u259-14uS).
In the present manuscript, we confirm the existence of such a barrier for the deep-sea mussel and argue that the observed positive correlation of w st with geographic distances at the mtCOI locus (significant Mantel Spearman tests) likely reflects the occurrence of a tension zone that originated at the Easter Microplate (as proposed by Johnson et al. [24]), but which subsequently moved northward and became captured by a second area of restricted larval exchanges (responsible for the cline), the 7u259S-14uS barrier to gene flow. As for many genetic barriers it is semipermeable and affects loci differentially.

Genetic differentiation of B. thermophilus is triggered by the recent formation of a series of transform faults at 7u259S-14uS
The Monmonier analysis and coalescence analyses statistically detected the presence of a barrier to gene flow at 7u259S-14uS using nuclear and COI genes. Sulfo1 and SAHH loci both showed two sets of divergent alleles distributed in the 9u509N-7u259S and 14uS-21u339S regions, with an almost reciprocal monophyly. The gene encoding Lyso exhibited two sets of divergent alleles that were geographically interspersed over the whole EPR and PAR, indicating that they were able to cross both the Easter Microplate and the 7u259S-14uS barriers. However, the 1-bp deletion allele found at a high frequency in populations from 14uS-21u339S was not observed from 9u509S to 7u259S nor within individuals from 38uS. The sample size (2 individuals) is far too low to exclude the absence of this deletion at 38uS. The absence of the deletion in the 9u509N and 7u259S individuals, however, strongly suggests that the spread of this new deletion has been blocked since its first occurrence in the southern populations (at least across the 7u259S-14uS barrier). This deletion is however old enough to provide recombinant alleles at a high frequency between the two Lyso divergent clades, as the deletion is present at the tip of both divergent clades (see Fig. 2). The geographic distribution of this deletion on the Lyso gene suggests that this newly-formed barrier is now impermeable or only sporadically permeable to gene flow. The multigenic estimation of migration rates using a Bayesian approach and the mitochondrial and nuclear sequences (IMa2) also fit the relative isolation of mussels located at 7u259S or 9u509N when compared to the more southern populations. The IMa2 analysis indeed showed the absence or the very limited number of mussel migrants across the 7u259S-14uS barrier. Although sequences of B. thermophilus from 11uS were not available for the other genes, we hypothesize this barrier is located, more precisely, between 11uS and 14uS based on haplotype frequencies of the mtCOI gene (w st ).
The barrier was however not sufficient to create an allele frequency shift at allozyme markers, highlighting an apparent discrepancy among marker types (allozymes vs. DNA sequences). Discrepancies in the genetic differentiation observed with allozyme and DNA markers have regularly been suggested to result from selection at allozyme loci, sometimes claimed to be under balancing selection [50], [51] and sometimes under disruptive selection [52], [53]. However, subsequent analyses often refute the hypothesis of selection on allozymes by demonstrating insufficient sampling [54], [55]. Here, EF1a is sufficient to show that a low level of differentiation is observed at a DNA marker in the area sampled for allozyme analysis. Furthermore, Johnson et al. [24] did not observe any genetic differentiation north of the Easter Microplate at the three DNA markers that are not in common between the two surveys. We conclude that there is no real discrepancy between the two categories of markers and that the few allozyme loci analyzed simply fell in genomic regions unlinked to isolation genes. These loci are thus able to cross the barrier freely and to organize themselves according to geographic distance, producing an IBD pattern. Moreover, the barrier may not be detectable on some of the allozyme loci (Mpi, Odh, Mdh 2 and Hk) because of their low level of polymorphism.
Geologic and hydrodynamic conditions from the equator to 15uS are consistent with limited gene flow observed for multiple species between north and south EPR [3], [15], [16], [17]. In terms of geology, dispersal in this region can be impeded by the equatorial triple junction between the EPR and the Galapagos Rift, as well as the Gofar/Discovery multiple transform fault complex near 4uS [56], [57], [58]. When looking at local hydrodynamism in the region, surveys of He-3 plumes produced by the venting activity along the EPR also indicated the occurrence of a westward flow centered at 15uS [59], [60]. This flow, possibly linked to the anticyclonic circulation of water masses in the eastern Pacific [61], creates strong cross-axis currents able to produce a hydrodynamic barrier at these latitudes. The cooccurrence of these geologic and hydrodynamic features represents ideal conditions for the establishment of new barriers to gene flow.  A dynamic effect of a couple of physical barriers to gene flow Based on tectonic plate history, the Easter Microplate barrier is older than the 7u259S-14uS barrier. The Easter Microplate originated from the progressive offsetting of the Pacific-Antarctic Ridge (PAR) along with the 'old' EPR (Nazca Pacific spreading centre) about 3.88 Mya (anomaly 3). This offsetting progressively expanded with the formation of an overlapping spreading center (OSC: see [62]). Ridge offsetting is known to have a profound impact on gene flow [17], [63], [64] and to isolate populations from each other. In the specific case of the Easter Microplate, species divergence time was thus expected to coincide with (or to be slightly older than) anomaly 3 (i.e., 3.88 Mya) at the Easter Microplate barrier. In comparison, transform faults responsible for offsetting the ridge axis between 9uN and 17uS initiated 1-2 Mya [56], [57], [58]. If the two geographical barriers were both impermeable to dispersal since their occurrence, one would expect to find reciprocal monophylies for nearly all genes, with a greater divergence, at the older barrier (i.e., the Easter Microplate barrier) and incomplete lineage sorting at the most recent barrier (i.e., the 7u259S-14uS barrier).
Contrasting results among genes across the Easter Microplate are in accordance with the presence of a tension zone at the latitude of Easter Island as proposed by Johnson et al. [24] who studied the PAC/EPR population connectivity between 21uS to 38uS in more detail. The present results on mtCOI integrating both Won et al. [15] and Plouviez et al. [17] datasets are consistent with geological evidence, indicating that gene flow is either blocked or greatly impeded by both the Easter Microplate and the 7u259S-14uS barriers. MtCOI haplotypes indeed show a reciprocal monophyly with a 2.1-4.3 Mya divergence time across the Easter Microplate [15], while populations display only haplotype frequency differences across the 7u259S and 14uS barrier (significant Q st , [17] and Table 2). Absence of detectable differentiation on the EF1a gene across the more recent 7u259S-14uS barrier ( Table 2) while populations were divergent across the older Easter Microplate barrier was expected considering nuclear genes generally evolve slower than mitochondrial genes. In contrast, the other nuclear gene structure did not fit this expectation; showing no differentiation across the Easter Microplate, whereas two of them (SAHH, Sulfo1) displayed clear but unexpected patterns of reciprocal monophyly across the more recent barrier (7u259S-14uS). Such an inconsistency among genes might be explained by the presence of a tension zone extending further north toward 7u259S associated with the relaxation of the Easter Microplate barrier.
The hypothesis of a relaxation of the barrier at the Easter Microplate is also supported by the geological formation of the microplate itself. If initiated by the offsetting of the PAR and 'old' EPR, the plate formation was achieved when the ends of the two expanding ridges joined together following a clockwise rotation with the setting up of insular volcanic arcs (Orongo and Pito rifts) about 2.5-1.75 Mya (end of anomaly 2a: see [65]). This junction between the PAR and EPR is likely to have provided steppingstones for episodic bursts of larvae between the PAR and the southern EPR (using vents located on active seamounts: e.g., the Pito Seamount [66]).
We propose that the progressive ridge-offsetting at 7u259-14uS played an important role in enhancing barriers to dispersal by efficiently trapping endogenously-induced genetic clines associated with the reconnection of B. thermophilus and B. antarcticus, which could have subsequently escaped from the relaxing barrier at the Easter Microplate. Divergence times between clades for the SAHH and Sulfo1 nuclear genes across the 7u259S-14uS barrier (4-6.6 Mya) were similar to those estimated at the Easter Microplate boundary between the PAC and SEPR clades of mtCOI and EF1a. This convergence in divergence times found at two distinct spatial locations, as well as the fact that nuclear mutation rates should be lower than that of the mtCOI (as found in the Atlantic Bathymodiolus sp. [23]) fits this ''dual relaxing/ enhancing barriers'' hypothesis. Indeed, the fixation of alleles in the two observed-divergent clades associated with the SAHH and Sulfo1 loci (and possibly alleles from the Lyso gene as well) might have been initiated by the separation of populations during the  Easter Microplate formation before migrating northward and being trapped by the emerging 7u259S-14uS barrier. To conclude, we have shown that the clinal distribution of mtCOI haplotypes along the EPR was in fact due to highly dynamic historical/geological processes. A system of two genetic barriers to gene flow, in which one would be progressively relaxed (i.e., Easter Microplate) and the other enforced (i.e., 7u259S-14uS transform faults), can explain the discordant distribution of genetic clines. A volcanic arc joining PAR and EPR probably played an important role in setting up of a secondary contact zone (and subsequent northward spreading of alleles/genetic clines) between individuals across the Easter Microplate. The proposed relaxation/enforcement scenario of barriers is consistent with the occurrence of a biogeographic transition zone along the southern EPR (resulting from the overlap of the North-EPR and the South-Easter Microplate biogeographic provinces) to explain the high diversity of vent fauna species observed between 17u259S and 21u339S [67]. Studies on motion of hybrid/tension zones are in their infancy (but see e.g., [12], [68], [69], [70]). Tension zone movement is expected theoretically although they are expected to move toward areas of lower population density [5], [71]. Evidence for the movement of hybrid zones has been obtained directly following known tension zones over time (e.g., [68]) or indirectly using molecular and/or phenotypic traits (e.g., [70], review in [69]). Theoretical analyses have demonstrated the ability of lowmigration areas (i.e., lack of suitable habitats/larval transport disruption) to trap genetic incompatibilities previously generated by other historical/ecological/biogeographical processes (see: [12]) and might apply to the case of B. thermophilus, which, in turn represents one of the first empirical observation of such theoretical predictions. The theory of genetic barriers and hybrid zones explains semi-permeable barriers to gene flow and the structure observed, possibly indicating selection against hybrids. Of course one could also speculate that different loci respond to different ecological pressures at the two positions. However, not only are the possible ecological differences unidentified but getting by chance a series of outlier loci (or hichhicked introns) under differential ecological selection out of four genes is highly unlikely. This therefore suggests that allele frequency shifts are linked to tension zones and are thus able to move. Together with what we know of the geological history of the EPR, the hypothesis of an Figure 5. Localisation of individuals with both north and south type alleles in SAHH using RFLP analyses. N/S individuals, individuals that present one allele from the northern clade (blue clade in Figure 2) and one from the southern clade (yellow clade in Figure 2). Left Y-axis is the number of N/S individuals. Right Y-axis is the total number of individuals. doi:10.1371/journal.pone.0081555.g005 interchange of clines from the southern barrier to the northern barrier explains the geographic discordance of cline positions.