Genetic diversity, genetic structure and demographic history of Cycas simplicipinna(Cycadaceae) assessed by DNA sequences and SSR markers

Study species

A total of 118 individual samples were collected from seven populations of C. simplicipinna (four populations were sampled in Yunnan Province, China and three populations were sampled in Laos). Of the 118 samples, 86 individuals from the seven populations were used for chloroplast and nuclear DNA sequencing. The population known as BOL was eliminated from SSR analysis because there were only 3 individuals. A total of 115 individuals from six populations were used for the microsatellite study. Information on each sampling location and the number of individuals from each population that were used in DNA sequences and SSR analyses are presented in Table 1 and Figure 1, respectively.

Population code	Population location	Latitude (N°)	Longitude (E°)	Altitude (m)	Individuals for DNA sequences/SSR (n)	cpDNA			nrDNA
						Haplotypes (No.)	Hd	Pi × 103	Haplotypes (No.)	Hd	Pi × 103
BOL	Bolikhamxay, Laos	18.456	103.802	140	3/0	Hap A (3)	0	0	Hap 1 (3)	0	0
LUA	LuangPrabang, Laos	19.897	102.161	300	20/25	Hap F (20)	0	0	Hap 5 (20)	0	0
LU	LuangPrabang, Laos	19.831	102.144	280	20/23	Hap E (20)	0	0	Hap 5 (20)	0	0
MM	Mengman town, Yunnan province	21.128	101.315	640	20/20	Hap B (20)	0	0	Hap 2 (20)	0	0
NZD	Nuozhadu Hydropower Station, Yunnan province	22.690	100.419	780	12/12	Hap C(5) Hap D(7)	0.530	0.37	Hap 3 (12) Hap 4 (9)	0.514	0.95
ML	Menglun town, Yunnan province	21.932	101.253	550	11/11	Hap G (11)	0	0	Hap 2 (20)	0	0
NBH	Nature reserve of Nabanhe, Yunnan province	22.166	100.663	694	10/24	Hap H (10)	0	0	Hap 3 (10)	0	0
Total	7				96/115		0.864	2.59		0.723	8.00

Molecular procedures

Young and healthy leaves were collected and dried immediately in silica gel for DNA extraction. Genomic DNA was extracted from dried leaves using the modified CTAB method [23]. After preliminary screening of 21–28 samples (representing approximately 3–4 individuals from each population) with universal chloroplast and nuclear primers, we chose two cpDNA intergenic spacers, psbA-trnH [24] and trnL-trnF [25], and one nrDNA internal transcribed spacer, ITS4-ITS5 [26], for complete analysis. The three pairs of fragments were amplified for the most polymorphic sites of the 86 individuals. PCR amplification was carried out in 40 μL reactions. For cpDNA, the PCR reactions contained 20 ng DNA, 2.0 μL MgCl₂ (25 mM), 2.0 μL dNTPs (10 mM), 4.0 μL 10 × PCR buffer, 0.6 μL of each primer, 0.6 μL Taq DNA polymerase (5 U/μL) (Takara, Shiga, Japan) and 26 μL double-distilled water. For nrDNA, the PCR reactions contained 40 ng DNA, 2.4 μL MgCl₂ (25 mM), 2.0 μL dNTPs (10 mM), 2.0 DMSO, 4.0 μL 10 × PCR buffer, 0.7 μL of each primer, 0.7 μL Taq DNA polymerase (5 U/μL) (Takara, Shiga, Japan) and 24.6 μL double-distilled water. PCR amplifications were performed in a thermocycler under the following conditions: an initial 5 min denaturation at 80°C, followed by 29 cycles of 1 min at 95°C, 1 min annealing at 50°C, and a 1.5 min extension at 65°C, and a final extension for 5 min at 65°C for cpDNA intergenic spacers. For nrDNA sequences we used an initial 4 min denaturation at 94°C, which was followed by 29 cycles of 45 s at 94°C, 1 min annealing at 50°C, and a 1.5 min extension at 72°C, and a final extension for 9 min at 72°C. All PCR products were sequenced in both directions with the same primers for the amplification reactions, using an ABI 3770 automated sequencer at Shanghai Sangon Biological Engineering Technology & Services Company Ltd. For nrDNA, we cloned individuals which had one or more heterozygous sites in the first sequencing round. Six to ten clones were randomly selected and sequenced until the heterozygous site split into two alleles.

