Distribution range of C. pinnata is concentrated in Beijing and adjacent areas in Tianjin and Hebei Province, only with a few records in Liaoning Province (explained in Taxonomic Treatment). From 2017 to 2019, we surveyed C. pinnata over its entire distribution range and found multiple distribution sites in Beijing and Liaoning Province (Table 1). In all its distribution sites, both of its putative parental species are present (Table 2 and Supplementary Figure S1). We recorded GPS coordinates and altitude of each population, and collected leaf tissues as well as specimens for flow cytometry experiments, phylogenomic studies, and morphological analyses. Individuals of C. pinnata are rare, often with a few plants in each site. So, we collected all the individuals in each site. In total, we gathered 27 C. pinnata individuals and one dubious plant from eight collecting sites in Beijing and Liaoning Province. A dubious individual was collected in Baihuashan (BHS), Beijing. It was not in flower when collected. The plant was creeping with pinnate leaves, which was similar to C. pinnata. However, the leaflets were smaller than typical C. pinnata and very similar to those of C. brevicaudata. We included this dubious individual in our phylogenomic study to determine its identity.

From all the eight collecting sites, we gathered 94 individuals of C. brevicaudata. We also collected 91 individuals of C. tubulosa from five of the eight sites, and 49 individuals of C. heracleifolia from the rest of the three sites. The individuals were randomly collected in each site, and the distance between individuals is more than 20 meters. Furthermore, five other Clematis species (C. hexapetala Pall., C. acerifolia Maxim., C. ochotensis (Pall.) Poir, C. intricata Bunge, and C. fruticosa Turcz.) commonly found in northern China were also included in the phylogenomic analysis as outgroups. Voucher specimens were deposited in the Herbarium of Beijing Forestry University (BJFC, Supplementary Table S1).

We used a flow cytometry (FCM) method (Doležel et al., 2007; Bourge et al., 2018) to determine whether C. pinnata is homoploid with its putative parents. We used samples from Jiufeng forest park to check if C. pinnata is homoploid with the putative parents in the same site. Because only C. brevicaudata and C. tubulosa are present in Jiufeng population, we checked C. heracleifolia using samples from Yanqing district.

The FCM measurement was carried out in the Key Laboratory of Photobiology of Institute of Botany, the Chinese Academy of Sciences. Preparation of dried leaf samples of the four taxa basically followed a standard two-step method to isolate plant cell nuclei (Qu et al., 2018). For each sample, about 20–30 mg of dried leaf was chopped in 1 mL of LB01 lysis buffer (15 mM Tris, 2 mM EDTA-Na₂, 0.5 mM spermine tetrahydrochloride, 80 mM KCl, 20 mM NaCl, 0.1% Triton X-100, 15 mM β-mercaptoethanol, pH 7.5) to release cell contents. The resulting culture was gently pipetted and filtered through a 400-mesh screen to remove cell debris. The samples were then stained with 50 μL 1 mg⋅mL^–1 propidium iodide (PI) and 50 μL 1 mg⋅mL^–1 RNase A in an ice bath for 10 min before being analyzed using a MoFlo XDP flow cytometer (Beckman Coulter Inc., United States). We measured about 5000 nuclei from each sample (individual) and performed three replicates (three individuals) for each species. PI conjugated to Alexa Fluor 488 (Molecular Probes) was excited with a 488-nm argon-ion laser and the fluorescence was detected using a microscope equipped with a 625/26-nm HQ bandpass filter. Mean channel positions and coefficients of variation (CVs) of the G0/G1 peaks were calculated using Summit 5.2 software (Beckman Coulter Inc., United States). We used C. brevicaudata as the external standard reference for our flow cytometry measurements.

To examine the possible hybrid origin of C. pinnata, we generated transcriptomes from samples of the eight collecting sites. Fresh leaf materials of the Clematis samples were collected from field and quickly deposited in liquid nitrogen. For C. pinnata, leaf tissues of all the individuals (from 1 to 10 individuals per site, Table 2) in each site were mixed together, respectively. Population transcriptome data of C. brevicaudata and C. heracleifolia/C. tubulosa were represented by five individuals for each species in each site. Five other Clematis species (C. ochotensis, C. intricata, C. hexapetala, C. fruticosa, and C. acerifolia, one individual for each species) were included for transcriptomes analysis, and Anemoclema glaucifolium was chosen as the outgroup.

