Skip to main content
  • Research article
  • Open access
  • Published:

Genetic diversity and demographic instability in Riftia pachyptilatubeworms from eastern Pacific hydrothermal vents

Abstract

Background

Deep-sea hydrothermal vent animals occupy patchy and ephemeral habitats supported by chemosynthetic primary production. Volcanic and tectonic activities controlling the turnover of these habitats contribute to demographic instability that erodes genetic variation within and among colonies of these animals. We examined DNA sequences from one mitochondrial and three nuclear gene loci to assess genetic diversity in the siboglinid tubeworm, Riftia pachyptila, a widely distributed constituent of vents along the East Pacific Rise and Galápagos Rift.

Results

Genetic differentiation (F ST ) among populations increased with geographical distances, as expected under a linear stepping-stone model of dispersal. Low levels of DNA sequence diversity occurred at all four loci, allowing us to exclude the hypothesis that an idiosyncratic selective sweep eliminated mitochondrial diversity alone. Total gene diversity declined with tectonic spreading rates. The southernmost populations, which are subjected to superfast spreading rates and high probabilities of extinction, are relatively homogenous genetically.

Conclusions

Compared to other vent species, DNA sequence diversity is extremely low in R. pachyptila. Though its dispersal abilities appear to be effective, the low diversity, particularly in southern hemisphere populations, is consistent with frequent local extinction and (re)colonization events.

Background

Demographic instability influences the distribution of genetic diversity within and among discrete colonies of a widely distributed species [1, 2]. Local extinction and (re)colonization events are expected to reduce the overall effective size of a metapopulation thereby limiting its capacity to retain genetic variation [3]. Furthermore, genetic drift contributes significantly less to geographical differentiation if the average persistence time of colonies is less than the time it takes to fix neutral alleles [4]. The number of colonists, their sources (i.e. migrant vs. propagule pools), and ongoing rates of gene flow all affect the degree to which demographic instability affects geographical differentiation [57]. Yet, these parameters are rarely known; instead they are commonly inferred from statistical analyses of genealogical variation that is assumed to be selectively neutral [8]. For a variety of practical reasons, recent studies have relied mostly on information from mitochondrial DNA variation, although the assumption of selective neutrality has been debated [911]. Independent genealogies involving nuclear genes are needed to distinguish the effects of idiosyncratic processes such as mitochondrial selective sweeps from demographic processes that leave genome-wide signatures [12, 13].

The communities of marine animals that occupy deep-sea hydrothermal vents provide unusual opportunities to study dynamic metapopulation processes [1416]. Vent communities are supported by chemosynthetic bacteria that oxidize volcanic gases (e.g., H2S and CH4) concentrated in the hydrothermal effluents. Adults of the iconic tubeworm Riftia pachyptila (Polychaeta: Siboglinidae), for example, are nourished entirely by sulfur-oxidizing endosymbiotic bacteria, linking their lives to the tempo of disturbances that create and destroy their habitats (Figure 1A-B). Eastern Pacific vents (Figure 1C) that host R. pachyptila are particularly ephemeral, persisting for a few years to several decades before fluid conduits are blocked, magma supplies shift, or lava flows extirpate local communities [17]. Earthquakes can open fluid conduits, re-activating old vents, and magmatic eruptions spawn new vents [18]. Studies of vent community succession [18] revealed that R. pachyptila is among the first species to colonize a new vent once suitable conditions are established. Within two years its numbers can grow to several thousand adult individuals, but changes in vent flow, or overgrowth by mytilids can lead to its replacement as the dominant species -- a process that can take months to years depending on the location. The frequency of habitat turnover varies with tectonic spreading rates [19]. Tectonic spreading at Eastern Pacific vents (Figure 1C) varies from moderate rates (65 mm/yr) along the Galápagos Rift (GAR) to fast rates (85-116 mm/yr) along the northern East Pacific Rise (NEPR) and superfast rates (142-158 mm/yr) along the southern EPR (SEPR). Cycles of habitat disturbance associated with these events are expected to favor organisms with rapid individual growth rates, early reproduction, and effective dispersal capabilities [20].

Figure 1
figure 1

Habitat instability and genetic diversity in R. pachyptila. (A) A healthy patch of tubeworms at the N27 locality. (B) An adjacent senescent patch on a rust-colored sulfide mound covered with numerous scavengers, the galatheid squat lobster Munidopsis subsquamosa. (C) Riftia pachyptila samples: blue and red dots indicate northern and southern sample locations; gray dots indicate active vents known to host R. pachyptila but not included in present analyses; white dots denote active vents that did not host substantial R. pachyptila colonies during the times of our expeditions; and white diamonds denote vents that did not support R. pachyptila colonies. Tectonic spreading rates are indicated along arrows for each location. (D) Allelic frequencies at four loci. Colors are coded to adjacent haplotype networks (black wedges are singletons). (E) Structure plots for the mean probability of assignment of individuals (Q values) to the northern cluster (blue) versus the southern cluster (white). (F) Haplotype networks for six genetic markers. Area of circle is proportional to the frequency of each haplotype, and straight lines denote single nucleotide differences.

Active vents along the EPR and GAR are restricted to axial grabens (rift valleys) that are expected to function as dispersal corridors for the larvae of species like R. pachyptila [21]. The distances between active vents vary from a few kilometers along a ridge segment to 100s of kilometers between contiguous segments (Figure 1C). A one-dimensional stepping-stone model [22] provides a reasonable starting hypothesis for connectivity among linear habitats such as these, but population genetic studies provide limited support for this configuration, partly due to under-sampling of populations and insufficiently polymorphic genetic markers [23]. However, R. pachyptila was first among the few species that exhibited evidence for stepping-stone dispersal [24]. Estimates of gene flow (Nm), inferred from pairwise F ST 's for allozyme differentiation, declined with geographical distances among NEPR and GAR populations. Based on larval durations and predominant currents along the NEPR, R. pachyptila is expected to have mean larval dispersal distances of 100 to 200 km [17, 21]. A subsequent phylogeographic survey examined additional samples from the SEPR axis, but found no evidence for stepping-stone dispersal due to a near absence of sequence variation in mitochondrial Cytochrome oxidase I, (COI) [25]. A number of co-distributed species of annelids and mollusks have been similarly examined, and they all exhibit substantial COI variation [2527]. Their mitochondrial haplotype diversities range from a minimum of 0.331 in R. pachyptila to a maximum of 0.960 in another annelid, with a mean of 0.615 across 11 species [summarized in Figure five of reference 15]. Examinations of nuclear allozymes, AFLPs and DNA microsatellites from limited portions of R. pachyptila's range [24, 28, 29] suggests that its low COI variation may be anomalous, resulting perhaps from a mitochondrial selective sweep. Nevertheless, allozymes, AFLPs and microsatellites exhibit mutation rates that do not compare to mitochondrial DNA sequences, and they might experience different selective constraints. A study including nuclear DNA sequences with sufficient numbers of individuals sampled from R. pachyptila's known geographic range has yet to be performed.

