1. Introduction

Bignoniaceae, predominantly distributed in Neotropical regions, plays a crucial ecological role in Neotropical Forest ecosystems while also occurring widely across tropical forests in Africa, Madagascar, and Southeast Asia. Comprising 79 genera and 901 species, this family exhibits remarkable ecological and morphological diversity, encompassing lianas, woody shrubs, and tree populations, as well as three independently derived herbaceous lineages (e.g., Incarvillea, Argylia, Tourrettia) (Glade-Vargas et al., 2021; Rana et al., 2021; Jiang et al., 2025; Shaheen et al., 2025) adapted to high-elevation environments in the Himalayas and Andes. This broad adaptive range and striking morphological variation reflect the family's complex evolutionary history, particularly in the inference of deep phylogenetic relationships among its major clades, making it an exceptional system for studying early evolutionary processes and diversification patterns (Gentry 1980; Olmstead et al., 2009).

Current understanding of the deep phylogenetic relationships within Bignoniaceae has remained unclear (Gentry 1980, 1992; Fischer et al., 2004). Molecular studies have divided Bignoniaceae into eight clades: Jacarandeae, Tourrettieae, Tecomeae, Oroxyleae, Catalpeae, Bignonieae, the Tabebuia alliance, and the Paleotropical clade (Olmstead et al., 2009). The Jacarandeae is consistently recognized as the most basal lineage within Bignoniaceae, forming a sister group to all other members of Bignoniaceae. The remaining seven clades, collectively referred to as 'Core Bignoniaceae', exhibit poorly resolved phylogenetic relationships among themselves. Moreover, molecular studies have revealed several persistent taxonomic uncertainties, particularly regarding the placement of Delostoma and Argylia, which were once classified within Tecomeae but whose phylogenetic position remains unresolved (Gentry 1980; Olmstead et al., 2009). Notably, significant topological conflicts have emerged between analytical approaches. In the parsimony strict consensus tree based on concatenated sequences of ndhF, trnL-F, and rbcL, Oroxyleae was recovered as sister to Catalpeae. However, this relationship was absent in the maximum likelihood (ML) tree inferred from only those taxa with complete three-gene data (Olmstead et al., 2009). Furthermore, the phylogeny of Olmstead et al. (2009) exhibited exceptionally short branch lengths across the deep nodes of major early-diverging lineages, with the exception of Jacarandeae, alongside low nodal support values. These patterns collectively suggest that Bignoniaceae likely underwent a complex early evolutionary radiation characterized by rapid diversification events.

In lineages that have undergone rapid diversification or extensive gene flow, clarifying deep-branching phylogenetic relationships is often confounded by factors such as ancient hybridization/introgression, and ILS (Feng et al., 2022; Bjornson et al., 2024; Wang et al., 2025; Yan et al., 2025). ILS arises from the stochastic retention of ancestral polymorphisms during rapid, successive speciation events, causing individual gene tree topologies to fail to fully capture the species tree divergence history (Guo et al., 2023; Xu et al., 2024; Xia et al., 2025). In contrast, gene flow typically occurs post-speciation, where hybridization introduces adaptive genetic variation into recipient populations, thereby augmenting genetic diversity (Lin et al., 2025).

For clades that underwent rapid early diversification, the difficulty lies not only in studying the respective impacts of gene flow and ILS on phylogenetic inference but also in the fact that signals from these processes can become intricately intertwined due to subsequent complex evolutionary events, making them challenging to accurately disentangle (Cai et al., 2021; Cui et al., 2025). To quantify these confounding factors, researchers have developed various phylogenomic analytical tools, such as PhyloNet for inferring hybridization/gene flow history (Than et al., 2008), QuIBL specifically designed to quantify introgression (Edelman et al., 2019), and the recently introduced Phytop for phylogenetic network reconstruction (Shang et al., 2025).

In this study, we aimed to reconstruct the deep-branching relationships within Bignoniaceae at the tribe level, investigate the causes of deep phylogenetic conflicts, and clarify the roles of hybridization, and incomplete lineage sorting in shaping the family's evolutionary history. To reconstruct Bignoniaceae phylogenetic relationships, we utilized extensive datasets of single/low-copy nuclear genes and chloroplast genomes from 88 samples covering 59 genera, all eight recognized tribes, and the phylogenetically critical genera Delostoma and Argylia. We then used phylogenomic approaches (i.e., PhyloNet, QuIBL, Phytop) to determine whether conflicts in deep evolutionary relationships within the Bignoniaceae phylogeny are signs of hybridization/gene flow or introgression.

2. Materials and method 2.1. Sample collection

