Evidence for continual hybridization rather than hybrid speciation between Ligularia duciformis and L. paradoxa (Asteraceae)

Sequencing of the three chloroplast intergenic spacers and one nrITS region

Three chloroplast DNA regions were amplified for all the sampled individuals using universal primer pairs for trnL-rpl32 (Shaw et al., 2007), trnQ-5′rps16 (Shaw et al., 2007) and psbA-trnH (Sang, Crawford & Stuessy, 1997; Tate & Simpson, 2003). Polymerase chain reaction (PCR) was performed using a reaction volume of 20 µL, containing 20–60 ng of template DNA, 2.0 µL of 10× PCR buffer (Mg²⁺ free), 1.0 µL of MgCl₂ (25 mM), 1.0 µL of dNTP (10 µM each), 1.0 µL of BSA (20 mg/mL), 0.3 µL of each primer and 0.15 µL of 1.5 units of Taq polymerase (Takara, Dalian, China). Amplification proceeded as follows: an initial 5 min at 95 °C for denaturation; 30 cycles of 45 s at 94 °C; 45 s at 53 °C; and 50 s at 65 °C; followed by a final extension of 7 min at 65 °C. The nrITS sequences of all sampled individuals (107) were amplified using the primers ITS4 and ITS5 (White et al., 1990) using the PCR procedure of Yu, Kuroda & Gong (2014). PCR products were purified by electrophoresis in a 1.2% agarose gel, from which the products were recovered with an E.Z.N.A. Gel Extraction Kit (Omega, Guangzhou, China). All accessions were subjected to sequencing with the amplification primers in an ABI 3700 DNA sequencer with a BigDye Terminator Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA).

For the nrITS sequences, direct sequencing was successful for 51.42% of L. duciformis individuals (18/35), 85% of L. paradoxa (34/40) individuals and 100% of L. lamarum individuals. For individuals with chimeric or unreadable peaks in the chromatograms, including each parent-specific and putative hybrid-specific ITS sequence, positive clones with accurate inserts were confirmed using colony PCR (Yu, Kuroda & Gong, 2014). Two to eleven clones were sequenced using ITS4 and ITS5 primers. Both strands of the ITS clones were sequenced. The number of cloned sequences for each accession was listed in Table S2. The DNA sequences generated in this study are available in GenBank with accession numbers: KY307306–KY307785.

Network analyses based on nrITS and cpDNA regions

All sequences were assembled and aligned using SeqMan (DNAstar 7.1, DNASTAR Inc., Madison, WI, USA) and adjusted manually with BioEdit 7.0.4.1 (Hall, 1999). Variable sites and haplotypes were obtained using the DnaSP 5.0 program (Rozas et al., 2003). Haplotype networks for the nrITS sequence and chloroplast DNA fragments were constructed using TCS 1.21 (Clement, Posada & Crandall, 2000) with parsimony probability set to 98%.

SSR genotyping

We selected 11 nuclear SSR markers from previous studies (Ahrens & James, 2013; Liu et al., 2004; Mao, Li & Wang, 2009) that could be successfully applied to sampled individuals in this study, among which loci with variable flanking regions were excluded to avoid null alleles. Moreover, individuals with more than two missing loci were removed from further analyses and a total of 89 individuals (L. duciformis, L. paradoxa and L. ×maoniushanensis) from two locations were used for SSR genotyping. PCR was optimally performed in a reaction volume of 25 µL containing 20–60 ng of DNA, 2.5 µL of 10× PCR buffer, 2 µL of MgCl₂ (25 mM), 1 µL of dNTPs (10 mM), 0.5 µL of each primer, 1.0 µL of BSA (20 mg/mL), and 0.15 µL of Taq DNA polymerase (5 unit/µL; Takara, Shiga, Japan). PCR was performed in a thermocycler under the following conditions: 1 cycle at 94 °C for 5 min, which was followed by 35 cycles at 94 °C for 60 s, 53–57 °C for 60 s and 72 °C for 90 s, followed by 1 cycle at 72 °C for 8 min. The PCR products were electrophoretically detected on a 2% agarose gel. We then performed preliminary screening for the SSR loci. An individual was considered null (no amplification) and treated as missing data if amplification failed more than two times at a locus. Individuals containing three or more missing data were excluded from the subsequent analyses.

