The origin and biogeographic diversification of fishes in the family Poeciliidae

The fish subfamily Poeciliinae (sensu Parenti, 1981) is widely distributed across the Western Hemisphere and a dominant component of the fish communities of Central America. Poeciliids have figured prominently in previous studies on the roles of dispersal and vicariance in shaping current geographic distributions. Most recently, Hrbek et al. combined a DNA-based phylogeny of the family with geological models to provide a biogeographic perspective that emphasized the role of both vicariance and dispersal. Here we expand on that effort with a database enlarged in the quantity of sequence represented per species, in the number of species included, and in an enlarged and more balanced representation of the order Cyprinodontiformes. We combine a robust timetree based upon multiple fossil calibrations with enhanced biogeographic analyses that include ancestral area reconstructions to provide a detailed biogeographic history of this clade. Key features of our results are that the family originated in South America, but its major diversification dates to a later colonization of Central America. We also resolve additional colonizations among South, Central and North America and the Caribbean and consider how this reconstruction contributes to our understanding of the mechanisms of dispersal.


Introduction
Species in the fish subfamily Poeciliinae (sensu [1]) are prominent features of fresh-and brackish water environments from Argentina in the south to New Jersey and Illinois in the north and on many islands throughout the Caribbean. The subfamily contains approximately 275 species in 27 genera [2,3]. Members of this subfamily have played important roles in scientific endeavors that range from cancer biology [4] to diverse facets of ecology and evolutionary biology [5]. Their value as research organisms is a function of their short generation times, hardiness in a laboratory environment, and their remarkable biological diversity. One noteworthy feature of the subfamily is variation in maternal provisioning, which ranges from egg laying in one species (Tomeurus gracilis) to livebearing with little or no post-fertilization maternal provisioning to having the functional equivalent of a mammalian placenta [6][7][8]. a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Other forms of variation that are foci of much of the research effort devoted to this subfamily are traits associated with sexual selection and conflict. The males of some species have highly elaborate ornamentation, including swords, enlarged dorsal fins, laterally compressed bodies and sexual dichromatism associated with elaborate courtship and male-male aggression [9]. Mapping these features on the phylogenetic tree of the subfamily shows that they have all evolved multiple times on different branches of the tree. This presence of multiple, independent origins of biologically significant traits makes the unravelling of the evolutionary history of this subfamily a goal of particular interest because that history can help define the circumstances that promoted the evolution of this diversity.
Rosen and Bailey [10] proposed that the Poeciliinae (sensu [1,11] originated in Central America because Central America harbors the greatest biodiversity. An alternative hypothesis is that the subfamily originated in South America [12]. This subfamily figured prominently in the development of vicariance biogeography, where parallel patterns of distribution were noted for some poeciliids, reptiles and plants [13,14]. This parallelism was interpreted as a common response to tectonic activity that successively fractured and rearranged their distributions, creating common patterns of geographic isolation and common pathways for dispersal. The advent of molecular phylogenetics has made it possible to add the specificity gained by integrating well-resolved family trees with estimates of the timing of diversification events and geological events to reconstruct evolutionary history in a more formal fashion. This integration has been applied to individual genera, subgenera or species (e.g. Brachyraphis [15]-Gambusia- [16,17]; Poeciliopsis- [18,19]; Girardinus- [20]; Pseudoxiphophorus- [21]; Poecilia- [22,23]; Limia- [24]; Mollienesia- [25]) . Hrbek et al. [12] presented the most recent comprehensive family-wide analysis.
Hrbek et al. [12] resolved the debate about the origins and history of dispersal of this group into three models, then deployed molecular phylogenetics to discriminate among them. One model was that the subfamily originated during the Cretaceous and had an early distribution that spanned South, Central and North America and hence that there was an earlier land bridge that joined the continents. A second alternative was that the subfamily originated in Central America after the split of North and South America and did not colonize South American until after the formation of the Panamanian land bridge. A third alternative was again that the subfamily originated during the Cretaceous, had an early distribution from Middle to South America, colonized the Greater Antilles from South America during the Eocene, then dispersed in both directions between South and Central America after the closure of the Isthmus of Panama.
Hrbek et al. [12] concluded that the Poeciliinae originated in South America during the Cretaceous. There results suggest a connection between Central and South America during the Cretaceous followed by bidirectional dispersal between the two regions, early colonization of the Greater Antilles, then independent radiations within all three regions. They argue that the colonization of Middle and North America was during the Eocene-Oligocene transition and was most likely from the Greater Antilles. Here we revisit Hrbek et al.'s [12] analysis with the benefit of: 1) the enhancement of the data on which the analysis is based, both in terms of the number and representation of species and in the number of loci per species, 2) a time tree based upon 16 primary fossil calibrations, and 3) the application of analytical techniques that enabled us to reconstruct formal ancestral area reconstructions and hence make inferences about the direction and timing of major dispersion events and subsequent radiations.

Materials and methods
The UCR IACUC committee approved of the research contained in this manuscript. All of this research was conducted under the auspices of an Animal Use Protocol administered by the IACUC at U. C. Riverside. UCR has Animal Welfare Assurance #A3439-01. The associated animal use protocol number is A-20140003.

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][30][31][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.
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] [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]. 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].
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.
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.
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][85][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).  Table 1 2. The full timetree is available in newick format (S1 Table). for understanding the biogeographic history of Poeciliinae (see below) are also provided in Table 1.

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.  [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].  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 sistertaxa 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.
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.
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  Table 3. Probabilities for ancestral ranges based on DEC + J analyses.
Node Number (Fig 2) Probabilities  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.

