Phylogenetic Systematics of the Water Toad (Bufo stejnegeri) Elucidates the Evolution of Semi-aquatic Toad Ecology and Pleistocene Glacial Refugia

Introduction

The discipline of phylogenetic systematics, and its goal to understand character evolution and speciation, has shaped our understanding of the patterns and processes of evolution (Wiley and Lieberman, 2011). Inter- and intra-specific phylogenies have been applied to elucidate several aspects of species evolution, such as adaptation, biogeography, and diversification. What begins as a simple study to understand the evolutionary genetics of a species can provide insight into broader evolutionary questions. Such is the case with our study of the water toad (Bufo stejnegeri).

The water toad (B. stejnegeri) is endemic to Northeastern Asia (South Korea, North Korea, and China). As the common name implies, B. stejnegeri has unique ecology for a toad by being semi-aquatic and breeding in lotic environments. Between November and January, males and females enter streams, immediately begin amplexus, and stay in amplexus for 3–6 months in frozen-over streams until spring spawning (March–April) (Kang and Yoon, 1975; Lee and Park, 2009; pers. obs.). Besides these basics, much of the biology of B. stejnegeri remains a mystery. Studying the genetics of B. stejnegeri provides a unique opportunity to understand the evolution of semi-aquatic toad ecology.

A majority of true toads in the family Bufonidae are terrestrial and characterized by dry, leathery skin, short legs, and parotoid glands (Pough et al., 2004). There are several exceptions to this (see Van Bocxlaer et al., 2010), with one being semi-aquatic stream species in the genus Ansonia, and three species in the genus Bufo (B. stejnegeri, B. andrewsi, B. torrenticola). All Ansonia species have similar adult (slender body, long legs) and tadpole (streamlined, sucker-type oral disk) morphologies, which are clear adaptations of this genus to live in fast-flowing, rocky streams (Inger, 1960). The other stream-dwelling toad species, however, are members of the largely terrestrial Bufo genus. It is unclear whether a semi-aquatic, stream lifestyle has independently evolved in these three species or was influenced by their phylogenetic relationships. Inferring the phylogenetic placement of B. stejnegeri and the relationship between these three species can provide insight into the evolution of a semi-aquatic lifestyle in Bufo.

This study also allows us to investigate the impact of Pleistocene glacial cycles on the biodiversity of Northeast Asia. The Pleistocene was a period characterized by climatic oscillations (glacial-interglacial cycles) that impacted the living organisms on Earth (Hewitt, 1996, 2004). In extreme cases, the glacial cycles caused species extinction (e.g., Picea critchfieldii; Jackson and Weng, 1999), while in other cases species ranges contracted to glacial refugia (e.g., Chorthippus parallelus, Erinaceus europaeus/concolor, Ursus arctos; Hewitt, 1999). In the latter case, the range of a species (in the Northern hemisphere) would retreat southward during a glacial period followed by expanding northward during an interglacial period. Genetically, what results is a pattern of southern richness and northern purity (Hewitt, 2004). One region with a controversial Pleistocene history is Northeast Asia (northeastern China, the Korean Peninsula, Japan, and Russian Far East). At least six glacial advances occurred during the Pleistocene (Hewitt, 1996, 1999; Avise, 2000), and it is unclear how biodiversity in the area responded.

Two hypotheses have been proposed for Northeast Asia–a multiple refugia hypothesis and a Korea refugium hypothesis. The multiple refugia hypothesis suggests that ≥2 refugia existed during the Last Glacial Maximum, preserving and promoting species diversity through allopatry (Qian and Ricklefs, 2000). Phylogeographic studies on several species support this multiple refugia hypothesis: Manchurian walnut (Juglans mandshurica, Bai et al., 2010), Asian spruce (Picea jezoensis, Aizawa et al., 2007), Siberian chipmunk (Tamias sibericus, Lee et al., 2008), Korean field mouse (Apodemus peninsulae, Serizawa et al., 2002), and Japanese treefrog (Dryophytes japonicus, Dufresnes et al., 2016). In contrast, the Korea refugium hypothesis posits that glaciation pushed temperate forests below 30°N, resulting in a single refugium on the Korean Peninsula (Kong, 2000; Harrison et al., 2001). Studies supporting this hypothesis were conducted on the Chinese black-spotted frog (Rana nigromaculata, Zhang et al., 2008) and short-tailed pit viper (Gloydius brevicaudus, Ding et al., 2011). A phylogeographic study of B. stejnegeri, restricted to Northeast Asia, provides another opportunity to evaluate these hypotheses and understand how glacial cycles impacted the biodiversity of Northeast Asia.

