1. Introduction

Natural hybridization plays an important role in speciation, genetic exchange, and adaptive evolution (Rieseberg and Willis 2007; Abbott et al., 2013; Ciezarek et al., 2024). The formation of a hybrid zone typically begins with two parental species that are isolated by spatial or ecological niche differentiation (Barton and Hewitt 1985; Harrison and Larson 2016). Upon contact, the pre- and post-reproductive isolating mechanisms between these species may break down, facilitating the formation of hybrids (Abbott et al., 2013; Abbott 2017; Nachychko and Sosnovsky, 2021; Mitchell et al., 2022). In some cases, hybrid zones are established (Arnold et al., 2012; Flowers et al., 2019). If the F₁ hybrids exhibit even slight fertility (Arnold et al., 1999), complex population genetic structure may arise within these hybrid zones, characterized by multiple hybrid generations that possess relatively high fitness and exhibit weak reproductive isolation from the parental species (Akopyan et al., 2020; Chiou et al., 2021).

The evolutionary consequences of hybridization can yield several outcomes, including the reinforcement of species boundaries (Abbott et al., 2013), hybrid speciation (Rieseberg, 2007; Todesco et al., 2016), introgression (DeVos et al., 2023; Groha and Coopa, 2024), and the possibility of extinction of rare parental species (Todesco et al., 2016). These outcomes are largely influenced by the population genetic structure of hybrid zones (Harrison and Larson 2016; Zieliński et al., 2019; Westram et al., 2021), which can manifest in various forms, such as tension zones, bounded hybrid superiority zones, mosaic hybrid zones, and F₁-dominated hybrid zones (Arnold et al., 1999; Abbott and Brennan, 2014). The population genetic structure of hybrid zones, in turn, is shaped by several factors, including the strength of the reproductive isolation barriers between the hybridizing taxa (Abbott et al., 2013), relative hybrid fitness (DeVos et al., 2023), and the proportions and distribution of the hybridizing taxa involved (Abbott et al., 2013; Mandeville et al., 2019; DeVos et al., 2023). Additionally, hybrid zones may occupy habitats that are distinct from those of one or both parental taxa, a phenomenon facilitated by niche partitioning and local adaptation resulting from novel combinations of genes and traits (Rieseberg et al., 2003; Arnold and Martin 2010; Arnold et al., 2012).

The population genetic structure of a hybrid zone can vary substantially even among hybrid zones formed by the same taxa. This variation is often shaped by a combination of environmental and genetic factors, including spatial distribution, abundance, habitat characteristics, and gene flow (Milne et al., 2003; Burgess et al., 2005; Abbott and Brennan 2014). Such factors have been investigated in genera such as Populus, Ophrys, and Lissotriton (Cotrim et al., 2016; Zieliński et al., 2019; Bolte et al., 2024). To gain a comprehensive understanding of how these factors influence the fate of the hybrid zones, it is essential to compare multiple hybrid zones formed by the same parental species. However, such comparative studies of this nature remain relatively limited (Mandeville et al., 2019; Zieliński et al., 2019).

Rhododendron L. (Ericaceae) is one of the largest and most taxonomically complex genera in the northern hemisphere (Chamberlain et al., 1996; Mo et al., 2022), with one of its centres of diversity located in the Himalaya-Hengduan Mountains (Shrestha et al., 2018; Fu et al., 2022). Natural hybridization is frequent within Rhododendron (Milne 2004; Zhang et al., 2007, 2020; Ma et al., 2022) and has played a crucial role in its extensive adaptive radiation and reticulate evolution (Soza et al., 2019; Ma et al., 2022; Mo et al., 2022; Xia et al., 2022). For instance, asymmetric and bidirectional hybridization between Rhododendron spiciferum Franch. and R. spinuliferum Franch. has been found to produce the natural hybrid R. × duclouxii H. Lévl. (Yan et al., 2017). The hybrid (R. × duclouxii) and its two parental species are diploid and long-lived shrubs. R. spinuliferum features red, tubular flowers conforming to a bird-pollination syndrome, while R. spiciferum has pink, open corollas indicative of bee pollination, with the hybrid flowers displaying varying traits between the two parents (Yan et al., 2013, 2019). This hybridization event has led to the establishment of at least 15 natural hybrid zones in Yunnan Province, which typically encompass multiple classes of offspring (e.g., F₁, F₂, backcrosses), although some sites may contain few F₁ individuals (Yan et al., 2017). Currently, the processes underlying the formation, variation in genetic structure, and potential evolutionary outcomes of these hybrid zones remain poorly understood. However, the observable structural variation among R. × duclouxii hybrid zones makes it a valuable system for investigating the factors that determine the trajectory of hybridization events.

