1. Introduction

Gene flow between sympatrically occurring plants is very common when their reproductive isolation has not been completely formed (Abbott et al., 2013; Ma et al., 2010a, 2010b; Marczewski et al., 2015; Milne et al., 1999, 2003; Zha et al., 2008,2010). Gene flow is generally associated with successful pollen flow between species that have overlapping flowering times, mediated by wind in anemophilous plants, or by shared pollinators. Importantly, gene flow from a common species may accelerate the threat of extinction for rare species. For instance, pollen flow from other species can cause pollen and/or ovule discounting, and gene flow can lead to genetic swamping, by which partially fertile and viable hybrids replace threatened genotypes (Ma et al., 2019). Due to weak reproductive barriers among close relatives, hybridization has been frequently reported in sympatric distribution Rhododendron species (i.e. Milne and Abbott, 2008; Ma et al., 2010b; Zha et al., 2008; Yan et al., 2015). Indeed, natural hybridization is typically considered deleterious in plant conservation genetics, especially when rare species hybridize with widespread species of related taxa (e.g., Gilman and Behm, 2011; Burgess and Husband, 2006; Ma et al., 2016).

Rhododendron pubicostatum T. L. Ming is listed as an endangered species in the Red List of Rhododendrons (Gibbs et al., 2011). Only one population (Jiaozi Mountain) of Rh. pubicostatum was recorded from specimen records, and the total mature individuals of the species are less than 5000 based on our present investigations. It is consistent with the concept of Plant Species with Extremely Small Populations (PSESP), which focuses on species that face an elevated risk of extinction, characterized by small remaining populations in restricted habitats and exposure to serious human disturbances (Ma et al., 2013; Sun, 2016; Sun et al., 2019a, 2019b). During our fieldwork in 2017 , we found Rh. pubicostatum, Rh. bureavii Franchet and Rh. sikangense var. exquisitum T. L. Ming sympatric distribution in Jiaozi Mountain. The overlapping flowering periods of the three taxa and intermediate morphological characters of Rh. pubicostatum compared with Rh. bureavii and Rh. sikangense var. exquisitum suggested that Rh. pubicostatum might be hybrid origin. For instance, sparse hairs underside the leaf can be found in Rh. pubicostatum whereas nearly glandular in Rh. sikangense var. exquisitum and dense hairs in Rh. bureavii. Also Rh. pubicostatumhad numbers with 5-7 per inflorescence comparing with Rh. bureavii having 10- 20 flowers and Rh. sikangense var. exquisitum 3-6 flowers per inflorescence. In addition, there has been little research on the population structure and genetic diversity of Rh. pubicostatum, although this information is important for the development of meaningful conservation management strategies for this species.

In this study, we developed highly polymorphic microsatellites (SSRs) in Rh. pubicostatum and cross-amplified these SSRs in two sympatric Rhododendron species. By comparing genetic diversity and genetic compositions of the three taxa for clarifying the species status of Rh. pubicostatum. Our aim was to develop meaningful conservation strategies for Rh. pubicostatum.

2. Materials and methods 2.1. Plant materials and study site

Rhododendron pubicostatum (subsect. Taliensia), was originally collected by P. I. Mao in 1952 from the Wumeng Mountains in Luquan County, Yunnan province at an altitude of 3650 m, but was not formally described as a new species until 1981 (Ming, 1981). There are no reports about Rh. pubicostatum until it was recalled as an endemic species in 2005 . Jiaozi mountain nature reserve (N 26.09°, E 102.84°) in the Wumeng Mountains is the only population of Rh. pubicostatum distributed, and the reserve is characterized by steep terrain, steep mountain and suspended altitude difference. There are abundant Rhododendron resources (e.g., Rh. pubicostatum, Rh. bureavii and Rh. sikangense var. exquisitum) in this area. Due the morphological traits of Rh. pubicostatum have the intermediate traits between Rh. bureavii and Rh. sikangense var. exquisitum. We suspect that Rh. pubicostatum is hybrid between Rh. sikangense var. exquisitum (subsect. Maculifera) and Rh. bureavii (subsect. Taliensia). In order to verify Rh. pubicostatum is a natural hybridization origin, we collected molecular experimental materials of Rh. pubicostatum, Rh. bureavii and Rh. sikangense var. exquisitum for the study.

