1. Introduction

The perennial herb Aristolochia delavayi Franch. (Aristolochiaceae) is an endemic species in the dry-warm valleys along the Jinsha River in southwestern China (Chen et al., 2018; Yu et al., 2020). The plant species has a short single-flower flowering period (~3 days), but the duration of flowering within a population is approximate three months (mostly from June to September). It mainly reproduces in a sexual way, which relies on cross-pollination by pollinators, but the plant itself tends to increase breeding opportunities through tillering (unpublished data). The fruits of A. delavayi are capsules, and its seeds are lacking elaiosomes. A previous study showed that seed dispersal mainly depends on gravity (Chen et al., 2015). Leaves of this plant species are locally used as a spice (Zhou et al., 1995; Sun et al., 2008) and also for medicinal purposes e.g. to strengthen the stomach, to increase the appetite and treat flatulence and malaria (Zhou and Yang, 1995). This plant species is also an important food source for the endangered butterfly species Byasa daemonius Alphéraky (Chen et al., 2015).

Aristolochia delavayi grows mainly in sparse vegetation in river valleys at an altitude of 1220–2250 m (Yang et al., 2014; He et al., 2017). The habitat in this area is fragile with high erosion rate, which makes vegetation restoration difficult (Zhong, 2000; Guan et al., 2013). Due to the habitat destruction by various human activities (e.g. hydropower construction, agricultural activities), and also the long-term overharvesting of this plant species, the populations of A. delavayi decreased gradually, so the species is on the verge of extinction (Yang et al., 2014; Chen et al., 2015). At present, A. delavayi is listed as a wild endangered species (EN) on the red list of higher plants in China and the red list of IUCN (Chen et al., 2015; Qin et al., 2017). It is also treated as a plant species with extremely small populations (PSESP) (Ma et al., 2013). To protect this species and to prevent extinction, protection measures are required.

Generally, before implementation of conservation, studying the genetic structure and diversity of a given species can help to develop scientific conservation strategies (Guan et al., 2013). The precipitous terrain and severe habitat fragmentation in the studied area may hinder gene flow among plant populations. Relevant studies revealed, the endemic taxa in this area generally show high genetic differentiation among populations, but the level of genetic diversity within populations is either high or low. This for example was shown on Cycas panzhihuaensis L. Zhou & S.Y. Yang (Xiao et al., 2020), Musella lasiocarpa (Franch.) C.Y. Wu (Ma et al., 2019), Trailliaedoxa gracilis W.W. Sm. & G. Forrest (Jia et al., 2016) and Munronia delavayi Franch. (Jia et al., 2014). By comparing the genetic differences between A. delavayi and other endemic plant species in this region, it is helpful to clarify the population history, its dynamics, and the effect of external environmental factors on the genetic pattern of the existing populations. Recent genetic analyses on A. delavayi was undertaken by Yang et al. (2014). They analyzed four populations using eight ISSR markers, which showed a high genetic diversity (PPB = 84.71%). However, during our comprehensive field surveys in recent years, we found that the above-mentioned four populations are only located in the upstream part of the distribution area of A. delavayi, whilst the populations in the middle and lower part of the distribution area were not involved. Considering that A. delavayi is vulnerable, it is necessary to use a higher number of molecular markers to re-study and evaluate its population genetics comprehensively.

At present, molecular genetic marker methodologies have been well developed which now allow performing sophisticated analyses. In particular, microsatellite markers have widely been used in the genetic research of wild natural populations and endangered species (Balloux and Lugon-Moulin, 2002; Chen et al., 2009; Tang et al., 2008; Wang et al., 2006; Yang et al., 2018). This kind of markers are rich in polymorphism, good in repeatability, mostly co-dominant, widely distributed in the genome, and reveal a large amount of information even within a small sample size (Jarne and Lagoda, 1996; Zhang and Hou, 2004; Habel et al., 2010; Zhang et al., 2019). In this study, SSR primers for A. delavayi were developed for the first time by using genome skimming technology. Based on these primers, we were able to analyze the genetic structure and its diversity of existing wild populations of this species. The obtained results allowed us to analyze the internal causes of the species' endangerment and propose reasonable conservation strategies.

2. Materials and methods 2.1. Plant materials

