Does the critically endangered Rhododendron amesiae deserve top priority for conservation?
Yi-Shan Ao1,2, Yu-Hang Chang1,2, De-Tuan Liu1,2, Yong-Bo Liu3,**, Yong-Peng Ma4,*     
1. Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, 650201, China;
2. University of the Chinese Academy of Sciences, Beijing, 100049, China;
3. State Key Laboratory of Environmental Criteria and Risk Assessment, Chinese Research Academy of Environmental Sciences, Beijing, 100012, China;
4. Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany, Chinese Academy of Sciences, Kunming, 650201, China

The genus Rhododendron is the largest genus of woody plants in China, having extremely important horticultural value (Chamberlain et al., 1996). However, efforts have only recently been made to assess the conservation status of the genus. According to conservation category assessments on the Red List of Rhododendrons (Gibbs et al., 2011) and the Threatened Species List of China's Higher Plants (Qin et al., 2017), 32.03% and 17.13% of rhododendrons are evaluated to be threatened. However, conservation research and actions have only been implemented in a limited number of species (e.g., Rhododendron protistum var. giganteum, Li et al., 2018; R. cyanocarpum, Liu et al., 2020; R. griersonianum, Ma et al., 2021; R. meddianum, Zhang et al., 2021; R. longipedicellatum, Cao et al., 2022).

In general, the species is the basic unit for conservation. If a threatened plant species is inaccurately delimited, conservation efforts might be unwarranted. This is likely to be particularly true for the taxonomically complicated Rhododendron, and especially for some endangered species with few diagnostic characters described from a small number of herbarium specimens (Marczewski et al., 2016; Li et al., 2018).

Rhododendron amesiae Rehd. has been designated as a critically endangered species in both the Red List of Rhododendrons and the Threatened Species List of China's Higher Plants. In addition, preliminary investigations have shown that it conforms to the characteristics of Plant Species with Extremely Small Populations (PSESP) (Sun, 2013; Ma et al., 2013). Specifically, it was found only in a few localities in Sichuan Province, China and each population has fewer than 1000 mature individuals. However, its close relative Rhododendron concinnum is widespread, and ranges across several provinces in China (Fig. S1). The main character separating the two species is the hispid petiole of R. amesiae. In the present study, we made a comprehensive comparison of R. amesiae and R. concinnum, employing both morphological and population genomic data, with the aim of providing insight into the conservation status of R. amesiae.

We investigated Rhododendron amesiae populations in the field from 2019 to 2021, and eventually found a total of 8 populations: Laoyulin (LYL), Pengtaxiang (PTX), Caiyuanzigou (CY), Guobayangou (GY), Zhonggu village (ZC), Mulonggou (ML), Siguniang Mountain (SG) and Jiajin Mountain (JJS), all in western Sichuan, China. We observed that hispid and non-hispid petioles can occur at the same time within the same individual (Fig. 1A). To evaluate the distribution pattern of hispid petioles within individuals, between individuals and within and among populations, we investigated five populations of R. amesiae and two populations of R. concinnum. In each population, we randomly selected 30 individuals from which five mature leaves were collected per individual.

Fig. 1 A, petioles with or without hispid hairs occurring within the same individual of Rhododendron amesiae; B, percentage of individuals with hispid hairs on all petioles (green), hispid hairs absent on all petioles (blue) and both hispid and non-hispid petioles present on the same individual (brown) in R. amesiae.

For double digest restriction site-associated DNA sequencing (ddRAD-seq), we collected 45 samples from 5 populations of Rhododendron amesiae and 12 samples from 2 populations of R. concinnum (Table S1). Total genomic DNA was extracted from silica-dried leaf tissue using the CTAB method (Doyle and Doyle, 1990). A ddRAD sequencing library was constructed following Peterson et al. (2012), using the restriction enzymes EcoR I and Mse I. Paired-end reads (150 bp) were generated on an Illumina Hiseq X-Ten platform (Illumina Inc., San Diego, CA, United States), generating 0.5 Gb data per sample.

Paired-end raw reads were demultiplexed, trimmed, and filtered using the process_radtags program in STACKS v.2.4, with the len_limit set to 140 bp to trim low-quality reads, and retain_header -t was set to 135. A reference-based assembly method was used to improve the accuracy of SNP calling and we used the BWA v.07.17 (Li and Durbin, 2009) MEM module with the default parameters to map the clean reads to the Rhododendron delavayi reference genome (Zhang et al., 2017). SAMtools v.1.3.1 was used to convert sequence alignment map (SAM) format files to sorted, indexed binary alignment map (BAM) files and the mapping rate of each sample to the reference genome was calculated using the -flagstat command (Li et al., 2009; Danecek et al., 2021). The ref_map.pl script, including the core programs of gstacks and populations in STACKS, was run to identify SNPs and generate population-level summary statistics. To remove the repeated sequences resulting from PCR amplification in the results of mapping, the --rm-pcr-duplicates parameter was employed and to mitigate the impact of missing data, we used the filtering parameter –r 0.7, which only included SNPs shared by at least 70% of individuals in a population. The -write-single-snp parameter was set to exclude linked loci and to keep only one SNP per locus.

