Recent and Rapid Radiation of the Highly Endangered Harlequin Frogs (Atelopus) into Central America Inferred from Mitochondrial DNA Sequences

Populations of amphibians are experiencing severe declines worldwide. One group with the most catastrophic declines is the Neotropical genus Atelopus (Anura: Bufonidae). Many species of Atelopus have not been seen for decades and all eight Central American species are considered “Critically Endangered”, three of them very likely extinct. Nonetheless, the taxonomy, phylogeny, and biogeographic history of Central American Atelopus are still poorly known. In this study, the phylogenetic relationships among seven of the eight described species in Central America were inferred based on mitochondrial DNA sequences from 103 individuals, including decades-old museum samples and two likely extinct species, plus ten South American species. Among Central American samples, we discovered two candidate species that should be incorporated into conservation programs. Phylogenetic inference revealed a ladderized topology, placing species geographically furthest from South America more nested in the tree. Model-based ancestral area estimation supported either one or two colonization events from South America. Relaxed-clock analysis of divergence times indicated that Atelopus colonized Central America prior to 4 million years ago (Ma), supporting a slightly older than traditional date for the closure of the Isthmus. This study highlights the invaluable role of museum collections in documenting past biodiversity, and these results could guide future conservation efforts. An abstract in Spanish (Resumen) is available as supplementary material.


Introduction
Amphibians are experiencing a global conservation crisis, with an estimated 41 to 50% of species suffering population declines [1][2][3]. These declines are likely the result of interactions among several factors of primarily anthropic origin, including overexploitation, habitat destruction, pollution, and the effects of emerging epizootic pathogens [4]. The most alarming cause of global declines and extinctions geological completion of the Isthmus of Panama is a controversial topic, because the traditional date of roughly 3 Ma [32,33] has been challenged recently by Montes et al. [34] on geological evidence, and by Bacon et al. [35] on DNA sequence data, who propose a much older closing of the Isthmus around 15-10 Ma. In turn, O'Dea et al. [36] charged both studies with biased data collection and erroneous interpretation, and the controversy continues. Nonetheless, DNA sequence studies on amphibians and reptiles are often supportive of a closure date older than 3 Ma (e.g., [37][38][39]). In this work, the phylogenetic relationships of the species of Atelopus of Central America are estimated from mitochondrial DNA sequence data in order to evaluate the current taxonomic status of populations and species, as well as to investigate the number and timing of colonization events between South America and Central America.

Sampling
The character data generated for this study consisted of mitochondrial DNA sequences from 103 individuals of Central American Atelopus (Figure 1) plus congeneric samples of individuals from Colombia, representing A. elegans, and A. laetissimus. Samples included toe-clips and tissue samples from vouchered specimens, including many tissue loans from natural history museum collections. DNA sequence data from an additional eight South American species of Atelopus were obtained from GenBank to evaluate the possible monophyly of Central American samples, including data from A.

Molecular Laboratory Protocols
Genomic DNA was extracted using a standard phenol-chloroform procedure.

