How the Central American Seaway and an Ancient Northern Passage Affected Flatfish Diversification

How the Central American Seaway and an Ancient Northern Passage Affected Flatfish Diversification Abstract While the natural history of flatfish has been debated for decades, the mode of diversification of this biologically and economically important group has never been elucidated. To address this question, we assembled the largest molecular data set to date, covering > 300 species (out of ca. 800 extant), from 13 of the 14 known families over nine genes, and employed relaxed molecular clocks to uncover their patterns of diversification. As the fossil record of flatfish is contentious, we used sister species distributed on both sides of the American continent to calibrate clock models based on the closure of the Central American Seaway (CAS), and on their current species range. We show that flatfish diversified in two bouts, as species that are today distributed around the equator diverged during the closure of CAS, whereas those with a northern range diverged after this, hereby suggesting the existence of a postCAS closure dispersal for these northern species, most likely along a trans-Arctic northern route, a hypothesis fully compatible with paleogeographic reconstructions. Pleuronectiformes, Bayesian dating, trans-Arctic migration, vicariance, Central American Seaway, Isthmus of Panama Introduction The Pleuronectiformes, or flatfishes, are a speciose group of ray-finned fish, containing 14 families and over 800 known species (Munroe 2015). Flatfish begin life in the pelagic zone, but undergo a larval metamorphosis in which one eye, either left or right, depending on the species, migrates to the other side of the cranium. The adult fish then adopts a benthic lifestyle. Flatfish have asymmetric, laterally compressed bodies, and have lost their swim bladders during transformation. With eyes facing upwards, flatfish are also capable of protruding them. This singular morphology long puzzled taxonomists (Norman 1934), and the phylogeny of this group remains poorly resolved. At the highest taxonomic level, flatfishes are generally considered to be monophyletic, based on both morphological (Chapleau 1993), and molecular evidence (Berendzen et al. 2002; Pardo et al. 2005; Azevedo et al. 2008; Betancur-R et al. 2013; Harrington et al. 2016). All these studies also support the monophyletic status of most families within the order, to the exception of the Paralichthyidae. As all molecular studies to date have essentially focused on the monophyletic status of the order, they were based on as many representative species of each order. As a result, intra-ordinal relationships, among genera and even families, are still debated. It can therefore be expected that taking advantage of both species- and gene-rich evidence, while incorporating paleontological and/or geological data in the framework of molecular clocks, should help clarify not only the phylogenetic status of this family (dos Reis et al. 2016), but also—and more critically here—their evolutionary dynamics. However, very few flatfish fossils are known (Chanet 1997; Friedman 2012), and placed with confidence on the flatfish evolutionary tree (Parham et al. 2012; Campbell, López, et al. 2014; Harrington et al. 2016). As a result, calibrating a molecular clock becomes challenging. On the other hand, a dense species sampling may enable us to take advantage of a singular feature of flatfishes: some extant species are found both in the Pacific and Atlantic oceans. Furthermore, the existence of geminate species pairs of flatfishes, where sister taxa have one member in each ocean, suggests a speciation event predating the formation of the Isthmus of Panama, which occurred ∼12–3 Ma (Haug and Tiedemann 1998; O’Dea et al. 2016). Also possible is the role played by the opening of the Bering Strait around 5.5–5.4 Ma (Gladenkov et al. 2002), which would have permitted the trans-Arctic migration of ancestral populations from one ocean to the other, as documented in marine invertebrates (e.g., Durham and MacNeil 1967; Reid 1990) or vertebrates (e.g., Grant 1986), but to date, never in flatfish. Our driving hypothesis is then that the sole information relative to the formation of the Isthmus of Panama can be used to calibrate molecular clocks, and allows us to unravel the timing of flatfish evolution, as how rapidly they diversified remains an unsolved question. We show here that the diversification of flatfish in the seas surrounding the Americas followed a complex pattern, related to not only the closure of the Isthmus of Panama, but also to a warming event that opened up a trans-Arctic northern route. Results To test how the evolutionary dynamics of flatfish were affected by a major geological event, the closure of the Central American Seaway (CAS), we reconstructed dated Bayesian phylogenetic trees from a large data matrix under four relaxed molecular clock models, each one of them being based on a different calibration scheme (fig. 1). Under the first model, no calibration priors were placed on internal nodes. This initial tree, with the rogue sequences removed (data on GitHub: see Materials and Methods), was used to identify pairs of sister taxa that are split between the two oceans, with one species in the Atlantic and the other in the Pacific. This led us to single out twelve pairs of such species. These sister species happened not to be evenly distributed on the estimated phylogeny (fig. 2, top), but they are the only geminate species included in GenBank (as of August 2016). On each pair, we placed calibration date priors corresponding to the closure of the CAS (the “ALL” model). An examination of the posterior distributions of their divergence times suggests that some species have a very narrow speciation window, where all the mass of the posterior distribution is between 5 and 3 Ma, whereas others have a wider distribution (fig. 3A). Closer inspection of these distributions further reveals that most of the species with narrow posterior distributions have a northern range (figs. 2 and 3A, in blue), whereas those with the wider posterior distributions have a “southern” distribution, closer to the Isthmus of Panama (figs. 2 and 3A, in red). To further assess this observation, we first went back to the original clock model, removed all priors initially placed on sister taxa (the “NONE” model), and were able to validate that even in this case: northern and southern species showed, to one exception each (Hippoglossus hippoglossus and H. stenolepis in the north, and Poecilopsetta natalensis and P. hawaiiensis in the south; fig. 2), shifted posterior distributions (fig. 3B). The former pair was actually not estimated as being sister species in any of the four clock models (fig. 1), while P. natalensis and P. hawaiiensis, although inhabiting the Western Indian and Eastern Pacific oceans, occupy ranges that do not extend to the coasts of the Americas as with the other identified sister taxa. Models with priors placed only on northern (the “NORTH” model: fig. 3C) or southern (the “SOUTH” model: fig. 3D) species also showed a similar temporal shift. This shift suggested that southern species diverged early, before the complete closing of the CAS, whereas northern species diverged later, at or possibly after the isthmus was completed. Averaging these posterior distributions for the northern and southern species, to the exception of the two outliers noted above, showed these results more clearly (fig. 3E–H). To assess the robustness of this result to the inclusion of the P. natalensis and P. hawaiiensis sister species, we conducted an additional set of analyses removing all prior information from this pair. The ALL and SOUTH models, which are the only ones where this prior occurs, show that the general pattern that we found above holds: southern species diverged early, before the complete closing of the CAS, and northern species diverged later, at or possibly after the isthmus was completed (supplementary fig. S1, Supplementary Material online). Fig. 1. View largeDownload slide Phylogenetic trees reconstructed based on relaxed clock models. Four models were employed, representing different specifications of prior distributions set on sister taxa. (A) No priors were set on sister taxa. (B) Priors were set on all pairs of taxa. (C) Priors were set only on sister taxa with a current northern range. (D) Priors were set only on sister taxa with a current southern range. Horizontal scale is in million years ago (MYA). The closure of the CAS (12–3 Ma) is shown between vertical gray broken lines. Fig. 1. View largeDownload slide Phylogenetic trees reconstructed based on relaxed clock models. Four models were employed, representing different specifications of prior distributions set on sister taxa. (A) No priors were set on sister taxa. (B) Priors were set on all pairs of taxa. (C) Priors were set only on sister taxa with a current northern range. (D) Priors were set only on sister taxa with a current southern range. Horizontal scale is in million years ago (MYA). The closure of the CAS (12–3 Ma) is shown between vertical gray broken lines. Fig. 2. View largeDownload slide Distribution of the geminate species used in this study. Geographic distribution of the twelve pairs of sister species of flatfish used to calibrate the relaxed molecular clock models. Original data come from GBIF (www.gbif.org; accessed November 9, 2017). The top panel shows the phylogenetic distribution of these species based on the clock model including prior ages on all twelve pairs of species; scale bar: 5 Ma. Fig. 2. View largeDownload slide Distribution of the geminate species used in this study. Geographic distribution of the twelve pairs of sister species of flatfish used to calibrate the relaxed molecular clock models. Original data come from GBIF (www.gbif.org; accessed November 9, 2017). The top panel shows the phylogenetic distribution of these species based on the clock model including prior ages on all twelve pairs of species; scale bar: 5 Ma. Fig. 3. View largeDownload slide Posterior densities of divergence times for sister taxa used as calibration points in the relaxed molecular clock analyses, under four different calibration schemes. NONE: no priors were placed on sister taxa; ALL: priors were placed on all pairs of sister taxa; NORTH: priors were placed only on pairs of sister taxa with a northern distribution; SOUTH: priors were placed only on pairs of sister taxa with a southern distribution. Top row (A–D) shows all 17 distributions, whereas bottom row (E–H) shows range-averaged distributions (solid) to the exception of outlier pairs (dashed lines). Densities are color-coded for species with northern (blue) and southern (red) range. Dashed vertical gray lines indicate the closure of the CAS (21–3 Ma). Fig. 3. View largeDownload slide Posterior densities of divergence times for sister taxa used as calibration points in the relaxed molecular clock analyses, under four different calibration schemes. NONE: no priors were placed on sister taxa; ALL: priors were placed on all pairs of sister taxa; NORTH: priors were placed only on pairs of sister taxa with a northern distribution; SOUTH: priors were placed only on pairs of sister taxa with a southern distribution. Top row (A–D) shows all 17 distributions, whereas bottom row (E–H) shows range-averaged distributions (solid) to the exception of outlier pairs (dashed lines). Densities are color-coded for species with northern (blue) and southern (red) range. Dashed vertical gray lines indicate the closure of the CAS (21–3 Ma). All these results were obtained assuming that Psettodidae is an appropriate outgroup for these analyses. However, the position of Psettodidae in flatfishes is still debated (Betancur-R and Orti 2014; Campbell, López, et al. 2014; Campbell, Chen, et al. 2014; Harrington et al. 2016), so that their inclusion might affect our results. To assess this possibility, we reran all four models without Psettodidae, letting the molecular clock root the tree (Aris-Brosou and Rodrigue 2012). Supplementary figure S2, Supplementary Material online shows the exact same results as above (fig. 3), so that our dating results are also robust to the inclusion, or not, of Psettodidae. Our dating results are also robust to the inclusion of internal calibration priors as used in Harrington et al. (2016), and based on the Scophthalmus rhombus/Bothus pantherinus and the Crossorhombus kobensis/B. pantherinus divergences, dated at 29.62 (Vandenberghe et al. 2012) and 11.06 Ma (Carnevale et al. 2006), respectively (supplementary fig. S3, Supplementary Material online). In an attempt to tease out these models (including Psettodidae and without the two internal calibration priors) and their predictions about the exact timing of divergence between northern and southern species, we assessed model fit by means of the modified Akaike’s Information Criterion (AICM; Baele et al. 2013), which accounts for uncertainty in the MCMC sampling (Raftery et al. 2007). Even if model ranking based on this measure is known to be unstable (Baele et al. 2013), it is clear that the models with priors only on the northern or on the southern species perform significantly (>200 AIC units) worse than the two other models, which may be difficult to tease apart (table 1). The predictions of the best models suggest that H. hippoglossus and H. stenolepis, both northern species, consistently diverged before the complete closure of the seaway, in tandem with the average southern species, and that the average northern species diverged in tandem with P. natalensis and P. hawaiiensis (fig. 3E and F). Table 1. AICM Values for the Phylogenetic Analyses Using Four Different Calibration Schemes. Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Note.—Models are ranked from the best AICM model (top) to the worst model. SE: standard error. Table 1. AICM Values for the Phylogenetic Analyses Using Four Different Calibration Schemes. Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Note.—Models are ranked from the best AICM model (top) to the worst model. SE: standard error. One intriguing result from these analyses is that the northern sister species seem to come from a single clade (fig. 2), which would suggest that it is possible to identify the direction of the trans-Arctic migration. For this, we retrieved the current species range of all flatfish analyzed here (supplementary fig. S4, Supplementary Material online), and indicated their average latitude on the phylogenetic trees (supplementary fig. S5, Supplementary Material online). This northern clade, highly supported (posterior probability > 0.9; see asterisk in supplementary fig. S5, Supplementary Material online) even in the analyses including the two additional internal priors (supplementary fig. S6, Supplementary Material online), is mostly comprised of species currently found in the Pacific (e.g., Reinhardtius evermanni), suggesting that the northern sister species that we identified originated from the Pacific, and that their trans-Arctic migration was from the Pacific into the Atlantic in a northern direction. Discussion Our results show that flatfish underwent a first speciation at the time when the CAS closed, which led to the species that, today, have an equatorial range (fig. 2), as the formation of the Isthmus of Panama resulted in a barrier to gene flow leading to their speciation. Our results also show that species that today have a northern range (fig. 2) either emerged at the closure of the seaway, or after its closure. These timing estimates imply that this second bout of speciation was not caused by gene flow impeded by the closure of the seaway, but demand an interpretation involving trans-Arctic gene flow, through a northern route, before being interrupted, most likely this time by a climatic event. Critically, these conclusions are robust to the inclusion of outgroup species (Psettodidae), of outliers species (P. natalensis and P. hawaiiensis), and are consistently found across the four clock models implemented here, including the one (NONE) relying on no other fossil information than the 47.8–42.1 Ma ancestor to this order (Chanet 1997; Friedman 2012)—hereby indicating that our data are highly informative with respect to the timing of the divergence of these sister species (and rejecting our initial hypothesis) following their trans-CAS and trans-Arctic migrations. Trans-Arctic migrations from the Pacific to the Atlantic are not new, and have been described in both invertebrates (e.g., Durham and MacNeil 1967) and vertebrates (e.g., Grant 1986), but never using a relaxed molecular clock relying on as little information as a single fossil and a large multigene data set covering an entire order (Pleuronectiformes). While seminal work was usually based on paleontological (Durham and MacNeil 1967), or on morphological evidence (Reid 1990), subsequent work started using allozymes and other proteins (Grant 1986; Zaslavskaya et al. 1992), followed by genetic evidence, mostly at the population level using mitochondrial DNA (Rawson and Harper 2009). However, most of the recent work that relies on molecular clocks either makes strong assumptions on substitution rates (Rawson and Harper 2009; Liu et al. 2011), uses divergences that are vicariance-motivated but not supported by any fossil information (McCusker and Bentzen 2010), or posit in which direction the trans-Arctic migration took place (Dodson et al. 2007) to infer trans-Arctic migrations. Our work is the first to use sister species to demonstrate the existence and elucidate the timing of both the trans-CAS and the trans-Arctic events, while also inferring in which direction the latter took place, without making too many assumptions. Strikingly, the geological evidence is directly in line with our date estimates. The fossil record suggests that the first aquatic connection between the Pacific and Arctic (and Atlantic) oceans through the Bering Strait occurred ∼5.5–5.4 Ma (Gladenkov et al. 2002) due to a rise in sea levels, linked to tectonic activity (Marincovich and Gladenkov 2001). This would have permitted the trans-Arctic migration of populations ancestral to today’s northern species from one ocean to the other through this ancient “northern passage.” This global warming event, between the late Miocene to early Pleiocene, was then followed by a significant period of cooling during the Pleiocene into the Pleistocene (Zachos et al. 2001), leading to periods of repeated glaciations and a subsequent ice age. These cold events would have resulted in the closure of this ancient “northern passage,” hereby stopping the trans-Arctic gene flow between the two oceans, and leading to the more recent speciation of the northern taxa. For approximately a million years after the first opening of the Bering Strait, water flowed through the strait in a southern direction, from the Arctic to the Pacific ocean, until the formation of the Isthmus of Panama occurred close to the equator (Berta 2012). With the formation of the Isthmus, and hence the closing of the CAS, the ocean currents reversed due to a change in global ocean circulation (Haug and Tiedemann 1998; De Schepper et al. 2015), and have since flowed in a northern direction through the Bering Strait, from the Pacific to the Arctic (Marincovich 2000). This again is absolutely consistent with our results (supplementary fig. S4, Supplementary Material online), with dispersal or migration from the Pacific into the Atlantic in a northern direction through this strait, which is known as the trans-Arctic interchange (Vermeij 1991). Fossil data also show that the Bering land bridge has been exposed and submerged on multiple occasions since the Pleistocene (Gladenkov 2004). These openings and closings of the Bering Strait could have provided a mechanism for divergence and the evolution of sister taxa (Taylor and Dodson 1994; Väinölä 2003). Our results also have implication at the family level of flatfish. Further significant global cooling during the Pleistocene resulted in major glaciation events (Zachos et al. 2008) that could be responsible for creating barriers that isolated populations. All of the remaining sister taxa in our analysis, who have divergence estimates of <2 Ma in our study belong to the family Pleuronectidae (sensuCooper and Chapleau 1998). The Pleuronectidae are the predominant flatfish family found in cold and temperate seas of the northern hemisphere (Norman 1934; Cooper and Chapleau 1998). There are far more Pleuronectidae species in the Pacific Ocean, most of them endemic to the north Pacific Ocean off of North America and Asia in the region extending from the Bering Strait to the gulf of California (Norman 1934). None of the arctic species is restricted solely to arctic waters (Munroe 2015). Munroe also noted that Cooper (in an unpublished manuscript) identified areas of endemism among the current distribution of the Pleuronectidae. It has been shown that during the trans-Arctic interchange, there was a far higher number of species (up to eight times higher) migrating to the Atlantic than to the Pacific (Vermeij 1991). Fossil and phylogenetic evidence suggest the Pacific Ocean as the origin for diversification of the Pleuronectidae (Munroe 2015), and our phylogenetic results are highly congruent with this hypothesis. It is possible that the outliers, the Atlantic and Pacific halibut (H. hippoglossus and H. stenolepis, respectively), diverged during one of the first openings of the Bering Strait. The estimated dates from the molecular clock analysis are ∼5.5 Ma, in accordance with the hypothesized dates for the first aquatic connection (Marincovich and Gladenkov 2001). The remaining sister taxa have a younger age estimate of ∼1–2 Ma, corresponding with global cooling during the Pleistocene and repeated glaciations (Zachos et al. 2001), which likely formed more barriers to genetic flow. In the second pair of outliers, P. hawaiiensis inhabits waters of the Eastern Central Pacific to the Hawaiian Islands, whereas P. natalensis inhabits coastal waters of Eastern Africa (fig. 2). Because of the far reaching range of P. natalensis, and the relatively younger age estimates for divergence (1–2 Ma), these results beg for future research into other vicariant hypotheses, or dispersal routes, as these two members of Poecilopsetta diverged long after the Isthmus had closed. One candidate would be the Indo-Pacific barrier, which formed an almost continuous land bridge between Australia and Asia, and could have hence limited dispersal between the Pacific and the Indian oceans, hereby promoting speciation (Gaither et al. 2011). This barrier formed during the Pleistocene glacial cycles (2.58–0.01 Ma; Voris 2000), which is, again, fully consistent with our date estimates (fig. 3, supplementary fig. S1, Supplementary Material online). On the basis of the most extensive multigene sequence alignment available to date across all flatfish species, we have showed here that the evolutionary dynamics of sister species that are distributed across the two oceans strongly supported the existence of two bouts of speciation: one triggered by the closure of the CAS 12–3 Ma, and a second one due to the closure of an ancient northern passage 5–0.01 Ma. Our work adds to a growing body of evidence pointing not only to the role of geological processes in shaping biotic turnovers (Bacon et al. 2015), but also on the consequences of climate change that can either trigger speciation events (Reid 1990), or promote faunal interchanges whose ecological and economic impacts are difficult to predict (Wisz et al. 2015). Materials and Methods Data Retrieval and Alignment Prior to retrieving sequence data, GenBank was surveyed to identify all the genes belonging to species of Pleuronectiform families (as of August 2016), based on GenBank’s taxonomy browser. DNA sequences for a total of 332 flatfish species (out of over 800 species in the order) were identified and downloaded for five nuclear genes (kiaa1239, myh6, ripk4, rag1, sh3px3), and four mitochondrial genes (12s, 16s, cox1, and cytb). These represented all the taxa having at least one of these six gene sequences in GenBank; see supplementary table S1, Supplementary Material online for the corresponding accession numbers. Diversity was richly sampled as species from 13 of the 14 families in the order Pleuronectiformes were included in our catchment. In line with current consensus in flatfish systematics, the family Psettodidae was chosen as the outgroup to all other taxa (Chapleau 1993). These sequences were aligned using MUSCLE ver. 3.8.31 (Edgar 2004) on a gene-by-gene basis. Each alignment was visually inspected with AliView ver. 1.18 (Larsson 2014), and was manually edited where necessary. In particular, large indels (> 10 bp) were removed prior to all phylogenetic analyses. The 5′ and 3′ ends of sequences were also trimmed. The aligned sequences were then concatenated into a single alignment using a custom R script. Data Preprocessing To gauge the phylogenetic content of our data set, we performed a first series of molecular clock analyses on the concatenated data matrix, with all the nine gene sequences obtained above (12S, 16S, COI, Cyt-b, KIAA1239, MYH6, RIPK4, SH3PX3, and RAG1), and all the 332 taxa representing all families of flatfish. To account for rate heterogeneity along this alignment, partitions were first identified with PartitionFinder 2 (Lanfear et al. 2017), selecting the best-fit model of nucleotide substitution based on Akaike’s Information Criterion (AIC). This optimal partition scheme was then employed in a first Bayesian phylogenetic analysis, conducted with BEAST ver. 1.8.3 (Drummond and Rambaut 2007). This program implements a Markov chain Monte Carlo (MCMC) sampler, which co-estimates both tree topology and divergence times. As determined by PartionFinder, a GTR + I + Γ model was applied to each data partition. These models of evolution across partitions were assumed to be independent (unlinked, in BEAST parlance), whereas both clock model and tree model partitions were shared (linked) by all partitions. An uncorrelated relaxed clock was employed with a lognormal distribution prior on rates, and a Yule speciation prior (Drummond et al. 2006). Because of the paucity of reliably placed fossils on the flatfish tree, a unique calibration point was placed on the most recent common ancestor (MRCA) of the ingroup, as a mean-one exponential prior, with an offset of 40 My reflecting the age of 47.8–42.1 Ma (Chanet 1997; Friedman 2012). To stabilize the estimate, a lognormal ln(0, 1.5) prior with a 40 Ma offset was also placed on the root of the tree (root height). Two separate MCMC samplers were run, each for 10,000,000 generations. Trees were sampled every 5,000 generations, and convergence was checked visually using Tracer ver. 1.6 (available at: http://tree.bio.ed.ac.uk/software/tracer/; last accessed May 24, 2018). Tree log files from each run were combined in LogCombiner, after conservatively removing 10% of each run as burn-in. The resulting maximum a posteriori (MAP) tree was generated with TreeAnnotator (Drummond and Rambaut 2007). As the topology of this resulting MAP tree was unconventional, we suspected the presence of rogue taxa, that is, species evolving either much faster or much slower than the majority, which can contribute negatively to consensus tree support (Aberer et al. 2013). Rogue taxa were identified using the RogueNaRok (Aberer et al. 2013) webserver (http://rnr.h-its.org; last accessed May 24, 2018). The consensus trees from the preliminary analysis using 10,000 iterations were used as the tree set for rogue taxon analysis. A total of 22 rogues were identified and pruned from the analysis, leaving 310 species. Molecular Dating To assess the impact of the closure of the CAS on flatfish evolutionary dynamics, a second set of partitioned relaxed molecular clock analyses was performed (without the rogue sequences). The timing of the closure of the CAS is estimated to have occurred between 12 and 3 Ma (Duque-Caro 1990; Coates et al. 1992; Haug and Tiedemann 1998; O'Dea et al. 2016), and we used this time window as a prior to inform the relaxed molecular clock-based phylogenetic reconstructions. The analyses were performed on the same concatenated data set, with the same single fossil calibration, but we also placed a lognormal prior ln(3, 1.5), that has most of its mass on the 12–3 Ma window, on the MRCA of each pair of sister taxa (the “ALL” model). From the initial BEAST analyses, twelve pairs of taxa were selected based on two criteria: 1) being sister species on that initial tree, 2) with one species distributed in the Pacific and one in the Atlantic ocean (fig. 2). Again, two independent MCMC samplers were run each for 100 million iterations, with samples taken every 5000 step. Because these pairs of sister species show a contrasted geographic distribution, having either a southern (equatorial) or a northern range (fig. 2), two additional sets of analyses were performed. In a first set, calibration priors (ln(3, 1.5), as per above) were placed only sister taxa that had a geographic range in the southern hemisphere (the “SOUTH” model), while in a second set, identical calibration priors were placed only on sister species with a northern range (the “NORTH” model). Finally, a set of analyses was performed using no sister taxa calibrations at all (the “NONE” model). For each analysis, results from the two MCMC runs were combined using LogCombiner after removing an even more conservative burn-in of 50%. The final MAP tree was generated with TreeAnnotator. In an attempt to rank these different models (priors on all sister taxa; only on southern taxa; only on southern taxa; no “CAS” priors), the AICM was computed for each model (Baele et al. 2013). These computations were performed in Tracer for each of the four different calibration models, based on 100 replicates. Species Distribution Data In order to map the distribution of each species of flatfish, we resorted to the Global Biodiversity Information Facility (GBIF: https://www.gbif.org/; last accessed May 24, 2018) database, accessed with the R library rgbif (Chamberlain 2017). Up to 1,000 records were retrieved for each species, and plotted on a geographic map based on GBIF observational records. Where needed, species locations were summarized by mean latitude and longitude (supplementary fig. S4, Supplementary Material online). Approximate attribution to a particular ocean was based on species-specific mean longitudes (Atlantic: between 20° and 120°; Indian: 120° and −100°; Pacific otherwise). Availability of Computer Code and Data Accession numbers, sequence alignments, BEAST models (as xml files), estimated phylogenetic trees (MAP trees from BEAST), and the R scripts used to plot the figures in this study are available from https://github.com/sarisbro, last accessed May 24, 2018. Supplementary Material Supplementary data are available at Molecular Biology and Evolution online. Acknowledgments This work was supported by the University of Ottawa (L.B.), and the Natural Sciences Research Council of Canada (F.C., S.A.B.). We are grateful to Jonathan Dench for discussions and comments, to two anonymous reviewers for very constructive comments, and to Compute Canada and Ontario’s Center for Advanced Computing for providing us access to their resource. This work was completed while S.A.B. was hosted by Prof. Yutaka Watanuki, at the University of Hokkaido in Hakodate, thanks to an Invitational Fellowship from the Japanese Society for the Promotion of Science. References Aberer AJ , Krompass D , Stamatakis A. 2013 . Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice . Syst Biol . 62 1 : 162 – 166 . Google Scholar CrossRef Search ADS PubMed Aris-Brosou S , Rodrigue N. 2012 . The essentials of computational molecular evolution. In evolutionary genomics: statistical and computational methods. In: Anisimova M , editor. Methods in molecular biology (methods and protocols) , Vol. 855 . Totowa (NJ ): Humana Press . Azevedo MFC , Oliveira C , Pardo BG , Martínez P , Foresti F. 2008 . Phylogenetic analysis of the order Pleuronectiformes (Teleostei) based on sequences of 12S and 16S mitochondrial genes . Genet Mol Biol . 31(1 Suppl) : 284 – 292 . Google Scholar CrossRef Search ADS Bacon CD , Silvestro D , Jaramillo C , Smith BT , Chakrabarty P , Antonelli A. 2015 . Biological evidence supports an early and complex emergence of the Isthmus of Panama . Proc Natl Acad Sci U S A . 112 19 : 6110 – 6115 . Google Scholar CrossRef Search ADS PubMed Baele G , Li WLS , Drummond AJ , Suchard MA , Lemey P. 2013 . Accurate model selection of relaxed molecular clocks in Bayesian phylogenetics . Mol Biol Evol . 30 2 : 239 – 243 . Google Scholar CrossRef Search ADS PubMed Berendzen PB , Dimmick WW , McEachran JD. 2002 . Phylogenetic relationships of Pleuronectiformes based on molecular evidence . Copeia 2002 : 642 – 652 . Google Scholar CrossRef Search ADS Berta A. 2012 . Return to the sea: the life and evolutionary times of marine mammals . Berkeley (CA ): University of California Press . Google Scholar CrossRef Search ADS Betancur-R R , Ortí G. 2014 . Molecular evidence for the monophyly of flatfishes (Carangimorpharia: pleuronectiformes) . Mol Phylogenet Evol . 73 : 18 – 22 . Google Scholar CrossRef Search ADS PubMed Betancur-R R , Li C , Munroe TA , Ballesteros JA , Ortí G. 2013 . Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: pleuronectiformes) . Syst Biol . 62 5 : 763 – 785 . Google Scholar CrossRef Search ADS PubMed Campbell J , López JA , Satoh TP , Chen WJ , Miya M. 2014 . Mitochondrial genomic investigation of flatfish monophyly . Gene 551 2 : 176 – 182 . Google Scholar CrossRef Search ADS PubMed Campbell MA , Chen WJ , López JA. 2014 . Molecular data do not provide unambiguous support for the monophyly of flatfishes (Pleuronectiformes): a reply to Betancur-R and Ortí . Mol Phylogenet Evol . 75 : 149 – 153 . Google Scholar CrossRef Search ADS PubMed Carnevale G , Bannikov AF , Landini W , Sorbini C. 2006 . Volhynian (early Sarmatian sensu lato) fishes from Tsurevsky, north Caucasus, Russia . J Paleontol . 80 : 684 – 699 . Google Scholar CrossRef Search ADS Chamberlain S. 2017 . rgbif: Interface to the Global ‘Biodiversity’ Information Facility API. R package version 0.9.9. https://CRAN.R-project.org/package=rgbif Chanet B. 1997 . A cladistic reappraisal of the fossil flatfishes record: consequences on the phylogeny of the Pleuronectiformes (Osteichthyes: teleostei) . Ann Des Sci Nat Zool Paris . 18 : 105 – 117 . Chapleau F. 1993 . Pleuronectiform relationships: a cladistic reassessment . Bull Mar Sci . 52 : 516 – 540 . Coates AG , Jackson JBC , Collins LS , Cronin TM , Dowsett J , Bybell LM , Jung P , Obando JA. 1992 . Closure of the Isthmus of Panama: the near-shore marine record of Costa Rica and western Panama . Geol Soc Am Bull . 104 7 : 814 – 828 . Google Scholar CrossRef Search ADS Cooper JA , Chapleau F. 1998 . Monophyly and intrarelationships of the family Pleuronectidae (Pleuronectiformes), with a revised classification . Fish Bull . 96 : 686 – 726 . De Schepper S , Schreck M , Beck KM , Matthiessen J , Fahl K , Mangerud G. 2015 . Early Pliocene onset of modern Nordic Seas circulation related to ocean gateway changes . Nat Commun . 6 : 8659. Google Scholar CrossRef Search ADS PubMed Dodson JJ , Tremblay S , Colombani F , Carscadden JE , Lecomte F. 2007 . Trans-Arctic dispersals and the evolution of a circumpolar marine fish species complex, the capelin (Mallotus villosus) . Mol Ecol . 16 23 : 5030 – 5043 . Google Scholar CrossRef Search ADS PubMed dos Reis M , Donoghue PCJ , Yang Z. 2016 . Bayesian molecular clock dating of species divergences in the genomics era . Nat Rev Genet . 17 2 : 71 – 80 . Google Scholar CrossRef Search ADS PubMed Drummond AJ , Ho SYW , Phillips MJ , Rambaut A. 2006 . Relaxed phylogenetics and dating with confidence . PLoS Biol . 4 5 : e88 – 710 . Google Scholar CrossRef Search ADS PubMed Drummond AJ , Rambaut A. 2007 . BEAST: Bayesian evolutionary analysis by sampling trees . BMC Evol Biol . 7 : 214. Google Scholar CrossRef Search ADS PubMed Duque-Caro H. 1990 . Neogene stratigraphy, paleoceanography and paleobiology in northwest South America and the evolution of the Panama Seaway . Palaeogeogr Palaeoclimatol Palaeoecol . 77 ( 3–4 ): 203 – 234 . Google Scholar CrossRef Search ADS Durham JW , MacNeil FS. 1967 . Cenozoic migrations of marine invertebrates through the Bering Strait region. In: Hopkins DM , editor. The Bering land bridge. Stanford (CA ), Stanford University Press . p. 326 – 349 . Edgar R. 2004 . MUSCLE: multiple sequence alignment with high accuracy and high throughput . Nucleic Acids Res . 32 5 : 1792 – 1797 . Google Scholar CrossRef Search ADS PubMed Friedman M. 2012 . Osteology of †Heteronectes chaneti (Acanthomorpha, Pleuronectiformes), an Eocene stem flatfish, with a discussion of flatfish sister-group relationships . J Vert Paleontol . 32 4 : 735 – 756 . Google Scholar CrossRef Search ADS Gaither MR , Bowen BW , Bordenave TR , Rocha LA , Newman SJ , Gomez JA , van Herwerden L , Craig MT. 2011 . Phylogeography of the reef fish Cephalopholis argus (Epinephelidae) indicates Pleistocene isolation across the indo-pacific barrier with contemporary overlap in the coral triangle . BMC Evol Biol . 11 : 189. Google Scholar CrossRef Search ADS PubMed Gladenkov AY , Oleinik AE , Marincovich L , Barinov KB. 2002 . A refined age for the earliest opening of Bering Strait . Palaeogeogr Palaeoclimatol Palaeoecol . 183 ( 3–4 ): 321 – 328 . Google Scholar CrossRef Search ADS Gladenkov AY. 2004 . Onset of connections between the Pacific and Arctic Oceans through the Bering Strait in the Neogene . Stratigr Geol Correl . 12 : 175 – 187 . Grant WS. 1986 . Biochemical genetic divergence between Atlantic, Clupea harengus, and Pacific, C. pallasi, Herring . Copeia 1986 3 : 714 – 719 . Google Scholar CrossRef Search ADS Harrington RC , Faircloth BC , Eytan RI , Smith WL , Near TJ , Alfaro ME , Friedman M. 2016 . Phylogenomic analysis of carangimorph fishes reveals flatfish asymmetry arose in a blink of the evolutionary eye . BMC Evol Biol . 16 1 : 224. Google Scholar CrossRef Search ADS PubMed Haug GH , Tiedemann R. 1998 . Effect of the formation of the Isthmus of Panama on Atlantic Ocean thermohaline circulation . Nature 393 6686 : 673 – 676 . Google Scholar CrossRef Search ADS Lanfear R , Frandsen PB , Wright AM , Senfeld T , Calcott B . 2017 . PartitionFinder 2: New Methods for Selecting Partitioned Models of Evolution for Molecular and Morphological Phylogenetic Analyses . Mol Biol Evol. 34 3 : 772 – 773 . Google Scholar PubMed Larsson A. 2014 . AliView: a fast and lightweight alignment viewer and editor for large datasets . Bioinformatics 30 22 : 3276 – 3278 . Google Scholar CrossRef Search ADS PubMed Liu JX , Tatarenkov A , Beacham TD , Gorbachev V , Wildes S , Avise JC. 2011 . Effects of Pleistocene climatic fluctuations on the phylogeographic and demographic histories of Pacific herring (Clupea pallasii) . Mol Ecol . 20 18 : 3879 – 3893 . Google Scholar CrossRef Search ADS PubMed Marincovich L , Gladenkov AY. 2001 . New evidence for the age of Bering Strait . Quat Sci Rev . 20 ( 1–3 ): 329 – 335 . Google Scholar CrossRef Search ADS Marincovich L. 2000 . Central American paleogeography controlled Pliocene Arctic Ocean molluscan migrations . Geology 28 6 : 551 – 554 . Google Scholar CrossRef Search ADS McCusker MR , Bentzen P. 2010 . Phylogeography of 3 North Atlantic Wolffish species (Anarhichas spp.) with phylogenetic relationships within the family Anarhichadidae . J Heredity . 101 5 : 591 – 601 . Google Scholar CrossRef Search ADS Munroe TA. 2015 . Systematic diversity of the flatfishes. In: Gibson RN , Nash RDM , Geffen AJ , van der Veer HW , editor. Flatfishes: biology and exploitation. Chichester, West Sussex : John Wiley & Sons . p. 13 – 44 . Norman JR. 1934 . A systematic monograph of the flatfishes (Heterosomata). London : Printed by Order of the British Museum . Google Scholar CrossRef Search ADS O'Dea A , Lessios HA , Coates AG , Eytan RI , Restrepo-Moreno SA , Cione AL , Collins LS , de Queiroz A , Farris DW , Norris RD. 2016 . Formation of the Isthmus of Panama . Sci Adv . 2 8 : e1600883. Google Scholar CrossRef Search ADS PubMed Pardo BG , Machordom A , Foresti F , Porto-Foresti F , Azevedo MFC , Bañon R , Sánchez L , Martínez P. 2005 . Phylogenetic analysis of flatfish (Order Pleuronectiformes) based on mitochondrial 16S rDNA sequences . Sci Mar . 69 4 : 531 – 543 . Google Scholar CrossRef Search ADS Parham JF , Donoghue PCJ , Bell CJ , Calway TD , Head JJ , Holroyd PA , Inoue JG , Irmis RB , Joyce WG , Ksepka DT , et al. . 2012 . Best practices for justifying fossil calibrations . Syst Biol . 61 2 : 346 – 359 . Google Scholar CrossRef Search ADS PubMed Raftery AE , Newton MA , Satagopan JM , Krivitsky PN. 2007 . Estimating the integrated likelihood via posterior simulation using the Harmonic Mean Identity . Bayesian Stat . 8 : 1 – 45 . Rawson PD , Harper FM. 2009 . Colonization of the northwest Atlantic by the blue mussel, Mytilus trossulus postdates the last glacial maximum . Mar Biol . 156 9 : 1857 – 1868 . Google Scholar CrossRef Search ADS Reid DG. 1990 . Trans-arctic migration and speciation induced by climatic change: the biogeography of Littorina (Mollusca: gastropoda) . Bull Mar Sci . 47 : 35 – 49 . Taylor EB , Dodson JJ. 1994 . A molecular analysis biogeography within (genus Osmevus) of relationships and a species complex of Holarctic fish . Mol Ecol . 3 3 : 235 – 248 . Google Scholar CrossRef Search ADS PubMed Väinölä R. 2003 . Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex . Mar Biol . 143 5 : 935 – 946 . Google Scholar CrossRef Search ADS Vandenberghe N , Hilgen FJ , Speijer RP , Ogg JG , Gradstein FM , Hammer O , Hollis CJ , Hooker JJ. 2012 . The Paleogene period. In: Gradstein FM , Ogg JG , Schmitz MD , Ogg GM , editors. The geologic timescale . Amsterdam (The Netherlands ), Elsevier . p. 855 – 921 . Vermeij GJ. 1991 . Anatomy of an invasion: the trans-Arctic interchange . Paleobiology 17 03 : 281 – 307 . Google Scholar CrossRef Search ADS Voris HK. 2000 . Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations . J Biogeogr . 27 5 : 1153 – 1167 . Google Scholar CrossRef Search ADS Wisz MS , Broennimann O , Grønkjær P , Møller PR , Olsen SM , Swingedouw D , Hedeholm RB , Nielsen EE , Guisan A , Pellissier L. 2015 . Arctic warming will promote Atlantic-Pacific fish interchange . Nat Climate Change . 5 9 :790–265. Zachos J , Pagani M , Sloan L , Thomas E , Billups K. 2001 . Trends, rhythms, and aberrations in global climate 65 Ma to present . Science 292 5517 : 686 – 693 . Google Scholar CrossRef Search ADS PubMed Zachos JC , Dickens GR , Zeebe RE. 2008 . An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics . Nature 451 7176 : 279 – 283 . Google Scholar CrossRef Search ADS PubMed Zaslavskaya NI , Sergievsky SO , Tatarenkov AN. 1992 . Allozyme similarity of Atlantic and Pacific species of Littorina (Gastropoda: littorinidae) . J Moll Stud . 58 4 : 377 – 384 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Molecular Biology and Evolution Oxford University Press

How the Central American Seaway and an Ancient Northern Passage Affected Flatfish Diversification

Loading next page...
 
/lp/ou_press/how-the-central-american-seaway-and-an-ancient-northern-passage-AGdrVVMd1G
Publisher
Oxford University Press
Copyright
© The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com
ISSN
0737-4038
eISSN
1537-1719
D.O.I.
10.1093/molbev/msy104
Publisher site
See Article on Publisher Site

Abstract

Abstract While the natural history of flatfish has been debated for decades, the mode of diversification of this biologically and economically important group has never been elucidated. To address this question, we assembled the largest molecular data set to date, covering > 300 species (out of ca. 800 extant), from 13 of the 14 known families over nine genes, and employed relaxed molecular clocks to uncover their patterns of diversification. As the fossil record of flatfish is contentious, we used sister species distributed on both sides of the American continent to calibrate clock models based on the closure of the Central American Seaway (CAS), and on their current species range. We show that flatfish diversified in two bouts, as species that are today distributed around the equator diverged during the closure of CAS, whereas those with a northern range diverged after this, hereby suggesting the existence of a postCAS closure dispersal for these northern species, most likely along a trans-Arctic northern route, a hypothesis fully compatible with paleogeographic reconstructions. Pleuronectiformes, Bayesian dating, trans-Arctic migration, vicariance, Central American Seaway, Isthmus of Panama Introduction The Pleuronectiformes, or flatfishes, are a speciose group of ray-finned fish, containing 14 families and over 800 known species (Munroe 2015). Flatfish begin life in the pelagic zone, but undergo a larval metamorphosis in which one eye, either left or right, depending on the species, migrates to the other side of the cranium. The adult fish then adopts a benthic lifestyle. Flatfish have asymmetric, laterally compressed bodies, and have lost their swim bladders during transformation. With eyes facing upwards, flatfish are also capable of protruding them. This singular morphology long puzzled taxonomists (Norman 1934), and the phylogeny of this group remains poorly resolved. At the highest taxonomic level, flatfishes are generally considered to be monophyletic, based on both morphological (Chapleau 1993), and molecular evidence (Berendzen et al. 2002; Pardo et al. 2005; Azevedo et al. 2008; Betancur-R et al. 2013; Harrington et al. 2016). All these studies also support the monophyletic status of most families within the order, to the exception of the Paralichthyidae. As all molecular studies to date have essentially focused on the monophyletic status of the order, they were based on as many representative species of each order. As a result, intra-ordinal relationships, among genera and even families, are still debated. It can therefore be expected that taking advantage of both species- and gene-rich evidence, while incorporating paleontological and/or geological data in the framework of molecular clocks, should help clarify not only the phylogenetic status of this family (dos Reis et al. 2016), but also—and more critically here—their evolutionary dynamics. However, very few flatfish fossils are known (Chanet 1997; Friedman 2012), and placed with confidence on the flatfish evolutionary tree (Parham et al. 2012; Campbell, López, et al. 2014; Harrington et al. 2016). As a result, calibrating a molecular clock becomes challenging. On the other hand, a dense species sampling may enable us to take advantage of a singular feature of flatfishes: some extant species are found both in the Pacific and Atlantic oceans. Furthermore, the existence of geminate species pairs of flatfishes, where sister taxa have one member in each ocean, suggests a speciation event predating the formation of the Isthmus of Panama, which occurred ∼12–3 Ma (Haug and Tiedemann 1998; O’Dea et al. 2016). Also possible is the role played by the opening of the Bering Strait around 5.5–5.4 Ma (Gladenkov et al. 2002), which would have permitted the trans-Arctic migration of ancestral populations from one ocean to the other, as documented in marine invertebrates (e.g., Durham and MacNeil 1967; Reid 1990) or vertebrates (e.g., Grant 1986), but to date, never in flatfish. Our driving hypothesis is then that the sole information relative to the formation of the Isthmus of Panama can be used to calibrate molecular clocks, and allows us to unravel the timing of flatfish evolution, as how rapidly they diversified remains an unsolved question. We show here that the diversification of flatfish in the seas surrounding the Americas followed a complex pattern, related to not only the closure of the Isthmus of Panama, but also to a warming event that opened up a trans-Arctic northern route. Results To test how the evolutionary dynamics of flatfish were affected by a major geological event, the closure of the Central American Seaway (CAS), we reconstructed dated Bayesian phylogenetic trees from a large data matrix under four relaxed molecular clock models, each one of them being based on a different calibration scheme (fig. 1). Under the first model, no calibration priors were placed on internal nodes. This initial tree, with the rogue sequences removed (data on GitHub: see Materials and Methods), was used to identify pairs of sister taxa that are split between the two oceans, with one species in the Atlantic and the other in the Pacific. This led us to single out twelve pairs of such species. These sister species happened not to be evenly distributed on the estimated phylogeny (fig. 2, top), but they are the only geminate species included in GenBank (as of August 2016). On each pair, we placed calibration date priors corresponding to the closure of the CAS (the “ALL” model). An examination of the posterior distributions of their divergence times suggests that some species have a very narrow speciation window, where all the mass of the posterior distribution is between 5 and 3 Ma, whereas others have a wider distribution (fig. 3A). Closer inspection of these distributions further reveals that most of the species with narrow posterior distributions have a northern range (figs. 2 and 3A, in blue), whereas those with the wider posterior distributions have a “southern” distribution, closer to the Isthmus of Panama (figs. 2 and 3A, in red). To further assess this observation, we first went back to the original clock model, removed all priors initially placed on sister taxa (the “NONE” model), and were able to validate that even in this case: northern and southern species showed, to one exception each (Hippoglossus hippoglossus and H. stenolepis in the north, and Poecilopsetta natalensis and P. hawaiiensis in the south; fig. 2), shifted posterior distributions (fig. 3B). The former pair was actually not estimated as being sister species in any of the four clock models (fig. 1), while P. natalensis and P. hawaiiensis, although inhabiting the Western Indian and Eastern Pacific oceans, occupy ranges that do not extend to the coasts of the Americas as with the other identified sister taxa. Models with priors placed only on northern (the “NORTH” model: fig. 3C) or southern (the “SOUTH” model: fig. 3D) species also showed a similar temporal shift. This shift suggested that southern species diverged early, before the complete closing of the CAS, whereas northern species diverged later, at or possibly after the isthmus was completed. Averaging these posterior distributions for the northern and southern species, to the exception of the two outliers noted above, showed these results more clearly (fig. 3E–H). To assess the robustness of this result to the inclusion of the P. natalensis and P. hawaiiensis sister species, we conducted an additional set of analyses removing all prior information from this pair. The ALL and SOUTH models, which are the only ones where this prior occurs, show that the general pattern that we found above holds: southern species diverged early, before the complete closing of the CAS, and northern species diverged later, at or possibly after the isthmus was completed (supplementary fig. S1, Supplementary Material online). Fig. 1. View largeDownload slide Phylogenetic trees reconstructed based on relaxed clock models. Four models were employed, representing different specifications of prior distributions set on sister taxa. (A) No priors were set on sister taxa. (B) Priors were set on all pairs of taxa. (C) Priors were set only on sister taxa with a current northern range. (D) Priors were set only on sister taxa with a current southern range. Horizontal scale is in million years ago (MYA). The closure of the CAS (12–3 Ma) is shown between vertical gray broken lines. Fig. 1. View largeDownload slide Phylogenetic trees reconstructed based on relaxed clock models. Four models were employed, representing different specifications of prior distributions set on sister taxa. (A) No priors were set on sister taxa. (B) Priors were set on all pairs of taxa. (C) Priors were set only on sister taxa with a current northern range. (D) Priors were set only on sister taxa with a current southern range. Horizontal scale is in million years ago (MYA). The closure of the CAS (12–3 Ma) is shown between vertical gray broken lines. Fig. 2. View largeDownload slide Distribution of the geminate species used in this study. Geographic distribution of the twelve pairs of sister species of flatfish used to calibrate the relaxed molecular clock models. Original data come from GBIF (www.gbif.org; accessed November 9, 2017). The top panel shows the phylogenetic distribution of these species based on the clock model including prior ages on all twelve pairs of species; scale bar: 5 Ma. Fig. 2. View largeDownload slide Distribution of the geminate species used in this study. Geographic distribution of the twelve pairs of sister species of flatfish used to calibrate the relaxed molecular clock models. Original data come from GBIF (www.gbif.org; accessed November 9, 2017). The top panel shows the phylogenetic distribution of these species based on the clock model including prior ages on all twelve pairs of species; scale bar: 5 Ma. Fig. 3. View largeDownload slide Posterior densities of divergence times for sister taxa used as calibration points in the relaxed molecular clock analyses, under four different calibration schemes. NONE: no priors were placed on sister taxa; ALL: priors were placed on all pairs of sister taxa; NORTH: priors were placed only on pairs of sister taxa with a northern distribution; SOUTH: priors were placed only on pairs of sister taxa with a southern distribution. Top row (A–D) shows all 17 distributions, whereas bottom row (E–H) shows range-averaged distributions (solid) to the exception of outlier pairs (dashed lines). Densities are color-coded for species with northern (blue) and southern (red) range. Dashed vertical gray lines indicate the closure of the CAS (21–3 Ma). Fig. 3. View largeDownload slide Posterior densities of divergence times for sister taxa used as calibration points in the relaxed molecular clock analyses, under four different calibration schemes. NONE: no priors were placed on sister taxa; ALL: priors were placed on all pairs of sister taxa; NORTH: priors were placed only on pairs of sister taxa with a northern distribution; SOUTH: priors were placed only on pairs of sister taxa with a southern distribution. Top row (A–D) shows all 17 distributions, whereas bottom row (E–H) shows range-averaged distributions (solid) to the exception of outlier pairs (dashed lines). Densities are color-coded for species with northern (blue) and southern (red) range. Dashed vertical gray lines indicate the closure of the CAS (21–3 Ma). All these results were obtained assuming that Psettodidae is an appropriate outgroup for these analyses. However, the position of Psettodidae in flatfishes is still debated (Betancur-R and Orti 2014; Campbell, López, et al. 2014; Campbell, Chen, et al. 2014; Harrington et al. 2016), so that their inclusion might affect our results. To assess this possibility, we reran all four models without Psettodidae, letting the molecular clock root the tree (Aris-Brosou and Rodrigue 2012). Supplementary figure S2, Supplementary Material online shows the exact same results as above (fig. 3), so that our dating results are also robust to the inclusion, or not, of Psettodidae. Our dating results are also robust to the inclusion of internal calibration priors as used in Harrington et al. (2016), and based on the Scophthalmus rhombus/Bothus pantherinus and the Crossorhombus kobensis/B. pantherinus divergences, dated at 29.62 (Vandenberghe et al. 2012) and 11.06 Ma (Carnevale et al. 2006), respectively (supplementary fig. S3, Supplementary Material online). In an attempt to tease out these models (including Psettodidae and without the two internal calibration priors) and their predictions about the exact timing of divergence between northern and southern species, we assessed model fit by means of the modified Akaike’s Information Criterion (AICM; Baele et al. 2013), which accounts for uncertainty in the MCMC sampling (Raftery et al. 2007). Even if model ranking based on this measure is known to be unstable (Baele et al. 2013), it is clear that the models with priors only on the northern or on the southern species perform significantly (>200 AIC units) worse than the two other models, which may be difficult to tease apart (table 1). The predictions of the best models suggest that H. hippoglossus and H. stenolepis, both northern species, consistently diverged before the complete closure of the seaway, in tandem with the average southern species, and that the average northern species diverged in tandem with P. natalensis and P. hawaiiensis (fig. 3E and F). Table 1. AICM Values for the Phylogenetic Analyses Using Four Different Calibration Schemes. Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Note.—Models are ranked from the best AICM model (top) to the worst model. SE: standard error. Table 1. AICM Values for the Phylogenetic Analyses Using Four Different Calibration Schemes. Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Calibration Model AICM Value SE ALL 374,082.199 ±1.579 NONE 374,110.488 ±3.789 SOUTH 374,367.749 ±0.696 NORTH 374,746.628 ±3.789 Note.—Models are ranked from the best AICM model (top) to the worst model. SE: standard error. One intriguing result from these analyses is that the northern sister species seem to come from a single clade (fig. 2), which would suggest that it is possible to identify the direction of the trans-Arctic migration. For this, we retrieved the current species range of all flatfish analyzed here (supplementary fig. S4, Supplementary Material online), and indicated their average latitude on the phylogenetic trees (supplementary fig. S5, Supplementary Material online). This northern clade, highly supported (posterior probability > 0.9; see asterisk in supplementary fig. S5, Supplementary Material online) even in the analyses including the two additional internal priors (supplementary fig. S6, Supplementary Material online), is mostly comprised of species currently found in the Pacific (e.g., Reinhardtius evermanni), suggesting that the northern sister species that we identified originated from the Pacific, and that their trans-Arctic migration was from the Pacific into the Atlantic in a northern direction. Discussion Our results show that flatfish underwent a first speciation at the time when the CAS closed, which led to the species that, today, have an equatorial range (fig. 2), as the formation of the Isthmus of Panama resulted in a barrier to gene flow leading to their speciation. Our results also show that species that today have a northern range (fig. 2) either emerged at the closure of the seaway, or after its closure. These timing estimates imply that this second bout of speciation was not caused by gene flow impeded by the closure of the seaway, but demand an interpretation involving trans-Arctic gene flow, through a northern route, before being interrupted, most likely this time by a climatic event. Critically, these conclusions are robust to the inclusion of outgroup species (Psettodidae), of outliers species (P. natalensis and P. hawaiiensis), and are consistently found across the four clock models implemented here, including the one (NONE) relying on no other fossil information than the 47.8–42.1 Ma ancestor to this order (Chanet 1997; Friedman 2012)—hereby indicating that our data are highly informative with respect to the timing of the divergence of these sister species (and rejecting our initial hypothesis) following their trans-CAS and trans-Arctic migrations. Trans-Arctic migrations from the Pacific to the Atlantic are not new, and have been described in both invertebrates (e.g., Durham and MacNeil 1967) and vertebrates (e.g., Grant 1986), but never using a relaxed molecular clock relying on as little information as a single fossil and a large multigene data set covering an entire order (Pleuronectiformes). While seminal work was usually based on paleontological (Durham and MacNeil 1967), or on morphological evidence (Reid 1990), subsequent work started using allozymes and other proteins (Grant 1986; Zaslavskaya et al. 1992), followed by genetic evidence, mostly at the population level using mitochondrial DNA (Rawson and Harper 2009). However, most of the recent work that relies on molecular clocks either makes strong assumptions on substitution rates (Rawson and Harper 2009; Liu et al. 2011), uses divergences that are vicariance-motivated but not supported by any fossil information (McCusker and Bentzen 2010), or posit in which direction the trans-Arctic migration took place (Dodson et al. 2007) to infer trans-Arctic migrations. Our work is the first to use sister species to demonstrate the existence and elucidate the timing of both the trans-CAS and the trans-Arctic events, while also inferring in which direction the latter took place, without making too many assumptions. Strikingly, the geological evidence is directly in line with our date estimates. The fossil record suggests that the first aquatic connection between the Pacific and Arctic (and Atlantic) oceans through the Bering Strait occurred ∼5.5–5.4 Ma (Gladenkov et al. 2002) due to a rise in sea levels, linked to tectonic activity (Marincovich and Gladenkov 2001). This would have permitted the trans-Arctic migration of populations ancestral to today’s northern species from one ocean to the other through this ancient “northern passage.” This global warming event, between the late Miocene to early Pleiocene, was then followed by a significant period of cooling during the Pleiocene into the Pleistocene (Zachos et al. 2001), leading to periods of repeated glaciations and a subsequent ice age. These cold events would have resulted in the closure of this ancient “northern passage,” hereby stopping the trans-Arctic gene flow between the two oceans, and leading to the more recent speciation of the northern taxa. For approximately a million years after the first opening of the Bering Strait, water flowed through the strait in a southern direction, from the Arctic to the Pacific ocean, until the formation of the Isthmus of Panama occurred close to the equator (Berta 2012). With the formation of the Isthmus, and hence the closing of the CAS, the ocean currents reversed due to a change in global ocean circulation (Haug and Tiedemann 1998; De Schepper et al. 2015), and have since flowed in a northern direction through the Bering Strait, from the Pacific to the Arctic (Marincovich 2000). This again is absolutely consistent with our results (supplementary fig. S4, Supplementary Material online), with dispersal or migration from the Pacific into the Atlantic in a northern direction through this strait, which is known as the trans-Arctic interchange (Vermeij 1991). Fossil data also show that the Bering land bridge has been exposed and submerged on multiple occasions since the Pleistocene (Gladenkov 2004). These openings and closings of the Bering Strait could have provided a mechanism for divergence and the evolution of sister taxa (Taylor and Dodson 1994; Väinölä 2003). Our results also have implication at the family level of flatfish. Further significant global cooling during the Pleistocene resulted in major glaciation events (Zachos et al. 2008) that could be responsible for creating barriers that isolated populations. All of the remaining sister taxa in our analysis, who have divergence estimates of <2 Ma in our study belong to the family Pleuronectidae (sensuCooper and Chapleau 1998). The Pleuronectidae are the predominant flatfish family found in cold and temperate seas of the northern hemisphere (Norman 1934; Cooper and Chapleau 1998). There are far more Pleuronectidae species in the Pacific Ocean, most of them endemic to the north Pacific Ocean off of North America and Asia in the region extending from the Bering Strait to the gulf of California (Norman 1934). None of the arctic species is restricted solely to arctic waters (Munroe 2015). Munroe also noted that Cooper (in an unpublished manuscript) identified areas of endemism among the current distribution of the Pleuronectidae. It has been shown that during the trans-Arctic interchange, there was a far higher number of species (up to eight times higher) migrating to the Atlantic than to the Pacific (Vermeij 1991). Fossil and phylogenetic evidence suggest the Pacific Ocean as the origin for diversification of the Pleuronectidae (Munroe 2015), and our phylogenetic results are highly congruent with this hypothesis. It is possible that the outliers, the Atlantic and Pacific halibut (H. hippoglossus and H. stenolepis, respectively), diverged during one of the first openings of the Bering Strait. The estimated dates from the molecular clock analysis are ∼5.5 Ma, in accordance with the hypothesized dates for the first aquatic connection (Marincovich and Gladenkov 2001). The remaining sister taxa have a younger age estimate of ∼1–2 Ma, corresponding with global cooling during the Pleistocene and repeated glaciations (Zachos et al. 2001), which likely formed more barriers to genetic flow. In the second pair of outliers, P. hawaiiensis inhabits waters of the Eastern Central Pacific to the Hawaiian Islands, whereas P. natalensis inhabits coastal waters of Eastern Africa (fig. 2). Because of the far reaching range of P. natalensis, and the relatively younger age estimates for divergence (1–2 Ma), these results beg for future research into other vicariant hypotheses, or dispersal routes, as these two members of Poecilopsetta diverged long after the Isthmus had closed. One candidate would be the Indo-Pacific barrier, which formed an almost continuous land bridge between Australia and Asia, and could have hence limited dispersal between the Pacific and the Indian oceans, hereby promoting speciation (Gaither et al. 2011). This barrier formed during the Pleistocene glacial cycles (2.58–0.01 Ma; Voris 2000), which is, again, fully consistent with our date estimates (fig. 3, supplementary fig. S1, Supplementary Material online). On the basis of the most extensive multigene sequence alignment available to date across all flatfish species, we have showed here that the evolutionary dynamics of sister species that are distributed across the two oceans strongly supported the existence of two bouts of speciation: one triggered by the closure of the CAS 12–3 Ma, and a second one due to the closure of an ancient northern passage 5–0.01 Ma. Our work adds to a growing body of evidence pointing not only to the role of geological processes in shaping biotic turnovers (Bacon et al. 2015), but also on the consequences of climate change that can either trigger speciation events (Reid 1990), or promote faunal interchanges whose ecological and economic impacts are difficult to predict (Wisz et al. 2015). Materials and Methods Data Retrieval and Alignment Prior to retrieving sequence data, GenBank was surveyed to identify all the genes belonging to species of Pleuronectiform families (as of August 2016), based on GenBank’s taxonomy browser. DNA sequences for a total of 332 flatfish species (out of over 800 species in the order) were identified and downloaded for five nuclear genes (kiaa1239, myh6, ripk4, rag1, sh3px3), and four mitochondrial genes (12s, 16s, cox1, and cytb). These represented all the taxa having at least one of these six gene sequences in GenBank; see supplementary table S1, Supplementary Material online for the corresponding accession numbers. Diversity was richly sampled as species from 13 of the 14 families in the order Pleuronectiformes were included in our catchment. In line with current consensus in flatfish systematics, the family Psettodidae was chosen as the outgroup to all other taxa (Chapleau 1993). These sequences were aligned using MUSCLE ver. 3.8.31 (Edgar 2004) on a gene-by-gene basis. Each alignment was visually inspected with AliView ver. 1.18 (Larsson 2014), and was manually edited where necessary. In particular, large indels (> 10 bp) were removed prior to all phylogenetic analyses. The 5′ and 3′ ends of sequences were also trimmed. The aligned sequences were then concatenated into a single alignment using a custom R script. Data Preprocessing To gauge the phylogenetic content of our data set, we performed a first series of molecular clock analyses on the concatenated data matrix, with all the nine gene sequences obtained above (12S, 16S, COI, Cyt-b, KIAA1239, MYH6, RIPK4, SH3PX3, and RAG1), and all the 332 taxa representing all families of flatfish. To account for rate heterogeneity along this alignment, partitions were first identified with PartitionFinder 2 (Lanfear et al. 2017), selecting the best-fit model of nucleotide substitution based on Akaike’s Information Criterion (AIC). This optimal partition scheme was then employed in a first Bayesian phylogenetic analysis, conducted with BEAST ver. 1.8.3 (Drummond and Rambaut 2007). This program implements a Markov chain Monte Carlo (MCMC) sampler, which co-estimates both tree topology and divergence times. As determined by PartionFinder, a GTR + I + Γ model was applied to each data partition. These models of evolution across partitions were assumed to be independent (unlinked, in BEAST parlance), whereas both clock model and tree model partitions were shared (linked) by all partitions. An uncorrelated relaxed clock was employed with a lognormal distribution prior on rates, and a Yule speciation prior (Drummond et al. 2006). Because of the paucity of reliably placed fossils on the flatfish tree, a unique calibration point was placed on the most recent common ancestor (MRCA) of the ingroup, as a mean-one exponential prior, with an offset of 40 My reflecting the age of 47.8–42.1 Ma (Chanet 1997; Friedman 2012). To stabilize the estimate, a lognormal ln(0, 1.5) prior with a 40 Ma offset was also placed on the root of the tree (root height). Two separate MCMC samplers were run, each for 10,000,000 generations. Trees were sampled every 5,000 generations, and convergence was checked visually using Tracer ver. 1.6 (available at: http://tree.bio.ed.ac.uk/software/tracer/; last accessed May 24, 2018). Tree log files from each run were combined in LogCombiner, after conservatively removing 10% of each run as burn-in. The resulting maximum a posteriori (MAP) tree was generated with TreeAnnotator (Drummond and Rambaut 2007). As the topology of this resulting MAP tree was unconventional, we suspected the presence of rogue taxa, that is, species evolving either much faster or much slower than the majority, which can contribute negatively to consensus tree support (Aberer et al. 2013). Rogue taxa were identified using the RogueNaRok (Aberer et al. 2013) webserver (http://rnr.h-its.org; last accessed May 24, 2018). The consensus trees from the preliminary analysis using 10,000 iterations were used as the tree set for rogue taxon analysis. A total of 22 rogues were identified and pruned from the analysis, leaving 310 species. Molecular Dating To assess the impact of the closure of the CAS on flatfish evolutionary dynamics, a second set of partitioned relaxed molecular clock analyses was performed (without the rogue sequences). The timing of the closure of the CAS is estimated to have occurred between 12 and 3 Ma (Duque-Caro 1990; Coates et al. 1992; Haug and Tiedemann 1998; O'Dea et al. 2016), and we used this time window as a prior to inform the relaxed molecular clock-based phylogenetic reconstructions. The analyses were performed on the same concatenated data set, with the same single fossil calibration, but we also placed a lognormal prior ln(3, 1.5), that has most of its mass on the 12–3 Ma window, on the MRCA of each pair of sister taxa (the “ALL” model). From the initial BEAST analyses, twelve pairs of taxa were selected based on two criteria: 1) being sister species on that initial tree, 2) with one species distributed in the Pacific and one in the Atlantic ocean (fig. 2). Again, two independent MCMC samplers were run each for 100 million iterations, with samples taken every 5000 step. Because these pairs of sister species show a contrasted geographic distribution, having either a southern (equatorial) or a northern range (fig. 2), two additional sets of analyses were performed. In a first set, calibration priors (ln(3, 1.5), as per above) were placed only sister taxa that had a geographic range in the southern hemisphere (the “SOUTH” model), while in a second set, identical calibration priors were placed only on sister species with a northern range (the “NORTH” model). Finally, a set of analyses was performed using no sister taxa calibrations at all (the “NONE” model). For each analysis, results from the two MCMC runs were combined using LogCombiner after removing an even more conservative burn-in of 50%. The final MAP tree was generated with TreeAnnotator. In an attempt to rank these different models (priors on all sister taxa; only on southern taxa; only on southern taxa; no “CAS” priors), the AICM was computed for each model (Baele et al. 2013). These computations were performed in Tracer for each of the four different calibration models, based on 100 replicates. Species Distribution Data In order to map the distribution of each species of flatfish, we resorted to the Global Biodiversity Information Facility (GBIF: https://www.gbif.org/; last accessed May 24, 2018) database, accessed with the R library rgbif (Chamberlain 2017). Up to 1,000 records were retrieved for each species, and plotted on a geographic map based on GBIF observational records. Where needed, species locations were summarized by mean latitude and longitude (supplementary fig. S4, Supplementary Material online). Approximate attribution to a particular ocean was based on species-specific mean longitudes (Atlantic: between 20° and 120°; Indian: 120° and −100°; Pacific otherwise). Availability of Computer Code and Data Accession numbers, sequence alignments, BEAST models (as xml files), estimated phylogenetic trees (MAP trees from BEAST), and the R scripts used to plot the figures in this study are available from https://github.com/sarisbro, last accessed May 24, 2018. Supplementary Material Supplementary data are available at Molecular Biology and Evolution online. Acknowledgments This work was supported by the University of Ottawa (L.B.), and the Natural Sciences Research Council of Canada (F.C., S.A.B.). We are grateful to Jonathan Dench for discussions and comments, to two anonymous reviewers for very constructive comments, and to Compute Canada and Ontario’s Center for Advanced Computing for providing us access to their resource. This work was completed while S.A.B. was hosted by Prof. Yutaka Watanuki, at the University of Hokkaido in Hakodate, thanks to an Invitational Fellowship from the Japanese Society for the Promotion of Science. References Aberer AJ , Krompass D , Stamatakis A. 2013 . Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice . Syst Biol . 62 1 : 162 – 166 . Google Scholar CrossRef Search ADS PubMed Aris-Brosou S , Rodrigue N. 2012 . The essentials of computational molecular evolution. In evolutionary genomics: statistical and computational methods. In: Anisimova M , editor. Methods in molecular biology (methods and protocols) , Vol. 855 . Totowa (NJ ): Humana Press . Azevedo MFC , Oliveira C , Pardo BG , Martínez P , Foresti F. 2008 . Phylogenetic analysis of the order Pleuronectiformes (Teleostei) based on sequences of 12S and 16S mitochondrial genes . Genet Mol Biol . 31(1 Suppl) : 284 – 292 . Google Scholar CrossRef Search ADS Bacon CD , Silvestro D , Jaramillo C , Smith BT , Chakrabarty P , Antonelli A. 2015 . Biological evidence supports an early and complex emergence of the Isthmus of Panama . Proc Natl Acad Sci U S A . 112 19 : 6110 – 6115 . Google Scholar CrossRef Search ADS PubMed Baele G , Li WLS , Drummond AJ , Suchard MA , Lemey P. 2013 . Accurate model selection of relaxed molecular clocks in Bayesian phylogenetics . Mol Biol Evol . 30 2 : 239 – 243 . Google Scholar CrossRef Search ADS PubMed Berendzen PB , Dimmick WW , McEachran JD. 2002 . Phylogenetic relationships of Pleuronectiformes based on molecular evidence . Copeia 2002 : 642 – 652 . Google Scholar CrossRef Search ADS Berta A. 2012 . Return to the sea: the life and evolutionary times of marine mammals . Berkeley (CA ): University of California Press . Google Scholar CrossRef Search ADS Betancur-R R , Ortí G. 2014 . Molecular evidence for the monophyly of flatfishes (Carangimorpharia: pleuronectiformes) . Mol Phylogenet Evol . 73 : 18 – 22 . Google Scholar CrossRef Search ADS PubMed Betancur-R R , Li C , Munroe TA , Ballesteros JA , Ortí G. 2013 . Addressing gene tree discordance and non-stationarity to resolve a multi-locus phylogeny of the flatfishes (Teleostei: pleuronectiformes) . Syst Biol . 62 5 : 763 – 785 . Google Scholar CrossRef Search ADS PubMed Campbell J , López JA , Satoh TP , Chen WJ , Miya M. 2014 . Mitochondrial genomic investigation of flatfish monophyly . Gene 551 2 : 176 – 182 . Google Scholar CrossRef Search ADS PubMed Campbell MA , Chen WJ , López JA. 2014 . Molecular data do not provide unambiguous support for the monophyly of flatfishes (Pleuronectiformes): a reply to Betancur-R and Ortí . Mol Phylogenet Evol . 75 : 149 – 153 . Google Scholar CrossRef Search ADS PubMed Carnevale G , Bannikov AF , Landini W , Sorbini C. 2006 . Volhynian (early Sarmatian sensu lato) fishes from Tsurevsky, north Caucasus, Russia . J Paleontol . 80 : 684 – 699 . Google Scholar CrossRef Search ADS Chamberlain S. 2017 . rgbif: Interface to the Global ‘Biodiversity’ Information Facility API. R package version 0.9.9. https://CRAN.R-project.org/package=rgbif Chanet B. 1997 . A cladistic reappraisal of the fossil flatfishes record: consequences on the phylogeny of the Pleuronectiformes (Osteichthyes: teleostei) . Ann Des Sci Nat Zool Paris . 18 : 105 – 117 . Chapleau F. 1993 . Pleuronectiform relationships: a cladistic reassessment . Bull Mar Sci . 52 : 516 – 540 . Coates AG , Jackson JBC , Collins LS , Cronin TM , Dowsett J , Bybell LM , Jung P , Obando JA. 1992 . Closure of the Isthmus of Panama: the near-shore marine record of Costa Rica and western Panama . Geol Soc Am Bull . 104 7 : 814 – 828 . Google Scholar CrossRef Search ADS Cooper JA , Chapleau F. 1998 . Monophyly and intrarelationships of the family Pleuronectidae (Pleuronectiformes), with a revised classification . Fish Bull . 96 : 686 – 726 . De Schepper S , Schreck M , Beck KM , Matthiessen J , Fahl K , Mangerud G. 2015 . Early Pliocene onset of modern Nordic Seas circulation related to ocean gateway changes . Nat Commun . 6 : 8659. Google Scholar CrossRef Search ADS PubMed Dodson JJ , Tremblay S , Colombani F , Carscadden JE , Lecomte F. 2007 . Trans-Arctic dispersals and the evolution of a circumpolar marine fish species complex, the capelin (Mallotus villosus) . Mol Ecol . 16 23 : 5030 – 5043 . Google Scholar CrossRef Search ADS PubMed dos Reis M , Donoghue PCJ , Yang Z. 2016 . Bayesian molecular clock dating of species divergences in the genomics era . Nat Rev Genet . 17 2 : 71 – 80 . Google Scholar CrossRef Search ADS PubMed Drummond AJ , Ho SYW , Phillips MJ , Rambaut A. 2006 . Relaxed phylogenetics and dating with confidence . PLoS Biol . 4 5 : e88 – 710 . Google Scholar CrossRef Search ADS PubMed Drummond AJ , Rambaut A. 2007 . BEAST: Bayesian evolutionary analysis by sampling trees . BMC Evol Biol . 7 : 214. Google Scholar CrossRef Search ADS PubMed Duque-Caro H. 1990 . Neogene stratigraphy, paleoceanography and paleobiology in northwest South America and the evolution of the Panama Seaway . Palaeogeogr Palaeoclimatol Palaeoecol . 77 ( 3–4 ): 203 – 234 . Google Scholar CrossRef Search ADS Durham JW , MacNeil FS. 1967 . Cenozoic migrations of marine invertebrates through the Bering Strait region. In: Hopkins DM , editor. The Bering land bridge. Stanford (CA ), Stanford University Press . p. 326 – 349 . Edgar R. 2004 . MUSCLE: multiple sequence alignment with high accuracy and high throughput . Nucleic Acids Res . 32 5 : 1792 – 1797 . Google Scholar CrossRef Search ADS PubMed Friedman M. 2012 . Osteology of †Heteronectes chaneti (Acanthomorpha, Pleuronectiformes), an Eocene stem flatfish, with a discussion of flatfish sister-group relationships . J Vert Paleontol . 32 4 : 735 – 756 . Google Scholar CrossRef Search ADS Gaither MR , Bowen BW , Bordenave TR , Rocha LA , Newman SJ , Gomez JA , van Herwerden L , Craig MT. 2011 . Phylogeography of the reef fish Cephalopholis argus (Epinephelidae) indicates Pleistocene isolation across the indo-pacific barrier with contemporary overlap in the coral triangle . BMC Evol Biol . 11 : 189. Google Scholar CrossRef Search ADS PubMed Gladenkov AY , Oleinik AE , Marincovich L , Barinov KB. 2002 . A refined age for the earliest opening of Bering Strait . Palaeogeogr Palaeoclimatol Palaeoecol . 183 ( 3–4 ): 321 – 328 . Google Scholar CrossRef Search ADS Gladenkov AY. 2004 . Onset of connections between the Pacific and Arctic Oceans through the Bering Strait in the Neogene . Stratigr Geol Correl . 12 : 175 – 187 . Grant WS. 1986 . Biochemical genetic divergence between Atlantic, Clupea harengus, and Pacific, C. pallasi, Herring . Copeia 1986 3 : 714 – 719 . Google Scholar CrossRef Search ADS Harrington RC , Faircloth BC , Eytan RI , Smith WL , Near TJ , Alfaro ME , Friedman M. 2016 . Phylogenomic analysis of carangimorph fishes reveals flatfish asymmetry arose in a blink of the evolutionary eye . BMC Evol Biol . 16 1 : 224. Google Scholar CrossRef Search ADS PubMed Haug GH , Tiedemann R. 1998 . Effect of the formation of the Isthmus of Panama on Atlantic Ocean thermohaline circulation . Nature 393 6686 : 673 – 676 . Google Scholar CrossRef Search ADS Lanfear R , Frandsen PB , Wright AM , Senfeld T , Calcott B . 2017 . PartitionFinder 2: New Methods for Selecting Partitioned Models of Evolution for Molecular and Morphological Phylogenetic Analyses . Mol Biol Evol. 34 3 : 772 – 773 . Google Scholar PubMed Larsson A. 2014 . AliView: a fast and lightweight alignment viewer and editor for large datasets . Bioinformatics 30 22 : 3276 – 3278 . Google Scholar CrossRef Search ADS PubMed Liu JX , Tatarenkov A , Beacham TD , Gorbachev V , Wildes S , Avise JC. 2011 . Effects of Pleistocene climatic fluctuations on the phylogeographic and demographic histories of Pacific herring (Clupea pallasii) . Mol Ecol . 20 18 : 3879 – 3893 . Google Scholar CrossRef Search ADS PubMed Marincovich L , Gladenkov AY. 2001 . New evidence for the age of Bering Strait . Quat Sci Rev . 20 ( 1–3 ): 329 – 335 . Google Scholar CrossRef Search ADS Marincovich L. 2000 . Central American paleogeography controlled Pliocene Arctic Ocean molluscan migrations . Geology 28 6 : 551 – 554 . Google Scholar CrossRef Search ADS McCusker MR , Bentzen P. 2010 . Phylogeography of 3 North Atlantic Wolffish species (Anarhichas spp.) with phylogenetic relationships within the family Anarhichadidae . J Heredity . 101 5 : 591 – 601 . Google Scholar CrossRef Search ADS Munroe TA. 2015 . Systematic diversity of the flatfishes. In: Gibson RN , Nash RDM , Geffen AJ , van der Veer HW , editor. Flatfishes: biology and exploitation. Chichester, West Sussex : John Wiley & Sons . p. 13 – 44 . Norman JR. 1934 . A systematic monograph of the flatfishes (Heterosomata). London : Printed by Order of the British Museum . Google Scholar CrossRef Search ADS O'Dea A , Lessios HA , Coates AG , Eytan RI , Restrepo-Moreno SA , Cione AL , Collins LS , de Queiroz A , Farris DW , Norris RD. 2016 . Formation of the Isthmus of Panama . Sci Adv . 2 8 : e1600883. Google Scholar CrossRef Search ADS PubMed Pardo BG , Machordom A , Foresti F , Porto-Foresti F , Azevedo MFC , Bañon R , Sánchez L , Martínez P. 2005 . Phylogenetic analysis of flatfish (Order Pleuronectiformes) based on mitochondrial 16S rDNA sequences . Sci Mar . 69 4 : 531 – 543 . Google Scholar CrossRef Search ADS Parham JF , Donoghue PCJ , Bell CJ , Calway TD , Head JJ , Holroyd PA , Inoue JG , Irmis RB , Joyce WG , Ksepka DT , et al. . 2012 . Best practices for justifying fossil calibrations . Syst Biol . 61 2 : 346 – 359 . Google Scholar CrossRef Search ADS PubMed Raftery AE , Newton MA , Satagopan JM , Krivitsky PN. 2007 . Estimating the integrated likelihood via posterior simulation using the Harmonic Mean Identity . Bayesian Stat . 8 : 1 – 45 . Rawson PD , Harper FM. 2009 . Colonization of the northwest Atlantic by the blue mussel, Mytilus trossulus postdates the last glacial maximum . Mar Biol . 156 9 : 1857 – 1868 . Google Scholar CrossRef Search ADS Reid DG. 1990 . Trans-arctic migration and speciation induced by climatic change: the biogeography of Littorina (Mollusca: gastropoda) . Bull Mar Sci . 47 : 35 – 49 . Taylor EB , Dodson JJ. 1994 . A molecular analysis biogeography within (genus Osmevus) of relationships and a species complex of Holarctic fish . Mol Ecol . 3 3 : 235 – 248 . Google Scholar CrossRef Search ADS PubMed Väinölä R. 2003 . Repeated trans-Arctic invasions in littoral bivalves: molecular zoogeography of the Macoma balthica complex . Mar Biol . 143 5 : 935 – 946 . Google Scholar CrossRef Search ADS Vandenberghe N , Hilgen FJ , Speijer RP , Ogg JG , Gradstein FM , Hammer O , Hollis CJ , Hooker JJ. 2012 . The Paleogene period. In: Gradstein FM , Ogg JG , Schmitz MD , Ogg GM , editors. The geologic timescale . Amsterdam (The Netherlands ), Elsevier . p. 855 – 921 . Vermeij GJ. 1991 . Anatomy of an invasion: the trans-Arctic interchange . Paleobiology 17 03 : 281 – 307 . Google Scholar CrossRef Search ADS Voris HK. 2000 . Maps of Pleistocene sea levels in Southeast Asia: shorelines, river systems and time durations . J Biogeogr . 27 5 : 1153 – 1167 . Google Scholar CrossRef Search ADS Wisz MS , Broennimann O , Grønkjær P , Møller PR , Olsen SM , Swingedouw D , Hedeholm RB , Nielsen EE , Guisan A , Pellissier L. 2015 . Arctic warming will promote Atlantic-Pacific fish interchange . Nat Climate Change . 5 9 :790–265. Zachos J , Pagani M , Sloan L , Thomas E , Billups K. 2001 . Trends, rhythms, and aberrations in global climate 65 Ma to present . Science 292 5517 : 686 – 693 . Google Scholar CrossRef Search ADS PubMed Zachos JC , Dickens GR , Zeebe RE. 2008 . An early Cenozoic perspective on greenhouse warming and carbon-cycle dynamics . Nature 451 7176 : 279 – 283 . Google Scholar CrossRef Search ADS PubMed Zaslavskaya NI , Sergievsky SO , Tatarenkov AN. 1992 . Allozyme similarity of Atlantic and Pacific species of Littorina (Gastropoda: littorinidae) . J Moll Stud . 58 4 : 377 – 384 . Google Scholar CrossRef Search ADS © The Author(s) 2018. Published by Oxford University Press on behalf of the Society for Molecular Biology and Evolution. All rights reserved. For permissions, please e-mail: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/about_us/legal/notices)

Journal

Molecular Biology and EvolutionOxford University Press

Published: Aug 1, 2018

There are no references for this article.

You’re reading a free preview. Subscribe to read the entire article.


DeepDyve is your
personal research library

It’s your single place to instantly
discover and read the research
that matters to you.

Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month

Explore the DeepDyve Library

Search

Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly

Organize

Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.

Access

Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals.

Your journals are on DeepDyve

Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more.

All the latest content is available, no embargo periods.

See the journals in your area

DeepDyve

Freelancer

DeepDyve

Pro

Price

FREE

$49/month
$360/year

Save searches from
Google Scholar,
PubMed

Create lists to
organize your research

Export lists, citations

Read DeepDyve articles

Abstract access only

Unlimited access to over
18 million full-text articles

Print

20 pages / month

PDF Discount

20% off