Phylogeography of Pulsatilla cernua (Ranunculaceae), a grassland species, in Japan

Abstract The genetic diversity and structure of Pulsatilla cernua, a continental‐grassland relict, were investigated using variations in chloroplast DNA (cpDNA) and microsatellites of nuclear DNA. In the analyses of three cpDNA regions, 17 haplotypes were found in 24 populations of P. cernua from Japan, Korea, and Russia. Although the route and time of migration between the continent of Asia and Japan could not be well resolved, the cpDNA haplotype network suggests the existence of several ancient lineages in Japan and a recent secondary migration from Japan to the continent. Microsatellite analyses did not indicate genetic structure among the Japanese populations, indicating the existence of gene flow across the distribution area until recently. These results indicate that the present fragmentation of P. cernua in Japan may reflect a rapid, recent reduction from a previously large, continuous distribution.


| INTRODUC TI ON
Recent phylogeographic studies have explained historical changes in species distribution based on the genetic structure of extant species (Avise, 2000(Avise, , 2004. Global climatic oscillations during the Quaternary caused major changes in the distribution of many species (Hewitt, 2000(Hewitt, , 2004Listl, Poschlod, & Reisch, 2017;. Even though East Asia was primarily free of ice sheets during the last glacial period (from approximately 2.6 million to 11,700 years ago), climatic oscillations during the Quaternary influenced the distribution of vegetation (Harrison, Yu, Takahara, & Prentice, 2001;Qiu, Fu, & Comes, 2011).
The Japanese Archipelago extends along a 3,000 km line from the northeast to the southwest, covering a wide range of climatic zones, from subarctic to subtropical. The distribution ranges of most plants in the Japanese Archipelago have repeatedly shifted following climatic changes during the Quaternary (Ikeda, Carlsen, Fujii, Brochmann, & Setoguchi, 2012;Kubota, Kusumoto, Shiono, & Tanaka, 2017;Yoshida, Kudo, Shimada, Hashizume, & Ono, 2016). Because their range shifts in the archipelago interacted with the neighboring continent of Asia as well as associated islands (Chiang et al., 2014;Fujii, Ueda, Watano, & Shimizu, 1997;Ikeda, Higashi, Yakubov, Barkalov, & Setoguchi, 2014;Lee, Lee, & Choi, 2013;Nakamura et al., 2014;Sakaguchi et al., 2012), elucidating the historical interaction with the continent and surrounding islands is necessary to understand the origin of Japanese flora. Thus, investigating the processes that have led to changes in their distribution may help us to understand their origin.
In the Japanese flora, there is a group of temperate grassland plants that are distributed in northeastern China, Far East Russia, the Korean Peninsula, and Japan. Koizumi (1931 ) named them "Mansen plants," after the geographical names of the continental regions mentioned above (Man-shu and Cho-sen in Japanese). On the continent, they commonly occur in meadows and constitute temperate grassland vegetation that is widely spread across northeastern China. Their continental range is likely the place of origin of the Mansen plants (Murata, 1988;Tabata, 1997). In Japan, most of these plants occur in the temperate southwestern part of the archipelago and are not found on the northernmost large island of Hokkaido (Hotta, 1974;Koizumi, 1931;Murata, 1988). Based on their present distribution in Japan, it is hypothesized that these plants migrated to Japan via the Korean Peninsula (Hotta, 1974;Kitamura, 1957;Murata, 1988;Tabata, 1997). Their distribution areas would have been expanded under the cool and dry climate of the glacial age in Pleistocene, followed by reduction in the postglacial period. Ushimaru, Uchida, and Suka (2018) called these plants "continental-grassland relicts." Kitamura (1957) considered that many of the Mansen plants immigrated to Japan approximately 150,000 years ago through the land bridge between the Korean Peninsula and Kyushu; however, he noted the possibility of much older immigration through a northern route of ancient land bridges in several species. Hotta (1974) also postulated multiple immigrations at different periods, based on the existence of various plant groups, such as warm temperate species (e.g., Potentilla discolor Bunge) and cool temperate species (e.g., Ribes maximowiczianum Kom.).
Today, natural temperate grasslands in Japan are restricted to places influenced by periodic disturbances, such as fires, floods, and volcanic eruptions. The outlines of their distribution at present occupy a rather large area in Japan, but most Mansen plants are found only in a few small isolated populations. Warming temperatures after the last glacial period probably caused habitat fragmentation due to encroaching forests (Suka, 2012;Tabata, 1997).
Pulsatilla cernua (Thunb.) Bercht. et C. Presl., a Mansen plant (Hotta, 1974;Murata, 1977Murata, , 1988, is a perennial herb that grows in the sunny grasslands of low mountains and river floodplains of Honshu, Shikoku, and Kyushu. This species has a relatively wide distribution in Japan, from Kyushu to the northern part of Honshu. They could have immigrated to Japan from the west and expanded eastward to the northern end of Honshu. In that case, we expect closer relationship F I G U R E 1 (a) Distribution of the cpDNA haplotypes of Pulsatilla cernua, (b) cpDNA haplotype network based on statistical parsimony. Circle size is proportional to the number of populations in the haplotype. Haplotypes with one-step distance from haplotype A are colored with white between the western populations and continental ones. Otherwise, they might have moved into Japan from the north and expanded westward to Kyushu, followed by extinction in Hokkaido. To understand migration history of Mansen plants in Japan, the genetic diversity and structure among populations of P. cernua were investigated.

| Plant materials and DNA extraction
One hundred eighty-nine individuals from 24 populations of P. cernua were collected in Japan, Korea, and Russia ( Figure 1a, Table 1).
The voucher specimens were deposited at Kumamoto University (KUMA). The leaves were dried in silica gel, and total genomic DNA was extracted using a modified cetyltrimethylammonium bromide (CTAB) method (Doyle & Doyle, 1987).

| Chloroplast DNA sequencing and microsatellite genotyping
In a preliminary investigation, six noncoding regions of cpDNA were sequenced for eight populations using one plant from each population.
The simple indel coding method (Simmons & Ochoterena, 2000) was employed for gap coding.

| Haplotype network
The three noncoding cpDNA regions were subsequently concatenated and analyzed. The concatenated sequences were used to construct an unrooted haplotype network, including the 24 samples (one from each population) of P. cernua using statistical parsimony software TCS 1.21 (Clement, Posada, & Crandall, 2000). Insertions and deletions were all nonoverlapping and were included as singlegap characters for statistical parsimony analysis.

| Genetic diversity
One hundred eighty-nine samples (eight from each population, except for five from the Yamanashi population) were used for the  (Weir & Cockerham, 1984), and the fixation index (F IS ) (Weir & Cockerham, 1984) were calculated using GenAlEx 6.5 (Peakall & Smouse, 2012) and averaged over all loci in each population. Linkage disequilibria were tested by the Markov chain algorithm for all pairwise combinations of the loci with sequential Bonferroni corrections (Rice, 1989) with GENEPOP version 4.0 (Rousset, 2008). Any significant deviation of F IS from zero was evaluated with 1,000 randomizations using FSTAT version 2.9.3.2 (Goudet, 1995).

| Population structure based on microsatellite data
To examine population structure, an individual-based assignment approach was used, assuming correlated allele frequencies and admixed ancestry, and implemented in STRUCTURE 2.3 (Pritchard, Stephens, & Donnelly, 2000). The analysis assumes K clusters and assigns individuals to one or more clusters through Markov chain Monte Carlo (MCMC) simulation. STRUCTURE was run 10 times for K = 1 to K = 24 clusters, and the length of the burn-in was set to 100,000, followed by 300,000 MCMC iterations. The most likely K value was determined by calculating ∆K (Evanno, Regnaut, & Goudet, 2005). The ∆K calculations were performed using the online version of STRUCTURE HARVESTER version 0.6.91 (Earl & von Holdt, 2012). To compare the likelihood between K = 1 and K = 2, genetic differentiation of haplotypes between the two groups, designated as K = 2, was estimated by the genetic differentiation among subpopulations (F st ) using DnaSP v6 (Rozas et al., 2017).
To consider genetic relationships among eastern Japan (pops.

| Correlation between genetic distance and geographic distance
We calculated the correlation coefficient between genetic distance and geographic distance for all the populations based on microsatellite data. The correlation coefficient was also calculated for the populations excluding the Russian populations, to avoid the influence of geographic barrier that is formed by the Sea of Japan.

| Correlation between genetic diversity and geographic latitude/longitude
We calculated the correlation coefficient between genetic diversity and geographic latitude/longitude for Japanese populations based on microsatellite data.

| Chloroplast DNA sequencing analysis
The concatenated and aligned cpDNA sequence was 1,620 bp.
Seventeen haplotypes were identified among the 24 individuals sampled from the 24 populations of P. cernua. Four haplotypes (A, B, K, and P) were found in multiple populations (Figure 1a). Figure 1b shows an unrooted haplotype network using statistical parsimony.
Two haplotypes from Russia (M and Q) were genetically distant from those from Japan, and these two haplotypes, and another haplotype from Russia (I), did not form a cluster in the network.

| Genetic diversity within and between populations
None of the linkage disequilibria among the six microsatellite loci were significant after the sequential Bonferroni correction in all populations examined (p > 0.05, data not shown). Thus, the six loci were sufficiently independent to apply Bayesian clustering using the admixture, instead of the linkage model (Falush, Stephens, & Pritchard, 2003

| Population structure
STRUCTURE analyses showed that the most likely number of clusters using Bayesian cluster analysis was K = 2. Based on K = 2, we recognized two groups: red and green (Figure 2). To elucidate the main genetic structure, individuals that were assigned to a single cluster with more than 80% similarity were defined as red or green groups. However, even when analyzing these representative individuals, genetic differentiation between the red and the green groups, based on the haplotype sequences, was not significant (F st = 0.039,  Figure 4 shows the relationship between the genetic distance and geographic distance based on the microsatellite data. The genetic distance (F st ) between populations is ranged from 0.040 to 0.442

| Correlation between genetic distance and geographic distance
with an average 0.158 (Table 5). There was a positive correlation between geographical and genetic distance across populations, including or excluding Russian populations (Figure 4a,b). The correlation coefficient for all the populations was R 2 = 0.126 (p = 0.0001) and for the populations excluding the Russian populations was R 2 = 0.095 (p = 0.0001). Figure 5 shows the relationship between the genetic diversity and geographic latitude/longitude for Japanese populations based on the microsatellite data. There found weak correlation in each combination: between latitude and N A (R 2 = 0.1489, p = 0.09), latitude and T A (R 2 = 0.0947, p = 0.09), latitude and H E (R 2 = 0.0735, p = 0.25), longitude and N A (R 2 = 0.0936), longitude and T A (R 2 = 0.0828), and longitude and H E (R 2 = 0.0526, p = 0.50). However, all of the correlations were not significant.

| Genetic diversity and genetic structure
The mean N A was 1.52 (range: 1.2-2.1), and the mean H E was 0.256 (range: 0.092-0.458) for the six microsatellite loci from the 24 populations (Table 4). In comparison, another congeneric species, P. patens (L.) Mill., a widely distributed circumboreal species, exhibited higher genetic diversity: N A = 3.74 and H E = 0.541 (Szczecińska et al., 2013). The diversity of another central For genetic structure among the populations, STRUCTURE analyses based on microsatellite indicated K = 2. However, there was no significant difference between the two groups based on the haplotype sequences (F st = 0.039, p = 0.449). The result of PCA analysis on microsatellite also did not support two groups (Figure 3).
On the other hand, the genetic distance is correlated with the geographic distance (Figure 4). These results mean low genetic differentiation between populations and suggest that the present fragmentation of P. cernua may reflect a rapid, recent reduction  from a large, continuous distribution. It is noted that P. cernua was a common species in Japan before the era of high economic growth (Murata, 1988;Suka, 2012). During the 1960s, Japanese grasslands declined rapidly. Nakahama, Uchida, Ushimaru, and Isagi (2018) determined that the recent decline in the genetic diversity and effective population size of Melitaea ambigua, a grassland butterfly species, between the 1980s and the 2010s was due to the rapid loss of seminatural grasslands. The rapid loss of seminatural grasslands has also resulted in the significant reduction of genetic diversity in P. cernua. In Russia, three haplotypes (I, M, and Q) were found in three populations of Primorsky, but they were distantly related to each other. Although these haplotypes were located at derived positions from haplotype A in the network (Figure 1b)

| Migration and history of P. cernua, a Mansen plant
In general, the Mansen plants, namely continental-grassland relicts in Japan, are considered to have originated in the temperate grasslands on the continent and migrated to Japan, eastward via the Korean Peninsula, under the cold and dry climate of the glacial age (Hotta, 1974;Koidzumi, 1931;Murata, 1988). The haplotype network shows that the widely distributed Japanese haplotype A is connected to the Korean haplotype E, by the intermediary haplotype D in Kyushu ( Figure 1b). This concurs with the abovementioned western-route hypothesis. But it is not sure because only a few populations from the continent were included in this study. Also, an immigration route from the north is not denied yet. The genetic diversity of the populations of northeastern Japan is slightly higher than that of western Japan ( Figure 5). Although the correlation is not significant, it is possible that the higher northeastern genetic diversity reflects their older age.
In Japan, the large distribution area of haplotype A and the existence of its satellite haplotypes (haplotypes B, C, D, F, G, H, I, K, N, and O) imply that the rapid expansion of haplotype A happened rather recently, followed by diversification of the satellites. It is congruent with the consideration that grassland expanded on the Japanese archipelago most largely in the last glacier period (Suka, 2012;Tabata, 1997). In addition to the satellite haplotypes of the haplotype A, it is noted that the haplotypes J and L, found near the Sea of Japan, are genetically distant from other haplotypes. It is possible that there was an unknown common ancestor of the haplotypes A, J, and L before the expansion of haplotype A. On the other hand, the occurrence of haplotype I, closely related to haplotype A, in Russia suggests a recent migration from Japan to the continent. Although the genetic distance between the haplotypes of Japan and continent, microsatellite analysis did not indicate genetic differentiation among them ( Figure 3). These results suggest a relatively long and complicated migration history for P. cernua in Japan, so that multiple distribution range shifts before the last glacial period must be taken into consideration. To elucidate the migration route and the range fluctuation of P. cernua, a Mansen plant, it is necessary to do further investigation using more genetic markers and samples from the continent.

ACK N OWLED G M ENTS
Sample collection would not have been possible without the help of many people, to whom we express our special thanks ( Sawa, S. Sei, and K. Tanaka. We are grateful for the public institutions that granted us permission to collect samples. We also extend our appreciation for the members of our laboratory in supporting our research. This study was supported by a grant from the Institute of Plant Science and Resources, Okayama University for Soejima (2754 & 2835). We would like to thank Editage (www.edita ge.jp) for English language editing. Our cordial thanks are going to the referee who kindly gave us advices to revise our manuscript.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
A.S. conceived the ideas; A.T., A.E.K, Z.V.K., N.F., and A.S. made fieldworks and collected the data; and A.S., A.T., and H.I. analyzed the data; and A.S. led the writing.

DATA ACCE SS I B I LIT Y
All chloroplast sequences and microsatellite genotypes are submitted to the Dryad database on 13 May 2019. https ://doi.org/10.5061/ dryad.3mf353m.