Natural hybridization among three Rhododendron species (Ericaceae) revealed by morphological and genomic evidence
BMC Plant Biology volume 21, Article number: 529 (2021)
Natural hybridization can influence the adaptive response to selection and accelerate species diversification. Understanding the composition and structure of hybrid zones may elucidate patterns of hybridization processes that are important to the formation and maintenance of species, especially for taxa that have experienced rapidly adaptive radiation. Here, we used morphological traits, ddRAD-seq and plastid DNA sequence data to investigate the structure of a Rhododendron hybrid zone and uncover the hybridization patterns among three sympatric and closely related species.
Our results show that the hybrid zone is complex, where bi-directional hybridization takes place among the three sympatric parental species: R. spinuliferum, R. scabrifolium, and R. spiciferum. Hybrids between R. spinuliferum and R. spiciferum (R. ×duclouxii) comprise multiple hybrid classes and a high proportion of F1 generation hybrids, while a novel hybrid taxon between R. spinuliferum and R. scabrifolium dominated the F2 generation, but no backcross individuals were detected. The hybrid zone showed basically coincident patterns of population structure between genomic and morphological data.
Natural hybridization exists among the three Rhododendron species in the hybrid zone, although patterns of hybrid formation vary between hybrid taxa, which may result in different evolutionary outcomes. This study represents a unique opportunity to dissect the ecological and evolutionary mechanisms associated with adaptive radiation of Rhododendron species in a biodiversity hotspot.
Natural hybridization is frequent in plants and plays a crucial role in the formation and maintenance of species. In recent decades, natural hybridization has been well-documented for numerous herbaceous (e.g., Helianthus [1, 2], Iris , Senecio , Viola , Gagea , Brassica , Mimulus ) and some woody plant taxa (e.g., Pinus [9, 10], Ostryopsis [11, 12], Rhododendron [13,14,15,16]). For the most part, studies reveal that when two or more closely related species are in sympatry, hybridization frequently occurs and natural hybrid zones are likely to form (e.g. [15, 17, 18]).
Hybrid zone formation is ultimately a response to selection and dispersal mechanisms acting on hybrids and parental species, where parental genomes are combined and functionally selected . In general, there are three main types of hybrid zones that have been documented in plant taxa; these include: tension zones, bounded hybrid superiority zones and mosaic hybrid zones [20, 21]. Within tension zones, hybrids are continually formed, but are selected against due to their low fitness relative to parental species (e.g., Populus fremontii × P. angustifolia [22, 23]; Senecio chrysanthemifolius × S. aethnensis [4, 24]). In bounded hybrid superiority zones, hybrids show higher fitness than both parental species in intermediate habitats, but lower fitness in parental habitats (e.g., Picea glauca × P. engelmannii ; Artemisia tridentata ssp. tridentata × A. tridentata ssp. vaseyana [26,27,28]). For mosaic hybrid zones, hybrids have higher fitness than parental taxa across a patchwork of multiple local habitats where the parental species overlap in range (e.g., Aquilegia formosa × A. pubescens ; Senecio ovatus × S. hercynicus ). Ultimately, the evolutionary outcome of natural hybridization between plant taxa is likely to depend on the type of hybrid zone that is formed, where understanding the genetic composition and structure of hybrid zones is an important first step in revealing the evolutionary mechanisms (e.g., selection, dispersal, reproductive isolation) and outcomes associated with natural plant hybridization [17, 31,32,33].
The genetic composition and structure of plant hybrid zones is often complex (i.e., multiple generations of filial and backcross hybrids), and depending on the strength of pre- and post-reproductive isolation among the parental taxa and their hybrids, hybrid taxa are able to persist, and can even form new species [34, 35]. In many cases, however, the formation and persistence of hybrids depends on the direction and symmetry of hybridization among parental taxa. For example, numerous studies have detected hybridization patterns among parental taxa that are asymmetric [8, 17, 36,37,38]. In such a case, one likely outcome is asymmetric introgression, where the genome of one parental taxon can become incorporated into the genetic background of another, due in part, to the preferential backcross of hybrids and the numerical asymmetries associated with parental taxa and their hybrids. The detection of such patterns can have profound implications for understanding the formation and persistence of potential hybrid species within a hybrid zone; hybrids, and even some parental taxa, that are not reproductively isolated, are unlikely to persist under such a scenario [39,40,41,42]. Alternatively, adaptive introgression can play a central role in the formation and persistence of novel hybrid taxa, and may even lead to the formation of new species. Under this scenario, habitat selection can favor hybrids that incorporate parental genes, which enable them to colonize and persist in novel habitats associated with the hybrid zone (i.e., hybridized habitats, which may be due, in part, to human disturbance ). Here, the adaptive introgression of parental traits in hybrids can buffer the deleterious effects of asymmetrical introgression, resulting in the formation and persistence of a complex hybrid zone that is composed of both filial and backcross generations [41, 44,45,46]. Detection of these unique gene × environment interactions (i.e., the composition of hybrids), in addition to an assessment of hybrid zone structure, may reflect plausible pathways for the formation and persistence of novel hybrid species [21, 40, 47].
Rhododendron L. (Ericaceae) is a large and taxonomically complex genus of woody plants that includes more than one thousand species, where more than half of which are endemic to China [48,49,50], and many species are of ornamental value. The species diversity of Rhododendron is rich in the Mountains of Southwest China, which is a global biodiversity hotspot [51, 52]. The morphology of Rhododendron species varies considerably, due in part, to substantial adaptive radiation and reticulate evolution of species within the genus . In addition, natural hybridization occurs frequent within Rhododendron [13, 16, 54,55,56,57]. Previous studies have verified that R. ×duclouxii H. Lévl., a hybrid taxon between R. spiciferum Franch. and R. spinuliferum Franch. that has intermediate leaf and flower morphology , exists in at least 15 natural hybrid zones in Yunnan, southwest China [16, 58]. A third species, R. scabrifolium Franch., also exists within some of these hybrid zones and it remains unknown as to the extent that this species contributes to hybrid formation and hybrid zone structure. It’s worth noting that all three species, in addition to R. ×duclouxii, are diploid and that R. scabrifolium and R. spinuliferum are sister taxa .
Recently, a putative novel hybrid taxon is suspected to exist in a region where R. spiciferum, R. spinuliferum and R. scabrifolium have overlapping ranges. These novel hybrids manifest a suite of morphological traits that not only differ from the three parental taxa at the site, but also differ from the known hybrid taxon R. ×duclouxii (e.g., flower size, flower color, corolla shape; Fig. 1) and may represent a new hybrid taxon. Here, we use genomic (ddRAD-seq) and morphological data to assess the genetic composition and structure of this Rhododendron hybrid zone in Yunnan Province, SW China (Fig. 2). We aim to address the following questions: (1) What is the origin of the novel hybrid taxon? (2) What is the genetic composition and structure of this hybrid zone? (3) Do patterns of hybridization differ between the novel hybrid taxon and R. ×duclouxii?
The distribution of variance in the PCoA based on morphological traits of 97 individuals distributed across the five taxa, was 69.7 % 10.0 % and 7.8 % for the first three axes, respectively, resulting in a cumulative value of 87.5 %. Individuals of the three parental species (R. spiciferum (SC), R. spinuliferum (SN), R. scabrifolium (SA)) form distinct clusters (Fig. 3a). For the two hybrid taxa, individuals of R. ×duclouxii (SN×SC) clustered between the three parent species, while individuals of the putative novel hybrid (SN×SA) were mainly distributed between R. scabrifolium and R. spinuliferum; hybrid individuals (R. ×duclouxii and SN×SA) showed some degree of overlap (Fig. 3a). When we considered variation among parental and hybrid individuals for the 14 quantitative traits separately, some interesting patterns were detected. Firstly, values for individuals of R. spinuliferum are larger than both R. spiciferum and R. scabrifolium for six traits (leaf length, leaf width, leaf area, flower tube length, style length, filaments length) revealing an inverted ‘V’ pattern (R. spinuliferum > hybrid > R. spiciferum/R. scabrifolium) (Fig. 3b, Fig. S1), while other traits (flower tube width, stigma width, ovary width) show a pattern that resembles an inverted ‘U’ (R. spinuliferum≈hybrid > R. spiciferum/R. scabrifolium) (Fig. 3c, Fig. S1). In addition, the flower width of SN×SA individuals were similar to that of R. scabrifolium individuals, while R. ×duclouxii individuals have a flower width that is more similar to individuals of R. spiciferum (Fig. 3d). Individuals of R. scabrifolium and R. spinuliferum had similar petiole length (Fig. 3e), and the two types of hybrid individuals had longer corolla lobes when compared to parental taxa (Fig. 3f). For some flower traits (e.g., style length, stigma width, filaments length), the magnitude of variation among SN×SA individuals was greater than those of R. ×duclouxii individuals (Fig. S1). The leaf width and pedicel length of SN×SA were significantly larger than R. ×duclouxii (Fig S1a, d). For the other quantitative traits, the differences between SN×SA and R. ×duclouxii are not significant.
The genetic structure of the hybrid zone
Approximately 40 Gb ddRAD-seq raw data was generated for all 45 individuals. After mapping to the Rhododendron williamsianum reference genome , on average, 617,833 loci (RAD tags) were assembled and the average coverage depth per locus was 27.7×. Furthermore, 30,346 SNPs were retained based on the populations pipeline analysis in Stacks. Finally, 11,011 SNPs were identified for Dataset A after Plink filtering. Admixture analysis based on the SNPs in Dataset A, showed that when K=3, the three parental species (R. spiciferum, R. spinuliferum, R. scabrifolium) formed three distinct genetic groups (Fig. 4a), and was the best-fit model. Individuals of SN×SA were a genetic admixture of R. scabrifolium and R. spinuliferum. Of which five individuals contained approximately 50 % of the genetic component from each parental species, another eight individuals contained a greater genetic proportion of R. scabrifolium than R. spinuliferum (Fig. 4a). The individuals of R. ×duclouxii revealed complex genetic profiles composed of varying amounts of R. spinuliferum and R. spiciferum genetic admixture (Fig. 4a). No R. scabrifolium × R. spiciferum hybrid individuals were detected in the hybrid zone. When K=4 and K=5, neither SN×SA nor R. ×duclouxii individuals formed a unique genetic cluster (Fig. S2).
The results of the PCoA analysis based on the SNPs Dataset A were in accordance with the results of the Admixture analysis. The first two principal axes explained 12.7 % and 10.1 % of the variance, respectively (Fig. 4b). The three parental species, and the two types of hybrid taxa (R. ×duclouxii and SN×SA), formed five respective well-differentiated clusters (Fig. 4b), which corresponded relatively close to those identified in the PCoA based on morphology. Rhododendron ×duclouxii individuals were intermediate between R. spinuliferum and R. spiciferum, while the individuals of SN×SA were intermediate between R. spinuliferum and R. scabrifolium (Fig. 4b). There were obvious genetic differences between individuals of SN×SA and R. ×duclouxii; the individuals of SN×SA were distinctly separated from the two parents, while individuals of R. ×duclouxii were placed on a continuum from R. spiciferum to R. spinuliferum, although more individuals were closer to R. spiciferum (Fig. 4b), which is consistent with the result of the Admixture analysis. Furthermore, no hybrid individuals between R. scabrifolium and R. spiciferum were detected. Minimum Spanning Network (MSN) analysis showed a similar result to the PCoA, where the three species separated well and the two types of hybrids fell between their parents, although two individuals of R. ×duclouxii clustered much closer to R. spinuliferum (Fig. 4c).
Extent and direction of hybridization
NewHybrids results indicated that nearly all the hybrid individuals were divided into a specific class except for one individual of R. ×duclouxii (Fig. 4a). For individuals of SN×SA, two were assigned to F1 class, while others were assigned to F2 class, corresponding to the results of the SNPs Admixture analysis (Fig. 4a). The R. ×duclouxii individuals were assigned to multiple classes: four were F1, three were F2, one was backcross to R. spinuliferum, three were backcrosses to R. spiciferum, and one was unidentified (possible F3 or later generation). One individual of R. spiciferum was assigned as backcross with R. spiciferum, suggestive of a possible hybrid individual, corresponding well to the results of Admixture and PCoA (Fig. 4a).
The sister groups R. scabrifolium and R. spinuliferum shared the same trnL-F haplotype (R. spinuliferum haplotype) as shown in Yan et al. . In this study, we found all the individuals of R. scabrifolium and SN×SA share the R. spinuliferum trnL-F haplotype (Fig. 4a). For the R. ×duclouxii individuals, four individuals contained the R. spiciferum haplotype, while eight comprised the R. spinuliferum haplotype. One individual of R. spiciferum had the R. spinuliferum trnL-F haplotype, but the nDNA genetic structure indicated it was R. spiciferum (Fig. 4a). The discrepancy between cpDNA and nDNA for this individual is probably due to hybridization followed by repeated backcrossing with R. spiciferum.
Confirmation of a novel hybrid taxon
In this study, we screened the origin of a putative novel hybrid taxon (SN×SA) for 13 individual plants collected from our Rhododendron hybrid zone. Morphologically, the majority of these individuals were distributed between R. scabrifolium and R. spinuliferum, while a few samples were distributed between R. spiciferum and R. spinuliferum (Figs. 1e and 3a). In addition, the Bayesian genetic clustering results indicated that the genetic composition of SN×SA individuals is composed of R. scabrifolium and R. spinuliferum (Fig. 4a), and the trnL-F haplotype analysis showed they share the same maternal chloroplast haplotype as the two parental species . Given that all three parental species (R. spiciferum, R. spinuliferum R. scabrifolium) formed independent clusters based on either morphology or ddRAD data, collectively our results indicate a novel hybrid taxon (SN×SA) is present at our study site and is of R. scabrifolium × R. spinuliferum origin. This result represents a first step to understand the potential evolutionary effects of natural hybridization among Rhododendron taxa within our hybrid zone, where further study on the ecological and genetic mechanisms (e.g., fitness differences, dispersal ability, reproductive barriers) governing formation and persistence of these novel hybrid taxa is needed [21, 33, 59, 60].
Complex hybridization among three closely related species
Our morphological and genetic results indicate that the hybrid zone is a complex hybrid swarm from two types of hybrid taxa among three closely related species. Our study site is located in one of the species diversification and distribution centers of Rhododendron in China, where most of the species were formed in Pliocene and Quaternary (< 5 Ma) due to adaptive radiation and subsequent rapid evolution [49, 53, 61,62,63]. In addition, this region contains many natural hybrid zones of Rhododendron species [13, 58, 64,65,66,67,68]. At our study site, which intersects the distribution range of our three parental species [16, 58], it would appear that gene flow among R. scabrifolium, R. spiciferum and R. spinuliferum is possible, although we found no evidence for the formation of hybrids between the former two species. Our previous studies show that post-zygotic reproductive barriers between R. spinuliferum and R. spiciferum are weak , and given that R. scabrifolium is the sister taxon of R. spinuliferum , it seems reasonable that the three species can hybridize with each other, although only two types of hybrid taxa were detected at the study site.
Based on field observations, the lack of hybrids produced between R. scabrifolium and R. spiciferum may be due to strong pre-zygotic barriers to reproduction. At our study site, R. scabrifolium and R. spiciferum appear to occupy different micro-environments: R. scabrifolium grows in lush evergreen broad-leaf forest, whereas R. spiciferum prefers to grow in more open disturbed habitats (Fig. 2; Table S1). In addition, the two species were often found in different plots within the hybrid zone, which may also reduce pollination opportunities, and subsequent gene flow, between R. scabrifolium and R. spiciferum at our study site. Thus, ecological niche divergence may play an important role for our study species, and this has also been found in other species (e.g., Silene , Senecio ). Based on our observations, the flowering phenology is also different between R. scabrifolium and R. spiciferum in the hybrid zone; peak anthesis for R. scabrifolium is in mid-February, whereas R. spiciferum typically flowers mid-March. Taken together, these factors may form a strong reproductive barrier between R. scabrifolium and R. spiciferum. Although we acknowledge that another possible reason for the lack of hybrids produced between R. scabrifolium and R. spiciferum may due to sampling bias (i.e., more sampling may increase the chances of finding such hybrids), these plausible barriers to gene flow could also explain why we didn’t find hybrid individuals between R. scabrifolium and R. spiciferum at our study site. Pollination experiments and further studies are needed to assess the magnitude of pre- and post-zygotic barriers to hybrid formation among the parental taxa.
Hybrid zone structure
Results of the PCoA and ANOVA analysis of morphological traits indicated that R. ×duclouxii and SN×SA share some similar morphological characters (Figs. 1 and 3, S1), which may shed light on the genetic architecture of shared phenotypic traits among the hybrid taxa present in the hybrid zone . However, our genetic data revealed that the genetic structure of these two types of hybrids is quite different. Rhododendron ×duclouxii individuals consisted of multiple kinds of offspring, which include not only F1 individuals, but also backcrosses and F2 individuals, whereas SN×SA individuals were only composed of F1 and F2 offspring, and no backcross individuals were detected. In addition, we found more R. ×duclouxii individuals have the R. spinuliferum trnL-F sequence compared to those with R. spiciferum. This result indicates that for the formation of R. ×duclouxii hybrids is bidirectional but asymmetric between R. spinuliferum and R. spiciferum, which was also reported in our previous studies [58, 69]. For the formation of the novel hybrid taxon (SN×SA), we could not ascertain the symmetry or direction of hybridization because the sequence for the trnL-F region, and even for the whole plastid genome (unpublished data), of R. scabrifolium is same as that for R. spinuliferum . Taken together, it seems likely that our hybrid zone represents a mosaic model, where different classes and generations of hybrids are selected by different microhabitats in areas of sympatry among the parental species at our study site, although we acknowledge that more environmental data, and fitness assessments, are necessary to verify this hypothesis. While this model is common in herbaceous perennial plants, e.g. Aquilegia [29, 73], Senecio , Petunia , its occurrence in hybridizing woody plant taxa is rare and requires additional evidence.
Although determining the structure of hybrid zones is essential for discerning the mechanisms and potential evolutionary outcomes of hybridization [33, 75, 76], systems that involve three or more parental taxa can make this process particularly challenging, especially when comparing results among studies. For example, three species of Ligularia (L. duciformis, L. yunnanensis and L. cyathiceps) are known to hybridize and form two groups of hybrids, a result that is similar to our findings, but unlike our study, both groups of the Ligularia hybrid taxa are restricted to F1 with no history of gene introgression . Although in our study the two groups of hybrid taxa shared R. spinuliferum as one of their parental species, they display a different pattern in population structure than that found in Ligularia . When compared to our study, the degree of reproductive isolation between parental taxa has been shown to be weak in both this study and our previous studies [58, 69], it seems likely that, R. ×duclouxii is in close proximity, may contribute to the formation of complex hybrid swarm in the future. Other hybrids zones that include three hybridizing parental species have shown that recurrent gene flow can induce adaptive introgression [43, 77,78,79]. Our study system represents a novel opportunity to explore mechanisms of potential adaptive introgression in Rhododendron, where further empirical study may determine whether the presence of novel hybrid taxa at our study site is the result of reproductive or ecological isolation resulting in unique genetic architectures that influence the adaption of hybrids in novel habitats. Further study, with increased sampling within our hybrid zone, may enable us to uncover the ecological and evolutionary mechanisms that ultimately contribute to hybrid speciation in the genus Rhododendron.
Our results uncover the presence of a complex hybrid zone in Yunnan, SW China that is comprised of three parental Rhododendron species (R. scabrifolium, R. spinuliferum, R. spiciferum) and two types of hybrid taxa, including a novel hybrid taxon of R. scabrifolium × R. spinuliferum recognized for the first time. Although the two types of hybrid taxa found at our site share R. spinuliferum as one of the parents, hybrid genetic admixture varies among them; individuals of R. ×duclouxii are an admixture of F1, F2, and backcrosses, while the novel hybrids are dominated by F2, with a low proportion of F1 and no backcross individuals. This study represents a unique opportunity to dissect ecological and evolutionary mechanisms associated with adaptive introgression, where hybridization among three parental taxa can result in the formation of complex hybrid swarms that consist of hybrid taxa that vary in their evolutionary potentials.
Study site and plant sampling
The hybrid zone is located in Nanhua County, Chuxiong Prefecture, Yunnan Province, southwest China (Fig. 2). Three closely related species in the Rhododendron subsect. Scabrifolia (R. scabrifolium, R. spiciferum, and R. spinuliferum) and their hybrids were found in sympatry within this hybrid zone. Voucher specimens (two to three duplicates per individual) and DNA samples (healthy leaves dried immediately with silica-gel) were collected from the wild. A total of 45 individuals were sampled: nine individuals of R. spiciferum, six of R. spinuliferum, five of R. scabrifolium, 12 individuals of R. ×duclouxii, and 13 individuals of the putative novel hybrid taxon (SN×SA) (Table S1). Total genomic DNA was extracted using the modified CTAB method . All specimens were identified based on morphology, and were deposited in the herbarium of Kunming Institute of Botany (KUN), Chinese Academy of Sciences, Kunming, China.
ddRAD-seq and cpDNA sequencing
Library construction of double digest restriction-site-associated DNA sequencing (ddRAD-seq) was according to the protocol of Yang et al. . Genomic DNA of each individual was digested with two restriction enzymes (AvaII, NEB Cat#: R0153S; MspI, NEB Cat#: R0106S), and then size-selected to a range of 500-700 bp. Library sequencing was performed on an Illumina HiSeq X Ten (San Diego, CA, USA) with 150 bp paired-end reads by Cloud Health (Shanghai, China). Libraries were pooled to a target of approximately 1.2 Gb raw data per individual. To determine cpDNA haplotypes, the chloroplast DNA region trnL-F was sequenced for all 45 sampled individuals from the hybrid zone. The protocols of PCR and sequencing followed Yan et al. .
Morphometric data collection
To examine morphological differentiation among R. spiciferum, R. spinuliferum, R. scabrifolium, and the two types of hybrid taxa (R. ×duclouxii and SN×SA), we chose 41 of the 45 individuals to measure 23 morphological traits (Table S2), which include nine traits that are qualitative (calyx lobe conspicuous/inconspicuous; density of abaxial and adaxial leaf hairs; density of corolla scales; filament hair present/absent; style hair present/absent; flower color; stigma color; anther color) and 14 quantitative traits (leaf length, width, thickness, and area; petiole length; pedicel length; corolla lobe length; corolla tube width and length; flower width; style length; filament length; stigma width; ovary width). For each individual, a total of three leaves and three flowers, without obvious symptoms of pathogen or physical damage, were selected for trait measurement. Qualitative floral traits were assessed by visual inspection at the site, while quantitative floral traits were measured by vernier caliper; leaf area measurements were based on scanned images (CanonScan LiDE 220; Canon, Tokyo, Japan). In addition, the 23 morphological traits were measured for 19 individuals of R. spinuliferum, 20 individuals of R. scabrifolium, and 17 individuals of R. spiciferum from other allopatric populations (close to the hybrid zone) to supplement the dataset.
To assess the dispersion of morphological traits in multivariate trait space, Principal Coordinates Analysis (PCoA) for morphometric distance (i.e., Gower distance) was conducted in the R package ggfortify . Variations of each quantitative trait were shown in boxplots generated by the R program . ANOVA analysis was performed for each quantitative trait using SPSS v19 .
The ddRAD sequencing data was analyzed by Stacks v.2.04 . In the first step, raw reads were de-multiplexed and quality-filtered by using the pipeline process_radtags. FastQC v.0.11.4  was used to assess read quality and GC-content. All reads were then trimmed to 140 bp to remove low quality bases by using -t parameter in Stacks v.2.04. The genome of R. williamsianum  was used as a reference in the reference-based analyses pipeline. To obtain a large number of high-quality SNPs, we treated the three Rhododendron species, and the two hybrid taxa (R. ×duclouxii and SN×SA), as five different populations in the populations pipeline. Loci were filtered in two steps. Firstly, the number of SNPs were reduced based on the following parameters: write_single_snp (restrict data analysis to only the first SNP per locus), r = 0.7 (a locus must exist in at least 70 % of the individuals for each of the populations), p = 5 (a locus should be present in each of the five populations), and min-maf = 0.05 (minimum minor allele frequency required to process a nucleotide site at a locus) . After this reduction step, SNPs with data missing call rates exceeding 0.2 (geno=0.2), and a Hardy-Weinberg equilibrium exact test p-value below 0.01 (hwe=0.01), were further filtered with Plink v1.90 . All SNPs generated from these two filtering steps were used as Dataset A. To distinguish hybrids in the NewHybrids analysis, another dataset (Dataset_B) was screened by the top 400 highest species pairwise FST unlinked SNPs using the function getTopLoc() in the hybriddetective R package [88, 89].
To uncover the genetic structure of the hybrid zone, Dataset A was analyzed by a Bayesian genetic clustering approach using Admixture v.1.3.0  and the best K value (number of genetic clusters ranging from 2 to 7) was evaluated by the CV (cross-validation) value in the log file. To visualize the genetic similarities among the five taxa, Principal Coordinate Analysis (PCoA) was performed using the ape package in R . The MSN (Minimum Spanning Network) analysis for all individuals was conducted by poppr package in R  and the graphs were produced using the ggplot2 package v3.3.3 .
We used NewHybrids v.1.1 to calculate the probability of each individual belonging to a particular genotypic class (P1, P2, F1, F2, BC1, BC2) . Since we identified a priori two different types of hybrid taxa in this hybrid zone, NewHybrids analysis was based on two different SNPs Datasets (high FST): Dataset_B1 includes all individuals from R. spinuliferum, R. scabrifolium and SN×SA, and Dataset_B2 includes all individuals from R. spinuliferum, R. spiciferum and R. ×duclouxii. This analysis was run for 100,000 rounds after a 100,000 burn-in iteration. If an individual was assigned to a genotype class with a probability ≥ 0.9, this individual was regarded as belonging to that class.
To detect the trnL-F cpDNA haplotype, trnL-F sequences of all individuals within the hybrid zone were aligned in Geneious v.8.1 . The informative sites were summarized, and the haplotypes belonging to R. spiciferum and R. spinuliferum were then assigned following Yan et al. .
Availability of data and materials
All newly generated trnL-F sequences had been deposited in the NCBI with the GenBank accession numbers MZ493237 to MZ493279 (Table S1). Raw ddRAD-seq reads are available at the National Center for Biotechnology Information (NCBI) database (https://www.ncbi.nlm.nih.gov) under BioProject ID PRJNA746356 (can be viewed at https://dataview.ncbi.nlm.nih.gov/object/PRJNA746356?reviewer=vls9t8u2codbav11mo4pnugsna and will be released after publication). The morphological data are provided in Table S3.
double digest Restriction-site Associated DNA
R. spinuliferum × R. spiciferum
R. spinuliferum × R. scabrifolium
Principal Coordinate Analysis
Single Nucleotide Polymorphism
Minimum Spanning Network
Rieseberg LH, Van Fossen C, Desrochers AM. Hybrid speciation accompanied by genomic reorganization in wild sunflowers. Nature. 1995;375:13.
Rieseberg LH, Sinervo B, Linder CR, Ungerer MC, Arias DM. Role of gene interactions in hybrid speciation: evidence from ancient and experimental hybrids. Science. 1996;272:741–5.
Arnold ML. Iris nelsonii (Iridaceae): origin and genetic composition of a homoploid hybrid species. Am J Bot. 1993;80:577–83.
Brennan AC, Bridle JR, Wang AL, Hiscock SJ, Abbott RJ. Adaptation and selection in the Senecio (Asteraceae) hybrid zone on Mount Etna, Sicily. New Phytol. 2009;183:702–17.
Gil HY, Kim SC. Viola woosanensis, a recurrent spontaneous hybrid between V. ulleungdoensis and V. chaerophylloides (Violaceae) endemic to Ulleung Island, Korea. J Plant Res. 2016;129:807–22.
Peterson A, Harpke D, Levichev IG, Beisenova S, Schnittler M, Peterson J. Morphological and molecular investigations of Gagea (Liliaceae) in southeastern Kazakhstan with special reference to putative altitudinal hybrid zones. Plant Syst Evol. 2016;302:985–1007.
Mason AS, Snowdon RJ. Oilseed rape: learning about ancient and recent polyploid evolution from a recent crop species. Plant Biol. 2016;18:883–92.
Vallejo-Marín M, Cooley AM, Lee MY, Folmer M, McKain MR, Puzey JR. Strongly asymmetric hybridization barriers shape the origin of a new polyploid species and its hybrid ancestor. Am J Bot. 2016;103:1272–88.
Song BH, Wang XQ, Wang XR, Ding KY, Hong DY. Cytoplasmic composition in Pinus densata and population establishment of the diploid hybrid pine. Mol Ecol. 2003;12:2995–3001.
Gao J, Wang B, Mao JF, Ingvarsson P, Zeng QY, Wang XR. Demography and speciation history of the homoploid hybrid pine Pinus densata on the Tibetan Plateau. Mol Ecol. 2012;21:4811–27.
Liu BB, Abbott RJ, Lu ZQ, Tian B, Liu JQ. Diploid hybrid origin of Ostryopsis intermedia (Betulaceae) in the Qinghai-Tibet Plateau triggered by Quaternary climate change. Mol Ecol. 2014;23:3013–27.
Wang ZF, Jiang YZ, Bi H, Lu ZQ, Ma YZ, Yang XY, Chen NN, Tian B, Liu BB, Mao XX et al. Hybrid speciation via inheritance of alternate alleles of parental isolating genes. Mol Plant. 2021;14:208–22.
Zhang JL, Zhang CQ, Gao LM, Yang JB, Li HT. Natural hybridization origin of Rhododendron agastum (Ericaceae) in Yunnan, China: inferred from morphological and molecular evidence. J Plant Res. 2007;120:457–63.
Milne RI, Abbott RJ. Reproductive isolation among two interfertile Rhododendron species: low frequency of post-F1 hybrid genotypes in alpine hybrid zones. Mol Ecol. 2008;17:1108–21.
Ma YP, Milne RI, Zhang CQ, Yang JB. Unusual patterns of hybridization involving a narrow endemic Rhododendron species (Ericaceae) in Yunnan, China. Am J Bot. 2010;97:1749–57.
Yan LJ, Burgess KS, Milne R, Fu CN, Li DZ, Gao LM. Asymmetrical natural hybridization varies among hybrid swarms between two diploid Rhododendron species. Ann Bot. 2017;120:51–61.
Liao RL, Ma YP, Gong WC, Chen G, Sun WB, Zhou RC, Marczewski T. Natural hybridization and asymmetric introgression at the distribution margin of two Buddleja species with a large overlap. BMC Plant Biol. 2015;15:146.
Zhang NN, Ma YP, Folk RA, Yu JJ, Pan YZ, Gong X. Maintenance of species boundaries in three sympatric Ligularia (Senecioneae, Asteraceae) species. J Integr Plant Biol. 2018;60:986–99.
Zieliński P, Dudek K, Arntzen JW, Palomar G, Niedzicka M, Fijarczyk A, Liana M, Cogălniceanu D, Babik W. Differential introgression across newt hybrid zones: evidence from replicated transects. Mol Ecol. 2019;28:4811–24.
Arnold ML, Bulger MR, Burke JM, Hempel AL, Williams JH. Natural hybridization: how low can you go and still be important? Ecology. 1999;80:371–81.
Abbott RJ, Brennan AC. Altitudinal gradients, plant hybrid zones and evolutionary novelty. Philos Trans R Soc B-Biol Sci. 2014;369:20130346.
Whitham TG. Plant hybrid zones as sinks for pests. Science. 1989;244:1490–3.
Sshweitzer JA, Martinsen GD, Whitham TG. Cottonwood hybrids gain fitness of both parents: a mechanism for their long-term persistence? Am J Bot. 2002;89:981–90.
James JK, Abbott RJ. Recent, allopatric, homoploid hybrid speciation: the origin of Senecio squalidus [Asteraceae], in the British isles from a hybrid zone on mount Etna. Sicily Evolution. 2005;59:2533–47.
De La Torre AR, Wang TL, Jaquish B, Aitken SN. Adaptation and exogenous selection in a Picea glauca x Picea engelmannii hybrid zone: implications for forest management under climate change. New Phytol. 2014;201:687–99.
Freeman DC, Turner WA, McArthur ED, Graham JH. Characterization of a narrow hybrid zone between two subspecies of big sagebrush (Artemisia tridentata: Asteraceae). 1991;78:805–15.
Wang H, Mcarthur ED, Sanderson SC, Graham JH, Freeman DC. Narrow hybrid zone between two subspecies of big sagerush (Artemisia Tridentata: Asteraceae). IV. reciprocal transplant experiments. Evolution. 1997;51:95–102.
Miglia KJ, Mcarthur ED, Redman RS, Rodrigue RJ, Zak JC, Freeman DC. Genotype, soil type, and locale effects on reciprocal transplant vigor, endophyte growth, and microbial functopnal diversity of a narrow sagerush hybrid zone in Salt Creek Canyon,Utah. Am J Bot. 2007;94:425–36.
Hodges SA, Armold ML. Floral and ecological isolation between Aquilegia formosa and Aquilegia pubescens. Proc Nati Acad Sci USA 1994;91:2493–6.
Raudnitschka D, Hensen I, Oberprieler C. Introgressive hybridization of Senecio hercynicus and S. ovatus (Compositae, Senecioneae) along an altitudinal gradient in Harz National Park (Germany). Syst Biodivers. 2007;5:333-44.
Arnold ML. Natural hybridization and evolution. Place:Oxford University Press;1997.
Rieseberg LH. Hybrid origins of plant species. Annu Rev Ecol Syst. 1997;28:359–89.
Abbott RJ, Barton NH, Good JM. Genomics of hybridization and its evolutionary consequences. Mol Ecol. 2016;25:2325–32.
Arnold ML, Ballerini ES, Brothers AN. Hybrid fitness, adaptation and evolutionary diversification: lessons learned from Louisiana irises. Heredity. 2012;108:159–66.
Baiz MD, Tucker PK, Cortés-Ortiz L. Multiple forms of selection shape reproductive isolation in a primate hybrid zone. Mol Ecol. 2019;28:1056–69.
Smith KL, Hale JM, Kearney MR, Austin JJ, Melville J. Molecular patterns of introgression in a classic hybrid zone between the Australian tree frogs, Litoria ewingii and L. paraewingi: evidence of a tension zone. Mol Ecol. 2013;22:1869–83.
Hülber K, Sonnleitner M, Suda J, Krejčíková J, Schönswetter P, Schneeweiss GM, Winkler M. Ecological differentiation, lack of hybrids involving diploids, and asymmetric gene flow between polyploids in narrow contact zones of Senecio carniolicus (syn. Jacobaea carniolica, Asteraceae). Ecol Evol. 2015;5:1224–34.
Testo WL, Watkins JE, Jr., Barrington DS. Dynamics of asymmetrical hybridization in North American wood ferns: reconciling patterns of inheritance with gametophyte reproductive biology. New Phytol. 2015;206:785–95.
Rieseberg LH, Carney SE. Plant hybridization. New Phytol. 1998;140:599–624.
Abbott RJ. Plant speciation across environmental gradients and the occurrence and nature of hybrid zones. J Syst Evol. 2017;55:238–58.
Arntzen JW, Vries WD, Canestrelli D, Martínez-Solano I. Hybrid zone formation and contrasting outcomes of secondary contact over transects in common toads. Mol Ecol. 2017;26:5663–75.
Pickup M, Brandvain Y, Fraïsse C, Yakimowski S, Barton NH, Dixit T, Lexer C, Cereghetti E, Field DL. Mating system variation in hybrid zones: facilitation, barriers and asymmetries to gene flow. New Phytol. 2019;224:1035–47.
Anderson E. Hybridization of the habitat. Evolution. 1948;2:1–9.
Hall RJ. Hybridization helps colonizers become conquerors. Proc Natl Acad Sci U S A. 2016;113:9963–4.
Barrera-Guzmán AO, Aleixo A, Shawkey MD, Weir JT. Hybrid speciation leads to novel male secondary sexual ornamentation of an Amazonian bird. Proc Natl Acad Sci U S A. 2018;115:E218-25.
Leroy T, Louvet JM, Lalanne C, Le Provost G, Labadie K, Aury JM, Delzon S, Plomion C, Kremer A. Adaptive introgression as a driver of local adaptation to climate in European white oaks. New Phytol. 2020;226:1171–82.
Liu YF, Li DW, Zhang Q, Song C, Zhong CH, Zhang XD, Wang Y, Yao XH, Wang ZP, Zeng SH et al. Rapid radiations of both kiwifruit hybrid lineages and their parents shed light on a two-layer mode of species diversification. New Phytol. 2017;215:877–90.
Chamberlain D, Hyam R, Argent G, Fairweather G, Walter KS. The genus Rhododendron: its classification and synonymy.Place:Royal Botanic Garden Edinburgh;1996.
Milne RI. Phylogeny and biogeography of Rhododendron subsection Pontica, a group with a tertiary relict distribution. Mol Phylogen Evol. 2004;33:389–401.
Soza VL, Lindsley D, Waalkes A, Ramage E, Patwardhan RP, Burton JN, Adey A, Kumar A, Qiu RL, Shendure J et al. The Rhododendron genome and chromosomal organization provide insight into shared whole-genome duplications across the heath family (Ericaceae). Genome Biol Evol. 2019;11:3353–71.
Yan LJ, Liu J, Möller M, Zhang L, Zhang XM, Li DZ, Gao LM. DNA barcoding of Rhododendron (Ericaceae), the largest Chinese plant genus in biodiversity hotspots of the Himalaya-Hengduan Mountains. Mol Ecol Resour. 2015;15:932–44.
Zou JY, Luo YH, Burgess KS, Tan SL, Zheng W, Fu CN, Xu K, Gao LM. Joint effect of phylogenetic relatedness and trait selection on the elevational distribution of Rhododendron species. J Syst Evol. 2020. https://doi.org/10.1111/jse.12690.
Shrestha N, Su XY, Xu XT, Wang ZH. The drivers of high Rhododendron diversity in south-west China: Does seasonality matter? J Biogeogr. 2018;45:438–47.
Handel-Mazzetti H. Die Abnahme eines Teiles verpflichtet zur Abnanme des ganzen Werkes. Symbolae Sinicae. 1936;7:775.
Zha HG, Milne RI, Sun H. Morphological and molecular evidence of natural hybridization between two distantly related Rhododendron species from the Sino-Himalaya. Bot J Linn Soc. 2008;156:119–29.
Ma YP, Zhang CQ, Zhang JL, Yang JB. Natural hybridization between Rhododendron delavayi and R. cyanocarpum (Ericaceae), from morphological, molecular and reproductive evidence. J Integr Plant Biol. 2010;52:844-51.
Zhuang P. Progress on the fertility of Rhododendron. Biodiversity Science. 2019;27(3):327–38.
Yan LJ, Gao LM, Li DZ. Molecular evidence for natural hybridization between Rhododendron spiciferum and R. spinuliferum (Ericaceae). J Syst Evol. 2013;51:426–34.
Richards TJ, Walter GM, McGuigan K, Ortiz-Barrientos D. Divergent natural selection drives the evolution of reproductive isolation in an Australian wildflower. Evolution. 2016;70:1993–2003.
Seehausen O. Conditions when hybridization might predispose populations for adaptive radiation. J Evol Biol. 2013;26:279–81.
Milne RI, Davies C, Prickett R, Inns LH, Chamberlain DF. Phylogeny of Rhododendron subgenus Hymenanthes based on chloroplast DNA markers: between-lineage hybridisation during adaptive radiation? Plant Syst Evol. 2010;285:233–44.
Fu CN, Mo ZQ, Yang JB, Cai J, Ye LJ, Zou JY, Qin HT, Zheng W, Hollingsworth PM, Li DZ, Gao LM. Testing genome skimming for species discrimination in the large and taxonomically difficult genus Rhododendron. Mol Ecol Resour. 2021;00:1–11.
Ding WN, Ree RH, Spicer RA, Xing YW. Ancient orogenic and monsoon-driven assembly of the world’s richest temperate alpine flora. Science. 2020;369,578–581.
Zha HG, Milne RI, Sun H. Asymmetric hybridization in Rhododendron agastum: a hybrid taxon comprising mainly F1s in Yunnan, China. Ann Bot. 2010;105:89–100.
Marczewski T, Chamberlain DF, Milne RI. Hybridization in closely related Rhododendron species: half of all species-differentiating markers experience serious transmission ratio distortion. Ecol Evol. 2015;5:3003–22.
Ma YP, Xie WJ, Sun WB, Marczewski T. Strong reproductive isolation despite occasional hybridization between a widely distributed and a narrow endemic Rhododendron species. Sci Rep. 2016;6:19146.
Zhang JL, Ma YP, Wu ZK, Dong K, Zheng SL, Wang YY. Natural hybridization and introgression among sympatrically distributed Rhododendron species in Guizhou, China. Biochem Syst Ecol. 2017;70:268–73.
Zhang XM, Qin HT, Xie WJ, Ma YP, Sun WB. Comparative population genetic analyses suggest hybrid origin of Rhododendron pubicostatum, an endangered plant species with extremely small populations endemic to Yunnan, China. Plant Diversity. 2020;42:312–8.
Yan LJ, Burgess KS, Zheng W, Tao ZB, Li DZ, Gao LM. Incomplete reproductive isolation between Rhododendron taxa enables hybrid formation and persistence. J Integr Plant Biol. 2019;61:433–48.
Karrenberg S, Liu XD, Hallander E, Favre A, Herforth-Rahmé J, Widmer A. Ecological divergence plays an important role in strong but complex reproductive isolation in campions (Silene). Evolution. 2019;73:245–61.
Chapman MA, Hiscock SJ, Filatov DA. Genomic divergence during speciation driven by adaptation to altitude. Mol Biol Evol. 2013;30:2553–67.
Scopece G, Palma-Silva C, Cafasso D, Lexer C, Cozzolino S. Phenotypic expression of floral traits in hybrid zones provides insights into their genetic architecture. New Phytol. 2020;227:967–75.
Fulton M, Hodges SA. Floral isolation between Aquilegia formosa and Aquilegia pubescens. P Roy Soc B-Biol Sci. 1999;266:2247–52.
Giudicelli GC, Turchetto C, Teixeira MC, Freitas LB. Morphological and genetic characterization in putative hybrid zones of Petunia axillaris subsp. axillaris and subsp. parodii (Solanaceae). Bot J Linn Soc. 2019;191:353–64.
Barton NH. The role of hybridization in evolution. Mol Ecol. 2001;10:551–68.
Grant PR, Grant BR. Hybridization increases population variation during adaptive radiation. Proc Natl Acad Sci U S A. 2019;116:23216–24.
Chhatre VE, Evans LM, DiFazio SP, Keller SR. Adaptive introgression and maintenance of a trispecies hybrid complex in range-edge populations of Populus. Mol Ecol. 2018;27:4820–38.
Dudek K, Gaczorek TS, Zieliński P, Babik W. Massive introgression of major histocompatibility complex (MHC) genes in newt hybrid zones. Mol Ecol. 2019;28:4798–810.
Grant PR, Grant BR. Triad hybridization via a conduit species. Proc Natl Acad Sci U S A. 2020;117:7888–96.
Doyle JJ, Doyle JL. A rapid DNA isolation procedure for small quantities of fresh leaf tissue. Phytochemical Bulletin. 1987;19(1):11–5.
Yang GQ, Chen YM, Wang JP, Guo C, Zhao L, Wang XY, Guo Y, Li L, Li DZ, Guo ZH. Development of a universal and simplified ddRAD library preparation approach for SNP discovery and genotyping in angiosperm plants. Plant Methods. 2016;12:39.
Tang Y, Horikoshi M, Li WX. ggfortify: Unified interface to visualize statistical results of popular R packages. R Journal. 2016;8:474-85.
Team RC. Changes in R from version 3.4.2 to version 3.4.3. R Journal. 2017;9:568–70.
IBM Corp. Released 2010. IBM SPSS Statistics for Windows, Version 19.0. Armonk: IBM Corp.
Rochette NC, Rivera-Colón AG, Catchen JM. Stacks 2: Analytical methods for paired-end sequencing improve RADseq-based population genomics. Mol Ecol. 2019;28:4737–54.
Andrews S. FastQC a quality control tool for high throughput sequence data. http://www.bio.informatics.babraham.ac.uk/projects/fastqc/. 2014.
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MA, Bender D, Maller J, Sklar P, de Bakker PI, Daly MJ et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.
Jeffery NW, DiBacco C, Wringe BF, Stanley RRE, Hamilton LC, Ravindran PN, Bradbury IR. Genomic evidence of hybridization between two independent invasions of European green crab (Carcinus maenas) in the Northwest Atlantic. Heredity. 2017;119:154–65.
Wringe BF, Stanley RRE, Jeffery NW, Anderson EC, Bradbury IR. HYBRIDDETECTIVE: A workflow and package to facilitate the detection of hybridization using genomic data in R. Mol Ecol Resour. 2017;17:e275-84.
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19:1655–64.
Paradis E, Claude J, Strimmer K. APE: analyses of phylogenetics and evolution in R language. Bioinformatics. 2004;20:289–90.
Kamvar ZN, Tabima JF, Grünwald NJ. Poppr: an R package for genetic analysis of populations with clonal, partially clonal, and/or sexual reproduction. PeerJ. 2014;2:e281.
Ginestet C. ggplot2: Elegant graphics for data analysis. J R Stat Soc Ser A-Stat Soc. 2011;174:245–6.
Anderson EC, Thompson EA. A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002;160:1217–29.
Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, Buxton S, Cooper A, Markowitz S, Duran C et al. Geneious Basic: an integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28:1647–9.
We thank Drs. Jie Liu, Guoqian Yang, Miss. Zhiqiong Mo, Mr. Linjiang Ye for their assistance in the field and lab work and their constructive suggestions on an early version of the manuscript. Laboratory work was carried out at the Laboratory of Molecular Biology in the Germplasm Bank of Wild Species in Southwest China, Kunming Institute of Botany, Chinese Academy of Sciences.
This work was supported by the National Natural Science Foundation of China (31670213, 31700179), the Strategic Priority Research Program of Chinese Academy of Sciences (XDB31000000), the Large-scale Scientific Facilities of the Chinese Academy of Sciences (2017-LSFGBOWS-02), the Program of Science and Technology Talents Training of Yunnan Province, China (2016HA005, 2017HA014), Yunnan Fundamental Research Project (202001AT070062), Construction of International Flower Technology Innovation Center and Industrialization of achievements (2019ZG006) and CAS President’s International Fellowship Initiative (2017VBB0008). The funding body played no role in study design, data analysis, and manuscript preparation.
Ethics approval and consent to participate
All plant materials used in this study were obtained from the hybrid zone and surrounding areas, which do not belong to nature reserves. In addition, all research materials did not involve any endangered or protected species.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
About this article
Cite this article
Zheng, W., Yan, LJ., Burgess, K.S. et al. Natural hybridization among three Rhododendron species (Ericaceae) revealed by morphological and genomic evidence. BMC Plant Biol 21, 529 (2021). https://doi.org/10.1186/s12870-021-03312-y