During the population and habitat surveys of the past few years, molecular samples of 193 individuals from ten populations of Aristolochia delavayi were collected along the Jinsha River Basin in July 2018 and August 2019. Detailed information about the plant material is provided in Table 1. Most of the samples were collected in Yunnan Province, except the samples belonging to the population ML, which was collected from Sichuan province (Fig. 1). From each population were more than 20 individuals collected, with a few exceptions, such as DD (only two individuals) and WB (seven individuals) due to the rare occurrence of this plant species. Fresh and healthy leaves were dried with silica gel immediately after collection and stored at −20 ℃ prior DNA extraction. Voucher specimens of partial populations were deposited in the Herbarium of Kunming Institute of Botany (KUN), Chinese Academy of Sciences.

2.2. Laboratory protocols

Genomic DNA was extracted from the dried leaves using a modified CTAB method (Doyle and Doyle, 1987), and the quality and concentration of DNA were measured by using a spectrophotometer (Nanodrop - 1000). Using the genome skimming technology, the genomic DNA from two individuals in the DJ and LQ populations were used for genome skimming. The specific process is as follows: The sequencing was performed on the Miseq Benchtop Sequencer (Illumina X-Ten) using the 2 × 250 bp read mode. The obtained data were assembled using software SPAdes and the microsatellite sequences of two samples were compared by using the software QDD 2.1 Beta (Meglecz et al., 2010). MISA software was used to design SSR primers based on the flanking sequence of the SSR loci.

The overall PCR reaction system volume was 20 μL, including 1 μL DNA template of 50–60 ng/μL; 2 ×Taq PCR MasterMix 10 μL, 0.3 μL of each primer at a concentration of 10 μmol/L; and 8.4 μL ddH₂O. The amplification procedure was set as follows: Pre denaturation at 95 ℃ for 3 min; denaturation at 95 ℃ for 30 s, annealing at appropriate temperature for 30 s, extension at 72 ℃ for 30 s, a total of 32 cycles; and the final extension step at 72 ℃ for 5 min. The annealing temperature (Tm) was determined according to the reference value of primer synthesis and optimized in the experiment. The PCR reaction was performed on the DNA Thermal Cycler (Applied Biosystems).

For the primer screening, the genomic DNA of four individuals from four populations were selected for PCR amplification. The PCR products were detected using 8% non-denaturing polyacrylamide gel electrophoresis. The gel was stained with silver and developed after electrophoresis, and the amplification was observed. The amplification effects of primers were compared, and polymorphic primers were screened for genetic analysis of all samples. The PCR products were separated and visualized using QIAxcel of capillary gel electrophoresis system (ABI PRISM 3730 XL, USA). The final banding data were read by GeneMarker v.2.2.

2.3. Data analysis

To detect linkage disequilibrium (LD) between the loci, we used the program GENEPOP v.4.7 (Rousset, 2008). The following parameters of genetic diversity at population and species level were calculated by the GenAlEx v.6.41 (Peakall and Smouse, 2006): Population sample sizes (N), number of alleles (Na), effective number of alleles (Ne), Shannon's information index (I), observed (Ho) and expected (He) frequency of heterozygotes, percentage of polymorphic loci (PPL), and fixation index (F). The GenAIEx program was also used to calculate the inbreeding coefficient (F_IS) and gene flow (Nm) at population level.

The geographic distances between populations were calculated using the Franson CoordTrans Program v.2.3 Meanwhile, POPGENE v.1.32 was used to calculate the Nei's genetic distance between populations of Aristolochia delavayi. We established a geographic distance as well as a genetic distance matrix, and used Mantel test in GenAlEx v.6.41 to detect correlation between geographic and genetic distance (IBD) (Mantel, 1967; Smouse et al., 1986; Peakall and Smouse, 2006). The recent bottleneck effects were detected by using BOTTLENECK v.1.2.2. For this purpose, two evolutionary models, infinite alleles model (IAM) and two-phase model (TPM) (Wang et al., 2010), were selected, and 1000 iterations statistics were carried out by performing the Sign test and Wilcoxon sign-rank test, respectively.