Phylogenetic Analyses
Phylogenetic inference was performed on the concatenated 2-gene data set totaling 1353 bp using maximum parsimony (MP), maximum likelihood (ML), and Bayesian criteria. The software PAUP * , version 4.0a149 for Windows [44] was used for heuristic tree searching under MP inference, based on 5000 replicate searches (each from random starting trees) with both the MaxTrees command and the rearrangement limit on each tree set to 100,000. These search conditions were repeated five additional times to evaluate the completeness of each search by counting the number of novel most-parsimonious trees obtained per each additional run. In total, 97 equally parsimonious trees were obtained as a result, which differed only slightly in terms of branch lengths and phylogenetic position of some highly similar samples of Atelopus senex, A. varius, A. zeteki, and A. glyphus.
The software PartitionFinder 2 [45] was used to simultaneously select the best-fit partition scheme and corresponding substitution models according to the Bayesian Information Criterion (BIC) and using the greedy algorithm of Lanfear et al. [46]. We assumed a maximum of six possible partitions across the concatenated data matrix (i.e., by gene and by codon). The resulting partition scheme and models are presented in Table S2. Posteriorly, a total of 300 independent tree searches (i.e., three runs with hundred replicates each) were conducted using GARLI version 2.01 [47] on the Cyberinfrastructure for Phylogenetic Research (CIPRES) online platform [48], maintaining the search parameter values as default.
For MP and ML analyses we assessed statistical support of branches via 1000 non-parametric bootstrap replicates [49]. MP bootstrapping was done in PAUP * , with 30 tree searches performed on each character matrix generated by sampling with replacement, and with other search parameters as above. For ML bootstrapping in GARLI, 1000 pseudo-replicates distributed in twenty runs of 50 independent searches each were run on CIPRES, again maintaining the parameter values as default.
Bayesian phylogenetic analysis (BA) using reversible-jump Markov chain Monte Carlo (rjMCMC) was performed with BEAST 2 version 4.8 [50], as implemented in the CIPRES portal, and using the based substitution model implemented in the RBS package of BEAST 2 [50,51]. The posterior distribution of trees including dates and rates was estimated by assuming a relaxed normal clock model of evolution and a birth-death tree prior (a diversification model that considers both speciation and extinction [52]). The prior on substitution rate assumed a normal distribution with a mean of 0.00957 (i.e., 0.957%) and standard deviation of 0.0010 per lineage per million years. The mean was based on the mitochondrial substitution rate of Macey et al. [53], obtained from mitochondrial protein-coding genes in Eurasian toads of the genus Bufo, as corrected by Crawford [54]. The prior distribution of other parameters were not modified from their default values. Bayesian rjMCMC analyses were run for 100 million generations, saving one sampled tree per 1000 generations and with the first 10% of trees discarded as burn-in. Stationarity and mixing of the log-likelihood values and parameter estimates were evaluated with Tracer version 1.7 [55]. We confirmed that all parameters estimated from the post-burn-in set of trees had effective sample sizes > 200. The resulting set of trees was summarized with TreeAnnotator 2.4.8, part of the BEAST 2 package.

Biogeographic Analyses
To determine the number and direction of colonizations between South America and Central America, we inferred ancestral areas using a model-based approach assuming the BA tree trimmed to one sample per species (including unconfirmed candidate species, see below). The R-package 'BioGeoBEARS' [56,57] was used to evaluate the fit of six biogeographical models (DEC, DEC-J, DIVA, DIVA-J, BAYAREA-LIKE, BAYAREALIKE-J) to our data based on the corrected Akaike information criterion (AICc). The model recovered as having the best fit was then used for conducting the ancestral area estimation, again using 'BioGeoBEARS'. Species of Atelopus were coded as being from either Central America or South America, with the dividing line assumed to be the Uramita fault that runs parallel to the Atrato River in Chocó, Colombia [34].

Pairwise Genetic Distances
The genetic distances between sequences of each gene among all samples were estimated with the software MEGA version 7 [58], assuming a Kimura 2-parameter model of sequence evolution (K2P [59]). Default values were retained for all parameters. Sites with ambiguous or missing bases were removed only in pairwise comparisons. An unconfirmed candidate species was identified as a reciprocally monophyletic clade recovered in the BA consensus tree that had a genetic distance from its sister clades larger than the minimum genetic distances between other named sister species (2.08% or 2.44% in COI or cyt b, respectively). Even though the genetic distance cut-offs used herein may seem low in comparison to those used for other frogs and for different loci (e.g., 3% for the case of the 16S ribosomal RNA gene and 10% for the barcode fragment of COI [60,61]), they allow defining clades that mostly correspond to the species recognized by the current taxonomy of the genus in Central America.