Testing selective sweep versus population bottleneck hypotheses requires an examination of multiple independent genetic markers [13]. To that end, we examined DNA sequences from three independent nuclear loci and another mitochondrial locus, Cytochrome b (Cytb), from R. pachyptila. Samples were collected from known vents along the GAR, NEPR and SEPR axes (Figure 1C). First, we tested the selective sweep hypothesis by assessing whether mitochondrial diversity within and among populations is idiosyncratically low. Second, we examined whether population structure is consistent with a linear stepping stone model, or if more complex models of subdivision are justified. Finally, we assessed whether demographic instability associated with tectonic spreading is correlated with reduced levels of genetic diversity.

Methods

Samples

Expeditions conducted between 1998 and 2005 explored and sampled hydrothermal vents at 19 locations along the NEPR, SEPR and GAR axes (Figure 1C). The samples examined in this study were obtained from nine locations (Table 1). Worms were harvested with robotic manipulators on the human occupied vehicle (HOV) Alvin (Woods Hole Oceanographic Institution, WHOI) or the remotely operated vehicle (ROV) Tiburon (Monterey Bay Aquarium Research Institute, MBARI). Samples were placed in an insulated container with ambient seawater. Upon recovery at the surface, whole worms were stored in 2°C filtered seawater prior to tissue sampling. A small sample of vestimentum, a tissue that is normally devoid of endosymbiotic bacteria, was frozen at -70°C or preserved in 95% EtOH.

Table 1 Coordinates and depths for Riftia pachyptil a samples.

DNA purification, PCR conditions and DNA sequencing

Genomic DNA was extracted using the Qiagen DNeasy kit, following manufacturer's protocols (Qiagen, Inc. Valencia, CA). The primers, optimal annealing temperatures, and PCR conditions are available in Additional file 1. We used two PCR programs: (i) initial denaturation temperature of 94°C for 10 min followed by 35 cycles of 94°C for 1 min, 55°C for 1 min and 72°C for 1 min with a single final extension of 72°C for 7 min; and (ii) initial denaturation temperature of 94°C for 10 min followed by 35 cycles of 94°C for 2 min, 48°C for 2 min, 72°C for 2.5 min, and final extension of 72°C for 7 min. PCR products from all loci were purified by gel excision and cleaned with Montage filter units (Millipore Corp., Billerica, MA) or diluted in 50 μl sterile water and purified with Multiscreen HTS PCR 96 vacuum manifold system (Millipore Corp., Billerica, MA). Unless noted (Additional file 1), DNA sequencing employed the same PCR primers. Purified PCR products were sequenced bidirectionally on an ABI 3100 capillary sequencer using BigDye terminator chemistry (Applied Biosystems Inc., Foster, CA). Protein-coding sequences were checked for stop codons.

Mitochondrial sequences

PCR mixtures for Cytb included 2.5 mM MgCl2, 0.4 μM of each primer, 200 μM of each dNTP, 1.25 units of Taq polymerase (Applied Biosystems Inc., Foster, CA), 1x PCR buffer (supplied by manufacturer), and 1 μL of DNA template in a final volume of 25 μL. Cytb sequences were amplified with previously published primer pairs [30] and PCR program-i (Additional file 1). The 377 bp amplified portion of Cytb corresponds to positions 10732-11108 in the R. pachyptila mitochondrial genome [GenBank acc. no. AY741662, deposited by 31].

Single-copy nuclear (scn) DNAs

One of us (SAK) previously generated a series of anonymous scnDNA sequences from R. pachyptila (unpublished, GenBank acc. nos. U68732-U68750). We designed primers from these sequences following Karl et al. [methodological details in 32]. Briefly, purified total DNA from a single R. pachyptila individual from the GAR was digested with the restriction enzyme Sau3a. Fragments were sorted by size, cloned, and screened to determine genomic copy number. To design primers, we sequenced each end of the low copy clone-inserts. Nine clone-inserts were tested for their potential utility in population screening. Seven inserts generated products that could not be used because they included length variation (large indels) that confounded our scoring of individual genotypes, or they lacked sufficient polymorphism. The Rpt46.1 primers amplified predictable products, but internal primers were required for Rpt84.1 to improve its sequence quality (Additional file 1). PCR program-i was used for both loci. Amplification with the Rpt46.1 primers generated two products; the larger band was used as template for sequencing. After preliminary sequencing, many genotypes contained two polymorphic sites. To resolve the ambiguous haplotypes, we cloned the PCR products from heterozygous individuals using the standard manufacturer's protocol (Topo TA cloning kit; Invitrogen, Carlsbad, CA).

Nuclear intron

Of seven primer pairs developed by Jarman et al. [33], only the ATPSα primers amplified a polymorphic product that could be scored in R. pachyptila. To improve its reliability, we designed non-degenerate internal primers (Additional file 1). The intron amplified using the same PCR recipe as above and PCR program-ii. Initial screening of the 473 bp product from ten individuals from nine populations revealed a single nucleotide polymorphism (A/G) at position 319 of the amplicon. Because the SNP could be recognized by the restriction enzyme NspI, population screening was based on restriction fragment length polymorphisms (RFLPs). Ten μl of the amplicons were suspended in a cocktail containing 1 × NE buffer-2, 10 μg/μl bovine serum albumen, 3 units of NspI (New England Biolabs, Ipswich, MA) and sterile water to a final volume of 15 μl, and the mix was incubated at 37°C for 1.5-11 hours. Products were separated on 3% agarose gels at room temperature for 1.5-3 hrs at 65-80 volts. Positive controls of known genotypes (from prior sequencing) were included with RFLP digestions and all gel analyses.