In our study site, these two species are closely related both to each other and to Rh. pubicostatum. Although Rh. sikangense var. exquisitum in Jiaozi mountain nature reserve with the lowest elevation and Rh. bureavii in the highest, but they have much wider distribution than Rh. pubicostatum in the reserve. A total of 24 individuals randomly sampled from Rh. pubicostatum used to develop microsatellite markers were collected from the study site in Jiaozi mountain nature reserve, to assess polymorphisms of the developed microsatellite markers. The feasibility of the developed microsatellite markers was also assessed in 24 individuals from Rh. bureavii and Rh. sikangense var. exquisitum, respectively. Healthy young leaves were collected from randomly selected plants and immediately stored in silica-gel until DNA extraction. Voucher specimens were deposited in Yunnan Key Laboratory for Integrative Conservation of Plant Species with Extremely Small Populations, Kunming Institute of Botany under Chinese Academy of Sciences.

2.2. SSR primer development and amplification of DNA from the three sympatric Rhododendron species

Genomic DNA was extracted from dried leaves by using the cetyltrimethylammonium bromide (CTAB) method as described by Doyle (1991). Extracted DNA from one sample each of the three species was sequenced in individual lanes on a MiSeq Benchtop sequencer (Illumina, Inc., San Diego, California, USA) using 2 × 250-bp read length. The expected products and primer sizes ranged from 100 to 280 bp and 18 to 23 bp, respectively. A summary of the microsatellite markers we developed are presented in Table 1.

PCR amplifications were performed in a total reaction volume of 20 m L containing 25e50 ng of template DNA, 1 m L of each primer (0.6 mM), and 10 m L 2 × Taq PCR MasterMix [Tiangen (Tiangen Biotech, Beijing, China); 3 mM MgCl₂, 100 mM KCl, 0.5 mM of each dNTP, 20 mM TriseHCl (pH 8.3), 0.1 U Taq polymerase]. PCR amplification was performed under the following conditions: an initial denaturation step for 3 min at 95 °C, followed by 33 cycles of denaturation at 94 °C for 30 s, annealing at 60 °C or 61 °C according to the Ta value (optimized for each locus; Table 1) for 30 s, extension at 72 °C for 45 s, and a final extension at 72 °C for 7 min. PCR amplification was carried out on a Veriti® 96-Well Thermal Cycler (Applied Biosystems, Foster City, USA). The PCR products were separated and visualized using an ABI3730 XL automated DNA fluorescent sequencer (Applied Biosystems, Foster City, CA, USA). Only products that were successfully amplified to an appropriate size and gave clear peaks were used for further polymorphism verification. From the 72 individuals studied, we validated a total of 12 microsatellite loci as polymorphic.

2.3. Comparative analysis of population genetics of three Rhododendron species