Phylogenetic Analyses
The complete, concatenated DNA sequence alignment consisted of a total of 1353 bp from 121 individuals of Atelopus plus samples from six individuals of outgroup bufonids. The majority of ingroup samples had sequences of the COI fragment (unavailable for two samples, or 1.6% of the total), whereas six ingroup terminals (4.7% of the total) lack DNA sequence data from cyt b. Premature stop codons were not observed in any of the sequences. Taking into account the ingroup only, of the 714 sites in cyt b, 142 were parsimony-informative and 35 were singletons, and of the 639 aligned sites of COI, 136 were parsimony-informative and 10 were singletons (Table 1). The trees with the highest likelihood and parsimony scores (one of the 97 more parsimonious and very similar trees was chosen randomly) are shown in Figure 2 and Supplementary Figure S1, respectively, but outgroup non-Atelopus bufonid samples were removed to aid in visualization of the ingroup. The Bayesian consensus phylogeny of ingroup samples is shown in Figure 3, while the complete BA tree with the outgroup is presented in Supplementary Figure S3. The topologies obtained by the three inference methods (MP, ML and BA) were highly congruent with respect to the Central American samples, differing only in the phylogenetic position of A. senex and A. chiriquiensis. The South American Atelopus samples consistently formed three well supported clades and the phylogenetic relationships among the three were distinct among the different inference methods (see below). In general, the majority of nodes in the three phylogenies were highly supported (e.g., bootstrap values > 75% or posterior probabilities > 0.95), with the exception of most internal nodes within A. zeteki and A. varius, among some nodes within Central American Atelopus in the ML topology (see Figure 2), and regarding the phylogenetic relationships between A. senex and A. chiriquiensis.
In all topologies, the genus Atelopus was recovered as monophyletic. The South American congeners, with the exception of A. elegans (see below), comprised three clades, all with significant support: (1) the species of the Amazonian foothills of Peru and Ecuador, plus the Guianan lowlands (Atelopus cf. spumarius, A. pulcher, A. flavescens, and A. cf. hoogmoedi), corresponding to the  (3) A. laetissimus from the Sierra Nevada de Santa Marta of north-western Colombia. In the MP topology, clade 1 was the sister to all other Atelopus, whereas in the ML and BA topologies, A. laetissimus was the sister to all other sampled congeners. Clade 2 was the sister to our focal Central American lineage of Atelopus in both MP and BA topologies, with relatively high support in the BA consensus tree (posterior probability of 0.93; Figure 3 and Supplementary Figure S2). In the ML tree, clades 1 and 2 were sisters and formed the sister lineage to the Central American clade of Atelopus.
The remaining species of Atelopus included in this work comprised a monophyletic group containing eight clades of mostly Central American distribution. All but two of these clades were assignable to a previously described name, because we included samples from near the type localities of all known Central American species of the genus. The two exceptions were that we had no samples of A. chirripoensis, and the group comprised of Caribbean Coast samples from Puerto Obaldía, Panama, and Capurganá, Colombia may represent an unnamed species, referred to as A. aff. limosus in Flechas et al. [62] and Lewis et al. [23]. This latter clade plus A. elegans from the Pacific coast of Colombia ( Figure 4D5) together formed the sister lineage of all other Central American species in all analyses ( Figures 2 and 3). The first split within this clade of Central American endemics separated (from all remaining samples) two species from eastern Panama, A. certus and A. glyphus, which were recovered as sisters in all analyses. The next split separated from all remaining samples either A. senex plus A. chiriquiensis (BA) or a clade containing A. senex, A. chiriquiensis, plus the central Panama species, A. limosus, A. varius, and A. zeteki (ML and MP), although neither of these hypotheses received high statistical support. In trees inferred by all three methods, A. limosus of central Panama was the sister group to a clade containing samples identified as A. varius and A. zeteki. Atelopus varius (e.g., Figure 4B3,B4) ranges from central Costa Rica to central Panama, and A. zeteki (e.g., Figure 4B5,C1,C2) is found near and east of Penonomé in central Panama. Specimens of these two species were found in sympatry in Juan Lana, near San Miguel Abajo, Panama. The identification of the Atelopus from the Monteverde Cloud Forest Reserve in northern Costa Rica ( Figure 4A1,A2), as A. varius, was not supported by the phylogenetic analyses. Therefore, A. varius seems to have a smaller distribution range in Costa Rica, whose type locality is "Veragoa", western Panama [63,64].
For both mitochondrial genes, the genetic distances between the eight clades of Central American Atelopus were generally much larger than within any clade (Table 2), with the exception of the mtDNA clades designated as A. varius and A. senex. However, this pattern is driven by the individual with institutional code MVZ 164818 (from Monteverde, Costa Rica) which has genetic distances comparable to, or even larger than, those found among the eight clades previously mentioned ( Table 2). Each of the eight clades had genetic distances greater than 2.08 at COI (and greater than 2.44 at cyt b) from their sister clade (Table 2), except for the case of A. glyphus, for which its estimated COI genetic distance (1.43-2.07) from its sister species is lower than the cut-off value. below). In general, the majority of nodes in the three phylogenies were highly supported (e.g., bootstrap values >75% or posterior probabilities >0.95), with the exception of most internal nodes within A. zeteki and A. varius, among some nodes within Central American Atelopus in the ML topology (see Figure 2), and regarding the phylogenetic relationships between A. senex and A. chiriquiensis.    Figure 3. A timetree of Atelopus diversification from a relaxed-clock MCMC Bayesian analysis using the software BEAST 2, calibrated using a mitochondrial DNA substitution rates of mitochondrial genes of bufonid frogs of Macey et al. [53], as corrected by Crawford [54]. Outgroup samples of non-Atelopus bufonid genera were trimmed for improved visualization of the ingroup and are shown in Supplementary Figure S3. Numbers above branches are the posterior probabilities of adjacent nodes. 95% credibility interval of the age of each node is indicated by horizontal bars. The grey vertical shading shows an interval between 3.2 to 2.76 million years ago (Ma) indicating the estimated time when gene flow among shallow marine animals and movement of surface water across the Panamanian Isthmus ended, according to O'Dea et al. [36]. Note that the minimum age of node H, the oldest node reconstructed as being entirely Central American, is still older than 3.2 Ma. Shading in the circles on each node correspond to the ancestral area reconstructions under the DEC model (white refers to South America, pale grey to Central America, and black to ancestors recovered as inhabiting both continents), and the letters in each circle correspond to rows in Table 3. Specimens are indicated by their field or museum voucher number, or if not available, with an arbitrary number and the name of the locality where they were obtained (sample details in Supplementary Table S1). Specimens from the locality in which A. varius and A. zeteki are found in sympatry (Juan Lana, Panamá) are marked with asterisks. Figure 3. A timetree of Atelopus diversification from a relaxed-clock MCMC Bayesian analysis using the software BEAST 2, calibrated using a mitochondrial DNA substitution rates of mitochondrial genes of bufonid frogs of Macey et al. [53], as corrected by Crawford [54]. Outgroup samples of non-Atelopus bufonid genera were trimmed for improved visualization of the ingroup and are shown in Supplementary Figure S3. Numbers above branches are the posterior probabilities of adjacent nodes. 95% credibility interval of the age of each node is indicated by horizontal bars. The grey vertical shading shows an interval between 3.2 to 2.76 million years ago (Ma) indicating the estimated time when gene flow among shallow marine animals and movement of surface water across the Panamanian Isthmus ended, according to O'Dea et al. [36]. Note that the minimum age of node H, the oldest node reconstructed as being entirely Central American, is still older than 3.2 Ma. Shading in the circles on each node correspond to the ancestral area reconstructions under the DEC model (white refers to South America, pale grey to Central America, and black to ancestors recovered as inhabiting both continents), and the letters in each circle correspond to rows in Table 3. Specimens are indicated by their field or museum voucher number, or if not available, with an arbitrary number and the name of the locality where they were obtained (sample details in Supplementary Table S1). Specimens from the locality in which A. varius and A. zeteki are found in sympatry (Juan Lana, Panamá) are marked with asterisks.

