Timetree analyses

We performed molecular dating analyses for the subfamily Poeciliinae by integrating the phylogeny presented by Pollux et al. [9] with 16 primary fossil calibration points. We estimated divergence times for Pollux et al.’s [9] RAxML tree and data set, which includes 293 Cyprinodontiformes and outgroups. We used the mcmctree program in PAML 4.4c [26] with (1) autocorrelated rates and hard-bounded constraints and (2) independent rates and hard-bounded constraints. A complete summary of both time trees is presented in newick format as S1 Table. Mcmctree implements Rannala and Yang’s [27] MCMC algorithms. We used six data partitions, each of which was allowed to have its own GTR + Γ model of sequence evolution. We employed a small number of partitions (6) to make mcmctree computationally feasible. Partitions were designed to maximize taxonomic overlap among genes in each partition and were as follows: six loci (TBR1, PTR, ZIC1, PLAG12, RYR3, SRBE) from Li et al. [28], six loci (CCND1, D2, D8, D29, RAB27, T36) with strong taxonomic sampling from Xiphophorus, β-ACTIN sequences for Girardinus spp. that have minimal taxonomic overlap with other studies, seven nuclear genes (GLYT, X-SRC, ENC1, Rag1, RH, SH3PX3, MYH6) with extensive taxonomic representation in Pollux et al. [9], mitochondrial tRNAs and rRNAs, and mitochondrial protein-coding genes (CYTB, ND1, ND2, COI). We set one time unit = 100 million years. Analyses were performed with cleandata = 0. Shape (α) and scale (β) parameters for the gamma prior of the overall rate parameter μ (i.e., rgene_gamma in mcmctree) were 1 and 2.11, respectively. Calculations for the shape and scale parameters of the gamma prior for the rate-drift parameter assumed an age of 157.5 million years for the most recent common ancestor of Otocephala. We employed minimum and maximum constraints at 16 nodes including constraints at ten nodes from Furness et al. [29]. Minimum ages were based on the oldest crown fossils that are assignable to each clade. Maximum ages were based on stratigraphic bounding, phylogenetic bracketing, and phylogenetic uncertainty [29–32]. Chains were run for 100,000 generations after a burn-in of 10,000 generations, and were sampled every 20 generations.

Minimum and maximum constraints on 16 calibrated nodes were as follows:

