1. Introduction

Cypripedium L. (Orchidaceae), commonly known as slipper orchids or lady's slipper orchids, is a well-known genus comprising 58 species (Cribb and Green, 1997) mostly distributed in the north temperate zone and extending south to the Himalayan regions and Central America. A total of 36 out of 58 species of Cypripedium are found in China, of which 25 are endemic species (Chen and Cribb, 2009). Many of these species (32) are listed in the Red List of China as threatened (Wang and Xie, 2004). The major threat is harvesting, breakdown in ecological connections (pollinators, mycorrhiza and chlorophyllous hosts), habitat loss and fragmentation (Swarts and Dixon, 2009).

Cypripedium tibeticum King et Rolfe is no exception to this general trend of population decline. For a long time, this beautiful wild flower has been intensively collected from the field for its horticultural and medicinal value. Overharvesting, as well as deterioration of the habitat due to overgrazing has led to the fragmentation of populations and a great decrease in their abundance. C. tibeticum is listed as a vulnerable species in Chinese Species Red List (Wang and Xie, 2004).

Knowledge of extent and structure of species genetic diversity is essential for the establishment of an efficient conservation strategy because genetic factors contribute to species extinction risk through inbreeding depression, loss of genetic diversity and loss of evolutionary potential (Frankham, 2012, Frankham et al., 2014). Giving priority to particular populations in conservation decisions is usually based on their contribution to the total species genetic diversity. For this reason it is important to know a level of a population's divergence from other populations, and its withinpopulation variation. Role of genetic variation in evolutionary context is coined by the concept of evolutionary significant units (ESU) (Ryder, 1986), intraspecific genetic lineages that resulted from historical isolation due to persistent barriers to gene flow. These genetically unique units are high priority targets for conservation (Moritz and Faith, 1998; Crandall et al. 2000). On the other hand, high within-population genetic variation is another important criterion for prioritization of populations for protection. Beside its contribution to the total species genetic diversity, sufficiently high within-population genetic variation is crucial for a species population long-term survival and capability to respond and adapt to the environmental changing environments changes (Avise and Hamrick, 1996; Reed and Frankham, 2003; Neale, 2012).

In order to gain the genetic knowledge needed for population conservation prioritization in C. tibeticum, we studied nine populations representing the entire species distributional range, using two chloroplast regions. The results were used to design an effective conservation strategy for C. tibeticum.

2. Materials and methods 2.1. Study species

C. tibeticum King et Rolfe is a rare perennial terrestrial herb with short creeping rhizome. It is one of the most widespread species of Cypripedium and is native to the East Himalaya-Hengduan Mountains (Chen and Cribb, 2009). This species is usually clustered in sparse forests, forest margins, scrubby slopes and grassy slopes at high altitudes (2300-4200 m). The plant has only one nodding and big flower with purple or dark maroon color (Fig. 1). The flowers are self-compatible but insects are required to deposit pollen on the stigma (Li et al., 2006). The effective pollinators are bumble bee queens. The produced seeds are dust-like and numerous (Wang et al., 2013).

2.2. Plant material

A total of 157 individuals from 9 populations of C. tibeticum were sampled. Fresh young leaves were collected and dried in silica-gel during field expeditions conducted from 2010 to 2016, and maintained at -80 ℃ in the lab. The randoμly chosen sampled individuals were at least 20 m apart. The accession codes, accession numbers, localities and altitude for each of the populations are listed in Table 1.

2.3. DNA extraction, amplification, and sequencing

Total DNA was extracted from silica gel-dried leaves using the OMEGA SP Plant DNA Kit (OMEGA BIO-TEK, Norcross, GA, USA). Two chloroplast fragments, the rps16 intron (Oxelman et al., 1997), the trnL intron and the trnL-trnF intergeneric spacer (trnL-F) (Taberlet et al., 1991) were amplified and directly sequenced for all the individuals. The rps16 was amplified using the two primers (rps16-1F: 5'- CACGGGTCGCCCTCGTTCCG -3', rps16-1R 5'- TCGGGATCGAACATCAATTGCAAC -3'). The trnL-F was amplified using primer "c" (5'- CGA AAT CGG TAG ACG CTA CG -3') and "f" (5'- ATT TGA ACT GGT GAC ACG AG -3'). Polymerase chain reaction (PCR) was carried out in a total volume of 20 μl containing 2 μl of 10 × PCR reaction buffer (Takara, Japan), 1.6 μl of 25 mM MgCl₂, 1 μl of each primer (Sangon, Shanghai, Chana) at 50 ng/μl, 1.6 μl of 2.5 mM dNTP solution in equimolar ratio, 0.1 μl of Taq DNApolymerase (5 units/μl, Takara, Japan) and 2 μl of genomic DNA at 5 ng/μl. The PCR protocol for rps16 intron was as follows: an initial denaturation step at 80 ℃ for 5 min, then 32 cycles of denaturation at 95 ℃ for 1 min, annealing at 59.5 ℃ for 1 min, an extension at 65 ℃ for 4 min, and a final extension at 65 ℃ for 10 min. The PCR protocol for trnL-F was as follows: an initial denaturation step at 94 ℃ for 4 min, then 33 cycles of denaturation at 94 ℃ for 1 min, annealing at 55 ℃ for 1 min, an extension at 72 ℃ for 2 min, and a final extension at 72 ℃ for 10 min. The PCR products were purified and sequenced by Sangon Cooperation, Shanghai, China.