Microsatellite markers were selected from recently developed nuclear microsatellites in Cycas[27–33]. PCR amplification was carried out in a 20 μL reaction, containing 10 ng DNA, 1.5 μL MgCl₂ (25 mM), 1 μL dNTPs (10 mM), 1.5 μL 10 × PCR buffer, 0.6 μL of each primer, 0.16 μL Taq DNA polymerase (5 U/μL) (Takara, Shiga, Japan) and 12.14 μL double-distilled water. PCR amplifications were performed in a thermocycler under the following conditions: an initial 4 min denaturation at 94°C, which was followed by 29 cycles of 40 s each at 94°C, 25 s annealing at 48–60°C, and a 30 s extension at 72°C, and a final extension for 10 min at 72°C. PCR products were checked with 8% non-denaturing polyacrylamide gel electrophoresis. Then, we made preliminary screening microsatellite loci for C. simplicipinna. The selected microsatellite loci were stained with a fluorescent dye at the 5' end, their PCR products were separated and visualized using an ABI 3770 automated sequencer, and their profiles were read with the GeneMapper software. An individual was declared null (nonamplifying) at a locus and was treated as missing data after two or more amplification failures. Finally, we chose polymorphic microsatellite loci for C. simplicipinna after calculating polymorphism indices.

Data analysis

Data analysis of SSR markers

Dataset editing and formatting was performed in GenAlEx, version 6.3 [49]. We tested for evidence of preliminarily selection of our selected loci because our microsatellites had been derived from recently developed nuclear microsatellites of Cycas. We also used the Fst-outlier approach to test for signs of positive and balancing selection on those loci [50, 51] by LOSITAN [52]. The outlier loci were identified by the expected distribution of Wright’s inbreeding coefficient Fst compared with H_E[53]. As recommended by Antao [52], we ran LOSITAN to identify the loci under neutral selection by using the infinite allele model and 10,000 simulations. Twenty microsatellites were first selected after detecting the levels of genetic diversity in the sample of 115 individuals of C. simplicipinna in the six populations. The results of positive and balancing selection on the twenty microsatellites detected balancing selection on locus A16 and positive selection on four other loci (A3, A9, A13, and A14). However, locus A13 did not reach the significant level of an Fst-outlier (Figure 2). Therefore, four loci (A3, A9, A14, and A16) with significant levels as Fst-outliers were removed from further analysis. Finally, we selected sixteen microsatellites with high polymorphism, stability, and conformity with neutral selection for our research (Additional file 1: Table S1).

The number of alleles (N_A), private alleles (A_P), effective number of alleles (N_E), expected heterozygosity (H_E = 1-∑Pi², Pi, population allele frequencies), observed heterozygosity (H_O = No. of Hets/N), information index (I), and fixation index (F = 1-(H_O/H_E)) were calculated using GenAlEx, version 6.3 [49], and POPGENE, version 1.32 [54], with mutual correction. Allelic richness (A_R) was estimated with FSTAT, version 2.9.3 [55], and percentage of polymorphic loci (PPB) was calculated with GenAlEx, version 6.3 [49]. Differentiation between pairs of populations was computed using F_ST and tested with GenAlEx, version 6.3 [49]. Isolation by distance (IBD) was tested on SSR data by computing Mantel tests in GenAlEx, version 6.3 [49] using a correlation of F_ST/(1-F_ST) with geographic distance for all pairs of populations. F_ST/(1-F_ST) was caculated with Genepop, version 4.1.4 [56]. Gene flow between pairs of populations was estimated using Wright’s principles Nm = (1-F_ST)/4F_ST[57]. Hardy-Weinberg equilibrium (HWE) was tested for each locus and each population using Genepop, version 4.1.4 [56].