Divergence Times and Ancestral Area Estimation
The Dispersal and Local Extinction and Cladogenesis biogeographical model (DEC [67]) was found to have the best fit to the geographic distribution of species and the BA phylogeny (Figure 3, Supplementary Table S3). The ancestral area estimation was mostly well-supported with probabilities > 0.95 across nodes (Table 3 and Supplementary Table S4). Combining these results with the BA timetree suggested that the most recent common ancestor (MRCA) of the genus Atelopus existed in the early Miocene of South America, 21.1 Ma (95% posterior credibility interval, CI: 14.21 to 28.95 Ma). In addition, the most recent common ancestor (MRCA) of the Central American species of the genus (plus the Colombian A. elegans) was recovered as diverging from the other species of Atelopus during the late Miocene, (i.e., stem age at 10.6 Ma, CI: 7.02 to 14.64 Ma; Figure 3, node J), with a high probability (0.958) that it was distributed in both Central and South America at the time. Atelopus sp. "Puerto Obaldía-Capurganá" (aka, A. aff. limosus) plus A. elegans shared a MRCA 2.72 Ma (late Pliocene; CI: 1.8 to 3.69 Ma) which was inferred to inhabit both Central and South America, while all members of the sister lineage (Figure 3, node F) to this clade were recovered as being firmly of Central American origin. The timetree allowed us to bound the timing of dispersal into Central America as follows. The hypothesis of just one colonization of Central America (followed by one back dispersal event by Atelopus elegans) implies that the colonization took place between the stem age ( Figure 3, node J; 10.6 Ma) and crown age (node H at 6.01 Ma, CI: 4.11 to 8.1 Ma) of the MRCA of all Central America samples, i.e., between roughly 10.6 Ma and 6.01 Ma. The alternative scenario of two independent invasions implies the age of first and more successful (in terms of descendants) took place between the stem age (node H at 6.01 Ma) and crown age (node F at 4.79 Ma, CI: 3.29 to 6.47 Ma) of the oldest clade of exclusively (i.e., monophyletic) Central American samples (Figure 3). Table 3. Estimated crown ages (in millions of years ago, Ma) and results of the ancestral area estimation using the Dispersal-Extinction-Cladogenesis (DEC) model for the species of Atelopus of Central America and other selected nodes within the genus. Ages were estimated by concatenated MCMC Bayesian analyses using the software BEAST 2. See text for details. Node labels are indicated on the tree in Figure 3. CA = Central America, SA = South America.

