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 warm-temperate 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, 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.

FIGURE 1

Figure 1. (A) Distributions of the two Sinalliaria species showing locations of the populations of S. limprichtiana (red circles) and S. grandifolia (green triangles) sampled in this study. (B) Photograph of S. limprichtiana. (C) Photograph of S. grandifolia.

In the present study, we first examined the ploidy level, morphological differentiation, and strength of two species' pre- and 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.

Materials and Methods

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⁸ 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.

FIGURE 2

Figure 2. (A) Hair presence/absence in S. limprichtiana (left) and S. grandifolia (right). (B) Contrasting sizes of the basal leaves of S. limprichtiana and S. grandifolia (Significance level, **P < 0.01). (C) Morphological clusters based on Principal Coordinate Analysis of 15 traits of S. limprichtiana (red circles) and S. grandifolia (green triangles).

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' post-pollination 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.

TABLE 1

Table 1. Characteristics of 12 EST-SSR and nine nuclear low-copy gene primer pairs validated in a survey of Sinalliaria populations.

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™ 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 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 and Nielsen, 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⁵ iterations, followed by 1 × 10⁶ 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⁻⁸ 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⁻⁶ substitutions per locus per year (yr) to scale the demographic parameters.

Results

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).

FIGURE 3

Figure 3. (A) Flowering periods of cultivated individuals of the two Sinalliaria species in the common garden (red lines for individuals representing different populations of S. limprichtiana and green lines for individuals representing different populations of S. grandifolia). (B) Fruit sets of monitored female parents of two species (red and green for S. limprichtiana and S. grandifolia respectively), pollinated interspecifically (groups 1 and 3) or by members of different populations of the same species (groups 2 and 4). (a) indicates results of tests of differences between groups 1 and 2, and (b) indicates results of tests of differences between groups 3 and 4. Significance level: **0.001 ≤ P < 0.01; *0.01 ≤ P < 0.05.

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).

FIGURE 4

Figure 4. (A) Results of Principal Coordinate Analysis of genetic distances among nine populations of two species, based on 12 SSR loci (red and green spots for individuals of S. limprichtiana and S. grandifolia, respectively). (B,C). 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.

TABLE 2

Table 2. AMOVA analyses of genetic partitions between and within S. limprichtiana and S. grandifolia for 12 SSR loci and nine nuclear low-copy genes.

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.

FIGURE 5

Figure 5. Genealogies of the nine low-copy nuclear loci. Colors in the pie chart indicate the haplotype origin; red, S. limprichtiana; green, S. grandifolia. The size of the pie is proportional to the haplotype frequency found in two species. Marked branch lengths indicate the mutation steps.

TABLE 3

Table 3. Summary of the nucleotide polymorphisms and neutrality tests of nine nuclear genes for S. limprichtiana and S. grandifolia.

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 Neighbor-joining tree analyses of all samples identified two well-delimited genetic groups (Figures 4C,D). Similarly, AMOVA detected higher differentiation between species than between intraspecific populations (Table 2).

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.74, N₁ = 82,200 individuals, 95% HPD interval: c. 57,600–116,000 individuals) than for S. grandifolia (c. θ₂ = 0.25, N₂ = 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).

TABLE 4

Table 4. MLEs and the 95% HPD intervals of demographic parameters.

FIGURE 6

Figure 6. Marginal distributions of posterior probabilities of six demographic parameters estimated by the IM model. (A) Population size of S. limprichtiana (θ_L), S. grandifolia (θ_G) and ancestral species (θ_A); (B) Migration rate between S. limprichtiana and S. grandifolia; (C) Divergence time of S. limprichtiana and S. grandifolia.

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, 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, 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 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, 2012; Wang 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., 2013, 2017; Han 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.

Funding

This work was supported by grants from National Natural Science Foundation of China (31590821, 31700323, 41771055), National Key Project for Basic Research (2014CB954100) and the Fundamental Research Funds for the Central Universities (YJ201714).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2018.00077/full#supplementary-material

References

Al-Shehbaz, I. A., and Guang, Y. (2000). A revision of the Chinese endemic Orychophragmus (Brassicaceae). Novon 10, 349–353. doi: 10.2307/3392983

CrossRef Full Text | Google Scholar

