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.

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 S_TACKS 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 S_TACKS, 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 (H_e), observed heterozygosity (H_o), inbreeding coefficient (F_IS) and F statistics (F_ST) were calculated using the populations program in the S_TACKS pipeline. Population genetic structure analyses based on SNPs was performed using PCA and S_TRUCTURE v.2.3.4 (Pritchard et al., 2000). The optimal K was chosen using the delta-K method implemented in S_TRUCTURE H_ARVESTER (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 S_TACKS 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 (H_o) ranged from 0.037 to 0.049 (average 0.041) and expected heterozygosity (H_e) 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 (F_IS: 0.147–0.295, average 0.215) (Table 1).

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

Bayesian clustering analysis was performed in S_TRUCTURE, 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.

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 F₁ 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.