We extracted total RNA from sampled leaves by using TRIzol Reagent (Invitrogen, Thermo Fisher Scientific, Shanghai, China) and sent the RNA samples to Biomarker Technologies¹ for cDNA library preparation and Illumina 2 × 150 bp paired-end sequencing. About 6 Gb of raw reads were obtained for each sampled species. The raw reads for the 30 newly generated transcriptomes were deposited in the NCBI Sequence Read Archive under BioProject PRJNA657443 (Supplementary Table S1). Then, they were filtered and trimmed using Trimmomatic v. 0.39 (Bolger et al., 2014). The clean RNA-Seq reads were de novo assembled using Trinity v. 2.5.1 (Grabherr et al., 2011) with default parameters. Because the assembly had a large number of redundant transcripts, we kept the longest isoforms of the related contigs by using the “get_longest_isoform_seq_per_trinity_gene.pl” utility in Trinity to let each unigene represent a collection of expressed sequences that apparently came from the same transcription locus. We used CD-HIT (Li and Godzik, 2006; Fu et al., 2012) to remove redundant sequences from the unigene and TransDecoder v 5.0 to predict protein-coding regions (Haas et al., 2013).

Single-copy orthologous nuclear genes (SCOGs) can be determined and screened out by reference-based methods (Zhao et al., 2013; Zeng et al., 2014, 2017; Feng et al., 2019) and all-against-all alignment methods without references (Sveinsson et al., 2014; Nishikawa et al., 2015; Majeed et al., 2019). The second strategy takes much more computational time, but can find more complete single-copy orthologous dataset than the former strategy, especially when there are no close genomic references available (Nishikawa et al., 2015). In this study, our method followed that of Nishikawa et al. (2015), which was based on all-against-all alignment with BLASTP as implemented in Proteinortho v. 6 (Lechner et al., 2011). The definition of an orthologous relationship between proteins was: e ≤ 1 × 10⁵, identity ≥ 25%, and alignment coverage ≥ 50%. The extracted SCOGs were aligned using MAFFT (Katoh et al., 2005), and each alignment was processed with a Python script² to remove ambiguous aligned positions.

For the phylogenetic analysis of the SCOGs data, both coalescent- and concatenation-based methods were applied. The coalescent-based method was conducted by the Accurate Species Tree Algorithm v. 4.4.4 (ASTRAL; Mirarab et al., 2014). Single-gene ML tree was reconstructed in RAxML v. 8.1.17 (Stamatakis, 2014) under the GTR + G model as suggested in the manual. Then, we used TreeShrink (Mai and Mirarab, 2018) to test whether long branch attraction caused incorrect phylogenetic relationship and then deleted all the trees with suspicious patterns of branch length. All the remaining ML trees were then used as input trees for ASTRAL. For the concatenation method, the aligned SCOGs sequences were linked using Geneious Prime 2019 (Kearse et al., 2012). Then, we used a ML method for phylogeny reconstruction as described above. In view that the presence of hybrid taxa may cause problems in phylogenetic inference (Xu, 2000), we exclude putative hybrid taxa (C. pinnata and C. ochotensis inferred by Phylonet analysis, see below) and re-do the above mentioned phylogenetic analysis to confirm the backbone of Clematis phylogeny.

Because complex evolutionary scenarios such as hybridization events do not proceed in tree-like manners (Huson and Bryant, 2006), we used the NeighborNet method implemented in SplitsTree 4.11.3 (Huson and Bryant, 2006) to reconstruct phylogenetic networks for the concatenated dataset. We excluded insertions/deletions (indels) and used the K2P model (Kimura, 1980) for distance analysis and support values at each node were estimated by running 1000 bootstrap replicates.

Then, we analyzed potential hybridization of C. pinnata using pseudolikelihood approach with PhyloNet (Wen et al., 2018), a software that detects phylogenetic networks based on the multispecies coalescent model, to directly infer the hybridization process. Because short gene alignments could cause random errors in phylogenetic inference in PhyloNet analysis (McLean et al., 2019), we selected gene alignments that were at least 500 bp in length for analysis. We divided SCOGs alignments into eight population datasets because PhyloNet cannot test too many hybridization events in a phylogenetic tree due to computational restrictions (Wen et al., 2018; Morales-Briones et al., 2021). Each dataset contains nine taxa: C. pinnata, both its putative parents, five other Clematis species, and Anemoclema as the outgroup. All the separate SCOGs phylogenies were inferred using an ML method as described above. We adopted the maximum pseudolikelihood method to model incomplete lineage sorting and gene flow using individual gene trees with the command InferNetwork_MPL (Yu and Nakhleh, 2015). Each Network search was allowed 0–4 reticulations and the log likelihood score for each network was also inferred. The best number of hybridization events was selected by plotting the likelihood scores. A sharp likelihood score increase is expected until it reaches the best number, after which the score increases slowly. For the dubious individual in BHS population, we also conducted another PhyloNet analysis by excluding another putative hybrid taxa C. ochotensis for comparison.