The genetic structures of sampled populations and individuals were estimated by unweighted pair group mean analysis (UPGMA) using TEPGA, version 1.3 [58], with 5,000 of permutations. An individual-based principal coordinate analysis (PCO) was visualized by the program MVSP, version 3.12 [59], using genetic distances among SSR phenotypes. We also conducted a Bayesian analysis of population structure on the SSR data using STRUCTURE, version 2.2 [60]. Ten independent runs were performed for each set, with values of K ranging from 1 to 6, a burn-in of 1 × 10⁵ iterations and 1 × 10⁵ subsequent MCMC steps. The combination of an admixture and a correlated-allele frequencies model was used for the analysis. The second-order rate of change of the log probability of the data with respect to the number of clusters (ΔK) was used as an additional estimator of the most likely number of genetic clusters [61]. The best-fit number of grouping was evaluated using ΔK by STRUCTURE HARVESTER, version 0.6.8 [62]. Finally, we identified geographical locations where major genetic barriers among populations might occur with a barrier boundary analysis, using BARRIER, version 2.2 [63], based on genetic distance matrices.

We calculated the effective population sizes of each population to establish the degree of endangerment of the species. We used the program LDNe at three levels of the lowest allele frequency (=0.01, 0.02, 0.05) at a 95% confidence interval [64]. We tested the bottleneck statistic at the population level to explore the demographic history of populations by using different models and testing methods implemented in BOTTLENECK, version 1.2.02 [65]. The computation was performed under a stepwise mutation model (SMM) and a two-phased model (TPM). We did not use the standardized differences test in this study because the test was usually used at the condition of having at least twenty polymorphic loci. Two other methods (Sign tests and Wilcoxon tests) were applied to the two models. We also used a mode shift model [66] to test for bottlenecks in each population. These methods implemented in BOTTLENECK have low power unless the decline is greater than 90% [66]. They are most powerful when bottlenecks are severe and recent [67]. In addition, a genetic bottleneck was further investigated with the Garza-Williamsion index (also called M-ratio [68], the ratio of number of alleles to range in allele size). When seven or more loci are analyzed, the Garza-Williamsion index is lower than the critical Mc value of 0.68, a value obtained by simulations based on the empirical data in bottlenecked populations, suggesting a reduction in population size [40, 68]. The Garza-Williamsion index is more powerful to detect genetic bottlenecks if the bottleneck lasted several generations or if the population made a rapid demographic recovery [67]. The index was analyzed by Arlequin, version 3.11 [40].

DNA sequences

The combined length of cpDNA (psbA-trnH and trnL-trnF) varied from 1,408 to 1,438 bp and aligned with a 1,452 bp consensus length that contained 14 polymorphic sites and 16 indels (Additional file 2: Table S2). A total of eight chloroplast haplotypes was identified, and each population was fixed for one particular haplotype, except for population NZD, in which two unique haplotypes was detected (Table 1). The aligned nrDNA (ITS4-ITS5) matrix ranged from 1,079 to 1,087 bp with a consensus length of 1,100 bp that contained 32 polymorphic sites and 11 indels (Additional file 3: Table S3). A total of five nuclear haplotypes was derived. Population BOL had one unique haplotype (Hap 1), MM and ML shared haplotype 2, LUA and LU shared haplotype 5, and NZD had two haplotypes (one was unique and another shared with NBH) (Table 1).

A phylogeny of cpDNA and nrDNA haplotypes was constructed by both maximum parsimony (MP) and Bayesian methods, using C. diannanensis as an outgroup. Both analyses produced phylogenetic trees with consistent topologies (Figure 3). Eight cpDNA haplotypes appeared as a comb-like structure because they lacked enough information sites (Figure 3, a). Five nrDNA haplotypes were clustered into three clades, showing that Hap 2 is more closely related to Hap 5, and Hap 3 is more closely related to Hap 4 (Figure 3, b). The neighbor-joining trees (NJ) supported the congruent phylogenetic relationship of the cpDNA and nrDNA haplotypes (Figure 4). The haplotype network analysis of cpDNA and nrDNA also yielded the same topological relationships (Figure 5). Most haplotypes were distributed in the outside nodes of the reticulate evolutionary diagram, and many missing haplotypes, specifically between Hap 1 and Hap 2, were evident in the reticulate evolutionary diagram of the nrDNA haplotypes (Figure 5, b).

