Population demographic history of a temperate shrub, Rhododendron weyrichii (Ericaceae), on continental islands of Japan and South Korea

Abstract Continental islands provide opportunities for testing the effects of isolation and migration on genetic variation in plant populations. In characteristic of continental islands is that the geographic connections between these islands, which are currently distinguished by seaways, have experienced fluctuations caused by sea‐level changes due to climate oscillations during the Quaternary. Plant populations on the islands have migrated between these islands via the exposed seafloors or been isolated. Here, we examined the demographic history of a temperate shrub, Rhododendron weyrichii, which is distributed in the southwestern parts of the Japanese archipelago and on an island of South Korea, using statistical phylogeographic approaches based on the DNA sequences of two chloroplast and eight nuclear loci in samples analyzed from 18 populations on eight continental islands, and palaeodistribution modeling. Time estimates for four island populations indicate that the durations of vicariance history are different between these populations, and these events have continued since the last glacial or may have predated the last glacial. The constancy or expansion of population sizes on the Japanese islands, and in contrast a bottleneck in population size on the Korean island Jeju, suggests that these islands may have provided different conditions for sustaining populations. The result of palaeodistribution modeling indicates that the longitudinal range of the species as a whole has not changed greatly since the last glacial maximum. These results indicate that exposed seafloors during the glacial period formed both effective and ineffective migration corridors. These findings may shed light on the effects of seafloor exposure on the migration of plants distributed across continental islands.