We sampled 88 accessions, comprising 86 ingroup species of Bignoniaceae and two outgroup species. These samples were sourced from a combination of 28 herbarium specimens, 14 garden-collected samples (from XTBG, SCBG, RBGE, and Kew Gardens), three field-collected samples. Data from an additional 43 samples were downloaded from NCBI and ENA, including six samples with genome resequencing data and 37 samples with target sequence capture data (only containing 353 nuclear genes) (Fonseca et al., 2023; Zuntini et al., 2024). Detailed information is provided in Tables S1 and S2. Our sampling covered all eight tribes of Bignoniaceae (Olmstead et al., 2009), as well as two additional genera, Argylia and Delostoma, which were not assigned to any of these tribes in previous studies. Given the close relationship between Paulowniaceae and Bignoniaceae, we selected two Paulowniaceae species as outgroups (Hug et al., 2016).

2.2. DNA isolation and genome resequencing

Total DNA was extracted from silica gel-dried leaf tissues of herbarium and field-collected samples using a modified CTAB method (Li et al., 2013). The extracts were examined using 0.8% agarose gel electrophoresis, and the Gene Genius Bio imaging system was used for imaging. All samples underwent rigorous quality control prior to sequencing, including measurements of DNA concentration and purity, to ensure they met the requirements of sequencing platforms. DNA was sequenced on the T7 platform (Modi et al., 2021). Most samples yielded approximately 25 Gb of 150-bp paired-end reads. We selected ten species with sufficient sequencing depth (~50 Gb) to represent the ten clades of Bignoniaceae (Olmstead et al., 2009). Given the family's estimated genome size of ~1 Gb (https://cvalues.science.kew.org/), this provided approximately 50× coverage, ensuring sufficient depth for reliable identification of orthologous nuclear genes (Table S3).

High-throughput sequencing data were filtered for quality and adapters were trimmed in Trimmomatic v.0.39 (Bolger et al., 2014). We used a sliding window parameter of 4:15 to remove sequences that did not meet quality standards and ultimately discarded clean reads shorter than 36 base pairs or those lacking paired-end information.

2.3. Chloroplast genome assembly and annotation

Complete chloroplast genomes were assembled in GetOrganelle v.1.6.2 (Jin et al., 2020) with varying k-mer lengths (75, 85, 95, 105, 115, and 125). For genomes that could not be successfully assembled using GetOrganelle v.1.6.2, we followed the approach of Dong et al. (2022) and Zhang et al. (2025). Newly assembled chloroplast genomes were annotated in Geneious Prime v.2021.1.1 (Kearse et al., 2012).

2.4. Orthologous nuclear gene identification

To independently validate our phylogenetic results using a universally conserved gene set, we first mined the Angiosperms353 genes from our entire sample set (88 samples). We employed GeneMiner v.2.0 (Zhang et al., 2022; Xie et al., 2024), a tool designed for efficient, reference-based gene extraction from heterogeneous sequencing data. A reference file containing Angiosperms353 gene (Zuntini et al., 2024) sequences from closely related taxa was compiled. GeneMiner v.2.0 was used to directly mine the target loci from the quality-filtered reads of all samples, bypassing de novo assembly for these specific genes. The resulting sequences for each locus were aligned and processed to produce a concatenated dataset. This dataset was named dataset1.

For the ten species with deep sequencing, we assembled draft genomes using SPAdes v.4.2.0 and assessed their completeness with BUSCO v.6 (Tegenfeldt et al., 2025). Completeness ranged from 71.1% to 95.8% (Table S3). Because SPAdes outputs only raw genomic contigs, we predicted gene models using GlimmerHMM v.3.0.4 (Majoros et al., 2004). Based on these exon–intron annotations, we extracted both the predicted coding sequences (CDS) and their corresponding genomic loci using JCVI v.1.5.7 and BEDTools v.2.31.0 (Quinlan and Hall, 2010; Tang et al., 2024). Orthologous nuclear genes were identified using Proteinortho v.6.3.6 (Klemm et al., 2023), which is specifically designed for the comparative analysis of multiple genomes or large-scale transcriptomes. This tool employed a reciprocal best hit (RBH) strategy combined with a sophisticated graph-based algorithm to distinguish orthologs from paralogs across species. To ensure the identification of high-confidence single/low-copy orthologs, we applied stringent post-filtering criteria: only gene groups present as single/low-copy in all sampled taxa and exhibiting clear, full-length alignments without internal stop codons were retained. This rigorous process yielded a preliminary set of 1471 single/low-copy nuclear genes (SCNGs). A custom reference was then constructed from ten representative genomes by aligning coding sequences of single/low-copy orthologs with MAFFT v.7.526 (Katoh et al., 2019), refining alignments with Trimmomatic v.0.39, and generating consensus sequences. This sequence-to-tree workflow closely follows the reed2tree (R2T) framework, which provides an efficient and standardized approach for integrating orthology assignment, alignment optimization, and downstream phylogenetic inference (Dylus et al., 2024). The other samples with re-sequenced data were mapped to consensus sequences using Bowtie v.2.5.4 (Langmead and Salzberg, 2012). The nuclear genes were generated using samtools v.1.22.1 and bcftools v.1.8 (Danecek et al., 2021), followed by filtering to retain only genes present in > 70% of samples with length > 500 bp. This process resulted in a final set of 1275 high-quality SCNGs for phylogenetic analysis, designated as dataset2.

