1. Introduction

Southwest China experienced glacial-interglacial cycles in the Pleistocene approximately 2.4-0.01 million years ago (MYA) (Zhou et al., 2006; Royden et al., 2008). Therefore, the dry-hot valleys of Southwest China, such as Red River valley, are recognized as potential refugia for some plant species (Guan and Zhou, 1996; Wang et al., 1996).

Extant cycads are composed of the two families (Cycadaceae and Zamiaceae) with ten genera, and are mainly distributed in Asia, Australia, South and Central America and Africa. Approximately 200 (62%) cycads are threatened with extinction (Jian et al., 2006; Hoffmann et al., 2010). There are about 25 Cycas species (21%) distributed in China (Calonje et al., 2016). The range of these species is often limited by habitat destruction and fragmentation, commonly attributed to planting economic crops along with over-collection for food, medicine and ornamentals (Wang et al., 1996). The Convention on International Trade in Endangered Species (CITES) of Wild Fauna and Flora (WFF) gave all cycads 'First Grade' conservation status in China (Xiao et al., 2004). China harbors abundant Cycas diversity, especially in the drainage areas of the Red River, which is considered a secondary diversification center of Cycas (Hill, 2008).

There are more than 14 Cycas species, of which 10 are endemic to the basin region of the Red River (Hill, 2008). Cycas chenii X. Gong & W. Zhou (Section Stangerioides) is a recently described species (Zhou et al., 2015). This species is distributed in Chuxiong and Honghe of Yunnan Province, China, along the upstream drainage areas of the Red River (also called Yuanjiang in China). It occurs on a range of substrates from limestone to shale or schist, which is characteristic of steep slopes at altitudes ranging from 500 m to 1300 m (Zhou et al., 2015). Only six populations have been discovered, four at the northeast side of the Red River and two at the southwest side. To date, the total population size of C. chenii is estimated at less than 500 individuals across its geographic distribution, with all the know populations of the species being far from any protected areas (Zhou et al., 2015). This species is threatened by ongoing land-clearing and over-collection. Narrow range, small population sizes and presence of a potential barrier to gene flow (the Red River) necessitates understanding extent and structure of genetic variation as a prerequisite for working out the appropriate conservation strategy for this species. The organelle DNA of cycads is maternally inherited and is transmitted only by seeds while their nuclear DNA is biparentally inherited and is transmitted via both seeds and pollens (Huang et al., 2004; Zhong et al., 2011; Feng et al., 2014). Therefore, using these two genetic markers together allows a greater understanding of the role that seed vs. pollen flow play in spatial structuring of genetic variation. Thus, we employed both maternally inherited chloroplast DNA (cpDNA) and biparentally inherited nuclear DNA (nDNA) markers to investigate (ⅰ) the extent and structure of genetic variation in C. chenii, and (ⅱ) the role of the Red River in shaping this structure. Based on our results, we propose suitable conservation strategy for the species.

2. Materials and methods 2.1. Study species and population sampling

All C. chenii populations have less than 100 individuals. All currently known populations of C. chenii were sampled during August to September in 2012. Four populations (ZS, DT, WJ and SP) were sampled from the northeast of the Red River (NE group) and the other two (ML and LH) were from the southwest of the Red River (SW group) (Fig. 1). The distance between sampled individuals was at least 5 m, increasing the likelihood of sampling genetically unrelated individuals. Fresh leaves were dried in silica gel after collection and stored at room temperature until DNA extraction. Voucher specimens were stored in the herbarium of the Kunming Institute of Botany, Chinese Academy of Sciences (KUN).

2.2. DNA extraction, PCR amplification and DNA sequencing