In this study, we integrated genomic ddRAD-seq data, along with morphological and ecological analyses, to investigate the formation, population genetic structure, direction of hybridization, and potential evolutionary trajectory of hybrid zones involving Rhododendron × duclouxii. Specifically, we aimed to address the following questions: 1) Are the multiple hybrid zones independent? 2) What factors influence the genetic structure of each hybrid zone? and 3) What are the possible evolutionary outcomes of hybridization between R. spiciferum and R. spinuliferum? This study will deepen our understanding of the R. × duclouxii hybridization event and provides insights into the evolutionary outcomes of hybridization and the diversification of the genus Rhododendron.

2. Materials and methods 2.1. Study site and sampling

Seven hybrid zones involving Rhododendron × duclouxii (that included hybrids and both parental species, R. spinuliferum and R. spiciferum, in sympatry) and two allopatric populations (i.e., populations not in proximity to each other and any hybrid zone) each for R. spinuliferum and R. spiciferum were sampled from their core distribution ranges in central Yunnan Province, SW China (Fig. 1 and Table S1). A total of 258 individuals were sampled, comprising 71 R. spiciferum (SC), 88 R. spinuliferum (SN), 98 R. × duclouxii (HY), and one individual of R. racemosum Franch. (RA), which served as an outgroup for molecular analyses (Table S2). Among the samples, 18 individuals each of R. spiciferum and R. spinuliferum were collected from the four allopatric populations, while the remainder were collected from seven hybridizing populations (hybrid zones).

For each individual, fresh healthy leaves were collected and immediately dried with silica gel before DNA extraction. Voucher specimens were also collected from each of the 258 sampled individuals for morphological trait measurement. All specimens were deposited in the herbarium of Kunming Institute of Botany, Chinese Academy of Sciences (KUN).

2.2. Data collection

To assess morphological variations, 23 morphological traits were measured from voucher specimens of all sampled 257 individuals for R. spiciferum, R. spinuliferum and R. × duclouxii, with the method of traits data collection following Zheng et al. (2021). These traits included fourteen quantitative traits (leaf length, width, thickness, and area; petiole length; pedicel length; corolla lobe length; corolla tube width, and length; flower width; style length; filaments length; stigma width; ovary width) and nine qualitative traits (calyx lobe conspicuous/inconspicuous; abaxial and adaxial leaf hairs density; corolla scales density; presence/absence of filaments and style hair; flower color; stigma color; and anther color).

Additionally, nine environmental (soil physicochemical properties) and thirteen functional traits were collected from living plants in four hybrid zones (AN, AZY, SHB, WD). Five individuals were selected from 20 × 20 m plots separately dominated by R. spiciferum, R. spinuliferum, and R. × duclouxii in each hybrid zone. For each of these individuals, soil samples were collected from the area near their roots, and thirteen functional traits (leaf fresh weight, leaf dry weight, leaf dry matter content (LDMC), leaf area, style length, pistil length, stamen length, corolla length, specific leaf area (SLA), petal length, flower width, chlorophyll content, and leaf thickness) related to energy capture and pollination were measured from five flowers and five leaves per plant (Table S3). Length measurements were taken using a vernier calliper, and weight was measured using an analytical balance. Chlorophyll content was measured with a SPAD 502 Plus chlorophyll meter. Soil samples were sent to Yunnan Tongchuan Agricultural Analytical and Testing Technology Co., Ltd. (Kunming, Yunnan, China) to analyse nine physicochemical properties (Table S3).

2.3. ddRAD-seq and cpDNA trnL-F sequencing

Total genomic DNA was extracted using a modified CTAB method (Doyle and Doyle, 1987). Library construction for double digest restriction-site-associated DNA sequencing (ddRAD-seq) was performed following the protocol of Yang et al. (2016). Library preparation and sequencing were conducted following Zheng et al. (2021). To determine chloroplast DNA (cpDNA) haplotypes, the trnL-F cpDNA region was sequenced for all 221 individuals sampled from the hybrid zones (i.e., excluding individuals from allopatric populations and the outgroup Rhododendron racemosum). PCR amplification and sequencing were performed following Yan et al. (2013).

2.4. Read mapping and SNP calling

The ddRAD-seq data were filtered following Zheng et al. (2021). The genome of Rhododendron simsii Planch. (Yang et al., 2020) was used as the reference in the reference-based analysis pipeline. This reference genome was 529 Mb in size with 13 chromosomes. The N50 of contigs and scaffolds was 2.2 Mb and 36 Mb, respectively. Loci filtering was conducted in two phases. In the first phase, the number of single nucleotide polymorphisms (SNPs) was reduced based on the following parameters: write_single_snp, r = 0.7, p = 22, and min-maf = 0.05 (Zhang et al., 2018). In the second phase, SNPs were further filtered by Plink v.1.90 with parameters geno = 0.2 and hwe = 0.01 (Purcell et al., 2007). The SNPs retained after these two filtering steps were compiled into Dataset_A. Subsequently, BayeScan (Foll and Gaggiotti, 2008) was employed to identify and remove SNPs under selection, resulting in Dataset_B for population genetic analysis. To detect hybrids in the NewHybrids analysis, a third dataset (Dataset_C), comprising 400 SNPs, selected based on high pairwise F_ST values and low linkage disequilibrium (LD), was generated using the getTopLoc function in the hybriddetective R package (Jeffery et al., 2017; Wringe et al., 2017).

