Investigating the Genetic Diversity, Population Differentiation and Population Dynamics of Cycas segmentifida (Cycadaceae) Endemic to Southwest China by Multiple Molecular Markers

Climate change, species dispersal ability and habitat fragmentation are major factors influencing species distribution and genetic diversity, especially for the range-restricted and threatened taxa. Here, using four sequences of chloroplast DNAs (cpDNAs), three nuclear genes (nDNAs) and 12 nuclear microsatellites (SSRs), we investigated the genetic diversity, genetic structure, divergence time and population dynamics of Cycas segmentifida D. Y. Wang and C. Y. Deng, a threatened cycad species endemic to Southwest China. High levels of genetic diversity and genetic differentiation were revealed in C. segmentifida. Haplotypes of networks showed two evolutionary units in C. segmentifida, with the exception of the nuclear gene GTP network. Meanwhile, the UPGMA tree, structure and PCoA analyses suggested that 14 populations of C. segmentifida were divided into two clades. There was significant effect of isolation by distance (IBD) in this species. However, this species did not display a significant phylogeographic structure. The divergence time estimation suggested that its haplotypes diverged during the Middle Pleistocene. Additionally, the population dynamics inferred from different DNA sequences analyses were discordant. Bottleneck analysis showed that populations of C. segmentifida did not experience any recent bottleneck effect, but rather pointed to a contraction of its effective population size over time. Furthermore, our results suggested that the population BM which held an intact population structure and occupied undisturbed habitat was at the Hardy–Weinberg equilibrium, implying that this population is a free-mating system. These genetic features provide important information for the sustainable management of C. segmentifida.