We extracted genomic DNA from dried leaves using the modified CTAB method (Doyle, 1991). After preliminary screening of DNA fragments from universal chloroplast and nuclear primers, we chose two cpDNA intergenic spacers, psbA-trnH and trnL-trnF (Taberlet et al., 1991; Shaw et al., 2005), and two nuclear genes, the phytochrome P gene PHYP (Zhou et al., 2015) and the gene that encodes the largest subunit of RNA polymerase Ⅱ, RBP-1 (Liu et al., 2015). PCR amplification was carried out in 30 μL reactions. For cpDNA, the PCR reactions contained 10 ng of DNA, 3.0 μL of 10× PCR buffer, 1.5 μL of dNTPs (10 mM), 1.5 μL of MgCl₂ (25 mM), 0.45 μL of Taq DNA Polymerase (5 U/μL), 0.45 μL of each primer, 1.5 μL of DMSO (20 mg/mL) and 19.65 μL of double-distilled water. For nDNA, the PCR reactions contained 20 ng DNA, 3.0 μL 10× PCR buffer, 1.5 μL dNTPs (10 mM), 2.25 μL MgCl₂ (25 mM), 0.45 μL Taq DNA Polymerase (5 U/μL), 0.525 μL of each primer, 0.75 μL DMSO (20 mg/mL) and 18 μL double-distilled water. PCR amplification for cpDNA included an initial denaturation stage for 5 min at 80 ℃, followed by 30 cycles of 1 min at 95 ℃, 1 min annealing at 50 ℃, extension for 1.5 min at 65 ℃, and a final extension at 65 ℃ for 10 min. For nDNA: an initial denaturation stage at 94 ℃ for 5 min, was followed by 35 cycles at 95 ℃ for 1 min, annealing at 55 ℃ for 1 min, extension at 72 ℃ for 2 min, and a final extension for 10 min at 72 ℃. All PCR products were sequenced in both directions with the same primers for the amplification reactions, using an ABI 3770 automated sequencer at Shanghai Majorbio Bio-pharm Technology Company Ltd.

2.3. Data analysis

We edited and assembled sequences using SeqMan (Swindell and Plasterer, 1997). Multiple alignments of the DNA sequences were performed in Clustal X, version 1.83 (Thompson et al., 1997), then the DNA sequences were adjusted by Bioedit, version 7.0.4.1 (Hall, 1999). Two cpDNA regions were combined by PAUP* 4.0b10 (Swofford, 2002). The concatenated sequence was used in the following analyses. For the two nuclear genes, heterozygous sites were resolved by applying the PHASE algorithm of DnaSP version 5.0 (Rozas et al., 2003). This program was also used for identification of haplotypes from the aligned DNA sequences and for calculation of Nei's nucleotide diversity (Pi) and haplotype diversity (Hd). Diversity and differentiation parameters (within-population diversity, Hs; total diversity, H_T; differentiation for unordered and ordered alleles, G_ST and N_ST respectively), and a test whether N_ST is larger than G_ST, indicative of a phylogeographic structure (a situation when closely related haplotypes are more often found in the same area than less closely related haplotypes) were calculated with PERMUT (http://www.pierroton.intra.fr/genetics/labo/Software). Significance of the difference between N_ST and G_ST was assessed with 1000 random permutations following Burban et al. (1999). The hierarchical analysis of molecular variance (Excoffier et al., 1992) as implemented in Arlequin (Schneider et al., 2000) was used to estimate among-groups, among-populations within groups and within populations variance components. Isolation by distance (IBD) was tested between all pairs of populations as a correlation between genetic and geographic distance by computing Mantel test using GenAlEx version 6.3 (Peakall and Smouse, 2006). We calculated the ratio of pollen flow to seed flow following the formula pollen/seed migration ratio=[2(1/Φ_STc-1)-(1/ Φ_STn-1)]/(1-1/Φ_STc), where Φ_STn and Φ_STc are levels of amongpopulation differentiation calculated from nuclear and chloroplast markers, respectively (Mousadik and Petit, 1996).

We inferred phylogenetic relationships among cpDNA and nDNA haplotypes using Bayesian approach implemented in MrBayes, version 3.2.1 (Ronquist et al., 2012), in which four simultaneous runs with four Markov chains each were run for 105 generations and trees were sampled every 100 generations, with the first 25% trees from each run being discarded. The nucleotide substitution model used was GTR. Phylogeographic relationships of haplotypes were inferred by statistical parsimony separately for cpDNA and nDNA data using NETWORK 4.2.0.1 software (Forster et al., 2007). Indels were treated as single mutational events in the Network analysis.