Phylogenetic Systematics and Biogeography
The phylogenetic relationships inferred here among Atelopus are mostly in agreement with those found by Lötters et al. [10], a study also based solely on mitochondrial DNA sequences (12S and 16S genes), but which included only three of the Central American species, A. varius, A. zeteki, and A. chiriquiensis. The current study and Lötters et al. [10] found that Atelopus is monophyletic and of South American origin. Additionally, both studies recovered a clade of Amazonian-Guianan species which is sister to a clade containing the Andean, Chocoan, and Central American species. In contrast with Lötters et al. [10], however, in this work we included sequences of one species (A. laetissimus) from the Sierra Nevada de Santa Marta (SNSM), an isolated massif in north-western Colombia well known for its endemicity. The phylogenetic position of A. laetissimus is noteworthy, as it diverged very early in the history of the genus, being either sister to the rest of Atelopus (Figures 2 and 3), or the sister to a clade formed by Andean, Chocoan, and Central American species (Supplementary Figure S1). A similar pattern was found in the frogs of the clade Allocentroleniae, in which the basal split separates the Guianas and Amazonian regions from a clade containing a species from the SNSM that is sister to all other known species [68].
While a single colonization event can explain the presence of Atelopus in Central America, our data are also compatible with two dispersal events (Figure 3). Further sampling of other Chocoan Atelopus, such as A. spurrelli, A. longibrachius, A. balios, and the mainland populations of A. elegans, may be needed to resolve this ambiguity. Savage and Bolaños [16] proposed that the enigmatic A. chirripoensis of Costa Rica is more closely related to the A. ignescens species complex of Ecuador and southern Colombia than to other Central American species of Atelopus, which would imply yet another potential colonization event. However, testing this hypothesis is difficult given that A. chirripoensis has not been seen since 1980, and it is now probably extinct [16,18]. The remarkable external similarity between A. chirripoensis and the A. ignescens complex might also be explained by morphological convergence, given that both are found in similar páramo environments [16].
Our analyses suggest that the ancestral lineage of Atelopus reached Central America between 14.64 and 3.29 Ma with point estimates ranging from 10.6 to 4.79 Ma, depending on whether one postulates one or two invasions. In all cases, the estimated time frame is consistent with recent molecular biogeographic studies (e.g., [38,69]) supporting a non-traditional, older closure date for the Panamanian Isthmus by 10 Ma [34]. One could also argue that Atelopus dispersed over water prior to the completion of the Isthmus, although amphibians are sensitive to dehydration and to salt water [66,70,71]. Assuming an old date for the formation of the Isthmus, one recent hypothesis predicts that the Isthmus was forested during the Miocene and Pliocene epochs until 2.5 Ma when savanna-like environments appeared [72]. Our proposed time-frame for colonization by Atelopus of the Isthmus is temporally consistent with this prediction, as Atelopus species are affiliated with moist forests. Studies of additional groups of terrestrial organisms with low dispersal capabilities (such as insects, fossorial reptiles, or understory birds) or that are sensitive to salt water (such as amphibians) are needed to confirm whether there may be a tendency of inferring early crossings between continents by wet-forest lineages [31].