To further verify hybrid origin of C. pinnata and possible gene flow between the putative progenitors, we used an alternative method similar to D-statistic (or ABBA-BABA tests) analysis (Patterson et al., 2012) for comparing the results with the PhyloNet analysis. In this study, we used HyDe v. 0.4.1a (Blischak et al., 2018) to check for hybridization and gene flow. This method considered both hybridization and coalescence in a unified framework and can be used to fast detect both current and ancient hybridization events using genomic SNP data (Kubatko and Chifman, 2019). It assesses statistical significance of hybridization by testing for hybridization among all possible triplet combinations of the sampled species. Since HyDe uses single nucleotide polymorphism (SNP) data, the SCOG alignments were concatenated to construct a super SNPs (complete SNP data, hereafter) matrix using Phyutility v. 2.2 (Smith and Dunn, 2008) with removing the missing sites. We then used HyDe with a Python script, run_hyde.py (Blischak et al., 2018) to test possible hybrid origin of C. pinnata from each of the eight collection sites. The software assumed that each one of the three tested species (C. pinnata, C. brevicaudata, C. heracleifolia/C. tubulosa) is of a hybrid origin of the other two species. Thus, in total, 169 possible hybrid triplet combinations in each collecting site were tested.

Linkage disequilibrium (LD) may bias the results when analyzing genomic SNP data (Nielsen and Signorovitch, 2003; Malomane et al., 2018). Studies have shown that using LD-pruned data can effectively correct such biases (Hoeffding et al., 2017; Malomane et al., 2018; Davenport et al., 2020). However, accurately pruning LD sites from total SNP dataset needs high quality whole genomic reference (at the chromosome level). Because Clematis has huge genome size with the reported mean 1C value of 10.48 pg (Min: 6.90 pg and Max: 15.80 pg³), there are still no whole genome references published to date. For this reason, we used another approach to relieve the influence of LD. We randomly selected one SNP from each SCOG alignment using a python script⁴ to obtain a reduced SNP dataset, and then the above mentioned HyDe analysis was also conducted using this dataset.

A Z-statistical test was conducted and probabilities of genetic contributions of putative parents (γ and 1–γ) were calculated. According to Meng and Kubatko’s (2009) hybrid model, the hybrid taxa is either sister to “P1” with probability (γ) or sister to “P2” with probability (1–γ), The null hypothesis was that when hybridization was absent, the expected value of γ should be 0.

Because chloroplast (cp) genomes inherit maternally in Ranunculaceae (Corriveau and Coleman, 1988), we used genome skimming data to assemble complete cp genomes to determine the maternal parent of C. pinnata. For all the eight collecting sites, one individual of C. pinnata and its putative parental species was chosen for DNA extraction and library construction. Five other Clematis and one Anemoclema species (Zhang et al., 2015; Jiang et al., 2017) were also included in the analysis (Table 2).

We used genomic DNA extraction kits (Tiangen Biotech Co. Ltd., Beijing) to extract total genomic DNAs from silica-dried leaves. Extracted DNA was then sent to Biomarker Technologies (Beijing, China) for library construction and next-generation sequencing (NGS). Paired-end reads of 2 × 150 bp were generated on an Illumina Hiseq 4000 genome analyzer platform and raw reads were filtered using the FASTX-Toolkit⁵ to obtain high-quality data by deleting adaptors and low-quality reads.

We then used the Map to Reference option in Geneious Prime 2019 (Kearse et al., 2012) and reference sequences (MG675223.1 and MG675222) to filter out cp reads. Putative cp reads were used for de novo assembly using Geneious Prime 2019 (Kearse et al., 2012) with a low sensitivity setting to reconstruct the complete cp sequence. Gaps were bridged using 20 replicates in a FineTuning step in Geneious Prime 2019 (Kearse et al., 2012). When a contig containing a large single copy (LSC), a small single copy (SSC), and an inverted repeat (IR) region was assembled, the other IR region was determined and attached to the contig using the Repeat Finder function in Geneious Prime 2019 (Kearse et al., 2012) to construct the complete cp genome sequence. The assembled plastome sequences were annotated using Plastid Genome Annotator (Qu et al., 2019), and then checked manually in Geneious Prime 2019 (Kearse et al., 2012). We used the Organellar Genome DRAW tool to illustrate the newly sequenced cp genomes (Lohse et al., 2013) with accession numbers from MT796599 to MT706622 (Supplementary Table S2).