To estimate coalescent time between lineages, we used the evolutionary rates 1.01 × 10^-9 and 5.1-7.1 × 10^-9 mutation per site per year for synonymous sites for cpDNA and nDNA (Graur and Li, 2000). Estimation of the time of divergence was performed by BEAST, version 1.6.1 (Drummond and Rambaut, 2007) using the HKY + G and HKY model for cpDNA and two nDNA fragments, respectively, chosen by model-test in MEGA 6.06 (Tamura et al., 2013), and a strict molecular clock. The BEAST program was also used to perform a Bayesian skyline plot analysis to infer the historical demography of C. chenii. Posterior estimates of the mutation rate and time of divergence were obtained by Markov Chain Monte Carlo (MCMC) analysis. The analysis was run for 10⁷ iterations with a burn-in of 10⁴. Genealogies and model parameters were sampled every 10⁴ iterations. Convergence of parameters and mixing of chains were followed by visual inspection of parameter trend lines and checking of effective sampling size (ESS) values in three preruns. The ESS parameter was found to surpass 200, which suggested acceptable mixing and sufficient sampling. Adequate sampling and convergence to the stationary distribution were checked using TRACER, version 1.5 (Drummond and Rambaut, 2004). The pairwise mismatch distributions were examined in DnaSP. The sum-of-squared deviations (SSD) between the observed and expected mismatch distributions were computed, and P-values were calculated as the proportion of simulations producing a larger SSD than the observed SSD. We also used Arlequin, version 3.11 (Excoffier et al., 2005) to calculate the raggedness index and its significance to quantify the smoothness of the observed mismatch distribution. DnaSP was used to examine neutrality tests, Tajima's D (Tajima, 1989) and Fu & Li's F* (Fu, 1997), for detecting departures from population equilibrium.

3. Results 3.1. Genetic diversity and differentiation

Combined, the two cpDNA fragments, psbA-trnH and trnL-trnF, comprised 1227 positions, of which four were nucleotide substitutions and four were indels, resulting in five chloroplast haplotypes (Table 1). Of those, three haplotypes (HapC3, HapC4 and HapC5) were population-specific (ML, SP and WJ), whereas haplotype HapC1 and HapC2 were detected in two populations each (Table 1, Fig. 1a). The nuclear gene PHYP had a length of 883 bp with 18 nucleotide substitutions and one indel. Fifteen haplotypes were detected, of which two (HapP2 and HapP3) were shared by four populations, and three (HapP1, HapP5 and HapP10) were shared by two populations. The remaining ten haplotypes were unique (Table 1, Fig. 1b). The nuclear gene RBP-1 had a length of 918 bp with 30 nucleotide substitutions and one indel. Of the 18 detected haplotypes, HapR2 was the most widely distributed haplotype (Table 1, Fig. 1c).

The within-population variation, H_s, was close to zero for cpDNA data, while exceeding 0.5 for both nuclear genes (Table 2). U tests showed that N_ST was not significantly greater than G_ST (G_ST=0.936, N_ST=0.984, P > 0.05), which indicated that there was no correlation between haplotype similarities and their geographic distribution in C. chenii.

AMOVA revealed substantial differentiation between the two groups (NE and SW) for cpDNA (F_CT=0.76), but no group differentiation for either of the two nuclear genes (F_CT=0, for both nuclear genes). Mantel tests revealed no isolation by distance (IBD) (Fig. 2). The pollen/seed migration ratios were 111.7 and 197 for PHYP and RBP-1, respectively.

3.2. Phylogeny and divergent time of haplotypes

The cpDNA phylogenetic tree revealed two sub-clades corresponding to the two sides of the Red River, comprising three haplotypes (NE group), and two haplotypes (SW group), respectively (Fig. 3a). However, no clear clade structure was revealed from nDNA data.

