A hybrid zone between Bathymodiolus mussel lineages from eastern Pacific hydrothermal vents

Background The inhabitants of deep-sea hydrothermal vents occupy ephemeral island-like habitats distributed sporadically along tectonic spreading-centers, back-arc basins, and volcanically active seamounts. The majority of vent taxa undergo a pelagic larval phase, and thus varying degrees of geographical subdivision, ranging from no impedance of dispersal to complete isolation, often exist among taxa that span common geomorphological boundaries. Two lineages of Bathymodiolus mussels segregate on either side of the Easter Microplate, a boundary that separates the East Pacific Rise from spreading centers connected to the Pacific-Antarctic Ridge. Results A recent sample from the northwest flank of the Easter Microplate contained an admixture of northern and southern mitochondrial haplotypes and corresponding alleles at five nuclear gene loci. Genotypic frequencies in this sample did not fit random mating expectation. Significant heterozygote deficiencies at nuclear loci and gametic disequilibria between loci suggested that this transitional region might be a ‘Tension Zone’ maintained by immigration of parental types and possibly hybrid unfitness. An analysis of recombination history in the nuclear genes suggests a prolonged history of parapatric contact between the two mussel lineages. We hereby elevate the southern lineage to species status as Bathymodiolus antarcticus n. sp. and restrict the use of Bathymodiolus thermophilus to the northern lineage. Conclusions Because B. thermophilus s.s. exhibits no evidence for subdivision or isolation-by-distance across its 4000 km range along the EPR axis and Galápagos Rift, partial isolation of B. antarcticus n. sp. requires explanation. The time needed to produce the observed degree of mitochondrial differentiation is consistent with the age of the Easter Microplate (2.5 to 5.3 million years). The complex geomorphology of the Easter Microplate region forces strong cross-axis currents that might disrupt self-recruitment of mussels by removing planktotrophic larvae from the ridge axis. Furthermore, frequent local extinction events in this tectonically dynamic region might produce a demographic sink rather than a source for dispersing mussel larvae. Historical changes in tectonic rates and current patterns appear to permit intermittent contact and introgression between the two species.


