Speciation history of a species complex of Primulina eburnea (Gesneriaceae) from limestone karsts of southern China, a biodiversity hot spot

Abstract Limestone karsts in southern China are characterized by high edaphic and topographic heterogeneity and host high levels of species richness and endemism. However, the evolutionary mechanisms for generating such biodiversity remain poorly understood. Here, we performed species delimitation, population genetic analyses, simulations of gene flow and analyses of floral morphological traits to infer the geographic history of speciation in a species complex of Primulina eburnea from limestone karsts of southern China. Using Bayesian species delimitation, we determined that there are seven distinct species that correspond well to the putative morphological species. Species tree reconstruction, Structure and Neighbor‐Net analyses all recovered four lineages in agreement with currently species geographic boundaries. High levels of genetic differentiation were observed both within and among species. Isolation–migration coalescent analysis provides evidence for significant but low gene flow among species. Approximate Bayesian computation (ABC) analysis supports a scenario of historical gene flow rather than recent contemporary gene flow for most species divergences. Finally, we found no evidence of divergent selection contributing to population differentiation of a suite of flower traits. These results support the prevalence of allopatric speciation and highlight the role of geographic isolation in the diversification process. At small geographic scales, limited hybridization occurred in the past between proximate populations but did not eliminate species boundaries. We conclude that limited gene flow might have been the predominant evolutionary force in promoting population differentiation and speciation.