The two clades (NE and SW groups) split at about 1.175 MYA (0.22-2.678 MYA, 95% HPD), while haplotypes within the two groups diverged much later (from 0.092 MYA to 0.342 MYA) (Fig. 3a). These results imply that haplotypes of C. chenii diverged in the Pleistocene.

3.3. Demographic analysis

The Bayesian Skyline Plots produced from the cpDNA, PHYP and RBP-1 were in disagreement. For cpDNA, it showed a slow decline in the population size until approximate 1000 years ago, at which point a slight expansion occurred (Fig. 4a). For PHYP, C. chenii had a long history of constant population size, followed by a population expansion (about 25, 000) (Fig. 4b). And for RBP-1, it showed that the species population size experienced a long period of continuous growth from 400, 000-50, 000 years ago, followed bya decline from about 50, 000 years ago to the present (Fig. 4c).

The mismatch analysis of cpDNA data revealed a multimodal pattern (Fig. 5a) with significantly positive SSD and raggedness index (Table 2), which indicated that C. chenii did not undergo a recent population expansion. This conclusion was also supported by the positive values of Tajima's D and Fu and Li' F* (Table 2). The mismatch analysis of nDNA data also showed a multimodal pattern but non-significantly positive SSD and raggedness index. Tajima's D and Fu and Li' F* values were negative for PHYP but not for the RBP-1 gene. Together, the results of mismatch analysis and the neutrality test applied to nDNA data indicate that C. chenii did not experience a recent population expansion (Fig. 5b-c; Table 2).

4. Discussion 4.1. Genetic variation

The high level of total genetic diversity in C. chenii based on cpDNA data (H_T=0.92) was comparable with what is commonly observed in other Cycas species: mean H_Tof 0.67 deduced from 170 Cycas species (Petit et al., 2005), H_T=1.000 for Cycas simplicipinna (Feng et al., 2014), H_T=0.896 for Cycas multipinnata (Gong et al., 2015) and H_T=0.56 for Cycas debaoensis (Zhan et al., 2011). The total genetic diversity of C. chenii was also high for both nuclear markers (PHYP and RBP-1) (H_T=0.85 and 0.662, respectively) (Table 2). The within-population diversity of cpDNA was low for the organelle markers (H_s=0.059) and high for nuclear markers (H_s=0.52 and 0.47). The organelle DNA is maternally inherited in Cycas and dispersed only by seeds, whereas nuclear DNA is biparentally inherited and dispersed by seeds and pollens (Huang et al., 2004). The seeds of Cycas are usuallylarge and heavy, falling near the mother plant. Therefore, the limited seed dispersal capacity of Cycas may be the main cause of high population genetic structuring as shown by cpDNA (Hamrick and Godt, 1990). Our study conforms to generally observed low within-population variation and high variation in genetic differentiation among populations of cycads (Walters and Decker-Walters, 1991), e.g.Cycas pectinata (H_S=0.077, G_ST=0.387), C. debaoensis (H_S=0.179, G_ST=0.684), C. multipinnata (H_S=0.225, G_ST=0.749), Cycas dolichophylla(H_S=0.32, G_ST=0.678)(Yang and Meerow, 1996; Zhan et al., 2011; Gong et al., 2015; Zheng et al., 2016).

Analysis of cpDNA data revealed two distinct genetic clades in C. chenii that correspond to NE and SW sides of the Red River, but also no difference in structuring of genetic diversity for ordered vs. unordered alleles (i.e. N_ST ≈ G_ST), and no correlation between geographic and genetic distances. The latter indicate a limited gene flow among the last remaining and highly isolated populations of C. chenii. The physiographic pattern of mountains dissected by deep valleys in the Red River Fault (RRF) zone (Li et al., 2008) could limit seed dispersal in a northeast-southwest direction, thus promoting the population isolation and differentiation of the NE group and SW group in C. chenii. Another factor contributing to the observed structure of genetic variation in C. chenii appears to be high fragmentation of the species to a large extent due to destroyed habitat. In sum, the present study suggests limited gene flow in C. chenii that result from ⅰ) vicariance (the Red River) and ⅱ) habitat destruction and fragmentation. These two processes are responsible for the existence of the two genetic clusters, which, however, display a mosaic-like genetic structure within each of the two parts of the species range (NE and SW groups), high genetic diversity at the species level and low genetic diversity within populations.

