V. Divergence time estimation within B. gargarizans
We estimated the divergence times between the major clades of B. gargarizans . Divergence dating analyses were performed using BEAST v.1.7.5 (Drummond et al., 2012) on the cytb dataset, as data from other species of the Bufonidae were publicly available (Recuero et al., 2012). We were unable to include B. gargarizans sampling from Zhan and Fu (2011), as they did not collect cytb data. We focused on the divergence dating results for the major B. gargarizansclades (A, B, G, and H), and within Clade H (Korea, Northeast China, and Russia). Five species were included as outgroups from our own data (B. stejnegeri , B. japonicus ) and GenBank (B. tibetanus , B. melanostictus , Telmatobius bolivianus ) (Table S2 ).
For two calibration points based on fossils, we followed Recuero et al. (2012) by setting a prior distribution for the root of the B. bufo group (B. bufo , B. eichwaldi , B. spinosus , and B. verrucosissimus ) (lognormal distribution with an offset of 9.7 mega-annum (Ma), and 95% of the values between 10.1 and 22.2 Ma) and for B. verrucosissimus (lognormal distribution with an offset of 1.81 Ma and 95% of the values between 2 and 4.5 Ma). The birth–death process was specified for the tree prior, since it was well suited for a multi-species dataset with deep genetic divergence across clades and species. Two independent runs of 100 million generations, sampling every 10 thousand generations, were combined after checking for convergence and adequate effective sample sizes (ESSs) of parameters using the software Tracer v.1.7.1 (Rambaut et al., 2018).