*2.3. Phylogenetic Analyses*

For both large and small datasets, we performed both maximum likelihood (ML) and multispecies coalescent-consistent phylogenetic analyses. We used IQ-TREE v1.5.5 [67] to perform our ML analyses, using a general time-reversible (GTR) model and assessing support with 10,000 ultrafast bootstrap replicates [68]. ML analyses were performed on an unpartitioned concatenated matrix (large matrix: 787,199 characters; small matrix: 786,510 characters) to increase computational e fficiency. As UCEs are usually not protein-coding, it is unclear which partition schemes should be used for them, or whether they should be used at all [55].

We also inferred the dendrobatid species tree using ASTRAL-III v5.6.1 [69], a summary method consistent with the multispecies coalescent. ASTRAL-III accounts for incomplete lineage sorting by summarizing gene trees constructed separately for each locus, rather than e ffectively assuming the whole set of loci acts as a single gene, as in a concatenated ML analysis. We made individual gene trees for each UCE locus in IQ-TREE with a GTR model and 1000 ultrafast bootstrap replicates. We contracted near-zero branch lengths to polytomies in the gene trees with IQ-TREE's *-czb* option, an approach recommended by Persons et al. (2016) [70] to reduce downstream bias. We used the gene trees as input for ASTRAL-III. For the large dataset, we assigned each sample to one of 38 putative species in a mapping file used as input with ASTRAL's *-a* option. For the small dataset, we omitted the mapping file so that ASTRAL-III correctly assumed that each sample corresponded to its own species.

## *2.4. Divergence Time Estimation*

We performed divergence time estimation of our dendrobatid phylogeny using BEAST 2 v2.5.0 [71]. Because Bayesian methods are computationally intensive, we reduced our small dataset, consisting of 1706 loci for 37 samples, to four subsets of 50 random loci each for analysis in BEAST, amounting to a total matrix size of 92,742 characters. We analyzed each subset twice in order to ensure that convergence occurred. We concatenated the loci in each subset and did not partition the alignment in order to avoid the intractably long periods of time a partitioned Bayesian analysis can take to converge [72].

We used the utility BEAUti to specify our BEAST settings. We used an HKY model with 4 gamma rate categories, with base frequencies set to "empirical," avoiding the GTR model because it can lead to overparameterization and subsequently low ESS values [73]. We used a relaxed log-normal clock model with a clock rate prior of 1e-10, the same order of magnitude as an estimate of avian UCE substitution rates from Winker et al. [74]. To further reduce computational demands, we fixed the analysis to our small ASTRAL topology by setting the *subtreeSlide*, *narrowExchange*, *wideExchange*, and *wilsonBalding* operators to zero, an approach used by Hsiang et al. [75]. We used a Yule tree prior with other priors set to their default values.

For our divergence time calibration, we used an indirect calibration derived from the timetree provided by Santos et al. (2009) [44]. Divergence time estimation in dendrobatids is a difficult problem because they lack a fossil record. Santos et al. calibrated their dendrobatid tree by nesting it within a tree for the whole of Amphibia, which has a voluminous fossil record, and using a combination of paleogeographic, fossil, and molecular clock evidence to date that amphibian tree. For our study, we use Santos et al.'s estimation of the divergence between Dendrobatidae and Aromobatidae (i.e., the node separating *Allobates femoralis* from the rest of our samples). The average of Santos et al.'s three estimations of this node's age was 38.534 Ma with σ = 5.151 Ma. Our calibration assigned a normal distribution with these values to this node.

We ran each analysis with an MCMC chain length of 100,000,000 generations, with a log sampling frequency of 100,000 generations and a tree sampling frequency of 10,000 generations. We assessed convergence between runs of each subset and ESS values for each run with Tracer v1.7.1 [76]. We found that all parameters had ESS values over the popular threshold of 200, and that convergence was reached for each subset. We used LogCombiner v2.5.0 [71] to combine the posterior distributions of trees from each of the eight runs, accounting for a burn-in of 10%. We then used TreeAnnotator v2.5.0 [71] to summarize the combined tree files, targeting the maximum clade credibility tree with mean node heights.

## *2.5. Visualizing Evolutionary History, Conservation Risk, and Spatial Distributions*

To better understand how the evolutionary history and the geographic distribution of poison frogs impacts their conservation status and extinction risk, we downloaded IUCN Red List classifications for all surveyed species of Dendrobatidae and its sister family Aromobatidae (*n* = 285) [77]. For each species, we recorded its population status, Red List status, and its countries of occurrence. The relationships between these factors were visualized in two circular plots created in the R package *circlize* [78]. The phylogenetic results of this study and those of Grant et al. (2017) [41] characterized the genus-level relationships among all input species.