Axelrod, D., Al-Shehbaz, I., and Raven, P. (1998). “History of the modern flora of china. floristic characteristics and diversity of east asian plants,”in Proceedings of the First International Symposium of Floristic Characteristics and Diversity of East Asian plants, eds A. Zhang and S. Wu (Beijing: China Higher Education Press), 43–55.

Bandelt, H. J., Forster, P., and Röhl, A. (1999). Median-joining networks for inferring intraspecific phylogenies. Mol. Biol. Evol. 16, 37–48. doi: 10.1093/oxfordjournals.molbev.a026036

PubMed Abstract | CrossRef Full Text | Google Scholar

Bitner-Mathé, B. C., and Klaczko, L. B. (1999). Plasticity of Drosophila melanogaster wing morphology: effects of sex, temperature and density. Genetica 105, 203–210. doi: 10.1023/A:1003765106652

PubMed Abstract | CrossRef Full Text | Google Scholar

Bleeker, W., and Hurka, H. (2001). Introgressive hybridization in Rorippa (Brassicaceae): gene flow and its consequences in natural and anthropogenic habitats. Mol. Ecol. 10, 2013–2022. doi: 10.1046/j.1365-294X.2001.01341.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Calonje, M., Meerow, A. W., Knowles, L., Knowles, D., Griffith, M. P., Nakamura, K., et al. (2013). Cycad biodiversity in the Bahamas Archipelago and conservation genetics of the threatened Zamia lucayana (Zamiaceae). Oryx 47, 190–198. doi: 10.1017/S0030605312000129

CrossRef Full Text | Google Scholar

Carstens, B., and Knowles, L. L. (2007). Shifting distributions and speciation: species divergence during rapid climate change. Mol. Ecol. 16, 619–627. doi: 10.1111/j.1365-294X.2006.03167.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Chiang, Y. C., Huang, B. H., and Liao, P. C. (2012). Diversification, biogeographic pattern, and demographic history of Taiwanese Scutellaria species inferred from nuclear and chloroplast DNA. PLoS ONE 7:e50844. doi: 10.1371/journal.pone.0050844

PubMed Abstract | CrossRef Full Text | Google Scholar

Coyne, J. A., and Orr, H. A. (2004). Speciation. Sunderland, MA: Sinauer Associates.

Google Scholar

Day, R. W., and Quinn, G. P. (1989). Comparisons of treatments after an analysis of variance in ecology. Ecol. Monogr. 59, 433–463. doi: 10.2307/1943075

CrossRef Full Text | Google Scholar

Dolezel, J., and Bartos, J. A. N. (2005). Plant DNA flow cytometry and estimation of nuclear genome size. Ann. Bot. 95, 99–110. doi: 10.1093/aob/mci005

PubMed Abstract | CrossRef Full Text | Google Scholar

Doyle, J. J. (1987). A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochem. Bull. 19, 11–15.

Google Scholar

Engelhardt, B. E., and Stephens, M. (2010). Analysis of population structure: a analysis of population structure: a unifying framework and novel methods based on sparse factor analysis. PLoS Genet. 6:e1001117. doi: 10.1371/journal.pgen.1001117

PubMed Abstract | CrossRef Full Text | Google Scholar

Evanno, G., Regnaut, S., and Goudet, J. (2005). Detecting the number of clusters of individuals using the software STRUCTURE: a simulation study. Mol. Ecol. 14, 2611–2620. doi: 10.1111/j.1365-294X.2005.02553.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Excoffier, L., and Lischer, H. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Favre, A., Widmer, A., and Karrenberg, S. (2017). Differential adaptation drives ecological speciation in campions (Silene): evidence from a multi-site transplant experiment. New Phytol. 213, 1487–1499. doi: 10.1111/nph.14202

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y. X., and Li, W. H. (1993). Statistical tests of neutrality of mutations. Genetics 133, 693–709.

PubMed Abstract | Google Scholar

Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, Q., Higashi, H., Mitsui, Y., and Setoguchi, H. (2016). Lineage isolation in the face of active gene flow in the coastal plant wild radish is reinforced by differentiated vernalisation responses. BMC Evol. Biol. 16:84–96. doi: 10.1186/s12862-016-0655-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Helenurm, K. (2003). Genetic diversity in the rare, insular endemic Sibara filifolia (Brassicaceae). Madroño 50, 181–186. Available online at: http://www.jstor.org/stable/41425513

