- Research article
- Open Access
A strong ‘filter’ effect of the East China Sea land bridge for East Asia’s temperate plant species: inferences from molecular phylogeography and ecological niche modelling of Platycrater arguta(Hydrangeaceae)
BMC Evolutionary Biologyvolume 14, Article number: 41 (2014)
In East Asia, an increasing number of studies on temperate forest tree species find evidence for migration and gene exchange across the East China Sea (ECS) land bridge up until the last glacial maximum (LGM). However, it is less clear when and how lineages diverged in this region, whether in full isolation or in the face of post-divergence gene flow. Here, we investigate the effects of Quaternary changes in climate and sea level on the evolutionary and demographic history of Platycrater arguta, a rare temperate understorey shrub with disjunct distributions in East China (var. sinensis) and South Japan (var. arguta). Molecular data were obtained from 14 P. arguta populations to infer current patterns of molecular structure and diversity in relation to past (Last Interglacial and Last Glacial Maximum) and present distributions based on ecological niche modelling (ENM). A coalescent-based isolation-with-migration (IM) model was used to estimate lineage divergence times and population demographic parameters.
Combining information from nuclear/chloroplast sequence data with nuclear microsatellites, our IM analyses identify the two varieties as genetically distinct units that evolved in strict allopatry since the mid-Pleistocene, c. 0.89 (0.51–1.2) Ma. Together with Bayesian Skyeline Plots, our data further suggest that both lineages experienced post-divergence demographic growth, followed by refugial isolation, divergence, and in the case of var. arguta post-glacial admixture. However, past species distribution modelling indicates that the species’ overall distribution has not greatly changed over the last glacial cycles.
Our findings highlight the important influence of ancient sea-level changes on the diversification of East Asia’s temperate flora. Implicitly, they challenge the notion of general temperate forest expansion across the ECS land bridge, demonstrating instead its ‘filter’ effect owing to an unsuitable environment for certain species and their biological (e.g., recruitment) properties.
Phytogeographers have long recognised the Sino-Japanese Floristic Region (SJFR) of East Asia as the world’s major centre of temperate plant diversity and endemism . Much of the potential primary vegetation of this vast region is composed of warm-temperate deciduous (WTD) forest, as presently found scattered in mid-elevation subtropical (Central/South/East) China, predominant in low-elevation North China and the Korean Peninsula, and disjunctively distributed in the main islands of Japan [1, 2]. Fossil pollen analyses have previously indicated that during the Last Glacial Maximum [LGM: c. 21,000–18,000 yr before present (BP)], the habitat of East Asian WTD forests in the northern parts of their range (e.g., in North China and North Japan) contracted, mainly in response to increased aridification [3, 4]. However, palaeo-biome reconstructions suggest that these forests also expanded across the large expanses of continental shelf (c. 1 million km2) that emerged in the East China Sea (ECS) as a consequence of sea level lowering by c. 85–130/140 m during cold periods [5–7]. In consequence, it is now widely believed that a near continuous belt of WTD forests spanned the ECS continental shelf during the LGM (and possibly earlier cold periods), connecting presently disjunct populations in East China, South Japan, and the Korean Peninsula [2, 8]. If correct, this land bridge may have served as a ‘dispersal corridor’ [9, 10], allowing intermittent migration of most WTD forest-restricted plant species from the Asian mainland into Japan (or vice versa), and/or periodic secondary contact and gene flow among formerly isolated populations, possibly up until the last shelf submergence (c. 16,000–7,000 yr BP; ). However, there is accumulating but still scanty evidence that the ECS land bridge also acted as a ‘filter’ in selectively hampering or even preventing the dispersal of certain plant species, while allowing those able to tolerate the environmental conditions of this palaeo-landscape to disperse more freely (see below).
Support for the ‘recent dispersal corridor hypothesis’ comes from phylogeographic data of two wide-ranging and abundant WTD forest tree species of the SJFR (Cercidiphyllum japonicum and Kalopanax septemlobus). In both instances, there is a near lack of differentiation for chloroplast and/or nuclear DNA sequences across the ECS, consistent with ecological niche modelling (ENM) predicting large expanses of suitable habitat for each species on the ECS land bridge during the LGM. By contrast, a high level of genetic differentiation across the ECS has been identified in several rare understorey herbs and shrubs restricted to the mountainous WTD forests of East China and South Japan (Croomia japonica, Platycrater arguta/Kirengeshoma palmata[15, 16] and Ligularia hodgsonii). In all of these latter cases, the estimated times of trans-ECS divergence based on chloroplast (cp) DNA fall into the early-to-mid Pleistocene, suggesting that the ECS land bridge imposed a formidable barrier to dispersal during the last glacial cycle(s) (reviewed in ). However, all of these previous time estimates relied on a single locus, i.e., chloroplast (cp) DNA, which thus renders them subject to error from the inherent stochastic nature of the coalescent [19–21].
The aim of the present paper is to re-examine the evolutionary history of one of these understorey shrubs, Platycrater arguta Siebold & Zucc. (Hydrangeaceae), using a multi-locus approach. Platycrater arguta, the only species of this monotypic genus, is a small deciduous shrub endemic to the mountainous WTD forests of East China and South Japan, where respectively two varieties, var. sinensis and var. arguta, are recognized [22–24]. In the cpDNA study previously conducted , these taxa were found to comprise distinct phylogroups, whose likely vicariant origin was dated to the mid-Pleistocene (c. 0.89 Ma). Here we use a broader sampling of P. arguta, and present additional datasets of two nuclear DNA (nDNA) sequence markers (ITS, Tpi) and nuclear microsatellite (nSSR) loci to infer a more robust divergence and demographic history of this species. Specifically, our goals were (i) to provide a refined time-scale for the divergence of var. sinensis and var. arguta by fitting all four datasets (ITS, Tpi, nSSRs, cpDNA) to an ‘isolation with migration’ (IM) model [20, 25]; (ii) to quantify the amount of post-divergence gene flow between them while accounting for potential changes in effective population sizes; and (iii) to model the species’ potential distributions in response to past climatic changes, specifically from the Last Interglacial (LIG/Eemian: c. 130,000–114,000 yr BP)  through the LGM to the present day. We attempted to reconcile these distribution patterns with our genetic data to determine the role of the exposed ECS shelf as a ‘corridor’ or ‘filter’ during lineage divergence at the time of the last glacial cycle(s). Finally, together with ENM, the added value of three nuclear data sets allowed for more confidence in the interpretation of the supposedly contrasting population histories of each variety than was derived previously from a single dataset (cpDNA) .
Plant material and sampling design
We sampled seven populations of each of var. sinensis and var. arguta, (272 individuals in total) covering the species’ entire range in East China and South Japan (Additional file 1: Table S1). Ten of the populations (100 individuals) were analysed previously for cpDNA , while four were newly sampled (C3–C5 and J2). As P. arguta is nested within a paraphyletic Hydrangea L., but with yet undefined sister species [27, 28], we arbitrarily designated H. chinensis Maxim., collected from Yandang Mt. (China), as one of the outgroups together with other members of Hydrangeaceae (see below). Voucher specimens of this species and all sampled populations of P. arguta are stored at the Herbarium of Zhejiang University (HZU).
Total genomic DNA was extracted from silica-dried leaf tissue by the cetyltrimethyl ammonium bromide (CTAB) method . The entire internal transcribed spacer (ITS) region of nuclear ribosomal (nr) DNA was sequenced in 72 individuals of P. arguta (4 to 8 individuals per population; Additional file 1: Table S1) and one individual of H. chinensis, following the methods described in .
After preliminary screening of six low-copy number nuclear genes using primers developed by Strand et al.  (ADHX2F–4R, CAMX1F–2R, CHIX1F–4R, CHSX1F–2RN, GPDX7F–9R, TPIX4FN–6RN), we chose the triose-phosphate isomerase (Tpi) gene (TPIX4FN–6RN) for the full survey of 173 individuals (5–20 individuals per population, plus H. chinensis; Additional file 1: Table S1) because this gene region was single copy in most cases and proved to be sufficiently polymorphic . Direct sequencing of Tpi was carried out using the same procedure as described in Zou et al.. Chromatograms of ITS and Tpi with additive peaks were further analysed by inferring the identity of the two haplotypes within a heterozygote through haplotype subtraction [33, 34]. When the chromatogram quality did not permit this procedure, PCR products were cloned using a pMD18-T vector system (Takara Biotechnology, Dalian, Liaoning, China) according to the manufacturer’s protocol. Six to ten clones were sequenced per individual when sequences contained two or more polymorphic sites. The two sequences of a heterozygote were separated by comparing sequences of the PCR product and the cloned sequences. All haplotype sequences have been deposited in GenBank [accession numbers for P. arguta: JQ978221–JQ978253 (ITS) and JQ978254–JQ978284 (Tpi); for H. chinensis: KF559183 (ITS) and KC853063 (Tpi)].
All 272 P. arguta individuals were genotyped at seven nuclear dinucleotide-repeat microsatellite loci (Pa1–Pa3, Pa5–Pa8) (Additional file 2: Table S2) according to the methods described in Qi et al.. Fluorescently labelled PCR products (HEX or 6-FAM; Applied Biosystems, Foster City, CA, USA) were separated on a MegaBACE 1000 (GE Healthcare Biosciences, Pittsburgh, PA). Alleles were scored manually with the aid of GENETIC PROFILER 2.2 (GE Healthcare Biosciences) using the ET550-R size standard.
Nuclear DNA sequence analyses
Sequences of ITS and Tpi were aligned and edited in GENEIOUS 4.8.2 . For each gene region, we estimated haplotype diversity (h) and nucleotide diversity (π) for each population, for each variety, and for the whole data set. Tests for non-neutral evolution were conducted by computing Tajima’s D and Fu & Li’s D* . In addition, we estimated Fu’s Fs  and Ramos-Onsins and Rozas’ R2 to detect population growth. Critical values for these statistics were obtained using 105 coalescent simulations. Recombination was calculated as the minimum number of recombination events (Rm) using the four-gamete test . All of the above analyses were performed in DNASP 4.1  and ARLEQUIN 3.11 .
Phylogenetic relationships of ITS and Tpi haplotype sequences were inferred using maximum parsimony (MP) and maximum likelihood (ML), with gaps (indels) treated as missing data. Heuristic MP searches were performed in PAUP* 4.0b10  using the same settings as described in Qiu et al. . The ML analysis was conducted under the GTR + Γ substitution model using RAXML 7.2.8 . Node support was assessed using 1,000 ‘fast bootstrap’ replicates. The following species were chosen as outgroups: Hydrangea chinensis (ITS, Tpi) as well as H. anomala D. Don and Schizophragma hydrangeoides Siebold & Zucc. (only ITS; downloaded as GenBank accessions JF976651 and GU98303, respectively). The alignments and phylogenetic trees were deposited in TreeBASE (submission number 15354). In addition, we constructed haplotype networks for each dataset using the 95% statistical parsimony criterion implemented in TCS 1.21 . For ITS, however, we had to increase the TCS connection limit to 50 steps to link the divergent networks of the two varieties. Indels were treated as single mutation events, and coded as substitutions (A or T). Population differentiation for unordered (GST) and ordered (NST) haplotypes were obtained with PERMUT to test whether NST is significantly larger than GST (1,000 random permutations), indicating the presence of phylogeographic structure . Analyses of molecular variance (AMOVAs) were carried out in ARLEQUIN using F-statistics. Sequence variation was hierarchically partitioned between the two varieties, among populations within varieties, and within populations. The significance of all estimated fixation indices was tested using 10,000 permutations as described in Excoffier et al..
Nuclear microsatellite analyses
Measures of nSSR diversity were assessed for each population, and across all loci, by calculating in FSTAT 2.9.3  the total number of alleles (NA), allelic richness (RS; standardized for 5 individuals using rarefaction), expected gene diversity (HS), and the inbreeding coefficient (FIS). Differentiation between populations was computed in ARLEQUIN using RST, which assumes a stepwise mutation model (SMM; ), and thus may be more realistic for microsatellite data than other measures (e.g., FST) based on the infinite alleles model (IAM).
Overall population structure was examined using the Bayesian clustering algorithm implemented in STRUCTURE 2.3 . This program was run 10 times on individual multi-locus nSSR genotypes for a number of clusters (K), ranging from 1 to 14 (the number of localities), using a burn-in length of 50,000 and run length of 500,000 iterations. We used the admixture model without prior information on sample population membership and allowed allele frequencies to be correlated among clusters . We plotted the posterior probability of the data [lnP(D)] and the ad hoc statistic ΔK for determining the most likely K.
Hierarchical AMOVA was performed using R-statistics, and partitioned as described above. In order to detect genetic signatures of recent population bottlenecks in the nSSR dataset, we applied two tests implemented in BOTTLENECK 1.2.02 : (i) Wilcoxon’s sign-rank test for heterozygote excess; and (ii) the ‘mode-shift test’ for detecting a shifted, rather than an equilibrium, L-shaped distribution of alleles.
Divergence time, effective population sizes, and gene flow
We used the ‘isolation with migration’ (IM) coalescent model as implemented in IMA to estimate population rate parameters (Θ) and effective population sizes (Ne) of var. sinensis in China (ΘC), var. arguta in Japan (ΘJ) and their common ancestral population (ΘA; Θ = 4Neu), as well as bidirectional migration rates (mC–J and mJ–C; M = m/u) and divergence times (τ = tu) between the two varieties. All parameters in the IM model are scaled by the neutral mutation rate (u). For this analysis we jointly employed the present nuclear datasets (ITS, Tpi, nSSRs) and previously obtained sequences of two cpDNA regions (psbA–trnH, trnD–trnE; ). Since IMA assumes no recombination within loci, the nDNA sequence data were analysed by using only maximally informative and nonrecombining blocks per individual as identified by the program IMGC.
Both nDNA and cpDNA sequence data were analysed under the HKY model rather than the infinite sites model (ISM) as we regularly found some sites with more than two changes in these data sets (note that HKY and ISM are the only models of nucleotide evolution implemented in IMA, whereby only HKY allows multiple changes at a site). Corresponding values of u were calculated as u = μkg, where μ is the number of substitutions per site per year, k is the average sequence length under study in base pairs (excluding indels), and g is the generation time in years (i.e., age of first reproduction), which is five years in P. arguta as observed under cultivation . As substitution rates for this species are unknown, we assumed the following mean values and confidence intervals (CI) in substitutions per site per year: ITS (for woody perennials), 2.15 × 10−9 (0.38–7.83 × 10−9) ; Tpi, 7 × 10−9 (0.5–3.0 × 10−8) [58, 59]; and cpDNA, 1.52 × 10−9 (1.0–3.0 × 10−9) [58, 60]. The nSSR loci were assumed to fit a stepwise model of mutation (SMM), with a mean mutation rate (μ) of 10−5 mutations per locus per year. Although the assumptions underlying this mutation rate are debatable (e.g., [61, 62]), this rate falls close to the median of average values (c. 7.7 × 10−4) reported for a wider range of plant species [63–67], and has also been recently employed for other woody perennials [68, 69]. The geometric average mutation rate of the four marker sets was used to rescale the IMA parameter estimates from the combined analysis. The inheritance scalars were set to 1 for the nuclear markers, and 0.5 for cpDNA, the latter value reflecting the expected effective population size (Ne) of maternally inherited cpDNA in a hermaphroditic plant relative to an autosomal locus  (note that maternal inheritance of cpDNA has been demonstrated in Hydrangeaceae ). Markov chain Monte Carlo (MCMC) simulations were conducted for 107 steps by using a linear heating scheme and 10 Metropolis-coupled chains with a burn-in period of 105 steps. To verify convergence on the same parameter values, we ran this analysis three times with different random seeds. Only estimates whose posterior distribution dropped to zero within the prior intervals were trusted.
To further examine the historical demography of each variety, we estimated the shape of the population growth function through time by constructing Extended Bayesian Skyline Plots (EBSPs) as implemented in BEAST 1.7.5 . This coalescent-based, nonparametric Bayesian MCMC algorithm incorporates multi-locus data to reduce estimate errors associated with single genes and increases the power to detect demographic dynamics . Analyses were performed for each variety across the sequence data sets (ITS, Tpi, cpDNA) assuming a strict molecular clock. The mutation rate priors for each locus were identical to those used for IMA. Based on the Akaike Information Criterion (AIC) as implemented in JMODELTEST, we selected the GTR + Γ and HKY substitution models for the nuclear and cpDNA sequences, respectively. Bayesian MCMC chains were run for 10 million generations, with a sampling frequency of every 1,000 generations, whereby the first 5,000 samples were discarded as burn-in. Convergence of the MCMC chains was inspected using TRACER 1.5  by visually checking that effective sample size (ESS) for all relevant parameters were well above 200. Skyline plots were visualized using EXCEL.
Present and past distribution modelling
To infer distributional changes of P. arguta since the last interglacial, we produced ENMs for three time periods (present, LGM, LIG) using the maximum entropy method implemented in MAXENT 3.2.1 . In addition to the distribution records included in this study, records were sourced from the Chinese Virtual Herbarium (http://www.cvh.org.cn), the National Specimen Information Infrastructure of China (http://www.nsii.org.cn), the Kyoto University Herbarium (KYO), and Red Data Books for the Aichi and Mie Prefectures of Japan [77, 78]. Based on a total of 59 records (China: 20; Japan: 39), a current distribution model was developed using six bioclimatic data layers (bio1, 12, 16, 17, 18, 19; WorldClim dataset, [79, 13]) at 2.5 arc-min resolution. This model was then projected onto the set of climatic variables simulated by the Model for Interdisciplinary Research on Climate (MIROC) 3.2  to infer the extent of suitable habitat during the LGM and the LIG [81–83, 13]. The accuracy of each model prediction was then tested by calculating the area under the ‘Receiver Operating Characteristic (ROC) Curve’ (AUC; ). We note that the ENM for the LGM explicitly accounted for changes in palaeo-coastlines (−110 m than at present) and palaeo-climate surfaces on the exposed ECS continental shelf [7, 85].
Nuclear sequence characteristics
The ITS sequences of the 72 P. arguta individuals (14 populations) from East China (34) and South Japan (38) were aligned with a total length of 737 base pairs (bp), revealing 36 nucleotide substitutions and three 1-bp indels. Together, these 39 polymorphic sites yielded 33 ITS haplotypes (‘ribotypes’, H1–H33) (Additional file 3: Table S3). For the Tpi locus, only one or two distinct sequences were detected in each of the 173 individuals surveyed, 28 of which were found to be heterozygotes. In total, 31 haplotypes (T1–T31), ranging from 309 bp to 314 bp, were designated based on 39 substitutions and three small (≤ 5 bp) indels (Additional file 4: Table S4). None of the loci showed significant deviation from neutral expectations for Tajima’s D or Fu and Li’s D* at the species or variety levels (Additional file 5: Table S5). Demographic tests of Fu’s Fs revealed negative but nonsignificant values, while R2 was consistently positive and significant, suggesting demographic growth (Additional file 5: Table S5). The minimum number of recombination events (Rm) estimated for ITS and Tpi were 15 and 2, respectively.
Genetic diversity, haplotype relationships and genetic structure at ITS and Tpi
There were high and comparable levels of total haplotype (hT) and nucleotide (πT) diversity at the species-level (ITS/Tpi: hT = 0.96/0.92; πT = 0.03923/0.02133), and the same was found for each variety (region), whereby each of the seven populations from East China (ITS/Tpi: hT = 0.88/0.90; πT = 0.0151/0.0108) and South Japan (ITS/Tpi: hT = 0.95/0.87; πT = 0.0188/0.0139) harboured broadly similar levels of diversity (Additional file 1: Table S1). Concordant with the previous cpDNA data , there was no sharing of haplotypes between China and Japan for either ITS (Figure 1a) or Tpi (Figure 2a).
For each of these nuclear markers, the tree topologies recovered from MP and ML analyses were similar, and only the MP strict consensus trees are shown. According to the ITS tree (Figure 3), rooted with three outgroup species, P. arguta was monophyletic (with MP/ML bootstrap values of 100/100%) and consisted of two major clades corresponding to var. sinensis from East China (100/99%) and var. arguta from South Japan (84/69%). In addition to these lineages previously identified by cpDNA, the ITS marker revealed two novel subclades, 1 (98/89%) vs. 2 (99/98%), predominant in the southern vs. northern parts of South Japan (Kyushu, Shikoku vs. central Honshu) but with apparent range overlap in the Kii Peninsula of south-central Honshu (population J5; Figure 1a). As can be seen from the unrooted ITS haplotype network (Figure 1b), the Chinese and Japanese haplotypes were separated from each other by 28 steps, and the two Japanese subclades by 13 steps.
The Tpi phylogeny, rooted with Hydrangea chinensis (Figure 4), had lower resolution than the ITS tree but still recovered all Chinese haplotypes as monophyletic (85/87%), and the same applied to the majority of those from Japan (100/100%). Both clades, however, formed a trichotomy with a deviant clade (86/93%) of the two remaining haplotypes (T5, T15) from Japan. These latter haplotypes, which showed no distinct geographic distribution (Figure 2a), occupied an intermediate position in the unrooted Tpi haplotype network (Figure 2b), yet with somewhat closer relationships to China rather than Japan (four versus nine mutational steps apart from each group).
For both ITS and Tpi, significant phylogeographic structure was obvious at the species-level and within each of the two varieties (NST > GST, all P values < 0.01; Additional file 6: Table S6). Hierarchical AMOVA revealed that most of the total nuclear sequence variation was partitioned between them (ITS: 70.3%; Tpi: 66.7%; Table 1). Nevertheless, for each variety (sinensis/arguta), NST values were high at each nuclear sequence marker (ITS: 0.751/0.907; Tpi: 0.458/0.775) (Additional file 6: Table S6), reflecting that the majority of haplotypes were population specific (ITS: 32/33; Tpi: 26/31; see also Figures 1a and 2a).
In 272 individuals from 14 populations of P. arguta, we detected a total of 220 alleles across 7 nSSR loci, with 17 to 48 alleles per locus. Average allele number (NA) was higher in China (mean: 66) than in Japan (mean: 43), but allelic richness (RS) and expected gene diversity (HS) were very similar (China: RS = 4.69, HS = 0.734; Japan: RS = 4.68, HS = 0.740; Additional file 1: Table S1). Within-population FIS values ranged from −0.054 to 0.480, with an average of 0.246.
The STRUCTURE analyses provided strongest support for K = 7, both when considering the probability of the data LnP(D) and ΔK (Additional file 7: Figure S1). With K = 2, individuals from one population in China (C7) and Japanese populations formed a joint cluster, while the two varieties remained separate from K = 3 upwards (Figure 5). At K = 7, five clusters were recovered in China (var. sinensis) and two in Japan (var. arguta), where most populations (J2–J5) comprised individuals representing both clusters (Figure 5).
The overall RST was 0.340, reflecting strong genetic differentiation over all populations and within each variety (var. sinensis/arguta: 0.273/0.366). However, in contrast to ITS and Tpi, the total genetic variance partitioned between the two varieties was relatively low (13.5%; Table 1), suggesting pronounced allele sharing. Wilcoxon’s test and the complementary mode-shift test provided little evidence for recent population bottlenecks, except for three populations from Japan (J2, J4, J6) and partly depending on the test applied (Additional file 8: Table S7).
Isolation-with-migration (IM) analyses
For this analysis, largest nonrecombining blocks of nDNA (ITS, Tpi) sequences were employed together with the nSSR markers and the previously obtained cpDNA sequences . The maximum-likelihood estimates (MLEs) and the 90% highest probability density (HPD) intervals of the six IMA-derived parameters are shown in Table 2, and their marginal posterior probability (MPP) distributions are illustrated in Additional file 9: Figure S2. Based on the geometric average mutation rate calculated (5.79 × 10−6 mutations per locus per year), these parameter estimates were converted to absolute values of years or individuals (Table 2). The split between var. sinensis and var. arguta was dated to about 889,358 yr BP, with a 90% HPD interval ranging from 509,438 to 1,193,295 yr BP (Additional file 9: Figure S2a, Table 2). For the ancestral (ΘA) and descendant (ΘC and ΘJ) population rate parameters, both var. sinensis (NC = 1.13 × 105) and var. arguta (NJ = 6.00 × 104) experienced a marked increase in effective population size (Ne) relative to that of their common ancestor (NA = 2.73 × 104) (Additional file 9: Figure S2b, Table 2). Peak posterior estimates of post-divergence migration were effectively zero in both directions (mC–J = mJ–C = 0.005; Additional file 9: Figure S2c, Table 2).
Bayesian skyline plot analysis of historical demography
The EBSPs show that both varieties maintained low but stable effective population sizes (Ne) up to approximately 0.4–0.5 Ma, and appear to have experienced a constant increase since (Additional file 10: Figure S3). While depicting a demographic trend over time, this analysis, however, cannot precisely estimate Ne because of the very broad confidence intervals (see Additional file 10: Figure S3).
Present and past ecological niche models
The AUC value for the current potential distribution of P. arguta was 0.981, indicating a good predictive model performance. For the present (Figure 6a), these predicted areas mainly include the species’ known distribution ranges in East China (Wuyi/Yandang Mts.) and South Japan (Kyushu, Shikoku, Kii Peninsula/south-central Honshu, and the Pacific Ocean side of central Honshu). Further suitable habitat is predicted in central-eastern Taiwan, but where the species is not known to occur. During the LIG (Figure 6b), the species’ potential range in China was apparently reduced compared to the present. Also in the LGM (Figure 6c) only small and disjunct areas are predicted as suitable in the Wuyi/Yandang Mts. In contrast, in Japan, the species’ current range is more or less similar to that during the LIG (Figure 6b), while during the LGM suitable habitat apparently contracted to the south (Kyushu, Shikoku) and the more northerly located Kii Peninsula (Figure 6c). Most strikingly, our LGM distribution model indicates almost no hospitable areas on the exposed ECS continental shelf, excepting those in the far east (i.e., extending from Kyushu to the delta region of the palaeo-Yellow River).
Intraspecific divergence and post-divergence gene flow
Our results provide strong evidence that the two recognized varieties of the WTD understorey shrub P. arguta in East China (var. sinensis) and South Japan (var. arguta) comprise reciprocally exclusive haplotypes for both ITS and Tpi (Figures 1 and 2). This is concordant with patterns for cpDNA  and consistent with the hypothesis that Chinese and Japanese conspecifics diverged over multiple glacial periods without inter-population gene flow. However, relationships inferred from the STRUCTURE analysis of nSSRs with K = 2 seem to conflict with this interpretation, as one Chinese population (C7) clusters with those from Japan (Figure 5). Given the congruence between nuclear and chloroplast sequence data, the sharing of nSSR alleles most likely reflects homoplasy and/or incomplete lineage sorting rather than admixture through rare long-distance dispersal and/or migration across the ECS land bridge. In fact, our IMA analysis of combined multi-locus data (cpDNA, ITS, Tpi, nSSRs) recovered only close to zero signals of bidirectional post-divergence gene flow between the two varieties following their estimated divergence in the mid-Pleistocene, c. 0.89 Ma (90% HPD: 1.2–0.51 Ma). Although this timing perfectly matches the respective point estimate drawn from cpDNA alone , single-locus cpDNA analysis suffered from a lack of convergence of the lineage divergence parameter (t), a problem that has been overcome with the present multi-locus approach (Additional file 9: Figure S2).
Broad-scale biogeographic history of P. arguta
According to Kimura [86, 87], there were three main stages of land connections between the Eurasian mainland and the Japanese/Ryukyu Islands since the Late Miocene, with most of the ECS sea-floor exposed at about 7.0–5.0 Ma, 2.0–1.3 Ma, and 0.2–0.015 Ma, whereby the latter interval includes land bridge formations in both the penultimate (Riss) and last (Würm) glacial periods [88–90]. Our estimated divergence time and its HPD intervals fall into the period (1.3) 1.0–0.2 Ma, the so-called ‘Ryukyu Coral Sea Stage’, when sea level rose tremendously and most formerly exposed land area submerged . Based on these palaeo-data, and with no genetic indication of long-distance dispersal, the disjunct distribution of P. arguta across the ECS is implied to have resulted via vicariance, that is, the species was probably more widely distributed on the ECS land bridge in the early Pleistocene, around 2.0–1.3 Ma, while divergence was subsequently driven by population contraction and extinction through land-bridge submergence.
Unfortunately, it is not possible to test this ‘early-Pleistocene expansion’ hypothesis further using ENM, as no proxy-climate records are currently available for time periods earlier than the Last Interglacial (LIG) . However, loess deposit and marine δ18O records indicate that both aridification and the influence of cold and dry winter monsoons in the northern subtropical areas of East Asia were less extensive before the mid-Pleistocene climate transition (MTP), c. 0.85 Ma . Afterwards, global ice volume drastically increased along with a change in the dominant orbital cycle from 41,000 to 100,000 years, resulting in drier and colder climates during glacial periods . Given that P. arguta is a moisture-dependent and drought-sensitive shrub, early-Pleistocene environmental conditions may have still favoured a more continuous distribution on the ECS land bridge, whereas subsequent glacials would have prevented the migration of individuals between China and Japan. Indeed, our LGM distribution model for P. arguta essentially shows no suitable habitat areas on the exposed ECS basin, with the exception of the delta region of the palaeo-Yellow River (Figure 6c). We therefore conclude that the ECS land bridge acted as a formidable barrier to the dispersal of P. arguta during the LGM and earlier cold periods of the Late Pleistocene. However, additional factors, other than climate-related niche requirements, may have had a role in preventing, or at least hampering such dispersal. This is because similar (i.e., genetic and ENM) evidence suggests that the ECS land bridge served as an intermittent route of range expansion for most of the Late Pleistocene to counteract isolation between Chinese and Japanese populations of two widespread WTD forest tree species, namely Cercidiphyllum japonicum and Kalopanax septemlobus. Both are tall canopy trees with long generation times, large effective population sizes and high seed production [94, 95], and additional traits generally considered important for population recovery, such as a generalist (i.e., wind) pollination/dispersal system (C. japonicum) and vegetative reproduction (C. japonicum, K. septemlobus). Similar traits are largely absent or only moderately developed in P. arguta (except for its wind-dispersed seed); the same applies to other rare and narrowly distributed understorey plants of the same forest biome with likewise ancient (early-to-mid Pleistocene) genetic breaks across the ECS (Croomia japonica: ; Kirengeshoma palmata: ; Ligularia hodgsonii: ). Moreover, in forest habitats particularly, one should not dismiss the effects related to edges (e.g., abundance to animal pollinators and/or dispersers), which may either increase or decrease seed production and recruitment [96–98]. Hence, species in the forest interior, such as P. arguta and other understorey plants, should be negatively affected by fragmentation; by contrast, species using edge or transitional habitats, such as the riparian-dwelling C. japonicum and the semi-invasive K. semptemlobus, may be favoured by fragmentation [96, 97]. We therefore propose that (i) range expansion in response to the formation of the glacially exposed ECS land bridge is species-specific; and (ii) predictions on the effects of this ephemeral land bridge, as ‘corridor’ or ‘filter’, have to account not only for habitat preferences per se but also for other biological features of each species, especially its recruitment properties.
Inferences of historical demography and range dynamics in P. arguta
Our results provide strong evidence that both var. sinensis and var. arguta underwent past population growth following their inferred sundering in the mid-Pleistocene. This is supported by significantly positive R2 statistics for both ITS and Tpi (Additional file 5: Table S5), and is also evident from the joint analysis of the four genetic data sets (cpDNA, ITS, Tpi, nSSRs) with IMA, suggesting a somewhat larger effective population size (Ne) in each variety compared to their ancestral population (Additional file 9: Figure S2b, Table 2). We caution, however, that this estimate of ancestral Ne probably reflects post-vicariant conditions, and thus is biased downward, as complex population dynamics within the ancestral population, such as contractions and extinctions, are not accounted for by IMA. Nonetheless, analyses of the combined nuclear and plastid sequence data with EBSPs were consistent in showing a strong increase of Ne in each variety from about 0.4–0.5 Ma onwards (Additional file 10: Figure S3). This is very similar to the point estimates of expansion times inferred from mismatch analyses of cpDNA alone, that is, c. 0.43 and 0.45 Ma for var. sinensis and var. arguta, respectively . Hence, the present results support our previous notion that this near-synchronized population growth may have been triggered by climate change at the beginning of China’s ‘Penultimate Interglacial Period’ (c. 0.46–0.33 Ma; ). This also accords with palaeo-climate data, suggesting that during interglacials since the mid-late Pleistocene (c. 1.0–0.78 Ma) the warm, wet East Asian summer monsoons have intensified [101–103]. In consequence, this may have also led to an increase of suitable habitats for P. arguta throughout the warmer periods of the Late Quaternary. Partly consistent with this, the ENM shows larger and more contiguous distributions at least in northern South Japan [i.e. Kii Peninsula and (south-)central Honshu] at both the LIG and the present compared to the LGM (Figure 6). Although suitable habitats may have shrunk in East China at the LIG (Figure 6b), and in both regions during the LGM (Figure 6c), we found no direct (EBSP) evidence for population decline (Additional file 10: Figure S3) and recent population bottlenecks (Additional file 1: Tables S1; Additional file 5: Tables S5). Nonetheless, the ENM suggests that potential glacial refugia in East China were more strongly affected by small-scale fragmentation compared to the situation in South Japan (Figure 6c) in general, and its southern parts (Shikoku, Kyushu) in particular, where such refugia might have even existed in nearby shelf areas. This inter-regional difference in glacial habitat suitability most likely reflects differences in topography and (micro-)climate, and should have differing consequences for population genetic structure and evolutionary history (see below).
Contrasting Late Quaternary evolutionary histories of P. argutain China and Japan
In var. sinensis, there is a marked subdivision between populations in nuclear genes (ITS, Tpi, nSSRs) (Table 1), along with high haplotypic and allelic diversity (Additional file 1: Table S1), and the same holds true for cpDNA . Together with our STRUCTURE result of five principal nSSR gene pools in this variety (Figure 5), all of these data indicate long-term population persistence and isolation over multiple glacial/interglacial cycles. With no indication of latitudinal range shifts, past population growth in var. sinensis (see above) may also reflect, at least in part, repeated downhill shifts in elevation range during glacials, perhaps by tracking favourable humidity conditions as imposed by the East Asian monsoon in areas of high relief . Over climatic cycling, such contiguous but localized range shifts along elevation gradients would have minimized bottlenecking and loss of genetic diversity, ultimately resulting in strong population subdivision .
In Japan, patterns of cpDNA diversity in var. arguta were previously interpreted to indicate southward retreat during glacials (to Kyushu and Shikoku) followed by expansion northward (up to central Honshu) during inter-/post-glacials . However, the current nuclear diversity and ENM data cast doubt on the validity of the ‘expansion-contraction’ (EC) model for var. arguta. For example, nuclear genetic diversity within var. arguta populations is more or less evenly spread throughout the distribution (Additional file 1: Table S1), and neither haplotypic (hS) diversity (ITS, Tpi) or allelic richness (Rs; nSSRs) show a significantly negative relationship with latitude (P = 0.33–0.58; data not shown), as would be expected under a scenario of south-to-north (re-)colonization [106, 107]. In addition, the ENM for the LGM indicates suitable habitats not only in the south (Kyushu, Shikoku, and nearby shelf areas) but also, more disjunctively, in the Kii Peninsula of south-central Honshu (Figure 6c). By contrast, large areas of contiguous suitable habitats were likely present at the LIG (Figure 6b), suggesting that opportunities for admixture occurred repeatedly during inter-/postglacial times, unless hampered by the spread of evergreen forest [1, 4].
Both of these latter predictions are basically bone out by the current nuclear data, even though each marker suggests a somewhat different historical scenario. Thus, the identification of two separate ITS lineages in Kyushu/Shikoku vs. central Honshu, with an apparent overlap of range in the Kii Peninsula (pop. J5; Figures 1a and 3), is strongly suggestive of a two-refugia scenario and a narrow contact zone between the two lineages. Interestingly, this putative contact region coincides with a well-known inter- and intraspecific suture zone of Japan’s temperate flora and fauna [14, 108, 109, 12]. On the other hand, the Tpi data are more consistent with a multiple-refugia scenario, that is, haplotypes are largely restricted to particular (sub)regions and/or populations, while still showing a relatively small amount of shared polymorphisms, mainly (but not exclusively) among adjacent sites (Figure 2a). Finally, the nSSRs reveal extensive admixture in most populations (J2–5), excepting those in Kyushu (J1) and central Honshu (J6, J7; Figure 5), suggesting that, in marked contrast to ITS, the contact area is much wider towards the south.
This discordance among patterns of genetic structure observed in the four genetic data sets makes it difficult to specify exactly which historical processes occurred in var. arguta. This discordance (even among nuclear genes and relative to cpDNA) is not readily explicable by a single difference in marker properties [e.g., mutation rates, transmission genetics, effective population sizes, concerted evolution (only ITS)], but likely results from a combination of marker features, with both incomplete lineage sorting and secondary admixture acting at different temporal and spatial scales . Taken together, we suggest that the genetic and ENM patterns found in var. arguta reflect a relatively ancient north–south vicariant event during glacials (sundering populations in central Honshu), followed by more recent climate-induced cycles of range contraction and expansion/admixture, with the latter producing more shallow patterns of divergence in more southerly areas (Kyushu, Shikoku, Kii).
This is the first study providing firm evidence that the ECS land bridge acted as a ‘filter’ during the last glacials in selectively preventing the dispersal of certain plant species of WTD forest, such as P. arguta, even though a near continuous belt of this forest biome is thought to have covered this land bridge during the LGM . Our data emphasize the species-specific recruitment and range expansion response of WTD forest-dwellers to the formation of the ECS land bridge, and caution against an uncritical interpretation of reconstructed palaeo-forest biome data as guides of past range dynamics of constituent species populations .
In addition, this multi-locus study strongly supports the two varieties of P. arguta in East China (var. sinensis) and South Japan (var. arguta) as genetically distinct units that evolved in strict allopatry since the mid-Pleistocene, which concurs with previous findings inferred from cpDNA alone . However, the combination of genetic structures from both nuclear and cytoplasmic loci can better depict the history of populations, demonstrating that (i) each lineage has undergone refugial isolation and divergence; and (ii) var. arguta likely experienced post-glacial admixture across a well-known suture zone. Yet, the extent of the species’ overall distribution does not seem to have greatly changed over the last glacial-interglacial cycles.
Availability of supporting data
The sequence alignments and phylogenetic trees were deposited in TreeBASE (http://treebase.org/treebase-web/home.html) under the submission number 15354. Sequence data used in this study have been deposited in GenBank (JQ978221––JQ978284, KF559183, KC853063).
Wu ZY, Wu SG: A proposal for a new floristic kingdom (realm) – the E. Asiatic kingdom, its delimitation and characteristics. Proceedings of the first international symposium on floristic characteristics and diversity of east Asian plants. Edited by: Zhang AL, Wu SG. 1996, Beijing, China: Springer-Verlag, 3-42.
Harrison SP, Yu G, Takahara H, Prentice IC: Palaeovegetation (Communications arising): diversity of temperate plants in east Asia. Nature. 2001, 413: 129-130. 10.1038/35093166.
Yu G, Chen X, Ni J, Cheddadi R, Guiot J, Han H, Harrison SP, Huang C, Ke M, Kong Z, Li S, Li W, Liew P, Liu G, Liu J, Liu Q, Liu K, Prentice IC, Qui W, Ren G, Song C, Sugita S, Sun X, Tang L, Van Campo E, Xia Y, Xu Q, Yan S, Yang X, Zhao J, Zheng Z: Palaeovegetation of China: a pollen data-based synthesis for the mid-Holocene and last glacial maximum. J Biogeogr. 2000, 27: 635-664. 10.1046/j.1365-2699.2000.00431.x.
Gotanda K, Yasuda Y: Spatial biome changes in southwestern Japan since the Last Glacial Maximum. Quat Int. 2008, 184: 84-93. 10.1016/j.quaint.2007.09.029.
Millien-Parra V, Jaeger JJ: Island biogeography of the Japanese terrestrial mammal assemblage: an example of a relict fauna. J Biogeogr. 1999, 26: 959-972. 10.1046/j.1365-2699.1999.00346.x.
Wang P: Response of Western Pacific marginal seas to glacial cycles: paleoceanographic and sedimentological features. Mar Geol. 1999, 156: 5-39. 10.1016/S0025-3227(98)00172-8.
Siddall M, Rohling EJ, Almogi-Labin A, Hemleben C, Meischner D, Schmelzer I, Smeed DA: Sea-level fluctuations during the last glacial cycle. Nature. 2003, 423: 853-858. 10.1038/nature01690.
Qian H, Ricklefs RE: Palaeovegetation - Diversity of temperate plants in east Asia - Reply. Nature. 2001, 413: 130-130. 10.1038/35093169.
Simpson GG: Mammals and land bridges. J Washington Acad Sci. 1940, 30: 137-163.
Lomolino MV, Riddle BR, Brown JH: Biogeography. 2006, Sunderland, MA: Sinauer Associates, Inc, 3
Chung CH: Vegetation response to climate change on Jeju Island, South Korea, during the last deglaciation based on pollen record. Geosci J. 2007, 11: 147-155. 10.1007/BF02913928.
Qi XS, Chen C, Comes HP, Sakaguchi S, Liu YH, Tanaka N, Sakio H, Qiu YX: Molecular data and ecological niche modelling reveal a highly dynamic evolutionary history of the East Asian Tertiary relict Cercidiphyllum (Cercidiphyllaceae). New Phytol. 2012, 196: 617-630. 10.1111/j.1469-8137.2012.04242.x.
Sakaguchi S, Qiu YX, Liu YH, Qi XS, Kim SH, Han J, Takeuchi Y, Worth JR, Yamasaki M, Sakurai S, Isagi Y: Climate oscillation during the Quaternary associated with landscape heterogeneity promoted allopatric lineage divergence of a temperate tree Kalopanax septemlobus (Araliaceae) in East Asia. Mol Ecol. 2012, 21: 3823-3838. 10.1111/j.1365-294X.2012.05652.x.
Li EX, Sun Y, Qiu YX, Guo JT, Comes HP, Fu CX: Phylogeography of two East Asian species in Croomia (Stemonaceae) inferred from chloroplast DNA and ISSR fingerprinting variation. Mol Phylogenet Evol. 2008, 49: 702-714. 10.1016/j.ympev.2008.09.012.
Qiu YX, Qi XS, Jin XF, Tao XY, Fu CX, Naiki A, Comes HP: Population genetic structure, phylogeography, and demographic history of Platycrater arguta (Hydrangeaceae) endemic to East China and South Japan, inferred from chloroplast DNA sequence variation. Taxon. 2009, 58: 1226-1241.
Qiu YX, Sun Y, Zong M, Zhang XP, Lee J, Murata J, Fu CX, Comes HP: Molecular phylogeography of East Asian Kirengeshoma (Hydrangeaceae) in relation to Quaternary climate change and landbridge configurations. New Phytol. 2009, 183: 480-495. 10.1111/j.1469-8137.2009.02876.x.
Wang JF, Gong X, Chiang YC, Kuroda C: Phylogenetic patterns and disjunct distribution in Ligularia hodgsonii Hook. (Asteraceae). J Biogeogr. 2013, 40: 1741-1754. 10.1111/jbi.12114.
Qiu YX, Fu CX, Comes HP: Plant molecular phylogeography in China and adjacent regions: tracing the genetic imprints of Quaternary climate and environmental change in the world's most diverse temperate flora. Mol Phylogenet Evol. 2011, 59: 225-244. 10.1016/j.ympev.2011.01.012.
Hein J, Schierup M, Wiuf C: Gene Genealogies, Variation and Evolution: A Primer in Coalescent Theory. 2004, Oxford: Oxford University Press
Hey J, Nielsen R: Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. persimilis. Genetics. 2004, 167: 747-760. 10.1534/genetics.103.024182.
Kuo C, Avise J: Phylogeographic breaks in low-dispersal species: the emergence of concordance across gene trees. Genetica. 2005, 124: 179-186. 10.1007/s10709-005-2095-y.
Katsuyama T: Platycrater arguta Siebold & Zucc. var. sinensis H. Hara (Saxifragaceae) collected in Fujian Province, China. J Japanese Bot. 1999, 74: 317-319.
Ohba H: Platycrater. Flora of Japan, Volume IIb. Edited by: Iwatsuki K, Boufford DE, Ohba H. 2001, Tokyo: Kodansha, 96-
Wei ZF, Bartholomew B: Platycrater. Flora of China. Edited by: Wu ZY, Raven PH. 2001, Beijing, China: Science Press and St. Louis, USA: Missouri Botanical Garden Press, 407-8
Hey J, Nielsen R: Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. P Natl Acad Sci USA. 2007, 104: 2785-2790. 10.1073/pnas.0611164104.
Kukla GJ, Bender ML, de Beaulieu J-L, Bond G, Broecker WS, Cleveringa P, Gavin JE, Herbert TD, Imbrie J, Jouzel J, Keigwin LD, Knudsen K-L, McManus JF, Merkt J, Muhs DR, Muller H, Poore RZ, Porter SC, Seret G, Shackleton NJ, Turner C, Tzedakis PC, Winograd IJ: Last interglacial climates. Quaternary Res. 2002, 58: 2-13. 10.1006/qres.2001.2316.
Soltis DE, Xiang QY, Hufford L: Relationships and evolution of Hydrangeaceae based on rbcL sequence data. Am J Bot. 1995, 82: 504-514. 10.2307/2445698.
Hufford L, Moody ML, Soltis DE: A phylogenetic analysis of Hydrangeaceae based on sequences of the plastid gene matK and morphological data. Int J Plant Sci. 2001, 158: 652-672.
Doyle J: DNA protocols for plants – CTAB total DNA isolation. Molecular Techniques in Taxonomy. Edited by: Hewitt GM, Johnston A. 1991, Berlin: Springer, 283-293.
Strand AE, Leebens-Mack J, Milligan BG: Nuclear DNA-based markers for plant evolutionary biology. Mol Ecol. 1997, 6: 113-118. 10.1046/j.1365-294X.1997.00153.x.
Malcomber ST: Phylogeny of Gaertnera Lam. (Rubiaceae) based on multiple DNA markers: evidence of a rapid radiation in a widespread, morphologically diverse genus. Evolution. 2002, 56: 42-57.
Zou XH, Zhang FM, Zhang JG, Zhang LL, Tang L, Wang J, Sang T, Ge S: Analysis of 142 genes resolves the rapid diversification of the rice genus. Genome Biol. 2008, 9: R49-10.1186/gb-2008-9-3-r49.
Clark AG: Inference of haplotypes from PCR-amplified samples of diploid populations. Mol Biol Evol. 1990, 7: 111-122.
Liao PC, Kuo DC, Lin CC, Ho KC, Lin TP, Hwang SY: Historical spatial range expansion and a very recent bottleneck of Cinnamomum kanehirae Hay. (Lauraceae) in Taiwan inferred from nuclear genes. BMC Evol Biol. 2010, 10: 124-10.1186/1471-2148-10-124.
Drummond AJ, Ashton B, Cheung M, Heled J, Kearse M, Moir R, Stones-Havas S, Thierer T, Wilson A: Geneious v4. 8. 2009, http://www.geneious.com,
Qi XS, Yuan N, Qiu YX: Development of 12 microsatellite markers for Platycrater arguta (Hydrangeaceae) endemic to East Asia. Am J Bot. 2012, 99: e304-e306. 10.3732/ajb.1100582.
Tajima F: Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics. 1989, 123: 585-595.
Fu YX, Li WH: Statistical tests of neutrality of mutations. Genetics. 1993, 133: 693-709.
Fu YX: Statistical tests of neutrality of mutations against population growth, hitchhiking and background selection. Genetics. 1997, 147: 915-925.
Ramos-Onsins SE, Rozas J: Statistical properties of new neutrality tests against population growth. Mol Biol Evol. 2002, 19: 2092-2100. 10.1093/oxfordjournals.molbev.a004034.
Hudson RR, Kaplan N: Statistical properties of the number of recombination events in the history of a sample of DNA sequences. Genetics. 1985, 111: 147-165.
Rozas J, Sánchez-De I, Barrio JC, Messeguer X, Rozas R: DnaSP, DNA polymorphism analyses by the coalescent and other methods. Bioinformatics. 2003, 19: 2496-2497. 10.1093/bioinformatics/btg359.
Excoffier L, Laval G, Schneider S: Arlequin ver. 3.1: an integrated software package for population genetics data analysis. 2006, Switzerland: Computational and Molecular Population Genetics Laboratory, Institute of Zoology, University of Bern
Swofford DL: PAUP*: Phylogenetic Analysis Using Parsimony (* and Other Methods), Version 4.0b10. 2002, Sunderland, MA: Sinauer Associates
Stamatakis A, Hoover P, Rougemont J: A rapid bootstrap algorithm for the RAxML web-servers. Syst Biol. 2008, 75: 758-771.
Clement M, Posada D, Crandall KA: TCS: a computer program to estimate gene genealogies. Mol Ecol. 2000, 9: 1657-1659. 10.1046/j.1365-294x.2000.01020.x.
Pons O, Petit RJ: Measuring and testing genetic differentiation with ordered versus unordered alleles. Genetics. 1996, 144: 1237-1245.
Excoffier L, Smouse PE, Quattro JM: Analysis of molecular variance inferred from metric distances among DNA haplotypes: application to human mitochondrial DNA restriction data. Genetics. 1992, 131: 479-491.
Goudet J: FSTAT, version 184.108.40.206. A program to estimate and test gene diversities and fixation indices. 2001, http://www2.unil.ch/popgen/softwares/fstat.htm,
Slatkin M: A measure of population subdivision based on microsatellite allele frequencies. Genetics. 1995, 139: 457-462.
Kimura M, Ohta T: Stepwise mutation model and distribution of allelic frequencies in a finite population. P Natl Acad Sci USA. 1978, 75: 2868-2872. 10.1073/pnas.75.6.2868.
Pritchard JK, Wen X, Falush D: Documentation for Structure software: Version 2.3. 2009, http://pritch.bsd.uchicago.edu/structure_software/release_versions/v2.3.2/structure_doc.pdf,
Falush D, Stephens M, Pritchard JK: Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics. 2003, 164: 1567-1587.
Evanno G, Regnaut S, Goudet J: Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol Ecol. 2005, 14: 2611-2620. 10.1111/j.1365-294X.2005.02553.x.
Piry S, Luikart G, Cornuet JM: Bottleneck: a computer program for detecting recent reductions in the effective size using allele frequency data. J Hered. 1999, 90: 502-503. 10.1093/jhered/90.4.502.
Woerner AE, Cox MP, Hammer MF: Recombination-filtered genomic datasets by information maximization. Bioinformatics. 2007, 23: 1851-1853. 10.1093/bioinformatics/btm253.
Kay KM, Whittall JB, Hodges SA: A survey of nuclear ribosomal internal transcribed spacer substitution rates across angiosperms: an approximate molecular clock with life history effects. BMC Evol Biol. 2006, 6: 36-10.1186/1471-2148-6-36.
Wolfe KH, Li WH, Sharp PM: Rates of nucleotide substitution vary greatly among plant mitochondrial, chloroplast, and nuclear DNAs. P Natl Acad Sci USA. 1987, 84: 9054-9058. 10.1073/pnas.84.24.9054.
Ossowski S, Schneeberger K, Lucas-Lledo J, Warthmann N, Clark R, Shaw R, Weigel D, Lynch M: The rate and molecular spectrum of spontaneous mutations Arabidopsis thaliana. Science. 2010, 327: 92-94. 10.1126/science.1180677.
Muse SV: Examining rates and patterns of nucleotide substitution in plants. Plant Mol Biol. 2000, 42: 25-43. 10.1023/A:1006319803002.
Boys J, Cherry M, Dayanandan S: Microsatellite analysis reveals genetically distinct populations of red pine (Pinus resinosa, Pinaceae). Am J Bot. 2005, 92: 833-841. 10.3732/ajb.92.5.833.
Bulut Z, Mccormick CR, Gopurenko D, Williams RN, Bos DH, Dewoody JA: Microsatellite mutation rates in the eastern tiger salamander (Ambystoma tigrinum tigrinum) differ 10-fold across loci. Genetica. 2009, 136: 501-504. 10.1007/s10709-008-9341-z.
Udupa SM, Baum M: High mutation rate and mutational bias at (TAA)n microsatellite loci in chickpea (Cicer arietinum L.). Mol Genet Genomics. 2001, 265: 1097-1103. 10.1007/s004380100508.
Thuillet AC, Bru D, David J, Roumet P, Santoni S, Sourdille P, Bataillon T: Direct estimation of mutation rate for 10 microsatellite loci in durum wheat, Triticum turgidum (L.) Thell. ssp durum desf. Mol Biol Evol. 2002, 19: 122-125. 10.1093/oxfordjournals.molbev.a003977.
Vigouroux Y, Jaqueth JS, Matsuoka Y, Smith OS, Beavis WD, Smith JS, Doebley J: Rate and pattern of mutation at microsatellite loci in maize. Mol Biol Evol. 2002, 19: 1251-1260. 10.1093/oxfordjournals.molbev.a004186.
Ally D, Ritland K, Otto SP: Can clone size serve as a proxy for clone age? An exploration using microsatellite divergence in Populus tremuloides. Mol Ecol. 2008, 17: 4897-4911. 10.1111/j.1365-294X.2008.03962.x.
Marriage TN, Hudman S, Mort ME, Orive ME, Shaw RG, Kelly JK: Direct estimation of the mutation rate at dinucleotide microsatellite loci in Arabidopsis thaliana (Brassicaceae). Heredity. 2009, 103: 310-317. 10.1038/hdy.2009.67.
Kondo T, Crisp MD, Linde C, Bowman DMJS, Kawamura K, Kaneko S, Isagi Y: Not an ancient relic: the endemic Livistona palms of arid central Australia could have been introduced by humans. Proc R Soc B. 2012, 279: 2652-2661. 10.1098/rspb.2012.0103.
Delplancke M, Alvarez N, Benoit L, Espíndola A, I Joly H, Neuenschwander S, Arrigo N: Evolutionary history of almond tree domestication in the Mediterranean basin. Mol Ecol. 2013, 22: 1092-1104. 10.1111/mec.12129.
Schaal BA, Hayworth DA, Olsen KM, Rauscher JT, Smith WA: Phylogeographic studies in plants: problems and prospects. Mol Ecol. 1998, 7: 465-474. 10.1046/j.1365-294x.1998.00318.x.
Soltis DE, Soltis PS, Ness BD: Maternal inheritance of the chloroplast genome in Heuchera and Tolmiea (Saxifragaceae). J Hered. 1990, 81: 168-170.
Heled J, Drummond AJ: Bayesian inference of population size history from multiple loci. BMC Evol Biol. 2008, 8: 289-10.1186/1471-2148-8-289.
Ho SYW, Shapiro B: Skyline-plot methods for estimating demographic history from nucleotide sequences. Mol Ecol Resour. 2011, 11: 423-434. 10.1111/j.1755-0998.2011.02988.x.
Posada D: jModelTest: phylogenetic model averaging. Mol Biol Evol. 2008, 25: 1253-1256. 10.1093/molbev/msn083.
Rambaut A, Drummond AJ: Tracer v1.5. 2009, http://tree.bio.ed.ac.uk/software/tracer/,
Phillips SJ, Anderson RP, Schapire RE: Maximum entropy modeling of species geographic distributions. Ecol Model. 2006, 190: 231-259. 10.1016/j.ecolmodel.2005.03.026.
Prefecture M: Red Data Book Mie, Plant and Fungi. 2005, Japan, Mie, Mie Kankyo Hozen Jigyodan
Aichi Prefecture Government: Red Data Book Aichi. 2009, Japan, Nagoya: Aichi Prefecture
Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A: Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005, 25: 1965-1978. 10.1002/joc.1276.
Hasumi H, Emori S: K-1 coupled GCM (MIROC) description, K-1 Technical Report 1. 2004, Japan, Tokyo: Center for Climate System Research, University of Tokyo, http://ccsr.aori.u-tokyo.ac.jp/~hasumi/miroc_description.pdf,
Otto-Bliesner BL, Schneider R, Brady EC, Kucera M, Abe-Ouchi A, Bard E, Braconnot P, Crucifix M, Hewitt CD, Kageyama M, Marti O, Paul A: A comparison of PMIP2 model simulations and the MARGO proxy reconstruction for tropical sea surface temperatures at last glacial maximum. Clim Dynam. 2009, 32: 799-815. 10.1007/s00382-008-0509-0.
Flanders J, Wei L, Rossiter SJ, Zhang SY: Identifying the effects of the Pleistocene on the greater horseshoe bat, Rhinolophus ferrumequinum, in East Asia using ecological niche modelling and phylogenetic analyses. J Biogeogr. 2011, 38: 439-452. 10.1111/j.1365-2699.2010.02411.x.
Qu Y, Zhang R, Quan Q, Song G, Li SH, Lei F: Incomplete lineage sorting or secondary admixture: disentangling historical divergence from recent gene flow in the Vinous-throated parrotbill (Paradoxornis webbianus). Mol Ecol. 2012, 21: 6117-6133. 10.1111/mec.12080.
Fawcett T: An introduction to ROC analysis. Pattern Recognit Lett. 2006, 27: 861-874. 10.1016/j.patrec.2005.10.010.
Sakaguchi S, Sakurai S, Yamasaki M, Isagi Y: How did the exposed seafloor function in postglacial northward range expansion of Kalopanax septemlobus? Evidence from ecological niche modelling. Ecol Res. 2010, 25: 1183-1195. 10.1007/s11284-010-0743-x.
Kimura M: Quaternary paleogeography of the Ryukyu Arc. J Geogr. 1996, 105: 259-285. 10.5026/jgeography.105.3_259.
Kimura M: Land connections between Eurasian continent and Japanese Islands – related to human migration. Migration and Diffusion. 2003, 4: 14-33.
Minato M, Ijiri S: The Japanese Archipelago. 1984, Tokyo: Iwanami-shoten, 3
Dobson M, Kawamura Y: Origin of the Japanese land mammal fauna: allocation of extant species to historically-based categories. Quaternary Res. 1998, 37: 385-395. 10.4116/jaqua.37.385.
Keally CT: Japanese Pleistocene landbridges and the earliest watercraft. Japanese archaeology. 2005, http://www.t-net.ne.jp/~keally/MiddlePalaeol/landbridges.html,
North Greenland Ice Core Project Members: High resolution climate record of the northern hemisphere reaching into the Last Interglacial period. Nature. 2004, 431: 147-151. 10.1038/nature02805.
Hao QZ, Guo ZT, Qiao YS, Xu B, Oldfield F: Geochemical evidence for the provenance of middle Pleistocene loess deposits in southern China. Quat Sci Rev. 2010, 29: 23-24.
Ruddiman WF, Raymo ME, Martinson DG, Clement BC, Backman J: Pleistocene Evolution: Northern Hemisphere ice sheets and North Atlantic Ocean Climate. Paleoceanography. 1989, 4: 353-412. 10.1029/PA004i004p00353.
Iida S, Nakashizuka T: Spatial and temporal dispersal of Kalopanax pictus seeds in a temperate deciduous forest, central Japan. Plant Ecol. 1998, 135: 243-248. 10.1023/A:1009779921463.
Sakio H, Tamura T: Ecology of Riparian Forests in Japan: Disturbance, Life History and Regeneration. 2008, Tokyo: Springer
Eriksson O, Ehrlén J: Landscape fragmentation and the viability of plant populations. Integrating Ecology and Evolution in a Spatial Context. Edited by: Silvertown J, Antonovics J. 2001, Oxford: Blackwell Science, 157-175.
Suzuki W, Osumi K, Masaki T, Takahashi K, Daimaru H, Hoshizaki K: Disturbance regimes and community structures of a riparian and an adjacent stand in the Kanumazawa Riparian Research Forest, northern Japan. For Ecol Manage. 2002, 157: 285-301. 10.1016/S0378-1127(00)00667-8.
Rodríguez-Sánchez F, Hampe A, Jordano P, Arroyo J: Past tree range dynamics in the Iberian Peninsula inferred through phylogeography and palaeodistribution modelling: a review. Rev Palaeobot Palynol. 2010, 162: 507-521. 10.1016/j.revpalbo.2010.03.008.
Won Y, Sivasundar A, Wang Y, Hey J: On the origin of Lake Malawi cichlid species: a population genetic analysis of divergence. P Natl Acad Sci USA. 2010, 102: 6581-6586.
Yi CL, Cui ZJ, Xiong HG: Numerical periods of Quaternary glaciations in China. Quat Sci. 2010, 25: 609-619.
Zhang R, Liu XD: The effects of tectonic uplift on the evolution of Asian summer monsoon climate since Pliocene. Chinese J Geophys. 2010, 53: 948-960. 10.1002/cjg2.1565.
Ao H, Dekkers MJ, Xiao G, Yang X, Qin L, Liu X, Qiang X, Chang H, Zhao H: Different orbital rhythms in the Asian summer monsoon records from North and South China during the Pleistocene. Glob Planet Change. 2012, 80: 51-60.
Han WX, Fang XM, Berger A: Tibet forcing of mid-Pleistocene synchronous enhancement of East Asian winter and summer monsoons revealed by Chinese loess record. Quaternary Res. 2003, 78: 174-184.
Liu H, Xing Q, Ji Z, Xu L, Tian Y: An outline of Quaternary development of Fagus forest in China: palynological and ecological perspectives. Flora. 2003, 198: 249-259. 10.1078/0367-2530-00098.
Hewitt GM: Genetic consequences of climatic oscillations in the Quaternary. Philos Trans R Soc Lond B Biol Sci. 2004, 359: 183-195. 10.1098/rstb.2003.1388.
Excoffier L, Foll M, Petit RJ: Genetic consequences of range expansions. Annu Rev of Ecol Evol Syst. 2009, 40: 481-501. 10.1146/annurev.ecolsys.39.110707.173414.
Kerdelhué C, Zane L, Simonato M, Salvato P, Rousselet J, Roques A, Battisti A: Quaternary history and contemporary patterns in a currently expanding species. BMC Evol Biol. 2009, 9: 14-10.1186/1471-2148-9-14.
Aoki K, Kato M, Murakami N: Phylogeographical patterns of a generalist acorn weevil: insight into the biogeographical history of broadleaved deciduous and evergreen forests. BMC Evol Biol. 2009, 9: 103-10.1186/1471-2148-9-103.
Kikuchi R, Jae-Hong P, Takahashi H, Maki M: Disjunct distribution of chloroplast DNA haplotypes in the understory perennial Veratrum album ssp. oxysepalum (Melanthiaceae) in Japan as a result of ancient introgression. New Phytol. 2010, 188: 879-891. 10.1111/j.1469-8137.2010.03398.x.
This research was supported by the National Science Foundation of China (grant nos. 31370241, 31170200), the Open Research Foundation of LSEB (State Key Laboratory of Systematic and Evolutionary Botany, Institute of Botany, Chinese Academy of Sciences), and the Zhejiang Provincial Natural Science Foundation of China (Grant No. LR12C02001). We are grateful to Dr. Richard G. Hodel (University of Florida) for assistance with the molecular data analysis, and Prof. Richard J. Abbott (University of St Andrews), Dr. Kenneth M. Olsen (Washington University in St. Louis), and three anonymous reviewers for their insightful comments and suggestions to improve the manuscript.
The authors declare that they have no competing interests.
YXQ conceived the ideas; SS and XSQ contributed to the sampling; XSQ and YN collected and analysed the molecular data and drafted the manuscript. SS performed the ENM. The manuscript was written by HPC and YXQ. All authors read and approved the final manuscript.