Statistical analyses

We used the program Phase v. 2.1 [34, 35] to reconstruct allelic haplotypes for nuclear loci using the cloned haplotypes from the double heterozygotes as a baseline. We also used the software to estimate background recombination rates per gene segment. Recombination rates were low enough to assume no effective recombination within each segment. Tcs v. 1.21 [36] was used to construct statistical parsimony networks for each locus. Arlequin 3.5 [37] was also used to estimate molecular diversity, θ S [38], from the number of segregating sites, and Fu's Fs [39], and to estimate gene diversity (H), adjusted for sample size according to Eq. 8.4 of Nei [40]. Allelic richness (a') estimates were obtained by the rarefaction methods implemented in Hp-rare v. 1.0 [41]. Estimates of F IS and tests for Hardy Weinberg equilibrium (HWE) and linkage disequilibrium (LD) were conducted with Genepop v. 3.4 [42].

The multi-locus genotypic assignment method employed in Structure v. 2.2 [43] was used to estimate the probable number of discrete population clusters (K) represented by the data. We assumed correlated allele frequencies between population samples and admixed ancestries of populations. Parameter estimates for correlated allele frequencies and admixture were estimated by averaging values from preliminary simulations that varied K from 1 to 8. Simulations were repeated at least five times for each value of K (number of data partitions) and for two disparate values of λ (0.731, 2.204), a shape distribution parameter for allele frequencies, based on initial simulations. 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. All simulations employed a burn-in value of 105 with 106 subsequent steps in the Markov Chain Monte Carlo (MCMC) sampling process.

To assess evidence for isolation-by distance (IBD) we used Arlequin 3.5 to obtain pairwise estimates of F ST , and the GeoDistances procedure in Le Progiciel R [44] to estimate great arc distances from geographical coordinates for the sample locations (Table 1). The Mantel test module implemented in Arlequin 3.5 was used to assess correlations between the genetic (A) and geographical (B) distance matrices. To assess possible affects of subdivision, we generated an assignment matrix (C) with values that equaled zero (0) if populations belonged to a cluster identified by Structure and one (1) if populations belonged to different clusters. The partial correlation (rAB|C) of A and B after accounting for C was obtained according to the methods of Smouse et al. [45] as implemented in Arlequin 3.5.

The F ST -outlier method [46] implemented in the program Lositan [47] was used to test for evidence of selection. The program uses empirical estimates of genic heterozygosity to simulate neutral expectations for F ST under an island model of migration. An iterative bisection algorithm is used to simulate data with a mean F ST that matches the empirical data, even when the data do not strictly adhere to island model assumptions [i.e. F ST ≈ 1/(4Nm + 1)] for neutral genes [47]. Initial simulations were conducted with 10,000 replications, assuming an infinite allele mutation model. We then conducted simulations with 40,000 replications and used the neutral mean for F ST to estimate its 95% confidence intervals (CIs). Empirical values were compared to this distribution and those significantly above or below the 95% CIs might be influenced by selection. We attempted to use the program Sweep_Bott [13] to distinguish population bottlenecks from a selective sweeps, but the simple star phylogenies obtained for the present genes prevented its execution.

Geographical and geophysical methods

Seafloor spreading rates (Figure 1C) were estimated by the NUVEL-1A model implemented in the web-based calculator < http://www.ldeo.columbia.edu/~menke/plates2.html>. This kinematic plate-motion model uses tables of Euler vectors previously estimated from geomagnetic reversals [48]. Space-based geodetic methods were used to improve accuracy of the model [49]. The NUVEL-1A estimates of spreading rates are highly correlated (r = 0.982) with empirical measurements based on near-bottom surveys and space geodesy methods [50].

Results

Distribution and abundance

We explored 19 active vent fields between 27°N and 38°S latitude and obtained adequate samples of R. pachyptila from eight localities (Figure 1C, blue and red dots). Robust populations were observed and sampled from the NEPR localities. Northern SEPR localities S7, S11 and S17 also host robust R. pachyptila populations, but only one specimen was observed and sampled from S23 on the NW flank of the Easter Microplate. Nine worms were observed and sampled from S32 during four dives in 1999, and none were found when we returned there in 2005. We conducted five submersible dives at S38 and found a single individual living deeply within a recess on the vertical wall of an active chimney, but it could not be sampled before deteriorating weather conditions forced us to abandon this location.

Polymorphism

Twenty-five alleles were identified across four polymorphic loci (Additional file 2). Sequences associated with these alleles were deposited with GenBank (acc. Nos. HQ405793-6724). For each locus, parsimony analyses produced very simple star-like genealogies with mostly one-step mutational differences from a primary (most frequent) allele (Figure 1D and 1F). Widely distributed secondary (less frequent) alleles occurred at three loci. The assignments of haplotype phases in double heterozygotes were resolved with 97% probability for Rpt46.1 (5 individuals) and 89-96% probability for Rpt84.1 (6 individuals) depending on the sample location. Haplotypes inferred with the program Phase were used in all subsequent analyses. Sequence diversity was very low in the samples for all four loci, with π values ranging from 0 to 0.006 and θ S values ranging from 0 to 0.004 (population values not reported).