INTRODUCTION
Southwest China is seen as one of the world's biodiversity hotspots because of its complicated terrain, diversified climates and habitats (Myers et al., 2000). The region possesses an extremely high richness of species and concentrates numerous endemic and endangered plant species (Zhao and Gong, 2015). Climatic oscillations in the Pleistocene led to drastic environmental changes repeatedly, which also substantially influenced the species' distribution and evolution, even extinction (Hewitt, 1996(Hewitt, , 2000(Hewitt, , 2004Comes and Kadereit, 1998). Climate change, dispersal, and habitat fragmentation are considered to be major factors influencing the current distribution and genetic diversity of species (Hewitt, 2004;Broquet and Petit, 2009;Arenas et al., 2012Arenas et al., , 2014Mona et al., 2014). In the Pleistocene, most parts of China were not covered by large ice sheet and had a relatively warm climate (Weaver et al., 1998;Williams et al., 1998). But climatic oscillations during the Pleistocene have also had severe effects on the divergence and population dynamics of many extant species (Hewitt, 1996(Hewitt, , 2000(Hewitt, , 2004Avise, 2000). In recent years, there has been an increasing body of literature on the effect of climatic oscillations on species' divergence and population dynamics (Zhan et al., 2011;Liu et al., 2013;Feng et al., 2014Feng et al., , 2016bGong et al., 2015).
As an ancient lineage, cycads are ideal materials with which one can explore how plants have responded to historical climate oscillations. Cycads belong to gymnosperms and are considered the most primitive living seed plants. However, recent research has proposed that a synchronous global rediversification occurred in cycads during the late Miocene (Nagalingum et al., 2011). Although the cycad lineage is ancient, the extant cycad species have evolved recently and are not older than 12 million years based on a fossil-calibrated molecular phylogeny (Nagalingum et al., 2011). They are distributed in tropical and subtropical regions and comprise two families (Cycadaceae, Zamiaceae) with 10 genera (Christenhusz et al., 2011). There is only one cycad genus, Cycas (Cycadaceae), in China (Hill et al., 2004). A study combining ancestral area reconstructions with fossil evidence revealed that South China is the origin of Cycas (Xiao and Moller, 2015). In South China, Cycas species usually grow on low-altitude slopes of ridges and cliffs along river valleys. Cycas species in China are all facing potential endangerment challenges due to over-collection because of their edible stem as well as ornamental attributes and the destruction of their habitats for the cultivation of commercial plants. Due to their evolutionary importance, cycads have been studied in many fields, including phylogeny, population genetics, phylogeography and conservation (Chiang et al., 2009;Cibrian-Jaramillo et al., 2010;Zhan et al., 2011;Feng et al., 2014Feng et al., , 2016aGong et al., 2015;Liu et al., 2015;Xiao and Moller, 2015;Zheng et al., 2016;Yessoufou et al., 2017). However, there are still some Cycas species that have not yet been targeted for protection.
For example, C. segmentifida D. Y. Wang and C. Y. Deng (Wang and Deng, 1995) draws little attention. It is endemic to Southwest China, and occurs primarily in the valleys of the You River basin of the eastern Yunnan, southwestern Guizhou, and northwestern Guangxi provinces. This species is characterized by its blue-green petiole of young leaves and dichotomous or sometimes forked lateral spine of megasporophyll (Wang and Deng, 1995). After the description of C. segmentifida, several new species of Cycas in this region were published, based on one or more specific morphological feature Zhang et al., , 1999. These laterdescribed species are morphologically similar to C. segmentifida, giving rise to the long-controversial confusion on species classification (Chen and Stevenson, 1999;Wang, 2000;Huang, 2001;Whiteloek, 2002;Ma, 2005). Combining evidence from chloroplast and nuclear DNA sequences, microsatellite analysis, and the geographical distribution, these ambiguous species were delimited and included in C. segmentifida (Feng et al., 2016a). Cycas segmentifida occupies two types of habitat according to the soil matrix, sand and karst, on which it grows. The populations of the sand type are found under the forest canopy, while the populations of the karst type are scattered on isolated limestone hills with few shrubs. The two types of habitat with obviously different environments lead to morphological differences, such as the length and width of leaflets or acuminate apex of pinnae and shorter carpophylls.
This species is dioecious, allogamous, and insect pollinated. Its pollinators are mainly beetles. As an inland Cycas species, C. segmentifida is classified into Cycas section Stangerioides Smitinand which has seeds that are short on thick spongy tissue, always sink in water and contain virulent cycasin, precluding their dispersal by water or animals for long distances (Dehgan and Yuen, 1983;Schneider et al., 2002;Xiao and Moller, 2015). Like other Cycas species from sect. Stangerioides, its fertile seeds are large and heavy, usually falling and germinating near the mother plant, as found in field survey.

Sampling and Genotyping
A total of 238 individuals were obtained from 14 populations, representing the entire natural distribution area of C. segmentifida. Young and healthy leaves were dried in silica gel immediately after collection. Within the 238 samples, 7 to 10 individuals from each population were chosen for chloroplast and nuclear DNA sequencing, while all of the 238 individuals were used for the microsatellite study. Geographical information of the 14 populations and numbers of individuals used in the DNA sequencing and microsatellite analyses are presented in Supplementary Table S1 and Figure 1, respectively.
Microsatellite loci were screened from nuclear microsatellites which were developed for other Cycas plants ( Cibrian-Jaramillo et al., 2008;Wang et al., 2008;Yang et al., 2008;Li et al., 2009;Zhang et al., 2009Zhang et al., , 2010Ju et al., 2011). PCR amplification and microsatellite genotyping were conducted with the same protocol as used in C. simplicipinna (Feng et al., 2014). Finally, we chose 12 polymorphic microsatellite loci, the same used in our study on species delimitation of the C. segmentifida complex (Feng et al., 2016a), for C. segmentifida.

Data Analysis of DNA Sequences
SeqMen was used to edit and assemble sequences. Sequence multiple alignments were subsequently performed in Bioedit, version 7.0.4.1 (Hall, 1999). The four cpDNA regions were combined for a congruency test using PAUP * 4.0b10 (Swofford, 2003). For three nuclear genes, heterozygous sites were identified by overlapping peaks in chromatograms, and the nuclear sequences were resolved by applying the algorithms of PHASE (Stephens et al., 2001;Stephens and Donnelly, 2003) in the software package DnaSP, version 5.0 (Librado and Rozas, 2009). The combined cpDNA sequences and the phased nuclear sequences were used in analyses that followed. The mapping work was done using the ArcGIS 10.2 software (Esri Inc.).
First, we used DnaSP, version 5.0, to detect recombination in nuclear genes. We calculated indices to measure the level of genetic variation, haplotypes, Nei's nucleotide diversity (Pi) and haplotype diversity (Hd). Gene diversity in total populations (H T ), within-population gene diversity (H S ), (Nei, 1973) and two indices of genetic differentiation, N ST and G ST , were also calculated using Permut 1.0 (Pons and Petit, 1996). We compared G ST and N ST using the U test. An estimation of genetic variation that was assigned within and among populations was performed with an analysis of molecular variance (AMOVA) using Arlequin, version 3.11 (Excoffier et al., 2005). The pollen/seed migration ratio (r) was calculated using a modified equation: r = mp/ms = [(1/F ST (n) -1) -2(1/F ST (c) -1)]/(1/F ST (c) -1) (Ennos, 1994;Petit et al., 2005), where F ST values (rather than G ST ) taken as estimators of population differentiation are derived from AMOVA, mp is the pollen migration rate, ms is the seed migration rate, F ST (n) is the nuclear (nDNA) F ST and F ST (c) is the chloroplast (cpDNA) F ST . Here, r G = mp (GTP)/ms; r P = mp (PHYP)/ms; r R = mp (PPRC)/ms. Phylogenetic relationships of haplotypes were constructed using Bayesian methods implemented in MrBayes, version 3.1.2 (Ronquist and Huelsenbeck, 2003), with Cycas edentata and Cycas rumphii as outgroups. We also used Network, version 4.2.0.1 (Bandelt et al., 1999), to estimate the degree of relatedness among cpDNA and nDNA haplotypes with indels treated as single mutational events.
The evolutionary rates previously estimated for seed plants, 1.01 × 10 −9 and 5.1-7.0 × 10 −9 (Graur and Li, 2000) mutations per site per year for synonymous sites, were used to estimate the coalescent time of haplotypes for cpDNA and nDNA, respectively. We used BEAST, version 1.6.1 (Drummond and Rambaut, 2007), to estimate the time of divergence using a strict molecular clock and the HKY model that was determined by MEGA, version 5, based on the Akaike Information Criterion (Tamura et al., 2011). We also used the BEAST program to create a Bayesian Skyline Plot to infer the history demography for C. segmentifida. Time of divergence and the mutation rate posterior estimates were obtained by Markov Chain Monte Carlo (MCMC) analysis. The ESS parameter of each process was checked by TRACER, version 1.5 (Rambaut and Drummond, 2009), until a stable value exceeding 200 was reached, suggesting that there was acceptable mixing and sufficient sampling. The subsequent three stable running log files and tree files were combined into pairs. The combined log and tree files were used to construct a Bayesian skyline plot in TRACER, version 1.5. We used DnaSP, version 5.0, to examine a pairwise mismatch distribution and neutrality tests, including Tajima's D, Fu and Li's D * and F * and Fu's F S (Fu, 1997), to further investigate population dynamics of the species. The sum-of-squared deviations (SSD) and raggedness index as well as P-values were calculated with the software Arlequin, version 3.11 (Excoffier et al., 2005).

Data Analysis of Microsatellites
Dataset editing and formatting were performed in GenAlEx, version 6.3 (Peakall and Smouse, 2006). Genetic diversity indices, including the number of alleles (N A ), effective number of alleles (A E ), private alleles (A P ), expected heterozygosity (H E ), observed heterozygosity (H O ), information index (I), fixation index (F) and percentage of polymorphic loci (PPB), were calculated using GenAlEx, version 6.3, and POPGENE, version 1.32 (Yeh et al., 1997), with mutual correction. Allelic richness (A R ) was calculated in the software FSTAT, version 1.2 (Goudet, 1995). The differentiation index F ST between pairs of populations was computed with Arlequin, version 3.11 (Excoffier et al., 2005). Isolation by distance (IBD) was tested by performing Mantel tests in GenAlEx, version 6.3, on the correlation of genetic distance [F ST /(1 − F ST )] with geographic distance for all pairs of populations. Genepop, version 4.1.4 (Rousset, 2008) was used to calculated F ST /(1 − F ST ). The pollen/seed migration ratio (r) was also calculated. Tests for departure from Hardy-Weinberg equilibrium (HWE) were performed in each locus and each population as well as a globally unified population using Genepop, version 4.1.4 (Rousset, 2008).
To gain insight into the population genetic structure of C. segmentifida, multiple approaches were used in this study. Initially, the unweighted pair group mean analysis (UPGMA) was performed using TEPGA, version 1.3 (Miller, 1997), with 5,000 permutations. Next, we conducted a Bayesian analysis of population structure with STRUCTURE, version 2.2 (Pritchard et al., 2000). Number of clusters (K) was set from 1 to 20, and each K run 20 times under 1 × 10 5 subsequent MCMC samplings and a subsequent burn-in of 1 × 10 5 iterations. The combination of admixture and correlated-allele frequencies model was used in this analysis. We evaluated the most likely number of groupings using K and the log-likelihood value in the program STRUCTURE HARVESTER, version 0.6.8 (Earl and vonHoldt, 2012). Finally, based on Nei's (1973) genetic distances, an individual-based principal coordinate analysis (PCoA) was conducted using MVSP, version 3.12 (Kovach, 1999).
The effective population sizes of each population were estimated in the program LDNe at three levels of the lowest allele frequency (0.01, 0.02, 0.05) with a 95% confidence interval (Waples and Do, 2008). The bottleneck effect based on different models and methods was tested in BOTTLENECK, version 1.2.02 (Piry et al., 1999), aiming to explore population dynamics. In this analysis, we chose the stepwise mutation model (SMM) and the two-phased model (TPM). Under the two models, the standardized differences test was removed from this study because this test is typically used only when at least 20 polymorphic loci are available. Two other methods (Sign tests and Wilcoxon tests) were applied in the analyses. We also used a mode shift model to test for bottlenecks in each population. These methods implemented in BOTTLENECK are most powerful unless bottlenecks are severe and recent. Moreover, we further investigated a genetic bottleneck using the Garza-Williamson index (GWI, also called M-ratio, the ratio of number of alleles to range in allele size) (Garza and Williamson, 2001) which was calculated by Arlequin, version 3.11 (Excoffier et al., 2005). When seven or more loci are analyzed, the GWI is lower than the critical Mc value of 0.68, a value obtained from bottlenecked populations, which suggests population decline in history (Garza and Williamson, 2001;Excoffier et al., 2005).

DNA Sequence Variations
The combined cpDNA had a 3,165 bp consensus length with a significant rate of homogeneity (P = 1, >0.5) based on the congruency test, suggesting that there was a high degree of homogeneity among the four cpDNA regions. They contained 22 polymorphic sites and seven haplotypes (segH1-segH7) across the 137 individuals (14 populations, Supplementary Table S1) of C. segmentifida. Of these, only two haplotypes segH1 and segH5, were shared by multiple populations, whereas the other five haplotypes were specific to one population (Supplementary Table S2). Detection of recombination in nuclear genes showed that no recombination occurred in GTP. In contrast, PHYP and PPRC had one and two recombination events, respectively. The aligned nuclear genes GTP, PHYP, and PPRC had 561, 930, and 718 bp consensus lengths, include 7, 12, and 16 polymorphic sites, and 5, 10, and 14 haplotypes, respectively. In GTP, the haplotype segG1 was the most abundant and predominant in all 14 populations. Haplotypes segG4 and segG5 were private. For PHYP, haplotype segP1 was the most frequent and was widely shared by all 14 populations. Haplotypes segP3, segp6, segP9, and segP10 were specific to single population JZ, YX, PHG, and BB. In PPRC, segR2 was the most widely distributed haplotype. Information on cpDNA and three nDNA haplotypes and their distribution in populations are shown in Figure 1 and Supplementary Table S2, respectively. Genetic diversity indices Hd and Pi for each population, summarized in Supplementary Table S2, were highly variable among cpDNA and nDNA. In sum, except for population BB, which displayed both cpDNA nucleotide and haplotype diversity, the remaining 13 populations have very low genetic diversity based on cpDNA data. Three nuclear genes showed variable genetic diversity within 14 populations. The nuclear gene PHYP had the highest haplotype diversity in C. segmentifida, PPRC had the highest nucleotide diversity and GTP showed the lowest genetic diversity. Total genetic diversity H T was higher than the average intrapopulation diversity Hs, except for the nuclear gene GTP, which had basic equal values (Supplementary Table S3), meaning that there were high levels of genetic differentiation (except GTP, N ST ranging from 0.282 to 0.996; G ST ranging from 0.229 to 0.949). N ST was not significantly greater than G ST (P > 0.05), indicating that C. segmentifida shows no correspondence between haplotype similarities and their geographic distribution.
The AMOVA revealed that almost all variation (99.80%) was partitioned among populations based on cpDNA data, whereas higher variations (98.98, 71.38, and 57.41%) were revealed within populations than among populations based on nuclear genes GTP, PHYP, and PPRC, respectively. Except for the gene GTP, F ST values for cpDNA and the other two nuclear genes ranged from 0.286 to 0.998 (Table 1), indicating highly significant genetic differentiation among C. segmentifida populations. The pollen/seed migration ratios were calculated as r G = 49399, r P = 1243.755 and r R = 581.086, suggesting that pollen flow was significantly higher than seed flow in C. segmentifida.

Haplotype Relationships and Divergence Times
Bayesian inference revealed the phylogenetic relationships among cpDNA and nDNA haplotypes with C. edentata and C. rumphii as outgroups. The four phylogenetic trees were strongly supported with high Bayesian probabilities, revealing the monophyly for C. segmentifida haplotypes (Supplementary Figure S1). The cpDNA haplotype phylogenetic tree showed that haplotypes segH5 and segH7 grouped together, and the other five haplotypes were clustered into a clade, indicating that they were more closely related than others. In the GTP haplotype tree, segG5 was the first diverged haplotype, while the relationship of the remaining four haplotypes, which formed a clade, was not resolved. Two clades were also revealed in the PHYP haplotypes Bayesian phylogram, showing that haplotypes segP3, segP4, segP6, and segP9 were clustered into one clade, while the other six haplotypes were clustered into another clade, which revealed a closer relationship among them. Fourteen PPRC haplotypes were grouped into three clades, and segR14 was the earliest diverged haplotype in the PPRC phylogram. Of the other two clades, one contained haplotypes segR3, segR4, segR6, segR9, segR10, segR12, and segR13, and the other comprised  Table S2. the remaining six haplotypes, indicating that they shared closer relationships.
The haplotype network of cpDNA was concordant with the haplotype network of nuclear genes PHYP and PPRC, which displayed two centrally located nodes representing hypothetical ancestral haplotypes with a higher frequency (Figures 2A,C,D). The remaining haplotypes were linked to these central haplotypes by one to nine steps in a star-like network. Absence of some haplotypes caused reticulate evolutionary relationships in networks of cpDNA and PPRC. In the network of GTP (Figure 2B), haplotype segG1 occurred at the highest frequency and was in an internal node location, indicating it may be an ancestral haplotype. The remaining four haplotypes were linked to the central haplotype by one or four steps in a star-like network.
The BEAST-derived trees for the four DNA sequences revealed similar topologies (Figure 3). The seven cpDNA haplotypes clustered into two lineages and diverged at approximately 1.701 million years ago (MYA). Of the two lineages, one included five haplotypes with segH1 as the ancestral haplotype and the other included two haplotypes with segH5 as the ancestral haplotype (Figures 2, 3A). For the three nDNA haplotype trees, haplotypes were also divided into two lineages, coalescing at 0.524 MYA (GTP), 0.480 MYA (PHYP), and 0.684 MYA (PPRC), respectively (Figures 3B-D). In total, most tip haplotypes in the four BEASTderived trees diverged recently from their common ancestors. The divergence times estimated from the four DNA sequences entirely implied that haplotypes of C. segmentifida diverged in the Middle Pleistocene.

Neutrality Test, Mismatch Analysis, and the Bayesian Skyline Plot
Values for Tajima's D, Fu and Li's D * , and F * , Fu's Fs tests as well as SSD and raggedness statistics are presented in Supplementary  Table S4. All of the values were positive and significant for cpDNA and were positive and non-significant for PHYP, indicating that C. segmentifida had not undergone a recent population expansion. For the nuclear gene GTP, except for SSD and raggedness, all values were negative and significant, suggesting the presence of population growth. The nuclear gene PPRC showed negative values for Tajima's D as well as for Fu and Li's D * and F * , indicating C. segmentifida had undergone a recent population expansion. Mismatch distribution analyses displayed multimodal graphs for cpDNA, PHYP, and PPRC, but a unimodal graph for GTP (Supplementary Figure S2). A unimodal curve is indicative of recent population expansion, whereas a multimodal curve is indicative of a population at demographic equilibrium. More detailed evidence of population dynamics came from the Bayesian Skyline Plots. Based on cpDNA data, we found that the population size of C. segmentifida remained constant over a long period of time; the population experienced a contraction only recently, i.e., between 0.1 MYA and the present day ( Figure 4A). During the recent 0.01 million years, rapid population growth and subsequently stabilization were detected in C. segmentifida based on the gene GTP (Figure 4B), while the Bayesian Skyline Plot of the gene PHYP ( Figure 4C) showed this species has had a population contraction since approximately 0.02 MYA, with subsequent slight expansion in demography. In contrast, the Bayesian Skyline Plot of the gene PPRC showed that C. segmentifida experienced a constant population size over 0.27-0.07 MYA, before undergoing a decline from approximately 0.07 MYA, and then experienced a population expansion from 0.01 MYA to present ( Figure 4D).

Nuclear Microsatellite Genotyping
A total of 117 alleles were identified by the 12 microsatellites across the 238 individuals of C. segmentifida, and the number of alleles per locus was between 3 (Cy-TaiEST-SSR11) and 44 (cha-estssr01). The locus cha-estssr01 had the highest genetic diversity while Cy-TaiEST-SSR11 had the lowest by comparison of genetic parameters A R , N A , A E , and I (Supplementary Table S5). Diversity estimates from 12 microsatellites also varied among populations (Table 2), with the highest and lowest measures of diversity consistently found in populations BB and LK, respectively. Fixation indices (F) (Wright, 1978) were negative for populations JZ, BM, BD, and LLB, but positive for the other 10 populations, with a mean value F = 0.066, indicating outcrossing within those four populations but inbreeding within most populations. This inbreeding resulted in a deficiency of heterozygotes and significant deviations from HWE at most loci and in most populations (Supplementary Table S6). It was noteworthy that population BM is accorded with HWE. Effective population sizes at the lowest allele frequency (= 0.05) are shown in Table 2, revealing only three populations, BD, PH, and BB, whose Ne (81.4,72.3,71.5) were greater than 50. The AMOVA revealed that more variation (77.11%) was partitioned within populations than between populations, with the significant genetic differentiation coefficient F ST = 0.229 ( Table 1). The pollen/seed migration ratio (r) was evaluated as 1678.039, suggesting more pollen flow than seed flow. The correlation between genetic and geographic distances was significant (P = 0.001, <0.05), suggesting that C. segmentifida has significant effect of IBD (Figure 5).
The UPGMA clustering dendrogram showed that individuals belonging to populations BY, JZ, LK, LKA, BM, NZ, and BD clustered into one group (Clade I), while the remaining individuals clustered into a second group (Clade II) ( Figure 6A). Structure ( Figure 6B) and PcoA (Figure 6C) revealed the same results as UPGMA, indicating that 14 populations were grouped into two clusters.
No population had a significant excess of heterozygosity in the two methods under two models in the bottleneck analysis, indicating that 14 populations did not deviate from mutationdrift equilibrium (Table 3)

High Levels of Genetic Variations and Differentiation
Genetic diversity is the basis for a species' survival, development, and evolution. Generally, geographical distribution, population number and size, and breeding system all affect genetic diversity in plant species (Hamrick et al., 1992;Hamrick and Godt, 1996;Nybom, 2004). Cycads are dioecious and long-lived plants. They have experienced a long process of evolution (Gao and Thomas, 1989;Axsmith et al., 2003) and should own high genetic diversity. However, a recent research proposed that the extant cycads evolved recently (Nagalingum et al., 2011). In Asia, most Cycas species are narrowly distributed and have similar life history traits, such as dioecy, pollination syndromes and longer life span (Wu and Raven, 1999). Comparing the genetic diversity of C. segmentifida with other Asian Cycas species using the similar molecular markers is reasonable. In contrast, C. segmentifida had slightly lower genetic diversity than C. simplicipinna (Feng et al., 2014), C. multipinnata , C. guizhouensis (Feng et al., 2016b) and C. dolichophylla (Zheng et al., 2016), but had higher genetic diversity than C. debaoensis (Zhan et al., 2011;Gong and Gong, 2016) and C. diannanensis (Liu et al., 2015). H T estimated from C. segmentifida by four cpDNAs was 0.745, which was higher than 170 plant species' mean value of H T = 0.67 (Petit et al., 2005), implying that C. segmentifida had a high level of genetic diversity. Generally, as an ancient and woody gymnosperm species, cycads are considered to possess high genetic diversity within population and low level of genetic differentiation among populations (Hamrick et al., 1992). However, in this study, cpDNA data demonstrated significant population differentiation (F ST = 0.998) (Wright, 1978) within C. segmentifida, which was much greater than the average F ST value estimated from other seed plants based on maternally inherited markers (mean F ST = 0.670) (Petit et al., 2005). In contrast to the significant high value of F ST obtained with cpDNA, lower values were derived from nuclear genes, but they still indicated high genetic differentiation within C. segmentifida (Wright, 1978). The pollen-to-seed migration ratios (r) (Ennos, 1994) illustrated that high pollen flow but limited seed flow among populations is the most likely explanation for the higher genetic differentiation based on cpDNA than nuclear genes. It can also explain the AMOVA results that more genetic variation existed among populations and less within populations based on cpDNA, with the opposite based on nuclear genes and microsatellites. Most of the C. segmentifida populations deviated significantly from HWE, together with the fixation indices (F) (Wright, 1978) basically greater than zero ( Table 2), illustrating that C. segmentifida populations are notably deficient in heterozygosity and have experienced severe levels of inbreeding. The fact is that its large and heavy seeds germinate and grow near the mother plant, increasing the chances of inbreeding within the species. Consequently, this species presents a high level of genetic differentiation.

Obvious IBD, Two Genetic Groups, and Middle Pleistocene Divergence
Population genetic structure is mainly affected by some factors such as habitat, differentiation history, and gene flow. Usually, populations cluster according to habitat types. For example, populations of C. debaoensis were divided into two groups according to two types of habitats (Zhan et al., 2011). Similarly to C. debaoensis, C. segmentifida occupied two habitat types (sand and karst); however, it is not clear whether its populations can be differentiated into two different clusters. N ST was not significantly greater than G ST following U tests, which indicated that there was no distinct phylogeographic structure in C. segmentifida. However, a significant correlation between nuclear genetic and geographic distance in this species was observed, indicating this species was in accord with IBD. In addition, a clear genetic structure of two genetic groups was detected in C. segmentifida, but the two genetic groups did not completely comply with the division of sand and karst habitats. The two central populations (BA and YX, Supplementary Table S1) distributed in the sand habitat were grouped into the karst habitat based on the microsatellite data, perhaps indicating the pollen gene flow occurred between them and the vicinity populations. Our result confirmed higher pollen flow than seed flow in C. segmentifida, and this was also demonstrated in other studies . The two populations growed on sand matrix may migrate to the karst habitat and adapt to distinct habitat under the circumstance of ecological selection (Wang et al., 2013).
Parsimony network analysis (except for the nuclear gene GTP) suggested that there are two star-like evolutionary units separately dominating two lineages (Figures 2A,C,D). The fewer mutational steps indicated a recent rapid divergence. Additionally, the UPGMA, Structure analysis and PCoA clustering of microsatellite data showed the 14 populations comprise two clades (I and II) (Figure 6A), and with similar patterns in haplotype distributions of the gene PPRC ( Figure 1D), but had a little conflict with cpDNA haplotype distributions ( Figure 1A). The population BA situated at the border of the two clades geographically clustered into Clade II based on microsatellites, while clustered into another group with the seven populations belonging to Clade I based on cpDNA haplotype distributions. This little conflictual relationship may reflect ongoing pollen-dispersed gene flow among nearby populations, insufficient lineage sorting at nuclear loci or shared ancestry due to recent species divergence. Meanwhile, the BEAST-derived trees also revealed that haplotypes were separated into two groups with an ancestral haplotype at the inner node of each group, further indicating a recent divergence.
Based on the BEAST-derived trees, the estimated divergence times within C. segmentifida fell within the Middle Pleistocene, further certifying that this lineage was the product of a recent rapid divergence, which is consistent with the former conclusion that the extant cycad species have been evolving recently (Nagalingum et al., 2011). Climatic oscillations in the Pleistocene are often regarded as a major factor shaping genetic diversity of many extant species (Hewitt, 1996(Hewitt, , 2000(Hewitt, , 2004Avise, 2000;Yessoufou et al., 2014). Although the Pleistocene ice sheet was evident during major glaciations in Europe (Taberlet et al., 1998), a relatively mild unglaciated Quaternary climate occurred in China, except at higher elevations (Weaver et al., 1998). In China, the lower slopes or valleys were not affected during the cooler periods (Qu et al., 2011). Therefore, we propose that the relatively mild Pleistocene climate in Southwest China contributed to the survival and divergence of C. segmentifida.

Different DNA Markers Reveal Different Patterns of Population Dynamics
In this study, cpDNA, three nuclear genes and microsatellites revealed inconsistent population dynamics for C. segmentifida. Based on cpDNA data, results of neutrality test and mismatch analysis suggested that C. segmentifida has not recently experienced population expansion events, which may be caused by long-term geographical isolation or geographical division. However, the Bayesian Skyline Plot revealed that C. segmentifida was in a stable state for a long time, until approximately 0.1 MYA when its population started shrinking. Population dynamics revealed by three nuclear genes were also inconsistent in this study. Nuclear gene GTP revealed that C. segmentifida experienced a recent population expansion. Our results from the nuclear gene PHYP showed that C. segmentifida has had a recent population contraction. For the nuclear gene PPRC, with the exception of the result of mismatch analysis, the neutrality test and Bayesian Skyline Plot revealed C. segmentifida had a recent population expansion. From the results of DNA sequences, two nuclear genes, GTP and PPRC, showed that C. segmentifida experienced population expansion in Marine Isotope Stages 1 (MIS1) (Ogg et al., 2008), namely, the postglacial period while cpDNA and nuclear gene PHYP showed population contraction in MIS1. In addition, BOTTLENECK analysis based on microsatellites showed that populations of C. segmentifida have not experienced a recent bottleneck event, but GWI estimation revealed these populations experienced a historical reduction in population size. Normally, different genes and markers are subjected to different selective pressures, which may generate the above inconsistency.
Distinct population dynamics were detected in different plant taxa as a result of glacial and interglacial climate oscillations. Gymnosperm species, such as C. debaoensis (Zhan et al., 2011), C. simplicipinna (Feng et al., 2014) and C. multipinnata , have experienced population contractions during the most recent glacial period, while Taxus wallichiana (Liu et al., 2013), Cycas revoluta and Cycas taitungensis (Chiang et al., 2009) have experienced population expansions. Improved knowledge about the history dynamic of species will help us to predict how they will react to environmental fluctuations in the future and to propose conservation strategies for species (Porretta et al., 2007;Lyons et al., 2012;Wei et al., 2013).

Conservation Suggestion
In the wild, C. segmentifida is distributed only in the boundaries of Yunnan, Guizhou, and Guangxi provinces. The destruction of its habitats for planting commercial or medicinal crops and massive illegal digging of this species for trading or ornaments has led to a sharp reduction of its population size. Delaying in conservation decision-making would give rise to the danger of extinction for this species. Conservation management decisions must be made rapidly to prevent this endangered species from extinction.
The objective of conservation of threatened species is to maintain their contribution to overall genetic diversity (Montalvo et al., 1997). This study detected a relatively higher level of genetic diversity in C. segmentifida than in some other Cycas species. If an effective population size is greater than 100, it can prevent inbreeding depression. However, the estimated Ne in most populations of C. segmentifida was less than 50, and in many cases less than 10, such as populations BY and JZ ( Table 2). According to the wild field investigation, we found that the habitat of population BM is not destroyed, the regeneration ability of the population is very strong, it has a high population density, and different ages of individuals exist in the population. The HWE test showed that the population BM was in HWE, indicating that it is a free-mating group. It was habitat destruction and massive illegal digging of C. segmentifida that jeopardized the existing populations of this species and its effective population size. Therefore, we suggest that protection zones or plots in the distribution areas of C. segmentifida be established to protect the habitat for this species. In this study, two genetic clusters were detected in C. segmentifida. We proposed that these two genetic clusters could be managed as two evolutionary units and should be given the highest priority protection. The divergence of C. segmentifida and network of its haplotypes indicated that this lineage recently evolved rapidly. We suggest protecting the species in its natural habitat (in situ). Alternatively, ex situ conservation such as seeds or seedling collection from distinct populations that possess pivotal genetic components and reintroduction are needed. Meanwhile, wild introgression from other Cycas species should be avoided. In addition, more efforts should be made to raise local farmers' conservation awareness and a prohibition on deforestation in Cycas distribution areas should be implemented.

AUTHOR CONTRIBUTIONS
XF was in charge of finishing the molecular genetic studies, performing the data analysis and writing the manuscript. JL contributed to materials collection, DNA extraction, and manuscript revision. XG and Y-CC designed the study, collected research materials and drafted the manuscript. All authors read and approved the final manuscript.