selection across heterogeneous environments can drive population divergence and lead to speciation with or without gene flow (Coyne & Orr, 2004;Via, 2009). To understand the mechanisms underlying the processes of speciation, we need to assess the relative importance of selection, drift and gene flow in driving population and lineage divergence. Natural terrestrial islands, such as limestone karst landforms, provide a particularly interesting system for studying speciation process and the mechanisms underlying colonization, adaptation and diversification (Barbará, Martinelli, Fay, Mayo, & Lexer, 2007;Sarthou, Samadi, & Boisselier-Dubayle, 2001). Like oceanic islands, terrestrial islands generate a range of microhabitats where ecological conditions contrast moderately to extremely from the surrounding landscape and therefore possess unique and endemic species with small population sizes and clearly defined geographic boundaries (Barbará et al., 2007;Sarthou et al., 2001). Thus, terrestrial islands act as "natural experiments" from which to study evolutionary processes of species diversity, in particularly assess the role of genetic drift and selection in population divergence and speciation.
Limestone karsts are the result of weathering of limestone rocks over time, leaving behind a terrain of dramatic peaks and caves (Ford & Williams, 2007). Across Southeast Asia and southern China, karst landscapes are disproportionately threatened by climate change and anthropogenic activities (Clements, Sodhi, Schilthuizen, & Ng, 2006). This has lead the International Union for the Conservation of Nature (IUCN) to recognize karst landscapes worldwide as being in significant need of protection (Watson, Hamilton-Smith, Gillieson, & Kiernan, 1997). With a total area of about 400,000 km 2 , southern China is recognized as the world's type area for karst landform development in the humid tropics and subtropics (Day & Urich, 2000). Karst landforms here are typically dominated by steep cone-shaped karst towers, caves, sinkholes and cliffs, and separated from other outcrops by lowland areas composed of different soil types. Such landscapes are characterized by diverse and extreme environmental conditions, harboring high levels of species diversity and endemism, and have been recognized as a global center of plant diversity (Davis, Heywood, & Hamilton, 1995). Substantial progress has been made in characterizing the geographic patterns of species richness and endemism in southern China (López-Pujol, Zhang, Sun, Ying, & Ge, 2011a, 2011bYing, Zhang, & Boufford, 1993), yet the evolutionary mechanisms for generating the great richness of narrow endemic species in this hot spot remain poorly understood. Very few karst plant species have been the subject of phylogeographic and population genetic analyses.
The few available studies revealed that allopatric speciation via geographic isolation might be the predominant mechanism in karst landforms (Chung et al., 2014;Gao, Ai, Kong, Kang, & Huang, 2015). In the allopatric mode of speciation, nonadaptive divergence, operating via genetic drift due to geographic isolation and founder effects, is expected to play a significant role in generating patterns of species diversity (Chung et al., 2014;Gao et al., 2015). However, karst areas in southern China are characterized by a high edaphic heterogeneity, with contrasting local-scale mosaics of soil types derived from bedrock of differing lithology (e.g., granite) (Hao, Kuang, & Kang, 2015). This edaphic complexity may be a strong driver of diversification and speciation via local adaptation to specific edaphic microhabitats (i.e., specialization), as widely reported in plants (Anacker & Strauss, 2014;Anacker, Whittall, Goldberb, & Harrison, 2011;Schnitzler et al., 2011).
Thus, both genetic drift and divergent selection could be strong in karst island populations. Yet, to date, no studies have reported the relative role of genetic drift, gene flow and selection as potential factors driving speciation in the karst flora in China.
Primulina (Gesneriaceae) is one of the species-rich plant genera in limestone karsts of southern China, where around 170 species have been recorded, among which 150 are endemic (Möller, Wei, Wen, Clark, & Weber, 2016). Primulina is a monophyletic group of evergreen perennials that are widely distributed throughout the lowland karst regions of southern China and northern Vietnam. They represent a group of typical rock-dwelling plants that have adapted to remarkably diverse habitats and niches in limestone caves and crevices. However, this genus displays a high degree of edaphic specialization, with the majority of species occurring in calcareous soils originated from limestone bedrock (i.e., calciphiles; Hao et al., 2015). As karst landforms in southern China are generally scattered and isolated as small limestone hills (i.e., terrestrial islands), most species are microendemics with narrow distribution range, often limited to a single cave or karst limestone hill system (Wei, 2010).
We focus on six species that have been shown by recent phylogenetic analysis to form a monophyletic group (Kang et al., 2014) Ning, Li, Pan, and Kang (2015) described a new rare species from southern Hunan, P. rubribracteata Z.L. Ning & M. Kang, which phylogenetically belongs to the species complex, but which we did not include in this study. Of the six studied species, four (P. lutea, P. polycephala, P. xiziae and P. alutacea) are calcareous soil specialists (i.e., calciphiles) and P. suichuanensis displays edaphic specialization on Danxia soils (after Danxia Mountain in Guangdong Province, China; they are usually red or purple in color and derived from coarse clastic rocks; Peng, Ren, & Pan, 2015), while P. eburnea occurs in both karst and Danxia landforms. Primulina eburnea has the widest distribution range of its genus, whereas the five remaining species are all narrow endemics ( Figure 1). Primulina lutea is one of the rare yellow-flowering species in the genus (Liu & Wei, 2004). Primulina polycephala, P. alutacea and P. suichuanensis occur around Nanling Mountains. Primulina xiziae is a newly described species morphologically similar to be P. eburnea (Li, Wang, Hua, & Wen, 2012). Primulina xiziae is the easternmost species of the genus and is found on shaded rocks in limestone hills at 70-110 m in Zhejiang province, China. This species complex shows clear geographic division and presents variation in floral and habitat specialization, making it an ideal system to study speciation of the karst flora. Previous phylogeographic analysis of this complex based on cpDNA and small numbers of microsatellites revealed strong spatial patterns of genetic structure . However, the cpDNA haplotypes were not sufficiently resolved to infer the geographic history of speciation. A clear understanding of speciation and diversification in karst landscapes is particularly important, as these habitats are disproportionately threatened by climate change and anthropogenic deforestation (Clements et al., 2006).
In this study, we used a panel of nuclear genes and SNPs to (i) determine the evolutionary relationships within and between species and estimate species boundaries using a coalescent-based, Bayesian species delimitation approach, (ii) discover the patterns of genetic structure and levels of genetic diversity across species and populations, (iii) test the temporal and geographic extent of gene flow and compare evolutionary scenarios using approximate Bayesian computation (ABC) and (iv) examine whether the divergence in quantitative traits (Q ST ) significantly exceeded that one detected using neutral markers (F ST ) as signal of divergent selection by pairwise comparisons of F ST and Q ST estimates. Our main objectives are to investigate processes of speciation in the P. eburnea complex and evaluate the role of genetic drift, gene flow and natural selection on population divergence and speciation.

