1. Introduction

Recent advances in phylogenomics have shifted our understanding of evolutionary relationships, moving from the traditional bifurcating tree-like phylogeny (e.g., APG IV, 2016; PPG I, 2016; Li et al., 2024) to a more dynamic, web-like framework (Stull et al., 2023; Xu et al., 2023). This new approach offers a more accurate depiction of evolutionary history by incorporating complex reticulate processes such as incomplete lineage sorting (ILS), hybridization, polyploidization, horizontal gene transfer, and introgression, all of which contribute to phylogenetic discordance (Duan et al., 2023; Steenwyk et al., 2023). Advances in bioinformatics tools now allow for more detailed analyses of these processes, revealing how gene flow and evolutionary interactions complicate traditional phylogenetic reconstructions. This shift in perspective has important implications for taxonomy, prompting a reassessment of species and genus boundaries. By integrating reticulate evolutionary processes, the web-like framework offers a more nuanced and accurate understanding of evolutionary relationships, especially in groups with complex evolutionary histories, such as the tribe Campanuleae (Xu et al., 2023). This approach not only refines genus and species delimitations but also enhances our understanding of the evolutionary dynamics that shape biodiversity.

Fischer (1823) established the genus Adenophora Fisch. by distinguishing a group of species from Campanula L. based on their distinctive annular or tubular disk structures located between the stamens and pistils (Hong, 1983; Hong et al., 2011). Currently, Adenophora includes approximately 41–70 species distributed across eastern Asia, with significant diversity in China and neighboring regions such as the Russian Far East, Mongolia, North Korea, Japan, northeastern India, and Nepal (Lammers, 2007a, 2007b). Only one species (Adenophora liliifolia (L.) A. DC.) extends into Europe, reaching Central Asia and the Caucasus (Castroviejo et al., 2010). Despite these clear morphological and geographical distinctions, defining the boundaries of Adenophora and its related genera remains challenging (Tsoong, 1935). For example, Campanula delavayi Franch. was reclassified as Adenophora delavayi (Franch.) D.Y. Hong after the rediscovery of a previously overlooked disk at the base of the style (Hong, 2015), despite being very short (less than 1 mm). However, continuous variation in both vegetative and reproductive traits complicates the identification of consistent characteristics for accurate generic delimitation (Tsoong, 1935; Federov, 1957; Hong, 1983; Lammers, 2007a; Ge and Hong, 2010; Hong et al., 2011). These findings suggest that morphology-based taxonomic studies presented significant limitations in lineages with rapid radiation.

Previous studies have employed various techniques to investigate the phylogenetic relationships within Adenophora, including pollen morphology, isozyme markers, hybridization experiments, and cytological analyses (Qiu and Hong, 1993; Ge and Hong, 1994a, 1994b, 1995, 1998, 2010). However, these methods have been insufficient to fully resolve its intergeneric relationships. Recent molecular-based phylogenetic studies have provided new opportunities to clarify the relationships among Adenophora-related species. Nevertheless, most of these studies have focused on resolving the relationship of the polyphyletic genus Campanula (Mansion et al., 2012; Xu and Hong, 2021; Xu et al., 2023) or even higher taxonomic levels, such as the East Asian Campanulaceae (Yoo et al., 2018) or the entire Campanulaceae family (Eddie et al., 2003; Crowl et al., 2016). Due to limited taxon and marker sampling in Adenophora, these phylogenetic frameworks have not consistently supported the monophyly of the currently recognized Adenophora (Fig. 1). Some studies (Eddie et al., 2003; Crowl et al., 2016; Yoo et al., 2018) have identified Adenophora to be polyphyletic, with the North Korean endemic genus Hanabusaya Nakai, Campanula delavayi, and C. rigescens Pall. ex Schult. (=C. turczaninovii Fed.) nested within it. These members are sisters to C. lehmanniana Bunge and members of C. subg. Rapunculus as defined by Hong (2015). In contrast, other studies have supported the monophyly of Adenophora, positioning it as sister to all the above members (Mansion et al., 2012; Xu and Hong, 2021). For clarity, we refer to Adenophora and its closely related taxa collectively as "Adenophora and its allies" in this study. Using single-copy nuclear (SCN) genes and plastome-level data, Xu et al. (2023) further confirmed a monophyletic group, including "Adenophora and its allies" and Brachycodonia fastigiatus (Dufour ex A. DC.) Fed. In this study, we conducted extensive taxon and genomic sampling and further elucidated the reticulate evolution of Adenophora and its allies for proposing an updated generic classification.

Deep genome skimming (DGS) has significantly advanced phylogenomic studies by enabling the assembly of hundreds to thousands of SCN genes and complete plastomes (Liu et al., 2021). This extensive genomic data enhances the accuracy of phylogenetic inferences by providing numerous genetic markers across the genome, reducing biases from limited gene sampling. The development of bioinformatics pipelines (Liu et al., 2022; Xu et al., 2023; Jin et al., 2024a; Xie et al., 2025) has further streamlined the use of DGS in phylogenomics. These pipelines automate key processes such as data quality control, assembly, annotation, and phylogenetic analysis, making DGS more accessible to researchers with varying levels of bioinformatics expertise. DGS-based data have resolved numerous evolutionary questions across diverse plant lineages. In the Rosaceae family, studies by Liu et al. (2023) and Xue et al. (2024) clarified complex phylogenetic relationships and species diversification. Similarly, research in Poaceae by Hu et al. (2024) uncovered the hybrid origin of Pseudosasa gracilis S.L. Chen & G.Y. Sheng. DGS has also effectively addressed phylogenetic challenges like reticulate evolution and hybridization in other plant families. Ongoing advancements in sequencing technologies and bioinformatics tools are expected to enhance DGS-based phylogenomics further, enabling more detailed investigations into plant evolution. Overall, DGS-based analyses are pivotal for deepening our understanding of plant biodiversity and evolutionary history (Liu et al., 2021).

In this study, we aim to 1) explore the phylogenetic backbone and complex evolutionary relationships within Adenophora and its allies (Campanulaceae) by employing dense taxon sampling and utilizing multiple sources of genomic evidence, including 1506 SCN genes and 77 plastid coding sequences (plastid CDSs). Additionally, we seek to 2) propose an updated generic delimitation integrating morphological and phylogenomic evidence.

2. Materials and methods 2.1. Taxon sampling, DNA extraction, and sequencing

To construct a robust phylogenetic framework for Adenophora and its related species, our PhyloAI team conducted extensive sampling of taxa and genomes, focusing particularly on close relatives such as Hanabusaya, Hyssaria Kolak., and recently transferred A. delavayi from Campanula. We sampled 48 species (out of 72, varying by different taxonomic references) and 13 subspecies from Adenophora (Lammers, 2007a, 2007b; Hong et al., 2011; Hong, 2015). This diverse sampling represents significant phylogenetic and morphological variation within the genus.

We employed various strategies to resolve reticulation and refine taxonomy across hierarchical levels. To establish a robust phylogenetic backbone, we conducted DGS to assemble hundreds to thousands of nuclear and organellar genes. To evaluate the monophyly of widespread species, we included multiple individuals from each species, such as seven individuals from three subspecies of Adenophora capillaris Hemsl., 11 individuals from two subspecies of A. gmelinii (Biehler) Fisch., 14 individuals from two subspecies of A. petiolata Pax & K. Hoffm., ten individuals from four subspecies of A. stricta Miq., and five individuals of A. triphylla (Thunb.) A. DC. These samples were strategically collected across the entire geographic range of each species to ensure comprehensive representation. Additionally, we sampled 24 species from ten genera—Asyneuma Griseb. & Schenk, Brachycodonia Fed. ex Kolak., Campanula, Hesperocodon Eddie & Cupido, Jasione L., Melanocalyx (Fed.) Morin, Poolea Morin, Rotanthella Morin, Triodanis Raf., and Wahlenbergia Schrad. ex Roth—as the outgroup within the tribe Campanuleae. In total, we collected 165 samples to construct the phylogenetic backbone of Adenophora and its relatives.

Whole-genomic DNA was extracted from silica-dried leaves and herbarium specimens using a modified CTAB method (Li et al., 2013), with the DNA extraction conducted at the Institute of Botany, Chinese Academy of Sciences (IBCAS). DNA quality and concentration were evaluated by agarose gel electrophoresis, and samples meeting the quality criteria were sent to Novogene (Beijing, China) for library preparation and Next-Generation Sequencing. DNA libraries were prepared using the NEBNext® Ultra^TM Ⅱ DNA Library Prep Kit, and paired-end (150 bp) reads were generated on the BGISEQ-500 platform (Novogene, Beijing, China). All raw sequencing data are available in the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) under BioProject PRJNA1205918, with SRA accession numbers SRR32029907–SRR32030042. Detailed information on sample collection and sequencing data is provided in Table S1.