Discussion
Rosen and Bailey [10] opened their monograph with the observation that "The broad tolerance for brackish and saltwater environments of many of these chiefly freshwater species makes them of particular interest in zoogeographic studies and overwater dispersal. . ." In contrast, Rosen [13] emphasized the importance of vicariance and associated overland dispersal as key determinants of biogeographic distribution. Our results can be interpreted as support for both of these processes. Our results also confirm and extend those of Hrbek et al [12].
As inferred by Hrbek et al. [12], we found that the most recent common ancestor of Poeciliinae was found in South America, dating to the late Paleocene (56.5 Ma, IR) or early Eocene (53.4 Ma, AUTO). With the exception of Xenodexia, all early splitting lineages of the family (Tomeurus, Phalloptychys, Phalloceros, Cnesterodon, Poecilia) are exclusively found or are rooted in South America. We confirmed that the common ancestors of the taxa in the Greater Antilles (Girardinus, Quintana) and the large Central American radiations dispersed from South America approximately 45-40 Mya. These dates are consistent with dispersal over the Aves Land Bridge between South America and the nascent Greater Antilles [12]. Hrbek et al. [12] concluded that the common ancestor of the Central American radiation was derived from the Greater Antilles. Our analysis could not discriminate between the direct colonization of Central America from South America versus a two-step colonization process (first Greater Antilles, then Central America-Node 18 in Fig 2). Our and Hrbek et al'.s [12] analyses also agreed on the approximate timing of the dispersal of the common ancestor of the subgenus Mollienesia from South America to Central America (node 21), the subgenus Limia (save L. heterandria) from South America to the Greater Antilles, (Node 19), the dispersals out of Central America to North America by the genus Gambusia (Node 14) and to South America by Pseudopoecilia festae (Node 17) and Neoheterandria (Node 15). Our analysis revealed additional dispersals among North America, Central America and the Caribbean by members of the genus Gambusia, (Nodes 1-13). While arguments have been made for vicariance as a cause of the nodes between 14 and 25, at least some of those that lie between nodes 1 and 13 must represent dispersal through brackish water or marine barriers, made possible by the ability of many species of Gambusia to tolerate high salt concentrations. Palacios et al. [25] arrived at the same conclusion for the subgenus Mollienesia (Node 21) and the same could be true for Limia (Node 19). Species in both taxa can be found in brackish water and hence are tolerant of high salt concentrations [10,93]. Either dispersal or vicariance represent plausible explanations for the origins of these nodes.
Our analyses confirmed Palacios et al.'s [25] finding of an explicit link between South America and the Caribbean for the Poecilia subgenus Limia. Prior investigators had inferred that Limia was derived from South America [12,24] but their analyses did not include Limia heterandria. This inclusion reveals a monophyletic cluster of species that spans the distance between the presumed origin of the common ancestor of the genus and the Caribbean radiation. Palacios et al. [25] also show that a triad of Caribbean species (Poecilia dominicensis, P. hispaniolae, P. elegans-not included in our analyses) formerly classified with the subgenus Mollienesia are instead more closely allied to Limia. The fact that this clade of three species is basal to the entirety of Limia suggests that they might represent a second colonization of the Greater Antilles from South America. Alternatively, L. heterandria could represent a back dispersal from the Greater Antilles to South America.
Our work also reaffirms Hrbek et al's [12] finding that Heterandria formosa, in the Southeastern United States, is most closely allied with Priapichthys panamensis in Central America This is a distribution that is difficult to reconcile with any form of vicariance. This association is coherent from the perspective of life histories. Both species have superfetation, or the ability to carry multiple broods of young in discretely different stages of development. The distribution of the trait suggests they inherited it from a common ancestor. Heterandra formosa also has matrotrophy, or the functional equivalent of a mammalian placenta. This distribution suggests that matrotrophy evolved in the H. formosa lineage after the evolution of superfetation. This same sequence of events happened three times in the genus Poeciliopsis [8] and once within the genus Poecilia [31,94,95].

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] [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  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.
Supporting information S1 Table. Complete time trees which served as the basis for Fig 1 and Fig 2. (DOCX) S2