Google Scholar

Hendry, A. P. (2009). Ecological speciation! or the lack thereof? Can. J. Fish. Aquat. Sci. 66, 1383–1398. doi: 10.1139/F09-074

CrossRef Full Text | Google Scholar

Hey, J. (2010). Isolation with migration models for more than two populations. Mol. Biol. Evol. 27, 905–920. doi: 10.1093/molbev/msp296

PubMed Abstract | CrossRef Full Text | Google Scholar

Hey, J., and Nielsen, R. (2004). Multilocus methods for estimating population sizes, migration rates and divergence time, with applications to the divergence of Drosophila pseudoobscura and D. Persimilis. Genetics 167, 747–760. doi: 10.1534/genetics.103.024182

PubMed Abstract | CrossRef Full Text | Google Scholar

Hey, J., and Nielsen, R. (2007). Integration within the Felsenstein equation for improved Markov chain Monte Carlo methods in population genetics. Proc. Natl. Acad. Sci. U.S.A. 104, 2785–2790. doi: 10.1073/pnas.0611164104

PubMed Abstract | CrossRef Full Text | Google Scholar

Hu, H., Al-Shehbaz, I. A., Sun, Y., Hao, G., Wang, Q., and Liu, J. (2015). Species delimitation in Orychophragmus (Brassicaceae) based on chloroplast and nuclear DNA barcodes. Taxon 64, 714–726. doi: 10.12705/644.4

CrossRef Full Text | Google Scholar

Hu, H., Hu, Q., Al-Shehbaz, I. A., Luo, X., Zeng, T., Guo, X., et al. (2016). Species delimitation and interspecific relationships of the genus Orychophragmus (Brassicaceae) inferred from whole chloroplast genomes. Front. Plant Sci. 7:1826. doi: 10.3389/fpls.2016.01826

PubMed Abstract | CrossRef Full Text | Google Scholar

Hubisz, M. J., Falush, D., Stephens, M., and Pritchard, J. K. (2009). Inferring weak population structure with the assistance of sample group information. Mol. Ecol. Resour. 9, 1322–1332. doi: 10.1111/j.1755-0998.2009.02591.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeda, H., Carlsen, T., Fujii, N., Brochmann, C., and Setoguchi, H. (2012). Pleistocene climatic oscillations and the speciation history of an alpine endemic and a widespread arctic-alpine plant. New Phytol. 194, 583–594. doi: 10.1111/j.1469-8137.2012.04061.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Ikeda, H., Fujii, N., and Setoguchi, H. (2009). Application of the isolation with migration model demonstrates the Pleistocene origin of geographic differentiation in Cardamine nipponica (Brassicaceae), an endemic Japanese alpine plant. Mol. Biol. Evol. 26, 2207–2216. doi: 10.1093/molbev/msp128

PubMed Abstract | CrossRef Full Text | Google Scholar

Ito, T., Yu, C. C., Nakamura, K., Chung, K. F., Yang, Q., Fu, C. X., et al. (2017). Unique parallel radiations of high-mountainous species of the genus Sedum (Crassulaceae) on the continental island of Taiwan. Mol. Phylogenet. Evol. 113, 9–22. doi: 10.1016/j.ympev.2017.03.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Koch, M. A., Haubold, B., and Mitchell-Olds, T. (2000). Comparative evolutionary analysis of chalcone synthase and alcohol dehydrogenase loci in Arabidopsis, Arabis, and related genera (Brassicaceae). Mol. Biol. Evol. 17, 1483–1498. doi: 10.1093/oxfordjournals.molbev.a026248

PubMed Abstract | CrossRef Full Text | Google Scholar

Krzanowski, W. (2000). Principles of Multivariate Analysis. Oxford: Oxford University Press.

Google Scholar

Levsen, N. D., Tiffin, P., and Olson, M. S. (2012). Pleistocene speciation in the genus Populus (Salicaceae). Syst. Biol. 61, 401–412. doi: 10.1093/sysbio/syr120

PubMed Abstract | CrossRef Full Text | Google Scholar