Background
Hybrid zones, regions of intergradation between genetically distinct populations or species, have been documented for a number of terrestrial and aquatic organisms. The locations, shapes and ages of these zones, and the unique dispersal, behavioral and fitness characteristics of the participating organisms have motivated a number of models that invoke endogenous (i.e. genetic) or exogenous (i.e. environmental) forces reviewed by [1][2][3][4][5]. Intergradation can result from primary or secondary processes that are difficult to distinguish [6]. Primary intergradation, which occurs at the contact boundaries between distinct habitats or along steep environmental gradients, is associated with parapatric speciation. Allelic and phenotypic frequencies often exhibit step-clines shaped by disruptive selection across the contact zone. Secondary intergradation, resulting from contact between previously allopatric populations, can produce similar step-clines. Different histories of dispersal and selection can produce similar patterns of concordance in allelic frequency clines at multiple gene loci, though subtle and discriminating footprints may remain in the genetic architecture of intergrading lineages [1].
Most examples of hybrid zones come from terrestrial and freshwater environments that have experienced past climatic changesfor example, at suture zones where formerly isolated biotic assemblages attain secondary contact [7,8]. Hybrid zones are uncommon in the open ocean, because marine environments tend to be more homogeneous temporally and spatially, and many marine species have high dispersal capabilities [9]. Nonetheless, marine hybrid zones are found in coastal environments that experience steep environmental gradients or were impacted by past climatic changes [10]. Hybridization is prevalent among subtidal mussels, Mytilus edulis, M. galloprovincialis, and M. trossulus, members of the "blue mussel complex" [11,12]. These closely related species are known to readily hybridize wherever two blue mussels are sympatric [13][14][15]. A mid-ocean hybrid zone was first discovered in deep-sea hydrothermal vent mussels (Mollusca: Bivalvia: Mytilidae) living along the Mid-Atlantic Ridge [16]. Hydrothermal vents differ from other open oceanic habitats because they occur as discrete island-like habitats distributed along the global mid-ocean ridge system (Figure 1), in disjunct back-arc spreading centers, and on active seamounts [17]. Because vent communities rely almost entirely on reduced volcanic gases (primarily H 2 S and CH 4 ,) and chemosynthetic bacteria for primary productivity, habitat quality can vary greatly in time and space [18]. The larvae of most vent animals are dispersed by hydrographic circulation that can be disrupted by geomorphological structures, such as ridge offsets (e.g., transform faults), intervening seamounts, and other bathymetric features [19,20]. Consequently, genetic studies have often revealed geographic subdivision and species boundaries associated with such features [21,22].
The Mid-Atlantic Ridge (MAR) mussels, Bathymodiolus azoricus and B. puteoserpentis, mussels segregate latitudinally and hybridize at an intermediate location that is relatively inhospitable for mussels [16,23]. The hybrid population exhibits an excess of the parental cytonuclear genotypes, leading Won et al. [24] to hypothesize that the contact region is Tension Zone [sensu 1] maintained primarily by immigration from the parental zones and possibly hybrid unfitness. A coalescence analysis of the nuclear GF1B locus identified recombination events between branch tips belonging to the parental lineages [25]. The degree of mitochondrial divergence and the apparent absence of earlier inter-lineage recombination events at GF1B were interpreted as evidence for a period of allopatric differentiation beginning about 0.760 million years (Myr) ago, followed by prolonged secondary contact and introgression.
Discovery of the MAR hybrid zone stimulated our interests in other potential contact zones along the global mid-ocean ridge system. Our target was the Easter Microplate boundary (EMB, Figure 1), where vent fauna associated with the 7000-km long East Pacific Rise (EPR) contact a related fauna associated with a NE segment of the Pacific-Antarctic Ridge (PAR) [22]. The EMB separates geminate species of bythograeid crabs [26,27], lepetodrilid limpets [28], and Bathymodiolus mussels [29]. The EMB does not contribute to significant geographical partitioning of all lepetodrilid species [28] and annelid species [30,31], however. Instead, this region acts as a "variable dispersal filter" with affects that depend on the life history characteristics and dispersal modes of individual taxa [22]. To better define this contact zone, we explored and sampled vents that were previously reported to exist in the region between 23°S on the NW flank of the Easter Microplate and 38°S latitude along the PAR [32,33].
Herein, we report the discovery of a zone of intergradation between two mussel lineages that contact one another at the EMB. Bathymodiolus thermophilus Kenk & Wilson [34] is distributed on the Galápagos Rift (GAR) and East Pacific Rise (EPR) between 13°N and 21°S latitude. A related lineage, Bathymodiolus aff. thermophilus [29], was found at 31-32°S latitude on segments of the EPR axis connecting the Easter and the Juan de Fernandez microplates. The mussel lineages are reciprocally monophyletic for distinct mitochondrial COI haplotypes that are 4.6% divergent, and they exhibit significant shifts in allozyme frequencies at five polymorphic loci. These concordant differences at multiple, independent, gene loci led Won et al. [29] to suggest that the southern lineage might warrant consideration as a distinct species, but additional samples and analyses of independent nuclear gene loci were needed to assess this hypothesis. With these goals, we examined mitochondrial COI and portions of five nuclear genes from mussels sampled along a 5600-km range from 13°N to 38°S latitude ( Figure 1). Multilocus genotypes constructed from the DNA sequences were used to assess geographical structure and subdivision, to identify putative hybrids, and to test for gametic disequilibrium and recombination in the contact zone.

Samples
Expeditions conducted between 1990 and 2005 explored segments of the GAR, EPR and PAR ( Figure 1) to obtain specimens for multi-species phylogeographic studies reviewed in [22]. Exploration in the Easter Microplate region was guided by shipboard coordinates and towedvideo surveys reported in a PhD thesis [35]. Mussel samples were obtained from 14 localities, subsequently labeled with abbreviations that include the approximate latitude, e.g., N13, N11, . . . S31, S38 (Table 1). Samples collected with the human-occupied vehicle (HOV) Alvin were placed in an insulated 'biobox' containing ambient seawater at 2-4°C. Adult mussels (up to~15 cm in length) were sampled with scoops, nets or directly with Alvin's robotic manipulators. The S23 sample of juveniles (≤25 mm) was attached by byssal threads to a rock sampled for geological studies. Upon recovery of the submersible, all samples were stored temporarily in refrigerated seawater prior to dissection of gill and adductor muscle tissues that were frozen at −80°C or preserved in 95% ethanol.