The complete cp genome dataset was aligned using MAFFT v 6.833 (Katoh et al., 2005). Ambiguous alignments and sites with more than 80% missing data were deleted automatically using a Python script⁶. We used the maximum likelihood (ML) and Bayesian inference (BI) methods for phylogeny reconstruction because these methods are less sensitive to long-branch attraction artifacts than the parsimony method (Bergsten, 2005). The ML trees were generated using RAxML v 8.1.17 (Stamatakis, 2014) under the GTR + G model, with bootstrap percentages computed after 100 replicates.

Bayesian inference was performed with MrBayes v3.2.3 (Ronquist and Huelsenbeck, 2003). The models were tested and set for BI analysis according to our previous study (He et al., 2019). Two Markov chain Monte Carlo (MCMC) chains were run independently. Each consisted of three hot chains and one cold chain for 5,000,000 generations. The trees were sampled every 100 generations. The convergence of the Markov chain was tested by calculating the standard deviation value of split frequencies (less than 0.01). The first 20% of trees were removed as burn-in. The remaining trees were used to construct the consensus tree.

Both hybrids and hybrid species are expected to fall in between parental progenitors in morphology (Zhang et al., 2020). To test hybridization hypothesis of C. pinnata, 16 qualitative and ten quantitative morphological characters were chosen for analysis (Supplementary Tables S3, S4). In total, 278 specimens of the four tested taxa (C. pinnata: 29; C. brevicaudata: 124, C. heracleifolia: 56, and C. tubulosa: 69) from multiple herbaria were measured (information of all the specimens is presented in Supplementary Table S5). The selected qualitative characters are taxonomically important traits based on previous taxonomic study (Wang and Xie, 2007). We measured more than 90 leaflets, stems, and flowers of each species to retrieve the quantitative characters.

For the 10 quantitative characters, Kolmogorov–Smirnov test and F-test were used for normality test. Normal distribution data was used one-way ANOVA (analysis of variance) to analysis, and non-normal distribution data was used Kruskal–Wallis test to analysis via IBM SPSS statistics software v25 (SPSS Inc., Chicago, IL, United States), to test the significant differences among the four tested species. Significant differences between any two taxa were identified using a post hoc Tukey’s honest significant difference (HSD) test with false discovery rate (FDR) correction (Benjamini and Hochberg, 1995). Box charts and principal component analysis (PCA) were applied to visualize the differences between all the taxonomic units by R (R Core Team, 2018).

To determine whether C. pinnata has its own ecological niche independent of its parents, we used species distribution models (SDMs) with 278 distribution records (Supplementary Table S5) (C. pinnata: 29; C. brevicaudata: 124, C. heracleifolia: 56, and C. tubulosa: 69 from their entire distribution ranges) and 33 high-resolution environmental variables to project and predict the current distribution patterns of each taxon. The distribution records were obtained from our field investigation and the Chinese Virtual Herbarium database (CVH⁷). For CVH data, the identification of specimen was checked based on the previous taxonomic revision (Wang and Xie, 2007). When no exact GPS information was recorded for a specimen, the geographic coordinates were determined using Google Earth 7.0⁸. In order to avoid sampling bias, only one individual was retained in each 1.0 × 1.0 km square using the “spatially rarify occurrence data” tool in the SDM toolbox for the ArcGis software (Esri, Redlands, CA, United States).

Because many ecological factors, such as bioclimatic, sun light, vegetation, terrain, and soil, may influence the distribution of the four tested species, we obtained 33 environmental variables with 2.5′ spatial resolution (Supplementary Table S6) to conduct niche modeling. Those variables consist of 19 bioclimatic variables from WorldClim-Global Climate Data⁹ (Hijmans et al., 2005); six UV-B radiation variables from a global UV-B radiation dataset for macroecological studies¹⁰ (Beckmann et al., 2014); five vegetation and terrain variables from Harmonized World Soil Database v 1.2 of the Food and Agriculture Organization¹¹ (Food and Agriculture Organization [FAO] and International Institute for Applied Systems Analysis [IIASA], 2012), and three soil variables from the University of Wisconsin¹² (New et al., 1999).

Because collinearity between environmental variables can lead to wrong modeling results (Dormann et al., 2013), we applied the Pearson correlation coefficients using the “banding collection statistics” tool in ArcGis 10.2 and the Jackknife test implemented in MaxEnt to evaluate the contributions of each variable. All the environment variables were converted to ASCII format using an ArcGIS 10.2 conversion tool. Next, we obtained a map of Asia from the global administrative region database¹³. After removing overly correlated variables (r < 0.7) (Guisan and Thuiller, 2005), we performed predictive analysis for each species using MaxEnt 3.4.1¹⁴. We ran 10 bootstrap replicates, in which 25% was used for model testing and the other 75% of the presence data was randomly selected for model training to optimize the model. The model was parameterized with a maximum of 10,000 background points, a convergence threshold of 0.00001, and a maximum of 500 interactions (Lin et al., 2020).