The nucleotide diversity (π), expected heterozygosity (He), observed heterozygosity (Ho), inbreeding coefficient (FIS) and F statistics (FST) were calculated using the populations program in the STACKS pipeline. Population genetic structure analyses based on SNPs was performed using PCA and STRUCTURE v.2.3.4 (Pritchard et al., 2000). The optimal K was chosen using the delta-K method implemented in STRUCTURE HARVESTER (Earl and vonHoldt, 2012).

We did not find any individuals with hispid petioles in the two Rhododendron concinnum populations. For R. amesiae, all examined individuals from GY had hispid petioles. However, among the remaining 4 populations, i.e., LYL, PTX, KZ and SG, 30%–56.67% of the examined individuals in each population had both hispid and non-hispid petioles (Fig. 1B).

A total of 422 million reads were produced for all samples. After removing the NoRadTag reads and filtering out the low-quality reads, 413 million clean reads remained for processing (Table S2). The mapping rate of samples to the Rhododendron delavayi reference genome was 39.2% on average (Table S3). After applying the strict filtering parameters (see method above), the reference-based analysis performed using the STACKS pipeline produced 565 high-quality SNPs.

Analysis of genetic parameters showed that the number of private alleles ranged from 13 to 34, with an average of 22.67 per population; observed heterozygosity (Ho) ranged from 0.037 to 0.049 (average 0.041) and expected heterozygosity (He) from 0.098 to 0.124 (average 0.111) at the population level. High genetic diversity was seen in each population (π: 0.116–0.132, average 0.123); and inbreeding coefficients were positive in all populations (FIS: 0.147–0.295, average 0.215) (Table 1).

Table 1 Statistics summarizing the genetic parameters of different populations of Rhododendron amesiae and R. concinnum. (Private: Number of private alleles; Ho: Observed heterozygosity; He: Expected heterozygosity; π: genetic diversity; FIS: inbreeding coefficients).
Species name Pop. ID Private Ho He Π FIS
R. amesiae LYL 24 0.039 0.117 0.126 0.257
R. amesiae PTX 28 0.042 0.111 0.121 0.206
R. amesiae CY 19 0.037 0.104 0.116 0.190
R. amesiae GY 25 0.049 0.102 0.118 0.149
R. amesiae ML 34 0.041 0.123 0.130 0.266
R. amesiae SG 13 0.037 0.115 0.123 0.246
R. amesiae JJS 19 0.044 0.108 0.120 0.182
R. concinnum MNXY 14 0.040 0.098 0.117 0.147
R. concinnum LDXY 28 0.040 0.124 0.132 0.295
Mean 22.67 0.041 0.111 0.123 0.215

The pairwise FST between the two populations of Rhododendron concinnum was 0.088 while within R. amesiae, pairwise FST among populations ranged from 0.062 to 0.144 (Table S4). There were no significant differences between the pairwise FST in populations of the two species (all p > 0.05), except for the pairwise FST between JJS (R. amesiae) and LDXY (R. concinnum) populations (p = 0.049) (Table 2). The result of a Mantel test between FST and geographical distance (GD) revealed that these were not significantly correlated (R = 0.767, p = 0.589).

Table 2 The pairwise FST/p-value (above diagonal) and geographical distance (GD, below diagonal) between populations. The p-values describe whether the FST measure is statistically significant according to Fisher's Exact Test.
Latin name Pop. ID LYL PTX CY GY ML SG JJS MNXY LDXY
R. amesiae LYL 0.086/0.076 0.076/0.056 0.112/0.107 0.082/0.082 0.086/0.084 0.078/0.062 0.106/0.093 0.070/0.062
R. amesiae PTX 54.256 0.099/0.082 0.110/0.087 0.090/0.090 0.095/0.092 0.094/0.073 0.109/0.079 0.075/0.067
R. amesiae CY 14.111 40.896 0.137/0.112 0.088/0.081 0.094/0.085 0.101/0.074 0.128/0.095 0.077/0.066
R. amesiae GY 118.319 67.386 106.499 0.097/0.096 0.111/0.106 0.132/0.111 0.144/0.097 0.089/0.080
R. amesiae ML 142.613 89.584 128.605 60.489 0.062/0.050 0.084/0.078 0.103/0.098 0.070/0.067
R. amesiae SG 144.874 90.792 131.688 37.085 33.307 0.091/0.081 0.111/0.101 0.067/0.059
R. amesiae JJS 124.597 70.546 111.438 24.587 35.920 20.278 0.130/0.100 0.067/0.049
R. concinnum MNXY 135.024 175.341 147.056 214.778 261.245 250.673 232.413 0.088/0.073
R. concinnum LDXY 32.343 55.633 35.609 106.048 143.383 138.110 118.363 119.735