Librado, P., and Rozas, J. (2009). DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics 25, 1451–1452. doi: 10.1093/bioinformatics/btp187

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J. Q. (2016). “The integrative species concept” and “species on the speciation way.” Biodivers. Sci. 24, 1004–1008. doi: 10.17520/biods.2016222

CrossRef Full Text

Liu, J. Q., Sun, Y. S., Ge, X. J., Gao, L. M., and Qiu, Y. X. (2012). Phylogeographic studies of plants in China: advances in the past and directions in the future. J. Syst. Evol. 50, 267–275. doi: 10.1111/j.1759-6831.2012.00214.x

CrossRef Full Text | Google Scholar

Liu, K. B. (1988). Quaternary history of the temperate forests of China. Quat. Sci. Rev. 7, 1–20. doi: 10.1016/0277-3791(88)90089-3

CrossRef Full Text | Google Scholar

Martin-Bravo, S., Valcarcel, V., Vargas, P., and Luceno, M. (2010). Geographical speciation related to Pleistocene range shifts in the western Mediterranean mountains (Reseda sect. Glaucoreseda, Resedaceae). Taxon 59, 466–482. Available online at: http://www.jstor.org/stable/25677605

Google Scholar

Mitsui, Y., and Setoguchi, H. (2012). Recent origin and adaptive diversification of Ainsliaea (Asteraceae) in the Ryukyu Islands: molecular phylogenetic inference using nuclear microsatellite markers. Plant Syst. Evol. 298, 985–996. doi: 10.1007/s00606-012-0608-6

CrossRef Full Text | Google Scholar

Nakamura, K., Chung, K. F., Huang, C. J., Kono, Y., Kokubugata, G., and Peng, C. I. (2012). Extreme habitats that emerged in the Pleistocene triggered divergence of weedy Youngia (Asteraceae) in Taiwan. Mol. Phylogenet. Evol. 63, 486–499. doi: 10.1016/j.ympev.2012.01.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Nei, M. (1987). Molecular Evolutionary Genetics. New York, NY: Columbia University Press.

Nosil, P. (2012). Ecological Speciation. Oxford: Oxford University Press.

Google Scholar

Peakall, R., and Smouse, P. E. (2012). GenAlEx 6.5: genetic analysis in Excel. Population genetic software for teaching and research-an update. Bioinformatics 28, 2537–2539. doi: 10.1093/bioinformatics/bts460

PubMed Abstract | CrossRef Full Text

Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959. Available online at: http://www.genetics.org/content/155/2/945

PubMed Abstract | Google Scholar

Qi, W., Lin, F., Liu, Y., Huang, B., Cheng, J., Zhang, W., et al. (2016). High-throughput development of simple sequence repeat markers for genetic diversity research in Crambe abyssinica. BMC Plant Biol. 16:139. doi: 10.1186/s12870-016-0828-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, H., and Ricklefs, R. E. (2000). Large-scale processes and the Asian bias in species diversity of temperate plants. Nature 407, 180–182. doi: 10.1038/35025052

PubMed Abstract | CrossRef Full Text | Google Scholar

Qian, H., Song, J. S., Krestov, P., Guo, Q., Wu, Z., Shen, X., et al. (2003). Large-scale phytogeographical patterns in East Asia in relation to latitudinal and climatic gradients. J. Biogeogr. 30, 129–141. doi: 10.1046/j.1365-2699.2003.00807.x

CrossRef Full Text | Google Scholar

Qiu, Y. X., Fu, C. X., and Comes, H. P. (2011). 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. 59, 225–244. doi: 10.1016/j.ympev.2011.01.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Schemske, D. W. (2010). Adaptation and the origin of species. Am. Nat. 176, S4–S25. doi: 10.1086/657060

PubMed Abstract | CrossRef Full Text | Google Scholar

Schluter, D. (2009). Evidence for ecological speciation and its alternative. Science 323, 737–741. doi: 10.1126/science.1160006

PubMed Abstract | CrossRef Full Text | Google Scholar

Schürch, S., Pfunder, M., and Roy, B. A. (2000). Effects of ants on the reproductive success of Euphorbia cyparissias and associated pathogenic rust fungi. Oikos 88, 6–12. doi: 10.1034/j.1600-0706.2000.880102.x