In this study, we use a suite of genetic analyses to infer the phylogenetic placement, phylogeography, and population genetics of B. stejnegeri. We apply these results to understand the evolution of a semi-aquatic lifestyle in Bufo and the impact of Pleistocene glacial cycles on the biodiversity of Northeast Asia.

Materials and Methods

Sampling

We collected genetic data from 43 individuals across the range of B. stejnegeri in northeastern China and South Korea. Although B. stejnegeri is distributed in North Korea, we were unable to obtain samples from the region. Tissue samples (liver or toe clips) for genetic work were preserved in 95% ethanol.

Laboratory Protocols

We extracted genomic DNA from tissues using a standard salt extraction protocol (Sambrook and Russell, 2001). A total of 12 loci were targeted for PCR and DNA sequencing (four mitochondrial DNA [mtDNA] loci, eight nuclear DNA [nuDNA] loci). The four mtDNA loci are 12S ribosomal RNA + 16S ribosomal RNA (12S-16S), cytochrome b (cytb), NADH dehydrogenase subunit 2 (ND2), and NADH dehydrogenase subunit 3 and associated tRNAs (ND3), while the eight nuDNA loci are 28S ribosomal RNA (28S), C-X-C chemokine receptor type 4 (CXCR-4), histone 3A (H3A), pro-opiomelanocortin (POMC), recombination activating gene 1 (RAG-1), rhodopsin (RHO), activity-dependent neuroprotector homeobox (ADNP), and neurotrophin 3 (NTF3). The primers used for each locus are listed in Table S1. All PCR products were purified and sequenced in both directions using the PCR primers on an ABI3730XL machine (Applied Biosystems, USA) at the National Instrument Center for Environmental Management (Seoul National University). Owing to the length of PCR products, additional internal primers were used to sequence 12S-16S (12S-7, 2571H) and cytb (Bufo3Inner1) (Table S1).

Phylogenetic Analyses

Sequences for each locus of each individual were assembled and checked by eye in Geneious 11.1.5 (Biomatters, Auckland, New Zealand). Consensus sequences for each locus were aligned using MUSCLE (Edgar, 2004). These multiple sequences alignments were checked manually, and sequences of protein coding genes were translated into amino acids to check for erroneous stop codons.

We performed phylogenetic analyses to (1) determine the phylogenetic placement of B. stejnegeri amongst other Bufo species, and (2) uncover phylogeographic patterns for B. stejnegeri. Different datasets and taxon sampling were used for each goal. For the phylogenetic placement of B. stejnegeri, three loci were used (12S-16S, CXCR-4, POMC), chosen based on the available data on GenBank. The taxon sampling consisted of 19 individuals in the genus Bufo (one B. stejnegeri, 18 individuals of 13 Bufo species from GenBank). Four extant Bufo taxa were not included as their data were not available on GenBank. Three taxa were included as outgroups (Incilius alvarius, Rhinella marina, Duttaphrynus himalayanus). We performed separate analyses for the individual loci, as well as for a concatenated dataset. For the phylogeography analyses, 10 loci were used (12S-16S, cytb, ND2, ND3, ADNP, CXCR-4, H3A, NTF3, POMC, RAG-1). Taxon sampling for these analyses was 46 B. stejnegeri across its geographic range [43 new, three GenBank (Jeong et al., 2013; Dong and Yang, 2016)] and two B. gargarizans outgroups. Collection localities are plotted on the map in Figure 1, grouped geographically by mountain range (Qian, Hwaak, Taebaek, Chiak, Jiri). Analyses were performed on two datasets: concatenated mtDNA (12S-16S, cytb, ND2, ND3) and concatenated data of all 10 loci. Both datasets included the complete taxon sampling of 46 B. stejnegeri and two B. gargarizans outgroups.

