Sample collection

57 individuals of B. g. miyakonis were collected from 53 locations throughout the Miyako Islands (Fig. 1B, Appendix I). All animals were euthanized using CO₂ (Suzuki and Kurosawa, 2005) and were then stored at –20°C.

Laboratory procedures

Total DNA was extracted from skeletal muscles using DNeasy Blood and Tissue Kit (QIAGEN, Hilden, Germany). Sequences from the mitochondrial control region of B. g. miyakonis were amplified by polymerase chain reaction (PCR). PCR was performed using ExTaq Kit (TaKaRa Bio., Kusatsu, Japan) using the following primers: CytbA-L (5′-GAA TYG GRG GWC AAC CAG TAG AAG ACC C-3′) and ControlB-H (5′-GTC CAT TGG AGG TTA AGA TCT ACC A-3′) (Goebel et al., 1999). PCR amplification was performed on a MiniAmp Thermal Cycler (Applied Biosystems, Massachusetts, USA) using the following protocol: predenaturation at 95°C for 5 min.; 29–35 cycles of denaturation at 95°C for 30 s, annealing at 65.5–67.5°C for 30 s, extension at 72°C for 90 s; and final extension at 72°C for 5 min. The PCR products were cleaned up by polyethylene glycol precipitation before determination of sequence identity by cycle sequencing. Cycle sequence reactions used a BigDye Terminator Cycle Sequencing Kit (Applied Biosystems) and the PCR primers listed above along with the following protocol: 95°C predenaturation for 1 min; 30 cycles of denaturation at 95°C for 10 s, annealing at 65.5–67.5°C for 5 s, and extension at 60°C for 4 min. Reaction products were then cleaned up by ethanol precipitation. The Sequences were read on an ABI PRISM 3130xl Genetic Analyzer (Applied Biosystems). All fragments were sequenced both in the forward and reverse senses and were assembled using Gene-Studio v.2.2.0.0 (Genestudio Inc., Suwanee, USA).

Inference of the genealogy of the Bufo gargarizans species group

To infer the phylogenetic position of B. g. miyakonis within the B. gargarizans species group, we used sequences from the mitochondrial control region. Although the obtained sequences reached a maximum of 348 bp, we used a 341 bp sequence that was shared by all sequences of the B. gargarizans species group. The sequences of B. g. miyakonis were aligned by ClustalW, as implemented in MEGA X (Kumar et al., 2018), with sequences from continental populations of B. gargarizans (Fu et al., 2005; Cao et al., 2006; Hu et al., 2007; Borzée et al., 2017; Tong et al., 2017) and Taiwanese populations of B. bankorensis (Liu et al., 2000; Yu et al., 2014) obtained from GenBank (NCBI; Fig. 1A, Appendix II). Similar to the previous studies, the common toad B. bufo (GenBank Accession No: MN122891), which had diverged with B. gargarizans species group in the Late Miocene (Zhan and Fu, 2011; Othman et al., 2022), was chosen as an outgroup. In total 120 haplotypes were analyzed (57 sequences from B. g. miyakonis and 63 sequences from members of the B. gargarizans species group collected at 50 locations: Fig. 1A, Appendix II). The continental populations of B. gargarizans were treated as four geographic groups as follows: western China (W), central China (C), eastern China (E), and northeastern China (NE), following Fu et al. (2005) and Hu et al. (2007). The previous studies suggested that the groups from eastern Taiwan (TE), western Taiwan (TW), and Korean peninsula (K) were all monophyletic (Yu et al., 2014; Borzée et al., 2017), so only one haplotype from each of these groups was used for our analyses.

