1. Introduction

The genetic diversity of a species, which is the product of long-term evolution, is the basis of species survival, adaptation, and development (Meffe and Carroll, 1997). Accordingly, protecting the genetic diversity of a species ensures that species maintain the ability to adapt to change (Huang et al., 2016). Thus, it is critical for researchers to understand the genetic structure and evolutionary history of endangered species in order to provide guidance for ex situ conservation, introduction, and cultivation (Groom et al., 2006).

Rhododendron L. (Ericaceae), which consists of over 1000 species around the world, is the largest genus of angiosperms in China (Goetsch et al., 2005). Rhododendron meddianum Forrest is a perennial evergreen plant of Rhododendron, subg. Hymenanthes, sect. Ponticum, subsect. Thomsonia, distributed in southwestern China and northeastern Myanmar (Fig. 1; Fang et al., 2005). R. meddianum has high ornamental value, but is listed as critically endangered (Qin et al., 2017). Recently, our field surveys discovered three wild populations of R. meddianum with a narrow area of distribution in western Yunnan, China (Fig. 1a). In this narrow region, R. meddianum habitats may be severely disturbed by anthropogenic activities (e.g., wind power stations and roads) (Fig. 1b). This is consistent with studies that indicate human-induced environmental change and habitat fragmentation threaten about 60% of wild Rhododendrons (Gibbs et al., 2011). Thus, it is urgent to carry out studies on the conservation biology of R. meddianum. Importantly, understanding the population genetics (e.g., genetic diversity, population structure) of these populations will help researchers set appropriate conservation strategies.

Population genetics of plants have long been assessed using molecular markers, including simple sequence repeat (SSR), amplified fragment length polymorphism (AFLP), random amplified polymorphic DNA (RAPD), and expressed sequence tag SSR (Yilmaz, 2016; Stone et al., 2019; Everaert et al., 2020; Wagutu et al., 2020). Recently, single nucleotide polymorphism (SNPs) based on next-generation sequencing (NGS) have become the markers of choice for determining population structure because they are abundant, stable in the genome, and can be accurately scored (Deschamps et al., 2012). One sample, low cost, high throughput approach to discovering hundreds to thousands of SNPs in non-model species is restriction-site associated DNA sequencing (RAD-seq). RAD-seq methods have been used on plant communities to reveal their evolutionary history and local adaptability (Jiang et al., 2019), identify the genetic structure of populations (Yamashita et al., 2019), and determine how environmental and topographical elements affect genetic diversity and differentiation (Xiao et al., 2020). RAD-seq methods have also been used to identify genetic lineages of plants, which has rendered references for the conservation of rare species (Amor et al., 2020).

In this study, neutral loci generated by double digest restriction-site-associated DNA sequencing (ddRAD) were used to determine the genetic diversity, population structure and demographic history of Rhododendron meddianum in southwestern China. In addition, we elucidate the evolutionary history between R. meddianum and its close relative Rhododendron cyanocarpum (Franch.) W.W. Smith (a threatened species endemic to Dali Cangshan Mountain, Yunnan province). The findings of this study will contribute to a better understanding of the evolutionary history of R. meddianum and provide a scientific basis for its conservation and management.

2. Materials and methods 2.1. Experimental material

Experimental material was taken from 45 samples of three wild Rhododendron meddianum populations from Yunnan province, southwestern China (Table S1, Fig. 1). The distance between sampling plants within each population was greater than 50 m. After collection, the young, tender leaves were quickly dried and preserved with silica gel.

2.2. DNA extraction and ddRAD library construction