2.5. Phylogenetic tree construction

The chloroplast dataset only contained the protein-coding sequences of the chloroplast genome. These protein-coding sequences exhibit high sequence similarity across species and provide sufficient genetic variation to accurately resolve phylogenetic relationships (Delport et al., 2009). Protein-coding sequences from each chloroplast genome was extracted with Geneious Prime 2021.1.1 and aligned with MAFFT v.7.526. Alignments were manually adjusted to ensure correct reading frames and to exclude potential misalignments. All 72 aligned protein-coding sequences were concatenated into a single dataset. We then used IQ-TREE v.2.3.6 (Minh et al., 2020) under the ML method to infer the phylogenetic tree, employing ModelFinder to automatically select the best-fit evolutionary model.

For nuclear genes, we employed two approaches to construct phylogenetic trees: the concatenation method and the coalescent-based method. For the concatenation approach, we concatenated the aligned sequences of each single/low-copy nuclear gene into a supermatrix and performed ML analysis using IQ-TREE v.2.3.6 with ModelFinder to select the optimal evolutionary model (Nguyen et al., 2015; Kalyaanamoorthy et al., 2017). The species tree was constructed through a multi-step process. Maximum-likelihood gene trees were inferred using IQ-TREE 2.3.6 and the coalescent-based species tree was generated from the complete set of gene trees using ASTRAL-Ⅲ v.5.6.3 (Zhang et al., 2018). All trees were visualized with Chiplot v.2.6.1 (Xie et al., 2023).

2.6. Tree topology consistency analysis

To comprehensively assess phylogenetic discordance, we performed two complementary analyses. First, we conducted a topological comparison between the chloroplast genome tree and the species tree inferred from nuclear datasets using ggtree package (v.3.16.3) (Xu et al., 2022) in R v.4.4.3, focusing on conflicting placements of major clades or specific taxa. Second, to quantify gene tree conflict at each node of the species tree, we employed Phytop (Shang et al., 2025) to generate a conflict network diagram illustrating the distribution and proportion of supporting versus conflicting gene trees. Additionally, we calculated the quartet score (QS) and gene concordance factor (gCF) for each node to evaluate conflicts. QS and gCF values close to 1 and 100, respectively, indicate strong concordance with the bipartitions defined by the species tree. QS values were computed using ASTRAL, while gCF values were calculated using IQ-TREE v.2.3.6. Finally, PhyParts (Smith and Kolpashchikov, 2017) was used to visualize conflicts among gene trees.

2.7. Assessing phylogenetic conflict signals

We used PhyTop to evaluate the relative contributions of ILS and introgression/hybridization (IH) for phylogenetic discordance, facilitating direct comparisons of network topologies and hybridization signals. To further investigate phylogenetic conflict at deep nodes, we optimized our sampling to balance computational efficiency with analytical rigor (Bordewich and Semple, 2007). From the ten clades (the eight recognized tribes plus the Argylia and Delostoma clades), we selected one representative species per clade with high coverage, plus an outgroup (Table S4). This created a curated subset of 11 taxa specifically designed to clarify tribal-level relationships. For this subset, we identified 727 high-quality SCNGs that were present in all taxa and exceeded 1500 bp in length, which constituted dataset3.

Based on dataset3, we statistically evaluated alternative species relationships using the approximately unbiased (AU) test with multiple testing correction to ensure robust inference (Shimodaira, 2002). Given the moderate dataset size (n = 727 gene trees), we performed 10, 000 multiscale bootstrap replicates (-zb 10, 000) in IQ-TREE v.2.3.6 to mitigate potential biases from limited sampling. All candidate trees were evaluated under identical optimal substitution models to ensure likelihood comparability, with Bonferroni correction (α = 0.05) applied to all p-values to strictly control the family-wise error rate in topological comparisons. The six highest-frequency topologies were visualized.

2.8. Reticulate evolution tests

We assessed the relative contributions of gene tree estimation error, ILS, and gene flow to observe reticulation patterns based on dataset3 using the reticulation index framework (Cai et al., 2021). We first computed reticulation indices using the bootstrap-supported reference species tree, where low nodal support values in individual gene trees indicate greater topological conflict. To distinguish between ILS and hybridization, we analyzed all possible combinations of three ingroup taxa (referred to as a "triplet"), along with a designated outgroup. For each triplet, we compared the frequencies of the three possible discordant gene tree topologies. ILS predicts roughly equal frequencies among these topologies, whereas hybridization predicts a significant deviation in one topology (indicating excess allele sharing between two taxa). These patterns were statistically evaluated against null expectations simulated under a pure ILS model using Phybase v.2.0 R package (von Hardenberg and Gonzalez-Voyer, 2025), with significance assessed via a χ² test. We then quantified ILS influence through θ values (RAxML/ASTRAL branch length ratios) and estimated gene tree error. Gene tree error was estimated by simulating 100 GTR-based alignments in SeqGen v.1.3.4 (Rambaut and Grass, 1997), using empirical parameters from our RAxML v.8 (Stamatakis, 2014) analyses. The length of these simulated alignments was set to 1500 bp to approximate the length of the genes in dataset3, ensuring that the simulation conditions reflected the properties of our empirical data. Finally, we performed relative importance analysis using R package relaimpo to partition variance in gene tree discordance among these three factors, a comprehensive evaluation of their evolutionary impacts.

