Interspecific Divergence of Two Sinalliaria (Brassicaceae) Species in Eastern China

How endemic species originated in eastern Asia has interested botanists for a long time. In this study, we combined experimental and computational modeling approaches to examine the morphological and genetic divergence and reproductive isolation of two tentative species of Sinalliaria (Brassicaceae) endemic to eastern China, S. limprichtiana and S. grandifolia. Most of the examined morphological characters (including hairs of leaf blades and stems, corolla length and width, and flower stalk length) were well-delineated between two species at the same ploidy level, and there was clear evidence of reproductive isolation between them (mainly due to post-pollination barriers) in the common garden environment. There were also strong and consistent divergences in the population genetic data. Coalescent simulations based on sequence variation of the nuclear genes suggest that interspecific divergence began during the Pleistocene when the climate oscillated in eastern Asia. Gene flow between two species appears to have been very limited and asymmetrical. Our results suggested that both species are well-differentiated and that the fast divergence between them might have been together shaped by both stochastic processes and habitat selection pressures.


INTRODUCTION
For several decades, biogeographic researchers have considered the reasons for the high species diversity of plants and numerous endemics in eastern Asia, especially compared with eastern North America, which has a similar size, climate, and floristic components (e.g., Tiffney, 1985;Axelrod et al., 1998;Qian and Ricklefs, 2000). The high diversity has been attributed to a combination of two major processes (Liu, 1988;Qian and Ricklefs, 2000;Qian et al., 2003). First, numerous ancient genera and species may have persisted due to geographical and ecological heterogeneity that could buffered the organisms from small-scale extinctions caused by the climate oscillations in eastern Asia, in conjunction with the lack of Pleistocene glaciation and subsequent large-scale extinctions. Second, many young species (for example, four Dysosma species, D. pleiantha, D. versipellis, D. difformis, and D. majoensis) have originated through fast divergence probably driven by Pleistocene climate oscillations in conjunction with physiographic heterogeneity (Ikeda et al., 2012;Han et al., 2016;Wang et al., 2017). Both hypotheses have been partly tested based on biogeographic analyses of a few genera and phylogeographic examinations of some widely distributed species (see summaries by Qiu et al., 2011;Liu et al., 2012). The direct evidence for such hypotheses should derive from the times and processes of the well-delimitated sister species that may have occurred during the Pleistocene (Zanella et al., 2016;Wang et al., 2017). Although Ikeda et al. (2012), Han et al. (2016), and Wang et al. (2017) have revealed the Pleistocene origins of a pair of endemic alpine species, two coastal wild radish lineages and four understory herbs based on population genetic data, it remains unknown whether these species or lineages were morphologically differentiated with stable distinct characters and reproductively isolated as "good species" based on the integrate species concept (Liu, 2016).
Here we focus on the two only known annual species of the genus Sinalliaria, S. limprichtiana, and S. grandifolia (Brassicaceae). Both species are endemic to eastern Asia (Zhou et al., 2014;Hu et al., 2015), and narrowly distributed in warmtemperate evergreen forests or forest edges in eastern China (Figure 1) (Al-Shehbaz and Guang, 2000;Zhou et al., 2014). These two taxa were treated as two subspecies by some authors (Zhou et al., 2014) and they were found to have clear genetic distinctions based on sequence variations of the chloroplast DNAs and the limited number of the sampled individuals for each species (Hu et al., 2015(Hu et al., , 2016. Nevertheless, whether morphological differentiation remains stable and reproductive isolation has been well-established between them still need to be confirmed. Moreover, the divergence time and demographic history of two species remains to be elucidated using population genetic data at multiple nuclear loci. In the present study, we first examined the ploidy level, morphological differentiation, and strength of two species' preand post-pollination reproductive isolation in a common garden when ecological factors were totally excluded. Then we measured the genetic divergence at the population level between them, according to sequence variation of both 12 Simple Sequence Repeat (SSR) loci and nine unlinked, low-copy nuclear genes. Finally, we quantified changes in their effective population sizes and amounts of gene flow during their speciation and established whether the interspecific divergence was related to Pleistocene climate changes with the relatively young age.

Common Garden Experiments
To assess phenotypic differences between two species, we used three S. limprichtiana populations and two S. grandifolia populations (Figure 1, Table S1). For each population, we selected five individuals. For each individual, we used 10 mature seeds collected from the field. The seeds were kept at 4 • C for 3 weeks and then sown in Petri dishes containing water-soaked filter paper at 25 • C for 3-7 days. Five to 10 resulting seedlings for each field individual were transplanted to 42 mm-diameter pots and placed in a growth chamber at 25 • C during vegetative growth for 1 week. After forming three leaves, all seedlings were moved from pots to an experimental field in the Sichuan University campus in October 2015. We measured germination rates each month. All seedlings started to flower in March and set fruit in April 2016.

Determination of Ploidy Level by Flow Cytometry
We estimated the amounts of nuclear DNA in each species' genome following the protocol described by Dolezel and Bartos (2005). Briefly, we isolated the nuclei from fresh leaves of each three individuals, from three S. limprichtiana populations and two S. grandifolia populations mentioned above, using Otto I solution, and subsequently suspended them with a mixture of Otto I and Otto II solution (1:4) in phosphate/citric acid buffer, pH 7.3. After 30 min staining with propidium iodide, nuclear fluorescence was measured using a Partec CA-II flow cytometer, with Vigna radiata (2C = 1.40 pg DNA and 1 pg DNA = 9.78 × 10 8 base pairs) as a calibration standard (Dolezel and Bartos, 2005). We used the coefficient of variation (CV = standard deviation/peak mean × 100%) of DNA peaks to quantify the precision of individual measurements with the value of CV < 3% calculated by FlowJo software.

Morphological Differentiation
To minimize the morphological variations caused by ecological factors, we randomly selected five descendants cultivated from each of the five individuals of each population of S. limprichtiana and S. grandifolia in the common garden at the full flowering stage. Firstly, we recorded the presence or absence of hairs on their leaf blades and stems. Then we examined and/or measured the characters of five flowers and sets of three basal leaves, central-stem leaves, and upper-stem leaves produced by each individual. These traits included six floral characters (corolla length and width, flower stalk length, calyx length, number of inflorescence branches, and number of flowers; Figure S3), and nine leaf characters (length, width, and petiole lengths of each set of leaves) ( Figure 2B, Figure S2). We used SPSS version 23.0 software to calculate means and standard deviations of these variables, and apply one-way analysis of variance (ANOVA) to assess the variability of each measured trait (Day and Quinn, 1989;Schürch et al., 2000), deeming differences to be significant if p < 0.01. Finally, we used principal co-ordinate analysis (PCoA; Bitner-Mathé and Klaczko, 1999), implemented in the Multivariate Statistics Package (MVSP) version 3.2 (http://www.kovcomp.com/), to assess the distinctiveness of two species, based on all examined and measured traits.

Pre-and Post-pollination Reproductive Isolation
We evaluated the strength of pre-pollination isolation between the species by examining the overlap of their flowering periods, defined as the time (about 56 days) when an individual had at least one open flower, with pollen and/or a wet receptive stigma. We daily observed and recorded the number of these "active" flowers on five individuals cultivated from seeds collected from each of five individuals from each population.
We also evaluated the strength of the species' postpollination reproductive isolation by comparing fruit set arising from inter-specific crossing, intra-specific (but inter-population) crossing, and spontaneous self-pollination. For all crossing experiments, we removed anthers before flowering and then bagged flowers. At the flowering stage, pollen was removed from the anthers of the designated "male" flowers and transferred to the receptive stigmatic cavities of the designated "female" flowers using a tooth pick with marker tags. All artificially pollinated flowers were re-bagged. The flowers assigned to the self-pollination treatment were continuously bagged and not otherwise manipulated. For each individual, only one flower was chosen for each of three above-mentioned experiments. A month after the pollination treatment we collected siliques and recorded the fruit set (defined as the number of full fruits divided by sum of fruits and undeveloped fruits). Differences in fruit set between the crossing types (and selfing treatment) were evaluated by ANOVA, followed by Student's two-tailed t-test for paired comparisons, with SPSS version 23.0 software, deeming differences to be significant if p < 0.05.

Extraction of Total DNA
Using a modified CTAB method (Doyle, 1987), we extracted total genomic DNA from ∼100 mg portions of dry leaf material collected from five and four natural populations of S. grandifolia and S. limprichtiana, respectively (Table S2). One population of S. limprichtiana is distributed far away from the other four (Figure 1), indicating a wide disjunction. However, there is no specimen record and our field explorations also failed to find its occurrence in other regions.

"SSR" Primer Development
First, we extracted RNA from fresh leaves of one individual of each species using TRIZOL reagent (Sigma Aldrich). Then we sequenced them to develop primers to amplify target sequences. Trinity (http://trinityrnaseq.sourceforge. net/) was used to assemble the RNA-Seq data, which retrieves large fractions of transcripts, especially abundant transcripts (Grabherr et al., 2011). To avoid redundancy in the datasets, we used CD-HIT (http://weizhong-lab.ucsd.edu/cd-hit) to obtain unigene sequences. The Microsatellite Identification Tool (MISA, http://pgrc.ipk-gatersleben.de/misa/) was used to identify microsatellites in the unigenes. Twelve microsatellite loci containing three-nucleotide motifs with at least five repeats were selected (Table 1). Using Primer version 3.0 software, primers with complementary sequences to flanking regions of these SSRs were then designed, varying in length from 20 to 23 base pairs (bp; close to the reportedly optimal length of 20 bp), with GC contents varying between 45 and 65% (close to the reported optimal 50%) (Qi et al., 2016). We also developed primers (Table 1) for nine low-copy nuclear genes with Primer version 6.0 software, using homologous EST genes of two species (Wang et al., 2011). Functional annotation of the corresponding genes was performed by BLAST searches of the NCBI database.

Genotyping Polymorphisms of SSR Loci
In total, 38 individuals representing four populations of S. limprichtiana and 36 individuals representing five S. grandifolia populations were genotyped using the 12 pairs of EST-SSR primers ( Table 1). After a 5-min denaturation at 95 • C, PCR was performed for 35-36 cycles (95 • C for 45 s, 55 • C for 40 s, and 72 • C for 80 s), with a 7-10 min final extension at 72 • C. The resulting PCR products were differentiated and genotyped using an ABI 3730 xL automated sequencer.

Sequencing Low-Copy Nuclear Genes
We amplified nine low-copy nuclear genes (Table 1) in 19 individuals representing four S. limprichtiana populations and 20 individuals representing four S. grandifolia populations. Amplifications were done by initial denaturation (94 • C, 5 min), 35-36 cycles of denaturation (94 • C, 30 s), annealing (55 • C, 45 s) and extension (72 • C, 50-60 s), and final extension (72 • C, 7-10 min). The resulting products were purified using a TIANquick Midi Purification Kit (Tiangen Biotech, Beijing, China), then sequenced with primers covering the whole PCR segments using an ABI Prism Bigdye TM Terminator Cycle Sequencing Ready Reaction Kit, and visualized using an ABI 3730 xL sequencer. Sequences were aligned using Clustal X (Thompson et al., 1997) with default parameters and curated manually. The boundaries of each dataset were delimited by comparison with sequences of the sister species. The obtained sequences have been deposited in GenBank under accession numbers KY689753-KY689824.

Genetic Diversity and Divergence at SSR Loci
We calculated genetic diversity parameters for the 12 SSR loci in both species, including information index (I), expected heterozygosity (H e ), observed heterozygosity (H o ), unbiased expected heterozygosity (uH), fixation index (F) values and F-statistics using GenAlex version 6.5 (Peakall and Smouse, 2012). We also examined the genetic variance within and among populations and species by analysis of molecular variance (AMOVA), using ARLEQUIN version 3.5.1.2 (Excoffier and Lischer, 2010).
We examined the population structure of two species, based on variation at the 12 SSR loci, using STRUCTURE version 2.3 software (Hubisz et al., 2009). Twenty independent runs with the number of clusters (K) varied from 1 to 20, in each case with 500,000 iterations and a 50,000 burn-in period (based on the initial comparison) using prior population information and an admixture model assuming correlated allele frequencies among clusters were applied to infer population structure according to Pritchard et al. (2000). K values based on the rate of change in the log probability of data between successive K values were estimated and the most likely cluster was chosen (Evanno et al., 2005). Then, we used PCoA to explore the internal structure of the presence/absence matrix for the entire data set using GenAlEx version 6.5 (Krzanowski, 2000). In PCoA multivariate statistics are used to depict the genetic structure, without relying on Frontiers in Plant Science | www.frontiersin.org the population genetics assumptions underlying STRUCTURE (Engelhardt and Stephens, 2010).

Genetic Divergence at Nine Nuclear Low-Copy Genes
We used PHASE version 2.1.1 to determine the haplotype phase of the directly sequenced PCR products for each nuclear locus for each individual (Stephens et al., 2001). We aligned the multiple sequences for each gene by MEGA version 6.0, then constructed haplotype networks, representing unique alleles separated by mutational steps, using the median-joining method (Bandelt et al., 1999) implemented in NETWORK version 5.0 (http:// www.fluxus-technology.com/sharenet.htm). Next, we applied DnaSP version 5.10 to calculate basic parameters (Librado and Rozas, 2009) for each gene, including the number of synonymous substitutions, the number of haplotypes (N h ), haplotype diversity (H e ), minimum number of recombination events (R m ), average number of nucleotide polymorphisms per site (π) between two sequences in a sample (Nei, 1987), and θ w = 4N e µ, where N e is the effective population size and µ is the mutation rate per site per generation) based on the number of segregating sites (Watterson, 1975). We also calculated Tajima's D (Tajima, 1989), Fu and Li's D * and F * statistics, using DnaSP version 5.10 (Fu and Li, 1993), to test the neutrality of each gene's mutations. Moreover, we constructed their phylogenetic relationships using the neighbor-joining approach with 500 bootstrap replicates, as implemented in MEGA version 6.0. The same as the SSRs dataset, we also used ARLEQUIN version 3.5.1.2 and STRUCTURE version 2.3 to evaluate the partitioning of genetic diveristy and population structure based on the nuclear genes.

Coalescence-Based Analyses of Divergence Times
We examined the species' speciation history using the following demographic parameters derived from analyses of the low-copy nuclear genes: population-split time (T), bidirectional migration rate m (m GL and m LG represent the migration rate from the S. limprichtiana to the S. grandifolia and from the S. grandifolia to the S. limprichtiana, respectively), and effective sizes of the ancestral (θ A ) and descendant populations (θ L and θ G ). We calculated all parameters using the IM model implemented in the IMa2 program. The IM model assumes no selection, no recombination within loci, no substantial population structure, and random mating in ancestral and descendent populations (Hey andNielsen, 2004, 2007). These assumptions are unlikely to be completely met in real populations. However, parameter estimates of the IM model have been found to be robust to high levels of population structure and recombination (Carstens and Knowles, 2007). We used the program to estimate parameters of the IM model, and assessed the posterior probability densities for model parameters using Markov Chain Monte Carlo (MCMC) methods (Hey and Nielsen, 2007;Hey, 2010). After several preliminary runs to optimize prior boundaries for the six primary demographic parameters (T, m LG , m GL , θ A , θ L, and θ G ), we set the length of the burn-in period of 1 × 10 5 iterations, followed by 1 × 10 6 MCMC steps and Metropolis-coupling of 120 independent heated chains. We assessed the MCMC mixing properties based on the effective sample sizes (ESS), swapping rates between successive MCMC chains, and trend-line plots of the parameters. We repeated the well-mixed runs to obtain reproducible results when three independent runs produced similar posterior distributions (Ikeda et al., 2009;Wang et al., 2013). All demographic parameters have to be scaled using the geometric mean of the mutation rate per year per locus during IMa2 calculations, so selection of appropriate mutation rates is crucial for estimates of almost all parameters (Strasburg and Rieseberg, 2010). The chalcone synthase (CHS) gene in Brassicaceae reportedly has a synonymous mutation rate (µ CHS ) of 1.5 × 10 −8 substitutions per site per year (Koch et al., 2000). Therefore, following Ikeda et al. (2009) and Wang et al. (2013), we estimated the mutation rate (µ) at silent sites of each locus using the equation µ = µ CHS × K Total /K S × L, where K Total /K S is the ratio of the average total genetic divergence to the average synonymous genetic divergence, and L is the length of the focal locus. The estimated geometric average of K Total /K S over all loci was 0.8353. We therefore used the derived geometric mean of 5.2947 × 10 −6 substitutions per locus per year (yr) to scale the demographic parameters.

Interspecific Differentiation in Genome Sizes and Morphology
We estimated genome size and ploidy level of each species by flow cytometry (FCM). As shown in Table S3, the flow cytometry results yielded 2C DNA values for the three sampled individuals from different populations of S. limprichtiana ranging from 1.54 to 1.70 pg (mean and standard deviation: 1.604 ± 0.09 pg), and corresponding values for S. grandifolia are 1.55 to 1.71 pg (1.647 ± 0.08 pg) (Table S3). Both species are at the same ploidy level with very similar genome sizes ( Figure S1).
There were no clear differences in overall germination and survival rates between two species, but we observed clear morphological differences. First, hairs were present on calyxes, leaf blades and stems of S. limprichtiana, but not S. grandifolia (Figure 2A). Second, S. limprichtiana has shorter and thinner basal and stem leaves than S. grandifolia (Figure 2B, Figure S2). Third, S. limprichtiana produces more flowers and inflorescence branches than S. grandifolia (Figures S3E,F). Finally, S. limprichtiana has shorter flower corollas than S. grandifolia ( Figure S3A). PCoA analyses of the 15 examined traits clearly distinguished the species (Figure 2C).

Pre-and Post-pollination Reproductive Isolation
The flowering times of the artificially constructed populations of each species in the common garden experiment (intended to minimize the effects of environmental factors) strongly overlapped (Figure 3A), although S. limprichtiana started and stopped flowering somewhat earlier than S. grandifolia. All bagged flowers for selfing experiments failed to set fruit. Interspecific crosses resulted in significantly lower (p < 0.001) fruit set (8.1-33.3%) than crosses between populations of the same species (50.0-100%) in the reciprocal crossing experiments (Figure 3B).

Genetic Diversity and Differentiation at SSR Loci
Averaged across all nine populations, observed heterozygosity (Ho) per locus was 0.378, and expected heterozygosity (He) was 0.414 (Table S4). The fixation index averaged across all loci (average F ST = 0.372; Table S4) indicated a pronounced level of genetic differentiation among populations (Bleeker and Hurka, 2001;Helenurm, 2003). PCoA of data acquired from all individuals again identified two well-delimited groups ( Figure 4A), which were supported by STUCTURE analyses, with each group representing a different species (best K = 2, Delta K = 283.38; Figure 4B, Figure S4). AMOVA analyses also revealed clear differentiation between both species and intraspecific populations (F ST = 0.461, p < 0.001; Table 2).

Genetic Diversity and Divergence of Nine Low-Copy Nuclear Genes
Lengths of the nine sequenced nuclear genes extracted from individuals representing populations of two species ranged from 253 to 638 bp, with a total concatenated length of 3803 bp, excluding gaps and missing data. More haplotypes of each gene were detected in S. limprichtiana than in S. grandifolia (Figure 5). Haplotype networks for each gene suggested that two species have completely diverged, and private SNPs of both species were detected in genealogies of seven of the nine loci ( Figure 5). Diversity in S. limprichtiana is higher than in S. grandifolia (Table 3). Tajima's D, Fu and Li's D * and F * neutrality indicators varied widely across the nine genes ( Table 3). No significant Tajima's D values were obtained for any of the nine genes in either species (p < 0.05), but significant D * values were obtained for three loci (SINTR2771_02, MEF9 and SINATR1498_34) and significant F * values for one gene (SINTR2771_02) in S. limprichtiana. Since the three tests did not provide consistent significant indications of deviation from neutrality in any of the nine loci ( Table 3), all of the genes were used for subsequent analyses.
Levels of nucleotide polymorphism (θ w and π) differed strongly among the nine genes ( Table 3). For each gene, the average nucleotide diversity was higher in S. limprichtiana than in S. grandifolia ( Table 3). The minimum number of recombination events (R m ) ranged from 0 to 1 in both species across the nine genes ( Table 3). Both STUCTURE and Neighborjoining tree analyses of all samples identified two well-delimited  . Two species' genetic structure, based on 12 SSR loci and nine low-copy nuclear genes, respectively. (D) Neighbor-joining tree derived from phylogenetic analysis of S. limprichtiana and S. grandifolia (S.l. and S.g. respectively) based on nine low-copy nuclear genes.

Gene Flow Analysis by IMa2
We used IMa2 to examine marginal posterior distributions of the probabilities of the six demographic parameters (T, m GL , m LG , θ A , θ L , and θ G ) of two species under different models based on sequence variations of the nine nuclear loci ( Table 4). All independent runs in the M-model gave consistent maximum likelihood estimates (MLEs) (Table 4, Figure 6), suggesting that the model has high reliability for the considered material. The MLEs indicate that interspecific divergence involved limited migration (m LG = 0.068) from S. grandifolia to S. limprichtian, and more limited migration (m GL = 0.004) in the opposite direction. Therefore, gene flow between S. grandifolia and S. limprichtian appears to have been very limited and asymmetrical. The marginal posterior densities of the divergence parameter according to MLEs, t, showed a sharp peak at 2.825, corresponding to a divergence time of c. 534,000 years, with a 95% highest posterior density (HPD) interval of c. 249,000-1,282,000 years. A much larger effective population size was derived for S. limprichtiana (c. θ 1 = 1.74, N 1 = 82,200 individuals, 95% HPD interval: c. 57,600-116,000 individuals) than for S. grandifolia (c. θ 2 = 0.25, N 2 = 11,900 individuals, 95% HPD interval: c. 5,100-22,900 individuals). The estimated effective population size of two species' ancestral population is c. θ A = 1.60, N A = 75,000 individuals (95% HPD interval: c. 56,700-330,000 individuals), smaller than that of S. limprichtiana, but clearly larger than that of S. grandifolia (Table 4).

DISCUSSION
We quantified morphological divergence and reproductive isolations between S. limprichtiana and S. grandifolia. In addition to the morphological differences mainly in the presence/absence of hairs recognized previously (Zhou et al., 2014), we further found the stable morphological differentiation between two species in the leaf and flower sizes, and numbers of flowers. Reproductive isolations were well-established between two species. Two sets of population genetic data similarly recovered high interspecific differentiation and divergence. All of these findings suggested that both species are "good species" rather than partly differentiated "subspecies" (Zhou et al., 2014). However, the divergence between them were estimated to occur in the recent past, probably during the Pleistocene, when the climate oscillated and forests accordingly fragmented and coalesced in eastern Asia (Qian and Ricklefs, 2000;Yu et al., 2000;Qian et al., 2003). Both stochastic processes and habitat selection probably contributed to such a fast divergence and strong post-pollination isolation.

Interspecific Differentiation
Our common garden experiments (designed to exclude effects of environmental variables) suggest that S. limprichtiana and S. grandifolia have developed stable morphological differentiation, including differences in the presence/absence of hairs on their calyxes, leaf blades and stems, leaf and flower sizes, and numbers of flowers per individual (Figure 2, Figures S2, S3). The morphological differentiation between them is as clear as that between closely related species of other genera in the family Brassicaceae (Zhou et al., 2001). All analyses of population genetic data obtained from analyses of the SSRs and low-copy nuclear genes also clearly delimited two species (Figure 4), in accordance with previous studies based on fewer individuals but different genetic markers (Hu et al., 2015(Hu et al., , 2016. All of these lines of evidence indicate that S. limprichtiana and S. grandifolia are sufficiently well-differentiated and delimited to be acknowledged as separate "good species" (Hu et al., 2015) rather than subspecies (Zhou et al., 2014), although the distinction between subspecies and species is always arbitrary (Liu, 2016). In addition, as both species have the similar genome sizes ( Figure S1, Table S3), their interspecific divergence did not involve any change in ploidy.

Reproductive Isolation
We failed to recover the pre-pollination isolation between two species. Pre-pollination isolation mechanisms, for example, differences in flowering phenology, are believed to have played major roles in the divergence of closely related species, especially species co-occurring in the same sites (Coyne and Orr, 2004). Such pre-pollination isolation mechanisms may frequently arise from selective divergence upon secondary contact following the development of an initial reproductive barrier (Hendry, 2009;Schemske, 2010;Sobel et al., 2010;Nosil, 2012). However, in our common garden experiments, flowering periods of S. limprichtiana and S. grandifolia strongly overlapped (Figure 3A), although we cannot exclude the possibility that this and other pre-pollination reproductive barriers probably occur in natural habitats of two species. S. limprichtiana prefers dry rocky places, while S. grandifolia generally occurs in wet sites along streams (Zhou et al., 2014). Such habitat differentiation probably results in differences in closely related plants' local flowering periods and rates of germination and/or survival probability because of the differences in water availability (Favre et al., 2017). In contrast to the overlapping flowering periods, our artificial crossing experiments suggested the presence of strong post-pollination isolation between S. limprichtiana and S. grandifolia. Interspecific crossing using individuals of either species as the female parents and the other as pollen donors consistently resulted in substantially lower fruit set than intra-specific crossing ( Figure 3B). Such intrinsic barriers arise mainly from Bateson-Dobzhansky-Muller genetic incompatibility due to genetic drift and natural selection, which result in the fixation of different alleles, new mutations and chromosomal rearrangements at different loci (Schluter, 2009). Further comparative analysis of two species' genomes and other molecular functional analyses are also warranted, to identify genetic changes associated with the observed post-pollination.

Pleistocene Divergence
We found that two species have fixed differentiated alleles at seven of nine examined nuclear loci, similarly suggesting the fast divergence in most loci. The high divergence was also recovered based on the SSRs results, consistent with the results based on other markers reported before (Hu et al., 2015(Hu et al., , 2016. However, we dated two species' divergence to within the Pleistocene (Figure 6, Table 4) although such calculations should be treated cautiously because of the no-fossil calibrated mutation rate and  Tajima D' and Fu & Li's D* and F* (Fu and Li, 1993). Significant level: *0.01 ≤ P < 0.05; **0.001 ≤ P < 0.01.  the sequence variation from the limited loci. In addition, we detected extremely low gene flow between two species based on the sequence variations of nine nuclear loci, much lower than rates detected between other pairs of species that diverged recently, especially during the Pleistocene (Ikeda et al., 2009(Ikeda et al., , 2012Wang et al., 2013;Han et al., 2016). What drove the fast divergence, development of the strong reproductive isolation and low gene flow between these two species within such a relatively short time? Three non-exclusive factors may have contributed to the fast differentiation between them. First, evergreen forests during the Pleistocene climatic oscillations might have fragmented as temperatures decreased in eastern Asia (Qian and Ricklefs, 2000;Yu et al., 2000), and the ancestor of these two species probably started to diverge allopatrically. In fact, all specimen records and our field investigations suggest that these two species currently occur in allopatric locations. Second, S. limprichtiana prefers dry habitats while S. grandifolia occurs in the wet habitats along the stream. These contrasted habitats may have produced strong selection pressure on the differentiations of two species (Wang et al., 2017). Finally, the effective population sizes of both species are small. The small populations may have accelerated the fixation of the different alleles between two species although both species are outcrossing. However, elucidating the relative importance of genetic drift and natural selection in the fast divergence and development of post-pollination isolation between S. limprichtiana and S. grandifolia remains challenging.
Our case study partly supported the prediction, derived from previous floristic analyses (e.g., Qian and Ricklefs, 2000), that Pleistocene climate changes likely promoted the production of the numerous young endemic species in eastern Asia through allopatric divergence and demographic oscillations as suggested for other endemic species occurring there (e.g., Chiang et al., 2012;Ikeda et al., 2012;Mitsui and Setoguchi, 2012;Nakamura et al., 2012;Han et al., 2016;Ito et al., 2017;Wang et al., 2017). Although further detailed studies of endemic species' origins are needed, these climate changes may have strongly contributed to the high species diversity in eastern Asia (Calonje et al., 2013). This study also adds to the limited number of known examples of Pleistocene speciation in plants from different regions of the world (e.g., Martin-Bravo et al., 2010;Ikeda et al., 2012;Levsen et al., 2012;Wang et al., 2013Wang et al., , 2017Han et al., 2016).

CONCLUSION
S. limprichtiana and S. grandifolia are well-differentiated in morphological characters, including the presence/absence of hairs, the leaf and flower sizes, and numbers of flowers. We also found that reproductive isolations were well-established between them. Interspecific divergence is high based on both sets of population genetic data. This divergence was further dated within the Pleistocene. All of these findings add our understanding of the divergence and the origin of the endemic species in eastern Asia.

AUTHOR CONTRIBUTIONS
LZ, TZ, and QH: conceived and designed this study; TZ and HH: collected the samples, performed crosses, and morphometric analyses in common garden experiments; TZ, HH, LF, HZ, and LZ: acquired and analyzed data for the study; TZ, LF, HH, LZ, and QH wrote the manuscript.

ACKNOWLEDGMENTS
We are grateful for John Blackwell for editing the English of the initial and final manuscripts. We also thank Professor Jianquan Liu for supervising this work.