We derived the estimated time of divergence of C. simplicipinna with the Bayesian method, using BEAST, version 1.6.1 [46]. The estimated time of divergence ranged from 0.276 MYA to 2.682 MYA according to the cpDNA data and 0.135 MYA to 1.429 MYA according to the nrDNA data (Figure 4). The cpDNA haplotype G (Hap G) was the earliest to diverge. Its time of divergence was estimated to have been 2.682 MYA. The time of divergence of the clade comprising Hap A, E, F, and B and the clade comprising Hap H, C, and D was 1.090 MYA (Figure 4, a). The phylogenetic tree of nrDNA shows that Hap 1 was the earliest haplotype to diverge. Its time of divergence was 1.429 MYA. The time of divergence between the clade comprising Hap 2 and 5 and the clade comprising Hap 3 and 4 was 0.935 MYA (Figure 4, b). These results imply that the C. simplicipinna haplotypes were diverged during the Pleistocene (2.6 Ma to 11 ka).

Population dynamic analysis using cpDNA and nrDNA data showed that the population demography of C. simplicipinna was stable until approximately 50,000 years ago, at which time a contraction event occurred (Figure 6). The results of the mismatch analysis for all C. simplicipinna populations displayed a multimodal distribution pattern (Figure 7) with significant SSD and raggedness values (Table 4), which indicates that C. simplicipinna has not undergone a recent population expansion. This conclusion is also supported by the results of the Neutrality Test, Fu’s F _S , which yielded positive values (Table 4). Based on a Bayesian simulation, the skyline plot showed recent declines in population size of all populations of C. simplicipinna during Quaternary glaciations and no subsequent expansion (Figure 6).

SSR data

A total of 169 alleles were identified at the sixteen loci. Diversity estimates varied in different populations (Table 5). Allelic richness was lowest in population MM (A_R, 2.628) and highest in population LUA (A_R, 5.014). The number of alleles (N_A) ranged from 2.875 to 6.063, the number of private alleles (A_P) ranged from 1 to 14, the effective number of alleles (N_E) ranged from 1.925 to 3.521, the information index (I) ranged from 0.635 to 1.268, observed heterozygosity (H_O) ranged from 0.306 to 0.473, and expected heterozygosity (H_E) ranged from 0.353 to 0.603. These indices all showed a similar trend, with the lowest values in MM and the highest values in LUA. Fixation indices (F) were positive for all six populations, with a mean value F = 0.170, which suggests a high level of inbreeding within each population. The percentage of polymorphic loci (PPB) was high, ranging from 75% to 100%. Population MM had the lowest genetic diversity, and LUA had the highest. The genetic differentiation coefficient F_ST varied from 0.036 to 0.467, with a mean value 0.261. No significant effect of isolation by distance (IBD) was detected (Figure 8), as the correlation between genetic and geographic distances was non-significant (P > 0.05), which was supported by the result of Mantel test. Estimates of gene flow between each pair of the six populations are showed in Table 6. Population LUA had the most gene flow with the other populations, and MM had the least. Excesses of homozygotes caused five populations and nine loci to deviate from Hardy-Weinberg equilibrium (Table 5, Additional file 4: Table S4).

Population	N T	N P	A R	N A	N E	I	H O	H E	F	HWE ( P)	PPB(%)	Ne
LUA	97	14	5.014	6.063	3.521	1.268	0.473	0.603	0.197	0.000***	100	164.9
LU	91	6	4.835	5.688	3.402	1.220	0.459	0.589	0.228	0.000***	100	33.9
MM	46	1	2.628	2.875	1.925	0.635	0.306	0.353	0.156	0.000***	81.25	56.6
ML	53	6	3.312	3.313	2.310	0.773	0.330	0.418	0.190	0.002**	75.00	24.8
NZD	64	6	3.889	4.000	2.385	0.885	0.441	0.457	0.093	0.477ns	93.75	43.7
NBH	63	8	3.417	3.938	2.292	0.846	0.387	0.442	0.154	0.000***	93.75	106.1
Mean	69	6.83	3.849	4.313	2.639	0.938	0.394	0.447	0.170		90.63	71.7

Note: N_T, number of total alleles; N_P, number of private alleles; AR, allelic richness; N_A, number of alleles; N_E, effective number of alleles; I, information index; H_O, observed heterozygosity; H_E, expected heterozygosity; F, fixation index; HWE, Hardy-Weinberg equilibrium; PPB, percentage of polymorphic loci; Ne, effective population size.-, Monomorphic; *, P < 0.05; **, P < 0.01; ***, P < 0.001.