We used area under the curve (AUC) values to evaluate model prediction accuracies. An AUC value is the area enclosed by the receiver operating characteristic curve and the abscissa, and the closer the value is to 1, the more predictively accurate the model is (Fielding and Bell, 1997; Babar et al., 2012). We used the jackknife method to evaluate the weight of each environmental factor on the distribution area, and based on these values, we classified potentially suitable habitats into four categories (which was widely applied by other niche modeling studies, such as Tang et al., 2017; Shitara et al., 2018; Lin et al., 2020, and many others) in a final predictive map of species: unsuitable (<0.2), barely suitable (0.2–0.4), moderately suitable (0.4–0.6), highly suitable (>0.6).

Then, the ecological niche overlap between tested species was checked following the methods by White et al. (2018). We tested niche overlap of species pairs using the ENMtools 1.4.4 (Warren et al., 2010) by calculating Schoener’s D (Schoener, 1968) and Warren’s I statistic values, which 0 indicates no overlap and 1 indicates full overlap (Warren et al., 2008). We then used a niche equivalence test species pairs by comparing statistics D and I using 100 pseudo-replicates (Warren et al., 2008) to test whether the species pairs have identical ecological niche modeling under the null hypothesis. When the statistic value is smaller than the pseudo-replicates value, the result indicates that the two species did not occupy the same ecological niche.

Next, to assess the extent of niche overlap between species pairs, we applied environmental PCA (PCA-env) method (Broennimann et al., 2012) as implemented in the R package of ecospat (Di Cola et al., 2017; R Core Team, 2018)¹⁵. We transformed the multidimensional space of environment variables into two-dimensional space by means of PCA. Following White et al. (2018), we set the resolution to 100 with each grid corresponding to a unique environmental space as suggested by Broennimann et al. (2012), and then we calculated the smoothness of species occurrence by using kernel density to project species onto a grid. We used Schoener’s D (Schoener, 1968) and Warren’s I statistics (Warren et al., 2008) to calculate niche overlap and we accepted the null hypothesis (two species’ niches were equivalent) when p < 0.05.

Our FCM results showed that the ratio of the mean (G₀/G₁) of C. × pinnata, C. tubulosa, and C. heracleifolia to that of the reference (C. brevicaudata) ranged between 0.95 and 1.22. These results showed that all tested species were of the same ploidy level, i.e., 2n = 16 (Supplementary Tables S7, S8 and Supplementary Figure S2).

We generated 30 transcriptome datasets from the four tested species and other outgroup species in this study. Detailed information of the transcriptome data was presented in Supplementary Table S1. Because the sampled species are closely related, we filtered out 3198 SCOGs without missing data, using all-against-all alignment strategy. Among them, 28 SCOGs were tested to generate bad trees with suspicious branch length by TreeShrink. So, we kept 3170 SCOGs for phylogenetic reconstruction. Our results showed that the phylogenetic backbone without putative hybrid taxa (C. × pinnata and C. ochotensis, Supplementary Figure S3) is fully consistent with the results from phylogenetic analysis with inclusion of the hybrid taxa (Figure 2). Clematis × pinnata, C. brevicaudata, C. heracleifolia, and C. tubulosa formed a well-supported clade by both concatenated (ML BS = 100) and coalescent methods (Local posterior probabilities, LPP = 1, Figure 2). However, the resolution and support values within this clade are low, and samples of C. × pinnata were not tested to be monophyletic within the clade. The network analysis (Figure 3) showed that C. brevicaudata and C. heracleifolia/C. tubulosa are well separated. C. heracleifolia and C. tubulosa are closely related but formed two strains, respectively. Samples of C. × pinnata did not form a single strain in this analysis.

We selected 3000 SCOGs with aligned length longer than 500 bp for PhyloNet analysis. All the eight tested populations showed a gently increasing probability value when the maximum number of reticulations was set to two (Supplementary Figure S4). For this reason, we displayed and discussed the results based on the setting of two maximum hybridization events (Figure 4). The results showed that, except for the dubious sample in population BHS, samples of C. × pinnata from the other seven populations were tested to be hybrids between C. brevicaudata and C. heracleifolia/C. tubulosa. Another possible hybridization event occurred in C. ochotensis which belongs to C. sect. Atragene (Figure 4). For population BHS, the dubious sample was not tested to be a hybrid between C. brevicaudata and C. tubulosa but more closely related to C. brevicaudata in the PhyloNet analyses either with (Figure 4) or without C. ochotensis (Supplementary Figure S5).