DNA methods
DNA extraction and purification, general PCR conditions, amplicon purification and DNA sequencing used methods that were previously reported for mussels [37]. Six primer sets were used to amplify DNA targets from mussels ( Table 2). All PCR products were diluted in 40 μl of sterile water and purified with a Multiscreen HTS PCR 96 vacuum manifold system (Millipore Corp. Billerica, MA). PCR products were sequenced bi-directionally on an ABI3130 sequencer with BigDye Terminator v3.1 (Life Technologies Corp., Carlsbad, CA) chemistry and primers used in PCR.

Molecular statistics
Computer programs used to format the data, obtain estimates of population genetic parameters, and conduct statistical tests are listed along with procedural references (Table 3). We used PHASE v. 2.1.1 to ascertain alleles in individuals heterozygous for nuclear genes. Characterization of alleles in individuals heterozygous for insertions and deletions (indels) required cloning and sequencing of at least five clones per individual methods in [39]. Analyses of genic diversity, Hardy-Weinberg   equilibrium, gametic disequilibrium, Isolation-with-Migration (IMA2), and Isolation-by-Distance (IBD) were limited to the exon portions of four nuclear genes that contained introns. None of the exon regions in these sequences exhibited stop-codons. The states of a 9-bp indel in the EF1α intron were treated as two-alleles. Tests for intra-genic recombination were conducted with PHASE v. 2.1.1. Phylogenetic analyses of historical recombination within nuclear genes were conducted with BEAGLE on complete sequences, including introns and exons. Gene networks for the exon portions of each locus were constructed with NETWORK PUBLISHER v. 1.3.0.0. Inferred haplotypes were re-labeled as discrete alleles (a 1 , a 2 , . . ., a n ) and concatenated manually into a multilocus genotype (MLG) file (a 1 /a 3 , b 2 /b 2 , . . ., g 1 /g 1 ). Uninformative singleton alleles were lumped with their most closely related haplotype. The MLG file was used as an input file for GENEPOP v. 4.0.10. Exact tests for Hardy-Weinberg equilibrium (HWE) were conducted according to the method of Raymond and Rousset [49]. The MLG file was also used as an input file for ARLE-QUIN v. 3.5.1.3 to conduct pairwise likelihood ratio tests for gametic disequilibrium and Mantel tests for Isolation-by-Distance (IBD). Sequential Bonferroni corrections were used to adjust α-levels [58].
The MLG file was also used to conduct assignment tests. STRUCTURE v. 2.3.3 analyses were conducted with and without prior information about population samples and admixture, and with correlated and uncorrelated allelic frequencies. The first 10 3 generations (burn-in) were discarded and 10 7 generations were iterated for 3 to10 times per estimate of K, the number of genotypic clusters. Assuming a uniform prior on K, we estimated a Bayes factor for each value of K after averaging the ln Pr(X/K) values between simulations. NEWHYBRIDS was used to assign the multilocus genotypes to putative parental, F1, F2, or backcross categories. We used default genotypic classes with no prior information on allelic frequencies and included uniform and Jefferey's priors [59] for θ and π. Five separate analyses were conducted with 10 5 sweeps of burn-in and 10 7 sweeps of data collection. Threshold values (Tq) of 0.9 [56] and 0.75 [60] were assigned separately for the genotypic categories if q ≥ Tq, and they were left unassigned if q < Tq [61].
We used the Isolation-with-Migration method, IMA2 [50,62,63] to estimate the following demographic parameters: (τ) the time of population splitting; (θ A , θ N , θ S ) the sizes of ancestral and descendant populations; and (m N , m S ) immigration into the descendant populations. An HKY substitution model with back mutation was assumed. An inheritance scalar to adjust the relative expected effective population size was 0.25 for mitochondrial COI and 1.0 for nuclear loci. Analyses were conducted under a twopopulation model based on STRUCTURE results. The two-population model lumped the N13-S23 and GAR  [ 56,57] samples into a northern group, and the S31-S38 samples into a southern group. The S23 sample was treated as northern due to the preponderance of northern alleles. A second set of analyses was also run without the S23 population. Analyses involved at least 10 8 steps with the first 10 4 steps discarded as burn-in (−b) [64], with 50 attempts at chain swapping per step (−k), 50-80 chains (−n) with geometric heating (−f), and g1 and g2 values of 0.99 and 0.3, respectively. Analyses were conducted multiple times and the number of steps between saving genealogies (−d) was increased to 1000. Initially we set upper bounds for each parameter, as recommended in the IMA2 manual, then performed multiple runs using random number seeds with more appropriate bounds based on the initial output.
We assessed the convergence of MCMC samplings of population parameters in multiple ways by allowing analyses to run ensuring effective sample sizes were a minimum of 100 for all parameters, examining output graphs for unimodality, ensuring genealogies were updating, and checking that multiple runs gave similar estimates. Multiple runs based on the same priors were then combined in L-mode and log-likelihood ratio tests were performed to test for isolation.

