*3.1. 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 cy<sup>t</sup> *b*. Premature stop codons were not observed in any of the sequences. Taking into account the ingroup only, of the 714 sites in cy<sup>t</sup> *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).

**Table 1.** Number of parsimony-informative sites, singletons, and invariant sites (in base pairs) for each gene region among the ingroup samples.


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 congruen<sup>t</sup> 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 'Amazonian-Guianan clade' of Lötters et al. (2011); (2) species from the high Andes of Ecuador and Peru (*A. epikheistos*, *A. nanay*, *A. bomolochos*, and *A. peruensis*); and (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.* <sup>a</sup>ff. *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 cy<sup>t</sup> *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-o ff value.

**Figure 2.** Maximum likelihood (ML) tree inferred from DNA sequences of two mitochondrial genes (COI and cy<sup>t</sup> *b*) from the species of *Atelopus* of Central and South America. Outgroup samples consisting of non-*Atelopus* bufonids were removed for better visualization of the ingroup and are shown in Supplementary Figure S2. Branch support for all major clades is indicated as non-parametric bootstrap values under ML and maximum parsimony (MP) above and below each branch, respectively. However, to improve visualization of the tree, support values for only the earliest-diverging lineages of within-population samples are shown. A dash indicates MP bootstrap support below 50%. Specimens are indicated by their field or museum voucher number, or if not available, by the name of the locality where they were obtained. Specimens from the locality in which *A. varius* and *A. zeteki* are found in sympatry (Juan Lana, Panama) 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.

**Figure 4.** Geographic color pattern variation of some harlequin frogs from central to eastern Panama and adjacent Colombia. Rows are arranged from west (**1**) to east (**5**) for columns (**A**–**C**), and north (**1**) to south (**5**) for column (**D**). Columns represent regions from Central America to Colombia: (**A**) Costa Rica; (**B**) western to central Panama; (**C**) central to eastern Panama and adjacent Colombia; (**D**) eastern Panama and Darién Province, Panama (**1**–**4**) to adjacent Colombia (**5**). Unscaled photographic images of: (**A1**,**A2**) *Atelopus* sp. "Monteverde" from the Monteverde Cloud Forest Reserve, Punta Arenas-Alajuela Provinces, Costa Rica (Figure 1, locality 12); (**A3**) male and (**A4**) female *Atelopus senex* from Volcán Barba, Parque Nacional Braulio Carrillo, Heredia Province, Costa Rica; (**A5**) male *Atelopus senex,* MVZ 203769, from the Refugio Nacional Tapantí, Costa Rica (Figure 1, locality 13); (**B1**) male and (**B2**) female *Atelopus chiriquiensis* from Río Candelaria, Parque Internacional La Amistad, Panama; (**B3**) *Atelopus varius* from Cerro Tute, Parque Nacional Santa Fe, Panama (Figure 1, locality 27); (**B4**) *Atelopus varius* from Río Tigrero, Coclé Province, Panama (Figure 1, locality 20); (**B5**) *Atelopus zeteki*, sample K37, Quebrada Juan Lana, San Miguel Abajo, Coclé Province, Panama (Figure 1, locality 24); (**C1**) *Atelopus zeteki* from the headwaters of Río María, Sorá, Panamá Oeste Province, Panama (Figure 1, locality 37); (**C2**) *Atelopus zeteki* from Cerro Azul (Cerro La Victoria), Panamá Province, Panama (Figure 1, locality 31); (**C3**) *Atelopus limosus* from type locality, San Juan de Pequení, Panamá Province, Panama (Figure 1, locality 8); (**C4**) *Atelopus limosus* from Nusagandi, Comarca Guna Yala, Panama (Figure 1, locality 9); (**C5**) *Atelopus* sp. "Puerto Obaldía-Capurganá" from Capurganá, Chocó Department, Colombia (Figure 1, locality 15); (**D1**) *Atelopus glyphus*, sample no. 10, from north-western Serranía de Pirre, Panama (Figure 1, locality 5); (**D2**) *Atelopus glyphus*, CH 5613, from southern Serranía de Pirre, Panama (Figure 1, locality 7); (**D3**) *Atelopus certus* from Cerro Sapo, Darién Province, Panama; (**D4**) *Atelopus certus*, CH 4665, from Serranía de Jingurudó, Darién Province, Panama (Figure 1, locality 1); (**D5**) *Atelopus elegans* from Isla Gorgona, Cauca Department, Colombia. Some of Central American *Atelopus* have bright colors with contrasting patterns and are known to have potent skin toxins [65] and are, therefore, considered to be aposematic [66]. Photo credits: David Cannatella (**A5**), David M. Dennis (**A1**,**B1**,**B2**,**C2**,**C4**), Sandra V. Flechas (**C5**,**D5**), Michael and Patricia Fogden (**A2**,**A3**,**A4**), Marcos A. Guerra (**C3**,**D1**,**D2**), Roberto Ibáñez (**D4**), César A. Jaramillo (**D3**), Erik D. Lindquist (**B3**,**B4**), and Guido Sterkendries (**B5**,**C1**).


**Table 2.** Pairwise genetic distances obtained using a Kimura 2-parameter model of sequence evolution (i.e., K2P distance × 100) for the mitochondrial fragments of the COI and cyt *b* genes of the species of *Atelopus* from Central America. Ranges represent minimum and maximum observed distances. Intraspecific genetic

## *3.2. 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.* <sup>a</sup>ff. *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.

