Evolutionary history of the fish genus Astyanax Baird & Girard (1854) (Actinopterygii, Characidae) in Mesoamerica reveals multiple morphological homoplasies

Background Mesoamerica is one of the world's most complex biogeographical regions, mostly due to its complex geological history. This complexity has led to interesting biogeographical processes that have resulted in the current diversity and distribution of fauna in the region. The fish genus Astyanax represents a useful model to assess biogeographical hypotheses due to it being one of the most diverse and widely distributed freshwater fish species in the New World. We used mitochondrial and nuclear DNA to evaluate phylogenetic relationships within the genus in Mesoamerica, and to develop historical biogeographical hypotheses to explain its current distribution. Results Analysis of the entire mitochondrial cytochrome b (Cytb) gene in 208 individuals from 147 localities and of a subset of individuals for three mitochondrial genes (Cytb, 16 S, and COI) and a single nuclear gene (RAG1) yielded similar topologies, recovering six major groups with significant phylogeographic structure. Populations from North America and Upper Central America formed a monophyletic group, while Middle Central America showed evidence of rapid radiation with incompletely resolved relationships. Lower Central America lineages showed a fragmented structure, with geographically restricted taxa showing high levels of molecular divergence. All Bramocharax samples grouped with their sympatric Astyanax lineages (in some cases even with allopatric Astyanax populations), with less than 1% divergence between them. These results suggest a homoplasic nature to the trophic specializations associated with Bramocharax ecomorphs, which seem to have arisen independently in different Astyanax lineages. We observed higher taxonomic diversity compared to previous phylogenetic studies of the Astyanax genus. Colonization of Mesoamerica by Astyanax before the final closure of the Isthmus of Panama (3.3 Mya) explains the deep level of divergence detected in Lower Central America. The colonization of Upper Mesoamerica apparently occurred by two independent routes, with lineage turnover over a large part of the region. Conclusion Our results support multiple, independent origins of morphological traits in Astyanax, whereby the morphotype associated with Bramocharax represents a recurrent trophic adaptation. Molecular clock estimates indicate that Astyanax was present in Mesoamerica during the Miocene (~8 Mya), which implies the existence of an incipient land-bridge connecting South America and Central America before the final closure of the Isthmus of Panama (~3.3 Mya).