We then performed phylogenetic analyses by both maximum likelihood (ML) and Bayesian inference (BI) approaches. The most appropriate substitution models were determined using MEGA X based on the corrected Akaike Information Criterion (Akaike, 1974; Kumar et al., 2018) for ML and Bayesian Information Criterion (Schwarz, 1978) for BI. For the ML tree, the HKY+G+I (Hasegawa et al., 1985) model was chosen. MEGA X was used to construct the ML tree. The reliability of the ML tree was assessed by calculating bootstrap probability (BP) (Felsenstein, 1985) with 1,000 replications; all other settings were set to default. To estimate divergence time as well as to infer phylogeny, the BI tree was constructed using BEAST v.2.6.6 (Bouckaert et al., 2019) based on the HKY+G+I model using a coalescent constant model and a relaxed log-normal molecular clock (Drummond et al., 2006). As calibration points, the divergence time of the B. gargarizans species group and B. bufo was set as 14.5±1 Ma, and the time to the most common recent ancestor (TMCRA) of the B. gargarizans species group was set as 9.99±1 Ma (Othman et al., 2022). This tree was initiated with random starting trees, consisting of Markov chain Monte Carlo (MCMC) chains that ran for 50 million generations, with sampling every 1,000 generations. The first 20% of generations were excluded as burn-in. Convergence of the chains to a stationary distribution and the effective sample size (ESS) were assessed visually using Tracer v.1.7.2. (Rambaut et al., 2018). Furthermore, to infer the phylogenetic relationships among haplotypes, median-joining network tree was constructed using NETWORK v.10.2.0.0. (Bandelt et al., 1999; Fluxus Technology Ltd., 2021).

Genetic diversity and population structure

We analyzed the genetic diversity and population structure of B. g. miyakonis samples, using a 615 bp sequence shared by all samples of this subspecies, including the extended region of the ≤348 bp sequence described above. Haplotype diversity (h) and nucleotide diversity (π) were calculated using Arlequin v3.5.2.2 (Excoffier and Lischer, 2010). Populations of B. g. miyakonis were then divided into two sub-populations: the Miyako Island–Kurima Island (MY–KR) and Irabu Island–Shimoji Island (IR–SM) groups. To investigate the genetic structure of these groups, we conducted an analysis of molecular variance (AMOVA; Excoffier et al., 1992) with a permutation test of 1,000 random replications using Arlequin v.3.5.2.2.

Demographic history

To investigate the population demographic history of B. g. miyakonis, the mismatch distribution (Rogers and Harpending, 1992), Tajima's D (Tajima, 1989) and Fu's Fs (Fu, 1997) were computed using Arlequin v3.5.2.2. The statistical significance of these tests was then tested using parametric bootstrapping with 1,000 replicates. In addition, an extended Bayesian skyline plot (EBSP) analysis (Drummond et al., 2005) was used to infer the demographic history. The EBSP was constructed based on an HKY model using a coalescent extended Bayesian skyline model and a relaxed log-normal clock consisting of MCMC chains run for 100 million generations, with sampling every 1,000 generations. The first 20% of generations were excluded as burn-in, and convergence of the chains to a stationary distribution and ESS were assessed visually using Tracer v.1.7.2. The evolution rate of the Miyako population was set to 1.69% per Ma, as inferred from the previously mentioned divergence time estimation. Because reproductive generations overlap in bufonids (e.g., Okuno, 1986), we did not estimate the effective population size in a strict sense. Bayesian skyline plots were produced using plotEBSP.R, an R package, implemented in R v.4.2.0 (R Core Team, 2021).

Table 1.

Genetic variation and statistics of Bufo gargarizans miyakonis (615 bp). In the IR–SM population, the number of samples was small, so the mismatch distribution analysis did not converge. n: sample size, hn: number of haplotypes, h: haplotype diversity value, π: nucleotide diversity value, SD: standard deviation, SSD: sum of square deviations, Rag: Harpending's Raggedness index. Asterisks indicate levels of significance: double asterisks (**) for P≤0.01.

Sequence characteristics

Of the 57 B. g. miyakonis individuals studied, we identified 12 haplotypes using a 347 bp sequence of the mitochondrial control region (shared by all sequences of B. g. miyakonis in the first dataset for phylogenetic inference of the B. gargarizans species group) and 19 haplotypes using the 615 bp sequence (used for population analyses) (Table 1; Appendix I; GenBank Accession No.: LC683658–LC683676). Haplotypes, MYK01, 05, 07, 10, recognized in the 347 bp data were divided into more than one haplotypes in the 615 bp dataset.