capacities for dispersal from migrating among islands and/or the mainlands, thus leading to increased genetic differentiation. During the late Quaternary, glacial-interglacial climate cycles occurred repeatedly on a global scale (Lisiecki & Raymo, 2005;Petit et al., 1999), and these oscillations have produced drastic changes in the connectivity of isolated populations by reconfiguring the spatial arrangements of islands due to exposure of seafloors that connected continental islands.
Spatiotemporal changes in the continuity of species distributions during climatic oscillations could therefore have been a key factor affecting the evolutionary dynamics of species on continental islands by influencing isolation and migration history. In addition, it is important for understanding how the flora of the islands has been shaped.
The Japanese Archipelago, continental islands in East Asia, is a good region for investigating effects of spatiotemporal changes in the continuity of species distributions. It has been inferred that in the Japanese Archipelago, the distribution of temperate forests became shrunken and fragmented across small patches along the coastal regions during glacial periods (Gotanda & Yasuda, 2008;Harrison, Yu, Takahara, & Prentice, 2001;Tsukada, 1985). It is increasingly recognized that the current genetic structure of plant populations has been shaped by descent from populations that sheltered in refugia (Aoki, Suzuki, Hsu, & Murakami, 2004;Iwasaki, Aoki, Seo, & Murakami, 2012;Worth, Sakaguchi, Tanaka, Yamasaki, & Isagi, 2013). However, processes of isolation and migration among refugium populations on islands via exposed seafloor are not clear and have received less attention (Burridge et al., 2013;Duncan, Worth, Jordan, Jones, & Vaillancourt, 2016). In addition, although some studies discussed effects of the last glacial maximum (LGM, Clark et al., 2009) to current genetic variation, these were based on fossil records without focal species, not based on time estimations from genetic data (Iwasaki et al., 2012;Worth et al., 2013).
Understanding the complexity of population demography (population expansion, population decline, divergence, admixture, migration, etc.) on the islands presents several difficulties. Firstly, most studies to date have examined only a few DNA sequences (e.g., regions of the chloroplast DNA [cpDNA] and/or nuclear ribosomal DNA), and stochastic variance between gene genealogies in such a limited number of genomic regions has made it difficult to obtain reliable estimates of population demographies. Secondly, if population genetic analyses are not conducted in the framework of an appropriate time scale, it is not possible to make any biogeographically meaningful inferences about the drivers of population divergence. This is particularly problematic in phylogeographic analyses of continental island species, because here population divergence can be triggered by interglacial geographic isolation and/or glacial habitat isolation. In the last decade, demographic modeling of natural populations based on the coalescent approach has become commonplace (Hey & Nielsen, 2004;Nielsen & Beaumont, 2009). This approach can statistically estimate demographic parameters such as migration rate, divergence time, and changes in population size within a given time scale. Multilocus DNA sequences offer the potential for inferring species' demographic histories in great detail because they provide many genetically informative segregating sites and independent genealogies (Wang et al., 2014).
Ecological niche modeling can also provide useful information about historical range shifts by reconstructing palaeodistributions during key periods, for example, the LGM. These approaches enable us to integrate phylogeographic inferences and information about spatial range shifts, and as a result, they can lead to deeper insights into species' biogeographic histories during the late Quaternary (Phillips, Anderson, & Schapire, 2006;Svenning, Fløjgaard, Marske, Nógues-Bravo, & Normand, 2011).
Rhododendron weyrichii Maxim. (Ericaceae) is a deciduous shrub species, occurring mainly in temperate regions that experience high rainfall during summer in southwestern parts of the Japanese Archipelago including the Kii Peninsula of Honshu, Shikoku, and Kyushu islands and small islands surrounding them, and in Jeju Island in the southwest of the Korean Peninsula, disjunctly (Figure 1; Chamberlain & Rae, 1990). These regions are thought to have been important glacial refugia for temperate plant species (Gotanda & Yasuda, 2008;Harrison et al., 2001). Rhododendron weyrichii produces numerous small seeds (30-40 mg/100 seeds; Yoichi W., personal observation). Its pollinators are butterflies and bumblebees, and thus, exchange of pollen between island populations is infrequent under present-day interglacial conditions. However, the dispersal barriers formed by seaways among neighboring islands disappeared during the LGM (ca. 21 kya) due to the shallow depth of the sea (<120 m, Clark et al., 2009;Park, Kimura, & Taira, 1996;Siddall et al., 2003), and this could have opened up effective dispersal corridors for the species.
As R. weyrichii is widespread across a number of islands and shows a moderate level of genetic differentiation among populations (Yoichi & Tomaru, 2014), and it provides a suitable study system for analyzing demographic history in the islands.
The configuration of the landscape in these island regions is believed to have affected the distribution and demographic history of plant populations. In this study, we examined the demographic history of R. weyrichii populations using statistical phylogeographic approaches based on chloroplast and nuclear DNA (nDNA) sequences.
We also examined the distribution of the species during the LGM using ecological niche modeling. The specific aims of the study were (1) to estimate population genetic divergence between continental islands and to address the isolation-migration histories of these populations; (2) to evaluate differences in genetic diversity and historical changes in population size among islands; and (3) to make inferences about population survival on each island during the LGM.
PCR was performed with an initial denaturation for 4 min at 94°C followed by 35 cycles of denaturation for 60 s at 94°C, annealing for 60 s at 55 or 60°C and extension for 60 s at 72°C, and a final extension for 7 min at 72°C, using AmpliTaq Gold Master Mix (Applied Biosystems).
After PEG precipitation (Hartley & Bowen, 1996), the PCR products were sequenced directly using the standard method provided with the BigDye Terminator Cycle Sequencing kit v. 3.1 (Applied Biosystems) and separated by electrophoresis on an ABI 3100 Genetic Analyzer (Applied Biosystems).
Mononucleotide repeats in the sequences were omitted from subsequent analyses to avoid the possibility of homoplasy. The number of tri-nucleotide repeats ([ATT] n ) found in rpl36-rps8 was counted, but these repeats were excluded from alignments when determining cpDNA haplotypes. Genealogical relationships for all cpDNA haplotypes were constructed as a parsimony network using TCS v. 1.06 (Clement, Posada, & Crandall, 2000).
F I G U R E 1 Species range and geographic distribution of (a) three haplotypes, showing the number of trinucleotide repeats (ATT) at chloroplast DNA loci and (b) two and (c) four genetic clusters detected by STRUCTURE analysis based on haplotypes at nuclear DNA loci in Rhododendron weyrichii populations. (a) Small thin letters indicate population codes corresponding to those in Table 1; numbers in pie charts indicate the numbers of tri-nucleotide repeats; a parsimony network superimposed on the map shows the relationships among the haplotypes using Rhododendron amagianum as an out-group. (b) and (c) A neighbor-joining tree superimposed on each map shows the relationships among the clusters based on net nucleotide distances; F ST of each cluster indicates the extent of genetic divergence from the ancestral population