Results
Adult mussels were abundant at most of the active vent fields between 13°N and 38°S latitude ( Figure 1; Table 1). The S23 locality on the NW flank of the Easter Microplate was an exception. During three dive-days (bottom-time = 14.1 hours), we observed a small vent field littered with the shells of dead bivalves. Post-expedition inspection of the dive videos identified one living adult that was not sampled; however, we inadvertently collected 21 juvenile mussels (≤ 25 mm length) that were attached to the underside of a sulfide rock sampled for geological studies. The S23 site, located in a deep basaltic crevice, was recently more active. Sulfide rubble from dead chimneys was scattered about, but we found only one active vent chimney with a high temperature of 296°C. The dominant animals were typical of the fauna found at the periphery of active ventse.g., serpulid annelids, caridean shrimp, anemones, enteropneusts, crinoids, brisingid sea stars, holothurians, sponges, and siphonophores [65]. Typical vent fauna, like alvinellid polychaetes were observed but not abundant. Two living Riftia pachyptila tubeworms were observed and one was collected [31]. Remnant tubes of Tevnia jerichonana occurred on sulfide rubble covering an area of diffuse low-temperature (≤ 8.5°C) venting, but no living worms were collected. Overall, the S23 locality exhibited hallmarks of a senescent vent field.

Geographical subdivision and admixture
Analyses were conducted with sequence polymorphisms in mitochondrial COI and the exon portions of four nuclear genes, SAHH, ANT, Cat and Col-1 (Table 4). A 9-bp indel in EF1α intron was treated as two discrete alleles. Haplotype networks for all six loci revealed differentiation between northern (N) and southern (S) mussel lineages ( Figure 2). Distinct N and S haplotypes for COI were reciprocally monophyletic and separated by at least 20 substitutions. Though less divergent, predominantly N and S alleles segregated at all five nuclear loci ( Figure 2). For each locus, the S23 sample (black pie slices, Figure 2) included N and S alleles. Latitudinal step-clines for five loci had common inflection points centered on the S23 locality ( Figure 3a). The SAHH cline was more complex (Figure 3b). Frequencies of the northern alleles, SAHH*1 and *2, declined to zero at S17, counterbalanced by increases of SAHH*7 and *8. The SAHH*4 allele only occurred in the intermediate region, S17 and S23. Indices of genic diversity differed among the genes ( Table 5). Diversity within localities (h) was greatest for COI and smallest for ANT, which was defined by one SNP, and Col-1, which was not polymorphic at the NEPR localities. Among the nuclear loci, SAHH locus had the highest diversity and richness. Expected heterozygosity and allelic richness peaked in the intergrade region, between S17 and S31, for each of the nuclear genes (Figure 4a and b). Heterozygote deficiencies (F IS ) also peaked in the intergrade zone ( Figure 4c). All statistically significant values for Fu's Fs were negative (Table 5). They were most pronounced for COI, where nine of 13 values were significantly negative. Negative Fs values reflect significant excesses of rare alleles in the simple star-like nodes for COI ( Figure 2).
Except for S23, the other samples did not exhibit evidence for non-random mating. Sporadically elevated values for F IS in these samples were not significant after sequential Bonferroni correction of the α-levels. Similarly, sporadic gametic disequilibria were not significant following correction for multiple tests. In contrast, violations of Hardy-Weinberg equilibrium were pronounced in the S23 sample. Three of the five nuclear loci exhibited significant heterozygote deficiencies (Table 6). Cytonuclear disequilibrium (CND) was tested in the S23 sample after recoding the mitochondrial haplotypes ( Figure 2) as N versus S "mitotypes". CND was significant for all five pairwise gametic combinations. Two of the tests remained significant following sequential Bonferroni correction ( Table 6). Significant gametic disequilibrium was evident in 10 pairwise comparisons of nuclear loci. Two contrasts remained significant following sequential Bonferroni correction (Table 6).
Pairwise sequence divergence within (d w ) and among (d a ) sample localities was estimated separately for each locus (Additional file 1: Table S1). The mean d w values were small for each locus: COI = 0.47%; SAHH = 0.46%; ANT = 0.01%; Cat = 0.38%; and Col-1 = 0.25%. The mean d a values were comparably low among localities within the northern (N13-S17) and southern (S31-S38) regions. Mean d a values between samples from the different regions were typically much greater than the d w values. The admixed S23 sample varied among loci in its affinities to the northern and southern regions.
These regional differences were also revealed by multilocus genotypic assignment methods. The most probable partitioning of the total sample resulted in K = 3 (Additional file 1: Table S2). The northern and southern partitions were self-evident from the previous analyses. An intermediate partition composed of the S17 and S23 was influenced by the SAHH locus, perhaps due to recombination or incomplete lineage sorting. Exclusion of SAHH from the analysis resulted in an estimate of K = 2, with S17 and S23 nested in the northern partition. The NEW-HYBRIDS analysis ( Figure 5) identified four of the S23 individuals with high posterior probabilities (PP > 0.95) as Table 4 Gene characteristics and GenBank accession numbers for six gene loci in eastern Pacific Bathymodiolus mussels    hybrids, and other potential hybrids (PP < 0.75) also were noted. Due to the limited number of loci, insufficient power existed to unequivocally assign the hybrid individuals to first-generation (F1), second-generation (F2) or specific backcross categories [60]. Several putative hybrids also occurred in the S32 sample. The S17 sample had a single individual that was assigned to the southern partition (PP > 0.95). Isolation-with-migration (IMA2) analyses were conducted with exon regions of genes that exhibited no evidence for recombination. Analyses were conducted under a two-population model, defined by the previous STRUC-TURE and NEWHYBRID analyses, a northern (N) population pooled the N13-S23 and GAR samples, and a southern (S) population pooled the S31-S38 samples. Maximum likelihood estimates (MLEs) for the effective population size scaled by mutation rate (θ = 2Nμ) were four-times larger for the northern population, as defined ( Table 7, Additional file 2: Figure S1). Pooling would limit our interpretation of these estimates, however, if significant structuring occurred within a partition. To assess this potential problem, we excluded the S23 sample from the IMA2 analysis. This exclusion reduced the estimates of effective population size by one-half (Table 7), but θ N still remained four-times greater than θ S . Decisions about pooling greatly influenced estimates of the immigration vectors (M i = m i /μ). Immigration into the northern population was greater than the reverse (M N > M S ) if the S23 sample was included in the northern population. Log likelihood ratio tests on migration rates showed exclusion of S23 from the analysis resulted in no significant bidirectional migration. Estimates of ancestral population sizes (θ A ) and population splitting time (τ = tμ) were unresolved, regardless of pooling. A sharp peak near zero in the highest posterior density (HPD) curve for τ is indicative of recent divergence; however, an escalating tail of the curve prevented an estimate for population splitting time from the multi-locus data.

