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).