Bayesian clustering analysis was performed in STRUCTURE, and the optimal K was calculated to be 2 (Fig. S2). No clear population genetic structure between Rhododendron amesiae and R. concinnum was found when K = 2 (Fig. 2A). Within R. amesiae, two genetic groups were roughly clarified. One group included the LYL, PTX, and CY populations, with all individuals dominated by a single genetic component (Fig. 2A, marked in orange). The other group included the ML and SG populations, with all individuals dominated by a different genetic component (marked in blue). The remaining populations included individuals with an admixture of the two genetic groups.

Fig. 2 Population genetic structure analysis showing nine populations of the two species (Rhododendron amesiae and R. concinnum). A, population structure analysis with delta K = 2 (burn-in: 100, 000, MCMC: 100, 000, K = 1–10, each K value repeated 20 times); B, PCA result.

In accordance with the STRUCTURE results, PCA analysis also showed no clear genetic separation between Rhododendron amesiae and R. concinnum, with low discrimination ability by both principal coordinate 1 (accounting for 5.12% of the total variance) and principal coordinate 2 (accounting for 5.12% of the total variance) (Fig. 2B). Furthermore, unlike the two rough genetic groups revealed in the STRUCTURE analysis, no clear genetic clustering was found in any population of either R. amesiae or R. concinnum. Overall, our population genetic results showed a high degree of similarity between the two species.

Our study incorporated a detailed morphological investigation and population genetic structure analysis of 9 populations of two closely related Rhododendron species sampled from western Sichuan, China. Our findings indicate that the species boundary between R. amesiae and R. concinnum is unclear. Consequently, the conservation status of R. amesiae should be reevaluated. Firstly, the key morphological characteristic supposedly distinguishing the two species can occur within populations or even within a single individual. Secondly, no clear genetic differences were found between populations of the two species. We recommend that R. amesiae be re-merged into R. concinnum, as R. concinnum is the older name (William, 1890; Rehder and Ernest, 1913). We also encourage the available conservation funding directed towards rescue of R. amesiae be assigned to other critically endangered Rhododendrons, particularly to those Critically Endangered (CR) plants that have been proven to be good species.