Genotypic frequencies at the nuclear loci revealed no evidence for deviations from Hardy-Weinberg equilibrium (F IS , Table 2). Multi-locus tests obtained by summing P-values for each locus (Fisher's method) provided no indication of significant heterozygote deficiencies within population samples (8 tests, Prange: 0.169-1.000). Furthermore, the grand total of these values obtained by summing across all loci and samples was not significant (P = 0.320). Tests for linkage disequilibrium between pairs of polymorphic nuclear loci provided no evidence for associations within populations (16 tests, Prange: 0.085-1.000). Furthermore no associations existed when P-values from the pairwise comparisons were summed across populations (3 tests, Prange: 0.484-0.971). Similarly, tests for cytonuclear disequilibrium involving Cytb and the three nuclear genes were non-significant within populations (13 tests, Prange: 0.158-1.000) and globally across populations (3 tests, Prange: 0.370-0.931).

Table 2 Indices of genic diversity in Riftia pachyptila samples.

Test of selection

Because overall differentiation, as estimated by F ST , varied among the four loci (0.163 for Rpt46.1; 0.188 for Rpt84.1; 0.269 for ATPSa; and 0.551 for Cytb), we used Lositan to test for evidence of natural selection. First, we determined the neutral expectation, F ST (nuc) = 0.138, for the nuclear loci alone. Though all the observed F ST (nuc) values were higher than the neutral mean, none fell outside of the 95% CIs. Assuming an even sex ratio and strictly maternal cytoplasmic inheritance, the neutral expectation for a mitochondrial gene is F ST (mt) = 4F ST (nuc)/(1 + 3F ST (nuc)) [51, 52]. The F ST (mt) = 0.551 for Cytb was only slightly greater than its neutral expectation (0.533); thus, no evidence existed for directional or balancing selection having shaped geographical differentiation for these genetic markers.

Geographical structure

We used Structure to estimate the number of population clusters (K) represented among 192 multilocus genotypes. For each combination of priors, K = 2 provided the greatest mean value for ln Pr(X|K). Bayes factors were 0.368 for K = 2 and zero for K = 1 and K = 3. We estimated Q, the probability of assignment of genotypes to the northern cluster (NEPR/GAR), and represented the mean of Q for each sample in Figure 1E. The means exhibited a clinal pattern between N27 and GAR, and are effectively homogeneous along the SEPR axis. A single individual sampled from S23 clustered with its neighboring samples. Hierarchical AMOVA of the eight population samples partitioned 18.8% of gene diversity between the northern and southern clusters, 9.3% within clusters, and 71.9% within samples.

Evidence also existed for isolation-by-distance, IBD (Figure 2). Pairwise F ST 's (Table 3) increased significantly with geographical distances (Mantel's r = 0.709, P = 0.003); however, the comparisons between northern and southern clusters greatly influenced this relationship (black dots in Figure 2). An IBD relationship remained after exclusion of SEPR samples from the analysis (r = 0.896, P = 0.167), but the correlation was not significant. To better assess how associations with the northern and southern clusters (C) affected the relationship between F ST (A) and distance (B), we estimated the partial correlation of A on B after accounting for C. A marginally significant relationship remained (rAB|C = 0.396, P = 0.069). It should be noted that the trans-Equatorial F ST between N9 and S7 (0.045) is smaller than the F ST 's between adjacent pairs of NEPR populations (range: 0.064-0.093). This trans-Equatorial F ST also is smaller than F ST 's between other population pairs that are about ~2000 km apart (N27 vs. N9; N9 vs. GAR; S7 vs. GAR) (Table 3).

Figure 2
figure 2

Correlation of genetic differentiation with geographic distances among eight sample localities. Negative values for F ST (Table 3) were set to zero prior to adjustment. Black dots denote contrasts between populations from the northern and southern clusters. White dots denote contrasts within clusters.

Table 3 Pairwise estimates of F S T (lower triangle) and geographic distances in km (upper triangle).

Environmental correlates of gene diversity

Heterozygosity declined strongly with tectonic spreading rate (Figure 3, r = -0.888), whether or not the off-axis GAR sample was included in the analysis (Table 4). Spreading rates are negatively correlated with latitude (r = -0.78) and the relationship is stronger with removal of the off-axis GAR sample (r = -0.98). Depth is relatively monotonous along the EPR axis and uncorrelated with the gene diversity indices. These correlations are clearly reflected in differences in the gene diversities between northern (NEPR/GAR) and southern (SEPR) populations, treated individually or as population clusters (Tables 2 and 5).

Figure 3
figure 3

Correlation of expected heterozygosity with tectonic spreading rates ( R2= 0.788, P = 0.0032). Allelic richness exhibits essentially the same relationship (R2= 0.730, P = 0.0069). Black squares are NEPR populations, white squares are SEPR samples and the white circle is the GAR population.

Table 4 Correlations between geographical and genetic measures for all eight samples including GAR (below diagonal) and for the seven EPR samples alone (above diagonal).
Table 5 Genic diversity in northern and southern population clusters.

Demographic tests were initially conducted separately for each locus at each locality, but limited polymorphism in each sample provided essentially no power to reject null hypotheses. Subsequent tests were conducted for the northern and southern population clusters (Table 5). Five out of six values for Fu's Fs were negative, though only one was statistically significant. Mismatch distributions were all unimodal, suggesting an expanding population size, but the tests of demographic stability had very limited power due to low polymorphism.

Discussion

