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.

We extracted total DNA using the modified CTAB method (Doyle, 1991) and then sequenced four cpDNA intergenic spacers, psbA-trnH, trnS-trnG, psbM-trnD, and trnL-trnT, and three nuclear genes, GTP, PHYP, and PPRC, for complete analysis after preliminary screening from universal chloroplast and nuclear primers. PCR amplification procedures were the same as those used in C. simplicipinna (Feng et al., 2014). The PCR products were sequenced on an ABI 3770 automated sequencer at Shanghai Major Biological Medicine and Technology Co Ltd in both directions. All sequences were deposited in GenBank with the accession numbers KU240448-KU240465, KU240469-KU240471, KU240474-KU240488, KU240490-KU240491 (Feng et al., 2016a); KT824880-KT824887, KT824890-KT824891 (Feng et al., 2016b); KY292488-KY292501.

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., 2009, 2010; Ju 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.

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

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⁵ subsequent MCMC samplings and a subsequent burn-in of 1 × 10⁵ 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).

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.

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 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 BEAST-derived 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.

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

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 mutation-drift equilibrium (Table 3). Meanwhile, Mode shift model tests revealed that all of the populations were under a normal L-shaped distribution, indicating that C. segmentifida had not experienced a severe bottleneck recently. However, Garza-Williamson indices of 14 populations were lower than the critical Mc value of 0.68 (Table 3), indicating that C. segmentifida had apparently experienced a population size decline (bottleneck) in history.

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 (Gong et al., 2015), 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.

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 (Gong and Gong, 2016). 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, 2000, 2004; Avise, 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.

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 (Gong et al., 2015), 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).

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.