| Population sampling, DNA sequencing and SNP genotyping
A total of 405 individuals from 43 populations (five to 14 per population) were sampled across the geographic range of six closely related Primulina species, P. eburnea, P. lutea, P. polycephala, P. alutacea, P. suichuanensis and P. xiziae ( Figure 1; Table S1). A previous study revealed that P. eburnea was divided into two deeply divergent phylogenetic lineages representing an eastern and a western distribution , so we treated these western and eastern populations as distinct species in the data analysis. In addition, one individual of P. tabacum was sampled as the out-group in phylogenetic analysis. Total genomic DNA was extracted from dried leaves using a modified cetyltrimethylammonium bromide (CTAB) method (Doyle & Doyle, 1987).
Six single-copy nuclear DNA (nDNA) loci developed from comparative transcriptomic data of multiple Primulina species  were amplified and sequenced (Table S2). Polymerase chain reactions (PCRs) were performed in total volume of 25 μl containing 10-40 ng DNA template, 50 mM Tris-HCl (Takara), 1.5 mM MgCl 2 (Takara), 0.5 mM dNTPs (Takara), 0.75 unit of Taq polymerase (ExTaq, Takara) and 2 mM of each primer in final concentration. All reactions were performed using the following temperature profile: 5 min at 94°C, 35 cycles of 30 s at 94°C, 30 s of annealing at 57°C, and 90 s at 72°C, with a final 10-min extension at 72°C. The PCR products were sequenced directly using the BigDye Terminator Cycle Sequencing kit according to the instructions of the manufacturer (Applied Biosystems, Foster City, CA, USA) on an ABI 3730XL automated sequencer. All sequences were deposited in GenBank with accession numbers KY239620-KY241086.  Table S1. Study region indicated by a rectangle in the lower right corner on the MassARRAY platform (Sequenom, San Diego, CA). The reliability and cross-species transferability of these loci have been validated in P. eburnea . The paired PCR and single extension primers were designed based on the flanking sequences using Sequenom Assay Design 3.1 software. The whole pipeline included PCR primer mix preparation, multiple PCR amplification, shrimp alkaline phosphatase reaction, extension reaction, resin mixture, sample spotting, mass spectrometry and genotype reading. The multiple PCR amplification was performed in a reaction volume of 5 μl containing 10 ng DNA template, 1 unit of Taq polymerase, 2 mM MgCl 2 , 0.5 mM dNTPs and 0.1 μM PCR primer mix in final concentration. The multiple PCR temperature profile was set as follows: 2 min at 95°C, 45 cycles of 30 s at 95°C, 30 s at 56°C, and 60 s at 72°C, with a final 5-min extension at 72°C.

| Genetic diversity
Heterozygous nuclear sequences were phased to alleles using Phase v2.2.1 (Stephens, Smith, & Donnelly, 2001). SeqPhase (Flot, 2010) was used to convert input and output formats for Phase. Phase computations were conducted with default settings (phase threshold = 90%, iteration steps = 100, thinning interval = 1, burn-in steps = 100) and were repeated three times to ensure consistency. The program DnasP 5.10.01 (Librado & Rozas, 2009) was used to calculate the number of segregating sites (S), the number of haplotypes (h), haplotype diversity (Hd) and nucleotide polymorphism (π) for each population and species. Because nonneutral data may bias the following analyses, we performed multiple neutrality tests to confirm whether the loci were neutral. We calculated Tajima's D (Tajima, 1989) to test for departure from the standard neutral model using DnasP. To discriminate between neutral and adaptive evolution, we performed the multilocus HKA test (Hudson, Kreitman, & Aguade, 1987) using the HKA package (https://bio.cst.temple.edu/~hey). Individual runs of the HKA test were performed for the contrast between each of the seven species and the out-group.
For SNP data, we trimmed the dataset by excluding loci with minor allele frequency (MAF) lower than 0.01 and individuals with more than 10% missing genotypes across all loci. This led to a total of 384 individuals from 43 populations with 414 loci in the data analysis. We estimated percentage of polymorphic loci (PPL), the number of different alleles (A), the number of effective alleles (A e ), observed heterozygosity (H O ), expected heterozygosity (H E ) and fixation index (F IS ) for each population and species using GenAlEx v6.502 (Peakall & Smouse, 2012). The population assignment, genetic differentiation, AMOVA and network clustering based on SNP data were similar as those based on nDNA sequence data.

