|Home | About | Journals | Submit | Contact Us | Français|
Reef fishes disperse primarily as oceanic “pelagic” larvae, and debate continues over the extent of this dispersal, with recent evidence for geographically restricted (closed) populations in some species. In contrast, moray eels have the longest pelagic larval stages among reef fishes, possibly providing opportunities to disperse over great distances. We test this prediction by measuring mitochondrial DNA (mtDNA) and nuclear DNA variation in 2 species of moray eels, Gymnothorax undulatus (N = 165) and G. flavimarginatus (N = 124), sampled at 14–15 locations across the Indo-Pacific. The mtDNA data comprise 632 bp of cytochrome b and 596 bp of cytochrome oxidase I. Nuclear markers include 2 recombination-activating loci (421 bp of RAG-1 and 754 bp of RAG-2). Analyses of molecular variance and Mantel tests indicate little or no genetic differentiation, and no isolation by distance, across 22000 km of the Indo-Pacific. We estimate that mitochondrial genetic variation coalesces within the past about 2.3 million years (My) for G. flavimarginatus and within the past about 5.9 My for G. undulatus. Permutation tests of geographic distance on the mitochondrial haplotype networks indicate recent range expansions for some younger haplotypes (estimated within ~600000 years) and episodic fragmentation of populations at times of low sea level. Our results support the predictions that the extended larval durations of moray eels enable ocean-wide genetic continuity of populations. This is the first phylogeographic survey of the moray eels, and morays are the first reef fishes known to be genetically homogeneous across the entire Indo-Pacific.
Coral reefs occur sporadically in a matrix of open ocean that spans 70% of the planet. Reef inhabitants that disperse only as adults can show population genetic structuring on scales of tens of kilometers (Bernardi 2000). However, a pelagic (oceanic) larval stage permits most reef fishes to disperse across open ocean. Larvae either swim actively or drift passively in surface currents, sometimes modifying their buoyancy to use both surface and deep currents (Leis 1991; Montgomery et al. 2001; Atema et al. 2002; Kingsford et al. 2002). This ability to disperse across unsuitable habitat prior to settling on a suitable reef enables dispersal among regions separated by thousands of kilometers.
The geographic distance over which larvae disperse varies across taxa and geographic regions and is probably influenced by the amount of time spent in the pelagic stage (the pelagic larval duration [PLD]; Ekman 1953; Scheltema 1968; Crisp 1978). The simplest expectation is that long PLDs yield proportionately greater gene flow among populations (Jones et al. 2005; Taylor and Hellberg 2005). For many reef fishes, this prediction is upheld by population genetic patterns (Craig et al. 2007). However, PLD is not necessarily predictive (Weersing and Toonen 2009), and many reef species maintain greater population genetic structure than expected from their PLD (Planes and Fauvelot 2002; Swearer et al. 2002; Bernardi et al. 2003; Taylor and Hellberg 2003; Thacker et al. 2007). These findings have prompted research into the types of geographic, oceanographic, and biological factors that mediate larval dispersal, including variation in effective population size, habitat preference, and demographic changes in response to major climatic and environmental shifts (Hellberg 2007; Leis 2007; Leis et al. 2007; Rocha et al. 2007; Floeter et al. 2008).
Genetic fragmentation among populations of reef fish species often coincides with major oceanic and geographic barriers (Briggs 1961; Springer 1982; Rocha et al. 2007). The 2 major phylogeographic breaks of the Indo-Pacific occur across the Eastern Pacific Barrier (EPB; 5000–7000 km of open ocean that separate eastern and central Pacific reefs) and the Sunda Shelf (SS; an area of shallow reefs, located just west of Indonesia, that are exposed during low sea levels and separate the Indian and Pacific oceans; Figure 1).
Some reef fishes show high genetic connectivity across thousands of kilometers (Muss et al. 2001; Rocha et al. 2002; Bay et al. 2004), indicating gene flow across the SS (Horne et al. 2008) or the EPB (Lessios and Robertson 2006). However, no reef fish with a nonpelagic adult stage has been shown to maintain connectivity across both the SS and the EPB. Fish species showing geographically widespread populations typically have PLDs exceeding 45 days, whereas most reef fishes have PLDs between 15 and 30 days (Lester and Ruttenberg 2005). Surgeonfishes of the genus Naso are among the most extreme examples of reef fishes maintaining geographically widespread gene flow. Naso vlamingii (Klanten et al. 2007) and 2 congeners (Horne et al. 2008) demonstrate similarly high levels of mitochondrial haplotype diversity and lack of phylogeographic structure across their Indo-Pacific ranges west of the EPB. These 3 species share extended PLDs and similar life histories.
Moray eels (Muraenidae) are solitary predators in reef ecosystems. They comprise about 200 species globally and approximately 150 species within the Indo-Pacific. Many moray species span the entire Indo-Pacific, occupying coral reefs and rocky ledges from the intertidal zone to 200+ m deep. Moray eels are poor swimmers as juveniles and adults and maintain high site fidelity to a few square meters of reef (Böhlke et al. 1989). In a pattern typical of reef fishes, dispersal occurs by pelagic larvae. However, the morays and elopomorph fishes (which includes all true eels, tarpons, tenpounders, and ladyfish) are unique in having a slender elongate larval form called a “leptocephalus” (Figure 1). These larvae are among the simplest, long-lived, and self-sustaining vertebrate forms. They are transparent except for eye pigmentation, and the body wall may be only a few cells thick. They have no identifiable digestive activity or nutrient stores; yet, leptocephalus larvae can persist up to 2 years in pelagic environments sustained by consuming dissolved organic carbon and the fecal pellets and waste products of zooplankton and other larvae (Mochioka and Iwamizu 1996; Bishop and Torres 1999; Ishikawa et al. 2001).
Moray eels are difficult to sample because of their solitary, aggressive, reclusive habits, and this study represents the first phylogeographic survey of the group. Previous authors have suggested that they have extensive oceanic dispersal, based on the extended PLD of the leptocephalus larvae and “the broad distributions of many species within the vast Indo-Pacific region” (Randall 2007). Here, we test this hypothesis for the first time with phylogeographic surveys and inferred coalescence times for haplotype variation within Gymnothorax flavimarginatus and G. undulatus. Both species are ecological generalists on reefs from 0 to more than 100 m deep, with distributions that traverse the SS and the EPB. We sample G. undulatus and G. flavimarginatus from reefs across the Indo-Pacific, from South Africa to Panama, using portions of the mitochondrial genes encoding cytochrome b (CYB) and cytochrome oxidase subunit I (COI), plus portions of 2 nuclear recombination-activating loci, RAG-1 and RAG-2. Our results reveal the first cases of high gene flow through larval dispersal across both the SS and the EPB.
Specimens were captured in fish traps and lobster traps, and by pole spears while snorkeling or scuba diving. Here, we report the collection of 289 specimens with only 2 injuries, a minor bite wound to the hand of author B.W.B. from G. undulatus and a bite to the foot of documentary filmmaker Ziggy Livnat. This is a notable accomplishment as moray eels are widely (and rightfully) feared for sudden and egregious acts of aggression (Riordan et al. 2004). In particular, G. undulatus is “especially vicious” (Hiatt and Strasburg 1960), and several bite victims are known to the authors (see also Randall 2007).
Fin clips or muscle samples were taken from live specimens and other tissue samples from tissue storage banks. Sampling includes 124 G. flavimarginatus specimens from 15 sites and 165 specimens of G. undulatus from 14 sites spanning the Indo-Pacific (Table 1, Figure 1). All samples were collected between 2001 and 2009. DNA was isolated using Viogene genomic DNA extraction column kits. Polymerase chain reactions (PCRs) were performed in 25-μl reactions of 5 μl of Promega (www.promega.com) 5× buffer, 2.5 μl of 25 mM MgCl2, 2.5 μl of 0.2 μM dNTPs, 2.5 μl of 0.2 μM of each primer, 0.125 μl (1 unit) of Promega GoTaq DNA polymerase, and 2 μl of template DNA at approximately 5 ng/μl. A 632-bp fragment of the gene encoding CYB was amplified using the primers L14725 (5′-GTG ACT TGA AAA ACC ACC GTT G-3′; Song et al. 1998) and H15573 (5′-AAT AGG AAG TAT CAT TCG GGT TTG ATG-3′; Taberlet et al. 1992) with an annealing temperature of 50°C. A 596-bp fragment of the gene encoding COI was amplified using primers FishF2 (5′-TCG ACT AAT CAT AAA GAT ATC GGC AC-3′) and FishR2 (5′-ACT TCA GGG TGA CCG AAG AAT CAG AA-3′; Ward et al. 2005) with an annealing temperature of 50°C. A 421-bp fragment of the nuclear recombination-activating gene RAG-1 was amplified using primers designed by author J.S.R. based on published marine fish RAG-1 sequences: RAG1-F3 (5′-GCC TCA GAA AAC ATG GTG CT-3′) and RAG1-R3 (5′-CCA CAC AGG TTT CAT CTG GA-3′) with an annealing temperature of 50°C. A 754-bp fragment of the nuclear recombination-activating gene RAG-2 was sequenced using primers RAG2-F3 (5′-AGG TGA CCC TTC GTT GTC AG-3′) and RAG2-R3 (5′-ATG AGG CTC CCT TCC AAA GT-3′), also designed by author J.S.R., with an annealing temperature of 52°C. The thermal profiles for PCR were as follows: 95°C for 3 min, followed by 35 cycles of 95°C for 1 min, annealing temperature for 40 s, and 72°C for 45 s, with a final elongation at 72°C for 7 min.
PCR products were visualized through agarose-gel electrophoresis and purified using Exo-Sap or Viogene gel purification kits (www.viogene.com). Sequences were generated on ABI 3130 and ABI 3330 automated DNA sequencers at the Washington University Genome Sequencing Center and the Smithsonian Museum Support Center. All loci were sequenced in both the forward and the reverse directions, and all data were generated and edited by the same researcher (J.S.R.). DNA sequences were manually edited using SEQUENCHER 4.8 (Ann Arbor, MI) and aligned by hand. For nuclear markers, heterozygous positions were identified as cases wherein a secondary peak in the electropherograms reached at least 25% of the intensity of the primary peak. Gametic phases of nuclear sequences with more than a single heterozygous site were estimated using a Bayesian approach implemented in the software program PHASE 2.1 (Stephens et al. 2001; Stephens and Donnelly 2003). All PHASE analyses were run through 5 iterations, each initiated using a different random-number seed and run for 1000 iterations with a single thinning interval and 100 burn-in iterations. Consistency of results was determined by examining allele frequencies and coalescent goodness-of-fit measures estimated for each of the 4 runs. In the case that the phase of a polymorphism could not be confidently estimated with 90% posterior probability, that site was coded as missing data. Less than 3% of all nucleotide characters were coded as missing data under this criterion. None of the markers for either species appear to be under selection according to Tajima's D (Tajima 1989) tests for selection (all P > 0.1) in ARLEQUIN 3.1 (Excoffier and Schneider 2005). Last, tests of linkage disequilibrium in ARLEQUIN between the nuclear loci RAG-1 and RAG-2 were nonsignificant (P > 0.1) for both species.
We constructed haplotype networks for the mitochondrial fragments from CYB and COI concatenated as a single genetic marker and separately for each of the nuclear markers RAG-1 and RAG-2. In all cases, haplotype networks were constructed with TCS 2.1 (Clement et al. 2000) under 95% statistical parsimony criteria. This approach permits a qualitative assessment of the geographic distributions of alleles. The standard molecular diversity measures of haplotype diversity (h) and nucleotide diversity (π) were calculated in ARLEQUIN accordingly to Nei (1987).
We tested for population structure and phylogeographic patterns with 3 independent analyses: 1) an analysis of molecular variance (AMOVA; Excoffier et al. 1992), 2) pairwise Φst and a Mantel test of geographic versus genetic distances (Φst), and 3) permutation tests of geographic distances on haplotype networks. Sample sizes less than N = 5 were excluded from any pairwise comparisons. After preliminary trials, the 2 mitochondrial DNA (mtDNA) fragments were analyzed as a single composite sequence. We subjected each marker (mtDNA, RAG-1, RAG-2) to AMOVA using ARLEQUIN based on haplotype frequencies and molecular distances among haplotypes. Individuals were grouped by localities and then by 3 larger subgroups comprising localities from the Indian Ocean, the western and central Pacific Ocean, and the eastern Pacific Ocean. These partitions correspond to the 3 oceanic regions demarcated by the SS and the EPB (Figure 1). Statistical significance of the corresponding ARLEQUIN statistic, Φsc, was determined by 10000 permutations.
We conducted Mantel tests in the program PASSAGE 2.0 (Rosenberg 2001) to measure correlations between geographic and genetic distances among populations. Genetic distances for all loci were computed in ARLEQUIN as pairwise Φst values based on the number of nucleotide differences between alleles and allele frequency differences between populations. We report only uncorrected differences because corrected pairwise differences under a variety of mutation models yielded congruent results. Significance of pairwise Φst values was computed using the B-H (Benjamini and Hochberg 1995; Misawa et al. 2008) and the more conservative Bonferroni (1936) corrections for multiple comparisons. The pairwise matrix of geographic distances between populations was obtained with ARCGIS 9.3 as the shortest straight-line distance that did not cross dry land. We could not incorporate oceanic currents in this analysis because most available data are for surface currents, whereas moray eel larvae use both surface currents and counterflowing deeper currents (Böhlke et al. 1989). Geographic and genetic distance (Φst) matrices were tested for correlation separately for each molecular maker (mtDNA, RAG-1, and RAG-2).
Associations of geographic and genetic distances were resolved using a permuted chi-square test of genetic distances among hierarchically nested groups of haplotypes and geographic areas covered by distributions of grouped haplotypes (Templeton 1998). This approach is the initial step in nested clade phylogeographic analysis (NCPA; Templeton 1998) and mirrors the Mantel test, but instead of using single populations as the sole level of sampling, it tests for associations between geographic distance and hierarchically nested groups from the haplotype network. Geographic distances were identical to those of the Mantel tests. Tests were performed separately on each locus (mtDNA, RAG-1, and RAG-2) in the program GEODIS 2.6 (which includes the Dunn–Sidak correction for multiple tests; Posada et al. 2000) using 10000 permutations. Inferences were made using the 2008 inference key available at http://darwin.uvigo.es/software/geodis.html. Loops denoting ambiguous groupings in haplotype networks were treated as follows: all possible nested combinations of resolutions to the loops were computed, and analyses were repeated for each possible network configuration. Alternative resolutions of loops did not alter inferences in this study.
There is little information available on the time since divergence from a most recent common ancestor of moray eels and their closest relatives, the genera Conger and Anguilla (Inoue et al. 2003). Without a strong fossil record or obvious cases of vicariance, one cannot derive a species-specific mutation rate for estimating coalescence times. For the mitochondrial gene CYB, we used a divergence rate of 2% per million year (My), which has been widely applied to phylogeographic studies of reef fishes (Bowen et al. 2001, 2006b) and is supported by measurements for goby (Gnatholepis) species separated by the Isthmus of Panama (1.95–2.17% divergence per My; Rocha et al. 2005). We used a mutation rate for the COI data of 1.2% sequence divergence per My as calibrated by Bermingham et al. (1997) for marine fishes. We estimated the time to coalescence for standing genetic variation in each of the 2 study species and for haplotype groups with significant inferences in NCPA, using the Bayesian Markov chain Monte Carlo (MCMC) approach as implemented in BEAST 1.5.2 (Drummond and Rambaut 2007). We used a congeneric outgroup and determined the target species (G. undulatus or G. flavimarginatus) to be monophyletic for haplotypic variation. Models of evolution were computed in MrModelTest (Posada and Buckley 2004), and the combined mtDNA data were partitioned by gene region by codon position, for a total of 6 partitions. A fixed mutation rate was estimated at 1.6% divergence per My (the weighted average of 2% and 1.2% for CYB and COI, respectively) and input as 0.008 substitutions per site per lineage per My. A coalescent tree prior of expansion growth was used, as recommended for estimating intraspecific coalescence times when population growth is inferred (see results of mismatch distributions; Drummond and Rambaut 2007). Simulations ran for 30 million generations, with sampling every 3000 generations. Two independent runs were computed for each species, and runs were then combined with a 10% burn-in using LogCombiner 1.5.2 (Drummond and Rambaut 2007). Adequate parameter estimation and convergence of chains in the MCMC run were computed in the program Tracer 1.4 (Rambaut and Drummond 2007). The effective Sample Size values for each parameter were more than 200, and 95% highest posterior density (HPD) intervals are reported as encompassing the uncertainty around coalescence time estimates within each species. This estimate of confidence intervals is directly proportional to and does not account for variation in the mutation rates used.
Frequency distributions of numbers of mutational differences between haplotypes sampled within a species (=mismatch distribution) were computed in ARLEQUIN using the least-square method (Schneider et al. 2000). To compute the confidence intervals around these values and the ability to reject a model of demographic expansion through time, we performed 100000 parametric bootstraps and calculated Harpending's raggedness index (HRI; Harpending 1994). We also computed Fu's (1997) FS statistics to identify cases of demographic expansion, which are indicated by significant negative values.
Nuclear and mitochondrial haplotype networks for G. flavimarginatus and G. undulatus are presented in Figure 2. Notably, haplotypes from the major geographic subdivisions of the Indo-Pacific are dispersed throughout the networks. Combining CYB and COI, we identify 104 mtDNA haplotypes for G. flavimarginatus (N = 124; GenBank accession number GU175445-GU175541) and 112 haplotypes for G. undulatus (N = 165; GU175576-GU175657). No mitochondrial haplotype occurs in samples of all 3 major oceanic regions (Indian Ocean, western/central Pacific Ocean, and eastern Pacific Ocean), but 4 haplotypes in G. flavimarginatus and 6 in G. undulatus occur in adjacent regions; all other haplotypes are observed in a single region. The nuclear markers RAG-1 and RAG-2 reveal less variation than the combined mitochondrial markers, as expected from their slower mutation rate and larger effective population size. RAG-1 sequences resolved 16 haplotypes for G. flavimarginatus (GU175542-GU175557) and 23 haplotypes for G. undulatus (GU175658-GU175680). RAG-2 sequences included 18 haplotypes for G. flavimarginatus (GU175558-GU175575) and 22 haplotypes for G. undulatus (GU175681-GU175702). The most common nuclear haplotypes occur in most localities, with haplotype sharing among the 3 major oceanic regions (Figure 2). Groupings of haplotypes endemic to 1 of the 3 major oceanic regions occur only in the mitochondrial haplotype networks. Molecular diversity indices for G. flavimarginatus at mtDNA were h = 0.997 and π = 0.012; at RAG-1, h = 0.706 and π = 0.0050; and at RAG-2, h = 0.912 and π = 0.0043. Molecular diversity indices for G. undulatus at mtDNA were h = 0.972 and π = 0.012; at RAG-1, h = 0.808 and π = 0.0053; and at RAG-2, h = 0.907 and π = 0.0067.
Results of a locus-by-locus AMOVA revealed no population structure among the 3 geographic groupings demarcated by the EPB and the SS. For G. flavimarginatus, the mtDNA yields a nonsignificant Φsc = 0.0088 (P = 0.29) and the nuclear markers a nonsignificant Φsc < 0.00001 (P = 0.89) for RAG-1 and Φsc < 0.00001 (P = 0.73) for RAG-2. The mtDNA results for G. undulatus reveal a low but significant Φsc = 0.039 (P = 0.007). The nuclear markers yield nonsignificant Φsc < 0.00001 (P = 0.30) for RAG-1 and Φsc < 0.00001 (P = 0.99) for RAG-2. Results for both species identify a vast majority of the variation (95–100%) occurring within populations.
The Mantel tests for correlations between geographic distances and Φst values are not significant for mtDNA, RAG-1, and RAG-2 (all R < 0.0001, P values >0.2). Tables 2 and and33 present mtDNA pairwise Φst values for G. flavimarginatus and G. undulatus, respectively. Structure in G. flavimarginatus was especially low, with 0 of 55 significant comparisons for mtDNA (Table 2) and no significant comparisons for RAG-1 or RAG-2. Similarly, G. undulatus showed low levels of structure, with 4 of 36 significant comparisons for mtDNA (Table 3) and none for RAG-1 or RAG-2.
The NCPA of the combined mtDNA sequences in G. flavimarginatus revealed only 3 nested groups that refute the hypothesis of panmixia. None of the resolutions of loops altered the chain of inference or the significance of pairwise comparisons; Figure 2 shows resolved loops only for clades with significant inferences. Group A (Figure 2, Table 4) supports an inference of contiguous range expansion, estimated to have occurred within the past 0.6 My (95% HPD 0.3–0.9 My). Tip haplotype groups in group A show expansion from the internal and mostly Hawaiian group to tip groups with haplotypes in the Seychelles, South Africa, and the eastern Pacific. Groups B and C yield inferences of restricted gene flow with isolation by distance coalescing within the past 0.5 My (95% HPD 0.2–0.8 My) and 0.4 My (95% HPD 0.3–0.7 My), respectively. Within group B, tip subgroup B-2 is observed only in Hawaii, whereas the more interior subgroup B-1 from which it derives includes haplotypes distributed across the Indo-Pacific. Both subgroups of group C include haplotypes from Hawaii, but they differ in having haplotypes from the Indian Ocean (subgroup C-2) versus the eastern Pacific Ocean (subgroup C-1).
The NCPA for G. undulatus includes 2 nesting levels (groups D and E) that refute the null model of panmixia, both involving range expansions of haplotypes (Figure 2, Table 4) that coalesce within 0.6 My (95% HPD 0.4–0.8 My) and 0.3 My (95% HPD 0.2–0.5 My), respectively. Range expansion in group D includes tip haplotypes sampled from Australia (subgroup D-1) or Fiji (subgroup D-2) from an inferred ancestral range in Hawaii, indicating some historical fragmentation followed by range expansion within the western Pacific. Group E shows a tip haplotype observed in Hawaii and the Indian Ocean (subgroup E-1) derived from an ancestral distribution in the Pacific Ocean (subgroup E-3).
The coalescence times for standing genetic variation based on mtDNA were 2.3 My (95% HPD 1.7–3.0 My) for G. flavimarginatus and 5.9 My (95% HPD 3.5–8.6 My) for G. undulatus. Estimates of coalescence times for G. flavimarginatus and G. undulatus are younger than those for other reef fishes with similar population structure (i.e., genus Naso; Klanten et al. 2007; Horne et al. 2008), but the confidence intervals surrounding those estimates broadly overlap.
We evaluated evidence for population expansion using mismatch distributions and Fu's FS. Mismatch distributions test the null hypothesis of demographic expansion, whereas Fu's FS tests the null hypothesis of constant population size, with significantly negative values indicating demographic expansion. Therefore, a nonsignificant result for the mismatch distribution test and a significantly negative Fu's FS both indicate population expansion. The mismatch distributions for mtDNA are shown in Figure 3. Overall, mismatch analyses indicate a pattern of population growth in both species for all markers. In G. flavimarginatus, mismatch analyses fail to reject the null hypothesis of demographic expansion for mtDNA (P = 0.78, HRI = 0.10; P = 0.72), RAG-1 (P = 0.10, HRI = 0.04; P = 0.55), and RAG-2 (P = 0.22, HRI = 0.28; P = 0.52). Similar results were obtained in mismatch analyses of G. undulatus for mtDNA (P = 0.42, HRI = 0.017; P = 0.51), RAG-1 (P = 0.5, HRI = 0.03; P = 0.68), and RAG-2 (P = 0.99, HRI = 0.004; P = 0.99). The results of the Fu's FS analyses support these conclusions, with significant negative values in G. flavimarginatus for mtDNA (−24.0; P < 0.001) and RAG-1 (−14.8; P < 0.001) but a nonsignificantly negative value for RAG-2 (−3.2; P = 0.15). Fu's FS analyses in G. undulatus indicated demographic expansion with significant negative values for mtDNA (−24.0; P < 0.001) and RAG-1 (−7.2; P = 0.019) but with a nonsignificant value for RAG-2 (0.57; P = 0.65). In summary, all analyses support conclusions of demographic expansion except Fu's FS for RAG-2 in both species. Peak pairwise differences of approximately 10 substitutions (G. undulatus) versus 17 substitutions (G. flavimarginatus) between paired haplotypes for the 1228 bp of mtDNA are consistent with transitory fragmentation of populations between 300000 and 600000 years before present (the estimated coalescent times for clades with significant inferences in NCPA). These events yielded divergent haplotype lineages that have persisted within each species through subsequent climatic cycles of approximately 100000 years (Paillard 2006).
Several reef fishes are known to have highly dispersive pelagic stages, including soldierfishes (genus Myripristis; Bowen et al. 2006b; Craig et al. 2007), pygmy angelfishes (genus Centropyge; Bowen et al. 2006a; Schultz et al. 2007), and unicornfishes (genus Naso; Klanten et al. 2007; Horne et al. 2008); however, even these species show genetic partitions (or distribution limits) at major biogeographic barriers. In the first moray species subject to phylogeographic analyses, we reveal through 3 independent tests a prevailing pattern of very little to no genetic structure among populations separated by the major biogeographic barriers of the eastern Pacific and SS. Leptocephalus larvae in general are extremely long lived in the pelagic environment and can delay metamorphosis until appropriate conditions are available (Schmidt 1923; Castonguay 1987; Tseng 1990; Crabtree et al. 1992), and it is likely that these capabilities have enabled the high levels of connectivity that we observe throughout the Indo-Pacific.
Haplotypes from the Indian, western/central Pacific, and eastern Pacific oceans are interspersed on the networks for all 3 loci (mtDNA, RAG-1, and RAG-2), with geographic clustering of haplotypes occurring only near some tips of the mtDNA network (Figure 2). Nuclear haplotypes are shared among localities separated by more than 22000 km and across 2 major biogeographic barriers. This high genetic continuity among populations with low adult vagility strongly implicates the leptocephalus larva (especially the long pelagic stage) as the dispersal engine.
The AMOVA results indicate a lack of persistent biogeographic barriers to gene flow. Only the mitochondrial locus of G. undulatus reveals a significant Φsc among regions delineated by the SS and EPB, with biogeographic barriers explaining only 3.9% of the overall variation. Mantel tests of genetic and geographic distances likewise fail to reject panmixia among geographic populations. Although the SS and EPB represent species boundaries in many taxa (Briggs 1961), and substantial obstacles to gene flow within species that cross them (Randall 1998), they appear porous to gene flow among populations of moray eels. Notably, the pairwise Φst comparisons show no consistent breaks across the SS, but G. undulatus has low/significant structure across the EPB in 2 comparisons with Hawaiian Islands. We conclude that the oceanic EPB can be a substantial barrier under contemporary conditions, whereas the SS seems to be important in a historical context but not a barrier under contemporary conditions. Leptocephalus dispersal capabilities quickly erase genetic separations developed across the SS during glacial maxima, whereas the signature patterns of partial isolation across the EPB are ongoing and detectable at low levels. Although sampling error should increase the variance on estimates of genetic divergence relative to actual divergence, our finding of widespread genetic homogeneity among populations is robust to the relatively small sample sizes within populations of the Indian and eastern Pacific oceans (Table 1).
Our NCPA results for the mitochondrial DNA indicate a pattern of transitory fragmentation among populations followed by genetic merging within 0.3–0.6 My. The NCPA is sensitive to historical events such as range expansion and previous fragmentation, and the mitochondrial locus is the only one variable enough to reveal such patterns. Three groupings of mitochondrial haplotypes from G. flavimarginatus, each group estimated to be less than 0.6 My old, show evidence of geographic fragmentation. Within group A (Figure 2, Table 4), mitochondrial haplotypes of G. flavimarginatus show evidence of a range expansion across the SS within the past approximately 0.6 My. Haplotypes in groups B and C also show some evidence of fragmentation on the timescale of 0.5–0.3 My, respectively. Within G. undulatus, mtDNA haplotype groups D and E show evidence of range expansion on this same timescale across the SS (group E) and between Australian and Hawaiian regions of the western/central Pacific Ocean. Results of the mismatch distributions and occurrence of 4 low but statistically nonzero pairwise Φst values for the mitochondrial markers of G. undulatus (Table 3) are consistent with this hypothesis of transitory geographic genetic structure within the past 0.6 My. Based on the estimated coalescence times within the past 2.3 My in G. flavimarginatus and 5.9 My in G. undulatus, the prevailing pattern seems to be genetic homogeneity of populations over the longer (My) timeframe, with transitory population fragmentation on the shorter timescale of 0.3–0.6 My.
Recent controversy regarding NCPA (Garrick et al. 2008; Knowles 2008; Petit 2008; Templeton 2008) concerns whether significant genetic associations of haplotype groupings with geographic distances justify inferences of isolation by distance or contiguous range expansion, as these are the 2 most common inferences from data sets generated under panmictic conditions by Panchal and Beaumont (2007). It should be noted that the limitations described by Petit (2008) and Panchal and Beaumont (2007) have been corrected (Templeton 2008, 2009). Our results from NCPA involve 3- and 4-step clades with significant nonrandom associations of genetic and geographic distances. Our interpretations from NCPA are consistent with the development of transient, geographically localized population genetic structure on the order of a few hundred thousand years, and these are supported in part by significant pairwise Φst comparisons and mismatch distributions. At broader scales and the entire cladogram level, NCPA, AMOVA, Mantel tests, and pairwise Φst all yield the same result: inability to reject panmixia. At best, NCPA has identified evidence for population genetic structure at a level that is lost to more reductionist approaches such as AMOVA. At worst, it has yielded a false inference, and all populations of G. undulatus and G. flavimarginatus surveyed here are panmictic. In either case, haplotypes likely to have arisen within the past few hundred thousands of years have not all acquired ocean-wide distributions, but genetic continuity of populations is nonetheless the prevailing pattern when results are averaged over the multi-My history covered by mtDNA markers.
Klanten et al. (2007) apply the term “temporal genetic partitioning” to Indo-Pacific reef fishes to describe a pattern of population genetic homogeneity on a My timescale, with transitory fragmentation and demographic expansion on shorter timescales. “Each expansion episode is likely to have been accompanied by re-colonization of distant reef habitats resulting in repeated and widespread secondary contacts amongst previously isolated, genetically divergent demes” (Klanten et al. 2007). Global climatic fluctuations provide the most obvious source of recurring cycles of fragmentation followed by expansion of populations for reef fishes having widespread larval dispersal. Randall (1998) hypothesizes 3–6 episodes within the past 700000 years during which sea level would have been low enough to preclude dispersal of marine species between the Indian and Pacific oceans. A review of geological climatic data by Paillard (2006) indicates sea-level fluctuation occurring on a cycle of approximately 100000 years throughout the past My. These results are highly congruent with our conclusion that the 2 Gymnothorax species are emerging from a transitory fragmentation event, with some remnants of fragmentation and recent range expansions over the past 0.3–0.6 My shown by some of the younger mitochondrial haplotypes.
This interpretation is consistent with our evidence for population growth in the form of mismatch distributions and significant negative Fu's FS values for the mitochondrial loci of each species and for RAG-1 of G. undulatus. As proposed for other reef fishes and invertebrates, this population expansion likely followed the last period of low sea levels that exposed portions of the SS (Barber et al. 2006; Klanten et al. 2007; Gaither et al. 2010).
The Hawaiian Archipelago has a special biogeographic significance for moray eels, and this is apparent in several facets of our data. First, morays are common and abundant on these islands, one of the most isolated tropical archipelagos in the world, whereas common reef predators such as groupers (Serranidae) and snappers (Lutjanidae) are nearly absent. Thus, the high dispersal evident in our survey has allowed the morays to colonize and to proliferate in Hawaii, to the exclusion of many west/south Pacific competitors (Randall 2007). Second, Hawaii may be the biogeographic crossroads between the tropical eastern Pacific and western Pacific. This is indicated by the distribution of shared haplotypes from Hawaii to the western Pacific and Indian oceans (subgroup C-2 and group B in G. flavimarginatus, subgroup E-1 in G. undulatus), and Hawaii to the eastern Pacific (subgroups A-1 and C-1 in G. flavimarginatus; Figure 2, Table 4). Based on the latter network configuration, NCPA reveals dispersal from Hawaii to the eastern Pacific in the past 0.3 My (directionality inferred from haplotype network rooting and historical biogeography; see Briggs 1999). Many tropical species show very shallow phylogeographic histories in the eastern Pacific, indicating that much of this fauna is extirpated during glacial conditions (Bowen and Karl 2007). Hence, Hawaii may be the essential source of moray diversity in the East Pacific. Third, the only significant mtDNA partitions on a local scale for either species are within the Hawaiian Islands (Table 3, populations number 10 and 14), and one of them is between the high islands of the southeastern Hawaiian Archipelago and the low islands to the northwest. The latter region is now protected as the Papahānaumokuākea Marine National Monument. Emerging data indicate that the jurisdictional and ecological differences between these regions are matched by population genetic differences in some other reef species (Rivera et al. 2004; Ramon et al. 2008; Eble et al. 2009).
Gymnothorax flavimarginatus and G. undulatus represent extreme cases among reef fishes in maintaining gene flow across two-thirds of the planet. Several other fishes found on or near coral reefs maintain high levels of connectivity across similar spatial scales (e.g., whale sharks, Castro et al. 2007; tunas, Theisen et al. 2008), but these all have highly vagile adult forms that are not restricted to coral reef habitat. Our favored hypothesis for the geographic–genetic homogeneity of these Gymnothorax species is that their extended pelagic larval phase permits widespread dispersal of individuals after the removal of transitory oceanic barriers to pelagic dispersal. The genetic connectivity of both Gymnothorax species across the EPB and SS indicates an ability to disperse across even the strongest obstacles to other reef fishes.
This work was supported by a United States National Science Foundation grant (OCE-0454873) and the Hawaii Institute of Marine Biology Northwestern Hawaiian Islands Coral Reef Research Partnership grant (NMSP MOA 2005-008/6682) to B.W.B.; National Geographic Young Explorers Award, the National Science Foundation Doctoral Dissertation Improvement Grant (DEB 0909756), the Society of Systematic Biologists Mini-Partnerships for Enhancing Expertise in Taxonomy Award, the DeepFin PEET Award, and the Professional Association of Diving Instructors Foundation grant to J.S.R; summer undergraduate research fellowships from the Howard Hughes Medical Institute to V.G. and K.J.
For logistic assistance, expertise, and sound advice we thank Randall Kosaki, Jeff Eble, Jane Ball, Jo-Ann Leong, Matt Iacchei, Yannis Papastamatiou, Carl Meyer, Luiz Rocha, Matthew Craig, David Smith, Michelle Gaither, Robert Moffitt, Joe O'Malley, Sheila Conant, Senifa Annandale, Toby Daly-Engel, Zoltan Szabo, Kelly Gleason, Robert Toonen, James Lakey, Malia Rivera, Elizabeth Keenan, Derek Skillings, Tim Clark, the officers and crew of the R/V Hi'ialakai, the Stock Assessment Program at the NOAA Pacific Islands Fisheries Science Center, and Matt Furtado of Koolau Pets. Thanks to Kevin Conway for the leptocephalus illustration in Figure 1. We thank the Papahānaumokuākea Marine National Monument, US Fish and Wildlife Services, and Hawaii Division of Aquatic Resources for coordinating research activities and permitting procedures for collections in the Northwestern Hawaiian Islands. This is HIMB contribution no. 1367 and the School of Earth Science and Technology contribution no. 7825. Collections made in the Northwestern Hawaiian Islands National Monument were made under permits PMNM 2007-32 to B.W.B. and PMNM 2008-31 to J.S.R. This study complied with current laws in the United States and was conducted in accordance with the regulations of the Washington University and the University of Hawaii Institutional Animal Care and Use Committee.