The STUCTURE analysis, using the ΔK method, showed that the optimal K value was K = 3 (Figure 9), which showed that the six populations were clustered into three groups. Populations LUA and LU were grouped into one cluster (Cluster I), MM and ML were grouped into another cluster (Cluster II), and NZD and NBH were grouped into a third cluster (Cluster III). The result of K = 6 was also present here to detect whether or not has further subdivision in the species. From the Figure 9 we can see that there is only further subdivision at K = 6 between the population LUA and LU. In contrast with K = 6, it is clear that K value was K = 3 is a better solution, because the existence of three groups was also supported by the PCO analysis (Figure 10). Two-dimensional PCO separated all individuals into three clusters along the two axes. The dendrogram (Additional file 5: Figure S1) obtained with the UPGMA clustering method showed that the six populations were separated into three clades with high bootstrap values (100). It is the same as STRUCTURE (K = 3) and PCO analysis. In the UPGMA clustering dendrogram, populations LUA, LU, MM, and ML were clustered into one large clade with a bootstrap value of 78.7. The BARRIER analysis showed that there was only one major genetic boundary (Barrier I), with a 52.7% mean bootstrap value, separating the six populations into two clusters (Figure 11).

Genetic variation and genetic structure

The genetic variation of a species is a product of its long-term evolution and represents its evolutionary potential for survival and development [69, 70]. Cycads, as ancient gymnosperms with millions of years of evolutionary history, a long life cycle, and overlapping generations, would be expected to have genomes that are responsive to different selective pressures. High levels of genetic variation would be expected to have accumulated during a long evolutionary history. As expected, we found that C. simplicipinna has high genetic diversity (Table 1, 2 and 5) at a species level compared with other species of Cycas by using similar markers e.g., an average value of H_T = 0.564 and Pi = 0.00132 were reported for two markers of type cpDNA in C. debaoensis[5], and an average value of H_O = 0.349 and H_E = 0.545 and the maximum value of Ap = 2.1, N_A = 5.8 were reported for 14 markers of type EST-microsatellites in C. micronesica[53]. Cycas simplicipinna also has higher genetic diversity than many conifers. Many individual conifer species show lower genetic diversity, e.g., an average value of H_T = 0.234 and Hs = 0.190 were reported for two markers of type cpDNA in Pinus tabulaeformis[71], an average value of π = 0.000573 and π = 0.006131 were reported for two markers of type cpDNA and one marker of type nDNA in Tsuga dumosa, respectively [72], and an average value of H_T = 0.77, Hs = 0.66, N_R = 3.98, H_E = 0.62 were reported for seven markers of type nuclear microsatellites in Taxus baccata[73]. The mean genetic diversity value of 170 plant species that was estimated from cpDNA-based studies was H_T = 0.67 [74]. However, at a population level, C. simplicipinna shows low genetic diversity; only population NZD has a relatively high genetic diversity.

The genetic diversity of C. simplicipinna among all populations (H_T = 1.000, 0.878 from cpDNA and nrDNA, respectively Table 2) is also higher than the average intra-population diversity (H_S = 0.076, 0.073 from cpDNA and nrDNA, respectively Table 2), which indicates that there are high levels of genetic differentiation among populations (G_ST = 0.924, 0.916, N_ST = 0.985, 0.992 from cpDNA and nrDNA, respectively Table 2). U tests showed that N_ST was not significantly greater than G_ST, suggesting that there is no distinct phylogeographical structure in C. simplicipinna. The F_ST value of C. simplicipinna (nSSR: F_ST = 0.261, G_ST = 0.246, Table 5) was higher than the mean value of outcrossing species (F_ST = 0.22) that was inferred from SSR [75]. Wright [76] had proposed that an F_ST value greater than 0.25 (C. simplicipinna: F_ST = 0.26 > 0.25) would indicate that there was significant genetic differentiation among populations. Additionally, according to the results of deviation from Hardy-Weinberg equilibrium test (Table 5, Additional file 4: Table S4), only population NZD was in Hardy-Weinberg equilibrium. The remaining five populations deviated significantly from Hardy-Weinberg equilibrium, and the fixation indices (F) were greater than zero. We therefore conclude that there is a notable deficit in heterozygosity and severe inbreeding in C. simplicipinna populations, resulting in a high among-population genetic differentiation as a whole.