Recombination
We reconstructed recombination histories from the complete sequence data (with introns) for four nuclear loci. Recombination could not be tested for ANT because it exhibited only one polymorphic site. Coalescence analyses conducted with BEAGLE [45] identified large numbers of recombination events in SAHH, Cat, Col-1 and EF1a (indicated with stars in Additional file 3: Figure S2). The northern and southern lineages were color-coded, as in Figure 1; however, homoplasy among the lineages obscured the identification of within-versus between-lineage recombination events. Nonetheless, the temporal distribution of recombination events suggests that the parental lineages may have had a long history of parapatric contact.

Isolation-by-distance
Previous studies reported that B. thermophilus populations distributed between 13°N and 17°S exhibited evidence for isolation-by-distance [29,36]. We re-tested this hypothesis with the present mitochondrial and nuclear polymorphisms. The combined pairwise F ST values for nuclear loci were not correlated with geographic distance (r = 0.008; P < 0.567, 1000 permutations) (Additional file 4: Figure  S3). In contrast, φ ST values estimated from mitochondrial sequences exhibited a significant correlation (r = 0.447; P < 0.019, 1000 permutations). Nonetheless, the φ ST correlation resulted from contrasts involving the S17 sample. Excluding S17 from the analysis eliminated the correlation. Based on the STRUCTURE and IMA2 analyses, the S17 sample appears to be influenced by introgressed southern alleles.