We visualized and quantified the SSR data using GENEMARKER v2.2.0 (Soft Genetics, State College, PA, USA). Based on the microsatellite data from all the sampled individuals, we then calculated certain population genetic diversity parameters, including number of alleles (N_A), effective number of alleles (N_E), Shannon information index (I), observed heterozygosity (H_O), and expected heterozygosity (H_E) using GenAlEx, version 6.502 (Peakall and Smouse, 2006). Deviation from HardyeWeinberg equilibrium (HWE) was calculated using the heterozygosity-based method (Nei's GIS) with 1000 permutations for each locus using the web version of GENEPOP, v4.0.10 (Raymond, 1995). Linkage disequilibrium (LD) for all loci pairs with probability adjusted using Bonferroni correction was also assessed using GENEPOP.

To further explore the genetic relationships among sampled individuals, principal coordinates analysis (PCoA) was performed in GenAlEx, version 6.502 (Peakall and Smouse, 2006). The pairwise genetic distances were based on Nei's standard genetic distance (Nei, 1978). Scatterplots were visualized in two-dimensional space for the first and second, as well as the first and third, principal coordinates.

To estimate the extent of admixture in the 72 study individuals, an analysis was conducted in STRUCTURE v.2.3.4 (Pritchard et al., 2000) using the Bayesian model-based clustering algorithm. Ten independent runs with a K value that varied from 1 to 10 were performed with a burn-in period of 1 × 10⁴ steps followed by 1 × 10⁵ Markov Chain Monte Carlo (MCMC) iterations, using an admixture model assuming correlated-allele frequencies. The bestfit number (K) of genetic clusters was estimated using STRUCTURE HARVESTER based on the maximum △K value (Earl, 2012; Evanno et al., 2005).

3. Results 3.1. Characteristics of SSRs

A total of 66, 229 microsatellites were identified from sequences obtained from the three study individuals. Dinucleotides were the most abundant repeat motifs (54.7%, 36, 241 sequences), followed by tetranucleotides (22.3%, 14, 791 sequences), pentanucleotide (11.5%, 7614 sequences), hexanucleotide (5.9%, 3937 sequences) and trinucleotides (5.6%, 3646 sequences). The dominant dinucleotide motifs were GA/CT repeats (37.7%) and the least common repeat was GC/CG (0.3%). We designed 100 primers and used them to evaluate PCR amplification efficiency and polymorphs: 17 primers did not generate clear microsatellite peaks or failed to amplify DNA from any of the three species; 76, 71 and 64 primers successfully amplified DNA samples from Rh. pubicostatum, Rh. bureavii and Rh. sikangense var. exquisitum, respectively. Only 58 primers generated clear microsatellite peaks from samples of all three species. Therefore, 58 primers were used for further evaluation, of which 12 showed polymorphisms. Primer sequences are given in Table 1.

Of the 12 polymorphic SSR loci used in our study, four were dinucleotides, whereas the remainder were trinucleotides. The genetic diversity parameters of these three species are shown in Table 2. The numbers of alleles per locus ranged from 2 to 11 (mean 6.583). In Rh. pubicostatum, the observed and expected heterozygosities ranged from 0.125 to 0.958 (mean 0.517) and 0.119 to 0.821 (mean 0.565), respectively. Three polymorphic microsatellite loci (RH16369, RH175950 and RH226095) significantly deviated from HardyeWeinberg equilibrium in Rh. pubicostatum. In total, RH119455, RH150193, RH16369, RH175950, RH178203, RH195227, RH19935, RH202706, and RH226095 polymorphic microsatellite loci were detected to be deviated from HardyeWeinberg equilibrium (Table 2). Only three polymorphic microsatellite loci (RH185982, RH19511 and RH209912) did not deviate from HardyeWeinberg equilibrium in three study species, indicating the possibility of null alleles or an insufficient number of samples (Table 2). After Bonferroni correction, no significant linkage disequilibrium was observed for any pair of loci.

3.2. Population genetics of the three Rhododendron species

Population genetic analysis showed that Rh. pubicostatum had the slightly higher number of effective alleles (3.326) than the two other species. Rh. pubicostatum observed, expected, and unbiased expected heterozygosities were 0.486, 0.566 and 0.576, respectively. Pairwise F_st values of the three species were highest (0.231) between Rh. bureavii and Rh. sikangense var. exquisitum, and lowest (0.081) between Rh. pubicostatum and Rh. sikangense var. exquisitum (Table 3 and 4) (Fig. 1).

The two-dimensional PCoA plot showed that the first principal coordinate accounted for 37.26% of total variation and separates Rh. sikangense var. exquisitum from Rh. bureavii. Rh. pubicostatum samples are intermediate between those of Rh. bureavii and Rh. sikangense var. exquisitum (Fig. 2). The second and third principal coordinates (17.26% and 13.15% of total variation) did not provide sufficient information for distinguishing the three taxa.

Based on STRUCTURE analysis, the value of △K was clearly highest (13.336) when K=2. Following ten independent STRUCTURE runs with K=2, 23 individuals morphologically identified as Rh. sikangense var. exquisitum were assigned to one cluster with high probability (all q > 90%). The remaining single individual (no.1) seemed to be backcrossed Rh. sikangense var. exquisitum with some genetic materials from Rh. bureavii. Unlike Rh. sikangense var. exquisitum, more individuals (10/24) of Rh. bureavii showed evidence of admixture using the threshold value of q > 90%. These Rh. bureavii-like individuals were probably intermediate or backcrossed genotypes. The Rh. pubicostatum individuals have genetic material from both Rh. bureavii and Rh. sikangense var. exquisitum, with proportions of Rh. bureavii genetic material ranging from 40% to 60%.

4. Discussion 4.1. High efficiency of microsatellite development using MiSeq

Although SSRs have been widely used in population genetics for numerous species, there are far fewer SSR polymorphic sites than in Miseq technology (Li et al., 2002; Glover et al., 2010; Guichoux et al., 2011). In the present study, multiple microsatellite loci were developed by using whole genome shotgun Illumina MiSeq sequencing, which is a useful tool for isolating new genetic markers in species with poor genomic resources and large genome sizes, while also circumventing the technical challenges of more traditional microsatellite enrichment protocols (Ekblom and Galindo, 2011; Ritchie et al., 2016). In high-throughput sequencing technology, RNA-seq technology is relatively lowcost and can result in large numbers of sites, but these tend to be non-neutral (Wang et al., 2009). Therefore, RNA-seq in not appropriate for obtaining the large numbers of neutral loci necessary for the genetic analysis of population structures. While RAD-seq technology and enzymatic digestion-based sequencing genotyping (GBS) can obtain more molecular marker sites than RNA-seq technology, the numbers of neutral sites revealed by these two technologies is still not comparable to that of Miseq technology (Mesak et al., 2014; Suchan et al., 2016). All 12 SSR primer pairs developed in this study gave an H_E > 0.01 across all three Rhododendron taxa studied, suggesting that the sequences they amplify are polymorphic across these taxa (Ott, 1992), which is consistent with the finding that the three species have 100 percentage of polymorphic loci (PPL). Additionally, in Rh. pubicostatum the primers RH119455, RH150193, RH16369, RH178203, RH19935 and RH226095 resulted in amplification products showing an H_E > 0.70 (Table 2), suggesting that these primers amplify highly polymorphic sequences (Ott, 1992). Miseq technology is indeed a relatively efficient and low-cost option for the development of a large number of neutral microsatellites.

4.2. Rhododendron pubicostatum maintains moderate levels of high genetic diversity

Many factors affect the genetic diversity of species, including geographic distribution, breeding system, and population size (Hamrick and Godt, 1989). Species with larger ranges tend to have higher genetic diversity than those with limited ranges (Hamrick and Godt, 1989 Hamrick et al., 1992). The widespread species Rh. simsii Planch (H^E=0.754), Rh. branchyacarpum G.Don (H_E=0.815) and Rh. ripense Makino (H_E=0.8) have higher genetic diversity than the narrowly endemic species Rh. pseudochrysanthum Hay. (H_E=0.42) and Rh. oldhamii Maximowicz (H_E=0.284) (Tan et al., 2009; Hirao, 2010; Kondo et al., 2009; Chen et al., 2014; Hsieh et al., 2013). In our study, Rh. sikangense var. exquisitum (H_E=0.549) and Rh. bureavii (H_E=0.508) showed lower levels of genetic diversity than were reported from the widespread species of Rh. simsii, Rh. branchyacarpum and Rh. ripense, which were analyzed with SSRs (Tan et al., 2009; Hirao, 2010; Kondo et al., 2009). Rh. bureavii and Rh. sikangense var. exquisitum showed quite low levels of genetic diversity for widespread species, but these figures are perhaps underestimates because we sampled only a single population of each taxon. Rh. pubicostatum (H_E=0.565) showed higher genetic diversity than some narrowly endemic species have, including Rh. bureavii and Rh. oldhamii (Chen et al., 2014; Hsieh et al., 2013). Nybom (2004) calculated that the average diversity of endemic species when assessed using SSRs was 0.42, whereas that of widespread species was 0.62. From these numbers, the Rh. pubicostatum in our study would fall between endemic and widespread species, although the taxon is described as a narrowly endemic species known only from a single locality. But, in our study, Rh. pubicostatum showed higher levels of genetic diversity than Rh. bureavii and Rh. sikangense var. exquisitum. This is consistent with previous findings that the genetic diversity of natural hybrid individuals of 'Rh. ⅹ duclouxii' (H_E=0.812) is higher than that of its parents Rh. spiciferum Franch. (H_E=0.765) and Rh. spinuliferum Franch. (H_E=0.728) (Yan et al., 2017). Hybridization or introgression between sympatric populations of species have been shown to increase genetic diversity (Dobes et al., 2004). In other words, these findings suggest that Rh. pubicostatum individuals are of a hybrid origin.

In the three sympatric taxa in our study, values of Ho are lower than H_E, suggesting the occurrence of inbreeding depression. Although many plants appear to minimize self-pollination, resulting in deposition of more outcrossed pollen (Jong et al., 1993 ), it is possible that geitonogamy plays a major role in Rhododendron. In Rh. cyanocarpum and other Rhododendron species with known breeding systems (Ma et al., 2015; Li et al., 2018), geitogamous pollination is likely to occur frequently, due to large floral displays (comprising thousands of flowers) from a single individual and pollinator behavior, where the pollinators move between nearby flowers. Inbreeding depression might result from geitonogamy, a pollination strategy frequently reported in several Rhododendron (e.g., Rh. cynocarpun and Rh. longipedicellatum) (Ma et al., 2015; Liu et al., 2019). It has been proposed that geitonogamy is favored by the foraging behavior of pollinators, as most single Rhododendron individuals often display a large number of flowers and pollinators usually minimize inter-flower travel by preferentially foraging in adjacent flowers (Ma et al., 2015).

4.3. The hybrid origin of Rhododendron pubicostatum and implications for conservation

Both predominant components analysis (PCoA) and Bayes STRUCTURE analysis showed distinct genetic clusters formed by the Rh. bureavii and Rh. sikangense var. exquisitum plants sampled in our study, although gene flow clearly occurs between the two taxa. However, the Rh. pubicostatum individuals showed admixed ancestry, having characteristics of both Rh. bureavii and Rh. sikangense var. exquisitum, suggesting that they are hybrids of these two taxa. This provides further evidence of natural hybridization in Rhododendron, which has been frequently reported (e.g., Kron et al., 1993; Milne et al., 1999; Zhang et al., 2007; Zha et al., 2008; Ma et al., 2010b; Yan et al., 2015). Furthermore, the sympatric distribution, close relationship and overlapping flowering times of Rh. bureavii and Rh. sikangense var. exquisitum all provide opportunities for hybridization. Moreover, habitat disturbance has also been suggested to promote hybridization (Bleeker and Hurka, 2001; Hasselman et al., 2014). Our study area is a tourism hotspot, especially in April and May, when Rhododendrons flower. The construction of roads in this area may result in habitats where hybrid plants may outcompete parent taxa. However, unlike other known hybrid zones in Rhododendron where the parental species are dominant and only limited numbers of hybrids are observed (e.g., Rh. cyanocarpum < 20), Rh. pubicostatum is dominant in our study site, and the two parental species mainly occupy marginal areas. This hybrid zone is, therefore, unimodal, which implies that it may have been maintained over a long period, possibly before the development of Jiaozi Mountain for tourism.

If the Rh. pubicostatum individuals are of hybrid origin, this would also explain the higher genetic diversity observed in these specimens than that of the two parental species. Our preliminary conclusion that Rh. pubicostatum is of hybrid origin may not provide insight into whether Rh. pubicostatum requires protection. If strong evidence supports very early generations of Rh. pubicostatum and the existence offrequent introgressionwith parents, Rh.pubicostatumshouldnotbe a conservation priority. However, if evolutionary potential (e.g., differentiated local adaptation, formation of strong reproductive barriers between Rh. pubicostatum and its parents) can be detected in Rh. pubicostatum, it should be made a conservation priority (Oro et al., 2004; Genovart et al., 2007; Genovart, 2009). Therefore, addressing the conservation status of Rh. pubicostatum requires further study by employing more markers, samples as well as detailed examination of reproductive barriers among the three taxa.

Author's contributions

Zhang X M performed the experiment and wrote the paper; Qin H T and Xie W J performed partly statistical analysis; Ma Y P and Sun W B designed the experiment and revised the paper.

Declaration of competing interest

None.

Acknowledgements

We are grateful to Junbo Yang and Zhirong Zhang to providing technical assistance of lab work and data analysis. Laboratory work was performed at the Laboratory of Molecular Biology at the Germplasm Bank of Wild Species, Kunming Institute of Botany, Chinese Academy of Sciences. This work was supported by Science and Technology Basic Resources Investigation Program of China (Grant No. 2017FY100100), the Second Tibetan Plateau Scientific Expedition and Research (STEP) Program (Grant No. 2019QZKK0502), the Young Academic and Technical Leader Raising Foundation of Yunnan Province (No. 2018HB066), Yunnan Innovation Team Program for Conservation and Utilization of PSESP (Plant Species with Extremely Small Populations) (Grant No.2019HC015) and Applied Basic Research Project of Yunnan Province (Grant No. 2018BB010).