SSR data analyses

Raw dataset editing and formatting were executed in GenAlEx 6.3 (Peakall & Smouse, 2006). A Bayesian approach was implemented in STRUCTURE 2.3.4 (Pritchard, Stephens & Donnelly, 2000) to identify population structure and potentially admixed individuals from multi-locus data. This Bayesian clustering algorithm identifies an optimal number K of genetic clusters of sampled individuals and simultaneously assigned the individuals to the genetic clusters by calculating the posterior probability of cluster membership. To determine the optimal number of clusters for putative parents and hybrids, 10 independent runs for each K value (K = 1–6) were executed with 10⁶ Markov Chain Monte Carlo iterations and a burn-in period of 100,000. The optimal K value was estimated using the mean value (ΔK) method (Evanno, Regnaut & Goudet, 2005). We adopted a threshold value of membership coefficient (q_i) to identify purebred parental individuals (0 < q_i < 0.10 or 0.90 < q_i < 1.0) and hybrids (0.10 < q_i < 0.90; Hoban et al., 2009; Vähä & Primmer, 2006). The software NewHybrids 1.1 (Anderson & Thompson, 2002) was used to assign sampled individuals to six genotype categories (pure species A, pure species B, F₁ hybrid, F₂ hybrid, F₁ backcross to pure species A, and F₁ backcross to pure species B) based on the 11 loci. The posterior probability (PP) of categorical membership for each individual was computed using a Bayesian approach. We use PP ≥ 0.90 to designate purebreds. For hybrids, PP ≥ 0.50 to assign a genotypic class; if PP ≤ 0.50, the individual was not assigned to any category.

In STRUCTURE and NewHybrids analyses, given that we wanted to ask whether the three taxa of Ligularia are part of one large panmictic species complex, we jointly analysed the SSR data of individuals from the two sampling sites to effectively exclude substantial genetic differentiation among the populations from different regions; to assess sensitivity of analyses to this procedure, samples from the two sampling sites were also analysed separately.

Introgressive patterns and genomic cline analyses of SSR loci were performed using the INTROGRESS package (Gompert & Buerkle, 2010) in R 3.1.2 (R Core Team, 2012). This package identifies each locus that deviates from expectations of neutral introgression using the log likelihood ratio given the observed data, from which a P value for the significance test is calculated. Significant deviations from neutral expectations for genomic clines were adjusted for multiple comparisons using the false discovery rate (FDR; Benjamini & Hochberg, 1995). We then used the clines.plot function to plot the genomic cline of each locus. Additionally, we estimated the hybrid index for the admixed individuals using functions est.h and mk.image to provide a clear visualization of the variation in introgression and ancestry across loci and for each of the admixed individuals (Gompert & Buerkle, 2010). Ligularia duciformis, L. paradoxa and L. ×maoniushanensis are recognized as homozygotes (Ad/Ad, Ap/Ap) or heterozygote (Ad/Ap) genotypes, respectively.

ISSR screen and analyses

We screened nine polymorphic ISSR marker primers (University of British Columbia: UBC807, UBC808, UBC811, UBC818, UBC828, UBC845, UBC849, UBC850, and UBC857) from a set of 100 ISSR for the sampled individuals of the putative parental species and hybrids. These markers were selected for further amplification. Final PCR amplifications were performed in a reaction volume of 20 µL, containing 30–50 ng of template DNA, 2.0 µL of 10× buffer, 1.6 µL of MgCl₂ (25 mmol/L), 1.2 µL of dNTPs (10 mmol/L), 0.3 µL of each primer (10 µmol/L), and 0.15 µL of Taq DNA polymerase (5 U/µL; Takara, Shiga, Japan). The PCR amplification conditions were as follows: 1 cycle at 94 °C for 7 min; 35 cycles at 94 °C for 45 s, 52/53 °C for 1 min and 72 °C for 2 min; followed by 1 cycle at 72 °C for 7 min. Because the possibility of L. lamarum participating in hybridization was excluded according to nrITS and cpDNA data (see Results), only 89 individuals (including 35 L. duciformis, 19 L. ×maoniushanensis, and 35 L. paradoxa individuals) were amplified in each PCR run. For each maker, if the PCR failed for more than five individuals, the marker was excluded from the subsequent analyses. The ISSR bands observed after electrophoresis on an agarose gel were manually recorded as binary characters (0 for absence, 1 for presence and “-” for failed amplification). Only polymorphic band data were used for downstream analysis, since monomorphic bands cannot discern relationships among individuals. The analyses of ISSR data were performed in NewHybrids 1.1 following Yu, Kuroda & Gong (2011).