Concerning the genetic structure, hierarchical analyses of molecular variance (AMOVA) were performed to assess the genetic structure within and between populations by using Arlequin v.3.11 (Excoffier et al., 2005). Principal coordinate analysis (PCoA) was conducted using GenAlEx v.6.41 to test the genetic similarity of all the individuals in the included population (Peakall and Smouse, 2006). UPGMA cluster analysis based on the Nei's genetic distance between populations was performed by using MEGA v.7.0.14. Bayesian clustering method in Structure v.2.3.4 was used to cluster all individuals individually and determine subsequently the number of genetic groups (Pritchard et al., 2000). The parameters for the calculation of the ΔK value and analyze possible genetic structure were set as follows: The K value was set to 10, 10 replicates, using the admixture model and assumed allele frequencies correlated, followed by a burn-in of 1 × 10⁶ iterations (Evanno et al., 2005).

3. Results 3.1. Development and screening of SSR primers

In this study, 3689 pairs of SSR primers were successfully discovered and designed based on the genome skimming technique, and 100 pairs of primers were randomly selected for polymorphism primer screening. Finally, 15 pairs of SSR primers with good amplification effects and high polymorphism were selected and used for genetic analysis. The details of these primers are shown in Table 2.

3.2. Genetic diversity of the Aristolochia delavayi

The GENEPOP test showed that only a few detectable loci were linked in a few populations, but none of them were linked with other loci in more than three populations. Due to this result, the above-mentioned 15 pairs of SSR primers were used for the genetic analysis of A. delavayi (Wang et al., 2013) and the results of these analyses are presented in Table 3, Table 4. Overall, a total of 78.50 alleles (Na) and 44.45 effective alleles (Ne) were detectable. Within these 15 pairs of SSR primers, 3–10 alleles and 2–6 effective alleles with an average of 5.23 alleles and 2.96 effective alleles could be detected (Table 3). The calculated observed heterozygosity (Ho) varied from 0.333 (DD) to 0.705 (SB), whilst the expected heterozygosity (He) varied from 0.233 (DD) to 0.677 (XZ). The calculated Shannon information index (I) ranges from 0.335 (DD) to 1.485 (SB). At population level, except the populations DD and WB, the percentages of polymorphic loci (PPL) in all other populations was 100% (Table 4). The calculated inbreeding coefficient (F_IS) value varied from −0.336 to 0.319, with an average value of 0.056. The average values of the genetic differentiation coefficient (F_ST) and the gene flow (Nm) among populations are 0.328 and 0.591, respectively.

In the performed Bottleneck test, under the IAM hypothesis, both methods exhibited, that the population YS deviated from the mutation drift equilibrium, showing a significant heterozygosity excess. Under the TPM hypothesis, the sign test revealed that the population LQ deviated from the mutation drift equilibrium, exhibiting a significant heterozygosity excess (Table 5). This test indicated, that the above-mentioned two populations have recently experienced a bottleneck effect.

3.3. Genetic differentiation and genetic structure

The AMOVA results (Table 6) indicated significant differences among populations of Aristolochia delavayi, because 31.62% of the variation existed among populations and 68.38% occurred within populations. The genetic differentiation coefficient (F_ST) at species level is 0.316, which is considered as significant, i.e. 0.25 < F_ST < 1. This indicates a significant genetic differentiation among populations (Wright, 1978). The performed Mantel test resulted in a significant positive correlation between the genetic distances of the studied populations and its geographic distance (P = 0.01 < 0.05 and r = 0.561 > 0.5) (Fig. 2), which also indicates the presence of IBD in the population structure. The result of principal coordinate analysis indicated, that the analyzed ten natural populations can be roughly grouped into three groups (Group A includes populations HTX, SB, DJ, XZ, ML, DD, Group B includes YS, WB, LQ and Group C includes HP) (Fig. 3). This result was supported by the performed UPGMA cluster analysis based on genetic distance (Fig. 4). Structural analysis showed, a ΔK is at the maximum, when K = 6. This indicated that all populations can be assembled in six groups. These results further refine the grouping results of the above-mentioned methods. Among all analyzed samples, six populations collected in Northwestern Yunnan consist of three groups; the three populations in Central Yunnan clustered in two groups, and the Huangping population was still remained separately (Fig. 5).

4. Discussion 4.1. Genetic diversity of Aristolochia delavayi