Background
Mesoamerica is one of the most complex biogeographical areas in the world [1][2][3][4][5]. This complexity reflects the confluence of Neotropical and Nearctic biotas and a long history of geological activity, stretching from the Miocene to the present, during which movements of the Cocos, North American, Pacific and Caribbean Plates [6,7] created barriers and land-bridges that have affected the distribution of freshwater fishes [8][9][10][11][12][13]. For example, the Pliocene (~3. 3 Mya) closure of the Panama Strait has been postulated to be one of the most important causes of faunal interchange between Neartic and Neotropical regions [14]. Climatic changes have also been invoked to explain the distribution of Mesoamerican fish fauna [15]. Distinguishing between climatic and geological effects requires information on phylogeny and species boundaries in a diversity of taxa [8,9,[16][17][18][19].
The absence of primary freshwater genera (e.g., Hypopomus, Pimelodella, Rhamdia and Roeboides) from Mesoamerica and the Antillean islands argues against an ancient colonization route through a proto-Antillean arc. Instead, it supports a colonization route through an incipient land-bridge formed during the gradual uplifting of the Panama Isthmus [8,12], over a time span of 3-20 Mya, combined with changes in sea level [25]. This is supported by molecular studies of arthropods [26], amphibians [27], and marine geminate species pairs on either side of the Panamanian Isthmus [28]. Support for a more recent colonization of Mesoamerica by primary freshwater fish through the Panama Strait comes from phylogeographical studies of Characids (e.g., Brycon, Bryconamericus, Eretmobrycon, and Cyphocharax). These studies indicate multiple waves of rapid expansion from South America during the Pliocene ~3. 3 Mya. [15].
The genus Astyanax provides an ideal model to investigate the relative importance of vicariance and dispersal on biogeographical patterns. This is partly because it is widely distributed across the region [29], and because its dispersal is confined to freshwater routes and dependent, therefore, on the formation of land-bridges.
Characiforms are generally assumed to have a Gondwanan (South American) origin [30][31][32], as supported by the fossil record [33], so the presence of Characidae in Northern America is viewed as a consequence of dispersal.
Astyanax comprises more than 107 recognized species and is, together with Hyphesobrycon (105 species), the largest and most diverse characiform genus [34,35]. Moreover, Astyanax has the widest distribution of American characids, being found from the Nearctic (Colorado River in Texas and New Mexico) to the Neotropics (Negro River in Patagonia) [4].
Previous phylogenetic studies [36,37] of the biogeography of Astyanax used a small number of samples from Mexico, Belize and Guatemala, and did not find geographical congruence for some of the groups recovered (i.e. Yucatan and Belizean populations were not the most closely related despite their geographical proximity). Furthermore, conspecific cave and nearest surface populations formed two separate lineages, in agreement with an earlier study of the genus [38]. This was attributed to at least two separate colonizations of Mesoamerica from South America during the Pleistocene. Estimated colonization times based on the cytochrome b gene were 1.8 and 4.5 Mya (3.1 Mya), with an estimated divergence rate of 1.5% per pairwise comparison per million years, which coincides with the closure of the Panama Strait (3.3 Mya) [37]. However, incomplete sampling (only few samples were included from upper Central America and Mexico) could lead to erroneous interpretations. This study provides a phylogeographical analysis based on a comprehensive distribution-wide sampling regime and more extensive sampling, and thus should provide new insights into the evolutionary history of the genus.
The genus Astyanax is characterised by high phenotypic plasticity and a capacity to adapt to diverse habitats [36,[38][39][40][41]. There is clear evidence of extremely rapid adaptations of fish to new habitats and environments, with ecological specialization and morphological differentiation, generally in accordance with genetic divergence [15,[42][43][44]. Considerable attention has been given to the evolution of developmental mechanisms and adaptation to cave environments [36,38,39], but less attention has been given to other habitat associated morphological plasticity. In this regard, Bramocharax, which is sympatric with Astyanax, is characterized by conspicuous trophic specializations, including differences in the number of premaxillary teeth, the presence of diastemas on the maxillary teeth, as well as differences in the shape and number of cuspids on the premaxillary, maxillary and dentary teeth, with some species (B. caballeroi and B.baileyi) having intermediate states between the morphotypes of Astyanax and Bramocharax species [45][46][47].
In this study we used mitochondrial and nuclear DNA sequences to develop a robust phylogenetic hypothesis for Astyanax and Bramocharax. This allows us to test biogeographical hypotheses for the Mesoamerican fish fauna, including the relative importance of historical geology and climatic factors.

Results
Three mitochondrial (Cytb, COI and 16 S) genes and one nuclear gene (RAG-1) were sequenced, giving a total of 3862 characters (2350 mitochondrial and 1512 nuclear).
RAG-1 was the most conservative of the genes analyzed (Table 1). Among the mitochondrial DNA genes, Cytb was the most variable, with COI exhibiting similarly high levels of variability and 16 s being the least variable. For the joint mitochondrial and nuclear analysis, 931 sites were variable, with 448 (~11%) being parsimony informative.
The topologies recovered by Maximum Parsimony (MP) and Bayesian Inference (BI) for the Cytb data set ( Figure 1) and combined data set were similar. In addition, Cytb and the combined data set analyses were also concordant, with discrepancies restricted to the tree topologies in Clades III and IV from Group I (the Maquinas population was grouped with Montebello in the Cytb topology but with Polochic-Grijalva-Usumacinta with the combined data set). The combined data matrix was useful to resolve the deeper nodes and recovered mostly higher support values, providing greater phylogenetic resolution. For this reason, description of the higher-level groups identified was based on the topology obtained with the combined data set.
All analyses supported the polyphyly of Bramocharax, with species of Bramocharax being sister groups to different clades of Astyanax (Figure 1), making Astyanax paraphyletic.
We identified six major phylogenetic groupings with high bootstrap support and significant posterior probabilities (Figures 1 and 2). Percentage divergences between groupings are given in Table 2. Groups V and VI correspond to the Chagres region (Panama) and Lagarto-Puntarenas basins of Costa Rica in Lower Central America, respectively. Groups II to IV are from Middle Central America and Group I is from Upper Central America and Mexico. These groups are non-overlapping geographically except for I with II and II with III (see phylogenetic Clades scheme in Additional file 1). Their inter-relationships were not resolved with either the Cytb or the combined data set (3.8 Kbp).

Geographic structuring within the major phylogenetic groupings GROUP I (Mexico and Upper Central America)
Four main clades were recovered from Group I. Clade I comprised most of the Mexican and Upper Central American (Guatemala and Belize) populations. Clades II-IV represented fewer populations with a patchy distribution over the range of Clade I.
Within Clade I we found geographical structure corresponding to the following four lineages: Lineage Id, from the Chiapas region of Mexico, was sister to a clade comprising Lineage Ie from the Candelaria region and a pair of sister lineages (1a and 1b) that occupy a wider region from the Yucatan Peninsula to the Bravo-Conchos basin.  [35,48]. We also found shared haplotypes or low divergence (see Additional file 2) between Bramocharax caballeroi and sympatric A. fasciatus from Lake Catemaco (see region C, Los Tuxtlas, Figure 4). Addi- Summarized Phylogenetic tree estimated with Bayesian Inference and maximum-parsimony methods using the Cytb gene   Lineage Ic, (region D, Figure 3) comprised populations from the Atlantic slope: Astyanax fasciatus (a synonym of A. mexicanus sensu Lima et al. [35]) from the type locality of A. petenensis (Peten -Itza Lake, Yucatan Peninsula) and A. fasciatus from the Candelaria karst region in Guatemala. This lineage included Bramocharax dorioni from its type locality (Semococh river). Similar to other instances of sympatry between Bramocharax and Astyanax, these two morphotypes showed low levels of genetic differentiation (less 0.5%, see Additional file 3).

Sampling sites
Lineage Id (Chiapas) grouped Pacific slope populations of A. fasciatus (Figure 1) Figure 3) with those from the upper Polochic basin (Cahabon river) and the Grijalva-Usumacinta basin.
Clade IV contained populations from the Montebello Lagoons (Montebello) in south-eastern Mexico. This group clustered with the Maquinas and La Palma populations in the analysis of Cytb alone, but this clustering was not supported by the combined nuclear and mitochondrial analysis (Figures 2 and 4).

GROUP II (Middle Central America)
Our analyses did not resolve relationships among the five main clades recovered within Group II. This group occurs widely over Middle Central America, ranging from Belize to Nicaragua and Costa Rica (E, Figure 3).    Subtree of Groups II to VI based on Cytb gene

A. nicaraguensis
Combined data set tree Clade IX included A. fasciatus from a very restricted area comprising the Ciruelas and Tempisque basins on the Pacific slope of Costa Rica ("Ciruelas-Tempisque", Figure  5).

Lower Central America
Four major groups (groups III-VI) in Lower Central America were characterized by more restricted ranges relative to Groups I and II. Group III comprised A. fasciatus populations from the Ciruelas and Chires basins on the Pacific Slope of Costa Rica (figures 1 and 2). Two highly divergent (D K81uf = 2.5%) haplotypes (corresponding to Groups II and III) were found in sympatry in the Ciruelas basin (figures 1 and 2).
Group IV included a single and well-differentiated population of A. orthodus from the Sixaola basin on the Atlantic slope of Costa Rica. Group V contained Pacific slope A. fasciatus populations from the Puntarenas basin (Costa Rica) to the Lagarto basin on the Panama-Costa Rica border. Group VI included A. fasciatus from the Chagres region on the Atlantic slope of Panama.

Systematics of the genera Astyanax and Bramocharax
Our analyses do not support the previously proposed monophyly of Bramocharax based on morphological analyses [45][46][47]. Moreover, Bramocharax specimens were present in two of the seven major Astyanax clades, with low levels of genetic differentiation when both morphotypes were found in sympatry (such as in lake Catemaco where it was possible to find haplotypes shared between individuals from both genera). Differentiation was equally reduced between allopatric populations of Bramocharax and Astyanax.
The genus Astyanax has been considered to be monophyletic in Mesoamerica [37], but polyphyletic in South America [30] on the basis of molecular analyses. Our results support the monophyly of Mesoamerican Astyanax only if we consider Bramocharax species to be morphotypes of Astyanax within the range of its phenotypic plasticity. This hypothesis is supported by the low genetic divergence between specimens of Bramocharax and Astyanax, and the evidence of recurrent evolution of the Bramocharax morphotype within Astyanax (Figure 3). This morphotype is associated with lacustrine habitats, suggesting that its recurrent evolution is a result of morphological convergence to similar ecological factors; similar patterns have been shown in other freshwater fishes [37,38,41,45,49]. If the "recurrent convergence" hypothesis is considered to be correct, then the taxonomy of Bramocharax needs to be revised and the evolutionary mechanisms giving rise to these morphological homoplasies need further investigation. Our analyses question the taxonomic utility of trophic characters (e.g., teeth shape or jaw modification), as previously done by Rosen [47]

on the basis of intermediate morphological states between Astyanax fasciatus and Bramocharax baileyi.
Further incidences of morphological convergence were found in troglobitic morphotypes of Astyanax jordani (this has been noted by previous authors [37,38,50]), providing further evidence of independent (at least two different times, see Figure 6) adaptation to troglobitic habitats. The presence of recurrent morphological convergence in Astyanax [50,51] makes the delimitation of species and genera difficult. Thus the absence of congruence between phylogenetic relationships uncovered in this study and previous taxonomic classifications for Astyanax and Bramocharax from Mesoamerica [35,37,52] is not surprising. In addition, our results are not in agreement with the idea that Astyanax (including samples from Mexico and Upper Central America) is a single species (i.e., A. fasciatus) as has previously proposed [37].
Although not a main goal of this study, we propose a provisional taxonomic nomenclature for Astyanax populations from Mesoamerica (see Additional file 3). The nomenclature proposed is based on well-defined monophyletic groups, high genetic divergences with Cytb (>2% K81uf), and in agreement with geographical distributions. In ascribing species names we gave priority to previous species descriptions and diagnostic morphological traits. Where monophyletic lineages could not be assigned to a valid species name, they were assigned to their own monophyletic group as Astyanax sp.

Genetic and Time Divergences
The penalized likelihood analyses performed in r8s for Cytb sequences was calibrated using the following events: 1) the Merida-Perija uplift about 8-12 Mya, 2) the presence of fossils of Colossoma macropomum in the Magdalena basin (from at least 15 Mya), and 3) the formation of the TMVB about 3-6 Mya ( Figure 6) [53]. The analysis gave an average divergence rate of 0.8% per million years with our in-group and the K81uf model of evolution ( Figure 6). While this is similar to divergence rates reported for other fishes [Cichlidae (0.7%), [22], Cobitidae (0.68%) [54] ( Table 3) and slightly lower than in cyprinid fishes (1.05%) [38], it is much lower than previous molecular clock rates (using fragments from the same gene) proposed for Astyanax (1.5% K2P divergences) [37]. This difference in estimated divergence rate is partially the cause of discrepancies between our study and previous historical biogeographical interpretations for Astyanax [37].

Biogeographical implications
We found a pattern of north-south phylogeographical structuring. The major phylogenetic groups were mostly Ultrametric tree based on Cytb topology using semi-parametric penalized likelihood Figure 6 Ultrametric tree based on Cytb topology using semi-parametric penalized likelihood. Ultrametric tree based on the topology obtained with the mitochondrial Cytb gene using semi-parametric penalized likelihood. The calibration points are indicated by arrows, the first, node 1, corresponds with the rising of the Sierra of Perija and Merida Mountains Díaz de Gamero [79], and the second, node 6, corresponds to the final closure of the Trans-Mexican Volcanic Belt Ferrari et al. [82]  non-overlapping, with the exception of Groups I and II, which overlap in the upper part of the Polochic basin of Guatemala, and Groups II and III, which overlap in the Ciruelas basin of Costa Rica. This north-south pattern is similar to that reported for other freshwater fishes [12,13,22]. We explain the observation of sympatric lineages of Astyanax in terms of niche overlap and lineage turnover, similar to that proposed in biogeographical models for other characins in Mesoamerica [15].
The lack of phylogenetic structuring in Astyanax in Middle Central America (even with the subset of data 3.8 Kbp) can be explained by a more recent colonization and rapid radiation about 6.9 Mya (Table 4); this pattern has been previously observed in other freshwater fishes of the region [12,13,15,22,55].

Dispersal hypothesis on the origin of genus Astyanax in Mesoamerica
We accept the widely held hypothesis of a South American origin for Astyanax and other Central American characids [15,30,33]. This is supported by the observation that Lower Central America lineages were most closely related to South American samples from Brazil and Argentina. Considering the widely used Cytb calibration rate for fish (1.09%/my HKY distances) [22,38,54] and our mean rate of 0.8%/my K81uf distances from R8s ( Figure 6 and Table  4), levels of divergence for populations of South and Central America imply a period of Mesoamerican colonization/expansion of Astyanax from South America about 7.8-8.1 Mya, before the final uplift of the Isthmus of Pan-ama ~3. 3 Mya [7,20,56]. The inclusion of more Astyanax samples from both sides of the Sierra de Perija and Merida Andes in further studies could improve this dating scenario.
The colonization of Central America prior to Late Cenozoic closure of the Panama Strait is incongruent with the geological data, and with other studies of characid genera (Brycon, Bryconamericus, Eretmobrycon and Cyphocharax) [15], including a previous study of Astyanax [37], all of which propose that closure of the strait ~3. 3 Mya provided the first opportunity for colonization of Central America from South America.
Divergence times similar to those found in this study have been reported in molecular studies of primary freshwater fauna [8,12]. For example, the Bermingham and Martin model [8] proposes a colonization of Mesoamerica  [12]. These estimates are in agreement with our study, and coincident with other primary and secondary freshwater fish fauna [12,13,22], as well as divergence times for invertebrates (9 Mya in pseudoscorpions) [26] and benthic foraminifera fossils (8 Mya) [58].
We found evidence of a unique biogeographical pattern involving multiple waves of expansion of Group I (Clades II -V) Astyanax in the upper part of the Polochic-Motagua fault (Mexico and Chiapas region). These clades have a restricted distribution overlapping that of Clade I, and in general occupy relatively stable ecological environments (springs or lakes), which can be less affected by climate change. Niche overlap and lineage turnover could explain this pattern, except in stable habitats where two lineages are found in sympatry (lineages Ia of Clade I and the lineage of Clade II in the Huasteca region, and in the Mezquital and Nazas-Aguanaval basin).

Main vicariant events in Astyanax populations from Mesoamerica
The vicariance events involving Astyanax in Mesoamerica occurred during the Plio-Miocene (4)(5)(6)(7)(8), occurring earlier in Lower Central America (Panama and mainly Costa Rica) than in Central and Upper Mesoamerica.
A pattern of restricted geographic ranges in Lower Central America (Groups IV-VI) supports pronounced geographical fragmentation as a consequence of tectonics movements [7,20], which eventually resulted in closure of the Panama Strait ~3. 3 Mya. Bermingham and Martin [8] have implicated multiple range fragmentation during the Miocene in patterns of diversity in other taxa of the primary freshwater fauna. With our data, five main vicariant events were identified for Lower Central America ( Figure  6). These are related to changes in eustatic sea level (5-8 Mya) [25] and to the formation of inter-oceanic biogeographical barriers [58] during the Middle-Late Miocene (8 Mya).
In Upper Central America, the main volcanic activity was produced by the Trans-Mexican Volcanic Belt (TMVB). This region was affected by periods of intense geological activity between 3 and 12 Mya, with some volcanic activity still occurring today [53,59,60]. The geographic structuring evident in Clade I of Astyanax indicates that the TMVB formed an effective geographic barrier during its development during the late Miocene 4 -6 Mya ( Figure   6). This date is in agreement with the geology of the region and previous studies of several groups of vertebrates [1, 5,10,11,61,62].

Other biogeographical patterns
Other biogeographical patterns were obtained in Lower and Middle Central America. These cannot be explained by geological barriers, but are in accordance with the main biogeographical regions proposed for other freshwater fishes [15,22,63].
While our results contrast somewhat with the Mesoamerican faunal regions recognized by Bussing [14], they do support his Isthmian Region, but with a fragmented pattern similar to that reported for other primary freshwater fishes [8,12,15]. Furthermore, we found that some Belize and Guatemala Astyanax clades (Clades VII and VIII in Fig Bussing's proposal [14], occurring from the Barbilla basin (Costa Rica) to Belize (Atlantic slope) and on the Pacific slope from the Rio Grande Basin to the Ciruelas basin in the Nicoya Gulf (Costa Rica). This pattern has been previously reported by other authors [22,63].
Finally, we dated the separation of Groups I and II to about 6 and 7.8 Mya (Table 4). These coincide geographically with Polochic-Motagua Fault [6], reported as a transition region for other freshwater fish groups [22,13]. In addition, we observed the presence of two well differentiated Lineages of Astyanax (Clade IV in Group I and Clade V in Group II) in sympatry in the Polochic basin, a finding that has not been reported for other characids [15]. We explain this pattern in terms of river capture whereby the Cahabon tributary was diverted to the Polochic river as a consequence of tectonic activity (Sierra de Chiapas [64]), while separating the Cahabon river from the Grijalva-Usumacinta (Clade IV in Group I). This has been proposed for Rhamdia [12].

Tissue collection and DNA extractions
A total of 208 specimens of the Astyanax and Bramocharax from 141 localities from Panama to the Mexican-USA border (Additional file 3; Figure 3) were analyzed, corresponding to 10 Table 1), with the exception of the RAG1 gene for which we followed the PCR conditions described in [67]. PCR products were run on 1.0% agarose gels to confirm amplification and purified with the EXOSAP-IT PCR Product Clean -Up (Usb) kit or by ethanol precipitation. Both strands were sequenced (see Table 5 for primers) and run on an ABI 3700 DNA automated sequencer (SECUGEN sequencing service).

Data analysis
Chromatograms and alignments were visually checked and verified. Saturation for transition and transversion substitutions was checked by plotting the absolute number of changes at each codon position against patristic distances for coding genes only.
Phylogenetic reconstruction was performed for Bayesian Inference (BI) using MrBayes version 3.1.2 [68]. We used Modeltest 3.07 [69] to find the best-fit model of evolution for each gene fragment using the Bayesian Information Criterion (BIC) [69] (see Additional file 4). BI was performed on two data sets as follows: (1) Cytb gene only with separate best-fit models for each codon position and (2) using the three mtDNA and RAG1 (separate best-fit models for each codon position were used) genes with a separate best-fit model of evolution for each gene fragment (partition). Analyses of best-fit models of evolution and BI were performed for a subset of the data.
Bayesian analyses were performed using two independent runs of four Metropolis-coupled chains of 10 million generations each to estimate the posterior probability distribution. The first 10,000 trees were discarded as burn-in. The program Tracer v1.4 [70] was used to assess run convergence and determine burn-in.
Sequence data were also analysed using maximum parsimony (MP) as implemented in PAUP* 4.0 b10 [71], NONA version 2.0 [72] and WINCLADA version 1.00.08 [73]. MP analyses in PAUP* and NONA/WINCLADA The fixed ages are shown in Bold. P = Uncorrected distances ML = K81uf distances. Minimum and maximum values are given in parenthesis.
Asterisks identify values with a normal distribution p  0.05.
were done under the same heuristic search strategy. Statistical support for recovered clades was assessed using bootstrap (1000 pseudo-replications). We applied different weights for transversions and transitions according to the empiric criterion obtained in PAUP* 4.0 b10 [71]. The two datasets run for BI were also run for MP.
Analysis 1 (for which most species and populations of the Astyanax and Bramocharax genera from Mesoamerica were represented) was used to infer phylogenetic relationships among populations. Analysis 2 (for which only a subset of the species/populations were available) was used to infer relationships among the main lineages identified in analysis 1.

Molecular clock and divergence time
Since we have a more complete data matrix, and in order to make our data comparable with previous studies in the region, we calibrated the molecular clock based on Cytb data alone rather than the combined data matrix (which included two more mitochondrial genes and one nuclear DNA gene).
Rate heterogeneity within the dataset was assessed using the likelihood-ratio test (LRT) [74,75]. LRTs were determined by comparing the log likelihood of the optimal topology recovered by maximum likelihood analysis (using appropriate models identified by Modeltest), while enforcing the molecular clock to the log likelihood of the optimal topology recovered by one that did not. The likelihood ratio statistic is twice the difference between the two log likelihoods. This statistic is compared to a  2 distribution with degrees of freedom equal to the number of terminals minus two following [76].
Because the molecular clock hypothesis was rejected, we conducted a non-parametric rate smoothing approach (NPRS) of divergence time estimation with the r8s package [75] to estimate divergence between the taxa. NPRS relaxes the molecular clock assumption by applying a least squares smoothing of estimates of substitutions rates.
Standard errors of divergence dates were estimated using the boot strapping procedure outlined in, and implemented by, Perl scripts in the r8s bootkit provided by Torsten Eriksson [77]. The first 100 bootstrapped datasets were created from the original Cytb dataset with the program Mesquite v. 2.01 [78]. Branch lengths were then reestimated for each bootstrapped dataset in PAUP* using the original ML parameters. The resulting trees, with branch lengths, were then imported into r8s. The TN (Truncated Newton) algorithm was implemented.