Taxonomic implications
Based on allozyme and mitochondrial DNA differences, Won et al. [29] suggested that the northern and southern B. thermophilus lineages might warrant recognition as distinct species. To assess this suggestion, we compared them with Bathymodiolus azoricus Cosel et al. [66] and B. puteoserpentis Cosel et al. [67], named geminatespecies that hybridize along the Mid Atlantic Ridge.
Although individuals of the MAR species can be discriminated with 95% confidence based on shell dimensions alone [24,68], the B. thermophilus lineages exhibit no comparable differences [29]. Also, sequence divergence for mitochondrial genes (COI and ND4) is greater between two MAR lineages, but differentiation at the nuclear loci is consistently less ( Table 8). The SAHH, Col-1 and Ef1α sequences examined in both species complexes are less between the MAR species (Table 8). Different numbers of allozyme loci were examined in the two complexes; so, we could only compare mean F ST values for the subsets of polymorphic genes. The proportion of allozyme variance that exists between the MAR species is half that of the B. thermophilus lineages. To summarize, shell dimensions and mitochondrial sequences indicate that the MAR species are more divergent than B. thermophilus lineages, but the opposite is true for nuclear genes. Because adaptively neutral divergence in protein-coding nuclear genes is generally believed to accumulate much more slowly than in mitochondrial genes, the consensus of evidence suggests that the EPR lineages split earlier than the MAR species. Assessing the species-status of allopatric and parapatric populations is an old problem in biology. Here we adopt the "metapopulation species concept" which equates species with separately evolving segments of the broader metapopulation [70]. In this sense, the widely distributed

Discussion
This is the second report of a mid-ocean hybrid zone involving hydrothermal vent mussels [16,29]. As in the previous case involving the Mid-Atlantic Ridge (MAR) mussels, the present example from the eastern Pacific Ocean also exhibits the properties of a Tension Zone: (1) coincident clines at multiple loci; (2) heterozygote deficiencies in the contact zone; and (3) significant gametic disequilibrium among nuclear and cytoplasmic genes (Figures 2, 3, 4, Table 6). Tension Zones are maintained by hybrid unfitness (i.e. endogenous selection) and recurrent immigration from the parental regions [1,75,76]. The abrupt step-clines at six independent gene loci suggest that introgression has been limited to adjacent localities.
We could not ascertain from the present genetic data whether the excess proportions of parental genotypes in the hybrid zone were products of selection against hybrids, or consequences of admixture (i.e. Wahlund effects). Comparison of juvenile and adult cohorts might help to identify ongoing selection against hybrids e.g., [77], but we are not aware of future efforts to explore this remote region with submersible vehicles. Simple admixture of immigrant recruits from the adjacent parental regions would generate the observed heterozygote deficiencies, gametic phase disequilibrium [78], and putative hybrid genotypes in the S23 sample ( Figure 5). Although the S23 vent field appeared to be in a senescent stage of the hydrothermal cycle a single active vent chimney supported small populations of chemosynthetically dependent animals. We observed one living adult mussel (unsampled) while reviewing the dive videos from 14 hrs of Alvin 'bottom-time'. Shell debris from recently deceased bivalves littered the floor of the  crevice that contained the active chimney. Other vents hosting more robust mussel populations might exist nearby, but they remain undiscovered. "Sweepstakes Reproductive Success" (SRS) as observed in some near-shore bivalves and other organisms reviewed in [42], might be invoked to explain the heterozygote deficiencies, but the hybrid zone juveniles were not the progeny of few females. The sample of 21 juveniles contained 11 distinct mitochondrial haplotypes (Table 5). Haplotypic diversity (h = 0.819) was not different from the mean diversity (h = 0.807 ± 0.077) for all other samples, which exhibited no evidence for deviations from random-mating expectations. The planktotrophic larvae of Bathymodiolus mussels appear to mix thoroughly before they settle at vents. Their capacity for long-distance dispersal is consistent with a lack of geographical subdivision among EPR (N13-S17) and GAR populations spanning thousands of kilometers [79] and this study. A previous report of subdivision across the Equator is an artifact of sampling gaps cf. [21,36]. Nonetheless, broad-scale homogeneity could coincide with small-scale patchiness in genotypic frequencies under the SRS model [42]. Hypervariable genetic markers might reveal localized patchiness due to episodic recruitment of cohorts from a limited number of parents in different source populations. For example, AFLP variation provides preliminary evidence for small-scale patchiness in the eastern Pacific vent annelid Riftia pachyptila [80], despite broad-scale homogeneity of nuclear and mitochondrial DNA markers [31]. In contrast, microsatellite DNAs provided no evidence for small-scale patchiness in the western Pacific vent snail, Ifremeria nautileii [81].