CrossRef Full Text | Google Scholar

Sobel, J. M., Chen, G. F., Watt, L. R., and Schemske, D. W. (2010). The biology of speciation. Evolution 64, 295–315. doi: 10.1111/j.1558-5646.2009.00877.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Stephens, M., Smith, N. J., and Donnelly, P. (2001). A new statistical method for haplotype reconstruction from population data. Am. J. Hum. Genet. 68, 978–989. doi: 10.1086/319501

PubMed Abstract | CrossRef Full Text | Google Scholar

Strasburg, J. L., and Rieseberg, L. H. (2010). How robust are “Isolation with Migration” analyses to violations of the IM model? A simulation study. Mol. Biol. Evol. 27, 297–310. doi: 10.1093/molbev/msp233

PubMed Abstract | CrossRef Full Text | Google Scholar

Tajima, F. (1989). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.

PubMed Abstract | Google Scholar

Thompson, J. D., Gibson, T. J., Plewniak, F., Jeanmougin, F., and Higgins, D. G. (1997). The CLUSTAL_X Windows Interface: flexible strategies for multiple sequence alignment aided by quality analysis tools. Nucleic Acids Res. 25, 4876–4882. doi: 10.1093/nar/25.24.4876

PubMed Abstract | CrossRef Full Text | Google Scholar

Tiffney, B. H. (1985). Perspectives on the origin of the floristic similarity between eastern Asia and eastern North America. J. Arnold Arbor. 66, 73–94. doi: 10.5962/bhl.part.13179

CrossRef Full Text | Google Scholar

Wang, Q., Abbott, R. J., Yu, Q. S., Lin, K., and Liu, J. Q. (2013). Pleistocene climate change and the origin of two desert plant species, Pugionium cornutum and Pugionium dolabratum (Brassicaceae), in northwest China. New Phytol. 199, 277–287. doi: 10.1111/nph.12241

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Q., Yu, Q. S., and Liu, J. Q. (2011). Are nuclear loci ideal for barcoding plants? A case study of genetic delimitation of two sister species using multiple loci and multiple intraspecific individuals. J. Syst. Evol. 49, 182–188. doi: 10.1111/j.1759-6831.2011.00135.x

CrossRef Full Text | Google Scholar

Wang, Y. H., Comes, H. P., Cao, Y. N., Guo, R., Mao, Y. R., and Qiu, Y. X. (2017). Quaternary climate change drives allo-peripatric speciation and refugial divergence in the Dysosma versipellis-pleiantha complex from different forest types in China. Sci. Rep. 7:40261. doi: 10.1038/srep40261

PubMed Abstract | CrossRef Full Text | Google Scholar

Watterson, G. A. (1975). On the number of segregating sites in genetical models without recombination. Theor. Popul. Biol. 7, 256–276. doi: 10.1016/0040-5809(75)90020-9

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Chen, X., Ni, J., Cheddadi, R., Guiot, J., Han, H., et al. (2000). Palaeovegetation of China: a pollen data-based synthesis for the mid-Holocene and last glacial maximum. J. Biogeogr. 27, 635–664. doi: 10.1046/j.1365-2699.2000.00431.x

CrossRef Full Text | Google Scholar

Zanella, C. M., Palma-Silva, C., Goetze, M., and Bered, F. (2016). Hybridization between two sister species of Bromeliaceae: Vriesea carinata and V. incurvata. Bot. J. Linn. Soc. 181, 491–504. doi: 10.1111/boj.12424

CrossRef Full Text | Google Scholar

Zhou, T. Y., Lu, L. L., Yang, G., and Al-Shehbaz, I. A. (2001). “Orychophragmus Bunge,” in Flora of China, Vol. 8, eds Z. Y. Wu and P. H. Raven (Beijing; St. Louis, MO: Science Press; Missouri Botanical Garden Press), 29–31.

Zhou, Y. Y., Zhang, H. W., Hu, J. Q., and Jin, X. F. (2014). Sinalliaria, a new genus of Brassicaceae from eastern China, based on morphological and molecular data. Phytotaxa 186, 188–198. doi: 10.11646/phytotaxa.186.4.2

CrossRef Full Text | Google Scholar