2.4. Molecular data analysis

Sequences were assembled with DNAStar (Gene Codes Corporation, USA). Variable sites in the data matrix were double-checked manually. Multiple alignments of the sequences were performed with ClustalX (Thompson et al., 1997) and subsequently adjusted in BioEdit (Hall, 1999). Haplotype diversity (H_d) (Nei and Tajima, 1983) and nucleotide diversity (Pi, π) (Nei, 1987) among populations were calculated with DnaSP (Librado and Rozas, 2009). Total chloroplast sequence diversity (H_T), sequence diversity within populations (H_S), and the coefficient of genetic differentiation (G_ST) were calculated with HAPLONST (Pons and Petit, 1996). The genetic variability was partitioned into within- and among-population components by an analysis of molecular variance (AMOVA) using Arlequin (Excoffier et al., 2005). Spatial genetic structure of chloroplast haplotypes was analyzed with SAMOVA (Dupanloup et al., 2002). The F_CT index of genetic differentiation among K groups was computed to obtain the best configuration of groups according to F_CT values. In this study, K ranged from 2 to 8, with each simulation starting from 100 random initial conditions and 1000 times permutation. Genealogical relationships between haplotypes were inferred from a maximum parsimony median-joining network using the program TCS (Clement et al., 2000). In order to rank populations in terms of their conservation priority, the contribution of each population to total gene diversity (CT) and total allelic richness (CTR) was computed according to Petit et al. (1998), using the CONTRIB software (available at https://www6.bordeaux-aquitaine.inra.fr/biogeco_eng/Scientific-Production/Computer-software/Contrib-Permut/Contrib). This method estimates contribution of each population to both within-population diversity (CS) and among-population differentiation (CD). Because the population sample sizes varied, we used rarefaction set to the smallest sample size (N = 10). Negative values indicate that the diversity or the differentiation of a population is lower than the mean of the whole dataset.

2.5. Species distribution modeling

We used species distribution modeling to predict the geographic distribution of suitable habitat for C. tibeticum under current climatic conditions. The 19 "Bioclim" variables (Hijmans et al., 2005) summarizing temperature and precipitation dimensions of the environment were obtained from WorldClim (Hijmans et al., 2005) with a resolution of 30" latitude/longitude (ca. 1 km² at the ground level). The MAXENT (Phillips et al., 2006; Phillips and Dudik, 2008) was used to generate an estimate of probability of presence of the species that varies from 0 to 1, where 0 being the lowest and 1 the highest probability. We used 116 occurrence records from GBIF (GBIF.org, 2018), Chinese herbaria and 9 sampled locations. In the analyses, we withheld 25% of the occurrence data for model evaluation, set number of iterations to 500 and used ten replicates under the 'crossvalidate' option. The accuracy of model predictions was tested by calculating the area under the 'Receiver Operating Characteristic (ROC) Curve' (AUC; Fawcett, 2006). The summary map was generated by averaging Maxent outputs.

3. Results 3.1. cpDNA sequence variation and haplotype distribution

The aligned length of the rps16 and trnL-F was 868 and 855 bp, respectively, of which two and seven nucleotides were variable and parsimony-informative (Suppl. Table 1). Indels were treated as missing data. In total 7 haplotypes (H1-H7) were identified in the 157 individuals from nine populations (Suppl. Table 1). All sequences were submitted to GenBank with accession numbers MF279233-MF279239 and MF279240-MF279246 for trnL-F and rps16, respectively.

Spatial distribution and relative frequencies of the haplotypes are shown in Fig. 2B. Haplotype H7 was the most abundant and widely distributed haplotype, found in 69 individuals from 5 populations (43.9% of the total sampled individuals). The haplotypes H5 and H6 were possessed by 38 and 25 individuals, respectively. The haplotype H5 was detected in populations GZ, DC and LJ, and haplotype H6 in populations XC, ZD and JDS. The haplotypes H2, H3 and H4 were "private" haplotypes, represented by 5 individual (36%) in CY, 6 individuals (50%) in DC, and 3 individuals (25%) in DC, respectively. The parsimony network indicated that H7, H5 and H6 are the ancestral relative to the other haplotypes (Fig. 2C).

3.2. Extent and structure of genetic diversity

The total diversity (H_T) and the within-population diversity (H_S) of C. tibeticum were 0.805 ± 0.078 and 0.293 ± 0.095, respectively, with large variation among the populations in the haplotype diversity (H_d) and the nucleotide diversity (π) (Table 1). The population DC had the highest level of haplotype diversity and population CY had the highest level of nucleotide diversity (Table 1). Four populations (GZ, LJ, ML and WC) were monomorphic and the other five populations had from two to three haplotypes.

Population differentiation was high (G_ST = 0.636, HAPLONST; Φ_ST = 0.86, AMOVA). Spatial genetic analysis of cpDNA haplotypes with SAMOVA failed to reveal clear phylogeographic structure.

Populations CY and DC, followed by GZ, LJ and XC contributed most to the diversity component of the total diversity (Fig. 3) due to their high haplotype diversity, while populations GZ, LJ and CY contributed most to the differentiation component of the total diversity. The high contributions of the former populations were due to presence of several rare alleles in each, while in the latter the high contributions were because these populations were composed of a single haplotype absent or rare in other populations. As to the contribution to allelic richness, populations CY and DC contributed most to the differentiation component.

3.3. Species distribution modeling

The accuracy of niche model prediction was high (AUC = 0.956). The highest contribution to niche prediction made variables Isothermality, which is a ration of Mean diurnal temperature range to Temperature annual range (* 100), Annual precipitation and Mean temperature of warmest quarter (0.42, 0.27 and 0.09, respectively). The predicted range of C. tibeticum with suitability ranging from 0 to 1 is shown in Fig. 2.

4. Discussion 4.1. Extent and structure of population genetic variation

Narrow distribution and low abundance make species more prone to genetic drift because small population sizes and their spatial isolation lead to the random fixation of alleles within populations while promoting among-population divergence (Barrett et al., 1991; Hamrick and Godt, 1996; Ross and Travers, 2016). Small range, small population sizes and low density are characteristic features of orchids including Cypripedium species. For example, Qian et al. (2014) studied six populations of Cypripedium japonicum from eastern and central China using inter-simple sequence repeats (ISSR) and found very low within-population genetic variation (ranging from 0.0297 to 0.0587) but high among-population genetic differentiation (G_ST = 0.671). Using ISSR, similar values were obtained for Dendrobium fimbriatum from south China (He = 0.087, G_ST = 0.744) (Ma and Yin, 2009), Cymbidium goeringii from central China (He = 0.194, G_ST = 0.244) (Yao et al., 2007), Gastrodia elata from Hubei province (China) (He = 0.176, G_ST = 0.256) (Wu et al., 2006), Tipularia discolor from eastern United States (He 0-0.069, Φ_ST = 0.415) (Smith et al., 2002), Piperia yadonii from California (He = 0.053-0.071, Φ_ST = 0.40) (George et al., 2009), Platanthera aquilonis, Platanthera dilatata, Platanthera huronensis from eastern and western North America (He = 0.084, 0.131, 0.119; Φ_ST = 0.70, 0.49, 0.36, respectively) (Wallace, 2004). Amplified fragment length polymorphism (AFLP), much less popular marker type than ISSR in studying orchid genetic variation, revealed a similar picture in Liparis loeselii from France and Great Britain (He = 0-0.063, Φ_ST = 0.382) (Pillon et al., 2007), and in Spiranthes romanzoffiana from the British Isles (He not reported, Φ_ST = 0.892) (Forrest et al., 2004).

In congeners of C. tibeticum for which genetic data are available, extent and structure of genetic variation was analyzed using allozymes. In C. candidum, Cypripedium fasciculatum, Cypripedium kentuckiense, Cypripedium macranthos and Cypripedium reginae species variation within populations (He) ranged from 0.030 to 0.185 (Case, 1994; Case et al., 1998; Aagaard et al., 1999; Chung et al., 2009). In Cypripedium calceolus, C. candidum, C. reginae and C. acaule from the northeastern United States and Canada Hs was 0.197, 0.50, 0.024 and 0.080, while G_ST was 0.194, 0.069, 0.164 and 0.349 (Case, 1994). In populations of C. calceolus from Poland He ranged 0.074-0.226 and F_ST was 0.137 (Brzosko et al., 2011). In C. kentuckiense the population average He was 0.042 and F_ST was 0.182 (Case et al., 1998). Very low within-population allozyme variation was also found in North American Cypripedium arietinum and C. acaule (Bornbusch et al., 1994). With a few exceptions, results of the above studies support a general trend in Cypripedium to have reduced genetic variation within populations with a strong difference in allele frequencies even among geographically close populations.

4.2. Causes of population genetic variation in C. tibeticum

The extent and structure of species genetic variation is determined by species life form, breeding system, reproductive characteristics, as well as effects of habitat fragmentation and population isolation (Loveless and Hamrick, 1984; Nybom and Bartish, 2000). We found low within-population variation but high genetic differentiation among populations of C. tibeticum with no clear structuring of this differentiation. Low within-population genetic variation apparently has several causes. One is self-compatibility and another one is low abundance, both reducing chances of a plant to be pollinated by pollen from another individual.

High genetic differentiation among populations of C. tibeticum can result from either limited pollen flow, localized seed dispersal or both. Flowers of C. tibeticum are pollinated by queens of Bombus lepidus, B. lucorum and B. hypnorum (Li et al., 2006), and theoretically, limited pollen movement among but not within populations of C. tibeticum can cause population differentiation in this species, However, the genetic markers used (chloroplast fragments) are maternally inherited in the study species, and can not be affected by the pattern of pollination.

Orchids have tiny seeds with light coat. Due to these features they can be can be easily carried by wind or water and travel significant distances (Summerhayes, 1951; Arditti and Ghani, 2000). However, despite the potential for dispersal over long distances, most of the seeds are dispersed close to the mother plant (Murren and Ellison, 1998; Machon et al., 2003; Jersáková and Malinová, 2007; Jacquemyn et al., 2007; Chung et al., 2009). The spatial genetic structure we found in C. tibeticum agrees with predominantly localized seed dispersal and occasional long-distance dispersal. The geographic isolation created by Hengduan Mountains by running in parallel high mountain ranges separated by deep river valleys appear to play minor, if any, role in C. tibeticum, as populations belonging to different mountain chains possess the same haplotypes (e.g. ZD and JDS, WC and ML).

This study investigated the species genetic variation using cpDNA, a marker type with very limited genomic sampling. Despite this limitation, some important conservation implications stem from the results, as discussed below.

4.3. Conservation implications

Although all wild species of Orchidaceae have been listed in Appendix Ⅱ of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) (UNEP-WCMC, 2014), collecting C. tibeticum in the field is observed occasionally. In recent years, human activities affected natural habitats of C. tibeticum, for example road construction in Demula, Balang and Maan mountains, and increased impact of tourist activities in Ganheba of Lijing County. It has been predicted that at least 20% reduction in number of populations will occur within the species next three generations (Wang and Xie, 2004).