| Population structure and genetic clustering
We used a Bayesian clustering approach implemented in Structure v2.3.4 (Pritchard, Stephens, & Donnelly, 2000) to infer the number of clusters and to assign individuals to these clusters. We ran 10 independent replicates for each K value between 1 and 10, with 100,000 burn-in MCMC steps followed by 150,000 recorded steps.
The admixture model and uncorrelated allele frequencies between populations were specified for each run. We used structure harvester (Earl & vonHoldt, 2012) to summarize the results including the identification of optimal K through the method of Evanno, Regnaut, and Goudet (2005). Results across replicate runs were combined using the program clumPP (Jakobsson & Rosenberg, 2007) and visualized with Distruct v1.1 (Rosenberg, 2004). We performed analysis of molecular variance (AMOVA) using arlequin v3.5.2.2 (Excoffier & Lischer, 2010) to calculate the respective proportion of genetic variation at three hierarchical levels: among groups, among populations within groups and within populations. The groups were defined according to species or the inferred clusters. The significances of variance components were estimated with 1,000 permutations. We also calculated genetic differentiation (F ST ) between population and species using arlequin v3.5.2.2 (Excoffier & Lischer, 2010).
We further performed a phylogeographic analysis using the Neighbor-Net method implemented in SplitsTree4 (Huson & Bryant, 2006) based on the uncorrected p-distance between individuals to reveal the overall genetic pattern of relationships among populations and individuals.