Genetic variation of species is the premise of local adaptation and evolution. It is also considered as an important parameter to determine the priority of population conservation in the protection of endangered plants (Schaal et al., 1998; Laikre et al., 2010; Zhao and Gong, 2015). So, understanding the genetic status of Aristolochia delavayi provides a rational basis for the evaluation of conservation work and formulating of effective protection measures. Young et al. (1996) showed that genetic variation of population eroded with reduced remnant population size. However, our analysis results showed that the expected and observed heterozygosity of A. delavayi at species level are 0.550 and 0.553, respectively. These values are higher than reported values from other endemic or endangered taxa in this region, such as Nouelia insignis Franch. (He = 0.149, Ho = 0.216) (Luan et al., 2006), Buddleja crispa Benth. (H = 0.314, I = 0.485) (Zhang et al., 2015), Cycas hongheensis S.Y. Yang & S.L. Yang ex D.Y. Wang (He = 0.435, Ho = 0.403) (Zhao and Gong, 2015), T. gracilis W.W. Smith et Forrest (Hs = 0.489) (Jia et al., 2016), Amorphophallus albus P.Y. Liu & J.F. Chen (He = 0.504, Ho = 0.528) (Tang et al., 2020), C. panzhihuaensis L. Zhou & S.Y. Yang (He = 0.328, Ho = 0.189) (Xiao et al., 2020). Our calculated result is consistent with the result reported by Yang et al. (2014), indicating that the species possess relatively high genetic diversity. Moreover, this indicates that the recent population size reduction only affects little the genetic variation of the species. According to Loveless and Hamrick (1984) and Nybom and Bartish (2000), the factors which affect the genetic variation of species generally include population history, population size, reproduction pattern, breeding system, genetic drift, gene flow, natural selection, geographical distribution, etc. The reduction of the population size cause genetic drift effects, but drift has only a little effect on genetic variation in limited generations (Young et al., 1996). The recent habitat fragmentation event may be below the fragmentation threshold, which probably not cause the loss of genetic variation (Prober and Brown, 1994). We assume that in the past few decades, with the rapid reduction of the population size, the genetic drift effects are not yet accumulated, and the species may retain its genetic variation by its sexual as well as asexual reproduction.

4.2. Genetic differentiation and structure

The results of AMOVA analysis showed that the genetic variation of Aristolochia delavayi occurred mainly within the population, but there was a significant genetic differentiation among the populations observable. The high percentage of genetic diversity within the population may be due to the retention of genetic resources prior the population reduction and to the outcrossing of sexual reproduction within the population. For genetic differentiation among populations, it is generally considered to be caused by restricted gene flow (Loveless and Hamrick, 1984; Slatkin, 1985). The results of this study showed, that the level of gene flow among the population is generally low (Nm = 0.591 < 1). This reduced gene flow may be related to the limited reproductive characteristics of the plant itself.

During the performed field work, it was observed that the sexual propagation strategy of Aristolochia delavayi relies on the participation of pollinators. The main pollinators are some small flies belonging to the families Ceratopogonidae and Chironomidae, which basically have weak flying abilities and low pollination efficiency. The seed dispersal only depends on gravity with a limited dispersal distance, and the seeds are prone to dry under this climatic circumstance. In addition, the active geological structure and erosion in the area caused major rivers in the region which formed valleys with a steep and rugged topographies. This became an obstacle for gene exchange among most plant species in the region as reported previously (Yue et al., 2012; Zhang et al., 2015; Tang et al., 2020). At the same time, the fragmentation of habitat caused by human activities became another obstacle for the gene flow of the studied plant species. Presumably, the restricted pollen and seed dispersal mechanisms of A. delavayi and its relative geographical isolation may be the main reasons for the observed genetic differentiation among populations of this species. Though the limited number of populations, the present outcrossing reproductive strategy may enable sufficient gene exchange among individuals, which maintain a high genetic diversity within the population.

The performed Mantel test detected isolation by distance (IBD) in the population structure of the studied species, which indicated that geographic barriers played an important role in the formation of the present genetic structure of populations, but there may be potential biases in further cluster analysis (Perez et al., 2018). Therefore, we have performed principal coordinate analysis. In this study, the results of principal coordinate analysis and UPGMA clustering are consistent, showing a clear geographic regionality. The ten populations of Aristolochia delavayi are assembled in three groups, which are clustered in Northwest Yunnan, Central Yunnan, and Huangping areas, respectively. Bayes cluster analysis further refined the above grouping results. The genetic clustering results of the three methods all disclosed the separation of the Huangping population. This indicates a distant genetic relationship with the other studied populations. The geographical distribution of this population is rather isolated, far away from the Jinsha River and its tributaries. Compared to other populations, there were no obvious differences in the main taxonomic characteristics (flower and leaf) of A. delavayi in Huangping population. But the different plant size (Fig. 6), leaf scent (Yu et al., 2020), and diverged genetic data suggest that this population could be a subspecies or variety of A. delavayi. But further efforts are required to validate this.