Geographic Color Pattern Variation
Most Central American species of Atelopus are strikingly variable in their color patterns, not only between but also within populations (e.g., [63]). In some species, such as A. chiriquiensis ( Figure 4B1,B2), A. senex ( Figure 4A3,A4) and populations of A. varius, a noticeable sexual dimorphism in coloration is sometimes present. On the other hand, color patterns can be very similar between species, resulting in species misidentification [28,73], especially in species that lack other clear diagnostic morphological characters. Our molecular data revealed the presence of a distinct species in the Monteverde Cloud Forest Reserve, northern Costa Rica, that has until now been confused with A. varius, very likely due to their similarity in color pattern.
In general, recently metamorphosed and juvenile individuals of A. certus, A. chiriquiensis, A. glyphus, A. limosus, A. varius, and A. zeteki have a barred dorsal color pattern, consisting of a series of dark and light colored chevrons, which usually changes or fades ontogenetically as they grow older ( [74], E.D.L. and R.I. pers. obs.). This barred pattern sometimes remains in adult individuals within certain populations, and may even be the predominant color pattern, as observed in some populations of A. glyphus ( Figure 4D1), A. limosus ( Figure 4C4), A. varius ( Figure 4B3,B4), and A. zeteki. Therefore, natural and/or sexual selection pressures might be behind this ontogenetic change in coloration, resulting in color pattern differences among geographic areas.
Geographic variation in color pattern has been documented for A. varius and A. zeteki [63,73], and we found distinct color pattern variants among populations of A. glyphus and A. limosus along their distribution ranges. In the northern region of the Serranía de Pirre (e.g., Figure 1, localities 5 and 6), adults of A. glyphus have a barred dorsal color pattern ( Figure 4D1), while in the highlands of the southern portion of Serranía de Pirre above Cana (e.g., Figure 1, locality 7), adults have a predominantly uniform brown dorsum with small, light markings ( Figure 4D2). In the western area of the distribution range of A. limosus (e.g., Figure 1, localities 8 and 10), adults have a uniform olive green dorsum ( [63,75]; Figure 4C3), whereas populations to the East (e.g., Figure 1, locality 9) have a barred color pattern on the dorsum ( Figure 4C4). The mtDNA data presented here supports our contention that these color pattern differences indeed represent geographic variation between conspecific populations. In contrast, the undescribed Atelopus sp. "Puerto Obaldía-Capurganá" from eastern Panama seems to exclusively have a barred dorsal color pattern ( Figure 4C5).