2.9. Phylogenetic network analysis

We employed the QuIBL (https://github.com/miriammiyagi/QuIBL) framework to investigate hybridization signals among deep nodes using dataset3 (Edelman et al., 2019). This Bayesian approach systematically evaluates phylogenetic discordance by comparing gene tree-species tree concordance under competing evolutionary scenarios: one considering only ILS and another incorporating both ILS and IH. The analysis requires two key inputs: (1) a reference species tree and (2) the complete gene tree set. QuIBL calculates Bayesian Information Criterion (BIC) values for both models, where relative support (ΔBIC = BIC_ILS - BIC_ILS + IH) determines the most plausible explanation for observed topological conflicts. Our threshold ΔBIC > 10 provides strong evidence for introgression events, while lower values favor pure ILS explanations. We randomly subsampled 500 trees for QuIBL analysis, repeating this process 100 times. After excluding topologies congruent with the species tree, we quantified genome-wide introgression proportions for the remaining two topologies using the formula:

\mathrm{H}=\text { mixpro} 2 \times(\text { count } / \text { total gene trees })

A higher H value indicate stronger inter-lineage gene flow. This dual approach combines hypothesis testing (BIC comparison) with quantitative estimation (introgression proportion) for comprehensive reticulation assessment (Tables S5–S13).

We implemented PhyloNet v.3.8.0 using maximum pseudo-likelihood (InferNetwork_MPL) (Than et al., 2008; Wen et al., 2018) to reconstruct reticulate relationships. Our analysis incorporated complete gene trees from all selected species as input. The network search protocol systematically evaluated potential hybridization scenarios: (1) an exhaustive search allowing 0–5 reticulation events; (2) five independent replicates per reticulation configuration to ensure identification of global likelihood optima. Convergence was achieved when successive iterations showed negligible improvement in log-likelihood scores (ΔlnL < 0.01), indicating network topology stability. The statistically best-supported phylogenetic network was then visualized and annotated using Dendroscope v.3.7.2 (Huson and Scornavacca, 2012) for clear interpretation of complex evolutionary relationships.

2.10. Time calibration of the phylogeny and trait data collection

We used the MCMCtree method from PAML v.4.10.9 (Yang 2007) to infer the divergence time of Bignoniaceae at the genus level. The species tree was calibrated using two well-documented fossil constraints from the Bignoniaceae family. The crown age was set to 56–66 Million years ago (Ma) according to the seed fossil (Pigg and Wehr, 2002; Lohmann et al., 2012). A Catalpa fossil from the Whitecap Knoll flora (38.4–39.17 Ma) was used to constrain the divergence between Catalpa and Delostoma (Zuntini et al., 2024).

Dating analysis was performed under an independent rates clock model with HKY substitution parameters. We initially generated branch rate priors using approximate likelihood computation (usedata = 3), then performed full dating using these priors (usedata = 2). The MCMC chain ran for 200, 000 generations following a 2000-generation burn-in, with sampling every 10 generations. Despite the moderate chain length, convergence was rigorously assessed using Tracer v.1.7.2 (Rambaut et al., 2018), which confirmed that all parameters had achieved effective sample sizes (ESS) well above the recommended threshold of 200, and trace plots showed good mixing and stationarity, indicating reliable parameter estimates.

We compiled data on four key traits for each genus of Bignoniaceae (i.e., life form, leaf type, fruit dehiscence mode, and seed wing presence) from Bignoniaceae monographs (Gentry, 1980, 1992; Fischer et al., 2004) and various regional floras (Leeuwenberg, 1973; Woodson et al., 1973; Gentry, 1980; Zhang and Santisuk, 1998; Stevens et al., 2001).

3. Results 3.1. Dataset information

This study employed four phylogenomic datasets, each designed to address specific analytical needs and collectively ensuring comprehensive phylogenetic resolution. Dataset1 (Angiosperms353 dataset) comprised 88 samples including 59 genera (74.68% of total 79 genera in Bignoniaceae). Dataset2 contained 51 samples representing 35 genera. This dataset included 1275 genes, ranging from 501 to 9093 bp. Dataset3 contained 11 samples (one per clade), comprising 727 genes longer than 1500 bp. This dataset was mainly used for computationally intensive analyses of deep-node evolutionary conflicts. Additionally, we successfully assembled 51 complete chloroplast genomes, with lengths ranging from 124, 193 bp to 182, 643 bp. Seventy-two protein-coding genes of the chloroplast genome were extracted forming a supermatrix (chloroplast dataset) with a total length of 69, 124 bp for the following analyses.

3.2. Phylogenetic relationships in Bignoniaceae

The phylogenetic tree based on the chloroplast dataset supported the division of Bignoniaceae into 10 clades (Fig. 1 and Fig. S7). Our phylogenetic tree identified the earliest-diverging lineage as Jacarandeae. The Delostoma clade was found to be sister to all other tribes except Jacarandeae. Bignonieae was sister to a clade containing Oroxyleae, the Tabebuia alliance, and the Paleotropical clade. Oroxyleae was sister to the Tabebuia alliance and the Paleotropical clade.

The species tree based on dataset1 was most consistent with the chloroplast tree (Fig. 1, Figs. S1 and S2), except for the placements of the Argylia and Delostoma clades. In the species tree, the Argylia clade was sister to a clade containing the Delostoma clade, Catalpeae, Bignonieae, Oroxyleae, the Tabebuia alliance, and the Paleotropical clade (Fig. 2a). Meanwhile, the Delostoma clade was resolved as sister to a group formed by Catalpeae, Bignonieae, Oroxyleae, the Tabebuia alliance, and the Paleotropical clade.

Both concatenation and coalescent-based analyses based on dataset2 produced congruent topologies that strongly supported the phylogenetic framework of Bignoniaceae (Fig. 2a and Figs. S4-S6). The results divided the family into 10 clades (Jacarandeae, Tourrettieae, Tecomeae, the Argylia clade, Oroxyleae, the Delostoma clade, Catalpeae, Bignonieae, the Tabebuia alliance, and the Paleotropical clade). Of these, Jacarandeae consistently formed the sister group to all other clades (Fig. 2a), matching previous findings (Olmstead et al., 2009). Tourrettieae and Tecomeae were sister clades. The Argylia clade was strongly supported as monophyletic and sister to the remaining six clades (Oroxyleae, the Delostoma clade, Catalpeae, Bignonieae, the Tabebuia alliance, and the Paleotropical clade), while Oroxyleae was sister to the Delostoma clade, Catalpeae, Bignonieae, the Tabebuia alliance, and the Paleotropical clade. The sister relationship between the Delostoma clade and Catalpeae was clearly resolved, addressing previous phylogenetic uncertainties. Bignonieae was sister to Catalpeae and the Delostoma clade, with these three clades collectively forming the sister to the Tabebuia-Paleotropical clade, which in turn were sister to each other.

While the species tree from dataset2 received maximal support (ASTRAL local posterior probability (LPP) = 1) at all deep nodes, the corresponding nodes in the dataset1 topology were less resolved (LPP < 0.9). This notable incongruence between the two primary phylogenetic hypotheses primarily involved the placements of Tourrettieae, the Delostoma clade, and Oroxyleae. In the dataset1 topology, Tourrettieae was not resolved as sister to Tecomeae, the Delostoma clade did not form a clade with Catalpeae, and Oroxyleae was sister to the Tabebuia-Paleotropical clade. In contrast, the dataset2 topology showed a strongly supported sister relationship between Tourrettieae and Tecomeae; a clear sister grouping of the Delostoma clade with Catalpeae; and Oroxyleae as sister to a larger clade comprising the Delostoma clade, Catalpeae, Bignonieae, the Tabebuia alliance, and the Paleotropical clade.

3.3. Evaluating conflicting phylogenetic signals in Bignoniaceae

To systematically investigate gene tree discordance, we performed a quantitative assessment of conflicting phylogenetic signals based on dataset2. QS and gCF analyses visualized early divergence patterns (Fig. 2a), showing low concordance values for major deep nodes (QS < 0.6; gCF < 0.4), including node C (Argylia clade; QS = 0.41, gCF = 0.11), the Oroxyleae-related node (QS = 0.53, gCF = 0.09), and the Delostoma clade-related node (QS = 0.54, gCF = 0.39). Other nodes with relatively low concordance values included B (QS = 0.45, gCF = 0.14), D-1 (QS = 0.39, gCF = 0.04), and D-2 (QS = 0.45, gCF = 0.14). These low QS and gCF concordance values provided quantitative support for complex evolutionary dynamics during early radiation, likely involving hybridization or ILS.

We further assessed relative contributions of ILS and IH to gene tree conflicts. The results showed 1) ILS indices (ILS-i) exceeded 50% for most deep nodes, with ILS explaining > 50% of conflicts at nodes C, D-1 and D-2, indicating ILS as the dominant factor shaping early evolutionary complexity (Fig. 2b); 2) IH had limited influence, with only node B showing a significant contribution (IH-e = 15.8%, IH-i = 38.9%), and only two nodes (nodes A and B) having IH-e > 10%. These findings quantitatively confirmed ILS as the primary driver of phylogenetic conflict during early radiation, with hybridization/gene flow playing only minor roles at specific nodes (Fig. 2c).