Divergence and intergradation
Reciprocally monophyletic differences for mitochondrial haplotypes and step-clines in allozyme frequencies suggest a relatively long history of isolation between B. antarcticus n. sp. and B. thermophilus s.s. Based on borrowed molecular clock rates for COI (1-2% divergence per million years), we estimated that the northern and southern mitochondrial lineages split 2.1-4.3 Myr ago, which temporally overlaps a span of dates encompassing orogeny of the Easter Microplate, 2.5-5.3 Myr ago [82,83]. An ongoing study that uses fossil data to calibrate major splitting events in the subfamily Bathymodiolinae estimates a similar time for the most recent common ancestor of the mitochondrial lineages (J. Lorien, personal communication). Nonetheless, mitochondrial 'gene trees' do not necessarily equal 'species trees' [84]. Different mutation rates among nuclear and mitochondrial genes, the amounts of ancestral polymorphism in various genes, and differing intensities of selection on independent characters can greatly complicate relationships between gene and organismic trees. Disagreements between these trees are further complicated by introgressive hybridization and reticulation [85]. Thus, it is not surprising that IMA2 analysis failed to resolve a splitting time (τ) from the multi-locus data for these mussels.
Whether intergradation between the northern and southern lineages results from primary or secondary processes remains uncertain. A coalescent analysis of the MAR mussel species provided evidence for "an initial period of allopatric differentiation during which recombination was blocked between lineages" followed by "a long history of gene flow" [25]. In contrast, our analysis of past recombination in the eastern Pacific mussels is consistent with a long history of parapatric contact. However, these contacts probably are episodic in this dynamic region of the global ridge system. Chaotic local extinctions and colonizations of vents as tectonic activities move up and down a ridge axis can create transient contact zones that allow gene flow [86]. The observed recombination history of these mussel lineages is consistent with episodic isolation and re-contact as tectonic conditions vary temporally.

Isolating mechanisms
Density troughs (regions of low population size) can function as powerful traps to stabilize hybrid zones and limit introgression [87][88][89]. The MAR hybrid zone coincides with a density trough at 'Broken Spur' [16,24], a geographically intermediate vent field that is relatively inhospitable for mussels due to precipitation of toxic polymetalic sulfides believed to interfere with respiration [23]. Mussel population densities were also low in the Easter microplate region during our 1999 and 2005 expeditions, but the causes are different. First, Won et al. [29] hypothesized that strong cross-axis currents will remove planktotrophic mussel larvae from the ridge axis, limiting local self-recruitment. These superfast-spreading centers have highly inflated axial calderas and essentially no lateral walls to constrain axial circulation [90]. The predominantly cross-axis currents in this region drive buoyant hydrothermal plumes to the west [91] presumably carrying larvae that rise above the axial calderas. Secondly, Coykendall et al. [31] hypothesized that vent habitats along the southern EPR are subject to high rates of disturbance. The habitats are frequently "repaved" with lava flows that drive local extinctions and create new vents [74]. We explored a number of hydrothermally active vent fields between 18.40°and 21.73°S latitude during our January 1999 expedition (Table 1), but none of them supported robust mussel populations. Soon after, Van Dover [92] visited the ' Animal Farm' vent field at 18.60°S and found it "clearly in a waning stage of the hydrothermal cycle' (p. 142). Venting water temperatures were only slightly elevated above ambient and most of the mussels were dead. The living individuals were smaller on average compared to dense mussel populations just to the north ( Table 1). The BIOSPEEDO expedition in 2004 sampled mussels at 18.40°and 21.42°S, but population densities were not reported [36]. Demographic fluctuations result in variance-effective population sizes that might be very small for some taxa, leading to losses of genetic variation [22]. For example, southern EPR populations of Riftia pachyptila are nearly devoid of nuclear and mitochondrial sequence variation [31], but other annelids are not similarly affected [30,36]. Together, cross-axis currents and metapopulation processes might reduce the genetically effective population size of mussel populations in this region, impeding along-axis dispersal and limiting introgression across the Easter Microplate.