Genomic DNA was extracted for each sample of Rhododendron meddianum by the CTAB method (Doyle and Doyle, 1987). DNA quality was assessed visually on 1.2% agarose gels and subsequently quantified using a Qubit 3.0 Fluorometer. Library construction was constructed following a previously reported ddRAD protocol (Peterson et al., 2012). Briefly, 100 ng of genomic DNA of each sample was digested with two restriction enzymes, EcoRI and MseI (New England Biolabs, Beverly, USA) at 37 ℃ for 5 h. Restriction enzymes were then inactivated by heating at 65 ℃ for 20 min. After ligation with individually barcoded EcoRI and MseI adapters with T4 DNA ligase (New England Biolabs, Beverly, USA) for each sample at 16 ℃ for 4 h, the reaction was stopped by heating at 65 ℃ for 20 min, and finally the product was maintained at 12 ℃. The fragments (300–500 bp) were recovered using a gel extraction kit (Omega Inc., USA). The purified product was amplified by PCR to the expected concentration, and then 150 bp paired-end sequencing was performed on the Illumina HiSeq X-ten platform. Sequencing was completed by Guangzhou Jierui Biotechnology Co., Ltd. Each sample generated about 0.5 Gb of data. The sequencing reads of R. meddianum in this study have been deposited in the NCBI sequence read archive (SRA) under bioproject accession PRJNA706564.

2.3. SNP mining