3.4. Quantitative analysis of ILS, gene flow, and gene tree estimation error

AU test based on dataset3 (11 samples) revealed that the top six competing topologies all exhibited supported values below 70% (ranging from 65.06% to 68.78%), with no statistically significant differences among them (p > 0.05; Fig. 3b). This pattern of low support and high conflict suggested substantial phylogenetic discordance, potentially reflecting significant deviation between gene trees and the species tree. Together with observed topological incongruence between chloroplast and nuclear phylogenies (Fig. 3a), this conflict may have been driven by ILS or ancient gene flow. Notably, the topology with the highest AU test support differed from the species tree only in the placement of Bignonieae, further confirming localized discordance. These results aligned with theoretical expectations for ILS or gene flow, where maximum topology support below 70% with non-significant differences among competing topologies typically indicated profound disturbance in coalescent history.

To further quantify the sources of conflict, we applied triplet frequency analysis to dataset3, evaluating the relative contributions of ILS, gene flow, and gene tree estimation error (Cai et al., 2021). The three factors collectively explained 39.96% of gene tree variation (R² = 39.96%), with ILS contributing 59.86%, gene flow 30.12%, and gene tree estimation error only 10.02% (Fig. 4). This quantitative partitioning showed ILS and gene flow as the two major evolutionary processes shaping gene tree conflicts in Bignoniaceae.