nrITS sequence analyses

Among 107 sampled individuals, 66 samples were directly sequenced and the remaining 41 samples (including 19 L. ×maoniushanensis, 16 L. duciformis and 6 L. paradoxa) were sequenced by cloning (Table S2). The aligned nrITS sequences had a length of 704 bp, of which 68 variable sites and 79 haplotypes were identified (Table 2). Sympatric L. lamarum, containing 18 specific variable sites, possessed specific haplotype H20 and could be distinguished clearly from other taxa. Thirty-six species-specific variable sites were used to identify the occurrence of hybridization between L. duciformis and L. paradoxa. As listed in Table 2, 65.6% sequences of L. duciformis have the following haplotypes: H1, H2, H3, H4, H5, H10, H14, H60, H62. Likewise, 84.8% sequences of L. paradoxa have haplotypes: H27, H43, H50, H54, H55, H57. For L. ×maoniushanensis, 63.2% sequences possessed haplotypes from L. duciformis and L. paradoxa (Table 2).

Haplotype networks for L. duciformis, L. paradoxa, L. ×maoniushanensis and L. lamarum were generated to illustrate their relationships (Fig. 2). The nrITS network contained three major subnetworks, where individuals of L. duciformis and L. paradoxa constituted two major subnetworks (here termed parts A and B), respectively, and sympatric L. lamarum individuals formed a third subnetwork (termed part C) located at the edge of the nrITS network (Fig. 2). Haplotypes derived from the cloned sequences of L. ×maoniushanensis were scattered among parts (A) and (B), which corresponded to L. duciformis and L. paradoxa, respectively, among which some haplotypes were unique to L. ×maoniushanensis. These unique haplotypes may have resulted from unsampled polymorphism in the parental species, mutation in hybrid individuals, and/or genetic recombination in the hybrids.

Analyses of the combined cpDNA sequences

The aligned and trimmed length of the amplified psbA-trnH, trnL-rpl32 and trnQ-5′rps16 regions were 453 bp, 916 bp and 905 bp, respectively. The concatenated length for three fragments was 2,275 bp, containing 19 variable sites in total (Table S3). Among loci, trnL-rpl32 had the most variable sites (52.6%). Ligularia ×maoniushanensis shared 31.6% variable sites with L. duciformis and 63.2% with L. paradoxa. Six haplotypes were identified among all the individuals (H1-H6; Table 2 and Table S3). A haplotype network was inferred based on the three combined cpDNA fragments of L. duciformis, L. paradoxa, L. ×maoniushanensis and L. lamarum (Fig. 3). Ligularia duciformis had two haplotypes (H1 and H2) and L. paradoxa exclusively had the H4 haplotype. Individuals of L. ×maoniushanensis at Mt. Maoniu shared haplotype H4 with L. paradoxa. However, at the Heihai Lake locality, 60% of L. ×maoniushanensis individuals shared haplotypes H1 with L. duciformis, 30% shared haplotype H4 with L. paradoxa, and 10% (sample HM5) had a unique H5 haplotype. Haplotypes H3 and H6 were unique to L. lamarum; each population of L. lamarum had a single haplotype (Fig. 3; Table 2).

SSR profiles

Eleven SSR loci were genotyped for L. duciformis, L. ×maoniushanensis, and L. paradoxa, with 69, 53, and 72 alleles identified, respectively. The average observed heterozygosity was 0.624, 0.696 and 0.581 for L. duciformis, L. ×maoniushanensis, and L. paradoxa, respectively.