4.3. Conservation suggestions

Aristolochia delavayi is an important species in the traditional medicinal and as spice, and its essential oil is essential for the development and application of antibacterial drugs (Li et al., 2013). At the same time, this species is an important host of the endangered butterfly Byasa daemonius. To sustain this plant species may effectively promote the gradual recovery of the B. daemonius, and may also avoid cascading effects (Koh, 2004; Brodie et al., 2014; Chen et al., 2015). Over the past few decades, the habitat destruction and overharvesting by human activities have led to a rapid decline in the population sizes. At present, only ten isolated populations of A. delavayi were found and most of them comprise of only few individuals. Therefore, the wild resource of the plant species needs effective protection urgently.

Because of the rarity of wild populations of Aristolochia delavayi and its significant genetic differentiation among the populations, it is necessary to protect the present populations. In particular, the two populations (YS and LQ), which have recently experienced genetic bottlenecks, would need special attention. For example, intensify the implementation of in situ protection may maintain the existing population size, and improve their population adaptability.

In view of the uniqueness of Huangping population from other populations, it is recommended to perform in situ conservation measures to sustain the current genetic resources for future research and utilization. For the Dadong population, due to the limited number of individuals found within this population, it is recommended to conduct first a thorough survey in order to assess the size of the existing population. Afterwards propagation of these plants, ex situ conservation should be implemented. Additionally, re-introduction of these cultivated individuals should be performed in order to increase the number of individuals within this population. If ex situ conservation measures are implemented, the pollination and seed dispersal processes should be taken into account, otherwise limited pollinators and seed dispersers may impact the reproduction systems negatively (Tang et al., 2019).

According to the genetic differences of the populations of Aristolochia delavayi in different regions, it is recommended to divide them into three protection units. Due to the low gene flow among populations, measures such as artificial pollination and seed dispersal may be taken into account. This would ensure and, as a consequence, increase the gene flow between populations. These measures should be carried out within the same protection unit in order to avoid outbreeding depression (Tallmon et al., 2004; Barmenttlo et al., 2018).

Remarkably, a large number of Aristolochia zhongdianensis J.S. Ma individuals are mixed with A. delavayi in the population SB. Our morphological data indicated that some suspected hybrids occurred in this population (unpublished data). Previous studies exhibited, that in two species showing hybridization phenomenon, gene introgression tends to dilute genetic variation within a small population (Levin et al., 1996; Wolf et al., 2001). Therefore, further research is needed to reveal possible gene introgression between A. delavayi and A. zhongdianensis in the population SB. This would clarify the direction of gene flow, and the trends of genetic variation of the population. Such studies may also provide a scientific guidance for reasonable and successful protection.

Author contributions

Yu-Long Yu and Gao Chen collected plant materials, and Zhi-Xiang Yu assisted in the field population survey. Yu-Long Yu independently performed the experiments, analyzed the data and wrote the manuscript. Hui-Chun Wang and Johann Schinnerl provided assistance in the language modification and polish of the manuscript. Rong Tang provided useful suggestions for the writing of the article. Gao Chen and Yu-Peng Geng designed and supervised the study, and also revised the manuscript. All authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that there are no conflicts of interest.

Acknowledgments

We thank Chencan Liao, Maoyao Yan and Qi Yu for their assistance in collecting samples. We appreciate Zhi Chen for picture optimization. Support for this study was provided by grants from the NSFC (National Natural Science Foundation of China)-Yunnan Joint fund to support key projects (No. U1602264) and Yunnan Ten Thousand Talents Plan Young & Elite Talent Project to G. Chen (YNWR-QNBJ-2018-017); and NSFC (National Natural Science Foundation of China) (No. 31660057) to Y.P. Geng; Science & Technology Basic Resources Investigation Program of China for Survey and Germplasm Conservation of PSESP in Southwest China (2017FY100100) to W.B. Sun.