Raw sequence data was analyzed using Stacks v.2.0 (Catchen et al., 2011, 2013). First, the module process_radtags was used to demultiplex and sort the raw data based on the barcodes used in each sample; low-quality reads with missing restriction sites or containing ambiguous barcodes were eliminated. Sequence quality was verified by Fastqc (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The ustacks module was used to cluster single sample sequences to form loci, with parameters were set as follows: minimum depth of coverage required to create a stack (-m), 2; maximum distance (in nucleotides) allowed between stacks (-M), 2; maximum distance allowed to align secondary reads to primary stacks, 4; max number of stacks allowed per de novo locus, 3; minimum alignment length, 0.8; alpha significance level for model, 0.01. Loci data of all of samples were merged into a catalog using cstacks module with default parameters. The sstacks module was run to compare the single sample loci to the catalog so as to obtain the SNPs, as well as allele and tag information of the sample. Finally, thepopulations module was used to screen SNPs, with parameters set as follows: minimum number of populations a locus must be present in to process a locus (-r), 6; minimum percentage of individuals in a population required to process a locus for that population (-p), 0.75; minimum minor allele frequency required to process a nucleotide site at a locus, 0.01; maximum observed heterozygosity required to process a nucleotide site at a locus, 0.6. One random SNP per locus was kept (-write_random_snp) to reduce linkage disequilibrium (Falush et al., 2003; Kaeuffer et al., 2007). The output result was exported in VCF format using the populations module.

2.4. Data analysis 2.4.1. Analysis of genetic diversity and structure

Tajima'D values were estimated to detect the selection effect (Tajima, 1989). We calculate Tajima'D using VCFtools v.0.1.15 (Danecek et al., 2011); the sliding window size of the operation was 3000 bp. The populations module (Catchen et al., 2013) of Stacks software was utilized to calculate genetic diversity parameters, including nucleotide diversity (π), the number of private alleles (Private), expected heterozygosity (H_E), observed heterozygosity (H_O) and inbreeding coefficient (F_IS). Outlier loci (False Discovery Rate (q) > 0.05) was detected by BayeScan v.2.1 (Foll and Gaggiotti, 2008) with default parameters to guarantee that neutral SNPs were generated for subsequent analysis. The format conversion was implemented using PGDSpider v.2.1.1.5 (Lischer and Excoffier, 2012).

Arlequin v.3.5.2 was used to estimate pairwise differentiation coefficient (F_ST) between populations (1000 permutations) (Excoffier et al., 2007) and to calculate the genetic variation between and within populations with analysis of molecular variance (AMOVA). STRUCTURE v.2.3.5 (Pritchard et al., 2000) was used to analyze the genetic relationship between samples via speculating the proportion of a certain hypothetical ancestor of the test samples, with the number of genetic clusters (K) set to 1–8. Each K was run with 10 independent iterations, each with a burn-in period of 20, 000 iterations and 100, 000 Markov chain Monte Carlo (MCMC) iterations by assuming the admixture model. The output result of LnP (D) under different K values was subsequently visualized by Structure Harvester v.0.6.94 (http://taylor0.biology.ucla.edu/struct_harvest/). Optimal K values were determined by maximum delta K and the mean LnP (K) value (Evanno et al., 2005). Principal components analysis (PCA) was performed via Plink v.1.90 (Purcell et al., 2007). The top 10 principal components (PCs) were extracted with 'pca' parameters. The top two PCs (PC1, PC2) were visualized in R v.3.3.1 (R Core Team, 2016) to confirm the subsection structure of all 45 samples. MEGA X (Kumar et al., 2018) software was used to construct Neighbor Joining (NJ) trees, and finally, EvolView online tool (http://www.evolgenius.info/evolview) was edited for visualization.

2.4.2. Historical demographic inference of Rhododendron meddianum

Stairway Plot v.2.1 was used to infer recent historical dynamics of effective population size of Rhododendron meddianum based on nuclear SNP data (Liu and Fu, 2015). One-dimensional unfolded SFS (1D-SFS) was generated from posterior probabilities of sample allele frequencies using the python script easySFS (https://github.com/isaacovercast/easySFS). To estimate the trend of effective population size (Ne) over time, 67% of all sites was randomly selected for training. The number of random breakpoints (nrand) in each test was set at 29–90. These break points may be thought of as points at which Ne changes in time. Using the median value of these estimates as final values, Morimoto et al. (2003) observed that Rhododendron reticulatum D. Don ex G. Don and Rhododendron macrosepalum Maxim. began to flower at about 10 years old. The Rhododendron ponticum L. began to flower at about 10 years (Cross, 1975). Yoichi et al. (2016) gave the mutation rate for Rhododendron weyrichii Maxim. as 1.581e-9. We therefore set the generation time and mutation rate to 10 years and 1.581e-9 for R. meddianum.

2.4.3. Evolutionary history of Rhododendron meddianum and R. cyanocarpum

Because we used ddRAD sequencing in both species (Liu et al., 2020a), we assumed SNPs were consistent with the 2.3 section analysis and estimated population structure the same as 2.4.1 analysis (BayeScan, STRUCTURE, PCA and NJ trees analysis). Raw sequence data of R. cyanocarpum was obtained from SRA in the NCBI (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA640883).

The results of our analysis of the population genetic structure (see the results section) of both Rhododendron meddianum and R. cyanocarpum indicated that these species are clearly divided into two genetic groups (Caojian (CJ), Pianma (PM) and Houqiao (HQ) for R. meddianum group; Huadianba (HDB), Xiaohuadianba (XHD), Dianshita (DST) and Dahuayuan (DHY) for R. cyanocarpum group, Fig. S1). To infer the divergence time and migration events of these two groups, we used Fastsimcoal2 v.2.6 (Excoffier et al., 2013) to test four demography models according to 1D-SFS (Model 1–4, Fig. 4). The mutation rate and generation time were set as 1.581e-9 and 10 years. Each model was run 50 times with 100, 000 coalescent simulations and 40 optimization cycles. Model comparison was based on the maximum likelihood value (Max log (Est_hood), Table 4) using the Akaike Information Criteria (AIC) and AIC weight of evidence. The maximum AIC weight value was selected for the optimal model (Excoffier et al., 2013). Finally, Parameter-confidence intervals (CIs) of the optimal model were calculated from 100 para-metric bootstrap replicates by simulating SFS from the maximum composite likelihood estimates (_maxL.par file), with 50 independent runs in each bootstrap.

3. Results 3.1. ddRAD sequencing and data processing

ddRAD sequencing was performed on 45 individuals of Rhododendron meddianum, a total of 268, 549, 498 raw reads were obtained (sample reads range: 1, 221, 942–1, 484, 088), and 265, 045, 933 reads were retained (per sample, 1, 210, 660–14, 707, 008) after quality control. A total of 3, 503, 562 reads were eliminated due to low quality and no rad tag. The average GC content per sample was 41.56% (40%–47%). Using ustacks module, 7, 252, 426 loci were retained. The mean coverage depth of locus of sample were different, which might be due to the degradation of sample DNA (Table S2). A total of 4, 204, 450 catalogs were obtained after merging sample loci. After mapping each sample loci to the generated catalogs, 3377 loci with 482, 194 variant loci were gained. After retaining one SNP per locus, 3377 SNPs were finally obtained for further analysis.

3.2. Genetic diversity and structure of Rhododendron meddianum

Tajima'D statistics were significantly negative for the three populations of Rhododendron meddianum (CJ, PM and HQ), with values ranging from −0.4151 to −0.5057 (Table 1), supporting the hypothesis that R. meddianum populations underwent a population expansion. Outlier loci with a significant elevated ofF_ST values (q < 0.05) were not found via BayeScan analysis, indicating that ddRAD-seq loci were not under strong selection (Fig. S2a).

The genetic diversity of 3377 neutral loci is shown in Table 1. The CJ population showed the highest genetic diversity (H_E = 0.0781 ± 0.0022, π = 0.0813 ± 0.0023). In contrast, the HQ population showed the lowest genetic diversity (H_E = 0.0646 ± 0.0021, π = 0.0671 ± 0.0022). The PM population had the largest number of private alleles (952), followed by the CJ population (946), and the HQ population had the fewest private alleles (902). The inbreeding coefficient (F_IS) of each population were positive (0.0620 ± 0.0172–0.1428 ± 0.0130), indicating a deficit of heterozygotes (Wright, 1922).

All the pairwise F_ST values between the populations were statistically significant (P < 0.01, Table 2). Moderate degrees of differentiation were found between the three populations (0.1582–0.2388). AMOVA indicated that 19.88% (P < 0.01) of the genetic variation occurred among populations, and 80.12% (P < 0.01) of the variation occurred within populations (Table 3). STRUCTURE analysis supported the best genetic cluster (K = 3) based on both Delta K and mean LnP (K) (Fig. 2a and b). In the case of K = 3, three wild populations of R. meddianum were clearly distinguished, forming three clusters. When it was assumed K = 2, the CJ and PM populations were grouped into the same cluster, indicating that they share a close genetic composition. The results of K = 4 was consistent with K = 3, i.e., the populations were clearly distinguished in accord with geographical distribution (Fig. 2c). Consistent with STRUCTURE results, both the PCA and NJ trees of R. meddianum samples revealed three major clusters corresponding to the CJ, PM and HQ populations (Fig. 2d and e).

3.3. Population demographic history

To reconstruct the demographic history of Rhododendron meddianum, we used a stairway plot to assess population size changes and obtained a demographic history from 150 to 100, 000 years ago (Fig. 3). The results revealed that a population bottleneck occurred about 70, 000 years ago, followed by a rapid expansion 60, 000–70, 000 years ago, and then the populations maintained a stable Ne from 150 to 60, 000 years ago, including the Last Glacial Maximum (LGM: 21, 000 years ago). However, the Ne declined 300 years ago (Fig. 3).

3.4. Phylogeographic patterns of Rhododendron meddianum and R. cyanocarpum

A total of 1495 loci with 207, 927 variant sites were retained after sequence data processing and quality control. After retaining one SNP per locus, 1495 SNPs were ultimately selected. Our BayeScan analysis did not identify any outlier loci with significantly elevated F_ST values (q < 0.05), which indicates that 1495 loci of the two species were neutral (Fig. S2b).

STRUCTURE analysis showed that the optimal K value was K = 2 with Delta K of 235.302 (Figs. S3a, b, c), which suggests that the populations of R. meddianum and R. cyanocarpum are grouped into different clusters. This result is accordance with the conclusions of PCA and NJ trees analysis (Figs. S3d and e).

3.5. Evolutionary history of Rhododendron meddianum and R. cyanocarpum

Model 4 (AIC weight = 0.9999; Table 4) estimated that the divergence time between Rhododendron meddianum and R. cyanocarpum from a common ancestor occurred in the late Pliocene, ~3, 048, 100 years ago (TDIV, 95%CI: 2, 213, 628–5, 026, 950) (Table 5, Fig. S4). The Ne of the two species ranged from 33, 827 to 35, 113, which is higher than that of the common ancestor (25, 546, 95%CI: 15, 435–37, 242). The gene flow/migration (Mmc) between the two species was extremely small (1.2477e-5, 95%CI: 1.1301e-5–1.7042e-5).

4. Discussion 4.1. Genetic diversity

Compared to widespread plant species, the genetic diversity of narrowly distributed species is likely to be lower because of inbreeding depression and genetic drift (Godt et al., 1997; Setoguchi et al., 2011). Unexpectedly, ddRAD sequencing indicated that nucleotide genetic diversity (0.0671–0.0813) of Rhododendron meddianum is higher than that of other endangered plant species, such as Viola uliginosa Muhl. (0.013–0.023), Cymbidium aloifolium (Linn.) Sw. (0.0136), and Geodorum densiflorum (Lam.) Schltr. (0.0359) (Roy et al., 2016; Lee et al., 2020). In addition, R. meddianum has higher genetic diversity than other widespread species, e.g. Taxus wallichiana Zucc. (0.0281–0.0433) (Qin, 2017). Yet, a similar level of genetic diversity was detected in R. cyanocarpum (0.0702) (Liu et al., 2020a). These results generally support the view that some rare and endangered species can maintain high-level genetic diversity even at small population sizes (Li et al., 2015; Stone et al., 2019; Liu et al., 2020a).

The genetic diversity of species is generally influenced by their distribution ranges, life span, breeding system, seed dispersal mechanism, and evolutionary history. In addition, genetic diversity is usually higher in outcrossing species than in selfing species (Hamrick, 1990; Nybom, 2004). Previous studies found that Rhododendron have adhesive pollen and require pollinators, indicating that Rhododendron may primarily rely on outcrossing (Ng and Corlett, 2000; Ma et al., 2015; Huang et al., 2017; Li et al., 2018). Our field surveys indicate that the primary pollination mode of R. meddianum is cross-pollination by insect (mainly Bombus spp.) and birds. The high genetic diversity of the three remaining populations of R. meddianum might also be derived from their ancestral populations, which has been the case in other endangered plants, e.g. Tricyrtis ishiiana (Kitag. et T. Koyama) Ohwi et Okuyama and Rhododendron protistum var. gigantum (Forrest ex Tagg) D.F. Chamb (Setoguchi et al., 2011; Wu et al., 2015).

4.2. Genetic differentiation and genetic structure

STRUCTURE analyses, PCA and NJ trees demonstrated that the three populations of Rhododendron meddianum (CJ, PM and HQ) are genetically separated in accordance with their geographical distributions. The F_ST of R. meddianum (0.1582–0.2388) revealed that a moderate genetic differentiation among the populations occurred (Wright, 1922). Furthermore, AMOVA showed that 19.88% of the genetic variation occurred amongR. meddianum populations. For rare and endangered plants, population genetic differentiation is usually affected by large geographical distances between sampled populations (Hamrick, 1990; Nybom, 2004). The distance among the populations in this study was more than 50 km, which hardly allows gene flow (pollen or seeds) between populations. Previous studies have indicated that the seed dispersal distance of Rhododendron species ranges from approximately 30–80 m (Ng and Corlett, 2000), and its pollen can be transmitted by bees and birds with dispersal distance of 3–10 km of (Ng and Corlett, 2000; Huang et al., 2017; Li et al., 2018). In addition, in line with our observation, seedlings grow sporadically around the parent tree. Hence, the differentiation among the populations could be attributed to distance-limited pollen flow and short-distance seed dispersal.

4.3. Population dynamic history

We used stairway plot to explore the demographic history of Rhododendron meddianum between 150 and 100, 000 years ago (Fig. 3). Our analyses showed that a population bottleneck occurred approximately 70, 000 years ago. This bottleneck event, which is consistent with that of other plant species (e.g. Cycas balansae Warb., see Feng et al., 2014; Rhododendron rex Lévl., see Zhang et al., 2020), may have been caused by climate oscillations. This result also indicates a population expansion around 0.06–0.07 Ma, which is supported by the significantly negative value of Tajima'D (−0.6740) (Tajima, 1989). However, the Ne of R. meddianum experienced a contraction 300 years ago, which may have been the result of human activity and climate change. According to our survey, the CJ population habitat was severely damaged and fragmented by the construction of wind power stations and roads. Moreover, the habitat of the other two populations are all located near farmlands and threatened by grazing. Endangered plants usually experience decreases in its genetic diversity after bottleneck events (Wu et al., 2017). In our study, genetic diversity of R. meddianum was higher than that of other species, suggesting that a historical bottleneck did not reduce its genetic diversity, similar to in other Rhododendron species, such as R. cyanocarpum, R. protistum and R. rex (Wu et al., 2017; Liu et al., 2020a; Zhang et al., 2020). Adaptive radiations often occur when a species expands its ecological niche to adapt to diverse habitats (Cronk, 1995). Thus, one explanation for the high-level genetic diversity of R. meddianum is that when climate oscillations led to a population bottleneck, this event may have initially reduced the genetic diversity of the population, but then generated new genetic variation. In addition, a rapid population expansion took place after the population bottleneck, which might have contributed to the accumulation of mutations.

4.4. Evolutionary history of Rhododendron meddianum and R. cyanocarpum

STRUCTRUE, PCA and NJ trees analysis indicated that genetic structure ofRhododendron meddianum and R. cyanocarpum populations differed clearly. These two species, which are distributed in Western Yunnan, are divided by rivers and mountains (elev. > 2500 m) (Fig. S1). Topographical uplift may have isolated populations of R. meddianum and R. cyanocarpum, shaping their patterns of genetic diversity.

We estimated that Rhododendron meddianum and R. cyanocarpum diverged around 3.05 (2.21–5.03) Ma. This date is consistent with the conclusion of Mline et al. (2004), who estimated that most species of the Rhododenron subgenus Hymenanthes (mainly distribute in Himalayas and southern China) diverged from one another around 3–5 Ma. Neogene climate reconstructions of Western Yunnan indicate that the climate during the late Pliocene was slightly warmer and much wetter than the present climate (Su et al., 2013a, b). The Asian winter monsoon became significantly stronger after the late Pliocene, which was characterized by drier winters in southwestern China (Su et al., 2013a; Wang, 2006). Our field surveys indicated that the climate conditions of R. meddianum populations are drier than those of R. cyanocarpum populations. This ecological difference may have promoted genetic differentiation and led to speciation (Spicer et al., 2020). Therefore, we hypothesize that the intensification of the Asian winter monsoon after the late Pliocene may have contributed to genetic divergence of R. meddianum and R. cyanocarpum.

4.5. Conservation implications

Our study on the genetic structure of Rhododendron meddianum provide important conservation implications for this narrowly distributed species. The most effective approach to conserving endangered species is in situ conservation (Wu et al., 2015). In wild populations of R. meddianum, seedlings and saplings are rarely found, and its habitat is disturbed and destroyed by various factors, including wind power and grazing. Therefore, it is imperative to establish conservation plots to protect its natural habitat. Meanwhile, seed collection from the all three R. meddianum populations (CJ, PM and HQ) should be conducted as soon as possible so they can be used for germplasm storage and ex situ conservation. Moreover, because moderate genetic differentiation and strong genetic structure were revealed among populations, germplasm resources from each population should be not mixed and instead should be used separately to avoid the risk of outbreeding depression (Liu et al., 2020b).

Author contributions

XJZ, ZHL, and DTL collected plant samples. XJZ, XFL and YRC performed the experiments. XJZ analyzed the data and wrote the manuscript. HM and YPM designed and supervised the study, and also wrote and revised the manuscript. All authors read and approved the final manuscript.

Declaration of competing interest

The authors declare that there are no conflict of interest.

Acknowledgements

This work was supported by the Fundamental Research Funds for the Central Non-profit Research Institution of the Chinese Academy of Forestry (Grant No. CAFYBB2019ZB007); the Biodiversity Survey and Assessment Project of the Ministry of Ecology and Environment, China (Grant No. 2019HJ2096001006); the Ten Thousand Talent Program of Yunnan Province (Grant No. YNWR-QNBJ-2019-010; YNWR-QNBJ-2018-174); and the Young and Middle-aged Academic and Technical Leader Reserve Talent Project of Yunnan Province (Grant No. 2018HB066).

Appendix A. Supplementary data

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