Efficient conservation of C. tibeticum, like any other threatened species, must include several components, complementing and enhancing each other (Volis, 2016a, b, 2018). From all the studied populations, only WC is located within a protected area (Wolong nature reserve). As a majority of C. tibeticum populations are unprotected, collection of seeds in as many populations as possible must be a necessary first step to prevent loss of unique alleles that exist in populations. Lack of geographic structuring of haplotypes implies that unique haplotypes may exist along the species range. The collected seeds should not only be stored in seed gene banks together with their mycorhiza, but also used for quasi in situ living collections. These collections, besides preserving species genetic diversity will serve as reliable sources of seeds for in situ actions. Among the existing in situ actions (enhancement, reintroduction, relocation), creation of new populations in nature reserves and national parks located with the predicted species range appears to be the most relevant for this species option. The results of species distribution modeling suggest that C. tibeticum can find suitable conditions in many existing protected areas (Fig. 2).

Among the studied populations, some appear to have higher priority than others for both seed collecting and protection in situ. The criteria for the selection of priority populations must include both the uniqueness and diversity level of their allelic composition (Petit et al., 1998). Populations CY and DC had the highest contribution to the total gene diversity (CT) as well as allelic richness (CTR). In addition, four rare haplotypes (H1, H2, H3, and H4) occurred in CY and DC populations privately. Therefore, these two populations should have the top priority in conservation planning and implementation.

We thank Zhi-Kun Wu and Zhi-Qiang Zhang for their help in collecting plant materials and two anonymous reviewers for constructive comments. This work was supported by the National Natural Science Foundation of China (grant nos. 31460050 and 31760048 to Yong-Hong Zhang).

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2018.12.002.