Traditional taxonomy is mainly based on morphological data. One problem with the use of morphology is that it is very difficult to capture the magnitude of intraspecific and interspecific variation (Marczewski, 2016). Such problems in Rhododendron are not uncommon. Li et al. (2018) demonstrated that the critically endangered R. protistum var. giganteum (the giant tree rhododendron) should be considered synonymous with the R. protistum var. protistum, which was evaluated to be Near Threatened (Gibbs et al., 2011). Additionally, natural hybrids with intermediate morphology between parental species may often be wrongly described as new species, after which they are designated as threatened because many hybrids are formed in very limited locations (e.g., sympatric distribution area) and numbers (particularly for early generation hybrids when reproductive barriers are still strong between parental species). For instance, the threatened plant Semiliquidambar cathayensis (Vulnerable status in the Threatened Species List of China's Higher Plants) was not a species but an F1 hybrid that originated from a natural hybridization between Altingia obovata and Liquidambar formosana (Wu et al., 2010).

Therefore, we propose that before initiating conservation actions, species delimitation should be confirmed by comparing the morphological and genetic differentiation between the target threatened species and a close relative, as was demonstrated in the present study (Li et al., 2018; Zhang et al., 2020). This applies particularly for conservation of endangered species belonging to taxonomically complex genera. The conservation issue revealed by this case study might be a general problem in conservation biology, especially for endangered taxa with few diagnostic characters (Marczewski et al., 2016; Li et al., 2019).

Author contributions

YSA drafted the manuscript, YHC and DTL collected samples, YBL and YPM designed experiment and revised the manuscript.

Data availability

Paired-end data are deposited in the NCBI SRA database under the BioProject accession PRJNA832044.

Declaration of competing interest

All authors declare that they have no competing interests and personal relationships and agree on the contents of the paper.

Acknowledgments

Thanks to Dr. Jane Marczewski from University of Edinburgh for many constructive suggestions. This study was supported by the Biodiversity Survey and Assessment Project of the Ministry of Ecology and Environment, China (2019HJ2096001006), CAS "Light of West China" Program (to Y. Ma), and the Ten Thousand Talent Program of Yunnan Province (YNWR-QNBJ-2018-174).

Appendix A. Supplementary data

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

References
Cao, Y.R., Ma, Y.P., Li, Z.H., et al., 2022. Genetic diversity and population structure of Rhododendron longipedicellatum, an endangered species. Trop. Conserv. Sci., 15: 19400829221078112.
Chamberlain, D.F., Hyam, R., Argent, G., et al., 1996. The Genus Rhododendron: its Classification and Synonymy. Royal Botanic Garden Edinburgh, Edinburgh, p. 181.
Danecek, P., Bonfield, J.K., Liddle, J., et al., 2021. Twelve years of SAMtools and BCFtools. GigaScience, 10: giab008. DOI:10.1093/gigascience/giab008
Doyle, J.J., Doyle, J.L., 1990. Isolation of plant DNA from fresh tissue. Focus, 12: 13-15.
Earl, D.A., vonHoldt, B.M., 2012. Structure HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Conserv. Genet. Resour., 4: 359-361. DOI:10.1007/s12686-011-9548-7
Gibbs, D., Chamberlain, D., Argent, G., 2011. The Red List of Rhododendrons. Botanic Gardens Conservation International, p. 16.
Li, H., Durbin, R., 2009. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25: 1754-1760. DOI:10.1093/bioinformatics/btp324
Li, H., Handsaker, B., Wysoker, A., et al., 2009. The sequence alignment/map format and SAMtools. Bioinformatics, 25: 2078-2079. DOI:10.1093/bioinformatics/btp352
Li, S.H., Sun, W.B., Ma, Y.P., 2018. Does the giant tree rhododendron need conversation priority?. Global Ecol. Conserv., 15: e00421. DOI:10.1016/j.gecco.2018.e00421
Liu, D.T., Zhang, L., Wang, J.H., et al., 2020. Conservation genomics of a threatened rhododendron: contrasting patterns of population structure revealed from neutral and selected SNPs. Front. Genet., 11: 757. DOI:10.1201/9781351187435-95
Ma, H., Liu, Y.B., Liu, D.T., et al., 2021. Chromosome-level genome assembly and population genetic analysis of a critically endangered rhododendron provide insights into its conservation. Plant J., 107: 1533-1545. DOI:10.1111/tpj.15399
Ma, Y.P., Chen, G., Grumbine, R.E., et al., 2013. Conserving plant species with extremely small populations (PSESP) in China. Biodivers. Conserv., 22: 803-809. DOI:10.1007/s10531-013-0434-3
Marczewski, T., Ma, Y.P., Zhang, X.M., et al., 2016. Why is population information crucial for taxonomy? A case study involving a hybrid swarm and related varieties. AOB Plants, 8: plw070. DOI:10.1093/aobpla/plw070
Peterson, B.K., Weber, J.N., Kay, E.H., et al., 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., Donnelly, P., 2000. Inference of population structure using multilocus genotype data. Genetics, 155: 945-959. DOI:10.1093/genetics/155.2.945
Qin, H.N., Yang, T., Dong, S.Y., et al., 2017. Threatened species list of China's higher plants. Biodivers. Sci., 25: 696-744. DOI:10.17520/biods.2017144
Rehder, A.W., Ernest, H., 1913. Rhododendron amesiae. Plantae Wilsonianae, 1: 523.
Sun, W.B., 2013. Conservation of Plant Species with Extremely Small Populations in Yunnan: Practice and Exploration. Yunnan Science and Technology Press, pp. 1-100.
William, B.H., 1890. Rhododendron concinnum. J. Linn. Soci. Bot. London xxvi: 21.
Wu, W., Zhou, R.C., Huang, Y.L., et al., 2010. Molecular evidence for natural intergeneric hybridization between Liquidambar and Altingia. J. Plant Res., 123: 231-239. DOI:10.1007/s10265-009-0275-z
Zhang, X.J., Liu, X.F., Liu, D.T., et al., 2021. Genetic diversity and structure of Rhododendron meddianum, a plant species with extremely small populations. Plant Divers., 43: 472-479. DOI:10.1016/j.pld.2021.05.005
Zhang, L., Xu, P.W., Cai, Y.F., et al., 2017. The draft genome assembly of Rhododendron delavayi Franch. var. delavayi. GigaScience, 6: 1-11.
Zhang, X.M., Qin, H.T., Xie, W.J., et al., 2020. Comparative population genetic analyses suggest hybrid origin of Rhododendron pubicostatum, an endangered plant species with extremely small populations endemic to Yunnan, China. Plant Divers., 42: 312-318. DOI:10.1016/j.pld.2020.06.012