STRUCTURE analyses

In the STRUCTURE analysis, SSR datasets from the two locations were analysed jointly. The highest ΔK value was obtained with K = 2 for the putative parents and hybrids (Table S4; Fig. S1), indicating the presence of two genetic clusters. For L. duciformis, all individuals (15) from Mt. Maoniu and most individuals (18/20) from Heihai Lake were assigned to one purebred cluster (q_i > 0.90; Fig. 4A). For L. paradoxa, all individuals (18) from Mt. Maoniu and most individuals (17/18) from Heihai Lake were assigned to the other purebred cluster (q_i > 0.90). However, a few individuals that were morphologically identified as L. duciformis or L. paradoxa were assigned to the mixed cluster with high probability (q_i (HD18) = 0.828, q_i (HD20) = 0.885, q_i (HP15) = 0.880). For L. ×maoniushanensis, eight individuals from Mt. Maoniu and nine from Heihai Lake were inferred to have a genetic composition deriving from both L. duciformis and L. paradoxa, and were unassigned to either cluster (0.10 < q_i < 0.90). The remaining two individuals of L. ×maoniushanensis were assigned to the L. paradoxa cluster (q_i > 0.90).

The STRUCTURE analyses for SSR data treating each of the two sampling sites individually generated very similar results. STRUCTURE analyses of the three taxa from Mt. Maoniu indicated that the optimal number of clusters was K = 2 (Table S5 and Fig. S2). Fifteen individuals of L. duciformis and 18 individuals of L. paradoxa were assigned to the two purebred clusters (q_i > 0.90), and nine individuals of L. ×maoniushanensis were assigned to the mixed cluster (0.10 < q_i < 0.90) (Fig. 4B). Two genetic clusters (K = 2) was also suggested for taxa from Heihai Lake (Table S6; Fig. S3). Twenty individuals of L. duciformis and 17 individuals of L. paradoxa were allocated to the two purebred clusters (Fig. 4C). Most individuals of L. ×maoniushanensis (8/10) showed an estimated membership (q_i) ranging from 0.16 to 0.84, and assigned to the admixed group, whereas the two remaining individuals were considered to be L. paradoxa.

In comparing the results of joint with separate analyses, we found some individuals with values near the critical value were assign to mixed class in joint analyses, e.g., q_i (HD18) = 0.828, q_i (HD20) = 0.885, however, they were assign to parental class (q_i (HD18) = 0.924, q_i (HD20) = 0.918, respectively) in separate analyses.

NewHybrids analyses

NewHybrids analyses were conducted for 89 individuals from the two sampling areas; the joint analyses (Fig. S4) were slightly different from the separate analyse for each sampling area (Fig. 5). For the Mt. Maoniu population, with the exception of one individual (MD15), all the individuals of L. duciformis and L. paradoxa were assigned to one parent with high PP (>0.94; Fig. 5A). Individual MD15 of L. duciformis may be a later hybrid progeny (PP < 0.90). Individuals of L. ×maoniushanensis were F₂ hybrids with relatively high PP (>0.80). For the Heihai Lake population, 18 out of 20 individuals of L. duciformis and 13 out of 17 individuals of L. paradoxa were assigned to their respective parental class with high PP (>0.90; Fig. 5B). The other two individuals of L. duciformis were not assigned to L. duciformis (PP < 0.90). Furthermore, the remaining four individuals of L. paradoxa were not assigned to any class (PP < 0.50), and these individuals might be backcrosses to L. paradoxa. For L. ×maoniushanensis, eight individuals were F₂ hybrids with a posterior probability higher than 0.50, whereas the other two individuals were not assigned to any class (PP < 0.50) which might result from backcrossing to L. paradoxa.

INTROGRESSION analyses

Figures of genomic clines (Figs. 6A–6K) were generated based on SSR data. This method unmasked marked heterogeneity in locus-specific patterns of introgression, and genomic clines of introgressed alleles can be identified. Eight loci (L38, L40, L41, L64, S4, S23, S24 and V40) showed significant deviations after FDR correction from expectation value based on a null model of neutral introgression (P < 0.021), indicating that introgression was detected in these loci (Fig. 6). For the six loci (L40, L64, S4, S23, S24, and V40), the homozygote and heterozygote genomic clines were steeper than predicted by the neutral model. Genotypic variation at locus L64 was also involved in steep transitions between parental homozygotes, but very few heterozygotes were observed.