Taxonomic and Conservation Implications
The phylogenetic results obtained here are broadly concordant with the current taxonomic arrangement of the species of Atelopus. Our data reveal two Central American lineages whose taxonomic status needs further attention, however. First, the specimens from Puerto Obaldía, Panama, and Capurganá, Colombia, were already of uncertain taxonomy. The names A. varius glyphus [76,77], A. varius [78] and A. aff. limosus [23,62] have been applied to this binational population, but our molecular phylogenetic results are not compatible with these proposed names. We cannot, however, reject the hypothesis that these Caribbean samples might correspond to a previously unreported population of A. spurrelli from the Pacific coast of Colombia, since we were unable to sample this latter taxon. More likely, however, these specimens represent an undescribed species.
Second, a specimen from the Monteverde Cloud Forest in Costa Rica (MVZ 164818) is most likely a member of an undescribed (and probably extinct) species, given its large genetic distance from all other specimens, including several sympatric individuals of its sister lineage, A. senex. The specimens from the Monteverde, Chompipe, and Tapantí areas in Costa Rica highlight a taxonomic problem that requires further study. Specimens from populations at Chompipe (the slopes of Volcán Barba) and Tapantí are traditionally considered to be A. senex ( Figure 4A3-A5), while the ones from around Monteverde are assigned to A. varius ( [15,63,79]; Figure 4A1,A2). The specimen (MVZ 164818) from the Monteverde Cloud Forest Reserve stood as a member of an undescribed (and probably extinct) species, although its phylogenetic position had no significant statistical support. This specimen was inferred to be the sister clade to either A. senex, A. chiriquiensis, or to a larger clade containing both. Richards and Knowles [28] also reported on the taxonomic confusion around Monteverde specimens and noted the low genetic divergence (i.e., 0.5-1.2%) between samples from Monteverde and A. senex specimens. However, we confirm that the specimen MVZ 164818 from Monteverde is certainly not A. varius, and it should likely be assigned to an undescribed species. Furthermore, the specimens from Monteverde MVZ 164816 and MVZ 164818 share the same morphology and color pattern based on our examination, but here we show they are genetically distinct at the level of species, as A. senex and A. sp. "Monteverde", respectively. Assuming no errors on the original source of the samples and DNA sequencing, possible explanations for the difference between morphological features and our mtDNA analysis are that genetic introgression between species has occurred or both species and/or their hybrids were present at Monteverde. To resolve this taxonomic problem, we suggest including more DNA data in future analyses, since further sampling of individuals would be difficult given the catastrophic population declines. Finally, the poorly known A. spurrelli, A. certus ( Figure 4D3,D4), and A. glyphus ( Figure 4D1,D2) are in need of further taxonomic study.
We confirm the results of Richards and Knowles [28] regarding the distinctiveness of A. zeteki ( Figure 4B5) from A. varius, and we uncovered the locality (Juan Lana, Panama) in which these species are found in sympatry, as was hinted by Zippel et al. [73]. As previously suggested [28,73], these two species could have been hybridizing in sympatry at this locality, a hypothesis that deserves further exploration given their morphological similarity and recent divergence (2.12 Ma, CI: 1.41 to 2.86 Ma). Additionally, we established the identity of individuals from a population in Cerro Azul (Cerro La Victoria), Panama ( Figure 4C2). Savage [63] suggested this to be an introduced population of A. varius that resembled those from El Valle de Antón and Cerro Campana. According to our phylogenetic analyses, the individuals from this population are A. zeteki.
Finding support for the traditional taxonomy of Central American Atelopus brings hope to the success of the ex situ conservation programs being undertaken with the described species, as it indicates that resources are being used efficiently in frog conservation, given that all surviving species are being included by these programs. In addition, since each named species appears to be a distinct genetic entity, species ranges are not being underestimated, and no species is being protected unnecessarily. However, this is not the case for the two candidate species uncovered in this study, as they are not currently part of any ex situ or in situ conservation program. Furthermore, finding that the original identification of most specimens (except some specimens of A. varius and A. zeteki, Supplementary Table S1) is concordant with the results of the phylogenetic analyses, suggests that cryptic diversity is low among Central American Atelopus, and that hybridization resulting from misidentifications of specimens in captive programs should be rare. At least in the case of A. limosus this has been shown to be the case [80]. Further study is needed using multilocus nuclear markers or genomic data, especially on the likely case of hybridization between A. varius and A. zeteki. Additional genotyping should be conducted on captive specimens currently housed in ex situ programs. Finally, recognizing the taxonomic distinctiveness of the candidate species found herein is urgently needed in order to initiate the corresponding conservation measures to guarantee their long-term survival.