Phylogenetic position of Bufo gargarizans miyakonis

Our phylogenetic trees were consistent with previous studies (e.g., Igawa et al., 2006, Zhan and Fu, 2011; Fig. 2A). The result of phylogenetic inference showed that B. g. miyakonis was included in the clade that consisted of many eastern populations of B. gargarizans (Fig. 2A). In contrast, the TE population, which was most geographically proximate, was relatively far from B. g. miyakonis in the tree. The monophyly of B. g. miyakonis was not supported (Fig. 2B). The monophyly of haplotypes, MYK01–10 (hereafter, haplogroup MYK), was supported by high BP and posteri- or probability (PP), but the remaining two haplotypes, MYK11 and MYK12, were located outside of the haplogroup MYK Clade along with a few the eastern populations of B. gargarizans (Fig. 2; E1 from Locality No. 90 and No. 93, and E2 from Locality No. 93, in Fig. 1A). In the network, B. g. miyakonis was not distinct from the continental B. gargarizans group. MYK11, MYK12 were placed between the other haplotypes from the Miyako Islands and those from the continental China (Fig. 2C).

Fig. 2.

Genetic trees of the Bufo gargarizans species group reconstructed from the mt control region (341 bp). Bootstrap probability (BP) for maximum likelihood (ML) and posterior probability (PP) values for Bayesian inference (BI) models are given in order at each node, although only when BP≥70% and PP≥0.80. A: ML trees of the B. gargarizans species group. The clade that includes B. g. miyakonis is indicated in gray. B: ML trees of the clade that includes B. g. miyakonis. Haplotypes of B. g. miyakonis are colored in gray. See the main text for information on how B. g. miyakonis samples were grouped. C: Median-Joining network trees of the clade that includes B. g. miyakonis. Gray circles indicate haplotypes of B. g. miyakonis, white rectangles indicate haplotypes of B. gargarizans or B. bankorensis, and black dots indicate missing haplotypes, respectively. When nucleotide substitution between haplotypes has occurred more than once, the number of times is shown as a number in parentheses.

Divergence time

The average evolution rate of the mitochondrial control region of B. gargarizans species group was estimated to be 1.37%±0.27% (ESS=5,188) per Ma. The divergence time between haplogroup MYK and other haplotypes was 0.75 Ma (95% HPD=0.289–1.02 Ma) in the Middle Pleistocene (Fig. 3), and the evolution rate of haplogroup MYK was 1.69% per Ma. In addition, the divergence times of MYK11 from E1 and MYK12 from E2 were estimated to be 0.62 Ma (95% HPD=0.072–0.749 Ma) and 0.54 Ma (95% HPD=0.022–0.736 Ma), respectively (Fig. 3).

Genetic diversity and population structure

For the 615 bp B. g. miyakonis population dataset, our estimate of h was 0.8972±0.0232, and our estimate of π was 0.004984±0.002919 (Table 2). When comparing the genetic diversity of the sub-populations, IR–SM samples showed approximately half the diversity (as measured by π) as MY–KR samples (Table 1).

Our AMOVA results revealed that there was significant genetic differentiation between sub-populations (Table 2).

Demographic history

For the whole B. g. miyakonis population, Tajima's D was negative (–1.13053) although not statistically significant (p=0.11). Moreover, Fu's Fs was also negative, and it was statistically significant (Table 1). Tajima's D was also not significant for both sub-populations. Fu's Fs was negative and significant only for the MY–KR group. The mismatch distribution of the population was a unimodal curve, with a peak around 3–4 (τ=3.129, θ₀=0.000, θ₁= 17.188, average pairwise difference=3.065). The sum of square variance (SSD) and Harpending's Raggedness index (Rag) were also not significant, suggesting the rapid expansion of this population (Table 1). The mismatch distribution of the MY–KR group was generally similar to that of the population as a whole. The distribution for the IR–SM group was not computed because of insufficient sample size. The EBSP of haplogroup MYK also indicated a rapid population expansion since 0.1 Ma ca (Fig. 4).

Fig. 3.