3.5. ILS and gene flow in early divergences of Bignoniaceae

QuIBL analyses further indicated that ILS was the main cause of deep-node conflicts (Fig. 5b). Topological conflicts at nodes A, B, and D were entirely ILS-driven. Node A primarily involved conflict between the Tourrettieae-Tecomeae clade and the Argylia-Oroxyleae-Delostoma-Catalpeae-Bignonieae-Tabebuia-Paleotropical clade, with Jacarandeae as the designated outgroup. Node B mainly reflected topological discordance between Tourrettieae and Tecomeae, relative to the larger Argylia-Oroxyleae-Delostoma-Catalpeae-Bignonieae-Tabebuia-Paleotropical clade. Node D principally showed conflicting signals for the Bignonieae-Catalpeae-Delostoma clade and the Oroxyleae-Tabebuia-Paleotropical clade, using the Argylia clade as the outgroup. Other nodes showed limited gene flow signatures and were predominantly influenced by ILS.

Gene flow was the second largest contributor to gene tree-species tree conflict (30.12%; Fig. 4). In QuIBL analysis, nodes D-3 and D-1 showed the highest proportions, both exceeding 10%, with D-3 reaching 22.41% (Fig. 5b). Node D-3 primarily resolved conflicts between Oroxyleae and the Tabebuia-Paleotropical clade, using the Bignonieae-Catalpeae-Delostoma clade as the outgroup. Node D-1 mainly involved discordant signals for Bignonieae versus the Catalpeae-Delostoma clade, with the Oroxyleae-Tabebuia-Paleotropical clade as the reference lineage. We further used PhyloNet on dataset3 to infer hybridization events (Fig. 5c). Model selection identified an optimal network configuration containing three hybridization events, revealing significant historical gene flow during early radiation. An ancient hybridization occurred between Tourrettieae (γ = 0.82) and the Delostoma clade (γ = 0.18), forming Tecomeae (Fig. 5c), though this contradicted QuIBL results showing minimal gene flow at node B (Fig. 5b). We also found that hybridization between the Delostoma-Catalpeae clade (γ = 0.35) and the Oroxyleae-Tabebuia-Paleotropical clade (γ = 0.65) produced Bignonieae (Fig. 5c), consistent with QuIBL's relatively high gene flow proportion at node D-1 (Fig. 5b). A third weak gene flow event occurred from the Argylia clade (γ = 0.22) into the common ancestor of the Oroxyleae-Delostoma-Catalpeae-Bignonieae-Tabebuia-Paleotropical clade (Fig. 5c), matching a weaker gene flow signal at node C in QuIBL (Fig. 5b). All three introgression events were supported by PhyTop validation analyses (Fig. 5a), confirming the role of gene flow in shaping key aspects of Bignoniaceae's complex evolutionary history.

3.6. Timescale and character states

We reconstructed the divergence times of Bignoniaceae (Fig. 6). The stem node of Bignoniaceae was dated to the Paleocene (61.39 Ma, 95% highest posterior density [HPD]: 56.29–66.12 Ma), and the crown node was estimated at 60.25 Ma (95% HPD: 55.08–65.55 Ma), during which the ancestor of Jacarandeae diverged. Our temporal framework revealed that the family underwent an exceptionally rapid diversification during the early Eocene, shortly after a period of pronounced global warming known as the Early Eocene Climatic Optimum (EECO, ~53–50 Ma) (Westerhold et al., 2020). The core radiation of major lineages occurred within a narrow window of less than 7 million years, from approximately 52.36 Ma to 45.51 Ma. This burst of diversification gave rise to all seven major tribal clades in quick succession, a pattern consistent with significant biological turnover and diversification events observed in other plant and animal groups around 51 Ma (Jaramillo et al., 2010).