Conclusions
Although B. thermophilus s.s. and B. antarcticus n. sp. cannot be distinguished with shell measurements, they differ in having reciprocally monophyletic mitochondrial lineages. Polymorphisms at nuclear loci (allozymes and DNA sequences) are not completely diagnostic but the two species comprise discrete genotypic clusters with predominantly non-overlapping geographical distributions. Some shared alleles might be remnants of incomplete lineage sorting between the metapopulation lineages, but these alleles are not widely distributed geographically, as expected for ancestral polymorphisms. Instead, they only penetrate localities that flank the Easter Microplate, a pattern expected for introgressive hybridization or disruptive selection along an environmental gradient [1].
Divergence between the two species probably initiated with orogeny of the Easter Microplate, 2.5-5.3 Myr ago. Superfast tectonic spreading rates in this dynamic part of the global mid-ocean ridge system result in frequent habitat destruction, leaving episodic gaps in the distribution of many vent-restricted species. Prolonged periods of partial to complete isolation would facilitate the accumulation of genetic differences in nuclear and mitochondrial genes exposed to genetic drift and potentially contrasting selection regimes. We suggest two hypotheses about possible adaptions that might result in disruptive selection across the Easter Microplate boundary. First, we hypothesize that novel life history traits emerged to accommodate different current regimes in the northern and southern regions. Northern B. thermophilus s.s. populations are subjected to a circulation pattern dominated by mesoscale eddies and along-axis currents [20,[93][94][95]. There, larval buoyancy and duration probably evolved to favor the retention of the planktotrophic prodissoconch stages. South of the Easter Microplate, B. antarcticus n. sp. is subjected to strong cross-axis vectors forced by the Antarctic Circumpolar Current [29]. Larval buoyancies and durations that function in the northern current regime would result in larval removal in this contrasting regime [96]. Though we currently have no developmental data on the larval life histories of these mussels, this hypothesis is testable. Early settling post-larvae retain their developmental history in prodissoconch sizes and morphologies that can provide significant clues regarding larval life history [97].
Secondly, we hypothesize that coadaptation between the mussel hosts and associated thiotrophic endosymbionts might result in disruptive selection. B. thermophilus s.s. and B. antarcticus n. sp. differ with respect to the symbiont strains they host in specialized cells (bacteriocytes) contained mostly in gill tissues [98]. Though we cannot at this time exclude some degree of vertical transmission, the endosymbionts are mostly acquired locally from the environment in which mussel larval settle [99]. An ongoing multi-gene analysis of symbiont diversity has only found the B. thermophilus s.s. symbiont strain in the hybrid zone mussels (Y-J Won, unpublished data), even in host individuals with predominantly southern genotypes. Subtle interactions involving these intracellular microbes might have affected evolution in mitochondrial and nuclear genes of the hosts. According to the "Red King" model, infectious horizontal symbionts and their hosts are expected to exhibit decelerated evolutionary rates, because the participants are subject to constraints that limit any form of change [100,101]. Episodic range expansions into the contact zone will result in secondary contacts that reestablish opportunities for gene flow and introgression in the hosts, and possible symbiont transfers. Yet, conflicts between divergent larval strategies, symbiont coadaptation, and other unknown processes would determine the strength of selection against intergrades. Disruptive processes involving the mussel hosts and their obligate endosymbionts might contribute to partial isolation and incipient speciation of these mussels.

Additional files
Additional file 1: Table S1. Mean percentage of sequence divergence (K2P corrected) within (d w , boldface on diagonal) and between (d a off-diagonal) eastern Pacific Bathymodiolus mussels. Pairwise comparisons involving the S23 sample are highlighted in blue. The d values between northern and southern regions are highlighted in gray. Table S2. Simulation results for the STRUCTURE analysis of eastern Pacific Bathymodiolus mussels under correlated and uncorrelated allelic frequency models. The number of possible genotypic clusters, K, varied from one to eleven. The natural log of the probability of the data for a given value of K was averaged across the number of simulations (runs) per K. Bayes factors were estimated based on these averages.