There was some heterogeneity among loci in the patterns of introgression between species, as well as some similarities, which was visually evident (Fig. 7A). On each locus, permutations of three genotypes, Ad/Ad, Ad/Ap and Ap/Ap, were observed. Genotypes with more Ad/Ap were mainly confined to the loci L40, L77, S2, S4, S23, S24, Sm083 and V40. Some individuals of L. duciformis shared genetic composition with L. paradoxa and L. ×maoniushanensis at most loci, except for L40, L64 and S24, and L. paradoxa shared genetic composition with most loci, except for L41. Moreover, genomic cline analysis revealed a gradually increasing pattern of introgression among L. ×maoniushanensis individuals in the proportion of L. paradoxa ancestry estimated (Fig. 7B).

ISSR profiles

Nine polymorphic ISSR markers were screened for the sampled individuals of L. duciformis, L. ×maoniushanensis, and L. paradoxa. Among these individuals, L. duciformis had 12 unique polymorphic loci and L. paradoxa had nine, whereas L. ×maoniushanensis did not contain any unique loci. Five loci that were shared among three Ligularia taxa were excluded.

NewHybrids analyses

NewHybrids analyses were separately implemented for each sampling site. For individuals from Mt. Maoniu, with the exception of one individual (MD3), all the individuals of the putative parents were assigned to L. duciformis and L. paradoxa with high PP values (>0.939). Moreover, all the putative hybrids (L. ×maoniushanensis) exhibited the highest PP (>0.749) for the F₂ hybrids compared to other genotype classes. For the Heihai Lake population, most individuals of L. duciformis (18/20) and L. paradoxa (12/17) were assigned to genotypes consistent with morphological identification (Fig. 8B). Two individuals of L. duciformis and one individual of L. paradoxa were not assigned to a parental class; this was with low PP (<0.90). The remaining four individuals of L. paradoxa were not assigned to any class (PP < 0.50). For L. ×maoniushanensis, nine individuals were assigned to F₂ hybrids with the PP > 0.50, and the remaining one was assigned to backcross to L. paradoxa.

Evidence for natural hybridization and introgression: nrITS and cpDNA markers

In this study, the nrITS data support the hypothesis that morphologically intermediate individuals are primarily hybrid descendants between L. duciformis and L. paradoxa and that the sympatric taxon L. lamrum does not participate in hybridization. In the nrITS analyse, some morphologically identified parental individuals showed double peaks of nrITS; the cloned nrITS sequences of these individuals suggest that these individuals may be the progenies of backcrossing. Backcrossing has also been suggested previously to cause frequent gene flow between Ligularia species, giving rise to complex relationships between them (Yu, Kuroda & Gong, 2011). In our results, some individuals had unique haplotypes; such variants might be caused by hybridization between parents containing different ITS sequences (Baldwin et al., 1995; Sang, Crawford & Stuessy, 1995; Li et al., 2014) or independent evolution of haplotypes through mutation exceeding the rate of concerted evolution or pseudogenization among ITS copies (Feliner & Rosselló, 2007; Kosnar et al., 2012).

Since cpDNA is reported to be maternally inherited in Ligularia (Zhang & Liu, 2003), analyses of chloroplast sequences can reveal the direction of hybridization. We observed that L. paradoxa was the only maternal species in Mt. Maoniu, whereas L. duciformis and L. paradoxa were both maternal species in Heihai Lake. Consequently, we deduce that the hybridization was primarily unidirectional (L. duciformis ♂ × L. paradoxa ♀) in Mt. Maoniu and bidirectional in Heihai Lake. Moreover, the hybrid combination of L. duciformis ♀ × L. paradoxa ♂ was predominant at Heihai Lake. These results congruent with those of Pan et al. (2008) in Mt. Maoniu where natural hybridization between L. duciformis and L. paradoxa was exclusively unidirectional; however, the results also suggest that hybridization in Heihai Lake is asymmetrical, with L. duciformis as the primary maternal parent.