2.5. Data analysis 2.5.1. Morphological trait differentiation

Principal Coordinates Analysis (PCoA) based on the 23 morphological traits (i.e., Gower distance) was conducted using the R package ggfortify to assess trait dispersion in a multivariate trait space (Tang et al., 2016). This analysis included 254 individuals, composing 88 R. spinuliferum, 98 R. × duclouxii, and 68 R. spiciferum individuals; the three remaining individuals were excluded due to insufficient morphological data. Box plots were generated to illustrate the variation within the 14 quantitative traits using R, and analysis of variance (ANOVA) for each quantitative trait was also performed in R.

2.5.2. Climate data and range modelling

To investigate the niche differentiation among the three taxa, distribution data were gathered from databases including the Chinese Virtual Herbarium (CVH: http://www.cvh.ac.cn/), the Global Biodiversity Information Facility (GBIF: https://www.gbif.org/), and additional sources, along with our field investigations. Duplicate records and occurrences within a 1 km range of one another were removed. Climate data used for model simulation were obtained from WorldClim (https://worldclim.org/), encompassing 19 climate variables derived in the ACCESS1.0 model at a 30 arc-second spatial resolution (Table S4). Climate variables that were highly correlated with each other (variance inflation factor (VIF) > 10) were excluded using the R package usdm. Finally, six climate variables (bio2, bio7, bio8, bio13, bio15, bio19) remained for data analysis.

We simulated the currently suitable ecological niche area for the three taxa ten times using MaxEnt v.3.4.1 with the default parameters, and the average results were used to map the existing suitable areas for R. spinuliferum, R. spiciferum and R. × duclouxii respectively in ArcGIS v.10.2. Additionally, the six climate variables for each distribution point were extracted in ArcGIS and then used for analyses of PCA (Principal Component Analysis) and DFA (discriminant function analysis) (Li et al., 2019).

2.5.3. Functional traits and soil physicochemical differentiation

PCA was conducted in R (using the ggfority package) with the 13 leaf and flower functional traits. ANOVA was used to evaluate the variation of each functional trait in R. The soil physicochemical properties were also analyzed via PCA and ANOVA.

2.5.4. Population genetic variation and structure

Population genetic differentiation (F_ST) was calculated using Arlequin v.3.5 (Excoffier and Lischer 2010) with Dataset_B. The R package gplots was used to plot the F_ST data with the heatmap.2 function. The genetic structure of each hybrid zone and allopatric parental populations was inferred using Admixture v.1.3.0 (Alexander et al., 2009) on Dataset_B, with K values ranging from 2 to 7 to investigate the genetic clustering of hybrid zones. The best K value was determined by cross-validation (CV). To visualize the genetic similarities among the three taxa, a PCoA was performed using the ape package in R (Paradis et al., 2004). A Minimum Spanning Network (MSN) analysis for all R. × duclouxii individuals was conducted using the R package poppr (Kamvar et al., 2014).

NewHybrids v.1.1 (Anderson and Thompson 2002) was employed to estimate the likelihood of individuals belonging to specific genotypic classes based on Dataset_C. This analysis recognized six classes and individuals were assigned to a class if their probability of belonging to a class was ≥0.9. The analysis was run for 100, 000 iterations after a burn-in of 100, 000. To explore genetic and morphological differences between classes, SNPs from Dataset_B and morphological trait data were integrated with NewHybrids results to perform PCoA in R.

To detect the direction of hybridization between parental species, the chloroplast trnL-F regions from hybrid zone individuals were aligned in Geneious v.8.1 (Kearse et al., 2012). Informative sites were summarized, and the haplotypes of R. spiciferum and R. spinuliferum were assigned according to Yan et al. (2013).

2.5.5. Gene flow and introgression

Gene flow between the parental species was detected by TreeMix with Dataset_B (Pickrell and Pritchard 2012), using one Rhododendron racemosum sample as the outgroup. F-branch analysis with Dsuite was also used to test the gene flow between populations (Malinsky et al., 2021). Barrier v.2.2 (Manni et al., 2004) was used to detect the geographic barriers between the populations of the parental species. Isolation by distance (IBD) analysis was performed using the R package ade4, with geographical distance and F_ST as the covariates.

The program BGC was used to perform introgression analysis with Dataset_A. This analysis was conducted for all individuals, and separately for each hybrid zone (Gompert et al., 2012). Results were plotted using the ClineHelpR package in R (Martin et al., 2021), and a Venn diagram was further generated using the OmicStudio tools (https://www.omicstudio.cn/tool). The SNPs associated with selection and introgression were annotated by Annovar (Wang et al., 2010).

3. Results 3.1. Morphological trait differentiation

In the PCoA analysis of the 23 morphological traits across 254 individuals, the first and second components accounted for 75.29% and 8.03% of the total variance, respectively (Fig. 2a). Individuals from the two parental species (Rhododendron spiciferum and R. spinuliferum) formed two distinct clusters, while R. × duclouxii individuals clustered between them, with separation primarily along the first component. Within these clusters, R. × duclouxii individuals exhibited greater morphological variability than either parental species, showing substantial overlap with R. spiciferum but more limited overlap with R. spinuliferum (Fig. 2a). Notably, R. × duclouxii individuals did not display significant differences among hybrid zones (Fig. S1). When considering only the first component, a comparison of individuals inside and outside the hybrid zones revealed that R. spinuliferum and R. spiciferum individuals outside the hybrid zones were relatively more different from R. × duclouxii than those within the hybrid zones in morphology (Fig. S2). Based on the offspring classification from the NewHybrids analysis (Table S2), the floral traits of F₁ (first generation) hybrids were more similar to R. spiciferum than to R. spinuliferum (Fig. 2b). Both F₁ and F₂ (second generation) hybrids clustered between the two parental species, but were more similar to R. spiciferum than to R. spinuliferum, with F₂ hybrids exhibiting greater variation than those of F₁ hybrids (Fig. 2c). The degree of morphological overlap between the hybrids and parental species varied across the seven hybrid zones (Fig. S3). ANOVA results indicated that 13 of the 14 quantitative morphological traits effectively distinguished the parental species, with R. × duclouxii generally displaying intermediate values. However, for flower width and corolla length, the R. × duclouxii plants showed higher mean values than either parental species (Fig. S4). Seven morphological traits showed significant differences across the hybrid zones for R. × duclouxii (Fig. 2d). Notably, in two hybrid zones (DZ and XD), the mean values for corolla tube width, corolla tube length and filaments length were lower than those observed in the four hybrid zones (HXL, WD, SHB and AZY).

3.2. Ecological niche differentiation

The results of the MaxEnt analysis revealed considerable overlap in distribution and ecological niche preferences between R. spinuliferum and R. spiciferum, particularly in central Yunnan Province, indicating that they occupy similar habitat types (Fig. 3a and Fig. S5). The overlap area between the parental species was almost entirely shared with R. × duclouxii, suggesting that hybrids do not occupy a distinct ecological niche (Fig. 3a). Neither PCA nor DFA analyses could differentiate the climate preferences of the three taxa (Fig. 3b and Fig. S6).

Of the 13 leaf and floral functional traits examined, 9 exhibited significant differences among the three taxa. Additionally, three traits (SLA, petal length, and chlorophyll content) distinguished R. spinuliferum from R. spiciferum, whereas R. × duclouxii did not show significant differentiation from its parents (P < 0.001). Leaf thickness was the only trait that did not differentiate any of the three taxa (Fig. S7). In the PCA of these traits, the first two axes accounted for 80.16% and 8.15% of the variance respectively, with R. × duclouxii situated between the two parental species along the first axis (Fig. 3c).

Significant differences in soil physicochemical properties were observed among the three taxa for four out of the nine variables analyzed (Table S5). Specifically, soil from the root zone of R. spinuliferum exhibited higher phosphorus (P) content compared to the other taxa. Conversely, soil associated with R. spiciferum showed lower organic carbon (OC) and nitrogen (N) contents. Additionally, soil from the root zone of R. spiciferum differed in pH from that of R. × duclouxii, while the pH of R. spinuliferum soil was intermediate and did not differ significantly from either taxon (Table S5). Importantly, the range of variation in most soil variables was much greater among sites than among the taxa (Table S5). The PCA of all soil physicochemical variables revealed a considerable overlap among the three taxa (Fig. 3d).

3.3. Population genetic variation

After mapping the reads to the reference genome of Rhododendron simsii (Yang et al., 2020), an average of 2, 888, 228 loci (RAD tags) were assembled for each of the 258 individuals. Following filtering criteria using Plink, a total of 5701 SNPs were identified in Dataset_A. Out of these, 388 SNPs were determined to be under selection according to BayeScan and were subsequently excluded from the population genetic analyses (Fig. S8), some of these SNPs were associated with functions related to plant resistance and growth (Table S6). The remaining 5313 SNPs were used to construct Dataset_B.

Population genetics analysis of Dataset_B revealed significant genetic differentiation between R. spinuliferum and R. spiciferum (F_ST = 0.15; Table S7), while no significant differentiation was detected between R. × duclouxii and either R. spinuliferum (F_ST = 0.05) or R. spiciferum (F_ST = 0.03). Among populations, the highest levels of genetic differentiation were observed between the allopatric populations of the two parental species (F_ST > 0.27), followed by the comparisons involving one allopatric parental population and the other parental species within hybrid zones (F_ST = 0.20–0.27), and subsequently between the two parental species within hybrid zones (F_ST = 0.02–0.13) (Fig. S9 and Table S8). Genetic differentiation among different populations of the same species was relatively low (F_ST = 0.02–0.16 for R. spinuliferum and F_ST = 0.03–0.10 for R. spiciferum), and that among populations of R. × duclouxii much lower (F_ST = 0.01–0.05). Notably, the greatest genetic differentiation for R. × duclouxii was observed between the DZ hybrid zone and all other hybrid zones.

3.4. Population genetic structure

The Admixture analysis showed that the best K value was 3 (Fig. S10). At K = 2, the two parental species (R. spiciferum and R. spinuliferum) formed distinct genetic groups, while all R. × duclouxii individuals exhibited genetic admixture between them (Fig. 4a). When K = 3, the R. spinuliferum populations were further divided into western and eastern clusters. Most R. × duclouxii individuals were found to be admixed between R. spiciferum and the eastern R. spinuliferum cluster, except for population DZ, which showed admixture between R. spiciferum and the western R. spinuliferum cluster. K = 4 did not identify R. × duclouxii as an independent group; instead, it consistently displayed a high degree of admixture with significant variation within and among populations. The degree of hybridization and the composition of hybrids varied across hybrid zones, with the western hybrid population DZ consistently exhibiting distinct characteristics compared to the other zones (Fig. 4a and Fig. S11).

The PCoA analysis corroborated the Admixture results (Fig. 4b). The first principal axis accounted for 8.48% of the variance, effectively separating Rhododendron spiciferum and R. spinuliferum into two well-differentiated clusters. R. × duclouxii individuals were positioned between R. spinuliferum and R. spiciferum, showing some overlap with R. spiciferum. The second axis explained 1.98% of the variance, further subdividing R. spinuliferum into eastern and western clusters, and R. spiciferum into a northern cluster that occurred in the hybrid zones and a southern allopatric cluster (Fig. 4b). Correspondingly, R. × duclouxii individuals also exhibited a distinct geographical pattern akin to that of R. spinuliferum.

When only Rhododendron × duclouxii individuals were examined across all hybrid zones, the first two principal axes explained 3.25% and 2.69% of the variance respectively (Fig. 4c). R. × duclouxii individuals from five hybrid zones (AN, HXL, SHB, WD and AZY) formed a compact cluster, whereas those from hybrid zone DZ formed a separate, distinct cluster. R. × duclouxii individuals from hybrid zone XD formed a loose cluster, with some individuals overlapping with the five-zone cluster, and others highly distinct from it. The MSN analysis showed a similar pattern among R. × duclouxii populations: the DZ zone formed a distinct cluster, while XD exhibited a mix of distinct and overlapping distribution (Fig. S12).

3.5. Hybrid identification

The NewHybrids analysis classified 32% of R. × duclouxii individuals as F₁ (first generation), 27% as F₂ (second generation), 18.4% as backcross to R. spiciferum, and 13.3% as backcross to R. spinuliferum (Fig. 4d). Eleven R. × duclouxii individuals could not be assigned to a specific generation, possibly representing F3 (third generation) or later-generation hybrids. Class composition varied across hybrid zones, with hybrid zones AN and HXLdominated by F₁s, whereas hybrid zone DZ comprised mostly F₁s and backcrosses to R. spiciferum and R. spinuliferum. In contrast, no F₁ individuals were detected in the hybrid zone XD. In the hybrid zones AZY and WD, F₂s were the most frequent class, while the hybrid zone SHB had a higher proportion of backcross to R. spiciferum than the other hybrid zones.

When integrating the NewHybrids results with SNPs data, the first principal component of the PCoA effectively distinguished the parental species and revealed larger genetic variation among F₂s compared to F₁s (excluding individuals from hybrid zone DZ) (Fig. S13). As expected, backcross individuals were positioned between F₁s and their respective parental species (Fig. S13). The combination of the NewHybrids and F_ST results indicated minimal genetic differentiation between F₁ and F₂ hybrids (Fig. S14).

3.6. Haplotype variation

The majority of R. × duclouxii individuals possessed the R. spiciferum chloroplast trnL-F haplotype (78%), whereas only 22% of individuals displayed the R. spinuliferum haplotype (Fig. 4d). The proportion of these haplotypes varied across hybrid zones, with hybrid zone XD exclusively containing the R. spiciferum haplotype, while hybrid zones AN and WD contained a higher proportion of the R. spinuliferum haplotype (Fig. 4d).

3.7. Gene flow

No evidence of gene flow was detected between different populations of R. spinuliferum, as the optimal migration parameter (m) derived from the Treemix analysis was 0. However, the optimal m value for the R. spiciferum populations was 2, indicating gene flow from the outgroup R. racemosum into populations AZY and DZ. No gene flow was detected between the R. spiciferum populations themselves (Fig. 5a), suggesting that the levels of gene flow may be too low to detect. When R. spinuliferum and R. spiciferum populations were analyzed together with Treemix, the optimal m value remained 2, indicating gene flow from R. spinuliferum to R. spiciferum (specifically in population DZ) and from the R. racemosum to population AZY of R. spiciferum, albeit with a weak signal for the latter.

The Treemix analysis was conducted for each hybrid zone to assess the direction and intensity of gene flow between parental species. Gene flow was detected in five hybrid zones, with four exhibiting weak gene flow from R. spinuliferum to R. spiciferum (Fig. S15), while the hybrid zone AZY showed gene flow in the opposite direction. No gene flow was detected within the hybrid zones AN and HXL. The results of the F-branch analysis were consistent with those obtained from Treemix, revealing multiple ancient gene flow events. We also identified some weak gene flow in R. spiciferum between population AN and other populations, as well as in R. spinuliferum between populations XD and DZ, DZ and SHB, and XD and WD (Fig. S16).

The Barrier analysis indicated that the most significant barrier to gene flow for R. spinuliferum occurred between population DZ and AN. Similarly, for R. spiciferum, the strongest barrier to gene flow was located between population DZ and AN, with the second strongest barrier situated between the southern allopatric populations (BFZ and JCZ) and the northern hybrid zone populations (Fig. 5b). The IBD analysis revealed a significant correlation between geographical distance and F_ST for both R. spinuliferum and R. spiciferum populations, indicating that genetic differentiation increased with increasing geographical distance (Fig. 5c).

3.8. Introgression

The introgression analysis showed distinct patterns of introgression across the hybrid zones examined. In the hybrid zones AN and HXL, the hybrid indices exhibited no significant excess of SNPs or outliers when compared to genome-wide expectations (Table 1). However, in the remaining five hybrid zones, only SNPs classified as α estimate-type were found in excess, whereas β estimate-type SNPs were not (Table 1). The five hybrid zones displayed varying numbers of introgressed SNPs. The hybrid zones SHB and WD had relatively few non-neutral introgressed SNPs, with a roughly equal distribution of α > 0 (introgression from R. spiciferum to R. spinuliferum) and α < 0 (introgression from R. spinuliferum to R. spiciferum) SNPs. In contrast, hybrid zone AZY exhibited a higher number of α > 0 SNPs, while the hybrid zones XD and DZ displayed a predominance of α < 0 SNPs. Notably, the loci associated with introgression were largely specific to each hybrid zone, with minimal overlap observed among populations (Fig. S17).

The annotation of the identified SNPs indicated that the introgressed SNPs were associated with 111 genes, several of which are linked to critical physiological functions. For instance, transcription factor MYB and genes such as AUX1 (AUXIN1) and CDH1 (CADHERIN1) were identified, suggesting potential adaptive introgression (see Table S9). The introgressed SNPs shared among the hybrid zones XD, AZY, and SHB were associated with the gene Rhsim01G0163000 (Table S10), which is implicated in multiple functions, including the regulation of growth (GO: 0040008) and positive regulation of transcription, DNA-templated transcription (GO: 0045893). Additionally, hybrid zones DZ and WD shared SNPs associated with the gene SMG7, which is involved in DNA damage stress response. This observation implies that these hybrid zones may have adapted to similar environmental conditions through this shared genetic characteristic.

4. Discussion

Our results, alongside findings from previous studies, demonstrate that hybridization between Rhododendron spiciferum and R. spinuliferum is bidirectional and asymmetric (Yan et al., 2017; Zheng et al., 2021). The variation in the proportion of F₁ individuals observed in this study compared to prior research (Yan et al., 2017) may be attributed to the more informative loci provided by ddRAD-seq data (Wang et al., 2014) and/or potential sampling biases. Overall, the floral morphology of the putative F₁ hybrids appears to be more similar to that of R. spiciferum than to R. spinuliferum. This observation suggests that pollination may occur more readily between F₁ hybrids and R. spiciferum, potentially contributing to the observed asymmetry in hybridization. Nevertheless, both genetic and morphological data imply that the introgression from R. spinuliferum into the genetic background of R. spiciferum may occur through genetic assimilation in certain hybrid zones. Our study significantly enhances the understanding of the evolutionary history, processes, and potential outcomes of natural hybridization.

4.1. Multiple independent hybrid zones

Our results presented here, along with previous research, confirmed that the formation of R. × duclouxii mainly occurs in its core distribution range located in the central region of Yunnan Province, southwest China, where R. spiciferum and R. spinuliferum are sympatric (Yan et al., 2013, 2017; Zheng et al., 2021). Our gene flow and barrier analyses indicated that gene flow among the sampled hybrid zones was notably limited. Moreover, the NewHybrids results revealed varying compositions and proportions of F₁, F_2, and bidirectional backcross hybrids across different hybrid zones, implying the existence of varying hybridization histories for each hybrid zone. The Kernel density distribution of morphological traits for parents and hybrids varied across hybrid zones, indicating that the parallel evolutionary hybridization events leading to R. × duclouxii occurred independently in different hybrid zones, with minimal influence from inter-zone gene flow. These multiple, independent natural hybrid zones serve as an ideal natural laboratory for exploring processes and evolutionary mechanisms, such as gene flow, natural selection, adaption, and the driving factors underlying natural hybridization (Hewitt, 1988; Abbott et al., 2014; Zieliński et al., 2019; Bolte et al., 2024).

4.2. Factors affecting the population genetic structure of hybrid zones

The overall patterns of morphological traits and population genetic structure of R. × duclouxii across different hybrid zones exhibited similarities to those observed in other reported hybrid taxa within the genera such as Pinus and Senecio (Gao et al., 2012; Walter et al., 2020). The intrinsic characteristics of the species, including genome architecture, gene function, gene interactions and selective pressures, may play significant roles in hybridization processes (Pickup et al., 2019; Scopece et al., 2020). Our results revealed a notable degree of variability in both morphological traits and the population genetic structure of R. × duclouxii among the hybrid zones studied. For instance, significant differences in corolla tube length and style length were observed across the hybrid zones. Furthermore, the composition and proportions of hybrid derivatives, including F₁, F₂, bidirectional backcrosses, and later generation hybrids, also varied across the hybrid zones. These variations in morphological and genetic characteristics provided a valuable opportunity to explore the factors shaping the population genetic structure of hybrid zones.

The interactions between the genetic and environmental factors may give rise to diverse hybrid zones, particularly when the genetic and environmental determinants of reproductive isolation or their effectiveness differ across locations (Seehausen et al., 1997; Taylor and McPhail 2000; Milne et al., 2003). The density and distribution of sympatric parental species can directly influence the history and population genetic structure of a hybrid zone (Burgess et al., 2005; Peters et al., 2017). For example, in hybrid zone AN, the prevalence of R. spinuliferum alongside the sporadic occurrence of R. spiciferum may have resulted in a relatively recently formed hybrid zone characterized by a higher proportion of F₁, a low proportion of F₂, and an absence of observed backcrossing hybrids. Conversely, the genetic composition of the parental species may also account for the variations in the genetic makeup of hybrids (Arnold et al., 2012; Araki and Sota 2021). For instance, hybrids from the western hybrid zone DZ exhibited a distinct genetic composition compared to those from the eastern hybrid zones, reflecting genetic differentiation between the eastern and western populations of R. spinuliferum due to geographical isolation. Furthermore, selection pressures or human activities may also play significant roles in generating differences among hybrid zones (Richards et al., 2016; Emel et al., 2021). In hybrid zone XD, which is surrounded by farmland, planted forests, and roads with intense human activity, the prevalence of R. spiciferum-like chloroplast haplotypes in hybrids, despite the low abundance of R. spiciferum individuals, likely reflects a historical impact of human disturbance on R. spiciferum populations. The absence of F₁ hybrid offspring further suggests a relatively long history of hybridization in this region. Collectively, our findings indicate that factors such as intrinsic characteristics, abundance, genetic composition of parental species, selection pressures, and human activities interact, influence the population genetic structure of hybrid zones.

4.3. Potential evolutionary outcomes for natural hybridization

Factors that influence the population genetic structure and morphology of hybrid zones can lead to a range of evolutionary outcomes, from extinction to hybrid speciation (Abbott et al., 2013; Harrison and Larson 2016; Mandeville et al., 2019; Westram et al., 2021). Hybrid speciation necessitates that hybrids establish an independent and stable genetic lineage while simultaneously achieving reproductive isolation from their parental species, although this isolation may be incomplete in certain instances (Abbott et al., 2013; Arntzen et al., 2017). Our previous research found that the reproductive isolation between R. × duclouxii and each parental species is weak (Yan et al., 2019). Furthermore, R. × duclouxii has not established itself as an independent genetic entity, exhibits morphological variability, and does not occupy a distinct ecological niche separate from its parental species. This suggests that R. × duclouxii is evolutionarily young and has not yet evolved into a distinct species. Instead, it may persist as various categories of hybrid offspring across multiple hybrid zones for an extended duration, similar to observations in Populus and Lissotriton (Zieliński et al., 2019; Bolte et al., 2024). In the present study, many individuals of one parental species within the hybrid zones contained genetic material from the other parental species, and the hybridization appeared to be asymmetric, favoring R. spiciferum. Consequently, the genetic composition and morphology of R. spiciferum may be more significantly altered than that of R. spinuliferum if these patterns persist over the long term (DeVos et al., 2023; Groha and Coopa, 2024). However, we cannot ascertain whether this will affect the species boundary, which requires further investigation.

The potential outcomes of hybridization may vary among the different hybrid zones studied, driven by various factors (Christe et al., 2016; Mandeville et al., 2019; Koci et al., 2020; Zhang et al., 2022). Overall, hybrids in our study exhibited a higher proportion of nuclear and chloroplast genetic components from R. spiciferum compared to R. spinuliferum, although the direction and extent of gene introgression differed across hybrid zones. Notably, greater introgression of loci from R. spiciferum to R. spinuliferum was observed in the hybrid zones AZY and WD, while the reverse pattern occurred in the hybrid zones SHB, DZ, and XD. Most introgressed loci were found to be specific to individual hybrid zones, indicating that each hybrid zone may result in a unique evolutionary outcome. Moreover, variations in floral traits, such as corolla tube length and style length, may influence pollinator preferences, potentially leading to different directions and degrees of hybridization (Anderson et al., 2016; Soares et al., 2021). The seven hybrid zones examined here displayed distinct geographical genetic patterns, divided into an eastern and western group, which may also influence the geographic pattern of the hybridization outcomes (Arntzen et al., 2017; Mandeville et al., 2019). Furthermore, anthropogenic activities, particularly habitat destruction, have significantly impacted the fate of these hybrid zones, with four out of five zones experiencing varying degrees of degradation over a brief four-year period (Fig. S18). Consequently, the potential evolutionary outcomes of each hybrid zone may be influenced by interactions between biotic and abiotic factors (Abbott et al., 2013).

4.4. Significance of hybridization in Rhododendron

Hybridization has long been recognized as a significant driver of rapid speciation (Stryjewski and Sorenson, 2017; Malinsky et al., 2018; Ronco et al., 2021). Rhododendron has undergone rapid evolutionary radiations, giving rise to over 1, 000 species within a relatively short geological time frame (Mo et al., 2022; Xia et al., 2022). Natural hybridization is prevalent within this genus and likely contributes to its species diversification (Ma et al., 2022; Mo et al., 2022). Although direct evidence of hybrid speciation events in Rhododendron has not yet been documented, extensive hybridization and genomic admixture suggest that repeated cycles of isolation, expansion, secondary contact, and hybridization may have facilitated rapid species diversification within Rhododendron subg. Hymenanthes in the Hengduan Mountains (Ma et al., 2022).

Rhododendron species are predominantly outcrossing, displaying considerable diversity in floral morphology, pollinator interactions, and ecological preferences (Shrestha et al., 2018; Song et al., 2019; Ye et al., 2021; Zou et al., 2021). The diversity of corolla shapes and colors observed in Rhododendron likely represents adaptions to attract a variety of pollinators (Berry et al., 2017). The hybrid R. × duclouxii exhibits a range of corolla shapes and colors that differ from those of the parental species, potentially implying a mechanism for promoting the evolution of novel adaptive traits and enhancing trait diversity (Barrera-Guzmán et al., 2018; Yan et al., 2019; Grant and Grant 2020; Ma et al., 2022). Further empirical studies, including pollination experiments and network analyses, are needed to investigate potential shifts in pollinator preferences associated with varying corolla shape and color in hybrids.

5. Conclusions

In this study, we investigated multiple hybrid zones of Rhododendron × duclouxii by integrating morphological, environmental, climatic, and genetic data. Our results further confirmed that hybridization between R. spinuliferum and R. spiciferum is bidirectional and asymmetric, with varying proportions of parental, F₁, F₂, and backcross individuals across different hybrid zones. Significant variations in morphology and genetic composition were observed both between individuals and hybrid zones. Each of the seven hybrid zones appeared to have originated from distinct hybridization events, with minimal or no gene flow between them. Factors such as intrinsic paternal traits, the abundance of parental species, genetic composition, selection pressures, and human activities may have interacted to influence the genetic structure of these hybrid zones. The observed variability in morphology, along with a lack of ecological, morphological, and genetic differentiation from the parental species, indicated that R. × duclouxii consists of hybrid swarms rather than representing a distinct hybrid species. Nevertheless, evidence of both historical and contemporary gene flow into parental species was observed in these hybrid zones, although the direction of gene flow varied between zones. The unique combinations of biotic and abiotic factors specific to each hybrid zone may potentially lead to different evolutionary outcomes. This study enhances our understanding of the ongoing evolutionary processes driven by natural hybridization.

Acknowledgements

We thank Drs. Jie Liu, Guoqian Yang, Linjiang Ye, and Mr. Xinfeng Liao for their assistance during the field collection, lab work and their constructive suggestions on the 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 study was supported by the National Natural Science Foundation of China, China (32300200, 31670213, 32160240); Postdoctoral Directional Training Foundation of Yunnan Province, China (E33O31C261), the Key Basic Research Program of Yunnan Province, China (202101BC070003), and the CAS President's International Fellowship Initiative, China (2024PVA0087). The Royal Botanic Garden Edinburgh is supported by the Rural and Environment Science and Analytical Services Division of the Scottish Government, United Kingdom.

CRediT authorship contribution statement

Wei Zheng: Writing – original draft, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation. Li-Jun Yan: Writing – original draft, Investigation, Funding acquisition. Kevin S. Burgess: Writing – review & editing. Richard I. Milne: Writing – review & editing. Han-Tao Qin: Methodology. Shao-Lin Tan: Methodology. Ya-Huang Luo: Methodology. Jia-Yun Zou: Investigation. Zhi-Qiong Mo: Methodology. Michael Mӧller: Writing – review & editing. Chao-Nan Fu: Writing – review & editing. Lian-Ming Gao: Writing – review & editing, Project administration, Funding acquisition, Conceptualization.

Data accessibility

All newly generated trait data, trnL-F sequences and ddRAD-seq reads are available at the Science Data Bank (DOI: https://doi.org/10.57760/sciencedb.13897).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

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