The concentration of tribal diversification within this brief geological interval suggested that the warm, stable post-EECO climates may have created new ecological opportunities that drove the rapid diversification of Bignoniaceae. This period of rapid lineage divergence also corresponded to the emergence of distinctive morphological traits that characterized the family's major clades.

Jacarandeae retained the ancestral traits of woody habit, simple/pinnate leaves, loculicidal capsule dehiscence, and winged seeds. As lineages diversified, herbaceous habit emerged independently in three clades, Tourrettieae (unsampled Tourrettia in our study), Tecomeae, and the Argylia clade. Pinnate leaves were present in the majority of lineages, while simple leaves were consistently observed only in Catalpeae and the Delostoma clade.

Fruit dehiscence mode, one of the earliest-used characters for classifying Bignoniaceae, diversified early in the evolution of Bignoniaceae. Our analysis showed that loculicidal capsules were consistently present in Jacarandeae, Tecomeae, the Argylia clade, the Delostoma clade, and Catalpeae. In contrast, septicidal capsules characterized Tourrettieae, Oroxyleae, and Bignonieae. Fleshy indehiscent capsules occurred only in parts of the Tabebuia alliance and the Paleotropical clade.

Wingless seeds were consistently produced in the Argylia clade. Beyond this clade, seed wing presence was variable in lineages containing fleshy-fruited members: both winged and wingless seeds were observed in some taxa with fleshy indehiscent capsules, while wingless seeds also occurred in a few lineages with dehiscent capsules.

4. Discussion 4.1. Conflicts in deep phylogenetic relationships of Bignoniaceae

A fundamental dynamic in the early evolution of Bignoniaceae is reflected in a conspicuous phylogenetic discordance between the chloroplast-derived tree and the species tree reflected (Cui et al., 2025). This cytonuclear conflict is most strikingly exemplified by the divergent placements of Tourrettieae, the Argylia clade, and the Delostoma clades (Fig. S3). While the species tree robustly supports Tourrettieae as sister to Tecomeae, the chloroplast tree resolves them as more distantly related. This discordance, coupled with our detection of strong ancient gene flow from Tourrettieae to Tecomeae (Fig. 5c), suggests a possible case of chloroplast capture, where the ancestral Tecomeae lineage may have acquired its chloroplast genome via introgression from a lineage closely related to the ancestral Tourrettieae. Similarly, the conflicting phylogenetic positions of the Argylia and Delostoma clades between the nuclear and chloroplast genomes are consistent with their involvement in additional ancient hybridization events. These systematic conflicts suggest that chloroplast capture via historical hybridization was a key mechanism underlying the observed cytonuclear discordance.

In this study, nuclear datasets generated a highly supported species tree (Fig. 2a), which delineated ten major clades and provided robust justification for elevating Argylia and Delostoma to tribal status. However, behind this clear species tree lay widespread gene tree conflict, reflected in low gene concordance factors and multiple statistically indistinguishable alternative topologies (Fig. 3b). These conflicts were not mere noise but served as critical evidence for reconstructing the family's early evolutionary dynamics (Degnan and Rosenberg 2009; Xu et al., 2024). Our comprehensive analyses indicated that while both ILS and ancient hybridization played roles, ILS was likely the dominant force shaping this complex history.

The prevalence of ILS pointed to a rapid early diversification of the major lineages. This was strongly supported by our divergence time estimation (Fig. 6), which placed the tribal diversification within a narrow window of 52.36–45.51 Ma, immediately following the Early Eocene Climatic Optimum (Westerhold et al., 2020). Under such circumstances, successive speciation events occurred faster than the extensive ancestral polymorphisms in a large population could be sorted (Maddison and Knowles 2006; Feng et al., 2022). As a result, different genes retained conflicting phylogenetic signals. PhyTop analysis, which showed support values below 60% for all deep nodes, combined with the failure of the AU test to reject multiple alternative topologies, constituted typical statistical evidence of the prevalence of ILS during this period.

This evolutionary process offers a new perspective for understanding the perplexing morphological homoplasy observed in Bignoniaceae. For instance, the species tree revealed that the Delostoma clade, characterized by simple leaves, formed a sister group with Catalpeae (Gentry, 1980; Olmstead et al., 2009). However, this shared trait was unlikely the result of independent convergent evolution but was more plausibly explained as an ancestral polymorphism retained through ILS (Liu et al., 2025). In other words, simple leaves might have been a polymorphic trait in an early ancestor of Bignoniaceae, which was later randomly sorted and fixed in the lineages leading to Delostoma and Catalpeae, creating the illusion of close relatedness based on morphological similarity. Similarly, the remarkable diversity in anomalous phloem structures within the Bignonieae may also have reflected ancestral polymorphisms in developmental regulatory genes that were inherited and differentially fixed in various descendant lineages through ILS (Pace et al., 2015).