Riftia pachyptila exhibits the lowest mitochondrial COI diversity observed to date among 14 species of invertebrates sampled from eastern Pacific hydrothermal vents [reviewed in 15]. Yet, the present data were not consistent with the hypothesis that a selective sweep eliminated mitochondrial diversity alone. Levels of sequence diversity were similarly low for three nuclear genes and another mitochondrial gene, Cytb. The shallow star-like haplotype networks exhibited by all four loci suggest that demographic factors affecting the entire genome were responsible for the reduced variation. Gene diversity indices (H and a') for eight population samples were inversely correlated with tectonic spreading rates (Figure 3), which can be taken as surrogates for rates of habitat turnover. The lowest gene diversities occurred at southern (SEPR) localities, which are subjected to superfast spreading rates of 142-156 mm/yr. Active vents are closely spaced along the southern SEPR axis (S14-S23), but R. pachyptila colonies were spotty there and individuals were scarce. In contrast, higher gene diversities occurred among the northern (NEPR/GAR) localities, where spreading rates are slower, 65-116 mm/yr. Active vents are more distantly spaced along the NEPR axis, but dense R. pachyptila colonies occupied nearly all the active vents that we have explored along the NEPR during the past 22 years. These observations are consistent with evidence from other vent species that gene diversity is correlated with vent "occupancy" [14, 15]. Species that occupy most of the active vent habitats in a region have high gene diversities compared to species with spottier distributions. Here we have found the same relationship within a species. Northern R. pachyptila have higher vent occupancy and gene diversity, whereas southern populations have low occupancy and low diversity. Frequent cycles of extinction and (re)colonization relative to time it takes to fix new mutations are expected to produce spotty distributions and correspondingly low gene diversity [4]. Apparently, high rates of population turnover along the southern SEPR axis has led to rapid coalescence on a small number of founders, but the origin of these founders remains a mystery.

Population structure

Overall, differentiation among R. pachyptila populations meets expectations for a linear stepping-stone model. Pairwise estimates of differentiation, F ST /(1 - F ST ), increased with geographical distances among the vent localities (Figure 2), particularly among the NEPR/GAR populations as previously reported for allozymes [24]. Large distances (≥2000 km) between the northern (NEPR/GAR) and southern (SEPR) population clusters greatly influenced the overall pattern, however. Hierarchical AMOVA attributed 18.8% of total gene diversity to differences between the northern and southern clusters, but this evidence for subdivision may be an artifact of sampling gaps between the two regions. Computer simulations have revealed that sampling gaps in a linear stepping-stone model can generate inflated pairwise F ST 's that may lead to incorrect inferences about population subdivision [23]. Here, the trans-equatorial F ST = 0.045 between N9 and S7 (1957 km apart, Table 3) is smaller than the F ST = 0.271 between the N9 and N27 along the NEPR (2050 km apart), or between N9 and GAR (2234 km apart). Even the F ST 's between closer pairs of NEPR samples (N27-N21 and N21-N9) were greater than the trans-equatorial pair. Estimates of pairwise migration rates (Nm) between adjacent EPR samples were similar, ranging from 2.0 to 5.3 after excluding the three southern SEPR samples for which F ST ≈ 0. Nonetheless, the existence of regional alleles (blue and red in Figure 1D) suggests that dispersal across the Equator might not have achieved equilibrium with genetic drift (discussed below).

The NEPR, SEPR and GAR axes meet at the Galápagos Triple Junction (Figure 1C), which plunges to 5400 m in the 50-km wide Hess Deep. Active vents are not known from the Hess Deep proper, but video surveys of the adjacent Hess Deep Rift Valley have identified vents that host R. pachyptila, though they were not sampled (T. Shank 2009, pers. comm.). Further exploration of the equatorial region is warranted because we do not know why it creates a variable dispersal filter for several vent taxa [reviewed in 15]. Hurtado et al. [25] first noted that several annelid species exhibit variable impedance to dispersal across the Equator, and Plouviez et al. [26] reported that one or more historical vicariance events might have separated NEPR and SEPR populations. The palm worm, Alvinella pompejana, and limpets of the Lepetodrilus elevatus species complex exhibit evidence for isolation across the Equator, but the elevated φ ST 's reported by Plouviez et al. for other species may not be greater than expected from the sampling gaps between these axes [23]. Variable impedance across such regions results from interactions between taxon-specific dispersal modes and physical factors [15]. For example, strong deep-ocean currents sweep eastward across the EPR axis near the Equator [53]. The removal of dispersing larvae from axial troughs by these currents and the apparent scarcity of vent habitats may together impede the dispersal of some species.

Non-equilibrium processes

Inferences about isolation-by-distance or stepping-stone models assume the constituent populations have achieved equilibrium between gene flow and genetic drift, but non-equilibrium scenarios can generate similar patterns of population structure [54]. If the Equatorial region is a contact zone between partially isolated NEPR and SEPR populations, this too could produce differentiation that resembles the product of stepping-stone dispersal, but this hypothesis fails to explain the low gene diversity in SEPR populations of R. pachyptila (discussed below). Magma supplies and associated tectonic and hydrothermal activities shift along ridge axes; thus, populations may be variably interconnected and disconnected by such shifts [17, 55]. Numerous closely spaced vents occur between 14°and 21°S latitude (diamonds in Figure 1C), and they are very young as evidenced by recent lava flows [56, 57]. Many of these vents host colonies of polychaetes, gastropods and bivalve mollusks [26, 58, 59], but R. pachyptila was rare during the expeditions we mounted between 1999 and 2005 to vents between 14°and 21°S along the SEPR. Genetic homogeneity among SEPR samples spanning 2673 km suggests that genetic drift contributes very little to differentiation among these colonies, as expected if their average persistence time is shorter than the time it takes to fix neutral alleles [4]. These ephemeral vents may offer no more than a demographic sink for R. pachyptila. Perhaps this region has been repeatedly colonized from more stable R. pachyptila sources that remain undiscovered. We searched for this species as far as 38°S latitude along a contiguous extension of the Pacific-Antarctic Ridge (PAR, Figure 1C). This latitude near the "roaring 40's" is the farthest south anyone has conducted operations with a human occupied research submersible. During five dives, we observed only a single R. pachyptila individual that could not be sampled before foul weather drove us from the location. If stable source populations exist elsewhere in this southern region, they would be expected to seed ephemeral SEPR vents with genetically diverse propagules, reversing the observed cline in genic diversity. Alternatively, the existence of many small (and genetically homogeneous) source populations seems unlikely given the high probability of local extinctions. Although R. pachyptila might exist elsewhere along the Pacific Antarctic Ridge, the present genetic evidence is consistent with colonization of the SEPR from the NEPR or GAR axes.

Conclusion

Other studies of mitochondrial variation in vent species have reported star-like genealogies that suggest recent range expansions or demographic instability [reviewed in [15]]. Yet, the nuclear and mitochondrial gene networks from R. pachyptila are simpler and shallower than those found in all other vent annelids and mollusks examined to date [e.g., [25, 26, 6063]]. Estimates of nucleotide diversity (π) for mitochondrial COI sequences are smaller by an order of magnitude than the nearest estimates from the other vent taxa. Comparable estimates do not exist for the present nuclear markers, but the π values were of the same order of magnitude as mitochondrial Cytb (Table 5). Although our attempts to examine indicators of demographic expansion with mismatch distributions were frustrated by the limited polymorphism, the overall lack of nuclear and mitochondrial variation itself suggests that R. pachyptila may be extraordinarily prone to population turnover that leads to rapid coalescence. High rates of local extinction coupled with relatively high rates of dispersal (colonization) from a limited number of source populations (i.e. propagule pool) are expected to produce the observed pattern of reduced diversity within and among disjunct colonies of a metapopulation [6, 7]. Which demographic factors have helped to retain genetic diversity in the co-distributed vent species, however, remain mysteries.

Abbreviations

ROV:

remotely operated vehicle

HOV:

human occupied vehicle

COI cytochrome-c:

-oxidase subunit I, Cytb cytochrome b

scnDNA:

single-copy nuclear DNA

IBD:

isolation-by-distance

SEPR:

southern East Pacific Rise

NEPR:

northern East Pacific Rise

GAR:

Galápagos Rift

AFLP:

amplified fragment length polymorphism.

References

  1. Wright S: Breeding structure of populations in relation to speciation. Am Natur. 1940, 74: 232-248. 10.1086/280891.

    Article  Google Scholar 

  2. Slatkin M: Gene flow and genetic drift in a species subject to frequent local extinctions. Theor Popul Biol. 1977, 12: 253-262. 10.1016/0040-5809(77)90045-4.

    Article  CAS  PubMed  Google Scholar 

  3. Maruyama T, Kimura M: Genetic variability and effective population size when local extinction and recolonization of subpopulations are frequent. Proc Natl Acad Sci USA. 1980, 77: 6710-6714. 10.1073/pnas.77.11.6710.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Slatkin M: Gene flow and the geographic structure of natural populations. Science. 1987, 236: 787-792. 10.1126/science.3576198.

    Article  CAS  PubMed  Google Scholar 

  5. Wade MJ, McCauley DE: Extinction and recolonization: their effects on the genetic differentiation of local populations. Evolution. 1988, 42: 995-1005. 10.2307/2408915.

    Article  Google Scholar 

  6. Whitlock M, McCauley D: Some population genetic consequences of colony formation and extinction: genetic correlations within founding groups. Evolution. 1990, 44: 1717-1724. 10.2307/2409501.

    Article  Google Scholar 

  7. Pannell J, Charlesworth B: Neutral genetic diversity in a metapopulation with recurrent local extinction and recolonization. Evolution. 1999, 53: 664-676. 10.2307/2640708.

    Article  Google Scholar 

  8. Emerson BC, Paradis E, Thèbaud C: Revealing the demographic histories of species using DNA sequences. Trends Ecol Evol. 2001, 16: 707-716. 10.1016/S0169-5347(01)02305-9.

    Article  Google Scholar 

  9. Bazin E, Glemin S, Galtier N: Population size does not influence mitochondrial genetic diversity in animals. Science. 2006, 312: 570-572. 10.1126/science.1122033.

    Article  CAS  PubMed  Google Scholar 

  10. Mulligan CJ, Kitchen A, Miyamoto MM: Comment on "Population size does not influence mitochondrial genetic diversity in animals.". Science. 2006, 314: 1390a-10.1126/science.1132585.

    Article  CAS  Google Scholar 

  11. Wares JP, Barber PH, Ross-Ibarra J, Sotka EE, Toonen RJ: Mitochondrial DNA and population size. Science. 2006, 314: 1388-1389. 10.1126/science.314.5804.1388.

    Article  PubMed  Google Scholar 

  12. Vitalis R, Dawson K, Boursot P: Interpretation of variation across marker loci as evidence of selection. Genetics. 2001, 158: 1811-1823.

    CAS  PubMed  PubMed Central  Google Scholar 

  13. Galtier N, Depaulis F, Barton N: Detecting bottlenecks and selective sweeps from DNA sequence polymorphism. Genetics. 2000, 155: 981-987.

    CAS  PubMed  PubMed Central  Google Scholar 

  14. Vrijenhoek RC: Gene flow and genetic diversity in naturally fragmented metapopulations of deep-sea hydrothermal vent animals. J Hered. 1997, 88: 285-293.

    Article  CAS  PubMed  Google Scholar 

  15. Vrijenhoek RC: Genetic diversity and connectivity of deep-sea hydrothermal vent metapopulations. Mol Ecol. 2010, 19: 4391-4411. 10.1111/j.1365-294X.2010.04789.x.

    Article  PubMed  Google Scholar 

  16. Shea K, Metaxas A, Young CR, Fisher CR: Processes and interactions in macrofaunal assemblages at hydrothermal vents: a modeling perspective. Magma to Microbe: Modeling Hydrothermal Processes at Oceanic Spreading Centers. Edited by: Lowell RP, Perfit MR, Seewald J, Metaxas A. 2009, Washington, D.C.: Geophysical Monograph Series, 178: 259-274.

    Chapter  Google Scholar 

  17. Chevaldonné P, Jollivet D, Vangrieshiem A, Desbruyères D: Hydrothermal-vent alvinellid polychaete dispersal in the eastern Pacific. 1. Influence of vent site distribution, bottom currents, and biological patterns. Limnol Oceanogr. 1997, 42: 67-80.

    Article  Google Scholar 

  18. Shank TM, Fornari DJ, Von Damm KL, Lilley MD, Haymon RM, Lutz RA: Temporal and spatial patterns of biological community development at nascent deep-sea hydrothermal vents (9°50'N East Pacific Rise). Deep Sea Res II. 1998, 45: 465-515. 10.1016/S0967-0645(97)00089-1.

    Article  Google Scholar 

  19. Juniper SK, Tunnicliffe V: Crustal accretion and the hot vent ecosystem. Philos Trans Math Physic Engin Sci. 1997, 355: 459-474. 10.1098/rsta.1997.0017.

    Article  Google Scholar 

  20. Lutz RA, Shank TM, Fornari DJ, Haymon RM, Lilley MD, Von Damm KL, Desbruyères D: Rapid growth at deep-sea vents. Nature. 1994, 371: 663-664. 10.1038/371663a0.

    Article  Google Scholar 

  21. Marsh AG, Mullineaux LS, Young CM, Manahan DT: Larval dispersal potential of the tubeworm Riftia pachyptila at deep-sea hydrothermal vents. Nature. 2001, 411: 77-80. 10.1038/35075063.

    Article  CAS  PubMed  Google Scholar 

  22. Kimura M, Weiss WH: The stepping stone model of genetic structure and the decrease of genetic correlation with distance. Genetics. 1964, 49: 561-576.

    CAS  PubMed  PubMed Central  Google Scholar 

  23. Audzijonyte A, Vrijenhoek R: When gaps really are gaps: statistical phylogeography of hydrothermal vent invertebrates. Evolution. 2010, 64: 2369-2384.

    PubMed  Google Scholar 

  24. Black MB, Lutz RA, Vrijenhoek RC: Gene flow among vestimentiferan tube worm (Riftia pachyptila) populations from hydrothermal vents of the Eastern Pacific. Mar Biol. 1994, 120: 33-39.

    Google Scholar 

  25. Hurtado LA, Lutz RA, Vrijenhoek RC: Distinct patterns of genetic differentiation among annelids of eastern Pacific hydrothermal vents. Mol Ecol. 2004, 13: 2603-2615. 10.1111/j.1365-294X.2004.02287.x.

    Article  CAS  PubMed  Google Scholar 

  26. Plouviez S, Shank TM, Faure B, Daguin-thiebaut C, Viard F, Lallier FH, Jollivet D: Comparative phylogeography among hydrothermal vent species along the East Pacific Rise reveals vicariant processes and population expansion in the South. Mol Ecol. 2009, 18: 3903-3917. 10.1111/j.1365-294X.2009.04325.x.

    Article  CAS  PubMed  Google Scholar 

  27. Plouviez S, Le Guen D, Lecompte O, Lallier F, Jollivet D: Determining gene flow and the influence of selection across the equatorial barrier of the East Pacific Rise in the tube-dwelling polychaete Alvinella pompejana. BMC Evol Biol. 2010, 10: 220-10.1186/1471-2148-10-220.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Shank TM, Halanych KM: Toward a mechanistic understanding of larval dispersal: insights from genomic fingerprinting of the deep-sea hydrothermal vent tubeworm Riftia pachyptila. Mar Ecol. 2007, 28: 25-35. 10.1111/j.1439-0485.2007.00146.x.

    Article  CAS  Google Scholar 

  29. Fusaro AJ, Baco AR, Gerlach G, Shank TM: Development and characterization of 12 microsatellite markers from the deep-sea hydrothermal vent siboglinid Riftia pachyptila. Mol Ecol Resources. 2008, 8: 132-134. 10.1111/j.1471-8286.2007.01897.x.

    Article  CAS  Google Scholar 

  30. Boore J, Brown W: Complete sequence of the mitochondrial DNA of the annelid worm Lumbricus terrestris. Genetics. 1995, 141: 305-319.

    CAS  PubMed  PubMed Central  Google Scholar 

  31. Jennings RM, Halanych KM: Mitochondrial genomes of Clymenella torquata (Maldanidae) and Riftia pachyptila (Siboglinidae): Evidence for conserved gene order in Annelida. Mol Biol Evol. 2005, 22: 210-222. 10.1093/molbev/msi008.

    Article  CAS  PubMed  Google Scholar 

  32. Karl SA, Avise JC: PCR-based assays of Mendelian polymorphisms from anonymous single-copy nuclear DNA: techniques and applications for population genetics. Mol Biol Evol. 1993, 10: 342-374.

    CAS  PubMed  Google Scholar 

  33. Jarman SN, Ward RD, Elliott NG: Oligonucleotide primers for PCR amplification of coelomate introns. Mar Biotech. 2002, 4: 347-355. 10.1007/s10126-002-0029-6.

    Article  CAS  Google Scholar 

  34. Stephens M, Smith NJ, Donnelly P: A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001, 68: 978-989. 10.1086/319501.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  35. Stephens M, Donnelly P: A comparison of Bayesian methods for haplotype reconstruction from population genotype data. Am J Hum Genet. 2003, 73: 1162-1169. 10.1086/379378.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  36. Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 4: 331-346.

    Google Scholar 

  37. Excoffier , Lischer H: Arlequin ver 3.5: an integrated software package for population genetics data analysis. Swiss Institute of Bioinformatics. 2009, 1-174.

    Google Scholar 

  38. Watterson GA: On the number of segregating sites in genetical models without recombination. Theor Popul Biol. 1975, 7: 256-276. 10.1016/0040-5809(75)90020-9.

    Article  CAS  PubMed  Google Scholar 

  39. Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147: 915-925.

    CAS  PubMed  PubMed Central  Google Scholar 

  40. Nei M: Molecular Evolutionary Genetics. 1987, New York, NY, USA: Columbia University Press

    Google Scholar 

  41. Kalinowsky S: Hp-Rare 1.0: a computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes. 2005, 5: 187-189. 10.1111/j.1471-8286.2004.00845.x.

    Article  Google Scholar 

  42. Raymond M, Rousset F: Genepop (Ver. 1.2): population genetics software for exact tests and ecumenicism. J Hered. 1995, 86: 248-249.

    Google Scholar 

  43. Pritchard JK, Stephens M, Donnelly P: Inference of population structure using multilocus genotype data. Genetics. 2000, 155: 945-959.

    CAS  PubMed  PubMed Central  Google Scholar 

  44. Casgrain P, Legendre P: The R Package for Multivariate and Spatial Analysis. 2001, Département de sciences biologiques, Université de Montréal, v. 4.0 d6 - Users Manual., [http://www.bio.umontreal.ca/casgrain/en/labo/R/index.html]

    Google Scholar 

  45. Smouse PE, Long JC, Sokal RR: Multiple regression and correlation extensions of the Mantel test of matrix correspondence. Syst Zool. 1986, 35: 627-632. 10.2307/2413122.

    Article  Google Scholar 

  46. Beaumont MA, Nichols RA: Evaluating loci for use in the genetic analysis of population structure. Proc R Soc Lond B. 1996, 263: 1619-1626. 10.1098/rspb.1996.0237.

    Article  Google Scholar 

  47. Antao T, Lopes A, Lopes R, Beja-Pereira A: Lositan: A workbench to detect molecular adaptation based on a F st -outlier method. BMC Bioinform. 2008, 9: 323-10.1186/1471-2105-9-323.

    Article  Google Scholar 

  48. DeMets C, Gordon RG, Argus DF, Stein S: Current plate motions. Geophys J Int. 1990, 101: 425-478. 10.1111/j.1365-246X.1990.tb06579.x.

    Article  Google Scholar 

  49. DeMets C, Gordon RG, Argus DF, Stein S: Effect of recent revisions to the geomagnetic reversal time scale on estimates of current plate motions. Geophys Res Let. 1994, 21: 2191-2194. 10.1029/94GL02118.

    Article  Google Scholar 

  50. Soudarin L, Cazenave A: Large-scale tectonic plate motions measured with the DORIS space geodesy system. Geophys Res Let. 1995, 22: 469-472. 10.1029/94GL03382.

    Article  Google Scholar 

  51. Birky CWJ, Maruyama T, Fuerst P: An approach to population and evolutionary genetic theory for genes in mitochondria and chloroplasts, and some results. Genetics. 1983, 103: 513-527.

    PubMed  PubMed Central  Google Scholar 

  52. Crochet P: Genetic structure of avian populations--allozymes revisited. Mol Ecol. 2000, 9: 1463-1469. 10.1046/j.1365-294x.2000.01026.x.

    Article  CAS  PubMed  Google Scholar 

  53. Reid JL: On the total geostrophic circulation of the South Pacific Ocean: Flow patterns, tracers and transports. Progr Oceanogr. 1986, 16: 1-61. 10.1016/0079-6611(86)90036-4.

    Article  Google Scholar 

  54. Slatkin M: Isolation by distance in equilibrium and non-equilibrium populations. Evolution. 1993, 47: 264-279. 10.2307/2410134.

    Article  Google Scholar 

  55. Jollivet D, Chevaldonné P, Planque B: Hydrothermal-vent alvinellid polychaete dispersal in the Eastern Pacific. 2. A metapopulation model based on habitat shifts. Evolution. 1999, 53: 1143-1156. 10.2307/2640817.

    Article  Google Scholar 

  56. MacDonald KC: Linkages between faulting, volcanism, hydrothermal activity and segmentation on fast spreading centers. Faulting and Magmatism at Mid-Ocean Ridges. Edited by: Buck W, Delaney P, Karson J, Lagabrielle Y. 1998, Washington, DC: American Geophysical Union, 27-58.

    Chapter  Google Scholar 

  57. Hey R, Massoth G, Vrijenhoek R, Rona P, Lupton J, Butterfield D: Hydrothermal vent geology and biology at Earth's fastest spreading rates. Mar Geophys Res. 2006, 27: 137-153. 10.1007/s11001-005-1887-x.

    Article  Google Scholar 

  58. Van Dover CL: Community structure of mussel beds at deep-sea hydrothermal vents. Mar Ecol Prog Ser. 2002, 230: 137-158. 10.3354/meps230137.

    Article  Google Scholar 

  59. Jollivet D, Lallier FH, Barnay AS, Bienvenue N, Bonnivard N, Briand P, Cambon-Bonavita MA, Comtet T, Cosson R, Daguin C, et al: The BIOSPEEDO cruise: a new survey of hydrothermal bens along the South East Pacific Rise from 7°24'S to 21°33'S. InterRidge News. 2004, 13: 20-26.

    Google Scholar 

  60. Teixeira S, Cambon-Bonavita MA, Serrão E, Desbruyères D, Arnaud-Haond S: Recent population expansion and connectivity in the hydrothermal shrimp Rimicaris exoculata along the Mid-Atlantic Ridge. J Biogeogr. 2010, 38: 264-574.

    Google Scholar 

  61. Kyuno A, Shintaku M, Fujita Y, Matsumoto H, Utsumi M, Watanabe H, Fujiwara Y, Miyazaki JI: Dispersal and differentiation of deep-sea mussels of the genus Bathymodiolus (Mytilidae, Bathymodiolinae). J Mar Biol. 2009, Article ID625672:15 pages

    Google Scholar 

  62. Young CR, Fujio S, Vrijenhoek RC: Directional dispersal between mid-ocean ridges: deep-ocean circulation and gene flow in Ridgeia piscesae. Mol Ecol. 2008, 17: 1718-1731. 10.1111/j.1365-294X.2008.03609.x.

    Article  CAS  PubMed  Google Scholar 

  63. Johnson SB, Young CR, Jones WJ, Warén A, Vrijenhoek RC: Migration, Isolation, and Speciation of Hydrothermal Vent Limpets (Gastropoda; Lepetodrilidae) Across the Blanco Transform Fault. Biol Bull. 2006, 210: 140-157. 10.2307/4134603.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

We thank Dr. Yong Jin Won and two anonymous reviewers for constructive criticisms that improved the manuscript. The project was conducted with the expert help and enthusiasm of captains, crews and pilots assigned to the R/V Atlantis and DSV Alvin (Woods Hole Oceanographic Institute) and to the R/V Western Flyer and ROV Tiburon (Monterey Bay Aquarium Research Institute, MBARI). Funding was provided by the David and Lucile Packard Foundation (to MBARI) and the following NSF grants: OCE-9910799 and OCE-0241613 (to RCV); and OCE-0937371, OCE-9529819, OCE-0327353 and ESI-0087679 (to RAL).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Robert C Vrijenhoek.

Additional information

Authors' contributions

RAL and RCV conceived of the project and led oceanographic expeditions that collected study materials. SAK developed the scnDNA primers. DKC optimized DNA primers, and conducted the DNA sequencing. RCV, DKC, and SBJ completed the statistical analyses and prepared a final draft of the manuscript.

Electronic supplementary material

12862_2010_1713_MOESM1_ESM.DOC

Additional file 1: Non-degenerate and internal primer pairs (indicated with *) developed for population screening of Riftia pachyptila. GenBank accession numbers are indicated along with the PCR program (see text for details). A table of DNA primer pairs used in screening nuclear and mitochondrial genes in R. pachyptila. (DOC 116 KB)

12862_2010_1713_MOESM2_ESM.DOCX

Additional file 2: Allelic diversity in Riftia pachyptila. A table of allelic frequencies for nuclear and mitochondrial genes in R. pachyptila. DNA haplotypes for alleles are identified by their polymorphic sites. (DOCX 15 KB)

Authors’ original submitted files for images

Below are the links to the authors’ original submitted files for images.

Authors’ original file for figure 1

Authors’ original file for figure 2

Authors’ original file for figure 3

Rights and permissions

Open Access This article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License ( https://creativecommons.org/licenses/by/2.0 ), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and permissions

About this article

Cite this article

Coykendall, D.K., Johnson, S.B., Karl, S.A. et al. Genetic diversity and demographic instability in Riftia pachyptilatubeworms from eastern Pacific hydrothermal vents. BMC Evol Biol 11, 96 (2011). https://doi.org/10.1186/1471-2148-11-96

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/1471-2148-11-96

Keywords