Given potentially influential factors such as habitat, population size, and gametophytic-sporophytic interactions during fertilization or organelle-nuclear gene interactions, asymmetric hybridization barriers in plants are common (Tiffin, Olson & Moyle, 2001; Turelli & Moyle, 2007). The phenomenon of asymmetrical hybridization has been reported in many taxa, such as Rhododendron (Ma et al., 2010; Yan, Gao & Li, 2013; Zha, Milne & Sun, 2009) and Melastoma (Liu et al., 2014). According to our field observation, individuals of L. duciformis are more abundantly distributed than any other sympatric Ligularia species in both hybrid swarms, possibly swamping available conspecific individuals. In addition, Mt. Maoniu and Heihai Lake have different habitat conditions. Mount Maoniu, possessing abundant wet soil under moist forest, appears to be more mesic habitat than Heihai Lake, which has less water availability among arid rock crevices where the plants occur. In the more mesic habitats of Mt. Maoniu, L. duciformis serves as paternal species, but in the arid habitats at Heihai Lake, L. duciformis and L. paradoxa are both maternal species. These different population dynamics may result from habitat or climatic factors, yet little is known of habitat requirements in these species and further studies are needed to clarify the most important factors in the present case.

Evidence for natural hybridization and introgression: SSR and ISSR markers

Given the evidence for substantial interspecific gene flow accrued in this study, it could be alternatively hypothesized that these Ligularia species distributed at the two study sites comprise a single large, randomly breeding complex to be treated as a single genetically distinct species. This hypothesis can be directly tested by the STRUCTURE analyses of the SSR data of all individuals from the two locations. The random-breeding complex hypothesis would be supported if L. duciformis and L. paradoxa cluster by sampling sites rather than by species, since gene flow is expected among individuals in sympatry. STRUCTURE analyses indicated that individuals from the two sampling sites cluster by species rather than by sampling sites, indicating the presence of genetic differentiation among these species; hence considering the parents as separate species is justified given the loci examined here.

Although overall results were congruent among datasets and analyses, we found that a few individuals were assigned to different genotype categories for purebred and hybrids in STRUCTURE and NewHybrids. This can be explained as a result of marker choice; accurate identification of genetic categories using molecular methods depends on the markers selected and the degree of differentiation among species (Anderson & Thompson, 2002; Vähä & Primmer, 2006). We employed a default threshold (q_i, 0.90) for discriminating purebred and hybrid specimens in our study, which in combination with differing evolutionary patterns among markers may lead to discrepancies among individual data partitions.

We also tested the hybrid origin of L. ×maoniushanensis and its genetic structure. In STRUCTURE analyses, L. duciformis and L. paradoxa formed two distinct genetic clusters; many individuals of L. ×maoniushanensis showed genetic admixture between the two clusters, implying some F₁ hybrids. Based on the NewHybrids analyses of the SSR markers, most individuals of L. ×maoniushanensis were assigned to as F₂ hybrids or backcrosses. The thresholds used for assigning individuals to the different genetic categories are different between STRUCTURE and NewHybrids. It is possible that different hybrid categories were suggested by the two methods partly for this reason. In addition, these two methods have different strengths. STRUCTURE is more efficient for evaluating the presence of hybrids in wild populations (Marie, Bernatchez & Garant, 2011). However, given the assumption of two parental categories, NewHybrids assigns explicit hybrid categories, which might be likely to show higher assignment accuracy than those obtained using STRUCTURE (Marie, Bernatchez & Garant, 2011). The NewHybrids results indicate that F₂ hybrids and backcrosses are more prevalent than F₁ hybrids. The occurrence of F₂ hybrids and the preponderance of backcrosses to both parental species are consistent with recent hybridization and potentially ongoing hybridization between the two species and among hybrid derivatives. Although SSR and ISSR markers had different resolution for distinguishing among the genotypes of the examined samples, a similar result was obtained when SSR results were compared with the ISSR results in NewHybrids analyses. Therefore, these results provide compelling evidence for the hybrid origin of L. ×maoniushanensis.