4.2. Ancient gene flow events and morphological divergence

Although ILS was the dominant factor, our tests revealed that IH played a significant role in the early evolution of Bignoniaceae. The warm, stable climates of the Early Eocene Climatic Optimum likely facilitated geographical range expansions and secondary contacts between the rapidly diverging nascent lineages, which created ideal opportunities for the ancient hybridization events detected in our analyses (Jaramillo et al., 2010). These detected signals of ancient hybridization helped explain some perplexing morphological similarities and instances of cytonuclear discordance that span major clades (Stull et al., 2023; Wang et al., 2024).

In the Bignonieae, our analysis indicated that the ancestor of Bignonieae received gene flow from the common ancestor of the Oroxyleae, the Tabebuia alliance, and the Paleotropical clade (Fig. 5c). Our analysis suggests that the ancestor of Bignonieae experienced gene flow from the common ancestor of several related clades. This introgression may explain a notable morphological anomaly: within the otherwise strictly didynamous Bignonieae, certain genera (e.g., Tanaecium) unusually possess the pentandrous stamens that are characteristic of the Oroxyleae (Gentry, 1980; Lohmann et al., 2012). Furthermore, although the Tabebuia alliance is known for its indehiscent fruits, the residual seed wings in its members were actually homologous to the seed wings found in the dehiscent capsules of Bignonieae. This suggests that the genetic basis controlling seed wing development might have been shared via gene flow early in the divergence of these lineages (Costa et al., 2019).

Our analyses detected a strong signal of gene flow from Tourrettieae to Tecomeae. This ancient hybridization event may have profoundly influenced the morphological evolution of Tecomeae. Although Tecomeae is a cosmopolitan and morphologically diverse clade where many New World taxa possessed palmately compound leaves, the presence of pinnately compound leaves within the group might be explained by this introgression. It was noteworthy that its sister group of Tourrettieae, which was also distributed in the New World, consistently retained pinnately or trifoliolately compound leaves. Therefore, we hypothesize that the genetic material acquired by the Tecomeae ancestor from Tourrettieae through gene flow may have included the potential genetic basis for developing pinnately compound leaves. This provided crucial raw material for the subsequent global diversification of Tecomeae and the evolution of diverse leaf forms, including the typically pinnate leaves of Old World lineages (Olmstead et al., 2009; Ragsac et al., 2022).

A similar pattern was observed in a gene flow event involving the Argylia clade. Distributed in the arid regions of South America, the Argylia clade shared palmate leaves with the Tabebuia alliance, a leaf form also occasionally observed in Tecomeae. This pattern of morphological sharing across distinct clades may also have originated from an ancient gene flow event (Cui et al., 2025). This further suggests that early hybridization may have played a widespread and important role in shaping key innovative traits in multiple lineages of Bignoniaceae.

5. Conclusion

Our phylogenomic analyses revealed that the early diversification of Bignoniaceae was characterized by a rapid diversification, concentrated within a narrow Eocene window (52.36–45.51 Ma). Through quantitative dissection of phylogenetic discordance using extensive nuclear gene sets, we demonstrated that ILS was the primary source of pervasive gene tree conflict, while ancient gene flow played a significant but secondary role. This complex evolutionary history, dominated by ILS and punctuated by hybridization, not only explained the challenging phylogeny but also provided a coherent framework for reinterpreting the evolution of key morphological traits in the family.

Acknowledgments

This study was supported by the Fundamental Research Funds for the Central Universities, China (NO. ZZK202502) and Discipline Crossing Foundation of School of Ecology and Nature Conservation, Beijing Forestry University, China (BH2025-JX-01). The authors thank The Huntington, the DNA and tissue bank at Kew Gardens, the Herbario Carlos Casanova Lenti, the Institute of Botany of the Chinese Academy of Sciences, the Xishuangbanna Tropical Botanical Garden of the Chinese Academy of Sciences, and Stefan Burger, Jing Gao, Shi-Hao He for their help in providing leaf material of several taxa.

CRediT authorship contribution statement

Pengpeng Yan: Writing – review & editing, Writing – original draft, Visualization, Validation, Methodology, Investigation, Conceptualization. Chang Guo: Validation, Investigation. Xingyong Cui: Writing – review & editing, Validation, Investigation, Methodology. Enze Li: Writing – review & editing, Validation, Investigation, Methodology. Yuran Bai: Visualization, Validation. Manuel R. Roncal-Rabanal: Data curation. Gangmin Zhang: Writing – review & editing, Validation, Supervision, Funding acquisition, Conceptualization. Wenpan Dong: Writing – review & editing, Writing – original draft, Visualization, Supervision, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization.

Author’s statement

All authors have read the final draft and agree to the publication, ensuring that all individuals deserving authorship have been appropriately acknowledged in the authorship contribution statement.

Data availability statement

All the genome resequencing data are deposited in the NGDC under the BioProject PRJCA055149, and the detailed information for each sample is referred to Table S1.

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