The genetic structure of C. simplicipinna based on SSR markers showed that the six study populations were divided into three clades (I, II and III). BARRIER analysis (Figure 11) showed that only one barrier (BS > 50%) exists among the six populations, suggesting that clade I and clade II are genetically more closely related to each other than either is to clade III. Genetic structures of C. simplicipinna derived from organelle and nuclear markers are different. It displays high differentiation at cpDNA markers and high distinct structure at biparental markers. This is caused by different features of the markers, such as mutation rates and the biased migration between organelle (pollen migration) and nuclear markers (pollen migration and seed migration, seed migration is very little). Comparisons of genetic with paleoecological data of temperate woody species are known to reveal unique genetic lineages and⁄ or endemic haplotypes in separate refuge populations [77, 78]. Although C. simplicipinna is a tropical or subtropical woody species, the species similarly possesses unique genetic lineages and endemic cpDNA haplotypes in its separate refuge populations. Duminil [79] had proposed that the level of genetic structure in temperate trees and the potential to reflect historical population isolation are determined in part by life history. Specifically, species with large geographical ranges and wide-ranging seed dispersal display low differentiation at maternally inherited cpDNA markers, and long-lived outcrossing species display low structure at biparental markers. However, C. simplicipinna displays the opposite results. This outcome may be due to its having a limited geographic distribution and severe population inbreeding or to its being a tropical or subtropical woody species, which differs from temperate trees.

The correlation between genetic and geographic distances was non-significant (Figure 8), which indicates that there is no significant effect of isolation by distance (IBD). The species lacks of clear geographic structure. It does not automatically mean that the genetic information has no value in directing management. Yang and Meerow [80] estimated that gene flow distance among local populations in Cycas was 2–7 km. However, the distances between extant populations of C. simplicipinna are all greater than 7 km (the geographic distance between LUA and LU was the smallest of all population pairs (7.52 km)). Because effective gene flow higher than 1 is often regarded as high enough to prevent population differentiation due to genetic drift, while gene flow less than 1 may be a major reason caused the genetic differentiation among populations[81], many gene exchanges between populations less than 1 reflect a low level of gene flow (Table 6). Cycas simplicipinna is dioecious and pollenated by insect (Curculionoidea, weevils) [82]. Unlike the birds, weevils can’t spread pollen over a long distance. Its seeds are too large to disperse naturally over such a long distance. Most seeds disperse near the mother plant, which increases inbreeding. Inter-population could not exchange gene flow easily. We therefore conclude that the strong genetic differentiation and structure observed in C. simplicipinna is due mainly to its limited gene flow and severe inbreeding. Our analysis therefore suggests that C. simplicipinna has high genetic diversity at the species level, low genetic diversity within populations, high genetic differentiation among populations and a clear genetic structure. This conclusion is also supported by the AMOVA analysis (Table 3), which shows that almost all of the genetic variation exists among populations (DNA sequences). Our results support the conclusion that low genetic variation within populations is biologically typical for Cycas, unlike other gymnosperms [83].

Phylogeny and demographic history

Phylogenetic trees constructed on the basis of DNA sequences with different criteria in four different software systems all suggested a consistent systematic relationship of haplotypes (Figure 3 and 4). The comb-like structure of cpDNA haplotypes is likely to be the result of insufficient information site due to insufficient evolutionary time [1, 84]. The network analysis (Figure 5) showed that most of the haplotypes were distributed in the outside nodes of the reticulate evolutionary diagrams and there were many missing haplotypes among them. The reason for the observed haplotypes distribution pattern is that diversity within populations is extremely low and differentiation among populations is high. With the exception of the population NZD, the populations have no haplotype and nucleotide polymorphism. This is because that missing haplotypes in the network are due to the species’ long evolutionary history during which climate variations, geological activities, and human activities formed the several scattered populations that currently exist.