Meanwhile, based on the STRUCTURE results, a few parental individuals of L. duciformis and L. paradoxa were genetically assigned to L. ×maoniushanensis, including HD18 and HP15, and a few individuals of L. ×maoniushanensis were genetically assigned as L. paradoxa, such as MM8 (from Mt. Moniu) and HM9 (from Heihai Lake). These individuals, for which morphology and molecular markers were incongruent, may result from frequent backcrosses to L. paradoxa, or the insufficiency of SSR markers to identify the individuals. In the NewHybrids analyses, some parental individuals with PP <0.90 are morphologically identified as L. duciformis and L. paradoxa, but they may be later generation hybrids according to the results of NewHybrids. Some hybrids that morphologically resemble parental species may result from continual backcrosses to the parental individuals. For example, one L. duciformis individual and one L. paradoxa individual in Heihai Lake were assigned to F₂ or backcross hybrids; this indicates the existence of bidirectional introgressive hybridization between L. duciformis and L. paradoxa.

Results from INTROGRESS further confirmed the presence of extensive introgressive hybridization between L. duciformis and L. paradoxa. Genomic cline analyses revealed a diversity of introgression patterns among loci. Clines for the majority of markers were inconsistent with neutral introgression in hybrids (8/11; Fig. 6). Excess and deficits of the three genotypes (Fig. 7) presented different patterns among loci. This may be consistent with the action of selection at these loci; the potential presence of non-neutral introgression in these populations requires further research.

From infertile seeds to repetitive backcrosses or hybrids

Seed germination rates in L. ×maoniushanensis have been examined previously; L. ×maoniushanensis was found to completely lack viable seeds (Pan et al., 2008). However, this experiment does not completely exclude the existence of fertile seeds from hybrids. It is possible that rare fertile seeds are occasionally generated from repeated hybridization events. The lower fertility of the initial hybrids compared to parental individuals could imply that hybrid progenies play a minor role in the evolution of a given species complex (Arnold, 1997). Although the fertility of the hybrids may limit the production of hybrid individuals, as initial hybrids repeatedly form, the opportunity of producing later hybrid generations increases, potentially restoring fertility in late-generational hybrids. A classic case comes from a repetitive cross experiment between Helianthus annuus and H. petiolaris (Asteraceae) which resulted in an increase of fertile seed (Heiser et al., 1969). Although low levels of fertility and viability often occur in initial hybrid generations, the fertility for late-generational hybrids could result in a stabilized hybrid species (Rieseberg, 1991). Therefore, the production of some fertile hybrid progenies is possible during hybridization between L. duciformis and L. paradoxa given the frequent occurrence of hybrid individuals.

Is L. ×maoniushanensis a hybrid species or a hybrid swarm?

Given the individuals we sampled, the status of L. ×maoniushanensis is unclear. A true hybrid species (an independent lineage of hybrid origin) should show patterns of continuous breeding among its populations at least beyond the F₂ generation (Zhou, Shi & Wu, 2005). Based on our observations, the present number of L. ×maoniushanensis individuals in the Mt. Maoniu is smaller than a decade ago. This phenomenon could be explained by the low fertility of hybrids and a preponderance of early-generational hybrids. Frequent hybridization and repetitive formation of L. ×maoniushanensis in these locations does not appear to have significantly increase the population size. However, L. ×maoniushanensis maintains a sympatric distribution with the parental species rather than occupying a distinct habitat or niche, and frequent backcrossing may inhibit the formation of reproductive isolation with parental species. The evolutionary species concept (Simpson, 1951; Blackwelder, 1962; Soltis & Soltis, 2009; Wheeler, 2000; Wiley, 1978) delimits as species those lineages evolve separately and have their own evolutionary tendencies. Under this concept, we conclude that L. ×maoniushanensis is repeatedly generated from the parental species, and since it is not independently evolving from its parental taxa it should not be considered as a new species. However, the intrinsic factors driving the directionality and asymmetrical contributions of parental genomic material, a more precise determination of hybrid generations, and the affinity of assigned individuals should be the topics of further investigation.