The nDNA data revealed lower differentiation between the NE and SW groups than the cpDNA (Fig. 3b-c), and the migration ratio of pollen/seed was higher than 100 for both nuclear genes. These results suggest much higher importance of pollen than seed flow in C. chenii, indicating an existing pollen flow not only among the populations within, but also between the NE and SW groups.

4.2. Demographic history

Some Gymnosperm species experienced population expansion during the recent glacial periods, such as Cycas revoluta (Chiang et al., 2009) and Taxus wallichiana (Liu et al., 2013). The mismatch analysis and neutrality tests showed that C. chenii did not experience a recent population expansion, and the Bayesian Skyline Plot based on cpDNA showed that C. chenii had experienced a slow range reduction since approximately 50, 000 years ago, which is similar to the population dynamics of C. debaoensis, C. simplicipinna and C. multipinnata (Zhan et al., 2011; Feng et al., 2014; Gong et al., 2015). The cpDNA data revealed a long period of slow decline in C. chenii until approximately 1000 years ago, at which point a slight population expansion occurred (Fig. 4a). During the recent 1000 years, the established warmer climate was able to promote population growth in C. chenii. The population dynamics revealed by the two nuclear genes differed. The gene RBP-1 showed a population expansion during the Quaternary, followed by a decline over the Last Glacial Maximum (LGM) (Fig. 4c); in contrast, the PHYP gene showed an expansion in C. chenii population about 25, 000 years ago (Fig. 4b). This contradiction in population dynamics revealed by the cpDNA and nDNA genes could possibly be explained by the different inheritance of the two genomes.

Although details are not clear, C. chenii appears to have experienced slow population contraction during the glacial period. The latter can explain the observed low within-population genetic variation in C. chenii because the larger genetic loss results from the slower range contraction or shift (Arenas et al., 2012). More widely distributed in the basin of the Red River before the glacial epoch, this species was probably forced into several isolated dry-hot Red River valley refugia during glaciation.

The deduced divergence time of the two cpDNA lineages (NE and SW) mainly fall into the Calabrian (1.175 MYA, Fig. 3a). In the Quaternary, the climate oscillated repeatedly from 2.4 MYA to the present, and the dry-hot valley of southwest China could serve at that time as refugia for many plant species suffering range contraction due to glaciations (Guan and Zhou, 1996; Wang et al., 1996).

4.3. Implications for conservation

The present decline of C. chenii throughout its distribution range is known to be largely caused by over-collection and habitat destruction. Designing a suitable conservation plan requires knowledge of the extent and structure of species genetic variation and demographic history. The results show that C. chenii has relatively high genetic diversity at the species level, low genetic diversity within populations and high genetic differentiation among populations. The last two features appear to be the result of range contraction during the species' evolution rather than recent habitat loss. Presence of a clear phylogeographic structure, i.e. two haplotype clades separated by the Red River, implies that conservation efforts cannot focus on one part of the distribution range. In order to preserve present genetic diversity, at least two protected areas must be established, representing NE and SW groups. Because of high population genetic differentiation and the presence of many unique haplotypes, every existing population is important. Populations which can-not be protected in situ, and those harboring the highest diversity and unique haplotypes (such as ML, WJ and SP) must be the highest priority for ex situconservation. Representative seed collection must be done in all existing populations.

Acknowledgments

This research was supported by the United Fund of the NSFC and the Yunnan Natural Science Foundation (Grant No. U1136602 to X. G.). We thank Wei Zhou, Jian Liu and Meng-meng Guan for their assistance with field sampling.