Understanding a species’ demographic history aids in understanding its ancient evolutionary environment. Quaternary glaciers are known to have profoundly affected the distribution of plant species [3]. Previous studies have shown that different plant species had different responses to glacial and interglacial influences. Most plant taxa are believed to have shifted the latitude or elevation of their ranges in response to glaciation [8]. Some plant species, e.g., Taxus wallichiana[85], C. revoluta and C. taitungensis[6], experienced population expansion during the most recent glacial period, and others, e.g., C. debaoensis[5], showed population contraction. Although there is growing evidence for population demographic stability or expansion throughout the Last Glacial Maximum (LGM) in a range of different organisms [86–91], C. simplicipinna appears to have exhibited population demographic contraction during the LGM and no later expansion. We infer this from the Bayesian skyline plot (Figure 6), a divergence of the observed mismatch distribution from the expected distribution (Figure 7). The significant positive values of SSD and raggedness and non-significant positive Fu’s Fs (Table 4) also imply C. simplicipinna did not undergo a population expansion, so it maybe undergo a population contraction or stay a population dynamic equilibrium. Bottleneck analysis (Table 7) showed that populations of C. simplicipinna have not experienced recent bottleneck event but a reduction in population size, which was in accord with Bayesian skyline plot (Figure 6). The divergence times of C. simplicipinna haplotypes fall generally in the Pleistocene. We therefore conclude that C. simplicipinna was widely and continuously distributed before the glacials and contracted into several isolated surviving populations during the glacials. Refugium populations typically have relatively high genetic diversity and unique haplotypes [77, 92]. In our study, C. simplicipinna had high genetic diversity at the species level. Each of its population had a unique haplotype for the cpDNA data but low level genetic diversity. The existing C. simplicipinna populations are distributed mainly in the tropics or subtropics. Because they are adapted to warm temperatures, temperature is the main factor affecting their growth. The average annual temperature in Chinese subtropical and tropical areas during the last glacial maximum was 4-6°C lower than it is today, which caused change and migration in vegetation [93]. This could have had a strong influence on the distribution of C. simplicipinna. The low Ice Age temperatures were not suitable for C. simplicipinna, which led to decline or even local extinction of populations. Some cpDNA and nuclear haplotypes were lost in the process, leading us to deduce that the present areas of distribution are its Ice Age glacial refugia. Cycas simplicipinna migrated to the scatter refugia that form its haplotypes current distribution pattern. They did not migrate to a common refugium during the Quaternary glacial period but instead survived in their original locales.

We conclude that the reduction in effective population size and limited gene flow were the main factors promoting genetic differentiation among populations of C. simplicipinna, which in turn led to the current population structure and distribution pattern.

Conservation implications

One goal for the conservation of threatened plants is to maintain the genetic diversity of native plant species [94]. Ex situ collections are also important because they provide a setting for breeding for introduction, which is an important way to increase the genetic diversity of native plant populations. We found that three (NBH, BOL, ML) of our seven study populations were distributed in protection zones in which there has been no overexploitation and habitat destruction. However, the other four populations (NZD, LU, LUA, and MM) occur outside of protected areas and should be protected. Our study, which revealed a clear population structure of C. simplicipinna that show low genetic variation within populations and high genetic differentiation among populations, has significant implications for conservation of this species. The population structure and demographic history of C. simplicipinna imply that conservation efforts cannot focus on only one part of the species’ range. We suggest that habitat protection be strengthened immediately by establishing protection zones or plots in the distribution areas of C. simplicipinna to improve the conservation awareness of local farmers and to prohibit deforestation in Cycas distribution areas. Protection of populations with the highest genetic diversity, such as populations NZD and LUA, should be of highest priority. Priority should also be given to populations with unique haplotypes. Based on the analysis of the two cpDNA sequences in our study, six of the seven sampled populations of C. simplicipinna each had one unique cpDNA haplotype. Therefore, every C. simplicipinna population should be given maximum protection to prevent losing any haplotypes. The reduction of effective population size in the Ice Age has led to a small effective population size for the species as a whole, which jeopardizes the species’ current genetic diversity. Different populations and/or individuals should be moved from current areas of rich genetic diversity to secure remote areas for ex situ conservation. These measures taken together could protect and enrich the genetic diversity of C. simplicipinna.

Availability of supporting data

The data set of the DNA sequencing data in our study are deposited in GenBank under accession numbers KM065478-KM065496.