Detuan Liu
Lu Zhang2
Jihua Wang
Yongpeng Ma- 1Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, China
- 2Flower Research Institute, Yunnan Academy of Agriculture Sciences, Kunming, China
Though it is well-acknowledged that next generation sequencing (NGS) technologies can provide further insights into plant conservation management than traditional molecular markers, studies employing NGS to address conservation genomics and subsequent conservation strategies for threatened plants are still rare. Rhododendron is the largest genus of woody plants in China, and many species are threatened, however, to date there has been no conservation genetic research using NGS in this genus. In the present study, we investigated the conservation genetics of R. cyanocarpum, a threatened species endemic to the Cangshan Mountains in Yunnan, China, using a double digest restriction-site-associated DNA-sequencing (ddRAD-seq) approach. Due to the availability of sufficient SNPs, we were able to distinguish between neutral and putatively selected SNPs and were able to further investigate the genetic diversity, population structure, and differentiation in R. cyanocarpum, as well as make an estimation of its demographic history. A total of 6,584 SNPs were obtained, of which 5,729 were neutral (detected using Tajima’s D). In terms of the 5,729 neutral SNPs, R. cyanocarpum had a higher genetic diversity (π = 0.0702 ± 0.0017, He = 0.0675 ± 0.0016) than other plant species assessed using Rad-seq methods, while population differentiation (Fst from 0.0314 to 0.0452) was weak. Interestingly, contrasting patterns of population structure were revealed from all neutral and selected SNPs, with distinct genetic clusters forming for all SNPs and neutral SNPs, but no distinct subgroups for selected ones. Moreover, we were able to detect changes in effective population size (Ne) of R. cyanocarpum from 150,000 years ago, including a bottleneck event ca. 60,000 years ago, followed by recovery of Ne over a short period, and a subsequent gradual decline in Ne to date. Implications for conserving R. cyanocarpum based on these main results are then discussed.
Introduction
Biodiversity conservation is key in the field of conservation, and the essence of all biodiversity is genetic variation (Culver et al., 2011). Indeed, one can argue that the conservation of biodiversity is ultimately the conservation of genetic diversity. Genetic diversity strongly influences the ability of plant species to persist in the face of threats (Hendricks et al., 2017), and loss of genetic diversity has been considered as a crucial factor that results in inbreeding depression, reduced adaptation and fitness, and a decrease in long-term species survival (Bruni et al., 2012). Endangered species with small isolated populations are at elevated risk of losing adaptive variation due to genetic drift and the genetic costs of inbreeding (Allendorf et al., 2007). It is therefore important to pay attention to the genetic diversity, genetic structure, demographic history, and genetic background of an endangered species before developing protection measures (Wang and Hu, 1996; Wu et al., 2017; Barbosa et al., 2018; Li et al., 2018). Understanding the extent of intraspecific genetic diversity and population genetic structure is important for planning conservation strategies for endangered species, and it is often necessary to investigate the causes of low levels of genetic variation in populations, especially for plant species with extremely small populations.
Population genetics, when based on a relatively wide distribution, can provide a rich and mathematically rigorous framework for understanding evolutionary processes in natural populations. Common marker types, such as microsatellites or AFLP, can only generate limited markers per sample (Hodel et al., 2016). With high-throughput sequencing technological advances, generally called Next Generation Sequencing (NGS), population genomics can now address evolutionary processes at a genomic scale in natural populations with thousands of genetic markers rather than a few genetic loci (Hohenlohe et al., 2010). Restriction site-associated DNA sequencing (RAD-seq), based on Illumina NGS technology, is relatively cheap and flexible, and has fueled studies in conservation genomics (Andrews et al., 2016). This has led to the discovery and genotyping of thousands of polymorphic genetic markers in non-model species, and the technique has now become a popular tool for the study of the selection and genetics of adaptation in natural populations (Andrews et al., 2016; Catchen et al., 2017). In addition, in datasets comprising fewer individuals, RAD-seq tends to outperform microsatellites, giving more reliable inferences on population structure and a higher resolution (Lemopoulos et al., 2019).
The genus Rhododendron L. (Ericaceae) contains more than 1,000 species and has a global distribution (Chamberlain et al., 1996). To date, a total of ca. 600 species are found in China (Ma, 2015; Tian et al., 2019), of which 70% are endemic and many are threatened, so the conservation of Rhododendron is urgent (Ma et al., 2014; Qin et al., 2017). Although there have been some reports on the population genetics of Rhododendron, few of these studies focused on threatened species (Ma et al., 2017). Furthermore, we are unaware of any population genetic studies employing RAD-seq methods in Rhododendron.
Rhododendron cyanocarpum Franch. ex W. W. Smith, Subgen. Hymenanthes (Blume) K. Koch and Subsect. Thomsonii Sleumer, is a threatened species endemic to the Cangshan Mountains in Dali, Yunnan Province, Southwest China (Zhang et al., 1998). R. cyanocarpum has four populations growing along the ridges of the Cangshan Mountains, and is classified as Vulnerable in the IUCN red list (World Conservation and Monitoring Centre, 1998), in the Red list of Rhododenrons (Gibbs et al., 2011) and in the Threatened Species List of China’s Higher Plants (Qin et al., 2017). R. cyanocarpum grows at altitudes of 3,400–3,900 m; 12 individuals were found at a 900 m2 plot (Zhang et al., 1998). Its main effective pollinators are bumblebees (Ma et al., 2015a), with a 10-km potential range (Goulson and Stout, 2001) and at least 1.5 km foraging range from their colonies (Osborne et al., 2008). In addition, birds might be another pollinator in the late flowering of R. cyanocarpum (Ma et al., 2015a). Temperature is a crucial factor influencing the survival rate of its seedlings (Zhang et al., 1998). Furthermore, polymorphism of flower color (Ma et al., 2015b), reproductive isolation, and natural hybridization were also reported (Ma et al., 2010a, b, 2016; Ma, 2010). However, to date, no study has investigated the population genetics of R. cyanocarpum. Here we hypothesized that, due to very similar habitats where R. cyanocarpum occurs, adaptive selection could tend to be harmonized among populations. Therefore, the population structure of R. cyanocarpum inferred from selected loci (if some can be detected) should show more similar patterns than neutral ones. Moreover, we supposed that historical demography (e.g., bottleneck) could play an important role in maintaining the current distribution of R. cyanocarpum.
In our study, we assessed the genetic diversity of R. cyanocarpum using ddRAD-seq (Peterson et al., 2012). Due to the availability of sufficient SNPs, we were able to investigate several aspects of the population genetics of R. cyanocarpum using neutral and putatively selected SNPs, including genetic diversity, population structure, and differentiation, as well as estimation of demographic history. We believe that both the analytical process and the results presented here will be useful in the future, not only for the management and conservation of R. cyanocarpum, but also as an example of population genetic studies employing NGS for the conservation of endangered plants, including the new conservation action concept of plant species with extremely small populations (PSESP) in China (Ma et al., 2013; Sun, 2013; Sun et al., 2019).
Materials and Methods
Sample Collection and Extraction of DNA
Rhododendron cyanocarpum is endemic to the Cangshan Mountains, Yunnan, southwestern China. We sampled 59 individuals from four populations at locations around the Cangshan Mountains (Figure 1 and Table 1). Within each population, individuals spaced about 50 m apart were randomly selected. Leaf material was collected into silica gel in the field. Total genomic DNA was then extracted from the silica gel-dried leaf material using a CTAB procedure (Doyle and Doyle, 1990). In brief, gel-dried fragments of samples were put into liquid nitrogen and quickly ground to powder; the powder was then immediately transferred into a clean 1.5 ml microcentrifuge tube and mixed immediately with 700 μl 2 × CTAB extracting solution and incubated in a water bath at 65°C for 30 min. The supernatant was then removed, and the pellet was resuspended in chloroform-isoamyl alcohol (24:1). The samples were then centrifuged at 12,000 rpm, after which the supernatant was transferred to another clean 1.5 ml microcentrifuge tube and mixed with the same volume of isopropanol, then incubated at −20°C for about 1 h. The solution was then centrifuged at 12,000 rpm for 10 min, then the supernatant liquor was carefully discarded. The DNA was then cleaned twice with 75% ethanol and air dried at room temperature. The cleaned DNA was then dissolved in 50 μl TE, and 0.5 μl RNase was added to digest RNA. DNA quantity was assessed with 1.2% agarose gel electrophoresis and quantified with a qubit 3.0 fluorescent quantitative assay. DNA extraction and ddRAD-seq library preparation were performed by JieRui BioScience Co. Ltd. (Guangzhou, China).
Figure 1. Distribution and habitat of R. cyanocarpum. (A) Dali (red dot) in Yunnan Province (green), southwestern China; (B) population distribution of R. cyanocarpum on the Cangshan Mountains; (C) habitat, the white arrows indicate plants of R. cyanocarpum.
Table 1. Collection sites of Rhododendron cyanocarpum on the Cangshan Mountains, Yunnan, China, including population id, location name, coordinates, and sample size (number of individuals used in the analysis).
ddRAD Library Preparation and Sequencing
Library preparation was conducted following Peterson et al. (2012). The genomic DNA (100 ng 10 μl) was double digested using 10 μl of the restriction enzymes EcoR I and Mse I for 5 h at 37°C, then 20 min at 65°C, and final incubation at 12°C. The resulting digested fragments were cleaned and subsequently quantified using agarose gel electrophoresis. Digested fragments were ligated to EcoR I and Mse I adapters containing sample specific barcodes with T4 DNA ligase (NEB) for 4 h at 16°C, then 20 min at 65°C, and final incubation at 12°C. Individually barcoded samples were cleaned and size-selected (350–500 bp) using agarose gel (Omega kit). Each library was then PCR-amplified to the desired concentration and paired-end sequenced (0.5 G each sample) on an Illumina X-ten (Illumina) with PE 150 mode.
SNP Calling
Quality filtering and locus assembly were conducted using the Stacks software, version 2.4 (Catchen et al., 2013). RAD-tags were demultiplexed using process_radtags, the len_limit was set as 140 bp to trim low-quality reads, and retain_header -t was set to 135. We then used ustacks (parameter as follows: –min depth of coverage to create a stack (m): 2, –repeat removal algorithm: enabled, –max distance allowed between stacks (M): 2, –max distance allowed to align secondary reads: 4, –max number of stacks allowed per de novo locus: 3, –deleveraging algorithm: disabled, –gapped assembly: enabled, –minimum alignment length: 0.8, –model type: SNP, –alpha significance level for model: 0.05) to cluster and generate loci. All the loci were merged into the catalog using cstacks with the default parameters. The populations program (–min-populations: 4, –min-samples-per-pop: 0.8, –max-obs-het: 0.6, –write-single-snp) was used to call the SNPs across all ddRAD sites, the parameter “–write-single-snp” was used to restrict data analysis to only the first SNP per locus (Catchen et al., 2013; O’Leary et al., 2018).
Analysis of Genetic Diversity and Structure
Tajima’s D was calculated in vcftools v.0.1.16 (Danecek et al., 2011) with 95% confidence limit (−1.795 to 2.052) to test all the loci for neutrality (Tajima, 1989b; Anderson et al., 2016). Window size for vcftools in calculating the Tajima’s D statistic was set as 3,000 bp. The loci were then divided into two sets: neutral loci and non-neutral loci evolving under putative selection. BayeScan v.2.1 was also used to infer loci under selection in each of the four sampled populations (Foll and Gaggiotti, 2008). PGDSpider v.2.1.1.5 was used for format conversion for the subsequent analysis (Lischer and Excoffier, 2011).
Population genetic statistics, including the number of private alleles, heterozygosity (Ho), nucleotide diversity (π), and Wright’s F statistics Fst and FIS statistics, were calculated using the “populations” program in Stacks (Catchen et al., 2013), and were tested for significance among populations by SPSS 16.0 software. Normality and variance homogeneity were first tested in all loci (26,336 observations), neutral loci (22,916 observations), and selected loci (3,420 observations) independently; because they violated the assumption of One-way ANOVA, Nonparamatric Tests were performed among populations. Geographic distance was calculated between each pair of sample locations using scripts from the Movable Type Scripts webpage1.
To examine the population structure of R. cyanocarpum, a Bayesian-based analysis was performed using the software Structure v.2.3 (Pritchard et al., 2000), with a burn-in of 1,000 steps and 5,000 replicates for each value of K. The optimal K for each analysis was chosen using Harvester v.0.6942. Genetic relationships among the studied individuals were also assessed with a principal component analysis (PCA) performed using Plink v.1.9 (Purcell et al., 2007) and R v.3.6.1 (R Development Core Team, 2013) to identify the population structure and to show the first two major axes explaining genetic variation. For this analysis, a variant call format (VCF) file was also generated using the “populations” program in Stacks. An analysis of molecular variance (AMOVA) with 1,000 permutations for each population using all loci, neutral loci, and selected loci were performed separately in the program Arlequin v.3.5 (Excoffier and Lischer, 2010) to calculate pairwise Fst-values, with statistically significant differentiation being determined using a p-value of <0.05.
Demographic History
To investigate the recent demographic history of R. cyanocarpum, we used the python script easySFS3 to generate the folded site frequency spectrum (SFS) formatted file with VCF files, and the demographic history was then inferred using the program Stairway plot v.0.2 (Liu and Fu, 2015) based on neutral loci with the recommended 67% of sites for training. Stairway plot is based on the site frequency spectrum (SFS) that does not require whole-genome sequence data or a reference genome, and is more accurate for inferring recent population size changes compared to the PSMC or MSMC methods (Liu and Fu, 2015). However, this method can be affected by the estimation of mutation rate (Lapierre et al., 2017). Patton et al. (2019) used four methods to infer the history of population size from genomic datasets and compared their performance. Their results suggested that the Stairway Plot has the most accurate estimation of the shape of recent trends, but that it tends to underestimate contemporary population sizes.
Cross (1975) observed that R. ponticum began to flower at about 12 years old; Tamaki et al. (2017) gave the generation time for R. japonoheptamerum as 30 years, as they described that R. japonoheptamerum is a shrub with a very low growth rate. To infer the recent demographic history, we therefore set the generation time for R. cyanocarpum to 10, 20, and 30 years per generation. The mutation rate of R. weyrichii had been determined to be 1.581 × 10–9 per site per year by Yoichi et al. (2016), and we therefore set a mutation rate of 1.581 × 10–8 per site per generation (generation time = 10), 3.162 × 10–8 per site per generation (generation time = 20) and 4.743 × 10–8 per site per generation (generation time = 30) for R. cyanocarpum.
Results
Sequence Data Quality and Processing
After all quality filters, 3,337,583 reads with no rad tag and 12,793 reads with low quality were discarded. We retained a total of 344,001,104 reads from the initial 347,351,480 raw reads, with an average of 5,830,527 reads per sample. After trimming and clustering, we obtained 9,148,349 loci. Mean locus coverage across all samples was 23.91×, ranging from 17.77× to 58.41×, with an average length of 136 bp per loci and with 40–48% GC content. After the “cstacks” module process, we obtained 2,464,683 catalogs (Supplementary Table S1). Finally, 6,584 SNPs were retained.
Genetic Diversity and Population Structure
After filtering using Tajima’s D in Vcftools, a total of 855 SNPs evolving non-randomly (”selected”), and 5,729 SNPs evolving randomly (”neutral”) were retained. The Tajima’s D value for the whole dataset at the species level was –0.2748. Significantly negative Tajima’s D values were found, indicating that the haplotype frequencies deviated from the neutrality model (Tajima, 1989b). Selected loci and neutral loci were also tested based on Bayescan v.2.1. Only two loci were detected to be under selection (q < 0.05) with this method. Hence, our final datasets for further analysis were all loci (6,582 neutral SNPs identified by Bayescan), neutral loci (5,729), and selected loci (855) identified by Tajima’s D.
The selected loci dataset has the lowest number of private alleles in each population (157–221), neutral loci dataset is in the middle (911–1136), and all loci dataset has the highest number of private alleles (1081–1351). The genetic parameters He and π showed significant difference among the four populations for all loci (Ho: p = 0.004; He: p = 0.000; π: p = 0.000; FIS: p = 0.001), neutral loci (Ho: p = 0.049; He: p = 0.000; π: p = 0.000; FIS: p = 0.001), and selected loci (Ho: p = 0.001; He: p = 0.000; π: p = 0.000; FIS: p = 0.074). The inbreeding coefficients were positive in the four populations, reflecting a deficit of heterozygotes (inbreeding) (Boscari et al., 2019). The summary statistics of genetic diversity using all loci, neutral loci, and selected loci are given in Table 2.
Table 2. Population genetic statistics of R. cyanocarpum on the Cangshan Mountains, Yunnan, China: observed heterozygosity (Ho), expected heterozygosity (He), genetic diversity (π), inbreeding coefficients (FIS).
The pairwise Fst values between populations were less than 0.05 (Table 3), indicating that all populations of R. cyanocarpum were genetically similar and there was little differentiation among them (Ewens, 1969), which was also supported by the analysis of molecular variance (AMOVA, Table 4). The global AMOVA revealed that most genetic variation was found within populations, with very little among populations.
Table 3. Genetic distances (Fst values, above diagonal) and geographic distances (km, below diagonal) between R. cyanocarpum populations on the Cangshan Mountains, Yunnan, China.
Table 4. AMOVA results: evaluation of genetic differentiation within and among sampling sites of R. cyanocarpum on the Cangshan Mountains, Yunnan, China.
ΔK = 2 was best supported for all three types of loci (Supplementary Figure S1). Similar patterns of population structure were inferred from both all loci and neutral loci, while DST had fewer admixture characteristics deduced from all loci (an average of 88% green vs. 12% red genetic backgrounds) than neutral loci (an average of 60% green vs. 40% red genetic backgrounds). Two clear genetic clusters are visible from the Structure results using all loci (Figure 2A) and neutral loci (Figure 2B), with populations HDB and XHD belonging to cluster 1, and populations DST and DHY belonging to cluster 2. However, the pattern revealed from selected loci differed from the other two types of loci, as only one genetically homogeneous group was detected (Figure 2C). Similar patterns were suggested using principal coordinate analysis. Three clusters were detected using all loci (Figure 2D) and neutral loci (Figure 2E), one for HDB and XHD, and one for DST and one for DHY. Additionally, the four populations were clustered to one group using PCA based on selected loci (Figure 2F).
Figure 2. Population structure of R. cyanocarpum on the Cangshan Mountains, Yunnan, China. (A–C) Structure results (K = 2), with different colors showing different genetic backgrounds [(A) using all loci, (B) using neutral loci, (C) using selected loci]; (D–F) results of principal components analysis [(D) using all loci, (E) using neutral loci, (F) using selected loci].
Effective Population Size and Demographic History
Effective population size and demographic history of R. cyanocarpum were inferred based on SNP frequency spectra and displayed on a stairway plot. A population bottleneck occurred at ca. 0.06 Ma (million years ago), after which followed a population expansion and then a gradual reduction in population size to ca. 100 years ago (Figure 3). Despite the differences observed in the effective population size, we found similar demographic patterns when using generation times of 30, 20, and 10 years.
Figure 3. Reconstructing the demographic history of R. cyanocarpum on the Cangshan Mountains, Yunnan, China. Stairway plot showing historical changes in effective population size (Ne × 1000, y-axis) for R. cyanocarpum population with generation times of 10 years (red line), 20 years (blue line), and 30 years (green line).
Discussion
Genetic Diversity
It is generally believed that characterizing genetic diversity and population structure are essential for the effective conservation of threatened species. The genetic diversity of R. cyanocarpum (0.0634–0.0683 with all loci, 0.0682–0.0730 with neutral loci) was higher compared with previous studies on other plant species using ddRAD sequencing. For example, a much lower level was found in endangered Viola uliginosa (0.013–0.023) (Lee et al., 2020), Clermontia fauriei (0.0014), and Cyanea pilosa (0.0012) populations (Jennings et al., 2016); we are not aware of other studies on rhododendrons using ddRAD-seq. We can therefore conclude that genetic diversity in Rhododendron cyanocarpum is higher despite its small population size. The higher genetic diversity seen in R. cyanocarpum could be due to its population demography. As a woody species growing in an alpine habitat, R. cyanocarpum is likely to have a long lifespan; for example, the generation time of R. japonoheptamerum is about 30 years (Tamaki et al., 2017). It has been recorded that ca. 60 years ago, many R. cyanocarpum trees were cut down to aid the building of the Huadianba Medical Factory (Ma et al., 2010a), close to the HDB and XHD populations, which probably resulted in a recent reduction in population size of R. cyanocarpum. In such a scenario, however, we would expect that R. cyanocarpum would still maintain a similar amount of genetic diversity as before, due to the very limited role played by genetic drift over such a short time period (Luan et al., 2006).
Genetic Structure and Genetic Differentiation
Gene flow is a main factor affecting genetic differentiation between populations. The low genetic differentiation among populations of R. cyanocarpum could be explained by gene flow resulting from long-distance gene dispersal either as pollen or as seeds (Zhao et al., 2012). All of the study populations of R. cyanocarpum were close geographically, with a distance range of 6.9–23.9 km, which is beneficial for gene flow facilitated by pollinators and seed dispersal. Previous studies suggested that birds were another important pollinator in R. cyanocarpum (Ma et al., 2015a), which can favor pollen flow between populations because of birds’ long-distance flight. Rhododendron seeds, characterized as small, light and winged, are frequently dispersed by wind (Ng and Corlett, 2000).
For the results using Structure and PCA, similar patterns were detected based on selected loci and neutral loci in R. cyanocarpum. When using all loci or neutral loci, the results of Structure and PCA revealed two or three distinct groups in R. cyanocarpum populations, respectively. In contrast, the patterns of population structure resulting from the selected loci dataset showed no distinct subgroups, which could be explained by natural selection. At small spatial scales, local environmental conditions often shape the patterns of genetic structure, especially in heterogeneous and fragmented habitats (Young et al., 1996; Zhu et al., 2016). While there are distances of 6–23 km among the populations, a valley within the mountain range, as well as other topographic variation, may represent stabilizing selection for R. cyanocarpum adaptation to shape the patterns of genetic structure.
Effective Population Size and Demographic History
According to our stairway plot to illustrate the demographic history of R. cyanocarpum, a population bottleneck occurred at about 0.058 Ma, which was consistent with the late Pleistocene extinction (ca. 0.05 Ma) (Stuart, 1999; Van Der Kaars et al., 2017), and was also consistent with the process termed the early Dali Glaciation (0.058 Ma), when the estimated snow line was lower than it currently is, glacial and periglacial processes were taking place in areas above 3,500 m (Yang et al., 2006; Wan et al., 2011), and the temperatures were lower than at present. Extremely low temperatures may seriously damage and inhibit growth of R. cyanocarpum. During our field investigations, we saw that R. cyanocarpum grew up to altitudes of about 3,860 m.
The stairway plot also suggested that R. cyanocarpum underwent a population expansion in the period before and after the Last Glacial Maximum (LGM, 0.04–0.006 Ma) (Garot et al., 2019). That this population expansion event took place was also supported by the negative value of Tajima’s D (–0.2748 with 95% confidence level), because strong negative values of Tajima’s D can indicate population expansions (Tajima, 1989a). Recently, over the last 7,000 years, the Ne of R. cyanocarpum has been declining, which is consistent with environmental exploitation, deforestation, clearance of forested areas, and other intensive human activities. The first evidence for anthropogenic disturbance on the Cangshan Mountains is from ca. 6,370 years ago based on pollen records, when forests were being cleared for shifting agriculture (Shen et al., 2006). Together with climatic warming, this deforestation could have led to the decline in the R. cyanocarpum populations. Further evidence of disturbance caused by anthropogenic activities on the Cangshan Mountains stems from Neolithic ruins (about 3,100 years ago), the Bronze Age (from 2,000 years ago) (Wan et al., 2011), the inward migration and establishment of the Han Chinese (about 2,300–1,700 years ago), and the ancient city of Nanzhao and the Dali Kingdom (from 1,000 years ago) (Dearing et al., 2008). Taken together, frequent anthropogenic activities, such as agriculture, pastoralism, irrigation, vegetation clearance and deforestation, and burning and grazing, as well as the onset of warmer summers (Dearing et al., 2008), all posed great threats to R. cyanocarpum and resulted in reduction of effective population sizes.
Implications for Conservation and Management of R. cyanocarpum
R. cyanocarpum underwent a bottleneck, a subsequent expansion, and a gradual reduction. With intensive anthropogenic activity and global warming, the species has become fragmented and now has a highly restricted range near the summit of the Cangshan Mountains. Based on our main findings, we recommend the following conservation measures for R. cyanocarpum. (i) Although all populations are located in a national nature reserve, habitat disturbance from anthropogenic activity still occurs frequently, including tourism, path construction, and pasture and land reclamation. Absolute prohibition of anthropogenic activities is not possible and is also not necessary, although in the national nature reserve deforestation and large-scale land reclamation should be banned absolutely, especially near the HDB and XHD populations, where there is a big pharmaceutical farm. (ii) Because most individuals of R. cyanocarpum are situated on ridges of the Cangshan Mountains, and no distinct subgroups were detected when using selected loci, we suppose that the most serious threat to R. cyanocarpum will be not anthropogenic activities but instead be global warming. Given the topographic and landscape constraints, the species has limited migration potential, and ex situ conservation is therefore urgent and should be conducted immediately, including seed collection and conservation, and searching for new suitable habitats. (iii) The DHY population had the highest level of genetic diversity of the four populations, and we observed more connectivity and admixture between DHY and other populations than between any of the others, suggesting its importance in maintaining genetic polymorphism. This population is geographically located in the middle of other populations, which is beneficial for facilitating dispersal and connecting the other populations as a corridor. Therefore, we suggest that more conservation effort be given to the population DHY.
Data Availability Statement
The raw sequencing data has been successfully uploaded to SRA in NCBI, https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA640883.
Author Contributions
YM conceived and designed the experiments and revised the manuscript. LZ carried out investigation and collected materials. DL analyzed the data and wrote the manuscript. YM and JW acquired the funding. All authors contributed to the article and approved the submitted version.
Funding
The authors acknowledge financial supports from the General Program of National Natural Science Foundation of China (Grant No. 31770418), the Ministry of Science and Technology granted funding for the national key program (Grant No. 2017YFC0505200), the Construction of International Flower Technology Innovation Center and Industrialization of Achievements (Grant No. 2019ZG006), the Youth Innovation Promotion Association, Chinese Academy of Sciences (Grant No. 2018428), the Reserve Talents for Academic and Technical Leaders of Middle-aged and Young People in Yunnan Province (Grant No. 2018HB066), and the Program of Science and Technology Talents Training in Yunnan Province awarded to JW (Grant No. 2016HA005).
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We would like to thank Yubing Zhou and Dr. Hantao Qin for her help with analysis. We also thank the editor and three reviewers for their comments and good suggestion.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2020.00757/full#supplementary-material
Footnotes
- ^ http://www.movable-type.co.uk/scripts/latlong.html
- ^ http://taylor0.biology.ucla.edu/structureHarvester
- ^ https://github.com/isaacovercast/easySFS
References
Allendorf, F. W., Luikart, G., and Aitken, S. N. (2007). “Small populations and genetic Drift,” in Conservation and the Genetics of Populations, (Washington, D.C.: USA), 117–146.
Anderson, C., Tay, W. T., McGaughran, A., Gordon, K., and Walsh, T. K. (2016). Population structure and gene flow in the global pest, Helicoverpa armigera. Mol. Ecol. 25, 5296–5311. doi: 10.1111/mec.13841
Andrews, K. R., Good, J. M., Miller, M. R., Luikart, G., and Hohenlohe, P. A. (2016). Harnessing the power of RADseq for ecological and evolutionary genomics. Nat. Rev. Genet. 17, 81–92. doi: 10.1038/nrg.2015.28
Barbosa, S., Mestre, F., White, T. A., Pauperio, J., Alves, P. C., and Searle, J. B. (2018). Integrative approaches to guide conservation decisions: Using genomics to define conservation units and functional corridors. Mol. Ecol. Notes 27, 3452–3465. doi: 10.1111/mec.14806
Boscari, E., Abbiati, M., Badalamenti, F., Bavestrello, G., Benedetti-Cecchi, L., Cannas, R., et al. (2019). A population genomics insight by 2b-RAD reveals populations’ uniqueness along the Italian coastline in Leptopsammia pruvoti (Scleractinia, Dendrophylliidae). Divers. Distrib. 25, 1101–1117. doi: 10.1111/ddi.12918
Bruni, I., De Mattia, F., Labra, M., Grassi, F., Fluch, S., Berenyi, M., et al. (2012). Genetic variability of relict Rhododendron ferrugineum L. populations in the Northern Apennines with some inferences for a conservation strategy. Plant Biosys. 146, 24–32. doi: 10.1080/11263504.2011.557093
Catchen, J., Hohenlohe, P. A., Bassham, S., Amores, A., and Cresko, W. A. (2013). Stacks: an analysis tool set for population genomics. Mol. Ecol. 22, 3124–3140. doi: 10.1111/mec.12354
Catchen, J. M., Hohenlohe, P. A., Bernatchez, L., Funk, W. C., Andrews, K. R., and Allendorf, F. W. (2017). Unbroken: RADseq remains a powerful tool for understanding the genetics of adaptation in natural populations. Mol. Ecol. Resour. 17, 362–365. doi: 10.1111/1755-0998.12669
Chamberlain, D., Hyam, R., Argent, G., Fairweather, G., and Walter, K. S. (1996). The Genus Rhododendron: Its Classification and Synonymy. Edinburgh: Royal Botanic Garden.
Cross, J. (1975). Biological flora of the British Isles. Rhododendron ponticum L. J. Ecol. 63, 345–364.
Culver, M., Fitak, R., and Herrmann, H. W. (2011). “Genetic methods for biodiversity assessment,” in Biological Diversity: Frontiers in Measurement and Assessment, eds A. E. Magurran and B. J. McGill (Oxford: Oxford University Press), 208–218.
Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., et al. (2011). The variant call format and VCFtools. Bioinformatics 27, 2156–2158. doi: 10.1093/bioinformatics/btr330
Dearing, J. A., Jones, R., Shen, J., Yang, X., Boyle, J., Foster, G., et al. (2008). Using multiple archives to understand past and present climate-human-environment interactions: the lake Erhai catchment, Yunnan Province, China. J. Paleolimnol. 40, 3–31. doi: 10.1007/s10933-007-9182-2
Ewens, W. J. (1969). Evolution and the genetics of populations: Vol. 2. The theory of gene frequencies. Science 168, 722–723. doi: 10.1126/science.168.3932.722
Excoffier, L., and Lischer, H. E. (2010). Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol. Ecol. Resour. 10, 564–567. doi: 10.1111/j.1755-0998.2010.02847.x
Foll, M., and Gaggiotti, O. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180, 977–993. doi: 10.1534/genetics.108.092221
Garot, E., Joët, T., Combes, M. C., Severac, D., and Lashermes, P. (2019). Plant population dynamics on oceanic islands during the Late Quaternary climate changes: genetic evidence from a tree species (Coffea mauritiana) in Reunion Island. New Phytol. 224, 974–986. doi: 10.1111/nph.16052
Gibbs, D., Chamberlain, D., and Argent, G. (2011). “Rhododendron cyanocarpum (Franchet) Franchet,” in In: The Red List of Rhododendrons, ed. W. W. Smith (Richmond: Botanic Gardens Conservation International), 32.
Goulson, D., and Stout, J. C. (2001). Homing ability of the bumblebee Bombus terrestris (Hymenoptera: Apidae). Apidologie 32, 105–111. doi: 10.1051/apido:2001115
Hendricks, S., Epstein, B., Schönfeld, B., Wiench, C., Hamede, R., Jones, M., et al. (2017). Conservation implications of limited genetic diversity and population structure in Tasmanian devils (Sarcophilus harrisii). Conserv. Genet. 18, 977–982. doi: 10.1007/s10592-017-0939-5
Hodel, R. G., Segovia-Salcedo, M. C., Landis, J. B., Crowl, A. A., Sun, M., Liu, X., et al. (2016). The report of my death was an exaggeration: A review for researchers using microsatellites in the 21st century. Appl. Plant Sci. 4:1600025. doi: 10.3732/apps.1600025
Hohenlohe, P. A., Bassham, S., Etter, P. D., Stiffler, N., Johnson, E. A., and Cresko, W. A. (2010). Population genomics of parallel adaptation in threespine stickleback using sequenced RAD tags. PLoS Genetics 6:e1000862. doi: 10.1371/journal.pgen.1000862
Jennings, H., Wallin, K., Brennan, J., Valle, A. D., Guzman, A., Hein, D., et al. (2016). Inbreeding, low genetic diversity, and spatial genetic structure in the endemic Hawaiian lobeliads Clermontia fauriei and Cyanea pilosa ssp. longipedunculata. Conserv. Genet. 17, 497–502. doi: 10.1007/s10592-015-0785-2
Lapierre, M., Lambert, A., and Achaz, G. (2017). Accuracy of demographic inferences from the site frequency spectrum: the case of the yoruba population. Genetics 206, 439–449. doi: 10.1534/genetics.116.192708
Lee, K. M., Ranta, P., Saarikivi, J., Kutnar, L., Vreš, B., Dzhus, M., et al. (2020). Using genomic information for management planning of an endangered perennial, Viola uliginosa. Ecol. Evol. 10, 2638–2649. doi: 10.1002/ece3.6093
Lemopoulos, A., Prokkola, J. M., Uusi-Heikkilä, S., Vasemägi, A., Huusko, A., Hyvärinen, P., et al. (2019). Comparing RADseq and microsatellites for estimating genetic diversity and relatedness—Implications for brown trout conservation. Ecol. Evol. 9, 2106–2120. doi: 10.1002/ece3.4905
Li, X., Li, M., Hou, L., Zhang, Z. Y., Pang, X. M., and Li, Y. Y. (2018). De novo transcriptome assembly and population genetic analyses for an endangered Chinese endemic Acer miaotaiense (Aceraceae). Genes 9, 1–20. doi: 10.3390/genes9080378
Lischer, H. E., and Excoffier, L. (2011). PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. Bioinformatics 28, 298–299. doi: 10.1093/bioinformatics/btr642
Liu, X., and Fu, Y.-X. (2015). Exploring population size changes using SNP frequency spectra. Nature genetics 47, 555–559. doi: 10.1038/ng.3254
Luan, S. S., Chiang, T. Y., and Gong, X. (2006). High genetic diversity vs. low genetic differentiation in Nouelia insignis (Asteraceae), a narrowly distributed and endemic species in China, revealed by ISSR fingerprinting. Ann. Bot. 98, 583–589. doi: 10.1093/aob/mcl129
Ma, H., Li, T. Q., Liu, X. F., Wan, Y. M., Li, Y. Y., Liu, X. X., et al. (2017). Research progress in conservation biology of Rhododendron. World Forestry Research 30, 13–17.
Ma, Y. P. (2010). The study of natural hybridization between Rhododendron cyanocarpum Franch. ex W.W. and R. delavayi Franch. Beijing: Chinese Academy of Sciences.
Ma, Y. P. (2015). A new species of Rhododendron (Ericaceae) from Shangri-La. NW Yunnan, China. Phytotaxa 238, 293–297. doi: 10.11646/phytotaxa.238.3.10
Ma, Y. P., Chen, G., Grumbine, R. E., Dao, Z. L., Sun, W. B., and Guo, H. J. (2013). Conserving plant species with extremely small populations (PSESP) in China. Biodivers. Conserv. 22, 803–809. doi: 10.1007/s10531-013-0434-3
Ma, Y. P., Milne, R. I., Zhang, C. Q., and Yang, I. B. (2010a). Unusual patterns of hybridization involving a narrow endemic Rhododendron species (Ericaceae) in Yunnan, China. Am J Bot 97, 1749–1757. doi: 10.3732/ajb.1000018
Ma, Y. P., Zhang, C. Q., Zhang, J. L., and Yang, J. B. (2010b). Natural hybridization between Rhododendron delavayi and R. cyanocarpum (Ericaceae), from morphological, molecular and reproductive evidence. J. Integr. Plant Biol. 52, 844–851. doi: 10.1111/j.1744-7909.2010.00970.x
Ma, Y. P., Nielsen, J., Chamberlain, D. F., Li, X. Y., and Sun, W. B. (2014). The conservation of Rhododendrons is of greater urgency than has been previously acknowledged in China. Biodivers. Conserv. 23, 3149–3154. doi: 10.1007/s10531-014-0764-9
Ma, Y. P., Wu, Z. K., Dong, K., Sun, W. B., and Marczewski, T. (2015a). Pollination biology of Rhododendron cyanocarpum (Ericaceae): An alpine species endemic to NW Yunnan, China. J. Syst. Evol. 53, 63–71. doi: 10.1111/jse.12114
Ma, Y. P., Wu, Z. K., Zhang, C. Q., and Sun, W. B. (2015b). Flower color polymorphism in Rhododendron cyanocarpum (Ericaceae), an endangered alpine species endemic to NW Yunnan, China. Plant Diver. 37, 21–28. doi: 10.7677/ynzwyj201514068
Ma, Y. P., Xie, W. J., Sun, W. B., and Marczewski, T. (2016). Strong reproductive isolation despite occasional hybridization between a widely distributed and a narrow endemic Rhododendron species. Sci. Rep. 6, 19146. doi: 10.1038/srep19146
Ng, S. C., and Corlett, R. T. (2000). Genetic variation and structure in six Rhododendron species (Ericaceae) with contrasting local distribution patterns in Hong Kong, China. Mol. Ecol. 9, 959–969. doi: 10.1046/j.1365-294x.2000.00958.x
O’Leary, S. J., Puritz, J. B., Willis, S. C., Hollenbeck, C. M., and Portnoy, D. S. (2018). These aren’t the loci you’re looking for: Principles of effective SNP filtering for molecular ecologists. Mol. Ecol. 27, 1–14. doi: 10.1111/mec.14792
Osborne, J. L., Martin, A. P., Carreck, N. L., Swain, J. L., Knight, M. E., and Goulson, D. (2008). Bumblebee flight distances in relation to the forage landscape. J. Anim. Ecol. 77, 406–415. doi: 10.1111/j.1365-2656.2007.01333.x
Patton, A. H., Margres, M. J., Stahlke, A. R., Hendricks, S., Lewallen, K., Hamede, R. K., et al. (2019). Contemporary demographic reconstruction methods are robust to genome assembly quality: A case study in Tasmanian devils. Mol. Biol. Evol. 36, 2906–2921. doi: 10.1093/molbev/msz191
Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., and Hoekstra, H. E. (2012). Double digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One 7:e37135. doi: 10.1371/journal.pone.0037135
Pritchard, J. K., Stephens, M., and Donnelly, P. (2000). Inference of population structure using multilocus genotype data. Genetics 155, 945–959.
Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M. A., Bender, D., et al. (2007). PLINK: a tool set for whole-genome association and population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575. doi: 10.1086/519795
Qin, H. N., Yang, Y., Dong, S. Y., He, Q., Jia, Y., Zhao, L. N., et al. (2017). Threatened species list of China’s higher plants. Biodiversity science 25, 696–744.
R Development Core Team (2013). R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing.
Shen, J., Jones, R. T., Yang, X. D., Dearing, J. A., and Wang, S. M. (2006). The Holocene vegetation history of Lake Erhai, Yunnan province southwestern China: the role of climate and human forcings. Holocene 16, 265–276. doi: 10.1191/0959683606hl923rp
Stuart, A. J. (1999). “Late Pleistocene megafaunal extinctions,” in Extinctions in Near Time: Causes, Contexts, and Consequences. eds R. D. E. MacPhee and H. D. Sues (New York, NY: Springer Press), 257–269. doi: 10.1007/978-1-4757-5202-1_11
Sun, W. B. (2013). Conservation Practice and Exploration of Wild Plant Species with Extremely Small Populations in Yunnan Province, (Kunming: Yunnan Science and Technology Press), 1–7.
Sun, W. B., Ma, Y. P., and Blackmore, S. (2019). How a new conservation action concept has accelerated plant conservation in China. Trends Plant Sci. 24, 4–6. doi: 10.1016/j.tplants.2018.10.009
Tajima, F. (1989a). The effect of change in population size on DNA polymorphism. Genetics 123, 597–601.
Tajima, F. (1989b). Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123, 585–595.
Tamaki, I., Yoichi, W., Matsuki, Y., Suyama, Y., and Mizuno, M. (2017). Inconsistency between morphological traits and ancestry of individuals in the hybrid zone between two Rhododendron japonoheptamerum varieties revealed by a genotyping-by-sequencing approach. Tree Genet. Genom. 13, 4. doi: 10.1007/s11295-016-1084-x
Tian, X. L., Chang, Y. H., Neilsen, J., Wang, S. H., and Ma, Y. P. (2019). A new species of Rhododendron (Ericaceae) from northeastern Yunnan, China. Phytotaxa 395, 66–70. doi: 10.11646/phytotaxa.395.2.2
Van Der Kaars, S., Miller, G. H., Turney, C. S., Cook, E. J., Nürnberg, D., Schönfeld, J., et al. (2017). Humans rather than climate the primary cause of Pleistocene megafaunal extinction in Australia. Nat. Commun. 8, 1–7. doi: 10.1038/ncomms14142
Wan, Y., Fu, K., Liu, Y., and Shi, Z. (2011). Geomorphologic structure, characteristics and processes in the Cangshan mountains: Explanations for the formation and development of the Dali Glaciation. Internat. J. Geosci. 2, 155–163. doi: 10.4236/ijg.2011.22016
Wang, H. X., and Hu, Z. A. (1996). Plant breeding system, genetic structure and conservation of genetic diversity. Biodiv. Sci. 1996, 92–96.
World Conservation and Monitoring Centre (1998). Rhododendron cyanocarpum. IUCN Red List Threatened Species. 1998:e.T32448A9707763.
Wu, F. Q., Shen, S. K., Zhang, X., Yang, G. S., and Wang, Y. H. (2017). Inferences of genetic structure and demographic history of Rhododendron protistum var. giganteum—The world’s largest Rhododendron using microsatellite markers. Flora 233, 1–6. doi: 10.1016/j.flora.2017.04.009
Yang, J., Zhang, W., Cui, Z., Yi, C., Liu, K., Ju, Y., et al. (2006). Late Pleistocene glaciation of the Diancang and Gongwang Mountains, southeast margin of the Tibetan Plateau. Quat. Int. 154, 52–62. doi: 10.1016/j.quaint.2006.02.003
Yoichi, W., Tamaki, I., Sakaguchi, S., Song, J. S., Yamamoto, S. I., and Tomaru, N. (2016). Population demographic history of a temperate shrub, Rhododendron weyrichii (Ericaceae), on continental islands of Japan and South Korea. Ecol. Evol 6, 8800–8810. doi: 10.1002/ece3.2576
Young, A., Boyle, T., and Brown, T. (1996). The population genetic consequences of habitat fragmentation for plants. Trends Ecol. Evol. 11, 413–418. doi: 10.1016/0169-5347(96)10045-8
Zhang, C. Q., Feng, B. J., Lv, Y. L., Zhou, B., and Gao, L. M. (1998). Research on causes for endangerment of Rhododendron protistum var. giganteum and R. cyanocarpum. J. Nat. Resour. 13, 276–278.
Zhao, X., Ma, Y., Sun, W., Wen, X., and Milne, R. (2012). High genetic diversity and low differentiation of Michelia coriacea (Magnoliaceae), a critically endangered endemic in southeast Yunnan, China. Int. J. Mol. Sci. 13, 4396–4411. doi: 10.3390/ijms13044396
Keywords: Rhododendron cyanocarpum, ddRAD-seq, genetic diversity, population demography, conservation
Citation: Liu D, Zhang L, Wang J and Ma Y (2020) Conservation Genomics of a Threatened Rhododendron: Contrasting Patterns of Population Structure Revealed From Neutral and Selected SNPs. Front. Genet. 11:757. doi: 10.3389/fgene.2020.00757
Received: 13 April 2020; Accepted: 25 June 2020;
Published: 04 September 2020.
Edited by:
Genlou Sun, Saint Mary’s University, CanadaReviewed by:
Shih-Ying Hwang, National Taiwan Normal University, TaiwanCharles Roland Clement, National Institute of Amazonian Research (INPA), Brazil
Yuji Isagi, Kyoto University, Japan
Copyright © 2020 Liu, Zhang, Wang and Ma. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Jihua Wang, wjh0505@gmail.com; Yongpeng Ma, mayongpeng@mail.kib.ac.cn