| Phylogenetic reconstruction
Prior to phylogenetic reconstruction, the best substitution model for each nDNA locus (Table S2) was determined under the Akaike information criterion as implemented in jModelTest v2.1.3 (Darriba, Taboada, Doallo, & Posada, 2012). Gene trees of the identified haplotypes were reconstructed using Bayesian inference as implemented in mrbayes v3.2.6 (Ronquist et al., 2012). The Bayesian analyses were run three times with one cold and three heated chains for ten million Markov chain Monte Carlo generations. Trees were sampled every 1,000 generations, and the first 25% generations were discarded as burn-in. tracer v1.5 (Rambaut & Drummond, 2009) was used to check convergence and effective sample sizes of parameter estimation. We estimated the species tree of the six Primulina species and the population tree of the 38 populations using *beast v1.8.3 (Heled & Drummond 2010). We ran the MCMC analysis three times for 500 million generations with sampling every 50,000 generations. The relaxed clock model with an uncorrelated lognormal distribution, the constant population size and the Yule speciation tree prior were also specified in the settings. The final trees were visualized using both FigTree v1.4.2 (http://tree.bio.ed.ac.uk/software/figtree/) and DensiTree in the beast package.

| Species delimitation
Bayes factor delimitation is a useful method that compares the marginal likelihoods of competing species delimitation hypotheses (Grummer, Bryson, & Reeder, 2014). Based on a combination of the results from population clustering and phylogenetic analyses, four alternative species delimitation models were tested against the current delimitation model (seven species): (i) a six-species model where the P. eburnea (east) and P. xiziae were lumped together; (ii) a second six-species model where the eastern and western populations of P. eburnea were lumped together; (iii) a five-species model where the three species of P. polycephala, P. alutacea and P. suichuanensis were lumped together; and (iv) a four-species model where the P. eburnea (east) and P. xiziae were lumped together and simultaneously, the three species of P. polycephala, P. alutacea and P. suichuanensis were lumped together. We used *beast v1.8.3 to estimate species trees for the above species delimitation models based on the six nDNA loci dataset. Each locus was given its own partition with the preferred substitution model, and the uncorrelated relaxed lognormal clock was specified. *beast analyses were performed for 500 million generations, logging every 50,000 generations. For each hypothesis, three *beast replicates were conducted to check convergence and effective sample sizes using tracer v1.5. Marginal likelihoods were estimated using path sampling (PS) and stepping stone (SS) methods (Baele et al., 2012;Grummer et al., 2014), with 100 path steps, a chain length of 100,000 generations and likelihoods logged every 100 generations.
The four alternative models were compared against the current delimitation model, and the significantly better models were determined whether the Bayes factor value (twice the difference in marginal likelihood estimates) was larger than 10 (Kass & Raftery, 1995).

| Estimation of gene flow
Four approaches were used to estimate spatial and temporal pattern of gene flow. Firstly, to distinguish hybridization or introgression from incomplete lineage sorting in causing phylogenetic incongruence between gene trees and species trees, we used Jml v1.3.0 (Joly, 2012) to test whether the minimum distance between sequences of two species is smaller than expected under a scenario without hybridization.
To do so, the observed minimum distance between sequences is compared to a null distribution obtained by simulating a scenario without hybridization. The posterior distributions of species trees, population sizes and divergence times included in the output of the *beast analyses were used as input for Jml. The best substitution models for individual genes were specified in the control file, while the default values were used in other settings.
Secondly, we estimated population migration rate based on the isolation with migration (IM) model implemented in ima2p (Sethuraman & Hey, 2016). We performed IM estimation with three separate datasets: the four lineages including all species (four-species model; Figure 2), the PAS lineage including three species of P. polycephala, P. alutacea and P. suichuanensis (three-species model; Figure 2) and the EX lineage including P. eburnea (east) and P. xiziae (two-species model; Figure 2). For each IM model, we ran analyses of 20 Markov chains with heating terms of 0.96 (-ha) and 0.90 (-hb), for 300,000 MCMC steps after a burn-in of 300,000 steps. Of the nucleotide substitution models implemented in the program ima2p, HKY was the best model and thus was applied for all genes. The mutation rates were calculated based on the transcriptome-wide divergence  and molecular dating results . Finally, we used IMfig (https:// bio.cst.temple.edu/~hey/software) to generate the result figures of the IM models.
Thirdly, we used the approximate Bayesian computation approach as implemented in the program MSABC (Pavlidis, Laurent, & Stephan, 2010) to evaluate the posterior probabilities of four alternative evolutionary scenarios: zero, constant, historical or recent gene flow models ( Figure 2). In the first two models, lineages evolve without gene flow or with constant gene flow throughout the divergence process. In the other two models, gene flow only occurs in or following the initial divergence stages. We ran 10 6 simulations under each model based on our empirical priors estimated in this study or previous attempts Gao et al., 2015), including locus number, sample size, sequence length, population mutation rate and divergence time. A migration rate of 1 N e m (0.5 N e m in each direction) was used because F I G U R E 2 Four scenarios of gene flow (indicated by shading) simulated under the neutral coalescent and evaluated using approximate Bayesian computation from the four lineages among EubW, EX, Lut and PAS (four-species model); three species among Pol, Alu and Sui (threespecies model) and two species between EbuE and Xiz (two-species model). Lineage and species abbreviations are defined in Table 1 Four species Three species Two species this is the maximum migration rate allowed between separate species (Hey, 2009). For the ABC analysis, we chose F ST as the summary statistic to describe the degree of differentiation among lineages. We transformed observed sequence data into ms-like files using the fas2 ms perl script in the MSABC package and calculated the same summary statistic. Based on the comparison between simulated and observed statistics, we used the "postpr" function in R package "abc" (Csilléry, François, & Blum, 2012) to estimate posterior probabilities, adopting the rejection method and a tolerance of 0.001.
Finally, we performed coalescent simulations in conjunction with the genealogical sorting index (gsi; Cummings, Neel, & Shaw, 2008) to determine the temporal period of gene flow, i.e., whether gene flow occurred among ancestral or extant populations. The gsi is a statistic that estimates the degree of exclusive ancestry of each delimited taxon on a scale ranging from 0 (nonexclusive) to 1 (monophyletic).
Coalescent simulations were performed using mccoal (Rannala & Yang, 2003;Yang & Rannala, 2010). In our simulations, a symmetric migration matrix with the migration rate constant at 0 or 0.5 N e m in both directions was used, and divergence times and population sizes were derived from our empirical estimates . The four simulated scenarios included zero, historical, constant, and recent gene flow scenarios ( Figure 5). Based on 10,000 simulated trees under each scenario, a gsi value was calculated for each species within each tree using the R package genealogicalSorting v0.92 (Cummings et al., 2008).
An empirical gsi value was also calculated based on the inferred tree in mrbayes v3.2.6. The comparison of empirical gsi values to the simulated gsi values was used to determine the timing of migration events. for normality using Kolmogorov-Smirnov normality test, and the nonnormal distributions (i.e., PW1 and PW2) were log transformed.

| Analysis of floral morphological data
To assess whether divergence of measured floral traits was associated with genetic divergence, we calculated the degree of quantitative phenotypic divergence between populations as Q ST is the variance between populations and σ W 2 is the within-population variance (Spitze, 1993

| Genetic diversity and pairwise differentiation
The six nDNA regions varied in length from 669 to 1,197 bp and totaled 5,613 bp with 803 variable sites, 532 of which were parsimonyinformative (Table S2) containing either a single or two haplotypes (Table S3). Species-level haplotype and nucleotide diversity estimates from the nDNA dataset indicated fairly even diversity across species, with slightly lower diversity in P. xiziae, the easternmost species (Table 1). For the SNP dataset, similarly low levels of genetic diversity were found across populations, with mean percentage of polymorphic SNPs varying from 8.94% to 68.36% (Table S4). Again, similar levels of effective number of alleles and expected heterozygosity over all loci (A e = 1.2-1.4 and H E = .096-.256) were found across species, with slightly lower diversity in P. xiziae (Table 1).
T A B L E 1 Summary statistics of polymorphism for the seven taxa based on six single-copy nuclear genes and 414 SNP loci We found high genetic differentiation at both the population and species levels. An AMOVA across all populations shows that most of the variance (68% and 59.1% for nDNA and SNP datasets, respectively) is found among populations (Table 2). For the nDNA sequence data, pairwise estimates of F ST between populations ranged from 0.073 to 1 (mean = 0.712), and all comparisons were highly significant (P < .001; Table S5). This translated into a high level of genetic differentiation among species (F ST = .204-.751; Table 3). For the SNP dataset, pairwise F ST values between populations varied from 0 to 0.893 (mean = 0.539; Table S6). Pairwise F ST values between species were lower than that calculated from the nDNA dataset, but highly significant (Table 3).

| Phylogenetic reconstruction
As expected, individual gene trees were poorly resolved (Fig. S1

| Genetic clustering
The Structure analyses based on nDNA dataset inferred K = 5 as the optimal number of groups. The supported major groups are generally congruent with genetic clusters generated with Neighbor-Net networks, but the western P. eburnea cluster was further divided into two groups. On the other hand, the model with K = 8 was also highly supported under the Evanno method, with a delta K value (Fig. S2) similar to the preferred K = 5. At K = 8, the cluster of the eastern P. eburnea and P. xiziae were further divided into three groups: the eastern P. eburnea populations excluding the NY01, the NY01 and P. xiziae. However, the clustering at K = 8 is essentially identical to K = 7, because no individuals are completely assigned to a new group (Figure 4). All of the main genetic clusters were geographically coherent, and only populations GXQZ01, NY01 and JXGF01 showed admixture above a cutoff level of 22%.
T A B L E 2 Results of the analyses of molecular variance (AMOVA) based on nuclear DNA sequence (nDNA) and SNP datasets For the SNP dataset, Structure analyses revealed fewer distinct groups, with the inferred K = 2 and K = 4 as the optimal and sub-optimal numbers of groups, respectively. However, a population grouping of of western P. eburnea (Fig. S3).

| Species delimitation
Bayes factor delimitation marginal likelihood values for alternative species delimitation models are listed in Table 4. The favored result is a seven-species model that corresponds to the morphological species description, with exception for the recognition of western and eastern clades of P. eburnea as two distinct species. In fact, the "six-species" model that grouped western and eastern populations of P. eburnea into a single species was the lowest ranked model, indicating that this species needs taxonomic revision.

| Test of hybridization and gene flow
The Jml analyses demonstrated that nine comparisons exhibited genetic distances that were significantly lower than expected (P < .05; Table S7), suggesting that incomplete lineage sorting alone could not explain the discordance across datasets and potential interspecies hybridization might have occurred.
ima2 analyses with the three datasets are summarized in Table 5.  Table 6).
F I G U R E 3 Estimated species trees from six single-copy nuclear loci using the framework of the multispecies coalescent algorithm implemented in *beast: (a) DensiTree visualization of consensus; (b) consensus tree. Posterior probabilities are given at each node. The four recovered lineages were defined as Lut, EbuW, EX and PAS

| Differentiation of floral traits
We found pronounced differences in the scored flower traits among populations, but none of the Q ST estimated differed significantly from SNP-based F ST estimates, with the exception of corolla length (CL) in the eastern populations of P. eburnea (Table S8) in all comparisons (Figure 6), suggesting that neutral processes seem to be the most probable cause for divergence in the studied traits.

| DISCUSSION
We conducted a detailed assessment of speciation history of the P. eburnea complex from a global biodiversity hot spot through the coupling of multilocus data, phylogenetic reconstruction, coalescent methods and model testing. Our data support the prevalence of allopatric speciation and highlight the role of neutral processes such as genetic drift in the diversification process. At small geographic scales, limited hybridization occurred in the past between proximate populations but did not eliminate species boundaries. We found no evidence of divergent selection contributing to population differentiation of a suite of flower traits.

| Phylogenetic relationships, species boundary and genetic differentiation
A previous study based on cpDNA sequences was not able to resolve phylogenetic relationships among taxa of the P. eburnea complex   . Based on the coalescent species tree reconstruction of multilocus nuclear DNA sequences, the present study revealed that most recognized species receive strong statistical support. Both nDNA and SNP markers revealed high genetic differentiation among most population-and species pairs. Actually, the six nuclear genes revealed a much deeper pattern of population structure than that of 414 SNP markers, suggesting that relatively small numbers of nuclear DNA loci without linkage may be sufficient for detecting population genetic structure when it is pronounced. The finding of high genetic differentiation among populations was consistent with our hypothesis of restricted gene flow and long-term persistence in geographically isolated karst island landscapes . As a habitat specialist confined to island-like karst caves or outcrops, the geographic isolation and small size of populations is expected to lead to reductions in genetic variation as a result of random genetic drift and a decrease in gene flow. Accordingly, most populations show lowto-moderate genetic variation (Table S3 and

Path sampling
Stepping stone sampling The four models of gene flow are illustrated in Figure 2. The model with the highest posterior probability is in bold. Lineage and species abbreviations are defined in Table 1. genera are characterized by their high levels of local endemism and karst habitat specialization in southern China. Taken together, these results have highlighted the role of isolation and genetic drift in shaping population differentiation in the karst endemics.

| Spatial and temporal patterns of gene flow
Despite the pattern of high genetic differentiation observed for within-and among species, the structure and SplitsTree analyses provide evidence for considerable admixture in several western  Gao et al., 2015) between the two species. In fact, P. xiziae was formerly recognized as P. eburnea and recently classified as a separate species due to its distinctive morphology and phenology (Li et al., 2012 EbuW of these results should be cautious because the choice of summary statistics and model selection criteria might influence the conclusions that one draws (e.g., Csilléry, Blum, Gaggiotti, & Francois, 2010).
The lack of genetic differentiation among species of the PAS lineage may be due to their recent divergence or high levels of gene flow.
However, IMA2 analysis did not find evidence of significant gene flow among species. This result seems to suggest a scenario of allopatric speciation, where a decline in gene flow could be due simply to the physical isolation between recently diverged species. However, the ABC models support a scenario of historical gene flow rather than a scenario of isolation model with no gene flow. Therefore, this lack of differentiation is most likely a consequence of their recent origin and the high level of shared ancestral polymorphism considerably reducing the power of tests for finding postdivergence gene flow.
Nevertheless, the level of gene flow we detected in all comparisons was far below 0.5 N e m per generation (Table 5), a value assumed insufficient to counteract random drift (Wright, 1969). The limited gene flow estimated for Primulina is in line with the results from the previous cpDNA and microsatellites study, where most populations of the complex were fixed for a single private haplotype, which resulted in a high population differentiation (F ST = 0.944) .
The high genetic differentiation revealed by multiple nuclear and SNP markers in this study also provides evidence for limited gene flow.
Therefore, gene flow, if any, has only a minor impact on population divergence and speciation of Primulina.

| Speciation history of the P. eburnea complex
Based on the strong signature of phylogeographic structure and current nonoverlapping distributions, Gao et al. (2015) argued for an allopatric mode of speciation for the P. eburnea complex. In agreement with this viewpoint, the nonsignificant gene flow inferred among species of the lineage PAS supports allopatric divergence without subsequent genetic exchange. In contrast, the significant postdivergence gene flow within the four lineages and between P. xiziae and eastern P. eburnea is therefore inconsistent with a strictly allopatric speciation mode, which assumes complete termination of gene flow after initial geographic isolation between two species. Furthermore, apparent genetic admixture in several populations (Figure 4; Fig. S3) along range boundaries is evidence for gene flow between species.
However, our simulation results showed that interspecies gene flow is a historical event and is not currently ongoing. These results are essentially in congruence with a scenario of parapatric speciation, where genetic introgression occurred at some point in the past between populations with proximate distributions before complete reproduction isolation. As posited above, however, the inferred gene flow for all comparisons was far below 0.5 N e m per generation, suggesting a minor role of gene flow on speciation. Therefore, allopatric speciation could have been the predominant geographic modes of speciation. Similar to our results, allopatric speciation accompanied by hybridization has been argued as a major evolutionary driver for endemics from the Balkan Peninsula (Olšavská, Slovák, Marhold, Štubňová, & Kučera, 2016).
Few studies have tested the geographic mode of speciation in karst areas; it is believed that the highly heterogeneous topography of the mountain chains in southern China favored the emergence of new lineages (López-Pujol et al., 2011a, 2011b, which are generally thought to have evolved as a result of microallopatric speciation. Seeds of Primulina are small but not small enough for effective wind dispersal and do not have dispersal aids (Li & Wang, 2004). Therefore, we would argue that the geographic isolation and the low dispersal ability in Primulina are main driving forces that have promoted genetic isolation among populations, leading to allopatric speciation at local scales, thus representing a case of "terrestrial island" speciation.
Recently, allopatric speciation has been proposed as the primary mode responsible for species diversification of Begonia (Chung et al., 2014), a genus similar to Primulina in that they exhibit extremely high species richness and endemism confined to karst landscapes in southern China and Southeast Asia. A similar scenario has also been suggested for explaining the high species diversity and narrow endemism of several members of New World Gesneriaceae, associated with geographically isolated rock outcrops in southern and southeastern Brazil (Perret, Chautems, Spichiger, Barraclough, & Savolainen, 2007).
Despite speciation with or without gene flow, adaptive differentiation can occur if natural selection is strong enough to cause significant fitness differences across different environment conditions. Floral traits are commonly thought to be a key innovation involved in the diversification of angiosperms. Our field work indicates that the P. eburnea complex is pollinated by bees (i.e., Amegilla spp., unpublished data). Although we cannot determine whether pollinator shifts have occurred among species, Ellis and Johnson (2009) suggested that floral variation can evolve without pollinator shifts, potentially driven by pleiotropic effects, herbivores, abiotic agents and/or genetic drift in isolated and small populations. In the P. eburnea complex, lack of significant difference in the levels of differentiation between flower traits (Q ST ) and neutral genetic markers (F ST ) indicates that genetic drift may have played a major role in the observed pattern of differentiation among species. As a habitat specialist confined to island-like karst caves or outcrops, population sizes of Primulina are usually very small.
In fact, both theoretical and empirical studies (Barrett & Schluter, 2007;Fisher, 1930;Woolfit & Bromham, 2005) suggest selection may become less efficient relative to drift in small populations. However, it should be noted that the difference in quantitative traits in this study was only estimated from floral traits. Leaf traits are another group of impotently functional traits, which deserve future investigations.
The high genetic differentiation, low levels of gene flow and lacking evidence of divergent selection on flower traits support the idea that neutral processes, such as genetic drift, have played an important role in driving species diversification in the complex. Although the contribution of genetic drift in speciation has been the subject of keen debate in this century (e.g., Sobel, Chen, Watt, & Schemske, 2010;Turelli, Barton, & Coyne, 2001), several empirical studies demonstrated that divergent selection had likely played little role in population divergence, suggesting that genetic drift might have played a large role in speciation (Boucher, Zimmermann, & Conti, 2016;Comes, Tribsch, & Bittkau, 2008;Kozak, Weisrock, & Larson, 2006;Verboom, Bergh, Haiden, Hoffmann, & Britton, 2015). Drift is especially likely to be an important driver of speciation for species occurring in small and isolated (allopatric) populations from island or island-like systems, where the small effective population size and limited gene flow are expected to increase the relative importance of drift as an evolutionary process (Templeton, 2008). In Primulina, the low genetic diversity and strong population differentiation could be related to small population sizes in highly dissected landscapes and low dispersal abilities, which strengthen the influence of genetic drift.

| Conservation implications
Due to its unique natural landforms and the associated special biota, "South China Karst" has been listed as a World Natural Heritage site, which warrants high priority for conservation. Plants restricted to karst island habitats, because of their small population sizes and physical barriers, are particularly sensitive to habitat loss and climate change (Damschen, Harrison, Arckerly, Fernandez-Going, & Anacker, 2012). The results of this study have obvious conservation implications, not just for the P. eburnea species complex, but also other karst endemics in this region. The low levels of genetic variation within populations (c.30%) suggest a substantial amount of genetic variation would be lost following even a few local populations. The naturally small population sizes and lack of gene flow make the populations highly vulnerable to extinction due to both demographic and genetic factors. Fine-scale studies of population connectivity incorporating landscape information will further inform conservation efforts and population management.

ACKNOWLEDGEMENT
This work was supported by grant from the Natural Science Foundation of China-Guangdong Natural Science Foundation Joint Project (U1501211 and 31370366). We want to thank Bin Hu and Yong Gao for their help in laboratory work. We also thank Eric Bishop von Wettberg for his critical reading and correction of English of the manuscript as well as two anonymous referees for their useful comments on a previous version of the manuscript.