Bayesian divergence time estimates of the mt control region (341 bp) of the Bufo gargarizans species group. The haplotypes of B. g. miyakonis (the haplogroup MYK and MYK11–12) are colored in dark gray and gray, respectively. The black dot nodes indicate calibration points based on Othman et al. (2022). Gray bars on nodes indicate 95% credible intervals.

Table 2.

AMOVA results. Asterisks indicate levels of significance: double asterisks (**) for P≤0.01.

Fig. 4.

The EBSP of haplogroup MYK. The vertical axis and the horizontal axis are the logarithmic relative population size and the absolute time before the present (Ma), respectively.

Origin and formation history of Bufo gargarizans miyakonis

Phylogenetic inferences suggested that the origin of B. g. miyakonis were from the eastern continental populations. In contrast, the Taiwanese groups, that are geographically closest to the Miyako Islands populations, were not closely related to the Miyako populations genetically. These phylogenetic relationships were consistent with the results of biochemical (Nishioka et al., 1990) and morphological (Matsui, 1984) analyses. Our results also suggested that B. g. miyakonis is not monophyletic in the mtDNA phylogeny. The non-monophyly of B. g. miyakonis might be attributable to an incomplete lineage sorting. In this species group, some cases of what appear to be incomplete lineage sorting have been reported (Yu et al., 2014; Othman et al., 2022). To clarify the relationship between B. gargarizans species group and B. g. miyakonis, more markers such as nuDNA loci should be used in future studies.

The divergence time between haplogroup MYK and the continental haplotypes was estimated to be 0.75 Ma, sometime in the Middle Pleistocene. The divergence times of MYK11 and MYK 12 from respective continental counterparts were estimated to be 0.54 and 0.62 Ma approximately, respectively. These estimates were slightly younger than 0.9 Ma, the estimate given by Igawa et al. (2006) who used the Beijing population as a representative of B. gargarizans. Since our dataset included closer lineages (E1 and E2) to Miyako populations, which collected in sites closer to the Miyako Islands than Beijing (Locality No. 90 and No. 93 in Fig. 1A), our estimates were slightly younger.

The distribution of close populations suggests that the ancestral B. g. miyakonis population migrated over a land bridge connecting eastern Eurasia and the Miyako Islands. A land bridge between continental Eurasia and the Miyako Islands has been suggested by fossil evidence, including many Palearctic species, from the Pinza-Abu cave on Miyako Island (Hasegawa, 1985a; Oshiro and Nohara, 2000). Almost all of these Palearctic species such as Capreolus miyakoensis and Rattus miyakoensis are considered to be closely related to other species present in northeastern Eurasia (Hasegawa, 1985a, 1985b; Zhen and Hasegawa, 1985; Kawaguchi et al., 2009). In the Early Pleistocene, the Miyako Islands were last connected to neighboring islands (Kizaki and Oshiro, 1977, 1980; Ota, 1998; Kimura, 2003). At 0.12–0.2 Ma in the Middle Pleistocene, the area between the continent and the Miyako Islands was likely shallow and there was a temporal land bridge at the ebb (Kimura, 1990; Kimura et al., 2001). Kimura et al. (2001) suggested a probability that some of the above Palearctic species might have migrated to the Miyako Islands across this temporal land bridge during the Middle Pleistocene. Considering that the toads depend on fresh water, it is unlikely that B. g. miyakonis could migrate through a shallow sea like mammal species (e.g., C. miyakoensis) could. However, our results suggest that the migration route of B. g. miyakonis was similar to other Palearctic species. We considered that this land bridge had been a corridor for the ancestor of B. g. miyakonis. Although recent phylogenetic analyses have suggested that several frog species inhabiting the Ryukyu Archipelago dispersed overseas (Komaki et al., 2017; Yu et al., 2020), the oversea dispersal of B. g. miyakonis is unlikely to happen, because its expected migration route is perpendicular to the present-day and past Kuroshio Current (Gallagher et al., 2015).

