Phylogeographic and population insights of the Asian common toad (Bufo gargarizans) in Korea and China: population isolation and expansions as response to the ice ages

Field work

Bufo gargarizans is usually found between 20 and 800 m a.s.l., in coniferous, mixed and deciduous forests, as well as grasslands, in humid but not saturated habitat (Kuzmin et al., 2004). Due to permit restriction from the Ministry of Environment (permit numbers: Yesan-2016-10; Boeun-2016-01; Gangwha-2016-01; Jeonju-2016-01; Hampyeong-2016-37; Daegu-2016-01; Geumsan-2016-02; Nonsan-2016-01; Changwon-2016-02) we limited our sampling to 47 individuals from 10 sites (Fig. 1), with all samples originating from road kills. For samples to be assigned to a same localities, they had to be collected within 10 km of each other, representative of an average maximum dispersion distance for amphibians (see review by Smith & Green, 2005). To define the GPS coordinates of the localities for pooled samples (Table 1), we calculated the centre of gravity created by the polygon of sites considered, weighted by the number of samples.

Molecular work

DNA was extracted through a Qiagen DNeasy blood and tissue kit (Qiagen, Hilden, Germany) following the instructions of the manufacturer. The targeted genes were the mitochondrial ND2 (this fragment comprised 74 bp ND1, complete tRNA-Ile, tRNA-Gln, tRNA-Met, 510 bp ND2; referred throughout as ND2 fragment for simplicity) and CR fragments. The primers used were: ND2: L-int 5′-CGA GCA TCC TAC CCA CGA TTT CG-3′ (Fu et al., 2005) and H4980 5′-ATT TTT CGT AGT TGG GTT TGR TT-3′ (Macey et al., 1998) and CR BGarF 5′-TTGGACGATAGCAAGGAACACTC-3′ and CR BGarR 5′-CCTGACTTCTCTGAGGCCGCTTT-3′ (this study). Templates were sequenced on both strands and the complementary reads were used to resolve rare, ambiguous base-calls in Sequencher v.4.9. In order to assess the positioning of the Korean individuals within the B. gargarizans species complex, several alignments and analyses were prepared depending on localities and GPS coordinates in Genbank records. A larger alignment for the RAxML and MrBayes analyses consisted of all the B. gargarizans specimens available in Genbank and other closely related taxa that fell within the B. gargarizans species complex, not associated to additional data such as exact geographical localities (n = 181, length = 1,222 base pair = bp, Fig. 2). A slightly smaller alignment that included individuals assigned to B. gargarizans senso strito that either had GPS or accurate Chinese localities (n = 148, CR = 427 bp, ND2 = 795, bp = 1,222, Figs. 3 and 4) was used to perform statistical coalescent and phylogeographic analyses (see below). Lastly, we built a smaller alignment that only included all the Korean B. gargarizans (n = 47) form Korea and all closely related matches from nearby China (n = 15; Fig. 6; Figs. S1 and S2). The length of this alignment was 422 and 795 bp for the CR and ND2 fragments, respectively. The combined data set resulted in 1,217 bp in length.

Sequences were aligned in Seaview v.4.2.11 (Gouy, Guindon & Gascuel, 2010) under ClustalW2 (Larkin et al., 2007) with default settings. Genetic p-distances (pairwise deletion) and standard error (% ± SE) were calculated using Mega 7 (Kumar, Stecher & Tamura, 2016). The most appropriate substitution model for the Bayesian inference (BI) analysis was determined by the Bayesian information criterion (BIC) in PartitionFinder v.2 (Lanfear et al., 2012). MrBayes v.3.2.6 (Ronquist & Huelsenbeck, 2003) was used with default priors and Markov chain settings, and with a random starting tree. Each run consisted of four chains of 20,000,000 generations (small alignment, n = 62, Figs. S1 and S2), 100,000,000 (large alignment, n = 181, Fig. 2) sampled each 2,000 generations (small alignment) and 10,000 (large alignment), respectively. For the BI analyses, a plateau was reached after few generations with 25% of the trees resulting from the analyses discarded as burn in. Phylogenetic relationships among haplotypes for each locus were estimated using a maximum likelihood (ML) approach, as implemented in the software RAxML v7.0.4 (Silvestro & Michalak, 2012; Stamatakis, 2014), using the default settings. Node support was inferred by the bipartition method in RAxML with 10 random addition replicates. All analyses were performed through the CIPRES platform (Miller, Pfeiffer & Schwartz, 2010). Candidate outgroup species for the trees were B. japonicus, Bufo stejnegeri and Bufo tibetanus. However, divergence of B. japonicus and B. tibetanus were low and closely grouped within the B. gargarizans clades. Thus, we employed B. stejnegeri as the most suitable outgroup for the RaxML and MrBayes analyses of the larger data set (n = 181) and used B. tibetanus in the smaller data set (n = 63). Networks were built on Haploviewer (Salzburger, Ewing & Von Haeseler, 2011) under the best tree topology as inferred in RAxML for sequenced Korean individuals (Fig. 5). The molecular clock test was done in MEGA 7 (Kumar, Stecher & Tamura, 2016), using the entire small dataset (CR+ND2 concatenated 63 samples) and applying the GTR+G model (the same used in RAxML analyses). The Partition Homogeneity Test was performed in PAUP 4 (Swofford 2003; Tamura et al., 2013), applying the partition of each gene, CR+ND2.