Conclusions
In this work, the phylogenetic relationships among seven of the eight described species of Central American Atelopus, plus eight South American congeneric species, were inferred based on mitochondrial DNA sequence data. The phylogenetic analyses revealed a ladderized topology showing the Central American species as a monophyletic group, and placing the species geographically furthest from South America more nested in the tree. We detected two previously unrecognized candidate species, including an undescribed species from Costa Rica involved in a taxonomic confusion that requires further study. We showed that species in eastern Panama, also show geographic variation in their color patterns, and clarify the species identity of individuals from some of these populations. Biogeographic models supported either one or two colonization events from South America, indicating that Atelopus reached Central America prior to 4 million years ago (Ma), a timing slightly older than the traditional date estimated for the closure of the Isthmus. Furthermore, this study underscores the invaluable role of museum collections in documenting biodiversity, and the relevance of genetic analyses for guiding conservation efforts.
Supplementary Materials: The following are available online at http://www.mdpi.com/1424-2818/12/9/360/s1; Resumen (Abstract in Spanish); Figure S1: One of the 97 most-parsimonious trees (selected arbitrarily) inferred from DNA sequences of two mitochondrial genes (COI and cyt b) for samples of Atelopus of Central and South America; Figure S2: Maximum likelihood phylogenetic tree inferred from DNA sequences of two mitochondrial genes (COI and cyt b) from species of Atelopus of Central and South America, including the extended non-Atelopus bufonid outgroup; Figure S3: Timetree of Atelopus based on a relaxed-clock MCMC Bayesian analysis of DNA sequences of two mitochondrial genes, COI and cyt b, calibrated using the estimated substitution rate for bufonid frogs of Macey et al. [53], as corrected by Crawford [54], and including the extended non-Atelopus bufonid outgroup; Table S1: Specimens of Atelopus and other bufonid frogs used for the present study, including museum voucher, locality, geographic coordinates, elevation and GenBank accession numbers; Table S2: Best-fit partitioning scheme for the two protein-coding mitochondrial genes, cytochrome b (cyt b) and cytochrome oxidase subunit I (COI) and models of nucleotide substitution for each partition, based on the Bayesian information criterion (BIC) as implemented in PartitionFinder 2 [45]; Table S3: Comparison among biogeographic models evaluated with the software BioGeoBEARS; Table S4: Estimated crown ages (in millions of years ago, Ma) and results of the ancestral area estimates under the Dispersal-Extinction-Cladogenesis (DEC) model, assuming the MCMC Bayesian consensus tree obtained using the software BEAST 2 and concatenating the two mitochondrial genes.