| Estimation of genetic diversity and genetic structure
The nDNA sequences acquired were also assembled using DNA Baser and aligned using Muscle in MEGA; indel polymorphisms in these sequences were omitted from subsequent analyses. The haplotype phases of the aligned sequences were determined using PHASE v. 2.1 (Stephens & Scheet, 2005;Stephens, Smith, & Donnelly, 2001) implemented in DnaSP v. 5 (Librado & Rozas, 2009). Three independent runs employing the general model for recombination rate variation were conducted with a burn-in period of 1,000 followed by 10,000 iterations with a thinning interval of 100 steps. Those haplotypes with a posterior probability of >0.90 were used in subsequent analyses.
Homology search using blastx was conducted against polypeptide sequences from the NCBI and Arabidopsis thaliana (TAIR; http://www. arabidopsis.org) databases to identify a possible function for the product of each locus.
Number of segregating sites (S) and nucleotide diversity (π; Nei, 1987) were calculated for each locus. The minimum number of recombinations (R M ) at each locus was estimated based on the four-gamete test (Hudson & Kaplan, 1985). Neutrality at each locus was tested using Tajima's D (Tajima, 1989) and Fay and Wu's H (Fay & Wu, 2000). Fay and Wu's H was tested with out-group (R. sanctum) sequences, and significance was determined using 10,000 coalescent simulations.
In addition, an HKA test (Hudson, Kreitman, & Aguadé, 1987) was conducted to test for neutral evolution across loci between species (R. weyrichii and R. sanctum). All of these analyses were performed using DnaSP, with the exception of the HKA test which was carried out using the HKA program (http://genfaculty.rutgers.edu/hey/software#HKA).
Within-population genetic diversity was evaluated for each population.
Nucleotide diversity, Watterson's θ (Watterson, 1975), and Tajima's D were calculated over all loci within each population using DnaSP.™ Population genetic structure was then estimated as follows. The average number of nucleotide differences (D xy ;Nei, 1987) over all loci between populations was calculated. Genealogical relationships among all haplotypes at each locus were estimated by constructing a parsimony network using TCS. Population genetic structure based on genotypes (two haplotypes) at each locus for each individual was estimated by model-based Bayesian clustering analysis using STRUCTURE v. 2.3.3 (Pritchard, Stephens, & Donnelly, 2000). Admixture, correlated allele frequency and LOCPRIOR models were used (Falush, Stephens, & Pritchard, 2003;Hubisz, Falush, Stephens, & Pritchard, 2009). MCMC simulations were performed with 20 independent runs for each cluster (K = 1-10), with 3.0 × 10 4 iterations after a burn-in period of 2.0 × 10 4 .
The optimal value of K was estimated by plotting the log likelihood of the data, Ln Pr(X|K) (Pritchard et al., 2000), and the ΔK statistic, which was calculated from the second-order rates of changes of Ln Pr(X|K) (Evanno, Regnaut, & Goudet, 2005), against K. A neighbor-joining among-clusters tree was constructed based on the net nucleotide distances between pairs of clusters, and F ST values between the ancestral T A B L E 1 Locality and the numbers of individuals used in chloroplast and nuclear DNA analyses and genetic diversity and neutrality statistics over eight nuclear loci, for 18 populations from eight continental islands population and each cluster, which represent the extent of genetic drift undergone by each cluster after differentiation from the ancestral population, were calculated (Pritchard et al., 2000).

| Demographic analysis
Demographic history was firstly inferred using the isolation with migration (IM) model. To estimate the IM history among four major regional population groups (Jeju, Kyushu, Shikoku, and Kii), the population demographic parameters (Figure 3) population split time (t = Tμ), which is a product of the absolute time in years (T) and the mutation rate (μ), population size (θ = 4Nμ, where N is effective population size), migration rate per mutation event (m = M/μ, where M is migration rate), and population migration rate (2NM = 0.5θm) were estimated using IMa2 (Hey & Nielsen, 2004, 2007; see Appendix S1).
The approximate Bayesian computation (ABC) approach enables us to infer demographic history by estimating the values of parameters in a complex model without calculating its likelihood (Beaumont, 2010;Csillery, Blum, Gaggiotti, & Francois, 2010). In this study, four different demographic models related to changes in population size (

| Ecological niche modeling
To infer changes in distribution between the LGM and the present, we produced distribution models for R. weyrichii based on modern distribution records and bioclimatic variables using maximum entropy methods (Phillips et al., 2006

| Genetic diversity and structure of chloroplast DNA
The sequences of two cpDNA loci were obtained from 72 individuals across 18 populations of R. weyrichii (Table 1). The lengths of aligned sequences for the two loci were 557 bp in trnG intron and 455 bp in rpl36-rps8. There were three haplotypes with three segregating sites across the two loci ( Figure 1a). Tri-nucleotide repeats (ATT) that were identified in a rpl36-rps8 intergenic spacer showed a high level of variation in Kyushu. The geographic distribution of haplotypes showed a clear division between Shikoku and Kyushu, and a distinct haplotype was recognized in a population from Kii.

| Genetic diversity and neutrality of nuclear DNA
The sequences of eight nDNA loci were obtained from 142 individuals across 18 populations (Table 1). The length of aligned sequences for each locus ranged from 339 bp for PHYE to 470 bp for EST136 (Table 2), and the total length was 3,340 bp. The eight nDNA loci showed a high degree of polymorphism, with the number of segregating sites (S) ranging from 9 to 25 and the nucleotide diversity for all sites (π) ranging from 0.0018 to 0.0060. Recombination events were recognized for seven loci, all except PHYE, and the number of minimum recombination events (R M ) ranged from 0 to 9 (

| Population genetic structure of nuclear DNA
The pairwise genetic divergence (the average number of nucleotide differences, D xy ) between populations ranged from 0.0031 to 0.0069 (Table S2). Each nuclear locus had major haplotypes, and variable patterns of haplotype distribution were detected among loci (Figure 2).
In the Bayesian clustering analysis, K = 2 was the optimal number of clusters supported by the ΔK statistics, and K = 4 was the optimal number of clusters that showed a plateau in the log-likelihood value and low variance among runs (Fig. S1). In the case of K = 2, the distribution of clusters showed a division between Kyushu and Shikoku ( Figure 1b). In the case of K = 4, clusters 1 and 2 were recognized in the populations of the Japanese Archipelago, cluster 3 was recognized mainly in Kyushu and the adjacent islands (Koshiki, Amakusa, and Fukue), and cluster 4 predominated in Jeju and was also recognized in Fukue (Figure 1c). The F ST values for clusters 2 and 4 showed higher values than those for clusters 1 and 3; in addition, clusters 2 and 4 were distant from clusters 1 and 3. Clear genetic divergence between Jeju-Kyushu and Shikoku-Kii was also supported by the phylogenetic relationships between island populations (Fig. S2).

| Demographic history of four island populations
The posterior modes and 95% highest posterior density (HPD) intervals for the split time (t) between Jeju and Kyushu (t JJ+KY ), between Shikoku and Kii (t SK+KI ), and between Jeju-Kyushu and Shikoku-Kii In the ABC analysis of each island population, the standard neutral model (Model 1) was supported for Kii and Kyushu, whereas the exponential growth model (Model 2) was supported for Shikoku and the instantaneous size reduction model (Model 3) was supported in the case of Jeju (Table 3, Figs S3, S4) with good agreement between the posterior predictive simulations for the observed and simulated data sets (Fig. S5). It was estimated that expansion or maintenance of population size in Kii, Shikoku, and Kyushu continued over the last glacial period (Figure 4)  relationship of θ 0 among islands was the same in both analyses, including the finding that the lowest value was that for Jeju.

| Differences in species distribution between the last glacial maximum and the present time
The Maxent model showed a relatively high accuracy (area under curve, AUC = 0.940 ± 0.009). The predicted distribution based on current climatic conditions was similar to the actual distribution of the species with high probability (>0.50), although some areas along the Pacific Ocean side and on the western edge of Honshu were also predicted to be suitable (Figure 5a). The contribution made by environmental variables to the model prediction was higher for precipitation (annual precipitation = 71.8%, precipitation in the driest month = 12.0%) than for temperature (annual mean temperature = 3.8%, mean temperature in the driest season = 6.0%). Under the two LGM climate simulations (CCSM and MIROC), it was estimated that suitable environments had been maintained on three of the four islands (Figure 5b, c), but the areas on the map with high probabilities (>50%) had shrunk. In addition, the area with a suitable environment on Jeju showed low probabilities and restricted distributions in both models. It is noteworthy that although the exposed seafloor physically connected the major islands during the LGM, the ranges suitable for the species were somewhat scattered between neighboring islands.

| Isolation and migration among island populations during the late Quaternary
The estimates of the times of divergence between R. weyrichii populations on these four islands indicate that timing and duration of population divergence is different. In particular, the clear genetic divergence  between Jeju-Kyushu and Shikoku-Kii was supported by the genetic structure given by both nuclear and chloroplast DNA sequences, and the time estimate suggests that the vicariance may have continued since the last glacial or have predated the last glacial (Petit et al., 1999). While, the most recent divergence between Shikoku and Kii may indicate a role of exposed seafloor for migration corridor. However, because time estimates are seriously influenced by mutation rate, and by variations in mutation rate among sites depending, for example, on whether substitutions are synonymous or nonsynonymous (Wolf, Li, & Sharp, 1987), careful consideration is necessary when interpreting population history.
In addition, the results of a likelihood ratio test of the IM analysis suggest that a few migration events among island populations happened after their divergence. The distribution of genetic clusters detected by STRUCTURE analysis indicates that Fukue, which is located midway between Jeju and Kyushu, may have been a region of land bridges that were used like stepping stones. These evidences may mean that the exposed seafloors during the last glacial have acted as both effective and ineffective migration corridors between islands (Siddall et al., 2003).

| Differences in population survival among islands
The three island populations in the Japanese Archipelago (those on Kyushu, Shikoku, and Kii) showed large population sizes (θ 0 ) than those on Jeju according to ABC analysis; the pattern was similar for the result of IM analysis. Palynological studies have suggested that warm-temperate forests may have retreated to southern refugia and been replaced by cold-tolerant coniferous and broadleaf vegetation throughout most of the Japanese Archipelago during the LGM (Gotanda & Yasuda, 2008). In addition, previous genetic studies on many warm-temperate species have indicated that southwestern parts of the Japanese Archipelago provided important refugia (Aoki et al., 2004;Zhai, Comes, Nakamura, Yan, & Qiu, 2012). These results suggest that climate oscillations had negligible impacts at the southern tips of the merged "Japanese island" at that time, and the populations may therefore have been maintained in situ in each region throughout at least one cycle of climate change (Clark et al., 2009). Another minor factor that affects genetic diversity and population size of certain species is introgression from related species. A chloroplast haplotype unique among R. weyrichii was recognized in the north of Kii. This haplotype was shared with R. kiyosumense Makino (Yoichi unpublished data), which is distributed on the Pacific Ocean side of Honshu. An intraspecific taxon with a different flower color, which is recognized as R. weyrichii f. purpureum Hara, in the northern part of Kii is considered to have resulted from introgression events from the other species, as has also occurred in other closely related Rhododendron species (Morimoto et al., 2005;Tagane, Hiramatsu, & Okubo, 2008).
In contrast to the populations in the Japanese Archipelago, the estimated value of effective population size on Jeju was low, and this is likely to have been caused not by a recent bottleneck event. In addition, the effective population size on Jeju has not expanded after the bottleneck event. Although the time estimates given by IM and ABC analyses exhibited different values, this difference may be the result of an estimation error due to difference in examined models and model simplification. At least the accepted model and the time estimate suggested that the Jeju population is not likely to have been formed by recent migration from the Japanese island. The Jeju population is currently restricted along rivers in warm temperate forests characterized by trees of Carpinus or Quercus whereas fossil pollen records suggest that the island has been dominated by grassland with patches of cool-temperate deciduous broadleaved woodland and that conditions were dry during the LGM (Chung, 2007;Lee, Lee, Yoon, & Yoo, 2008). In addition, ecological niche modeling predicted that relatively small areas of the Jeju region were suitable, and the probability values for the distribution were not high. These lines of evidences indicated that Jeju was not suitable for R. weyrichii; however, Bayesian clustering and demographic analyses suggested that R. weyrichii may have survived in restricted and scattered habitats on or around the island during the LGM. The possibility of the invisible refugia is supported by the high genetic diversity of tree species on the island (Chung et al., 2013;Lee, Lee, & Choi, 2013;Lee, Lee, Choi, & Choi, 2014). These regional differences in population size between the islands suggest that environmental sustainability is a key factor in maintaining a population under conditions in which there are few migrations among populations (Frankham, 2005;Morjan & Rieseberg, 2004).

| Survivals of island populations from the last glacial maximum up to the present
When considering isolation and migration among islands, the LGM is an important period in history (Hewitt, 2000(Hewitt, , 2004. Under both the CCSM and the MIROC simulation, suitable areas were predicted on each of the four islands, but these areas appeared to lack spatial connections across the exposed seafloor. The existence of at least possible two regions on the islands containing refugia during the LGM, Kii-Shikoku and Kyushu-Jeju, was also supported by IM analysis. Thus, the latitudinal and/or longitudinal range of the species may not have changed greatly since the LGM as is also the case for other forest species on continental islands (Duncan et al., 2016;Qi et al., 2014;Qiu et al., 2009;Worth et al., 2013), and the current distribution may not have been shaped by northward range expansion such as has occurred for widespread species (Fujii et al., 2002;Magri et al., 2006;Sakaguchi, Takeuchi, Yamasaki, Sakurai, & Isagi, 2011).

| CONCLUSION
Statistical phylogeographic analyses of R. weyrichii populations using model-based approaches revealed that history of population isolation between islands is likely to have a considerable influence on population survival on each island, even if the islands were sometimes connected geographically during climatic oscillations (Qi et al., 2014).
Especially, differences in survival history due to the suitability of isolated habitats are reflected in the current effective population size on each island. The genetic structure of long-lived woody species would certainly have been influenced by environments during the last glacial.