Figure 1. Maximum likelihood phylogeny of Bufo stejnegeri inferred from a concatenated dataset of 10 loci (12S-16S, cytb, ND2, ND3, ADNP, CXCR-4, H3A, NTF3, POMC, RAG-1). Thicker branches indicate a bootstrap support value ≥ 85. The shapes next to each individual correspond to the geographic locations on the map to the right. Geographic locations have been grouped into five populations based on the mountain range. The geographic range of B. stegnegeri is indicated by the light-green shading, based on the Map of Life (Jetz et al., 2012). Although this species is expected to occur in North Korea, no data are available. Darker lines on the map (solid and dotted) indicate the location of the major mountain ranges (adapted from Park and Son, 2008).

We performed maximum likelihood (ML) analyses in RAxML v.8.2.4 (Stamatakis, 2014) and Bayesian (BI) analyses in MrBayes v.3.2.5 (Ronquist et al., 2012). Datasets were first partitioned by locus and codon (for protein-coding genes), whenever applicable, followed by finding the best partitioning scheme in PartitionFinder 2.1.1 (Lanfear et al., 2016) using the greedy algorithm (Lanfear et al., 2012). Each ML analysis was performed using the combined rapid bootstrap and search for the best-scoring tree, 1,000 bootstrap replicates, and the GTR+G model of nucleotide substitution for each partition. For BI analyses, the substitution model for each partition was selected in MrBayes using the reversible jumping model choice, with both rate variation and invariable sites. Two BI searches with random starting trees were run and compared using four Markov Chain Monte Carlo (MCMC chains), 10 million generations, and sampled every 1,000 generations. Stationarity was evaluated by ensuring the average standard deviation of split frequencies was low (<0.01) and checking to see the plot log-likelihood scores had plateaued in Tracer v.1.7 (Rambaut et al., 2018). Burn-in data were discarded before building the consensus tree.