In this study, we obtained 316,066 and 3170 SNPs for the complete and reduced SNP datasets, respectively, for HyDe analysis. In general, the two datasets generated similar results (Table 3 and Supplementary Table S9), and hybrid origin of C. × pinnata was confirmed by this analysis. The results from the reduced SNP data (as well as the complete SNP data) showed that the genetic contributions of the putative parents were nearly equal in six of the eight populations (Table 3). In WGY, the parental genetic contribution of C. brevicaudata was 0.33, while genetic contribution (1-γ) of C. heracleifolia was 0.67. For the dubious sample from population BHS, HyDe analysis for the complete SNP data detected weak possible hybrid origin of the dubious plant between C. brevicaudata and C. tubulosa (positive z-score = 3.74) with very low genetic contribution of C. tubulosa (1-γ = 0.01, Supplementary Table S9). However, analysis of the reduced SNP data did not support hybrid origin of the dubious plant (Table 3) which is consistent with the PhyloNet analysis.

We obtained about 4 Gb of clean NGS data from each sample and filtered 132,798–443,873 cp reads for C. × pinnata samples and its close relatives by de novo assembly. Chloroplast genome sizes of C. × pinnata, C. brevicaudata, C. heracleifolia, and C. tubulosa ranged from 159,597 bp (C. tubulosa from YMS) to 159,667 bp (C. brevicaudata from YMS) and the overall GC content of all the four species were around 38%. All the acquired plastome sequences consisted of a pair of IRs (31,038–31,087 bp), separated by an LSC (79,392–79,419 bp) and an SSC (18,093–18,187 bp) regions (Supplementary Figure S6). The cp genomes of all the four Clematis species encoded an identical set of 112 genes, including 18 genes with introns, 78 protein-coding, 29 tRNA and four rRNA genes, and 25 genes are in IR region. Structural variation of the newly sampled cp genomes, such as gene inversion/translocation and IR expansion, was similar to that previously reported for other Clematis species (Liu et al., 2018).

Chloroplast phylogenomic analysis showed that all the samples of C. × pinnata and its possible progenitors formed a strongly supported clade (ML BS = 100, PP = 1, Figure 5). Within this clade, the branch lengths of subclades and terminal branches are very short, and the resolution as well as support values within the clade are relatively low. Individuals of C. × pinnata from different sites were separated from each other. Some C. × pinnata samples grouped with C. brevicaudata and clustered into a subclade (ML BS = 99, PP = 1), whereas others grouped with C. tubulosa and C. heracleifolia and were paraphyletic to the C. brevicaudata-C. × pinnata subclade. Unlike the nuclear phylogeny, the cp genomic phylogeny showed that both C. heracleifolia and C. tubulosa are not monophyletic.

Qualitative characters suggested that C. × pinnata exhibited uniparental character states (e.g., hairs on both sides of lamina are the same with C. brevicaudata, hair of inside sepal is more similar to C. heracleifolia/C. tubulosa), intermediate character states (e.g., hair on stamens, sepal color, and spreading direction of the sepals), and new character states to its progenitors (e.g., leaf type) (Supplementary Table S3).

Eight of the ten quantitative characters showed that C. × pinnata measures fell in between C. brevicaudata and C. heracleifolia/C. tubulosa (Figures 6A–D,F,H–J). Whereas, the other two characters, pedicel and filament lengths, were larger in C. × pinnata than in its putative progenitors (Figures 6E,G). PCA results showed that C. × pinnata’s quantitative characters are intermediate to either C. brevicaudata and C. heracleifolia or to C. brevicaudata and C. tubulosa (Figures 6K,L).

After removing overly correlated variables (r < 0.7) (Guisan and Thuiller, 2005), we reserved nine ecological variables for niche modeling (Supplementary Table S6). The AUC values of the four tested species ranged between 0.986 (C. brevicaudata) and 0.998 (C. × pinnata), indicating that the models had high levels of accuracy and validity. The predicted suitable habitats of the four species are shown in Supplementary Figure S7. Jackknife results (Supplementary Table S10) revealed that UVB4 and BIO18 were important for all the species, whereas BIO8 and BIO19 were only important for C. brevicaudata; UVB1 was important for C. × pinnata, C. heracleifolia and C. tubulosa; BIO15 was a key variable for C. × pinnata and C. tubulosa; and Soil pH was important for C. heracleifolia. The top four key environmental variables for C. × pinnata were the same as those for C. tubulosa.