Stratigraphic studies have suggested that the Miyako Islands submerged (completely or partly) during the Early–Middle Pleistocene, 0.41–0.95 Ma (Iryu et al., 2006). Since this estimate is more recent than the divergence time between B. g. miyakonis and the continental B. gargarizans group, 0.54–0.75 Ma, the migration event of B. g. miyakonis may have occurred later than the genetic divergence. It is consistent with the inference that the ancestral population of B. g. miyakonis had migrated to the Miyako Islands at 0.12–0.2 Ma, when the land bridge had connected the continent and the Miyako Islands (see above; Kimura et al., 2001). Furthermore, the scale of crustal movements during the Late Pleistocene and the uniqueness of the fossils found on individual islands of the Ryukyu Archipelago both provide reasons to be skeptical of the existence of a Late Pleistocene land bridge in the Ryukyu Archipelago (Hasegawa, 2012; Furukawa and Fujitani, 2014; Yokoyama et al., 2018). Our results suggested that the ancestral population of B. g. miyakonis migrated to the Miyako Islands at the Middle Pleistocene, consistent with the absence of a land bridge in the Late Pleistocene. The curious distribution of the B. g. miyakonis may be affected by the formation history of the Miyako Islands. Future research should focus on seeking geographic or stratigraphic evidence for this land bridge.

The migration history of B. g. miyakonis is notably different from that inferred for other species from the Miyako Islands (land snails: Kameda et al., 2007; freshwater crab: Segawa, 2011; Asian keelback snake: Kaito and Toda, 2016; skinks: Kurita and Toda, 2017; narrow-mouthed frog: Tominaga et al., 2019). Moreover, on the Miyako Islands, almost all animal species that originated from northeastern Eurasia are considered extinct (Hasegawa, 2012), and there are few modern species that have their closest relatives on the continent like B. g. miyakonis (e.g., Miyako grass lizard: Lin et al., 2002; Ota et al., 2002). Therefore, our results may be important for understanding the formation process of present-day terrestrial fauna of the Miyako Islands.

Population genetics

Our population genetic analyses suggested rapid demographic expansion of B. g. miyakonis. For the B. g. miyakonis population, we measured a high value for h (0.8972±0.0232) and a low value for π (0.004984±0.002919) relative to other members of the B. gargarizans species group. Hu et al. (2007) gave h and π values for the continental populations of the B. gargarizans species group: the mean h was 0.947±0.006 with a range from 0.000–0.855± 0.085, and the mean π was 0.04739±0.00215 with a range from 0.0000–0.04438±0.00524 (Hu et al., 2007). Inferred from the relationship between a relatively high h and a low π values (Grant and Bowen, 1998; Avise, 2000), it was suggested that B. g. miyakonis was initially genetically uniform and subsequently accumulated mutations due to a rapid expansion of the population size. Our mismatch distribution analysis and Fu's Fs test results also supported this demographic trend. The EBSP also suggested a rapid expansion from 0.1 Ma until the present. Although our estimate suggests that this expansion started before the last glacial maximum (LGM), it might also have been caused by increases in land area following glacial regression. Because of the high diversity of fossil fauna and the existence of large mammals in the Late Pleistocene, the environment of the Miyako Islands is thought to have been significantly larger in both land area and carrying capacity before the LGM relative to the present (Ota and Takahashi, 2008). After the LGM, the land area and carrying capacity decreased, and large animals became extinct (Ota and Takahashi, 2008). However, it was estimated that the expansion of B. g. miyakonis population was maintained after the LGM.

We also found evidence of moderate genetic differentiation of B. g. miyakonis sub-populations. We found many unique haplotypes in the MY–KR sample set, which had a relatively large sample size. On the other hand, we found no unique haplotype from the IR–SM dataset, partially because of the sample size was smaller. Nevertheless, we found evidence of differentiation among these sub-populations. Interestingly, no such genetic differentiation has been reported for other species present on the Miyako Islands (smooth skinks: Koizumi et al., 2014; anchialine atyid shrimp: Weese et al., 2012, 2013). All collection sites used in this study were part of a continuous landmass at the LGM (Japan Association for Quaternary Research, 1987). After the LGM, geographic isolation due to the transgression of the sea level likely contributed to the genetic differentiation between subpopulations.