We did a post hoc analysis of phylogenetic informativeness for each of the 10 molecular markers using the online program PhyDesign (http://phydesign.townsend.yale.edu; Townsend, 2007). The ML tree was converted to an ultrametric tree using FigTree v.1.3.1 (Rambaut, 2010), and HyPhy (Pond et al., 2005) was used to calculate the substitution rates.

Ancestral State Reconstruction

We performed ancestral state reconstruction analyses to understand the evolution of a semi-aquatic lifestyle in Bufo using the packages phytools (Revell, 2012) and ape (Paradis and Schliep, 2018) in R 3.5.3 (R Core Team, 2019). The ML tree resulting from the Bufo phylogeny was used, which was topologically identical to the BI tree. Lifestyle was coded for all 17 Bufo as a discrete character–semi-aquatic for B. stejnegeri, B. torrenticola, and B. andrewsi; and terrestrial for all other species/individuals. The equal rates [ER] model of evolution was implemented in the “ace” function, which implicitly assumes a flat prior probability distribution at the global root of the tree. We used stochastic character mapping (Huelsenbeck et al., 2003) to sample character histories from their posterior probability distribution and 1,000 simulation replicates, and “countSimmap” to count the number of character changes in the simulated trees. We use the “densityMap” method to visualize the posterior probability of each state at all points along the tree.

Divergence Dating Analyses

A 114-sample, 7-locus (12S-16S, 28S, CXCR-4, H3A, POMC, RAG-1, RHO) Bufonidae dataset (primarily from Frost et al., 2006; Pramuk, 2006; Pramuk et al., 2008) was used for the divergence dating analyses. These loci were chosen based on the available data in GenBank. The amount of missing data varied for each taxon, with some taxa represented by data from a single locus. As missing data have been shown to have minor impacts on the accuracy of divergence dating (Zheng and Wiens, 2015), we proceeded with the analysis. Analyses were performed in BEAST 2.5.2 (Bouckaert et al., 2014) and run on the CIPRES Science Gateway (Miller et al., 2010). Analyses were run for 100 million generations, sampling every 10 thousand generations, with a Yule model of speciation as a tree prior. Although using a relaxed clock can accommodate the more restrictive case of a strict clock (Drummond et al., 2006), we check for the best clock model by comparing marginal likelihoods under an uncorrelated relaxed lognormal clock and a strict clock. We estimate marginal likelihoods using path sampling/stepping-stone sampling with 1 path, 100 million generations, sampling every 10 thousand generations. The marginal likelihoods of each model were compared with Bayes factors, interpreting the results following Kass and Raftery (1995).

The substitution models determined in PartitionFinder were GTR + G + I for 12S-16S, 28S, CXCR-4, H3A, POMC, and RAG-1, and HKY + G for RHO. Three calibration points were used, following Pramuk et al. (2008) and Portik and Papenfuss (2015): (1) oldest bufonid (Baez and Gasparini, 1979) exponential distribution, mean = 20, log(SD) = 1, offset = 57, producing a credibility interval of 58.03–116.9; (2) Central-North American split based on the fossil of B. praevis (Tihen, 1951), lognormal distribution, mean in real space = 10, log(SD) = 1, offset = 19, producing a credibility interval of 20.17–50.42; and (3) Origin of Rhinella marina group (Sanchíz, 1998), lognormal distribution, mean in real space = 10, log(SD) = 1, offset = 10, producing a credibility interval of 11.17–41.42. We used Tracer 1.7.1 (Rambaut et al., 2018) to determine if runs had reached stationarity, and a burn-in of 25% was discarded.

Population Genetic Analyses

All population genetics analyses were performed on a concatenated ND2-ND3 dataset, as this dataset had the least missing data. A total 42 of 46 individuals were retained for these analyses. Individuals were grouped into populations based on mountain range: Qian (3 individuals), Hwaak (16 individuals), Taebaek (10 individuals), and Chiak (13 individuals). Unfortunately, the single sample from Jiri (downloaded from GenBank) did not have available sequence for ND2-ND3 and therefore was not included. We attempted to analyse datasets including this sample, but theses analyses were not informative since the sequences were either short (568 bp of 16S) or had low taxon sampling (9 individuals for cytb).

We used PopART (Leigh and Bryant, 2015) to infer a haplotype network using the median-joining network algorithm (Bandelt et al., 1999), setting epsilon to zero. Additionally, with PopART, we performed an AMOVA to compare the variation among and within populations and calculated Φ_ST. As we are dealing with mtDNA sequences, we chose Φ_ST as the appropriate fixation index to calculate. DnaSP 6.11.01 (Librado and Rozas, 2009) was used to calculate genetic diversity indices for each population (number of haplotypes [Nh], haplotype diversity [h], and nucleotide diversity [π]).

Inferring Demographic History

We inferred the effective population size of B. stejnegeri using BEAST 2.5.2, using the same ND2-ND3 dataset as the population genetic analyses. Substitution models were set based on results from PartitionFinder. The tree prior was set as coalescent Bayesian skyline tree, while the prior for the root height was based on the result of the divergence dating analysis. MCMC chains were run for 500 million generations, sampling every 10 thousand generations. In Tracer 1.7.1, the plot of the posterior score was viewed to determine burn-in, followed by construction of the Bayesian skyline plot.

Results

Phylogenetic Analyses

All new sequences were deposited in GenBank (Table S2). The ML phylogenetic tree inferring the phylogenetic placement of B. stejnegeri within the genus Bufo is shown in Figure 2. The BI phylogeny only differed in the branch lengths and support values of nodes, with BI posterior probabilities generally higher. Bufo stejnegeri is inferred to be the sister species of a group including B. japonicus, and B. torrenticola (ML bootstrap: 44, BI posterior probability: 0.96). All analyses of the individual locus datasets had lower resolution compared to the concatenated dataset (not shown), but B. stejnegeri always clustered with the same Bufo species listed above.

Figure 2. (A) The phylogenetic placement of Bufo stejnegeri within the genus Bufo and ancestral state reconstruction of a semi-aquatic lifestyle. This phylogeny is from the partitioned maximum likelihood (ML) phylogenetic analysis performed on a concatenated dataset of three loci (12S-16S, CXCR-4, POMC). Numbers at nodes represent (ML bootstrap/Bayesian posterior probability). The Bayesian inference phylogeny only differed in branch lengths and support values of nodes. Colors of the branches represent the results from the ancestral state reconstruction analysis, while circles at the nodes represent inferred ancestor states nodes present the lifestyle state of each species. For branches, red represents a semi-aquatic lifestyle, blue represents a terrestrial lifestyle, while magenta-purple represent intermediate posterior probabilities. Plot probability distributions of the number of reconstructed changes across all simulated histories (B) terrestrial → aquatic and (C) aquatic → terrestrial.

The ML tree used for understanding B. stejnegeri phylogeography is shown in Figure 1, and the BI tree in Figure S1. Both phylogenies were largely concordant, with minor differences in support values and the presence of more polytomies in the BI tree. Results from the concatenated mtDNA dataset (not shown) were also largely concordant to the 10-locus dataset. Bufo stejnegeri is monophyletic and exhibits genetic variation. It is noteworthy that the single Jiri sample was the most basal on the tree, while the three Qian samples were monophyletic and in a nested position. However, it is unclear if these are true patterns or an artifact due to small samples sizes. Of the three remaining populations (Hwaak, Taebaek, and Chiak) there is no clear pattern relating this genetic variation to geography–individuals from each population are relatively evenly dispersed across the phylogeny.

For the phylogenetic informativeness analysis. A majority of the markers reached maximum phylogenetic informativeness closer to the base of the tree (Figure S2). The only exception was H3A, which showed to be maximally informative closer to the tips of the tree.

Ancestral State Reconstruction Analyses

Figure 2A depicts the results of the ancestral state reconstruction simulations on the ML tree The posterior probability of being in each state is illustrated by the colors of branches, with red being a 1.0 posterior probability of being semi-aquatic, blue being a zero posterior probability of being terrestrial, and purple-magenta representing intermediate posterior probabilities. For clarity, the true states of the tips were included. Figures 2B,C plot probability distributions of the number of reconstructed terrestrial-aquatic and aquatic-terrestrial changes across all simulated histories; our simulations most frequently inferred three terrestrial-aquatic and zero aquatic-terrestrial changes. The semi-aquatic ecology of B. stejnegeri, B. torrenticola, and B. andrewsi are inferred to have evolved independently of each other.

Divergence Dating Analyses

The marginal likelihoods of the analysis were compared under an uncorrelated relaxed lognormal (−68379.725) and strict (−68858.507) molecular clock. With 2log_e(B₁₀) = 12.46 Kass and Raftery (1995) interpret this to be very strong evidence to support the use of an uncorrelated relaxed lognormal molecular clock. All subsequent results are from the analysis using the relaxed clock model.

The results from the divergence dating analysis are found in Figure S3. Effective sample sizes of the majority of parameters were >200, providing us with confidence that analyses reached stationarity. The two divergence dates of interest for our study were (mean, 95% confidence interval): between the two Korean Bufo species (B. stejnegeri and B. gargarizans) (7.3 Ma [mega-annum], 4.9–10.2 Ma), and the origin of Bufo stejnegeri (4.3 Ma, 2.7–6.2 Ma). The relationships within the genus Bufo were largely consistent to the RAxML analysis, but with stronger posterior probability support. The two differences in topology were the placement of B. torrenticola (sister species to B. japonicus [RAxML] vs. sister species to B. j. formosus [BEAST]) and B. bufo (sister species to B. verrucosissimus [RAxML] vs. sister species to all other Bufo species [BEAST]).

Population Genetics Analyses

All population genetics analyses excluded the Jiri population (see explanation in the “Materials and Methods” section). The population genetics summary statistics for each of the remaining populations (separated by mountain range) are in Table 1. The three different population genetics analyses yielded similar results–populations designated a priori based on mountain range are not distinct and likely retain ancestral polymorphism and/or experience gene flow. First, the median-joining haplotype network recovers genetic clusters, but these clusters do not correspond to populations (Figure 3). Second, the AMOVA found the majority of variation to be within populations (81.5%; d.f = 38; sum of squares = 2394.265; σ² = 63.007) as opposed to among populations (18.5%; d.f = 3; sum of squares = 789.854; σ² = 14.306). Third, the fixation index value (Φ_ST = 0.185; p < 0.001), indicates there is significant genetic variation between populations, but the variation is not large.

Table 1. Population genetics summary statistics for populations of Bufo stejnegeri designated based on mountain range, based on a concatenated ND2-ND3 dataset.

Figure 3. Median-joining network of Bufo stejnegeri haplotypes based on a concatenated ND2-ND3 dataset constructed in PopART. Haplotypes are represented by circles, whose sizes are proportional to the number of individuals. Different colors represent different mountain ranges, and potentially different populations. The color scheme is the same as in Figure 1. Hatch marks correspond to nucleotide substitutions, while missing haplotypes are indicated by small black circles.

Demographic History

Figure 4 shows the Bayesian skyline graph inferring the effective population size of B. stejnegeri over time. Effective population size began to decrease during the Pliocene (~3.5 Ma) and continued to decline through the Pleistocene, reaching the lowest point at ~0.4 Ma. After this point, B. stejnegeri experienced a population expansion, where the population size recovered to or exceeded Pliocene levels.

Figure 4. Bayesian skyline plot inferring the effective population size of Bufo stejnegeri.

Discussion

Bufo stejnegeri is a toad species endemic to Northeastern Asia (China, South Korea, and North Korea) with a unique ecology in that it is semi-aquatic and breeds in lotic environments. Males and females enter streams in autumn and begin amplexus, and remain in amplexus during winter in frozen streams until spring spawning in March–April (Kang and Yoon, 1975; Lee and Park, 2009; pers. obs.). We undertook this study to understand both the higher-level and population-level genetic relationships of B. stejnegeri to gain insight into the evolution of B. stejnegeri's distinctive ecology and the impact of Pleistocene glacial cycles on the biodiversity of Northeast Asia.

Evolution of Semi-aquatic, Lotic Lifestyle

In the Bufo phylogeny, B. stejnegeri is relatively closely related to B. torrenticola and B. andrewsi, two species that have some similarities in ecology and reproductive biology. Bufo species are largely terrestrial and breed in lentic environments (Duellman and Trueb, 1994), but B. stejnegeri, B. torrenticola, and B. andrewsi are semi-aquatic and breed in lotic environments (Matsui, 1980; Tsuji and Kawamichi, 1996; Lee and Park, 2009). These three species have similar morphological characters that are likely adaptations to life in mountainous streams–indistinct or small tympanums, long and slender limbs, and less divergent short parotoids (Matsui, 1980). Additionally, the skin of B. stejnegeri and B. torrenticola becomes smooth during the breeding season (Tsuji and Kawamichi, 1998a; pers. obs.), and spawning pairs mate submerged underwater (Tsuji and Kawamichi, 1998a,b; pers. obs.).

These similar characters could be a result of convergent evolution to lotic environments, due to phylogenetic relatedness, or a combination of both. Our ancestral state reconstruction analyses infer independent evolution of this unique ecology in the three Bufo species (Figure 2). Upon closer investigation of the features in these three species, we find differences that support independent evolution of a semi-aquatic lifestyle. In terms of amplexus and egg-laying strategy, B. torrenticola and B. andrewsi exhibit behavior more typical Bufonidae species; amplexus occurs immediately before laying eggs, and eggs are laid in pools along the streams where the water flow is slower (Liu and Hiu, 1961; Matsui, 1976; Tsuji and Kawamichi, 1996, 1998a). In contrast, B. stejnegeri has extended amplexus that lasts 3–6 months, and eggs are primarily laid in fast-flowing streams under rocks (Matsui, 1980; Lee and Park, 2009; pers. obs.). These additional characters support the inference that B. stejnegeri has independently evolved a semi-aquatic, lotic-breeding lifestyle.

Phylogeography and Population Genetics of B. stejnegeri

Bufo stejnegeri lives in the riparian zone of mountainous forests (200–700 asl), coming down to streams and rivers to breed (Matsui and Zhao, 2004). For phylogeographic and population genetics analyses, individuals were grouped into a priori populations to test whether population structure corresponds to mountain range (Qian, Hwaak, Taebaek, Chiak, Jiri). Population structure may result from either geography acting as barriers to gene flow or B. stejnegeri exhibiting site fidelity. Site fidelity (returning to natal pond to breed) has been observed in several Bufonid species (e.g., Heusser, 1969; Sinsch, 2010).

For B. stejnegeri, phylogeographic structure does not correspond to mountain range, as individuals from Hwaak, Taebaek, and Chiak populations are interspersed across the phylogeny (Figure 1). However, when the two populations of B. stejnegeri at the northernmost (Qian) and southernmost (Jiri) extremes of the distribution are considered, a preliminary pattern appears. The southern Jiri sample is the most ancestral, while the northern Qian samples cluster in a nested position. This needs to be tested with additional sampling from the extremes of the distribution, but if the pattern holds, it would be similar to the “southern richness and northern purity” pattern seen in the northern hemisphere resulting from glacial-interglacial cycles of the Quaternary–contraction southward followed by northward expansion (Hewitt, 2004; Clark et al., 2009). If this is the case for B. stejnegeri, the population would have retreated to a single glacial refugium in the Jiri Mountain Range, followed by a northward population expansion along associated mountain ranges, supporting the Korea refugium hypothesis. Further studies are needed to make a conclusion, with denser sampling from the extremes of the distribution, and (based on our analysis of phylogenetic informativeness) more quickly-evolving markers to resolve shallower nodes in the phylogeny.

The geology of the Korean Peninsula and northeastern China would have facilitated population expansion because all of the mountain ranges in the region are connected and part of the larger Baekdudaegan Mountain Range. The Baekdudaegan is often referred to as the “spine” of the Korean Peninsula, running 1,625 km from Jiri Mountain in southern South Korea up to Baekdu Mountain at the China-North Korea border. The China portion of the mountain range is called the Changbai Mountains, extending across Heilongjang, Jilin, and Liaoning Provinces (where the Qian population is located). The formation of these mountain ranges began with volcanism in the Oligocene (28.4 Ma) and lasted until the Pleistocene (1.2 Ma), with most activity during the Miocene and Pliocene (Wang et al., 1989; Wan, 2012). Results from divergence dating and estimation of demographic history are consistent with the proposed scenario.

The date inferred for the origin of B. stejnegeri (4.3 Ma, 2.7–6.2 Ma) places the origin of B. stejnegeri after the majority of the Baekdudaegan had already formed and before glacial cycles of the Pleistocene. Consistent with this scenario, the Bayesian skyline plot infers that B. stejnegeri experienced a decrease in effective population size during the Pliocene-Pleistocene, followed by an expansion (Figure 4). We use these results to infer the demographic history of B. stejnegeri, although we recognize that effective population size is almost always smaller than census population size, and changes in effective population size do not necessarily indicate a corresponding change in census population size. Additionally, this scenario should be viewed as a hypothesis, as the sample sizes for Qian (n = 3) and Jiri (n = 1) were small in our study. Additional samples from these two regions are needed to verify whether there is a true pattern of southern richness and northern purity.

Korean Bufo

Two Bufo species are found on the Korean Peninsula: B. stejnegeri and B. gargarizans. Two previous studies based on mtDNA (Lee and Park, 1992) and allozymes (Yang et al., 2000) looked at intra- and inter-species relationships. Studying the divergence time between the two species, they inferred an age of 3.6 or 2.8 Ma, respectively. This corresponds to a Pliocene divergence. Our results imply a much older divergence date (7.3 Ma, 4.9–10.2 Ma), pushing the divergence back into the Late Miocene-Early Pliocene. The disparity between the results is likely due to the difference in analytical methods; Lee and Park (1992) and Yang et al. (2000) used a fixed divergence rate (2%/million years [mtDNA] and 20%/million years [allozymes]), while we used a relaxed clock and several fossil calibration points. We believe our methods, taking into account uncertainty in historical events and potential differences in evolutionary rates of lineages, provide a more realistic representation of the biological process and therefore a better estimation of the true divergence time.

Our results push back the timing of the split between Korean Bufo species closer to the age of another Korean frog, Bombina orientalis (14.80 Ma, 6.50–23.52 Ma) and its four subclades across Northeast Asia (mean ages 3.12–9.20 Ma) (Fong et al., 2016). These studies of frog species in Korea point toward the Late Miocene-Early Pliocene as being an important period of diversification. This was a period inferred to be climatically warm and stable (Liu and Leopold, 1994; Kong, 1996), but a time of major geologic events such as forming the major mountain ranges in the area (Wang et al., 1989; Park and Son, 2005, 2008; Wan, 2012) and opening of the East Sea/Sea of Japan (Otofuji et al., 1985, 2003). Comparative phylogeographic studies testing for simultaneous divergence of animal and plant species (Hickerson et al., 2006) will be helpful to identify the key geologic events influencing diversification on the Korean Peninsula.

Conclusion

In our study, genetic data of B. stejnegeri helped us understand the inter- and intra-specific relationships of this species. From a broader evolutionary view, these data suggest an independent evolution of a semi-aquatic, lotic lifestyle in Bufo, as well as support the presence of a single refugium in Northeast Asia during the Last Glacial Maximum.

Data Availability Statement

The datasets generated for this study can be found in GenBank with the following accession numbers: MK031411–MK031654.

Ethics Statement

The animal study was reviewed and approved by Seoul National University Ethics Committee.

Author Contributions

JF, B-TY, P-PL, BW, and M-SM designed the study and wrote the manuscript. JF, B-TY, and M-SM collected samples. JF and B-TY analyzed the data.

Funding

This research was supported by grants from the National Research Foundation of Korea (2012K1A2B1A03000496) funded by the Government of the Republic of Korea (MSIP), and the Research Institute of Basic Sciences, Seoul National University (to BW), and from the National Natural Science Foundation of China (31211140342) (to P-PL). We thank Profs. Wenbo Liao and Jinzhong Fu for information on the reproductive biology of B. andrewsi. We would like to thank members of the Min and Waldman labs for help with fieldwork, and thank Adam Leache for help with divergence dating analyses.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2019.00523/full#supplementary-material

Figure S1. Bayesian inference phylogeny of Bufo stejnegeri inferred from a concatenated dataset of 10 loci (12S-16S, cytb, ND2, ND3, ADNP, CXCR-4, H3A, NTF3, POMC, RAG-1). Thicker branches indicate a bootstrap support value ≥0.95. The shapes next to each individual correspond to the geographic locations on the map from Figure 1.

Figure S2. Calculation of phylogenetic informativeness of the 10 markers using PhyDesign.

Figure S3. Phylogeny and divergence dates inferred using BEAST. Numbers at nodes represent the mean divergence date (in Ma [mega-annum]), while bars at nodes represent the 95% confidence interval of the divergence date. Nodes with high support (posterior probability ≥ 0.95) are denoted with black, bold branches, while nodes with low support (posterior probability <0.95) are denoted with red, thin branches.

Table S1. Primers used in this study.

Table S2. Geographic localities and GenBank numbers of new sequences used in this study.

References