*2.1. 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.* cf. *spumarius*, *A. flavescens*, *A.* cf. *hoogmoedi*, *A. nanay*, *A. peruensis*, *A. pulcher*, *A. epikheistos*, and *A. bomolochos*. As outgroup to the genus *Atelopus*, we included additional GenBank sequence data representing the bufonids *Bufo gargarizans minshanicus*, *B. japonicus*, *B. stejnegeri*, *Duttaphrynus melanostictus*, *Anaxyrus boreas*, and *Bufotes viridis*. See Supplementary Table S1.

**Figure 1.** Map of Isthmian Central America showing the collecting localities of the samples used for this study. Each species is indicated with a di fferent color, as follows: *Atelopus chiriquiensis*, grey; *A. senex*, white; *A. varius*, yellow; *A. zeteki*, orange; *A. limosus*, green; *A. certus*, red; *A. glyphus*, brown; *A.* sp. "Puerto Obaldía-Capurganá", light blue; *A.* sp. "Monteverde", violet. The yellow|orange and violet|white circles indicate the localities where *A. varius* and *A. zeteki* or *A.* sp. "Monteverde" and *A. senex*, respectively, were found in sympatry.

## *2.2. Molecular Laboratory Protocols*

Genomic DNA was extracted using a standard phenol-chloroform procedure. A 714-base pair (bp) fragment of cytochrome *b* (cyt *b*) was amplified with the primers CB1 (5--CCATCCAACATCTCAGCATGATGAAA-3-) and CB3 (5--GGCGAATAGGAAGTATCATTC-3- [40]). A 639-bp fragment of the 3- end of cytochrome oxidase I (COI) was amplified using the primers

COIf (5--CCTGCAGGAGGAGGAGAYCC-3-) and COIa (5--AGTATAAGCGTCTGGGTAGTC-3- [41]). PCR, DNA sequencing, and contig assembly protocols followed those of Crawford et al. [42]. As no internal gaps were inferred, sequences were aligned by eye and no premature codons were found when translating inferred codons with the software Mesquite version 3.10 [43].

## *2.3. 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<sup>∗</sup>, 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 rearrangemen<sup>t</sup> 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 di ffered 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<sup>∗</sup>, 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 e ffective sample sizes > 200. The resulting set of trees was summarized with TreeAnnotator 2.4.8, part of the BEAST 2 package.

## *2.4. 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].

## *2.5. 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 cy<sup>t</sup> *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.