All subsequent analyses were conducted using the Ortho2Web pipeline (available at https://github.com/PhyloAI/Ortho2Web), a pipeline recently developed by the PhyloAI team (Xu et al., 2023). Ortho2Web is designed to disentangle the roles of hybridization and allopolyploidization in reticulation, as demonstrated in a case study of the family Campanulaceae. The pipeline offers an 11-step, user-friendly procedure, making it incredibly accessible to novice researchers.

2.2. Developing single-copy nuclear gene reference

MarkerMiner v.1.0 (Chamala et al., 2015) was used to obtain single-copy nuclear genes for subsequent phylogenetic analysis. Transcriptome data of two Campanulaceae species, Platycodon grandiflorus A. DC. (sample code: IHPC) and Lobelia siphilitica L. (sample code: IZLO), were obtained from the OneKP database (https://db.cngb.org/onekp). The SCN library of Arabidopsis Heynh. was selected as the reference database, and MarkerMiner v.1.0 (Chamala et al., 2015) was employed to align and filter the transcriptome data against the nuclear genes to identify SCNs. Default parameters were used, except for setting minTranscriptLen to 600 to obtain more candidate genes of sufficient length and ensure gene quality.

2.3. Reads processing and single-copy nuclear gene assembly

The raw data were processed to remove low-quality sequences and adapters using Trimmomatic v.0.39 (Bolger et al., 2014) with the following parameters: ILLUMINACLIP.fa: 2:30:10:1, LEADING: 5, TRAILING: 15, SLIDINGWINDOW: 4:15, and MINLEN: 40. FastQC v.0.12.1 (available at https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was then used to assess the quality of the processed data and ensure it met the standards for downstream analysis.

Although the hybpiper-nf pipeline (Jackson et al., 2023) was recently developed, its limited capacity to handle large datasets has significantly hindered its broader application in Tree of Life reconstruction. Given our substantial dataset size (up to 9.89 terabytes (TB) of raw data), we opted to use HybPiper v.2.1.6 (Johnson et al., 2016) for assembling SCN genes and plastid CDSs. This choice ensures efficient processing and reliable assembly despite the extensive volume of data.

HybPiper v.2.1.6 (Johnson et al., 2016) was utilized to assemble SCN genes from the cleaned data of each sample, using references designed by MarkerMiner v.1.0 (Chamala et al., 2015). Sequencing reads were sorted with BWA v.0.7.17 (Li and Durbin, 2009) and SAMtools v.1.19 (Li et al., 2009), and the corresponding reads were then mapped to the references. HybPiper v.2.1.6 was subsequently employed for the de novo assembly of each gene. Reads corresponding to each gene were assembled into contigs using SPAdes v.3.15.5 (Bankevich et al., 2012), with the "–cov_cutoff" parameter set to 5 to maximize contig retention.

To extract coding sequences from assembled contigs, including those containing partial introns, we merged multiple contigs aligning to the reference into stitched contigs using Exonerate v.2.4.0 (Slater and Birney, 2005). We then utilized various scripts to retrieve and visualize the assembled data. Assembly statistics for each sample were generated using the "hybpiper stats" script, and SCN sequence coverage was visualized as a heat map with "hybpiper_heatmap". To identify and retrieve potential paralogous genes, we applied the "hybpiper_paralog_retriever" script to the assembled SCN dataset, presenting the results both statistically and visually. Only non-chimeric sequences were retained for subsequent orthology inference.

2.4. Orthology inference for nuclear genes

Given the extensive genome duplication events within the tribe Campanuleae (Xu et al., 2023) and the numerous paralogs identified in our SCN gene assembly, we employed a tree-based orthology inference method initially proposed by Yang and Smith (2014). This method has been subsequently adapted for DGS datasets by Jin et al. (2024b) and then integrated into the Ortho2Web workflow (Xu et al., 2023).

The orthology inference process began with the construction of homologous gene trees. First, sequences for each SCN gene assembled by HybPiper were renamed, and paralogs within the same sample were appropriately marked. Multiple sequence alignments of the renamed SCN genes were performed using MAFFT v.7.520 (Nakamura et al., 2018) with default settings. Phyx v.1.1 (Brown et al., 2017) was then utilized to trim the alignments by removing sequences with more than 10% missing bases. Subsequently, gene trees for each SCN gene were inferred using RAxML v.8.2.13 (Stamatakis, 2014) with the GTRCAT model and 100 bootstrap replicates to assess tree robustness. In each gene tree, we reduced redundancy by retaining only the sequence with the highest number of characters from clades containing multiple sequences from the same sample, ensuring that each clade represented a single sequence per sample. Finally, TreeShrink v.1.3.9 (Mai and Mirarab, 2018) was applied to the gene trees to identify and remove abnormally long branches, resulting in the final set of homologous gene trees for downstream analyses.

The second step involved inferring orthologous genes from the homologous gene trees by further pruning them to ensure one sequence per species. We employed three orthologous gene inference methods proposed by Yang and Smith (2014), each utilizing distinct strategies to generate three datasets: one-to-one (1to1), Monophyletic Outgroup (MO), and Rooted Ingroup (RT). In generating these datasets, Wahlenbergia marginata (Thunb.) A. DC. was designated as the outgroup, and the minimum number of ingroups was set to 35 to maximize the number of orthologs. The MO method identified homologs with monophyletic and non-repeating outgroups from the homologous gene trees. It rooted the gene trees and sequentially pruned them from the roots to the tips, removing extra branches formed by paralogous genes. This process retained only the branch containing the largest non-repeating ingroup, ensuring it met the minimum ingroup size requirement. In contrast, the RT method did not require a monophyletic outgroup. Instead, it split branches containing non-repeating ingroups and removed outgroups, thereby maximizing the number of orthologous genes. To incorporate outgroups back into the RT dataset for downstream analysis, we retained only orthologs present in both the MO and RT datasets and then added the MO outgroups to the corresponding RT orthologs. The 1to1 method applied the most stringent criteria, retaining only orthologous genes that contained strictly non-repeating groups. This method ensured that each orthologous gene set included only one sequence per species without any duplication. Consequently, the 1to1 method produced the fewest orthologs, followed by the MO method, while the RT method yielded the largest number of orthologs.

2.5. Nuclear sequence cleaning and dataset generation

This step-by-step data-cleaning pipeline was initially proposed by Liu et al. (2022). Recently, multiple successive steps have been integrated into a single script (https://github.com/PhyloAI/Ortho2Web/wiki/5-Data-Cleaning) for seamless use within the Ortho2Web workflow. The MO, RT, and 1to1 datasets were subsequently filtered and cleaned using this enhanced pipeline to remove low-quality sequences. This ensures the accuracy and reliability of the data for downstream phylogenetic tree construction.

Orthologs were aligned using MAFFT v.7.520 (Nakamura et al., 2018) with the parameters "–localpair –maxiterate 1000". TrimAl v.1.4 (Capella-Gutiérrez et al., 2009) was then employed to trim the sequences using parameters "-gt 0.8 -st 0.001", which removed sequences with more than 20% gaps and low similarity. The trimmed sequences were concatenated into a supermatrix using AMAS v.1.0 (Borowiec, 2016). Next, Spruceup v.2022.2.4 (Borowiec, 2019) was utilized to remove abnormal bases from the supermatrix by applying the uncorrected p-distance calculation method. A window size of 50 characters with an overlap of 25 characters between sliding windows ensured accuracy. The fraction parameter was set to 1 to align each sample with all others, and the log-normal distribution criterion was used to identify outliers. This process resulted in a supermatrix with an abnormal sequence cutoff of 0.95. The supermatrix was then split into individual gene sequences based on partitions using AMAS v.1.0 (Borowiec, 2016). TrimAl v.1.4 was reapplied with the same trimming parameters to further refine the data. To retain sequences longer than 150 bp, the exclude_short_sequences.py script (Liu et al., 2022) was used. Orthologous gene trees were estimated using RAxML v.8.2.13 (Stamatakis, 2014) with the GTRGAMMA model and 200 bootstrap replicates to assess tree robustness. Finally, TreeShrink v.1.3.9 (Mai and Mirarab, 2018) was applied to remove abnormally long branches from the orthologous gene trees, simultaneously deleting the corresponding sequences. This data processing pipeline ultimately yielded three high-quality ortholog datasets (MO, RT, and 1to1) for subsequent phylogenetic analysis, and all these three nuclear datasets have been deposited in the Dryad [DOI: https://doi.org/10.5061/dryad.9zw3r22s2].

2.6. Plastid protein-coding gene assembly and dataset generation

Due to significant rearrangements in the plastome of Campanulaceae (Li et al., 2020), we focused on assembling plastid CDSs. To retrieve sequencing reads associated with plastome for each sample, we utilized GetOrganelle v.1.7.7.0 (Jin et al., 2020) with default parameters, given its effectiveness and popularity within the plant systematics community. We selected the complete plastome of Adenophora remotiflora (Siebold & Zuccs.) Miq. (Accession no. OP920648) as the reference. Subsequently, Bowtie 2 v.2.5.4 (Langmead and Salzberg, 2012; Langmead et al., 2019) was employed to recruit reads related to the target plastome. The initially recruited reads served as bait in a multiple extension iterations strategy, enabling the retrieval of nearly all plastome-related reads from the clean data. The resulting plastome-related reads were then used for plastid CDSs assembly.

In the subsequent step, we constructed references for plastid CDSs assembly. First, we downloaded six plastomes from GenBank: Adenophora divaricata Franch. & Sav. (Accession no. ncbi-n:KX462129), A. erecta S. Lee (Accession no. ncbi-n:KX462130), A. kayasanensis Kitam. (Accession no. MZ365443), A. racemosa J. Lee & S. Lee (Accession no. MT012303), A. remotiflora (Accession no. OP920648), and A. stricta (Accession no. ncbi-n:KX462131). Seventy-seven protein-coding genes were extracted from these plastomes using the gb2fasta.py script (available at https://github.com/wpwupingwp/rename), serving as references for plastid CDSs assembly. We then employed HybPiper v.2.1.6 (Johnson et al., 2016) to assemble the plastid CDSs for each sample. Gene recovery for each sample was calculated using the "hybpiper stats" script and visualized as a heat map with "hybpiper_recovery_heatmap". Assembled plastid CDSs were retrieved using the "hybpiper_retrieve_sequences" script. Finally, all assembled sequences were manually reviewed and verified for accuracy using Geneious Prime v.2023.2 (Kearse et al., 2012).

Finally, we generated the organellar datasets through multiple cleaning steps. Briefly, we first employed MAFFT v.7.520 (Nakamura et al., 2018) to align the plastid CDSs datasets using the parameters "–localpair –maxiterate 1000". The aligned sequences were subsequently trimmed with TrimAl v.1.4 (Capella-Gutiérrez et al., 2009) using the parameters "-gt 0.8 -st 0.001" to obtain high-quality alignments for each gene. Consequently, we generated the clean plastid CDSs dataset, and this data matrix has been deposited in Dryad [DOI: https://doi.org/10.5061/dryad.9zw3r22s2].

2.7. Concatenation- and coalescent-based phylogenetic inference

We employed both concatenation- and coalescent-based methods for nuclear phylogenetic inference. In contrast, as reviewed by Doyle (2022), plastid CDSs should be treated as single estimates of the underlying species phylogeny. Therefore, we performed only concatenation-based phylogenetic inference for the plastid CDSs dataset.

For the concatenation-based method, AMAS v.1.0 (Borowiec, 2016) was utilized to concatenate the sequences from each dataset into a supermatrix. PartitionFinder2 (Stamatakis, 2006; Lanfear et al., 2017) was then employed to determine the optimal partitioning scheme and nucleotide substitution model for each dataset, selecting linked branch lengths to optimize evolutionary models. All models were specified within the screening range, using the Akaike Information Criterion corrected (AICc) as the selection criterion. The rcluster algorithm was applied to the nuclear gene dataset (Lanfear et al., 2014), while the greedy algorithm was utilized for the plastid CDSs dataset (Lanfear et al., 2012). Based on the best partitioning scheme and evolutionary model identified by PartitionFinder2, Maximum Likelihood (ML) trees were inferred using both IQ-TREE2 v.2.2.2.7 (Minh et al., 2020) and RAxML v.8.2.13 (Stamatakis, 2014) for plastid CDSs ML tree, while IQ-TREE2 v.2.2.2.7 (Minh et al., 2020) was employed for nuclear ML trees. IQ-TREE2 was employed with the parameters "-B 1000 -alrt 1000" to evaluate tree robustness through 1000 rapid bootstrap replicates and approximate likelihood ratio tests (ALRT). RAxML was used to estimate an alternative tree based on the GTRGAMMA model, utilizing 200 bootstrap replicates.

For the coalescent-based method, RAxML v.8.2.13 (Stamatakis, 2014) was employed to infer gene trees for each SCN gene using the GTRGAMMA model with 200 bootstrap replicates. Phyx v.1.1 (Brown et al., 2017) was then utilized to collapse branches in each gene tree that had support values lower than 10. Following this, ASTRAL-Ⅲ (Zhang et al., 2018) was applied with default parameters to construct a species tree for Adenophora and its allies, thereby inferring the shared evolutionary history among the different gene trees.

2.8. Detecting and visualizing nuclear gene trees discordance

Due to the limited sequence length and informative sites for most plastid CDSs, which result in poorly resolved trees with weak support for groupings (Walker et al., 2019), we conducted phylogenomic conflict analyses. Notably, we employed two alternative methods to provide complementary evidence.

Initially, we utilized phyparts v.0.0.1 (Smith et al., 2015) to analyze topological conflicts among the gene trees. All gene trees and the target phylogenetic trees (ML and species trees) were rooted using Newick utilities v.1.6 (Junier and Zdobnov, 2010), with Wahlenbergia marginata designated as the outgroup. To differentiate the effects of missing taxa and uninformative sites at each node in the gene trees, we performed two separate phyparts analyses using distinct parameters. In the first analysis, we conducted a full concordance analysis by setting the parameter "-a" to 1. Nodes with bootstrap support (BS) below 50% were considered uninformative and subsequently collapsed in each gene tree. Each rooted gene tree was reconstructed and compiled into bipartitions, which were then mapped onto the target ML or species tree to calculate topological concordance and the number of conflicting gene trees. However, this approach assumes no missing taxa in the input gene trees. We performed a second phyparts analysis with the parameter "-a" set to 0 and no BS value cutoff for each node to address this limitation. This analysis comprehensively accounts for potential missing taxa in each gene tree. The results from both analyses were combined using an R script and visualized with a Python script, phypartspiecharts_missing_uninformative.py (available at https://bitbucket.org/dfmoralesb/target_enrichment_orthology). Additionally, we calculated the Internode Certainty All (ICA) value for each node to quantify the degree of inconsistency across the gene trees.

We employed Quarter Sampling (QS; Pease et al., 2018) to assess branch support by repeatedly sampling quartets, thereby distinguishing between topological conflicts and insufficient information at weakly supported nodes. For each clade, we conducted 100 replicates, setting the log-likelihood threshold to 2 while maintaining other parameters at their default values. In the target phylogenetic tree, each internal clade was partitioned into four non-overlapping species subsets. QS randomly sampled quartets from these subsets, and the resulting phylogenetic topologies were statistically evaluated based on the sequence data using likelihood scores. Finally, strong conflicts and weak support at each node of the target phylogenetic tree were assessed using three indicator scores: consistency (quartet concordance, QC), difference (quartet differential, QD), and informativeness (quartet informativeness, QI). The results were visualized using the plot_QC_ggtree.R scripts, available at https://github.com/ShuiyinLIU/QS_visualization.

2.9. Incomplete lineage sorting analyses

In this study, we employed two methods to investigate the role of ILS in the evolution of Adenophora and its allies. First, we utilized Coalescent Simulation (CoSi, https://github.com/PhyloAI/Coalescent-Simulation-Analysis) to examine ILS across the phylogeny globally, and this method has been successfully applied in various previous studies (Liu et al., 2022; Jin et al., 2024a). We extracted subtrees containing the exclusive Adenophora species with Wahlenbergia marginata as the outgroup from each gene tree and inferred species trees using ASTRAL-Ⅲ (Zhang et al., 2018). Subsequently, we simulated 10,000 gene trees based on the inferred species tree using Phybase v.2.0 (Liu and Yu, 2010) under the Multispecies Coalescent (MSC) model. DendroPy v.4.5.2 (Sukumaran and Holder, 2010) was then used to calculate the distances between the simulated gene trees and the species tree, as well as between the empirical gene trees and the species tree. Finally, we compared and visualized the differences in distance distributions between the simulated and empirical gene trees.

As an alternative approach, we performed Mutation Calculation based on Coalescent (MuCCo, https://github.com/PhyloAI/Mutation-Calculation-based-on-Coalescent-Model) analysis to evaluate the effect of ILS on the diversification of Adenophora and its allies. This method assesses ILS at specific nodes by calculating the population mutation parameter, theta (θ) (Cai et al., 2021; Nie et al., 2023; Jin et al., 2024a; Xie et al., 2025). We used RAxML v.8.2.13 (Stamatakis, 2014) to convert the branch lengths of the species tree from coalescent units to mutation units based on nucleotide sequences. For each node, theta value was calculated by dividing the branch length in mutation units by the corresponding branch length in coalescent units. A theta value greater than 0.1 indicates a significant influence of ILS.

2.10. Single-nucleotide polymorphisms calling and gene flow analyses

In this study, we downloaded the genome of Adenophora triphylla var. japonica (Regel) H. Hara (=A. triphylla; accession number: GCA_022495995.1) from NCBI, which then served as the reference genome for single-nucleotide polymorphism (SNP) calling. BWA v.0.7.17 (Li and Durbin, 2009) was then used to align the sequencing reads from each sample to the reference genome, mapping the reads to their corresponding positions. The alignment results were converted and sorted using SAMtools v.1.19 (Li et al., 2009). Duplicate reads were identified with the MarkDuplicates tool in Picard v.3.1.1 (https://github.com/broadinstitute/picard), and SAMtools v.1.19 was again employed to index the marked duplicate files. SNPs were detected using HaplotypeCaller program in the Picard v.3.1.1. The CombineGVCFs function of GATK4 v.4.5.0.0 (McKenna et al., 2010) was used to merge haplotypes from multiple samples, followed by the GenotypeGVCFs function to extract genotypes. Preliminary filtering of variants was performed using VariantFiltration function of GATK4 v.4.5.0.0 with the following parameters: QD < 2.0, FS > 60.0, MQ < 40.0, MQRankSum < −12.5, and ReadPosRankSum < −8.0. All SNPs were then extracted using the SelectVariants function of GATK4 v.4.5.0.0. Secondary filtering was conducted with VCFtools v.0.1.16 (Danecek et al., 2011) using the parameters "–max-missing 0.5 –mac 3 –minQ 30" to retain high-quality SNPs suitable for downstream analyses.

To assess gene flow between species, we employed the Dsuite package (Malinsky et al., 2021). The f₄-ratio statistical test was performed on the SNP data using the Dtrios program from Dsuite. The resulting species relationships were subjected to Benjamini-Hochberg correction for multiple testing using the R package (p.adjust). Finally, the f₄-ratio results were visualized using Ruby (https://github.com/mmatschiner/tutorials/tree/master/analysis_of_introgression_with_snp_data).

2.11. Network analyses

To visualize the net-like phylogenetic relationships among Adenophora and its allies, we constructed a global split network using SplitsTree v.4.19.2 (Huson and Bryant, 2006). This analysis was based on the nuclear 1to1 dataset and utilized default settings, including uncorrected P distance, the EqualAngle network construction algorithm, and the NeighborNet method.

To further explore the complex network evolution of Adenophora species and its allies, the PhyloNetworks toolkit (Solís-Lemus et al., 2017) was employed for phylogenetic network inference. Given its substantial computational requirements, thereby reducing computational pressure during runtime. The dataset includes a representative species from each of the Adenophora and its related groups (Ⅰ–Ⅴ), focusing on the relationships among Adenophora species and its allies. Wahlenbergia marginata was designated as the outgroup.

First, we utilized the raxml.pl script from the PhyloNetworks toolkit (Solís-Lemus et al., 2017) to construct gene trees for each nuclear gene using RAxML v.8.2.13 (Stamatakis, 2014). Based on these gene trees, species trees were constructed using ASTRAL-Ⅲ (Zhang et al., 2018). Subsequently, the readTrees2CF module was employed to count quartet concordance factors (CFs) for all nuclear gene trees. We then used Species Networks applying Quartets (SNaQ, Solís-Lemus and Ané, 2016) to infer hybridization events based on CFs and the species tree. Initially, we set the maximum number of hybridization events (h_max) to 0 to create the initial network; this generated phylogenetic network served as the starting point for building subsequent networks (h_max = 1), and so forth. In this study, we set h_max from 0 to 10 and conducted 11 phylogenetic network inferences. We then counted and visualized the pseudo-deviance scores of these 11 networks. The h_max value at which the pseudo-deviance score drops sharply and begins to stabilize was used as the basis for determining the optimal phylogenetic network. The inferred phylogenetic network was visualized using Dendroscope v.3.8.8.0 (Huson et al., 2007).

3. Results 3.1. Next-generation sequencing and single-copy nuclear gene development

In this study, our DGS sequencing generated 165 raw data, with clean paired-end reads (150 bp) ranging from 196,439,166 (27.44 gigabytes (GB)) to 833,475,842 (116.44 GB). The average number of high-quality reads was 428,256,594. Overall, we produced 9.89 TB of raw data. Using the genome size of Adenophora triphylla (2.9 GB; Kang et al., 2024) as a reference, the sequencing depth for the DGS samples ranged from 9 × (A. pulchra Kitam.) to 29 × (A. cordifolia D.Y. Hong), with an average depth of 15 ×. These sequencing depths are sufficient for high-quality nuclear gene assembly (Liu et al., 2021).

Through sequential analyses of the transcriptomes from Lobelia siphilitica and Platycodon grandiflorus, we identified and compiled 1506 SCN genes. These gene sequences serve as references for our subsequent phylogenomic analyses.

3.2. Assembly of plastid and nuclear orthologous gene datasets

The majority of samples successfully recovered all 1506 SCN genes (Fig. S1), with recovery rates ranging from 53% (799 genes) to 100% (1506 genes). Statistical analysis of the assembled genes revealed the presence of paralogous sequences in several genes within Adenophora, with 292 genes identified as paralogs (Fig. S2). To distinguish orthologous genes from paralogs, we applied orthology inference methods based on phylogenetic tree analysis. This approach resulted in three datasets: the 1to1 dataset containing 1243 orthologous genes, the MO dataset with 1499 orthologous genes, and the RT dataset with 1610 orthologous genes. Following a series of data filtering and cleaning steps, the concatenated alignments for these datasets yielded lengths of 1,253,449 bp (1to1), 1,520,430 bp (MO), and 1,647,358 bp (RT), respectively.

In addition to the nuclear dataset, we also assembled plastid CDSs for all samples. Most samples successfully recovered 77 plastid CDSs, with recovery rates ranging from 82% (63 genes) to 100% (77 genes) (Fig. S3). The final concatenated matrix for the plastid CDSs had an alignment length of 69,798 bp. All datasets are available in the Dryad Digital Repository [https://doi.org/10.5061/dryad.9zw3r22s2].

3.3. Nuclear phylogenetic relationship and gene tree discordance of Adenophora and its allies

Phylogenetic relationships of Adenophora were inferred using both concatenation- and coalescent-based methods across three datasets (1to1, MO, and RT), resulting in six distinct phylogenetic trees (Figs. S4–S9). Three species trees shared identical topologies, and the three gene trees also exhibited consistent topologies. However, significant discordances were observed between the gene trees and the species trees (Figs. S4–S9). We focused on the phylogenetic backbone from the 1to1 dataset (Figs. S4 and S5) for further analysis. All six nuclear trees supported the paraphyly of Adenophora (Group Ⅴ), with the monotypic genus Hanabusaya (Group Ⅳ) nested within it (Figs. 2 and S4–S9). In the species tree, Group Ⅳ was supported as a sister to Adenophora delavayi (Group Ⅴ₁; Local Posterior Probabilities (LPP) = 0.58), and together they were sister to the remaining Adenophora species (Group Ⅴ₂; LPP = 1.0; Figs. 2A and S4). Phyparts analysis showed that 16 out of 233 informative gene trees supported the sister relationship between Groups Ⅳ and Ⅴ₁, indicating weak topological concordance among the gene trees (ICA = 0.04, Figs. 2B, S10, and S11). Our QS analysis further indicated weak support for this topology (QS score = 0.21/0.97/0.99; Figs. 2B and S12). In contrast, the concatenation-based tree presented an alternative topology, where Groups Ⅳ and Ⅴ₂ formed a sister relationship, and together they were sister to Group Ⅴ₁. This group had weak support (SH-aLRT/ultrafast bootstrap support (SH-aLRT/UFBoot) = 12.6/80; Figs. 2B and S5), with only one out of 356 informative gene trees supporting it (ICA = −0.01, Figs. 2B, S13, and S14), and substantial conflict was noted in the QS analysis (QS score = 0.48/0.32/0.99; Figs. 2B and S15), indicating the possible presence of a secondary evolution history.

In the nuclear trees, closely related members of Adenophora formed three distinct groups (Groups Ⅰ, Ⅱ, and Ⅲ) rather than a monophyletic group (Figs. 2A and B, and S4–S9). This relationship was also supported by the results from SplitsTree (Figs. 2C and S16). The sister relationship between Adenophora (Groups Ⅴ + Ⅳ) and Group Ⅲ (Campanula rigescens) was consistently supported by all six nuclear trees (SH-aLRT/UFBoot = 100/100, LPP = 1.0; Figs. 2 and S4–S9), and QS analysis provided full support for this relationship (QS score = 1/−/1; Fig. S12). However, this relationship was supported by only 29 out of 267 informative gene trees (ICA = 0.04; Figs. S10 and S11). The monophyly of Group Ⅱ (Hyssaria) received extremely high support (SH-aLRT/UFBoot = 100/100, LPP = 1.0; Figs. 2 and S4–S9), and was supported by 369 out of 518 informative gene trees (ICA = 0.56; Figs. S10 and S11). Similarly, the monophyly of Group Ⅰ (C. aristata Wall., C. immodesta Lammers, C. calcicole W.W. Sm., C. crenulate Franch., and C. chrysosplenifolia Franch.) was robustly supported (SH-aLRT/UFBoot = 100/100, LPP = 1.0; Figs. 2 and S4–S9), with 415 out of 605 informative gene trees supporting this topology (ICA = 0.49; Figs. S10 and S11). QS analysis also provided the highest support for these three closely related groups (QS score = 1/−/1; Fig. S12).

3.4. Plastid phylogenetic relationship of Adenophora and its allies

Our plastid CDS trees (based on both RAxML and IQ-TREE2) consistently supported the inclusion of Hanabusaya (Group Ⅳ) within the currently circumscribed Adenophora (Group Ⅴ; Figs. 3, S17, and S18). This close relationship was further corroborated by 12 out of 26 informative gene trees in phyparts analysis (ICA = 0.14) and received full support from QS analysis (QS score = 1/−/1; Figs. 3 and S19–S21). Group Ⅳ was sister to Group Ⅴ₂, and together they formed a sister group to Group Ⅴ₁ (A. delavayi; Figs. 3A, S17, and S18). However, the sister relationship between Group Ⅳ and Group Ⅴ₂ was only weakly supported (SH-aLRT/UFBoot = 74.3/98, BS = 76), likely due to limited support from informative genes (five out of 28, ICA = 0.05; Figs. S19 and S20) in the phyparts analysis and the lower QS scores (0.38/0.63/0.84; Fig. S21) from QS analysis. The phylogenetic position of Groups Ⅰ, Ⅱ, and Ⅲ was consistent across the plastid and nuclear trees, with all groups receiving strong support (SH-aLRT/UFBoot > 95/100, BS > 95; Figs. S17 and S18) from the majority of gene trees in phyparts analysis (Figs. S19 and S20) and full support from QS analysis (QS score = 1/−/1; Fig. S21).

3.5. Negligible role of incomplete lineage sorting in explaining gene tree discordance

Each node was assessed using the population mutation parameter (theta, θ) to evaluate the influence of ILS. The node corresponding to the clade containing Adenophora delavayi (Group Ⅴ₁) and Hanabusaya asiatica (Nakai) Nakai (Group Ⅳ) showed a low θ value (θ = 0.0483; Figs. 2B and S22). Similarly, most nodes within the genus Adenophora exhibited low θ values (θ < 0.1; Fig. S22), suggesting the negligible role of ILS for explaining gene tree discordance throughout the phylogeny. Additionally, a global evaluation across the phylogeny of Adenophora and its allies showed that the distribution distances between the practical and simulated gene trees did not overlap, indicating that ILS alone could not account for the observed phylogenetic conflicts (Fig. 2D).

3.6. Rampant gene flow as a driver of reticulate evolution

In the phylogenetic network constructed from the Adenophora dataset and its closely related clades (14 taxa), the pseudo-deviance score gradually stabilized as the h_max value exceeded five, indicating the optimal network. Three strongly supported hybridization events were identified. One event is associated with the diversification of five major groups, specifically the hybrid origin of the most recent common ancestor (MRCA) of the combined group (Ⅲ + Ⅳ + Ⅴ). The other two events occurred within Groups Ⅰ and Ⅴ, respectively. Comparing cytonuclear topologies, we hypothesized that Group Ⅱ served as the maternal parent in the hybrid origin of the MRCA of Group Ⅲ + Ⅳ + Ⅴ, with a genetic contribution value (γ) of 0.609 (Figs. 2E and S23). Additionally, the ABBA-BABA test revealed significant gene flow between Group Ⅱ and species of Group Ⅲ + Ⅳ + Ⅴ (Fig. 4), consistent with the hybridization events detected in the SNaQ analysis. Notably, frequent gene flow within Adenophora suggests that gene introgression may explain the substantial conflicts among gene trees (Fig. 4), further highlighting the complex evolutionary history of Adenophora.

4. Discussion 4.1. Taxonomic reassessment: splitting the polyphyletic genus Campanula

Over the past two decades, extensive studies from plant systematists have consistently confirmed that Campanula is polyphyletic (Eddie et al., 2003; Mansion et al., 2012; Crowl et al., 2016; Jones et al., 2017; Yoo et al., 2018; Xu and Hong, 2021). Various genera, such as Adenophora, Asyneuma, Homocodon D.Y. Hong, Peracarpa Hook. f. & Thomson, and Triodanis, are nested within this polyphyletic group (Hong et al., 2011), and have long been recognized as separate entities. While merging these genera into a single large genus, encompassing 420 to 500 species (Contandriopoulos, 1984; Lammers, 2007a), may seem a plausible solution to the polyphyly of Campanula, this approach has not been embraced by the Campanulaceae community. The consensus within this community has instead been to accept the polyphyletic Campanula, alongside these nested genera, and to maintain the existing classification system.

However, taxonomy plays a critical role in advancing biodiversity exploration and evolutionary research, and in facilitating effective communication among plant scientists and the broader public (Crane et al., 2017; Wen et al., 2017; Jin et al., 2023). Therefore, reevaluating the classification of Campanula is necessary. To accommodate its polyphyly, the optimal solution is to split the genus, recognizing distinct lineages that reflect their evolutionary histories. In this context, we propose redefining the boundaries of Adenophora and its allies, which fall within subclade F of the subtribe Phyteumatinae (tribe Campanuleae) as outlined by Xu et al. (2023).

Elucidating the reticulation within these lineages is crucial to this reclassification. By doing so, we can establish a solid phylogenetic framework that will support the splitting of Campanula into monophyletic genera, each with a well-defined common ancestor. This approach ensures that taxonomy accurately reflects the evolutionary relationships within the group, providing a more coherent and useful classification system for both scientific and public purposes.

4.2. Best practices for generic delimitation in the context of reticulation

The unresolved phylogenetic relationships between Adenophora and its allies largely have been due to the limited markers used in previous studies, such as Mansion et al. (2012) and Xu and Hong (2021) (Fig. 1). Considering this limitation, we performed extensive taxon and genomic-level samplings, including 1506 pre-screened SCN genes and 77 plastid CDSs assembled from ca. 9.89 TB of raw data (165 taxa). However, the significant gene tree discordance greatly challenged our accurate phylogenetic inference, and our traditional tree-like phylogenetic inference yielded well-supported but conflicting topologies (Fig. 2, Fig. 3, S4–S9, S17, and S18). Our further reticulation analysis revealed the negligible role of ILS (Fig. 2B and D) in the evolution of Adenophora and its allies. Instead, hybridization (Fig. 2E) and introgression (Fig. 4) were found to play significant roles in promoting the early diversification of this lineage.

Although the principle of "monophyly" has traditionally guided taxonomic delimitations, the ladder-like topology observed among the major groups (Ⅰ–Ⅴ) in our study presents a significant challenge to redefining the circumscription of Adenophora (Fig. 2, Fig. 3, Fig. 5). In the following, we present multiple lines of evidence to support a redefined circumscription for Adenophora.

While phylogenetic inferences based on hundreds of SCN genes typically yield a well-supported topology, underlying gene tree discordance complicates the situation. This discordance, like a double-edged sword, can lead to alternative phylogenetic hypotheses depending on the inference method or gene sampling used (Smith et al., 2020; Stull et al., 2023). For example, our coalescent- and concatenation-based inferences produced distinct phylogenetic placements for Hanabusaya (Group Ⅳ) and Adenophora delavayi (Group Ⅴ₁), either as sister groups or as successive sister to the core Adenophora clade (Fig. 2). However, all these relationships were supported by very few genes (Fig. 2B), suggesting that the phylogenetic signal remains weak for these groups. Additionally, our gene flow analysis revealed introgression among the five groups, further complicating the interpretation of their evolutionary relationships (Fig. 2, Fig. 4). The significant gene tree heterogeneity observed in these groups likely contributed to the difficulty in achieving a clear, consistent phylogenetic inference (Duchêne et al., 2018). This gene flow, coupled with the lack of strong phylogenetic resolution, underscores the complexity of accurately defining the boundaries of Adenophora.

Morphologically, Adenophora was characterized by the presence of an annular or tubular disk between the stamens and pistils (Hong, 1983; Hong et al., 2011). Notably, the rediscovery of the annular disk in A. delavayi led one of our authors (De-Yuan Hong) to transfer this species to Adenophora (Hong, 2015). To preserve the monophyly of the redefined Adenophora, we considered two potential approaches. First, maintaining Hanabusaya as a separate genus would necessitate establishing A. delavayi as a new genus. However, our examination of herbarium specimens (PE) by De-Yuan Hong revealed that both Hanabusaya and A. delavayi possess disks between the stamens and pistils, albeit very short (less than 1 mm). This morphological similarity suggests a closer relationship between these taxa than previously recognized. Given these findings, merging Hanabusaya and A. delavayi into Adenophora presents a more cohesive solution. This approach ensures that all members of the revised Adenophora share the synapomorphic trait of a nectar disk, thereby maintaining the monophyly of the redefined Adenophora. Additionally, our investigation of Campanula rigescens (Group Ⅲ) confirmed the absence of a nectar disk in this species, further supporting the distinctiveness of the redefined Adenophora. Integrating our phylogenomic data with morphological evidence, we have redefined the genus Adenophora to include all currently recognized Adenophora species, the monotypic genus Hanabusaya, and A. delavayi (Fig. 2, Fig. 3). Furthermore, we propose that the other groups identified in our analysis (Groups Ⅰ, Ⅱ, and Ⅲ) be recognized as separate genera, which we have named Rosomala D.Y. Hong & B.B. Liu, Hyssaria, and Boreoasia D.Y. Hong & Chao Xu, respectively. This revised circumscription not only resolves previous taxonomic uncertainties but also aligns with both genetic and morphological data, providing a robust framework for the classification of these taxa within the Campanulaceae family.

The genus Adenophora exemplifies the persistent challenges in modern plant systematics, particularly in reconciling morphological and molecular evidence for taxonomic delineation. Establishing a robust infrageneric classification remains problematic, as discordances between phenotypic traits and phylogenetic relationships primarily stem from inadequate phenotypic analysis. Such insufficiency may arise from a limited scope of phenotypic observations, methodological deficiencies, and interference from complicated evolutionary mechanisms. Notably, phenomena such as converted or parallel evolution can decouple phenotypic similarity from true phylogenetic relationships, thereby hindering scientific taxonomic treatment unattainable (Lyu et al., 2015; Cheng et al., 2021; Xue et al., 2022). To address this challenge, future systematic studies of Adenophora must adopt a more integrative framework. Recent advances in deep learning (DL, Pyron et al., 2022; Karbstein et al., 2024) offer a promising avenue for synthesizing multi-modal data—including phenotypic traits, genomic information, ecological factors, and geographic distributions—ultimately enabling more efficient and objective species delimitation.

4.3. Proposal for two new genera and one reinstated genus

Our phylogenetic analysis consistently placed Group Ⅰ as a basal lineage, a result in agreement with previous studies that either identified it as the basal clade (Mansion et al., 2012; Yoo et al., 2018) or generated unresolved, comb-like phylogenies (Crowl et al., 2016; Xu and Hong, 2021) based on limited plastid or nuclear data. Hybridization analysis using SNaQ further supported a hybrid origin for Group Ⅰ (h_max = 3; Fig. S23), with the MRCA of a combined group (Ⅱ + Ⅲ + Ⅳ + Ⅴ) as the maternal parent (γ = 0.201) and the MRCA of Group Ⅱ (Hyssaria) as the paternal parent (γ = 0.799). However, we found no evidence of gene flow between Rosomala (Group Ⅰ) and Hyssaria (Group Ⅱ), which contradicts the hybrid origin hypothesis (h_max = 3; Fig. S23). The divergence between these two clades occurred ca. 15–9 Million years ago (Mya; Mansion et al., 2012; Crowl et al., 2016). We hypothesize that the uplift of the Qinghai-Tibet Plateau during the Miocene (ca. 23–5 Mya), which followed the initial divergence of these groups, acted as a major geographic barrier that disrupted migration routes and substantially reduced gene flow between them (He and Jiang, 2014; Liu et al., 2017). This geological event likely facilitated genetic isolation and allopatric speciation, offering a plausible explanation for the absence of gene flow between Group Ⅰ (Rosomala) and Group Ⅱ (Hyssaria) in our study.

Morphologically, all members of Group Ⅰ share distinctive traits, such as rosulate basal leaves that persist at anthesis, ranging from reniform to cordate or nearly orbicular in shape. The cauline leaves are typically basal, with the upper ones, if present, being linear and sessile or nearly so. Geographically, this clade is restricted to southwest China, particularly the Hengduan Mountain region and the southern edge of the Himalayas. Integrating the phylogenomic, hybridization, morphological, and geographic evidence, we propose recognizing Group Ⅰ as a new genus, Rosomala. This recognition is supported by its unique phylogenetic position and distinct morphological characteristics.

Group Ⅱ (Hyssaria) has a complex evolutionary history and plays a pivotal role in the hybrid origin of several other clades, notably the combined group (Ⅲ + Ⅳ + Ⅴ, i.e., Boreoasia + Hanabusaya + Adenophora; Fig. 2E). Our gene flow analysis (Fig. 4) reveals frequent introgression between Group Ⅱ and all the other members of this combined clade. Morphologically, Group Ⅱ is characterized by its smaller stature and diminutive flowers, setting it apart from the other groups (Ⅲ + Ⅳ + Ⅴ). As a monotypic genus, Hyssaria is geographically confined to the arid regions in Central Asia, further reinforcing the distinctiveness of this clade. Based on the phylogenetic, hybridization, morphological, and geographical evidence, we propose reinstating Hyssaria as a separate genus, distinct from Campanula. This reclassification reflects its unique evolutionary lineage and ecological traits.

Group Ⅲ, the closest related lineage to the redefined Adenophora, can be easily distinguished by the absence of a nectar disk, a key characteristic present in the latter. Based on this distinctive morphological feature, we propose recognizing this lineage as a new genus, Boreoasia.

In summary, these proposed taxonomic revisions—Rosomala for Group Ⅰ, reinstating Hyssaria, and establishing Boreoasia for Group Ⅲ—are supported by comprehensive phylogenetic, hybridization, morphological, and geographic evidence, and they reflect the distinct evolutionary histories and ecological specializations of these groups.

4.3.1. Key to Adenophora and its allied genera

1. Basal leaves rosulate, persistent at anthesis, reniform, cordate to nearly orbicular; cauline leaves mostly at base, upper ones, if present, usually linear, sessile, or nearly so; endemic to Southwest China and the adjacent regions................................................................................................................ Rosomala (Group Ⅰ)

1. Basal leaves rosulate, with the base covered by dry remnants of leaf petioles, lanceolate to nearly linear (or no basal leaves); cauline leaves many and evenly distributed, mostly alternate, rarely opposite or verticillate, upper ones usually linear, sessile....................................................................................................................... 2

2. Stems short, usually less than 10 cm; corolla is small, shorter than 2 cm; endemic to Central Asia............... . Hyssaria (Group Ⅱ)

2. Stems high, usually more than 10 cm.............................................. 3

3. Flowers with a conspicuous tubular (rarely annular) nectar disk; East Asia, extending to Europe through Central Asia.................................................................. Adenophora (Group Ⅳ + Ⅴ)

3. Flowers without a conspicuous nectar disk; endemic to Siberia to Russian Far East and Mongolia............. Boreoasia (Group Ⅲ)

4.3.2. Group Ⅰ: Rosomala D.Y. Hong & B.B. Liu, gen. nov. (Chinese name: 莲座风铃属)

Type: Rosomala aristata (Wall.) D.Y. Hong & B.B. Liu ≡ Campanula aristata Wall.

Diagnosis. Plants in this genus are characterized by basal rosette leaves that persist during flowering, with a shape ranging from reniform to cordate, ovate, or nearly orbicular, and typically long petioles. The cauline leaves exhibit considerable variation along the stem. Most are located at the lower part of the stem and are typically petiolate. However, the upper cauline leaves gradually become linear, with petioles nearly absent, except in Rosomala chrysosplenifolia. The corolla is relatively small, not exceeding 3 cm in length, and the base of the style lacks a nectar disk. This genus is endemic to Southwest China (Hengduan Mountain and southern edges of the Himalayas) and adjacent regions.

Description. Plants perennial, roots carrot-shaped, stems 2 to several, caespitose, sometimes slender, pubescent or glabrous, usually branching or not branching. Basal leaves rosulate, persistent at anthesis, reniform, cordate, ovate to nearly orbicular, with long petioles; cauline leaves usually showing a transitional form (except for Rosomala chrysosplenifolia), lower cauline leaves mostly ovate, lanceolate to broad-linear, often petiolate, upper cauline leaves gradually transitioning to linear, with or without petioles. Flowers usually solitary on main stems and branches, corolla pendulous or erect, calyx lobes filiform or linear, sometimes toothed along the margins, corolla blue or blue-purple, small, 5–30 mm long. Capsule cylindrical or inversely ovate-conical, with pore dehiscence above the middle.

Distribution. Species of this genus are found in Gansu, Qinghai, Shaanxi, Xizang, and Yunnan in China, with distribution extending to Bhutan, Nepal, Pakistan, Kashmir region, and India, growing at elevations of 2300–5000 m in rocky areas, grasslands, and shrublands.

Etymology. The name Rosomala reflects the morphological and geographical features of this newly proposed genus. The prefix "Roso-" highlights the distinctive rosulate basal leaves, a key characteristic that sets this genus apart within the Campanulaceae family. The suffix "-mala" is derived from "mālā," meaning garland, symbolizing its endemic distribution along the southern edges of Himalayan mountain ranges, extending across the Pan-Himalayan region, including Hengduan Mountain, Bhutan, Nepal, northern India, Kashmir, and North Pakistan.

Including five species: Rosomala aristata, R. calcicola, R. chrysosplenifolia, R. crenulata, and R. immodesta.

New combinations

Rosomala aristata (Wall.) D.Y. Hong & B.B. Liu, comb. nov. ≡ Campanula aristata Wall. in Roxburgh, Fl. Ind. 2: 98. 1824.

=Campanula aristata var. longisepala Marquand, J. Linn. Soc., Bot. 48: 196. 1929.

=Wahlenbergia cylindrica Pax & K. Hoffm., Repert. Spec. Nov. Regni Veg. Beih. 12: 501. 1922. ≡ Campanula cylindrica (Pax & K. Hoffm.) Nannf., Acta Horti Gothob. 5: 24. 1930.

Distribution: Afghanistan, China (Gansu, Qinghai, Shaanxi, Sichuan, Xizang, and Yunnan), India, Kashmir, Nepal, and Pakistan.

Rosomala calcicola (W.W. Sm.) D.Y. Hong & B.B. Liu, comb. nov. ≡ Campanula calcicola W.W.Sm., Notes Roy. Bot. Gard. Edinburgh 12: 196. 1920.

Distribution: China (Sichuan and Yunnan).

Rosomala chrysosplenifolia (Franch.) D.Y. Hong & B.B. Liu, comb. nov. ≡ Campanula chrysosplenifolia Franch., J. Bot. (Morot) 9: 364. 1895.

= Campanula leucotricha C.Y. Wu, Yunnan Trop. Subtrop. Fl. Res. Rep. 1: 58. 1965.

Distribution: China (Sichuan and Yunnan).

Rosomala crenulata (Franch.) D.Y. Hong & B.B. Liu, comb. nov. ≡ Campanula crenulata Franch., J. Bot. (Morot) 9: 365. 1895.

= Campanula nephrophylla C.Y. Wu, Yunnan Trop. Subtrop. Fl. Res. Rep. 1: 60. 1965.

Distribution: China (Sichuan and Yunnan).

Rosomala immodesta (Lammers) D.Y. Hong & B.B. Liu, comb. nov. ≡ Campanula immodesta Lammers, Novon 8: 34. 1998.

Distribution: China (Sichuan, Xizang, and Yunnan), India (Sikkim and Darjiling), and Nepal.

4.3.3. Group Ⅱ: Hyssaria Kolak., Soobshch. Akad. Nauk Gruz. SSR 103(1): 151. 1981. (Chinese name:天山风铃属)

Type: Hyssaria lehmanniana (Bunge) Kolak. ≡ Campanula lehmanniana Bunge.

Diagnosis. Plants of this genus are typically small, usually not exceeding 10 cm in height, with the base of the stems often covered by fibrous, dried leaf petioles. The basal and lower cauline leaves are oblong-elliptic to nearly linear, while the upper cauline leaves are linear and sessile. The flowers are upright before blooming and are relatively few. The corolla is small, generally less than 2 cm in length.

Distribution. Kazakhstan, Kyrgyzstan, Tajikistan, and Uzbekistan.

Including only one species and three subspecies, Hyssaria lehmanniana (subsp. lehmanniana, subsp. capusii, and subsp. pseudohissarica).

New combination

Hyssaria lehmanniana (Bunge) Kolak. subsp. pseudohissarica (Kamelin ex Rassulova) D.Y. Hong & B.B. Liu, comb. nov. ≡ Campanula lehmanniana subsp. pseudohissarica Kamelin ex Rassulova, Fl. Tadzhiksk. S.S.R. 9: 525, 156. 1988.

Distribution: Tadzhikistan.

4.3.4. Group Ⅲ: Boreoasia D.Y. Hong & Chao Xu, gen. nov. (Chinese name:北亚风铃属)

Type: Boreoasia rigescens (Pall. ex Roem. & Schult.) D.Y. Hong & Chao Xu (≡ Campanula rigescens Pall. ex Roem. & Schult.).

Diagnosis. Plants of this genus are typically tall, often exceeding 40 cm in height, with the base of the stem covered by fibrous, dried petioles. The basal and lower cauline leaves are lanceolate to nearly linear, while the upper cauline leaves are linear and sessile. Flowers are either solitary or arranged in a racemose inflorescence. The calyx lacks appendages, with linear or linear-lanceolate teeth. The corolla is relatively large, typically exceeding 3 cm in length.

Description. Perennial plant; root fusiform, slightly branching; stem 10–40 cm or taller, erect, generally few-branched, smooth or sparsely pubescent, stem covered at base with fibrous remnants of leaves. Basal and lower cauline leaves oblong-elliptic, lanceolate, sometimes sublinear, usually entire or denticulate, with petioles nearly equal to the lamina; upper cauline leaves lanceolate-linear, sessile, with an entire margin. Flowers solitary or 2–5 in racemose inflorescences; calyx without appendages, teeth linear or linear-lanceolate; corolla large, 3–3.5 cm long, blue, narrowly infundibular-campanulate, with slightly ciliate margins. Capsule opening by small holes at apex or sometimes about middle.

Distribution. This species is native to regions extending from Siberia to the Russian Far East and Mongolia, predominantly thriving in cold climates.

Etymology. The name Boreoasia reflects both the morphological and geographical essence of the genus. The prefix "Boreo-" highlights its adaptation to cold climates and northern distribution, while the suffix "-asia" references its geographical range across north Asia. Together, the name emphasizes its ecological specialization in cold regions of the Asian continent.

Including only one species, Boreoasia rigescens (Pall. ex Roem. & Schult.) D.Y. Hong & Chao Xu.

New combination

Boreoasia rigescens (Pall. ex Roem. & Schult.) D.Y. Hong & Chao Xu, comb. nov.

≡ Campanula rigescens Pall. ex Roem. & Schult., Syst. Veg., ed. 15 bis [Roemer & Schultes] 5: 102. 1819.

= Campanula infundibulum Rchb., Iconogr. Bot. Pl. Crit. 1: t. 75, f. 158. 1823, nom. illeg., non Vest (1819). ≡ Campanula silenifolia Fisch. ex A. DC., Monogr. Campan. 320. 1830, nom. illeg., non Host (1827). ≡ Campanula stevenii var. silenifolia (Fisch. ex A. DC.) Regel, Bull. Soc. Imp. Naturalistes Moscou 40(3): 187. 1868. ≡ Campanula simplex var. silenifolia (Regel) Trautv., Trudy Imp. S.-Peterburgsk. Bot. Sada 5: 540. 1877.

= Campanula baicalensis Pall., Prodr. 7(2): 479. 1839, pro syn.

= Campanula ciliata Patrin, Prodr. 7: 479. 1839, pro syn.

= Campanula stevenii var. dasycarpa Regel, Bull. Soc. Imp. Naturalistes Moscou 40(3): 188. 1868. ≡ Campanula simplex var. dasycarpa (Regel) Trautv., Trudy Imp. S.-Peterburgsk. Bot. Sada 6: 85. 1879.

= Campanula stevenii var. integerrima Regel, Bull. Soc. Imp. Naturalistes Moscou 40(3): 188. 1868.

= Campanula turczaninovii Fed., Komarov. Fl. URSS. 24: 304. 1957. ≡ Neocodon turczaninovii (Fed.) Kolak. & Serdyuk., Zametki Sist. Geogr. Rast. 40: 29. 1984. ≡ Campanula stevenii subsp. turczaninovii (Fed.) Victorov, Novosti Sist. Vyssh. Rast. 34: 230. 2002.

Distribution: Siberia (Altai, Buryatia, Chita, Irkutsk, Krasnoyarsk, Tuva, Yakutia); Russian Far East (Khabarovsk, Magadan); Mongolia.

4.3.5. Group Ⅳ + Ⅴ: Adenophora Fisch., Mém. Soc. Imp. Naturalistes Moscow 6: 165. 1823. (Chinese name:沙参属)

≡ Campanula sect. Adenophora (Fisch.) D. Dietr., Syn. Pl. 1: 755. 1839. ≡ Campanula subg. Adenophora (Fisch.) Borbas, Magyar Bot. Lapok 3: 190. 1904. Type: Adenophora verticillata Fisch. (= Adenophora triphylla (Thunb.) A. DC.)

= Hanabusaya Nakai, Bot. Mag. (Tokyo) 25: 161. 1911. ≡ Keumkangsania Kim, Fl. Coreana 6: 94. 1976, nom. illeg. superfl. Type: Hanabusaya asiatica (Nakai) Nakai ≡ Symphyandra asiatica Nakai (= Adenophora asiatica (Nakai) D.Y. Hong & Chao Xu).

= Floerkea Spreng., Anleit. ed. 2, 2: 523. 1818. nom. illeg. non Willd. (1801). ≡ Campanula sect. Floerkea Spreng., Syst. 1: 735. 1825. Type: Floerkea marsupiiflora Spreng. (= Adenophora stenanthina (Ledeb.) Kitag.)

Diagnosis. Annular or tubular disk structures located between the stamens and pistils.

Distribution. Mainly distributed in East Asia—namely China, the Korean Peninsula, Japan, and the Russian Far East—with a few species extending into Siberia, Central Asia, and Europe.

Including ca. 72 species: Adenophora amurica C.X. Fu & Y.M. Liu, A. asiatica (Nakai) D.Y. Hong & Chao Xu, A. biformifolia Y.Z. Zhao, A. biloba Y.Z. Zhao, A. borealis D.Y. Hong & Y.Z. Zhao, A. capillaris Hemsl. subsp. capillaris; A. capillaris subsp. dzukoensis A.A. Mao, N.Sarma & D.K. Roy, A. capillaris subsp. leptosepala (Diels) D.Y. Hong, A. capillaris subsp. paniculata (Nannf.) D.Y. Hong & S. Ge, A. changaica Gubanov & Kamelin, A. coelestis Diels, A. contracta (Kitag.) J.Z. Qiu & D.Y. Hong, A. cordifolia D.Y. Hong, A. daqingshanica Y.Z. Zhao & L.Q. Zhao, A. dawuensis D.Y. Hong, A. delavayi (Franch.) D.Y. Hong, A. divaricata Franch. & Sav., A. elata Nannf., A. erecta S. Lee, A. fusifolia Y.N. Lee, A. gmelinii (Biehler) Fisch. subsp. gmilinii, A. gmelinii subsp. hailinensis J.Z. Qiu & D.Y. Hong, A. gmelinii subsp. nystroemii J.Z. Qiu & D.Y. Hong, A. golubinzevaeana Reverd., A. grandiflora Nakai, A. hatsushimae Kitam., A. himalayana Feer subsp. himalayana, A. himalayana subsp. alpina (Nannf.) D.Y. Hong, A. hubeiensis D.Y. Hong, A. jacutica Fed., A. jasionifolia Franch., A. khasiana (Hook. f. & Thomson) Collett & Hemsl., A. koreana Kitam., A. lamarckii Fisch., A. liliifolia (L.) A. DC., A. liliifolioides Pax & K. Hoffm., A. linearifolia D.Y. Hong, A. lobophylla D.Y. Hong, A. longipedicellata D.Y. Hong, A. maximowicziana Makino, A. micrantha D.Y. Hong, A. morrisonensis Hayata subsp. morrisonensis, A. morrisonensis subsp. uehatae (Yamam.) Lammers, A. nikoensis Franch. & Sav., A. ningxianica D.Y. Hong ex S. Ge & D.Y. Hong, A. palustris Kom., A. pereskiifolia (Fisch. ex Schult.) G. Don, A. petiolata Pax & K. Hoffm. subsp. petiolata, A. petiolata subsp. huadungensis (D.Y. Hong) D.Y. Hong & S. Ge, A. petiolata subsp. hunanensis (Nannf.) D.Y. Hong & S. Ge, A. pinifolia Kitag., A. polyantha Nakai subsp. polyantha, A. polyantha subsp. scabricalyx (Kitag.) J.Z. Qiu & D.Y. Hong, A. potaninii Korsh. subsp. potaninii, A. potaninii subsp. wawreana (Zahlbr.) S. Ge & D.Y. Hong, A. probatovae A.E. Kozhevn., A. pulchra Kitam., A. remotidens Hemsl., A. remotiflora (Siebold & Zuccs.) Miq., A. rupestris Reverd., A. rupincola Hemsl., A. sajanensis Stepanov, A. sinensis A. DC., A. stenanthina (Ledeb.) Kitag. subsp. stenanthina, and A. stenanthina subsp. sylvatica D.Y. Hong, A. stenophylla Hemsl., A. stricta Miq. subsp. stricta, A. stricta subsp. aurita (Franch.) D.Y. Hong & S. Ge, A. stricta subsp. confusa (Nannf.) D.Y. Hong, A. stricta subsp. sessilifolia D.Y. Hong, A. subjenisseensis (Kurbatsky) A.V. Grebenjuk, A. sublata Kom., A. taiwaniana S.S. Ying, A. takedai Makino, A. tashiroi (Makino & Nakai) Makino & Nakai, A. taurica (Sukaczev) Juz., A. teramotoi Hurus., A. trachelioides Maxim. subsp. trachelioides, A. trachelioides subsp. giangsuensis D.Y. Hong, A. tricuspidata (Fisch. ex Schult.) A. DC., A. triphylla (Thunb.) A. DC., A. tuvinica Knjaz., A. uryuensis Miyabe & Tatew., A. wilsonii Nannf., A. wulingshanica D.Y. Hong, A. xiaoxiensis D.G. Zhang, D. Xie & X.Y. Yi, A. xifengensis (P.F. Tu & Y.S. Zhou) P.F. Tu & Y.S. Zhou, and A. ×izuensis H. Ohba & S. Watan.

New combination

Adenophora asiatica (Nakai) D.Y. Hong & Chao Xu, comb. nov. ≡ Symphyandra asiatica Nakai, Bot. Mag. (Tokyo) 23: 188. 1909. ≡ Hanabusaya asiatica (Nakai) Nakai, J. Coll. Sci. Imp. Univ. Tokyo 31: 62. 1911. ≡ Keumkangsania asiatica (Nakai) Kim, Fl. Coreana 6: 94. 1976.

= Hanabusaya latisepala Nakai, Bot. Mag. (Tokyo) 35: 147. 1921. ≡ Keumkangsania latisepala (Nakai) Kim, Fl. Coreana 6: 95. 1976. ≡ Hanabusaya asiatica var. latisepala (Nakai) W. Lee, Lineamenta Florae Koreae 1071. 1996.

Distribution: Korean Peninsula.

5. Conclusion

In this study, we conducted the most comprehensive taxon and genomic sampling to date, generating up to 9.89 TB of DGS data and encompassing most of the species in the core Adenophora and all members within its allies. Our reticulation analyses offer a detailed view of the evolutionary dynamics within this lineage. Despite these advancements, future research should focus on further elucidating the ecological and evolutionary mechanisms driving reticulation. In particular, studies should explore how gene flow and introgression shape these lineages, with an emphasis on the role of hybridization events and the influence of geographical and ecological factors. Incorporating additional genomic data, such as pangenomes, could offer a more nuanced understanding of the evolutionary history.

As sequencing technologies and bioinformatics tools continue to progress, enhancing the methodologies for handling large and complex datasets will be crucial. The Ortho2Web pipeline (Xu et al., 2023) used in this study provides a robust framework for investigating reticulation, and refining this tool for broader applications in plant systematics will be key to advancing future research. Our findings also carry important implications for conservation biology and plant genetics. Understanding the role of hybridization and gene flow not only clarifies species boundaries but also supports the conservation of genetically diverse and potentially cryptic lineages. Furthermore, insights into genomic conflict and reticulation can inform the management and sustainable utilization of plant genetic resources. Ultimately, expanding our understanding of reticulate evolution and improving genomic analysis tools will deepen our insight into the complex evolutionary history of Adenophora and its relatives, leading to more accurate and comprehensive taxonomic classifications that reflect the true diversity of the Campanulaceae family.

Acknowledgements

We thank Bing Liu and Dan Xie (IBCAS), as well as Li-Fang Wang (Yanshan University) for their valuable contributions in sample collections. All the phylogenomic analyses have been run on the PhyloAI supercomputer (https://doi.org/10.12282/PhyloAIHPC), managed by Bin–Bin Liu. All the molecular experiments were performed on the Plant DNA and Molecular Identification Platform (PDMIP). This work was supported by the National Natural Science Foundation of China (grant numbers 32270216 and 32000163 to B.B.L.), the Youth Innovation Promotion Association CAS (2023086 to B.B.L.), and the Biological Resources Programme, Chinese Academy of Sciences (CAS-TAX-24-013 to B.B.L.).

CRediT authorship contribution statement

Xiao-Hua Lin: Writing – original draft, Methodology, Formal analysis. Si-Yu Xie: Writing – review & editing, Writing – original draft. Dai-Kun Ma: Writing – review & editing, Writing – original draft. Shuai Liao: Investigation. Bin-Jie Ge: Methodology. Shi-Liang Zhou: Writing – review & editing, Conceptualization. Liang Zhao: Supervision. Chao Xu: Supervision. De-Yuan Hong: Supervision. Bin-Bin Liu: Supervision, Project administration.

Data accessibility statement

All the DGS data are deposited in the NCBI SRA under the BioProject PRJNA1205918, and the detailed information for each sample is referred to Table S1. The predesigned gene references, sequence alignments, phylogenetic tree files, and other data files generated in this study are deposited at Dryad Digital Repository [https://doi.org/10.5061/dryad.9zw3r22s2].

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.05.010.