1. Tetraodontidae. Minimum = 32.25 Ma; maximum = 56.0 Ma. The geologically oldest representative of crown Tetraodontidae is A. winterbottomi, which dates to the early Rupelian (34–32 Ma)[33]. The maximum age for Tetraodontidae allows for the possibility that the fossil genus Eotetraodon is one of the first two outgroups to crown Tetraodontidae [34, 35]. See Furness et al. [29] for a detailed discussion of minimum and maximum calibrations for this node. Betancur-R et al. [36] used similar minimum and maximum ages of 32.0 and 50.0 Ma, respectively.
2. Tetraodontiformes + Lophiiformes. Minimum = 96.9 Ma; maximum = 126.0 Ma. We used minimum and maximum dates from Furness et al. [29] for this node. Tetraodontiformed and Lophiiformes are sister taxa in Miya et al. [37]. The oldest fossils belonging to Lophiiformes (e.g., anglerfish, frogfish) include frogfish (i.e., Eophryne barbutii) and batfish (i.e., Tarkus) from early Eocene (Ypresian) deposits in Italy [38, 39]. The oldest stem Tetraodontiformes fossil is Plectocretacicus clarae from the Cenomanian of Lebanon [40]. P. clarae has a minimum age of 96.9 Ma [41]. The maximum age is based on stratigraphic bounding and is 126.0 Ma.
3. Moronidae + Lutjanidae. Minimum = 47 Ma; maximum = 84.2 Ma. We used minimum and maximum dates from Furness et al. [29]. The minimum age is based on fossil moronids from middle Eocene Messel deposits that are ~ 47 Ma [42]. Skeletal remains of moronids are also known from Lutetian deposits of the Green River Formation, Wyoming (e.g. Priscacara serrate, [43]). Lutjanid skeletal remains are known from Lutetian deposits of the Monte Bolca locality (e.g. Ottaviania leptacanthus, [44]). The maximum age for this node (83.5 +/- 0.7 = 84.2 Ma) is based on otoliths from the Late Campanian with putative moronid affinities [45]. Alfaro et al. [46] and Santini et al. [47] both used a minimum of 74 Ma for crown Moronidae base on fish otoliths (Morone sp.) from the Late Cretaceous Coffee Sand of Mississippi [48], but in a later publication Nolf [49] refers to the oldest putative moronid fossil from the Cretaceous as? Monoridae indet. We therefore used a younger minimum age as in Furness et al. [29].
4. Zoarcoidea + Gasterosteiformes. Minimum = 13.0 Ma; maximum = 84.2 Ma. Minimum and maximum ages are from Furness et al. [29]. The minimum age of 13.0 Ma is established by Gasterosteus aculeatus, which is known from the Monterey Formation, California (13–13.3 My; Serravallian; [50]) and is the oldest known non-indostomid gasterosteiform. The maximum age allows for the possibility that the Campanian Gasterorhamphosus has stem group affinities with Zoarcoidea + Gasterosteiformes.
5. Poeciliinae + Anablepidae. Minimum = 39.9 Ma; maximum = 71.2 Ma. Minimum and maximum ages are from Furness et al. [29]. The minimum age is based on undescribed poeciliid fossils that are known from the Lumbrera Formation in Argentina [51]. As noted by Furness et al. [29], geological data suggest an upper Ypresian-lower Lutetian age for the Faja Verde of the Lumbrera Formation. However, the minimum age in only 39.9 Ma based on U/Pb zircon dating of a tuff bed that overlies this formation [51]. The maximum is established by cf. Cyprinodontiformes fossils from the El Molino Formation in Bolivia of Maastrichtian to Paleocene age [52] that have uncertain relationships to living taxa. Santini et al. [47] employed a similar calibration (minimum = 55 Ma, maximum = 99 Ma) for Poeciliidae to Fundulidae.
6. Otocephala = Ostarioclupeomorpha (Ostariophysi + Clupeomorpha). Minimum = 149.85; maximum = 165.2 Ma. The oldest crown Ostarioclupeomorpha is Tischlingerichthys viohli from the upper Tithonian of Muhlheim, Bavaria, Germany [53, 54]. We followed Benton et al. [41] and Furness et al. [29] and employed a minimum age of 149.85 Ma based on T. viohli. The maximum age is based on stratigraphic bounding as in Furness et al [29]. Alfaro et al. [46] employed minimum and maximum ages of 149 and 152 Ma for this node.
7. Clupeocephala (= Otocephala + Euteleostei). Minimum = 149.85 Ma; maximum = 165.2 Ma. The oldest crown otocephalan is Tischlingerichthys viohli from the upper Tithonian of Muhlheim, Bavaria, Germany [53, 54]. The oldest euteleosts are Leptolepides and Orthogonikleithrus from the same deposits [41]. The maximum age is based on stratigraphic bounding. Our minimum and maximum agree with Benton et al. [55]. Subsequently, Benton et al. [55] suggested a much older maximum (235 Ma) for this clade while acknowledging that a maximum bound for crown Clupeocephala can only be defined arbitrarily.
8. Clupeomorpha. Minimum = 82.8 Ma; maximum = 149.5 Ma. Clupeomorph taxa in our phylogeny include one clupeid and two pristigasteroids. The minimum is based on Gasteroclupea branisai, which is found in Santonian deposits of Bolivia (Chaunaca Formation; [56]) and has been treated as Clupeidae [57], Pristigasteroidea [57], or Pristigasteridae [56]. The oldest stem clupeomorphs are Neocomian (i.e.. Scutatuspinosus itapagipensis, [58, 59]. The base of the Neocomian is 149.5 Ma.
9. Siluriformes + Characiformes. Minimum = 92.8 Ma; maximum = 126.0 Ma.Bertini et al. [60] described Siluriformes indet. from skeletal elements from the lower portion of Adamantina Formation (Locality 99, Brazil). The age of this formation is Turonian-Santonian [61]. Cione and Prasad [62] suggest the oldest catfish records (cf. Ariidae) are Campanian (see [52]). The oldest fossils assignable to the characiform are isolated teeth from the Wadi Milk Formation (Cenomanian), Sudan [63]. We based the minimum age for this node on these Cenomanian fossils (min = 93.6 +/- 0.8 = 92.8 Ma). The placement of the fossil species Santanichthys diasii is controversial. Historically, Santanichthys has been treated as Clupeocephala incertae sedis, Clupeomorpha, or Clupavidae (reviewed in Malabarba and Malabarba, [63]. Most recently, Filleul and Maisey [64] place Santanichthys in Characiformes (Otophysi) while Malabarba and Malabarba [63] and Diogo et al. [65] argue that Santanichthys is a stem otophysian. Santini et al. [47] treated Santanichthys diasii as the oldest stem characiform. We treat Santanichthys diasii as Otophysi incertae sedis. Stratigraphic bounding suggests a maximum age of 126.0 Ma.
10. Ostariophysi. Minimum = 133.9 Ma; maximum = 155.6 Ma. The oldest undisputed crown ostariophysian is the gonorynchiform Rubiesichthys gregalis (Valanginian-Berriasian, Serra del Montsec, Lerida Province, Spain) [66]. The minimum age for Ostariophysi is based on the upper boundary of the Valanginian (133.9 Ma). The maximum age based on phylogenetic bracketing/phylogenetic uncertainty is based on Tischlingerichthys viohli (Tithonian, Solnhofen Limestone Formation, Mühlhei, Bavaria, Germany; [53, 54]), which has been traditionally treated as a stem characiform [54], but more recently [65] has been hypothesized as a stem ostariophysian. In this study we treat Tischlingerichthys viohli as Ostariophysi incertae sedis. The maximum age based on stratigraphic bounding is two chronologic units below the Berriasian, i.e., 155.6 Ma.
11. Esociformes + Salmoniformes. Minimum = 70.0 Ma; maximum = 126.0 Ma. We used minimum and maximum dates from Furness et al. [29] for Esociformes + Salmoniformes. Estesesox foxi is the oldest esociform fossil and is from the Aquilan (early Campanian), Milk River Formation, Alberta [67, 68]. We use the end of the Campanian (70.6 + 0.6 = 70.0 Ma) to establish a minimum age for this clade. The maximum age is based on Helgolandichthys schmidi (early Aptian, Tock, Helgolang, Germany), which was originally described as a salmoniform [69] although we are unaware of any phylogenetic studies that have included this taxon.
12. Macrouridae to Gadidae (in Gadiformes). Minimum = 58.5 Ma; maximum = 131.5 Ma. Minimum and maximum dates are from Furness et al. [29]. The oldest crown taxa in this clade are macrourids and gadids that are known from otoliths of Selandian age [70, 71]. The definitive sister-taxon to Gadiformes is unclear, but the candidate lineages include other Neoteleostei, among which are Barremian Aulopiformes [72] that establish a maximum age of 131.5 Ma.
13. Cichlidae. Minimum = 40.2 Ma; maximum = 100.5 Ma. Our minimum and maximum dates are from Furness et al.[29]. The oldest crown cichlids are New World fossils (e.g. Proterocara argentina, Gymnogeophagus eocenicus) from the Lumbrera Formation of Argentina and are Ypresian-Lower Lutetian (Eocene) in age [73]. Phylogenetic analyses (e.g., [74]) recover Proterocara argentina and Gymnogeophagus eocenicus as crown members of the South American Cichlinae and establish a minimum age for Cichlidae [29]. The maximum age allows for phylogenetic uncertainty in the placement of the Plectocretacicus clarae from the Cenomanian, i.e., 99.6 +/- 0.9 = 100.5 Ma.
14. Neoteleostei. Minimum = 124 Ma; maximum = 154.8 Ma. Minimum and maximum dates are from Furness et al.[29]. The oldest neoteleostians (e.g. Atolvorator longipectoralis) belong to Aulopiformes and are Barremian in age [72]. These fossils establish a minimum age of 124 Ma for Neoteleostei. Pachyrhizodontoidea are stem euteleosts [75] or clupeocephalans incertae sedis [76] and may belong to the second outgroup to Neoteleostei. The oldest pachyrhizodontoid is an indeterminant Tithonian species from Banos del Flaco Formation, Termas del Flaco, Chile [77], which establishes a maximum age of 154.8 Ma.
15. Goodeidae. Minimum = 5.332 Ma; maximum = 23.03 Ma. Tapatia occidentalis (Jalisco, Mexico, late Miocene; [78]) and Empetrichthys erdidi (Posey Canyon shale, late Miocene, southern California;[79]) are the oldest known goodeids [80].[81, 82]. Given the uncertain position of Tapatia within Goodeidae we used this taxon as a minimum for crown Goodeidae. The maximum is based on stratigraphic bounding and is 23.03 Ma.
16. Empetrichthys to Crenichthys. Minimum = 5.332 Ma; maximum = 23.03 Ma. Empetrichthys erdisi is known from the Posey Canyon Shale. This member of the Peace Valley Formation was originally considered to be Pliocene in age (e.g., [79]), but more recent evidence [83] suggests a Miocene age. We used the top of the Miocene as the minimum age for the split between Empetrichthys and Crenichthys, and the base of the Miocene as the maximum based on stratigraphic bounding.

Biogeographical reconstructions

Ancestral areas were reconstructed with minimum area change (MAC) parsimony [32] and BioGeoBears [84–86] for 252 cyprinodontiform taxa from Pollux et al.’s [9] data matrix (area assignments for each taxa are summarized in S2 Table). We performed the BioGeoBears analysis with the timetree that was reconstructed with mcmctree using autocorrelated rates and hard-bounded constraints. MAC parsimony analyses were performed with Mesquite [87]. Taxa were coded for their occurrence in one or more of the following six areas: North America (excluding Mexico), Central America + Mexico, South America, West Indies, Africa, Europe + Asia (Supplementary Information). We allowed a maximum of three areas per node, which exceeds the range of extant taxa, all of which are confined to one or two areas. Analyses with BioGeoBears were performed in R [85] and employed the dispersal-extinction-cladogenesis (DEC) model of Ree and Smith [88] and the DEC + J model of Matzke [84, 86] that combines DEC with an allowance for jump dispersal. DEC and DEC + J analyses were performed with a maximum size of three areas. DEC and DEC + J models were compared using a likelihood ratio test.

Molecular dating analyses

Full timetrees (293 taxa from Pollux et al., [9]) based on autocorrelated and independent rates, including 95% credibility intervals, are provided in Supplementary Information (S1 Table). Fig 1 shows the full tree and the location of the 16 calibrated nodes while Fig 2 expands on the portion of the tree in the neighborhood of the Poeciliidae. Divergence times for select nodes are provided in Table 1 based on autocorrelated rates (AUTO) and independent rates (IR). Analyses based on AUTO and IR provide similar estimates for the most recent common ancestor of Clupeocephala (AUTO = 163.8 Ma, IR = 163.1 Ma) and other deep nodes within this clade (Table 1). The most recent common ancestor of Cyprinodontiformes was estimated as 78.0 Ma (AUTO) to 90.7 Ma (IR). Estimates for the most recent common ancestor of Poeciliinae ranged from 53.4 Ma (AUTO) to 56.5 Ma (Table 1). Dates for other nodes that are relevant for understanding the biogeographic history of Poeciliinae (see below) are also provided in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Timetree for 293 Cyprinodontiformes and outgroup taxa.

Divergence times were estimated for Pollux et al.’s [9]RAxML tree using the mcmctree program in PAML 4.4c [26] with autocorrelated rates and hard-bounded constraints. Numbered nodes correspond to the sixteen fossil calibrations described in Materials and Methods. The estimated age for Tetraodontiformes to Lophiiformes (85.0 Ma) is younger than the minimum calibration age (96.9 Ma) for this node. Family names (right) are indicated for species in the order Cyprinodontiformes. The clades Poeciliinae and Anablepidae (top) have been collapsed; an expanded view of these clades can be found in Fig 2. The full timetree is available in newick format (S1 Table).

https://doi.org/10.1371/journal.pone.0172546.g001

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Timetree and ancestral area reconstructions for Poeciliinae + Anablepidae.

Divergence times were estimated for Pollux et al.’s [9] RAxML tree using the mcmctree program in PAML 4.4c [26] with autocorrelated rates and hard-bounded constraints. Ancestral areas were reconstructed with minimum area change parsimony [32].

https://doi.org/10.1371/journal.pone.0172546.g002

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Divergence times (Ma) taken from this and previous studies.

Details on the calibrations used in each study can be found in the main text.

https://doi.org/10.1371/journal.pone.0172546.t001

Biogeographic analyses

Our tree topology replicates all of the 9 major clades within the Poeciliinae defined by Hrbek et al. [12], but adds details not present in the earlier analyses largely because of our broader representation of the order Cyprinodontiformes and additional species of Poeciliinae. One such detail is that our tree provides a robust definition of the sister taxa to Parenti’s Poeciliinae. Our phylogeny differs from Parenti’s [1, 11] by defining the Anablepidae as the sister taxa to Poeciliinae. Both of these taxa are exclusively New World clades. Parenti instead identified the subfamilies Fluviphylacinae (genus Fluviphylax–represented here by two species) and Aplocheilichtyinae (genus Aplocheilichthys, represented here by two species) as sister to the Poeciliinae. It was this sister relationship that motivated her to reclassify Rosen and Bailey’s [10] Poecillidae as one of three subfamilies of her newly defined Poeciliidae. In our tree, Aplocheilichthys is most closely allied with Valencidae. Both of these taxa are from the Old World and both are exclusively egg layers. We found that the Fluviphylacinae are sister to (Poeciliinae + Anablepidae). The Fluviphylacinae are composed of egg layers found in the Amazon basin. The closer affinity of Poeciliinae with Anablepidae and Fluviphylacinae, to the exclusion of Aplocheilichthyinae means that the radiation of these taxa likely took place on South America after it became isolated as an island continent. Our time tree analysis supports this conclusion, since the estimate for the node that represents the common ancestor of the Anablepidae, Poeciliinae and Fluviphylacinae dates to between 60 and 65 million years (Fig 2).

We applied additional scrutiny to our conclusion that Parenti’s [1] Poeciliinae and Anablepidae are sister taxa because this topology represents a departure from the conventional definition of the Poeciliidae [1, 11]. Statistical tests (Kishino-Hasegawa, Templeton, winning-sites) provide additional support for the hypothesis that Parenti’s [1, 11] Poeciliinae and Anablepidae are sister taxa and form a monophyletic clade to the exclusion of Aplocheilichthys and Fluviphylax (Table 2). All three topological resolutions of Parenti’s [1, 11] Poeciliidae required significantly more steps than Anablepidae + Poeciliinae, including tests that excluded the X-SRC gene that provided the basis for Meyer and Lydeard’s [89] original proposal of a sister-taxa relationship between the Poeciliinae and Anablepidae. Indeed, a tree patterned on Parenti’s hypothesis required 51 (with X-SRC) or 48 (without X-SRC) additional steps to constrain the monophyly of Parenti’s [1, 11] Poeciliidae relative to Meyer and Lydeard’s [89] hypothesis. We thus have strong confidence in the conclusion that the Anablepidae are sister to the Poeciliinae.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Summary of statistical tests regarding the sister group to Poeciliinae.

https://doi.org/10.1371/journal.pone.0172546.t002

The results of biogeographic analyses are reported in Fig 2 (MAC parsimony) and Table 3 (DEC +J). We also performed an analysis with DEC, but this model was rejected by a likelihood ratio test (lnL DEC = -203.3, lnL DEC + J = -162.8, p = 2.3x10^-19) and DEC results are not shown. Overall, MAC parsimony and DEC + J analyses yielded very similar results (Fig 2, Table 3), which strengthens our interpretation of the results reported below. Importantly, MAC parsimony and DEC + J both allow for dispersal to a new and unique area on a given branch. Biogeographic analyses strongly support a South American origin for the common ancestor of Poeciliinae and Anablepidae. Reconstructions of ancestral areas also confirm that Poeciliinae originated in the New World, most likely South America, and reveal details of the subsequent dispersal and radiation of this subfamily. Dispersal events are summarized in Table 4 and detailed in S3 Table. The alignment between our estimated divergence dates and those of Hrbek et al.[12] appear in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 3. Probabilities for ancestral ranges based on DEC + J analyses.

https://doi.org/10.1371/journal.pone.0172546.t003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 4. Dispersal events in Poeciliinae based on MAC parsimony and DEC + J^*.

https://doi.org/10.1371/journal.pone.0172546.t004

There were two dispersal events from South America to Central America. The first was the dispersal of the ancestor of Xenodexia ctenolepis (node 25 in Fig 2), dating to ≤53.4 Ma with AUTO and ≤ 56.4 Ma with IR. The second, more recent dispersal was the common ancestor of Limia + Pamphorichthys + Mollienesia (node 21) to the most recent common ancestor of subgenus Mollienesia (25.5–19.0 Ma, AUTO/20.1–14.3, IR).

There are two dispersal events associated with Node 18 that were ultimately derived from South America (node 24), but the path from South America through node 18 to its immediate descendant nodes cannot be resolved. Both events resulted in substantial radiations after dispersal. One dispersal event (40.8–37.2 Ma, AUTO/41.5–37.5 Ma, IR) leads to the most recent common ancestor of the Central American clade, which has radiated into the majority of extant species. The second dispersal event (40.8–25.9 Ma, AUTO/41.5–25.5 Ma, IR) gave rise to the common ancestor of the Girardinus/Quintana clade of the Greater Antilles. Our lack of resolution at node 18 means we cannot resolve whether the common ancestor of Girardinus/Quintana was derived from South or Central America. It also means we cannot resolve whether the common ancestor of the Central American radiation was derived from South America or the Greater Antilles. The timing of these colonization events is close to that reported by Hrbek et al. [12].

There was one dispersal from South America to the Caribbean which resulted in the radiation of the Poecilia subgenus Limia (Node 19, 22.8–16.2 Ma, AUTO/ 16.7–11.4 Ma, IR). One interesting feature of our tree pertains to the status of Poecilia (Limia) heterandria, a species whose subgeneric classification has varied over time. It has formerly been allied either with subgenus Poecilia or subgenus Pamphorichthys. Most recently, Poeser [91] classified it as its own genus (Pseudolimia). Our analysis shows that it is the basal lineage to the subgenus Limia, but is the only member of the subgenus found in South America.

There were nine later dispersal events out of Central America—two to South America, four to North America and three to the West Indies. The two dispersals from Central to South America resulted in Neoheterandria elegans (node 15, Ma ≤ 21.4 Ma, AUTO/≤ 21.7 Ma, IR) from Colombia and Pseudopoecilia festae (node 17, ≤ 24.4 Ma, AUTO/≤ 22.5 Ma, IR) from Colombia and Ecuador. These are the only representatives of the subfamily that have colonized the Pacific coast of South America. Ho et al. [23] report an additional dispersal event from Central to South American by species of Poecilia. This event is not present in our results because our data set lacked all Poecilia species derived from the ancestor that migrated from Central to South America.

The earliest colonization of North America is represented by Heterandria formosa, which is most closely related to Priapichthys panamensis (node 16) and arrived in North America ≤ 25.7 Ma (AUTO) or ≤ 22.8 Ma (IR). This inference, first made by Hrbek et al. [12], is of particular interest because Rosen and Baily [10] originally allied it with what is now classified as Pseudoxiphophorus, a genus represented by 12 species in Guatemala and Mexico [92].

The most recent dispersal events in our tree (nodes 1–14, <1.6 Ma– 11 Ma) all involve the genus Gambusia. These include three independent dispersals from Central America to North America and three from Central America to the Caribbean. The only dispersal events out of North America were three to Central America, again all by species in the genus Gambusia. There are also species in this genus that have ranges that span North America and the Caribbean.

Timetree comparisons with other studies

Some of the studies of individual genera in the subfamily Poeciliinae yield time tree estimates for some of the same nodes contained in our tree. These estimates sometimes differ by a wide margin. Such variation invites a consideration of the means by which these estimates are arrived at and their efficacy. Here we consider some of this variation and some of the possible causes.

Our taxon sampling differs substantially from other time tree studies of bony fishes and necessarily includes its own unique set of calibrated nodes. For example, Betancur-R et al. [36] included a diverse assemblage of 1410 bony fishes plus four tetrapod and two chondrichthyan outgroups in their analyses whereas our analysis was limited to 293 taxa belonging to Clupeocephala. Betancur-R et al.'s [36] data set included 13 species in Cyprinodontiformes whereas our data set included 242 Cyprinodontiformes. Our estimated timetree dates also show both similarities and differences with Betancur-R et al. [36] and other studies (e.g., [12, 15, 20–25, 47, 96] (Table 1).

Divergence dates for several nodes that are deep in our tree, including Clupeocephala, Otocephala, Ostariophysi, and Euteleostii, are consistently younger than dates from Betancur-R et al. [36], but are comparable to dates reported by Alfaro et al. [46] and Santini et al. [47] (Table 1). We suspect that our younger dates result from constraining the root of the tree to have a maximum date of 165.2 Ma following Benton et al. [41]. As noted by Benton et al. [55], a maximum bound for Clupeocephala is highly arbitrary and a much older maximum is also possible. Betancur-R et al. [36] did not calibrate Clupeocephala, but calibrated several deeper nodes in their tree including Gnathostomata (519–426 Ma), Osteichthyes (438–418), Actinopterygii (423–398 Ma), and Actinopterii (415–375 Ma). This calibration strategy allowed Betancur-R et al.’s [36] estimated age for Clupeocephala to be significantly older than our estimated age.

Given our younger dates for deeper nodes, it is perhaps surprising that we recovered older dates within Cyprinodontiformes than Betancur-R et al. [36] (Table 1). Our estimates for the most recent common ancestor of Cyprinodontiformes are 78.0 and 90.7 Ma with autocorrelated and independent rates, respectively. Betancur-R et al.'s [36] date for this node is slightly younger at 69.3 Ma. Betancur-R et al.'s [36] dates for nodes within Cyprinodontiformes exhibit even more disagreement with our estimated dates. For example, Betancur-R et al.[36] estimated a divergence date of 21.8 Ma for Cyprinodontidae (Floridichthys to Jordanella) whereas our estimated dates for this node are ~2.5X larger at 44.1 to 49.6 Ma. Similarly, Betancur-R et al. [36] estimated that the Poeciliinae genera Heterandria and Poecilia last shared a common ancestor 17.4 Ma whereas our analyses recovered divergence time estimates of 40.8 Ma (AUTO) and 41.5 Ma (IR). Unlike Betancur-R et al. [36], who did not include any calibrations inside of Cyprinodontiformes, we included three calibrations (Fig 1, Node 5—Poeciliinae to Anablepidae; Node 15 –the Goodeidae; Node 16 –Empetrichthys to Crenichthys) in this order. The inclusion of these cyprinodontiform calibrations in our analysis may explain, at least in part, the different dates that were recovered in our analyses. Our data set also included a greater proportion of mitochondrial genes than Betancur-R et al [36], which adds resolution to the characterization of more recent nodes. Finally, as noted above, our taxon sampling for Cyprinodontiformes is much more dense than Betancur-R et al.'s [36] taxon sampling within this order.

Our dates for Poeciliinae and its subclades based on a fossil-calibrated tree are in good agreement with Hrbek et al.'s [12] dates for several overlapping nodes (Table 4). The similarities between the two studies is striking because there is little overlap in the DNA sequence data on which the trees were based and no overlap in way the tree was calibrated. Whereas we employed fossil calibrations to date our tree, Hrbek et al.[12] employed a secondary calibration from Reznick et al. [8], who in turn calibrated their tree for Poeciliopsis based on a vicariant event. What we have in common is extensive, subfamily-wide coverage.

Other studies have obtained both older and younger dates for divergences within Poeciliinae. Alda et al. [22] calibrated their tree based on alternate dates (25–20 Ma, 17–14 Ma) for the separation of Hispaniola versus Cuba (Wind Passage), which they tied to the cladogenetic separation of Limia melanonotata and L. vittata. Based on these alternate calibrations, Alda et al.[22] recovered dates that in some cases are more than twice as old as our estimated ages (Table 4). For example, our point estimates for Mollienesia to Limia are 20.1 and 25.5 Ma with independent and autocorrelated rates, respectively. By contrast, Alda et al. [22] obtained a date of 56 Ma for this split based on their 25–20 Ma calibration (date with 17–14 Ma calibration not reported). Doadrio et al [20] and Weaver et al. [24] used the 17–14 Ma date for the separation of Hispaniola and Cuba for their calibration of the cladogenic separation of Limia vittata from the remainder of the genus. Doadrio et al. [20] arrived at a point estimate of 62 Ma for the cladogenic separation of the common ancestor of Girardinus/Quintana from the Central American clade (node 18). We instead arrived at a point estimate of 40.8/41.5 Ma (autocorrelated/independent). Weaver et al. [23, 24] arrived at a point estimate of 32.9 Ma for the cladogenic separation of the Poecilia subgenus Limia from the Poecilia subgenus Pamphorichthys. We arrived at point estimates of 24.1/19.1 Ma (autocorrelated/independent) for this same event.

At the other extreme, Ho et al.[23] obtained much younger dates for various nodes within Poeciliinae and criticized previous analyses that employed biogeographic calibrations. More specifically, Ho et al.[23] suggest that inferences based on biogeographic calibrations always incur circular reasoning and should be avoided. However, this assessment is only logical if the estimated dates for a node used as a calibration point is subsequently employed to infer the age of the same vicariant event. Similarly, estimated ages for fossil-calibrated nodes have the potential to become circular if the estimated ages from a molecular dating analysis are championed as supporting the same divergence time as suggested by said fossils. Perhaps more importantly, Ho et al.'s [23] calibrated their timetree analysis with a single secondary calibration taken directly from Betancur-R et al. [36] for the separation of Poecilia and Heterandria (14.3 Ma). As noted above, Betancur-R et al.'s [36] analysis included far fewer Cyprinodontiformes than were included in our analysis (13 versus 242) and lacked calibrations within this clade.

Palacios et al. [25] took a different approach to developing a time tree, which was to use the “universal MtDNA mutation rate” to calibrate their tree for the Poecilia subgenus Molliensia. Their calibrations yielded consistently younger divergence dates than ours. For example, their estimated date of origin of the genus Poecilia was 16.4 Ma, as compared to our point estimates of 30.5/26.4 Ma (autocorrelated/independent).

Given widespread discordance among published studies for divergence times within Poeciliinae, we suggest that the true divergence dates remain to be determined. Our contribution is of value for what it reveals about the patterns of colonization and divergence, but what it can tell us about vicariance versus dispersal is limited by the uncertainty in the estimation of the timing of events. Future fossil discoveries loom as critical for gaining increased confidence in the timing of the Poeciliinae radiation. We also suggest that the time is ripe to perform simulation studies that will illuminate how the choice of calibration points and the nature of the data that serve as a basis for defining a phylogeny can jointly influence the values generated for a time tree.