Both niche equivalency test and PCA-env results rejected the null hypothesis that ecological niches for C. brevicaudata vs. C. heracleifolia and C. brevicaudata vs. C. tubulosa were overlapping (p < 0.05). Whereas, the hypothesis of niches overlapping of C. × pinnata vs. C. heracleifolia, C. × pinnata vs. C. tubulosa, C. × pinnata vs. C. brevicaudata, and C. heracleifolia vs. C. tubulosa, was accepted by the two statistical analyses (Table 4).

The application of genomic data has become an effective means to validate hybridization events. Many natural hybridization events have been proposed and validated in numerous organisms (Roberts and Roalson, 2018; Cao et al., 2019; Glémin et al., 2019; Wang et al., 2019; Zhang et al., 2019; Yang et al., 2020). In this study, our nuclear and cp genome analyses confirmed hybrid origin of C. × pinnata between C. brevicaudata and C. heracleifolia/C. tubulosa (Figures 2–5). Hybrid origin of C. × pinnata was further validated by PhyloNet and HyDe analyses, which have been designed for inferring hybridization (Blischak et al., 2018; Wen et al., 2018). Using populational transcriptome data, we not only confirmed hybrid origin of C. × pinnata but also obtained its parental genetic contributions in all the sampling sites except BHS (Table 3). From HyDe analysis, the genetic contributions of the progenitors of most populations were nearly equal, thus strongly suggesting that most plants of C. × pinnata represent the F1 of the parents, because segregation may occur in F2 or later generations causing significant biased parental contributions (Soltis and Soltis, 2009). Meanwhile, for the dubious sample from BHS, PhyloNet analysis failed to detect hybridization of this sample and that it clustered closely with C. brevicaudata. This result was further confirmed by HyDe analysis using the reduced SNP data (Table 3). The sample from BHS may represent a morphologically variable individual of C. brevicaudata

There are some methodological limitations in this study that should be concerned. The complete SNP data may generate biased results due to linkage disequilibrium (LD). In this study, our strategy (randomly select one SNP from each SCOG) greatly relieved the influence of LD. However, because there are no high-quality whole genome references for Clematis, and some of the 3170 SCOGs may be still physically closely related, it is possible that the effects of LD still exist. Furthermore, separating hybridization signal from other sources of incongruence, such as incomplete lineage sorting (ILS) and ancestral polymorphism in population, is difficult (Maddison, 1997). The ABBA-BABA tests implemented in HyDe analysis (or D-statistic test, Patterson et al., 2012) may fail to distinguish signal between hybridization and ancient population structure using SNP data, because this test do not test for hybridization per se, but for deviations from the Wright–Fisher model of random mating which can occur as a result from many different evolutionary processes (Huynh et al., 2019).

Another issue is concerned with our RNA-seq sampling. In the eight collecting sites, we collected multiple C. × pinnata individuals from five sites (Table 2). These multiple individuals from each site were mixed together for RNA-seq to represent genetic features in each population. However, in view that each C. × pinnata individual may have evolved independently, mixed samples may obtain inaccurate or incorrect result, especially when backcrossing occurred in some individuals. Further analysis using transcriptome data from separate individuals should be conducted to review the impact of this mixed sampling strategy.

In recent years, complete cp genomes have been widely used in plant phylogenetic reconstructions (e.g., Moore et al., 2010; Li et al., 2019). Because cp genomes show uniparental inheritance in most angiosperm species, it is also a good molecular marker for inferring hybridization and introgression (Liu et al., 2020). In this study, the cp genome analysis demonstrated that multiple C. × pinnata samples do not form a single lineage (Figures 2, 3, 5). From cp genome phylogeny, we can clearly tell which one acted as the maternal parent of C. × pinnata (Figure 5). The result showed that all the three putative parental species may have contributed maternally to C. × pinnata. This results also depicted that different C. × pinnata individuals evolved from different hybrid events.

From our FCM and phylogenomic results, individuals of C. × pinnata are homoploid hybrids between C. brevicaudata and C. heracleifolia/C. tubulosa. Different hybridization events happened in different distribution sites, and plants of C. × pinnata may be predominantly F1 individuals. This conclusion is also supported by several lines of evidence. Our field investigations found that C. × pinnata’s individuals always occur in the places where both C. brevicaudata and C. heracleifolia/C. tubulosa are present. The flowering time of all four taxa overlap from July to September, thus creating the potential for cross-pollination. For this reason, if C. × pinnata is a species, it cannot easily survive from the parental introgression. In this condition, homoploid hybrid speciation is often accompanied by ecological isolation between daughter species and its parental species (Liu et al., 2014; Kadereit, 2015; White et al., 2018). However, the niche modeling and the niche equivalency test results demonstrated that C. × pinnata has not adapted to a new ecological niche independent of its parents’. From long-term observation, we found that C. × pinnata suffers from reduced fertility of a high proportion of pollen and fruit abortion, which was also reported by Shi (2003). For example, a C. × pinnata individual in Jiufeng, Beijing, noted only a few normal achenes per year, while most of them did not develop at all.