In order to time-calibrate the population tree, we fixed the mutation rate in the CR to 7.3 × 10⁻⁹ substitutions/site/year following (Stöck et al., 2006) and the ND2 rate estimated from the prior. This CR rate lies as a mean from previous estimates from European Bufo genus species and B. gargarizans, between 0.067 and 0.087 per lineage per million years (Stöck et al., 2006) and are established based on the last connection between North Africa and Sicily about 5.3 million years ago (M. y. a.). For population size, we used a prior normal distribution with mean = 7.3 × 10⁻⁹ and standard deviation = 0.5. In addition to the population tree, we co-estimated the dispersal history using a discrete phylogeographic (ancestral state reconstruction) model (Lemey et al., 2009) implemented in BEAST v1.8.2 (Drummond et al. 2012; http://beast.bio.ed.ac.uk). Given that our geographic sampling of populations is uneven and our state-space is low, we chose a symmetric continuous-time rate matrix. As priors for the rates, we selected the approximate reference (CTMC) prior (Ferreira & Suchard, 2008). The discrete BI analyses (alignment n = 62) was run for 10 million generations and sampling every 1,000 generations. In addition, the statistical dispersal-vicariance analysis (S-DIVA) was run using RASP v3.0 (Yu et al., 2015), allowing all possible area overlaps and keeping a maximum of three areas per reconstructed node. Results were saved in text and PNG format and then edited by plotting the maximum-clade-credibility (MCC; see below) tree in R using APE (Paradis, Claude & Strimmer, 2004; Fig. 6).

Secondly, we estimated the diffusion of the B. gargarizans CR and ND2 sequence data through time using the continuous Bayesian phylogeographic approach (Lemey et al., 2010). This approach estimates population range changes through time and ancestral population locations or origins. We used the Gamma RRW model (alignment n = 62) and Cauchy RRW model (alignment n = 148) with all individuals assigned to GPS coordinates. A random ‘jitter’ was added to each GPS coordinate with a window size of 0.5. We applied a fixed clock rate to the CR partition (see above) and estimated from the prior for the ND2. We used marginal likelihood estimations (MLE) and Bayes factors (BF) to select for the continuous trait model. MLE were calculated through path sampling (PS) and stepping stone (SS) analyses in BEAST (Baele et al., 2012). We tested for Brownian, Gauchy, Gamma and Lognormal RRW models under the strict and relaxed lognormal models running 300 million generations, with sampling every 30,000 generations. The larger data set (n = 148) revealed chain convergence problems (e.g. Gamma RRW) with MLE analyses and therefore models were tested only on the strict clock. The posterior distributions of parameter estimates were visually inspected in Tracer v1.6 (Rambaut & Drummond, 2007). For all the models tested, MLE analyses were run for 50 path steps and 100,000 generations with each step. The BF were calculated as two times the difference in MLE between different models, and the significance was determined if the BF value was >10 (Kass & Raftery, 1995). Spatial continuous space MCC trees (n = 62 and n = 148) were computed in SPREAD (Bielejec et al., 2011) and viewed in Google Earth (http://earth.google.com; Figs. 7 and 8).

Thirdly, we inferred changes in effective population size of the Korean individuals through time using a Bayesian skyline plot (BSP) model (Drummond et al., 2005) with strict clock for the CR data and same prior and a Lognormal relaxed clock for the ND2 fragment (Fig. 9). For both analyses, we ran two independent MCMC chains, each with 20 million states and sampling every 2,000 generations.

Independent runs were evaluated for convergence and mixing by observing and comparing traces of each statistic and parameter in Tracer v1.6 (http://beast.bio.ed.ac.uk/tracer). We considered effective sampling size (ESS) values >200 to be good indicators of parameter mixing. The first 10 states of the discrete, BSP and continuous Bayesian phylogenies were discarded as burnin. Samples were merged using LogCombiner v1.8.2. The resulting trees were summarised using TreeAnnotator v1.8.2, where a MCC tree with mean values was generated under heights = ca (Heled & Bouckaert, 2013). In the case of BSP, log and tree files were uploaded into Tracer, in which the plot was generated.

In addition, to examine whether the species showed any sign of historical population expansion, we estimated Tajima’s D (Tajima, 1989), Fu’s Fs (Fu, 1997) and a mismatch distribution analysis (Rogers & Harpending, 1992). Negative values of Tajima’s D can be interpreted as evidence of population expansions (Fu, 1997), and negative values of Fs indicate an excess of recent mutations and reject population stasis. In addition, we performed a Ramos-Onsins and Rozas (R²) analyses (Ramos-Onsins & Rozas, 2000). Diagrams of frequencies of pairwise genetic differences were drawn using DnaSP v.5.0 (Rozas et al., 2003). A thousand bootstrap replicates were used to generate an expected distribution using a model of sudden demographic expansion (Excoffier, Laval & Schneider, 2005). This mismatch distribution is usually multimodal in samples drawn from populations at demographic equilibrium, but is usually unimodal in populations following recent population demographic and range expansion (Slatkin & Hudson, 1991; Rogers & Harpending, 1992; Ray, Currat & Excoffier, 2003). Based on the lack of population differentiation between Korean (n = 47) and Chinese (n = 15) sequences were all considered as one population for demographic population analyses.

Discussion

Although numerous amphibian species see their range restricted to the Korean Peninsula, our results emphasise a broad distribution of B. gargarizans, supported by the presence of common haplotypes in Korea as well as North East Asia. The Korean Peninsula served as a refugium for several species during either glacial or inter-glacial periods for anurans (Zhang et al., 2008; Yoshikawa et al., 2008), as well as for other more efficient dispersers (Hosoda et al., 2000; Kim et al., 2013; Koh et al., 2014). Our findings suggest colonisation of the Korean Peninsula by individuals originating from the Chinese mainland between 2 and 0.8 M. y. a. followed by weak isolation or repeated gene exchange, as seen through the monophyly of a Korean clade within the widespread B. gargarizans species complex. This was followed by dispersal events originating from the Korean Peninsula, northwards over land and westwards over land bridges during sea level falls.

The population expansion seen on the Chinese mainland, and on then Korean Peninsula, for the last 2 m.y. started during the paeloperiod when the drainage basins of the Han, Amur, Liao and Yellow Rivers were joined into a single unit, flowing into the ocean south of the current Yellow Sea (Ryu et al., 2008). At the period in the early Quaternary, the totality of the Chinese Mainland, Taiwanese and Japanese Islands and the Korean Peninsula were connected by the continental shelf for the last time (see Fig. 1, He et al., 2015). This increased landscape connectivity (Xie, Jian & Zhao, 1995; Diekmann et al., 2008) enabled toads to breed, disperse, and colonise new areas unhindered. Besides, this period matches with the establishment and reinforcement of the monsoon weather (Wang & LinHo, 2002; Lee, Jhun & Park, 2005; He et al., 2015), which also contributed to an increase of habitat suitability for amphibians, and the proliferation of populations.

During the next glacial maximum of the next 2 m.y. the Korean Peninsula may have not been isolated during glacial maxima, but during inter-glacials, as a complementary hypothesis to the four options raised by Kim et al. (2013), and explaining the unexpected ‘expansion period (0.02 ± 0.02 [m.y.a]) corresponding to the LGM’ (Kim et al., 2013). The glacial maxima resulted in a larger continuous refugium connecting the Korean Peninsula and the Chinese mainland south of the ice line (Zhan & Fu, 2011; Yan et al., 2013), due to the absence of ice in the Korean Peninsula during the LGM (Kong, 2000; Walker et al., 2009; Yi & Kim, 2010), and the drop in sea levels, accumulated in glaciers, during glacial maxima. This will have allowed for the continuous gene exchange, and the later population expansions, followed by the northward dispersion events. The potential increased connectivity between Chinese mainland and Korean Peninsula during this period are the Mindel-Riss inter-glacial (MIS 11; in Europe: Kukla, 2005; in Russia: Galuskin et al., 2009), the Riss glacial period (MIS 6; Kazantseva glaciation, Vasil’ev, 2001; Müller, Pross & Bibus, 2003) and the glaciation period called the Wurm Glaciation in Europe (=MIS 2–4; Ehlers & Gibbard, 2008; Zyryanka in Asia: Astakhov, 1998).

The recovery of two distinct clades within the peninsula reflects likely episodes of dispersal events to/from the Korean Peninsula. The median time since the most recent common ancestor (TMRC) between the two clades generally matches with the interglacial period before the Gunz Glaciation (=Marine isotope stage (MIS) 16; Schweitzer, 2004), circa 700 K. y. a., when the Korean and the mainland became isolated due to sea level resurgence. Topographical and environmental heterogeneity is more often than not a predictor to patterns of increased genetic variability sustained by different levels of allopatric diversification. Within Korea, the lack of North–South genetic structure may have been foreseen due to the continuity of low plains facilitating dispersion of anurans, as visible through the current range of Dryophytes suweonensis (Borzée, Yu & Jang, 2016; Borzée et al., 2017) and P. chosenicus. In contrast, several studies have shown the role of East-West topographical barriers in the Korean Peninsula, contributing to allopatric speciation in amphibians, such as for Hynobius spp. (Baek et al., 2011; Min et al., 2016) and D. japonicus (Dufresnes et al., 2016). However, despite the presence of landscape barriers throughout the sampled area, such as the Baekdu Mountain Range and the sea channel between the mainland and Gangwha Island, genetic exchange was not constrained, as shown by the presence of common haplotypes on both sides of these candidate geographical barriers. This pattern was unexpected due to the topographical relief of Baekdu Mountain Range above the occurrence threshold of B. gargarizans. We interpret this apparent lack of genetic structure as a result of efficient dispersal capabilities in B. gargarizans, as seen in most toads (Vences et al., 2003), or dispersion over emerged lands during glacial maxima.

Bufo gargarizans in Chinese localities (HLJ, HB and GS in Table 1) recovered common haplotypes for the CR in Korean individuals, implying recent dispersal and gene flow between localities. The clustering of South Korean and Chinese haplotypes in a single clade subsequently to the last major interglacial period is historically coherent due to the regular fluctuations of the Yellow Sea level following minor ice ages (Jingtai & Pinxian, 1980; Oba et al., 1991; Liu et al., 2009). The Korean Peninsula was last connected to the Chinese mainland over the Yellow Sea during the late Pleistocene (Millien-Parra & Jaeger, 1999), and species were able to take advantage of land bridges to disperse until the end of the LGM, circa 15 K. y. a. (Chung, 2007; Lomolino et al., 2010). However, we need to address that while the presence of a Korean and a Korean + mainland clade may suggest possible expansions from different areas of potentially ancestral allopatric populations that reconnected in the Peninsula, this pattern can also be caused by incomplete lineage sorting, where some individuals still retain ancestral haplotypes (i.e. partial divergence with no secondary contact). Although the idea of two disjoint populations coming into secondary contact in the peninsula sounds appealing, our data does not lend to test hypotheses of admixture, thus testing apart secondary contact from lineage sorting would be difficult (Qu et al., 2012).

The case of Gangwa Island is peculiar as the two sampled individuals were breeding males, identifiable through the unmistakable forearm muscular development, but about a third smaller than average male size, 2.83 cm (± 0.12 cm), versus 8.08 cm (± 0.67 cm) in average for males B. gargarizans (Cheong, Sung & Park, 2008). This size variation may be explained by the island effect, where isolated individuals drift in average size, while being genetically closely related, i.e. same haplotype, to individuals on the mainland. This is, for instance, the case of D. japonicus between Jeju Island and the Korean Mainland (Jang et al., 2011). Further sampling in this and other satellite islands is likely to elucidate further information on the processes of dwarfism in the area.

The topographical effects caused by sea level changes throughout the ice ages in Asia and Europe are contrasting. In Europe, ice shelves lead to the southwards displacement of species resulting in isolation through allopatric speciation on southern peninsulas (i.e. Iberian, Italian, Greek, Balkans; Hewitt, 2000). In contrast, sea level recession in North East Asia at glacial maxima resulted in the drainage of the Yellow Sea (Jingtai & Pinxian, 1980; Haq, Hardenbol & Vail, 1987; Fairbanks, 1989; Park, Khim & Zhao, 1994; Kim & Kennett, 1998; Ijiri et al., 2005; Liu et al., 2009) and the near closing of the Tsushima Straight (Oba et al., 1991; Millien-Parra & Jaeger, 1999; He et al., 2015). Land bridge formations between the mainland, the Korean Peninsula (Haq, Hardenbol & Vail, 1987; Millien-Parra & Jaeger, 1999; Chung, 2007; Lomolino et al., 2010) and the then continuous formation of the Japanese Islands effectively allowed genetic admixture between populations as seen in this study, and blurring of species boundaries through introgressive extinction.