Except for several prerequisites for natural hybridization, e.g., close kinship, overlapping distribution, similar flowering period, shared pollinators, and same chromosome numbers (Ning et al., 2019), habitat disturbance has been often considered as one of the most important factors promoting hybridization (Arnold, 1997; Li et al., 2017; Liao et al., 2021). In the study area of Beijing and Liaoning Province, we found all the C. × pinnata individuals occurred along the mountain roadside. Whereas, C. brevicaudata and C. heracleifolia/C. tubulosa can occupy much larger distribution ranges than C. × pinnata. In autumn, the mountain roadside area will be cleared by forest workers for fire prevention in northern China. This provides opportunity for the parental species to contact each other, and open up habitat for the new hybrids (Li et al., 2017).

Recent studies have shown that F1 hybrids are common in angiosperms and that they can successfully impede gene flow and thus maintain species boundaries between parental species especially in areas where habitat disturbance is high (Liao et al., 2021). Morphologically, C. sect. Clematis and sect. Tubulosae have been considered distantly related in the genus due to their great morphological divergence. They were often placed into different subgenera due to their different floral characters (Tamura, 1995; Grey-Wilson, 2000). Only recently have molecular phylogenetic studies clarified their close relationship (Xie et al., 2011; Lehtonen et al., 2016; Yan et al., 2016; He et al., 2021). Representative species of these two sections, C. brevicaudata and C. heracleifolia/C. tubulosa shared part of their distribution areas in Beijing, Hebei, and Liaoning Provinces, which are exactly the distribution ranges of C. pinnata. The presence of C. × pinnata (F1 hybrid, with reduced fertility) can impede gene flow between parental species, and maintain species boundaries of C. brevicaudata and C. heracleifolia/C. tubulosa in their contact zones.

In Clematis, other species like C. ochotensis (sect. Atragene) also showed hybrid origin from other two different sections (sect. Montana and sect. Fruticella and/or sect. Meclatis, Figure 4). However, we neither have morphological evidence nor sufficient sampling to discuss about this hybridization event. This study demonstrates that interspecific hybridization between two morphologically highly diverged species can occur naturally in Clematis, and hybridization may play an important role in the evolution and diversification of the genus. Taxonomy of Clematis may have suffered from widely hybridization among morphologically diverged species.

This study provided opportunity to investigate how morphological characters of hybrids can be expressed in comparison with its morphologically highly diverged parents. Hybrids are often expected to be morphologically intermediate. However, morphological analyses of natural and artificial hybrids showed that characters of hybrids can be truly intermediate, or identical to those of either paternal or maternal parent, or even new traits (Soltis and Soltis, 2009). In this study, our morphological analysis showed that C. × pinnata exhibited all kinds of morphological outcomes, such as intermediate characters, uniparental characters, and new characters compared to its parental species (Supplementary Table S3 and Figure 6).

Clematis × pinnata shows extensive variation in morphological characters especially in its leaf types and shapes. Leaves of C. × pinnata are predominantly pinnate with five leaflets, but ternate, bi-ternate, or even simple leaves sometimes also occur in different development stages of the plants (Xie et al., 2005). Morphological characters often exhibit higher variability in hybrids than in hybrid species (Zhang et al., 2020). High morphological variation in C. × pinnata also supports the point that C. × pinnata represent early generation of hybrids rather than a hybrid species.

Our morphological analysis did not clearly distinguish C. heracleifolia from C. tubulosa. Except some qualitative characters, e.g., pedicel length, sepal shape, and pollen type (Wang and Xie, 2007), other morphological characters failed to distinguish the two species from one another. Furthermore, our genome analysis and niche equivalency test also did not clearly separate them. These results raise an interesting issue of species delimitation in C. sect. Tubulosae that need to be studied in the future.

All the previous taxonomic studies recognized C. × pinnata as a distinct species based solely on morphology (Fang, 1980; Xie et al., 2005; Wang and Xie, 2007). However, our study has clearly shown that plants of C. × pinnata has not formed a self-evolving lineage, and are generated by recurrent hybridization events between C. brevicaudata and C. heracleifolia/C. tubulosa in their overlapping zones. Clematis × pinnata cannot hold species status from our analysis, and we make a taxonomic treatment as below.