1. Introduction

Species delimitation is a fundamental cornerstone for conservation, management, and sustainable utilization of biodiversity (Hong, 2016; Wu et al., 2023). However, both conceptually and in practice, defining what constitutes a species remains highly challenging. Traditionally, species identification relied on morphological traits (Mishra et al., 2016; Shen et al., 2019). However, morphological analysis is often complicated by phenotypic plasticity—including both divergence within the same population and convergence across different species—resulting in morphological overlap (Fu et al., 2023). The advent of DNA barcoding using short-standardized sequences has significantly helped to overcome these constraints, offering a more reliable method for delineating species boundaries in complex groups (Hebert et al., 2003). Commonly used universal or standard plant barcodes include three plastid regions (matK, rbcL, and trnH-psbA) and the nuclear ITS (CBOL, 2009; Li et al., 2011). However, the discriminatory power of these fragments is rather limited for certain closely related taxa, particularly within recently diverged plant groups, mainly due to the lack of adequate genetic variation, potentially pervasive hybridization (Coissac et al., 2016; Zhang et al., 2023) and/or rapid radiation (Zan et al., 2023). Consequently, there is a pressing need for a more robust dataset of plant DNA barcodes to enhance deeper taxonomic resolution (Hollingsworth et al., 2016).

Urticaceae comprises over 2000 species and approximately 54 genera with a pantropical distribution (Plants of the World Online, 2024). Despite several phylogenetic analyses (Hadiah et al., 2008; Wu et al., 2013, 2015, 2018), some clades within the family remain poorly resolved due to the lower discriminatory power of limited DNA loci (Kim et al., 2015; Tseng et al., 2019). Studies have indicated that Debregeasia Gaudich (Fig. 1) is a monophyletic genus of the tribe Boehmerieae endemic to the Old World (Wang et al., 2020b, Wu et al., 2013, Wu et al., 2018) and distributed mainly in the Himalaya, East Asia, and Southeast Asia. Several Debregeasia species are valued for their fruits and stem fibers (Chen et al., 2003), as well as medicinal properties (Almubayedh and Ahmad, 2019). Available evidence indicates that Debregeasia originated from Eurasia and spread to Africa and Southeast Asia (Wu et al., 2018), with the current pantropical and rapid radiation likely occurring in Pliocene (Wang et al., 2020a, Wang et al., 2020b). However, these studies involved limited genome coverage and sampled few individuals per taxon. Hence, a deeper examination of the historical biogeography of Debregeasia is required. Since Debregeasia was described (Weddell, 1869), the number of accepted Debregeasia species has varied (see Table S1) among studies (Wilmot-Dear, 1988, 1994; Chen et al., 2003; Wilmot-Dear and Friis, 2012; Wang, 2016). Furthermore, the convergent evolution of morphological traits evident in Urticaceae (Wu et al., 2015) and likely in Debregeasia (Fig. S1) calls into question taxonomic assignments based on morphology alone. Hence DNA based analysis is particularly relevant in this genus, but analyses using few gene loci have met limited success in resolving relationships within the Boehmerieae tribe (Wu et al., 2013, 2018). This might be because the evolutionary trajectory of Debregeasia was markedly influenced by multiple reticulation events. Moreover, the role of biogeographical radiation in the speciation of the genus remains largely unexplored. Therefore, much greater discriminatory power marker is needed compared to previous work.

Next-generation sequencing (NGS) technologies have made it possible to generate low-coverage genomes of three organelles (plastome, mitochondrial, and nuclear) at a relatively low cost through genome skimming (Zhang et al., 2017). Plastomes and entire nuclear ribosomal DNA (nrDNA) have been proposed as suitable targets for next-generation barcodes in plants (Tong et al., 2022; Wang et al., 2022) and are often referred to as ultra-barcodes (Kane et al., 2012), organelle-scale barcodes (Yang et al., 2013), or super-barcodes (Ji et al., 2019). With more informative sites (Hollingsworth et al., 2016; Ji et al., 2019), these next-generation barcodes have improved species discrimination (Zhou et al., 2019; Guo et al., 2021; Jiang et al., 2024) and phylogenetic resolution in various plant taxa, such as the Millettioid/Phaseoloid clade from Fabaceae (Oyebanji et al., 2020), Rubiaceae (Amenu et al., 2022), the tribe Urticeae (Ogoma et al., 2022), and genus Triplostegia (Fu et al., 2023). In addition, plastome sequencing allows the detection and exploitation of hypervariable loci and lineage-specific indels, which has proven crucial in discriminating Taxus species (Liu et al., 2011, 2018). Furthermore, high-throughput sequences also offer a pathway for deciphering the organellar discordance in plant lineages, such as cyto-nuclear conflicts, by exploring 353 single-copy nuclear genes (Johnson et al., 2019). Unlike plastomes, which are maternally inherited, nuclear genomes are biparentally inherited, allowing for the detection of complex biological phenomena such as introgression, hybridization, and incomplete lineage sorting (ILS) (Ravinet et al., 2017). This capacity makes nuclear data particularly valuable for resolving relationships where plastome data may fall short. Despite the increasing application of ultra-barcodes (Ślipiko et al., 2020; Yu et al., 2022), their discriminatory power remains disputed. Also, the biological events leading to topological discordances are worth exploring. It is essential to continuously evaluate the utility of ultra-barcodes and biological complexities driving the evolutionary relationships for different plant groups (Li et al., 2021).

Here, standard DNA barcoding sampling was conducted for the nine recognized Debregeasia morphospecies, permitting comparative and integrative species delimitation and evolutionary analyses using plastome and nuclear genetic data. Our goals were to (ⅰ) characterize Debregeasia plastomes for potential informative structural variations, ⅱ) compare the discriminatory power of the three categories of barcode data using multiple species delimitation algorithms, ⅲ) explore trait evolution and identify key taxonomic diagnostic traits for monophyletic clades in Debregeasia populations and species based on the newly reconstructed phylogeny, and iv) investigate the reticulate evolution and biogeographic history of Debregeasia. The findings from this study will clarify the diversity and evolution of Debregeasia and contribute to a better understanding of evolutionary processes in Urticaceae lineages.

2. Material and methods 2.1. Sample collection

A total of 1963 individuals of Debregeasia were collected between 2018 and 2023. In addition, 60 newly sampled individuals representing the nine recognized species of Debregeasia (Fig. 1) were included for genomic analysis. Leaves from newly sampled individuals were dried and stored in silica gel for preservation. To provide a comparative framework for phylogenetic and evolutionary analysis, we also sampled 11 outgroup species (Table S2), representing six tribes within the Urticaceae family(Conn and Hadiah, 2009, Wu et al., 2013). Three datasets (Ⅰ, Ⅱ, and Ⅲ) were generated based on our research goals. Dataset Ⅰ, consisting of 51 individuals of Debregeasia and 11 outgroups, was used to characterize Debregeasia plastomes, assess barcode performance, and evaluate species boundaries in Debregeasia. For each species, 2 to 10 individuals were sampled from different locations, ensuring comprehensive coverage of their entire distribution range (Fig. 2 and Table S2). Dataset Ⅱ, comprised 270 Debregeasia individuals from 91 populations, was designed to study evolution of diagnostic morphological traits in Debregeasia (Table S3). Finally, dataset Ⅲ, which included all individuals from dataset Ⅰ, with nine additional Debregeasia individuals, was used to investigate the evolutionary relationships and biogeographic history of Debregeasia (Table S6).

2.2. DNA extraction, library preparation, and sequencing

A modified cetyltrimethylammonium ammonium bromide (CTAB) protocol (Doyle and Doyle, 1987) was used to extract total DNA from ~20 mg of the silica gel-dried leaves. Genomic DNA from each sample was then assessed for quality using a NanoDrop® ND-1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, USA), and 1% agarose gel electrophoresis was run before library preparation. The library was built using the NEBNext Ultra Ⅱ DNA Library Prep Kit for Illumina (New England BioLabs, Ipswich, Massachusetts, USA) according to the manufacturer's instructions. Sequencing was done on an Illumina DNBSEQ-T7 (Illumina, San Diego, CA, United States) platform, yielding 150 bp paired-end reads and ca. 8 Gb of clean data generated for each individual.

2.3. Data assembly and annotation

The GetOrganelle v.1.7.5.3 toolkit (Jin et al., 2020) with k-mer values 21, 45, 65, 85, 105, and 127 was used for de novo assembly of the newly generated plastome and nrDNA (18S-ITS1-5.8S-ITS2-26S) data. The complete plastome (MN189948) and nrDNA (OM892755) from Debregeasia elliptica accession De19 were used as reference genomes for both assembly and annotation. The assembled plastomes were visualized and filtered to generate a complete circular genome in Bandage v.0.8 (Wick et al., 2015). The inverted repeat (IR) regions were identified by automatically flanking the IR using the function "find repeat" in Geneious v.8.1 (Kearse et al., 2012). Further curation of the annotated plastomes was done to confirm the start and stop codon, and tRNAscan-SE v.1.21 (Schattner et al., 2005) was used in default settings to verify the tRNA genes. The newly assembled and annotated genomes were deposited at the China National GenBank for accession numbers (Tables S2 and S6).

2.4. Plastome structural variation, divergence, and mutational hotspot analyses

We characterized the features of 51 plastomes, including genome size and structural elements, gene content [protein-coding (PCG) genes, tRNA, rRNA], and GC content. Comparative analyses of IR expansion and contraction at the four junctions of the sampled plastomes were conducted in Geneious and plotted in IRscope (Amiryousefi et al., 2018). Further, gene rearrangements were investigated in Mauve v.2.4.0 (Darling et al., 2010) with default settings, while sequence divergence analysis was performed in mVISTA in Shuffle-LAGAN mode (Frazer et al., 2004) using D. elliptica De19 as the reference plastome. We further assessed the hypervariable regions based on nucleotide diversity and the variation of parsimony informative sites and nucleotide diversity of all Debregeasia plastomes as an additional tool for species barcode assessment. The coding and non-coding regions of the sampled plastomes were extracted and separately analyzed. Parsimony informative sites were estimated in PAUP* v.4.0b10 (Swofford, 2002) while the nucleotide diversity analysis was conducted using a sliding window analysis in DnaSP v.5 (Librado and Rozas, 2009), with a step size of 200 bp and a window length of 600 bp. The top ten sequences based on nucleotide diversity were selected as hypervariable regions. We generated a physical circular map of the plastome in Chloroplot (Zheng et al., 2020).

2.5. Species discrimination analysis

We aligned the newly generated plastomes and nrDNA data separately using MAFFT v.7.22 (Katoh and Standley, 2013) with the LINSI algorithm; from this, multiple data matrices were generated to explore species delimitation among the sampled individuals. Six regions were extracted from the plastome using Geneious: coding sequences (CDS), non-coding sequences (nCDS), the hypervariable ycf1 region, and three universal standard DNA barcode sequences—rbcL, matK, and trnH-psbA. Two nuclear datasets were extracted: ITS (ITS1-5.8S-ITS2), and nrDNA (18S-ITS1-5.8S-ITS2-26S). These datasets were combined to create 12 data matrices that were then classified into three barcode categories: ultra-barcodes (whole plastome, whole plastome + nrDNA, CDS, and nCDS), specific (nrDNA, ycf1, matK + rbcL + trnH-psbA, and matK + rbcL + trnH-psbA + ITS), and standard barcodes (matK, rbcL, trnH-psbA, and ITS). Each data matrix was aligned separately using MAFFT and manually adjusted in Geneious, while the most suitable nucleotide substitution model was determined in jModelTest 2 (Darriba et al., 2012). Sites that were parsimony-informative, constant, and variable were determined using PAUP* v.4.0b10 (Swofford, 2002), with gaps treated as missing data (Table 1).

Two types of broadly used species discrimination methods, i.e. tree-based and distance-based (see details in supplementary file 1), were employed in each of the three barcode data categories. Tree-based methods use a phylogenetic tree approach to calculate the variation along the branches to delimit species, and four were used here: Maximum likelihood (ML), Bayesian Inference (BI), Generalized Mixed Yule Coalescent (GMYC: Pons et al., 2006), and Poisson Tree Processes (PTP: Kapli et al., 2017). Conversely, distance-based methods mainly rely on the level of differences between sequences to calculate intraspecific and interspecific thresholds; of these, we employed three: Barcoding gap (BG), Automatic barcode gap discovery (ABGD; Puillandre et al., 2012) and assembly species by automatic partitioning (ASAP; Puillandre et al., 2021). All these methods use assembled data, and so a third type of method, skmer analysis (Sarmashghi et al., 2019), was performed on the unassembled reads so that its power could be compared to tree and distance-based methods. The species discrimination power of ML, BI, and skmer was determined by whether all representatives of a given species clustered into a monophyletic group, whereas in BG analysis, a species was successfully identified when the minimum interspecific distance was higher than its maximum intraspecific distance. Other analyses relied on matches, defined as occurring when morphospecies correspond to operational taxonomic units (OTUs, defined as comprising all individuals of one species but no others). Hence, GMYC, bPTP (Bayesian Poisson Tree Processes), mPTP (multi-rate Poisson Tree Processes), ABGD, and ASAP were evaluated based on the number of matches obtained, following Magoga et al. (2021).

2.6. Morphological trait analysis

We used dataset Ⅱ to comprehensively analyze 14 morphological traits that vary within Debregeasia, measuring two to three specimens in each sampled population (Table S3). All morphological terms used for both the qualitative (leaf base, leaf margin, leaf tip, leaf blade abaxial colour, leaf shape, perianth type, and pubescence of branch and petiole) and quantitative (number of secondary veins, inflorescence length, peduncle length, leaf area, leaf length, leaf width, and leaf weight) traits followed Harris and Harris (1994) and Flora of China (Chen et al., 2003). Quantitative traits were measured manually using a vernier calliper (SHAHE caliper, Wenzhou, Zhejiang, China). The leaf area was calculated using ImageJ software (Abràmoff et al., 2004), and leaf weight using ME-T analytical balance (Mettler-Toledo, Columbus, Ohio, USA). To account for within-individual variation, we measured leaf traits from three leaves of the middle stem, and the longest spreading inflorescences of the lower, middle, and upper stem parts. We calculated Spearman's rank correlation coefficients to define the relationship between characters and performed Principal Coordinates Analysis (PCoA) using 'ggfortify, ' 'ggrepel, ' and 'ggplot2' packages on the Gower distances from character trait data generated using the 'gowdis' function in 'FD' package to visualize the variations and clustering in the morphological traits of the Debregeasia species. All analyses were done in R (R Core Team, 2024).

Additionally, we performed ancestral state reconstruction (ASR) based on seven qualitative characters, with Debregeasia australis scored following Wilmot-Dear and Friis (2012). We used maximum parsimony (MP) in Mesquite v.2.75 (Maddison, 2007) on the newly reconstructed phylogeny based on the concatenated plastome and nrDNA data while treating character states as unordered. The character coding and data matrices are presented in Tables S4 and S5, respectively.

2.7. Single-copy nuclear genes assembly and alignment

We extracted the universal set of 353 single-copy nuclear genes from genome skimming data using Captus v.0.9.83 (Ortiz et al., 2023) following the recently published pipeline (https://github.com/edgardomortiz/Captus) based on dataset Ⅲ (Table S6). MEGAHIT v.1.2.9 (Li et al., 2015) was used with the "assemble" function to combine paired reads into contigs. We then extracted the target genes using the "extract" function and aligned the extracted markers in MAFFT v.7.22.

2.8. Phylogenetic inference and reticulate evolution analysis

Phylogenetic relationships were inferred using ML and the species tree for the single-copy nuclear genes matrix. For whole plastome data, we relied solely on the ML tree. This was based on the performance of this method and the barcode in species delimitation in prior analysis. ML followed the procedure described in section 2.5 (see details in supplementary file 1), while the species tree was inferred using the multispecies coalescence (MSC) model in ASTRAL-Ⅲ (Zhang et al., 2018). The MSC approach can detect biological events such as incomplete lineage sorting and hybridization causing topological conflicts or discordance. Thus, we compared the plastid and nuclear trees from ML and ASTRAL-Ⅲ for gene and species trees, respectively. These alternative tree topology frameworks were used to provide insights into the reticulate evolutionary histories of the genus Debregeasia. We quantified the gene concordance using PhyParts v.0.0.1 (Smith et al., 2015). Next, the quartets (SNaQ) method implemented in the Julia package PhyloNetworks v.0.6.0 (Bezanson et al., 2017; Solís-Lemus et al., 2017) was employed to infer phylogenetic networks from the species networks, for five separate runs with the maximum number of hybridization events varying from 0 to 4. The best phylogenetic network with the smallest pseudolikelihood value was identified as showing the most likely pattern of past hybridization events. A coalescent simulation analysis was performed to measure the goodness-of-fit of the coalescent model and whether ILS explains the gene tree discordance. If the simulated gene trees correspond well to the empirical tree, ILS may explain discordance (Yang et al., 2020; Morales-Briones et al., 2021).

2.9. Molecular dating and ancestral area reconstruction

We conducted molecular dating on the newly reconstructed ML tree of 352 single-copy nuclear genes based on a penalized likelihood (PL) in treePL (Smith and O'Meara, 2012). We utilized the ML tree from the single-copy nuclear genes matrix due to its higher support values and resolved relationships. The nodes were calibrated using two selected fossil points from the outgroups, Forsskaoleeae (33.9 Ma) and Girardinia (11.6 Ma). Additionally, one secondary calibration (crown divergence) for Urticaceae 68.7 Ma (56.2–81.7 Ma) was applied, as inferred from previous study based on the fossil (Wu et al., 2018). Following the empirical guide described by Maurin (2020), the smoothing parameter was determined using the cross-validation option, and priming was used to determine the best optimization scores. We used TreeAnnotator v.2.6.7 (Bouckaert et al., 2014) to summarize the dated bootstrapped trees and produce the confidence intervals (CI) for a consensus tree. The results were visualized in FigTree v.1.4.0 (Rambaut, 2009).

We conducted the biogeographical analyses using Bayesian Binary MCMC Analysis (BBM) implemented in RASP (Yu et al., 2020) based on the dated ML tree from single-copy nuclear genes data with the outgroup pruned. The number of maximum areas for ancestral nodes was set at two, representing the maximum number of areas observed in extant Debregeasia species. The distribution area was divided into seven regions: A, Southern China; B, Third Pole (Tibetan Plateau and adjacent mountain regions); C, Southeast Asian archipelago; D, Indian subcontinent; E, Australian archipelago; F, Yemen, Saudi Arabia, and Ethiopia; G, Taiwan (China) and Japan. These areas were extrapolated from the literature in Flora of China and plants of the world based on the geological and climatic history, distribution, and endemism of Debregeasia (Chen et al., 2003).

3. Results 3.1. Plastome structural variation, sequence divergence, and hypervariable regions

Plastomes of all 51 individuals of the nine Debregeasia species examined had a quadripartite structure (Fig. 3A) and exhibited limited intraspecific variation in size (Table S7). The lengths of the large single-copy (LSC), small single-copy (SSC), and IR regions had ranges of 84,904–90,893 bp, 18,532–19,217 bp, and 22,536–25,689 bp, respectively (Table S7). There was marginal variation in total GC content, ranging from 36.30 to 36.80%. However, the GC content for IR regions was higher (42.6–42.7%) than that of LSC (34.0–34.5%) and SSC (29.4–30.2%). Debregeasia plastomes are highly conserved, and encode between 125 and 130 genes, of which 80–85 are protein-coding, eight encode rRNA, and 37 tRNA genes.

The plastomes of Debregeasia species studied have experienced independent boundary shifts in the IR both between (Fig. 3C) and within species (Fig. S2). The LSC/IRb border is embedded in the rps19 gene (of which 27–43 bp is located within IRb), causing an independent boundary shift (Fig. S2). The IRb/SSC boundary is located within the ndhF gene (of which 21–36 bp is located within IRb), while the IRa/LSC boundary is within the rpl2–trnH-GUG region (Fig. S2). The boundary SSC/IRa is conserved throughout the studied plastomes and is located within the ycf1 coding gene (of which 970–1370 bp is located within IRa).

There were no significant evolutionary gene rearrangements or inversions detected in the studied plastomes (Fig. S3). Similarly, sequence variation analysis of Debregeasia plastomes revealed low sequence divergence within and between species (Fig. S4). Nucleotide diversity (π) ranged from 0 to 0.0422 across the plastomes, and the most hypervariable regions included rps16–trnQ-UUG, ycf4–cemA, trnT-GGU–psbD, rpl32–trnL-UAG, and ycf1 (Table S8 and Fig. 3B). Parsimony-informative sites for all regions ranged from 0 to 516 sites (Figs. S5 and S6; Tables S9 and S10). Because ycf1 exhibited both the highest π value and parsimony-informative sites, its species discriminatory power was further investigated.

3.2. Characteristics of data matrix

Of the 12 data matrices, an ultra-barcode (plastome + nrDNA) and standard barcode (ITS) had the longest (170,286 bp) and shortest (741 bp) alignment length, respectively. However, one standard barcode (ITS) had the highest percentage of variable (40.85%) and parsimony-informative (38.30%) sites, in contrast to another standard barcode (rbcL), which had the smallest percentage (2.62%) of both variable and parsimony-informative sites (Table 1).

3.3. Species discrimination 3.3.1. Tree-based species discrimination

BI and ML topologies derived from all 12 data matrices for the nine species were congruent regarding which species were monophyletic; in addition, these methods outperformed other tree-based methods used in this study (Table 2 and Figs. S7–S18). Six species were monophyletic in all analyses, whereas Debregeasia longifolia, D. hekouensis, and D. libera were non-monophyletic, and hence could not be discriminated by any barcode. Three data matrices from the ultra-barcodes (plastome + nrDNA, plastome, and CDS) and two specific barcodes (nrDNA and matK + rbcL + trnH-psbA + ITS) showed the highest discriminatory power (67%; six of nine species discriminated). The standard barcodes matK, rbcL, and ITS gave an average performance of 56% (five of nine species discriminated), while trnH-psbA had the least resolution, with only two out of nine species discriminated (22%; Table 2). Additionally, we noted significant topological incongruence between the plastome vs. nrDNA results (Fig. 4A and B).

We observed contrasting discriminatory powers of the data matrices with the different algorithm thresholds. By delimiting five species, the ycf1 specific barcode and CDS ultra-barcodes had the highest discriminatory power in bPTP and mPTP, respectively. In contrast, the nrDNA specific barcode and all the four standard barcodes (ITS, matK, rbcL, trnH-psbA) had the least resolution, delimiting no species in mPTP. In mGMYC, the plastome, plastome + nrDNA, and ITS performed slightly below-average, delimiting four out of nine species, whereas the CDS ultra-barcodes could not delimitate any species (0/9). In sGMYC, the standard barcode matK gave the highest resolution (5/9), whereas both the ultra-barcodes (plastome and plastome + nrDNA) and ITS (standard barcode) delimited 4/9, and the others even fewer. In putative OTU analysis (Tables S11–S14), plastome + nrDNA yielded similar ML OTUs in both sGMYC and mGMYC, whereas the plastome ultra-barcodes alone gave more in mGMYC compared to sGMYC (Table S11). The CDS ultra-barcodes stand out with the highest OTUs (38) in sGYMC, whereas ITS exhibited the lowest (1). Both matK and rbcL resulted in more ML clusters in sGMYC than mGMYC (Table S11). The bPTP analysis showed a high acceptance rate and maximum split for the specific barcode (matK + rbcL + trnH-psbA) (Table S12), and most OTUs in nrDNA specific barcode (22) and the fewest (2) in ITS standard and CDS ultra-barcodes. In contrast, mPTP analysis gave the most OTUs in the CDS ultra-barcodes (10) and the fewest in the nrDNA-specific barcode (1), exhibiting a significant discrepancy between the two methods (Table S13). Plastome ultra-barcodes gave the same number of OTUs (8) in both mPTP and bPTP, whereas the four standard barcodes (matK, rbcL, trnH-psbA, and ITS) gave fewer (between 2 and 6) in both mPTP and bPTP.

3.3.2. Distance-based species discrimination

The discriminatory power of the distance-based analyses varied considerably, identifying 1–6 species of the nine morphospecies. In barcode gap analysis, a specific barcode (matK + rbcL + trnH-psbA + ITS) had the highest discriminatory power (6/9 species; 67%); next best were the four ultra-barcodes (plastome + nrDNA, plastome, CDS, and nCDS) and a standard barcode (ITS), each of which identified five of the nine species (Fig. S19A). In ABGD analysis, five barcodes, including three standard barcodes (matK, rbcL, and ITS), one ultra-barcodes (plastome + nrDNA) and one specific barcode (matK + rbcL + trnH-psbA), performed best (6/9; 67%) whereas nrDNA performed the worst (1/9; 11%). Using the ASAP algorithm, the specific barcode (matK + rbcL + trnH-psbA) and standard barcode (rbcL) performed best (5/9; 56%), while other barcodes discriminated up to 4/9 taxa (Table 2). The number of generated OTUs varied across ABGD analysis with the different prior intraspecific divergence both in initial and recursive partitions (Tables S14 and S15), with the standard barcode (ITS) generating the highest number of OTUs (25) in recursive partition and the fewest in the ycf1 specific barcode (2) in initial partition, both at 0.001 prior intraspecific divergence. The three models, K80 Kimura, Jukes-Cantor (JC69), and simple distance, produced comparable OTUs for each dataset. The number of discriminated OTUs similarly varied across different ASAP score values (Table S16), with the nCDS specific barcode generating the most OTUs (33), while matK + rbcL + trnH-psbA, ycf1, and ITS generated the fewest (2). The skmer method equally discriminated up to six out of nine species (67%; Fig. 4C).

3.4. Comparison of species resolution for different analytical methods

Species discrimination rates varied across analysis method and barcodes (Fig. S19B). ML and BI outperformed other methods, with better resolution in six out of twelve barcodes and a maximum resolution of 67% (6 out of 9 species discriminated). In contrast, mPTP performed the worst, with no species discriminated by any of the standard barcodes (0% resolution). Among the distance-based analyses, the barcode gap method showed a better resolution power (> 56% in nine of twelve barcodes) than ASAP and ABGD. Equally, skmer analyses resulted in the best resolution (67%), similar to the maximum performance achieved in BI and ML methods.

3.5. Morphological evolution

Spearman's rank test showed varying relationships among the examined morphological characters, with eight characters recovered as uncorrelated (Table S17). PCoA of all 14 morphological traits showed a better clustering than when the qualitative and quantitative traits were separately analyzed, with four morphospecies (Debregeasia elliptica, D. libera, D. squamata, and D. wallichiana) displaying distinct clusters. A similar pattern was resolved when only the eight selected uncorrelated traits were analyzed (Fig. 5).

ASR revealed multiple cases of homoplasy involving key diagnostic morphological characters (Figs. S20A–B). Leaf shape barely defined half of the species, with the rest of the characters exhibiting rampant homoplasy. The leaf tip, base, perianth type, and margin showed high homoplasy with shared character states across species (Figs. S21–S25). Three synapomorphies were found in Debregeasia elliptica (rounded leaf base, finely denticulate leaf margin, and elliptic leaf shape), two in D. saeneb (thickly snow-white tomentose and sparsely spreading hirtellous), and D. squamata (cordate leaf base, fleshy spreading protuberances, and mixed appressed hair), and a single trait in D. wallichiana (shortly acuminate leaf tip).

3.6. Sequence characteristics of single-copy nuclear genes

We recovered 352 exons of the universal set of 353 targeted loci. The length of each locus after trimming ranged from 195 to 3657 bp, with GC content ranging from 39.0% to 68.7%. These genes were assembled into a concatenated matrix of aligned length 277,835 bp, with 153,277 constant sites, 5822 parsimony-uninformative, and 66,336 parsimony-informative sites. GC content was 46.8%.

3.6.1. Phylogenetic resolution and reticulation patterns

The nuclear tree was largely consistent with the plastome tree but with conflicting phylogenetic positions detected in the concatenated Debregeasia tree (Fig. S26), where D. hekouensis and D. libera placements varied significantly. In both data matrices (Figs. S27–S29), Debregeasia is supported as monophyletic (BS = 100%; PP = 1.00), excluding D. wallichiana, which clustered with Coussapoa contorta and Maoutia setosa with strong support (BS = 100%; PP = 1.00). The single-nuclear copy data matrix recovered six Debregeasia species as fully supported monophyletic clades (D. elliptica, D. saeneb, D. orientalis, D. squamata, D. longifolia, and D. australis), whereas the plastome recovered only five (with D. longifolia in two clades). The conflict analysis based on the nuclear gene trees showed that most gene trees (144/118) supported the monophyly of Debregeasia, suggesting potential evolutionary processes such as ILS and hybridization (Fig. 6A). In line with this, phylonet recovered two phylogenetic networks of possible hybridizations from the four runs with (Likelihood: −5463.62) (Fig. 6C) and (Likelihood: −5286.41) (Fig. 6D) where one hybridization event was identified as the best-fit model based on the smallest pseudolikelihood value (Fig. 6C). This network suggests that 68.3% of ancestral lineage in Debregeasia are from D. australis and 31.7% from another ancestral sister Urticaceae species (possibly the close relatives used as outgroups). Coalescent simulation analysis between the empirical and simulated gene-to-gene tree at the species level showed much overlap between the plastid and simulated gene tree (Fig. 6B). Although the simulations did not rule out ILS, SNaQ analysis suggests that other alternatives, including gene flow and plastid capture, might explain the conflict.

3.6.2. Molecular dating and ancestral area reconstruction

Our dating analysis suggests that Debregeasia started to diversify at ca. 12.82 Ma (crown-age; 95%HPD: 11.54−14.63Ma) (Fig. 7D) in the mid-Miocene era of the Neogene period. Our analyses further indicated that the biogeographical ancestral area of Debregeasia may be the Australian archipelago (Prob: 51.75%), and that the genus originated during the mid-Miocene. A variance event subsequently split Debregeasia into two lineages. The first lineage consists of the Australian species, D. australis, whereas the second lineage consists of the remaining species from Asia and Africa (Fig. 7A and B). At least six dispersals and one variance event occurred during the Miocene, which influenced the present-day distribution pattern of the genus (Fig. 7C).

4. Discussion 4.1. Plastome characteristics and implications for DNA barcoding

Plastomes of Debregeasia were found to be highly conserved and exhibited typical quadripartite structure, similar to other Urticaceae plastomes such as Oreocnide (Wu et al., 2022), Urtica (Li et al., 2022), members of Urticeae (Ogoma et al., 2022) and other photosynthetic angiosperms (Wang et al., 2020b; Amenu et al., 2022). No inversions were detected, but several species had independent IR-SC boundary shifts, as detected for other taxa (Park et al., 2018; Ogoma et al., 2022). In Debregeasia, the maximum possible species discrimination based on plastomes was 6/9, with D. libera, D. hekouensis, and D. longifolia similarly non-monophyletic in nrDNA (Fig. 4). Possible reasons for this are discussed below, but in practice, barcode performance could therefore not exceed six discerned species, or 67%. Hence, plastome barcodes achieved the best possible performance alone or combined with nrDNA and also as coding sequences only (CDS) (Fig. S19B). This confirms the reported utility of complete plastomes as barcodes (Fu et al., 2022; Xu et al., 2022; Zhang et al., 2023), especially where species are difficult to resolve, as also shown in Ficus (Bruun-Lund et al., 2017), Diospyros (Li et al., 2018), Panax (Ji et al., 2019), and Taxus (Fu et al., 2019). Although standard barcodes have performed well for Taxus (Liu et al., 2018), Alnus (Ren et al., 2010), Dendrobium (Singh et al., 2012), and Pedicularis (Yu et al., 2015), their performance may be suboptimal in recently diverged species or taxonomically complex plant taxa (Liu et al., 2022; Yu et al., 2022). Here, most standard individual barcodes performed relatively poorly, resolving a maximum of five Debregeasia species, indicating a need for specific or ultra-barcodes. However, combining the standard barcodes with others (matK + rbcL + trnH-psbA + ITS) offered an improved resolution, discriminating all six monophyletic species in barcoding gap analysis and both ML and BI. This indicates that discrimination success is improved by increasing the number of variable characters examined.

Taxon-specific markers offer a potential compromise between the expense involved in ultra-barcodes, such as the whole plastome, and the potential problem of limited genetic variation for standard barcodes. We identified five mutational hotspots for Debregeasia: rps16–trnQ-UUG, ycf4–cemA, trnT-GGU–psbD, rpl32–trnL-UAG, and ycf1. These mutational hotspots were similar to hotspots for other Urticaceae (Li et al., 2021; Ogoma et al., 2022). These regions had π values ranging from 0.0223 to 0.0421. The ycf1 hotspot offered similar discrimination power (56%) to three of the standard barcodes (matK, rbcL, and ITS) and was better compared to trnH-psbA (22%). Therefore, it can serve as an alternative where these standard barcodes fail. Considering the cost and time constraints of plastome sequencing, these gene regions are cost-effective and efficient for future population genetic studies on Debregeasia, contributing to the expanding dataset for taxon-specific barcodes.

4.2. Species delimitation of Debregeasia

Species delimitation success varied depending on the data and methods used (see details in the supplementary file 1, section 4.2.1), similar to previous studies (Liu et al., 2018; Magoga et al., 2021). Complete plastome data indicated that five of nine Debregeasia species are monophyletic, with three species, D. hekouensis, D. libera, and D. longifolia, polyphyletic; the remaining species, D. wallichiana should be removed from the genus. However, phylogenies based on the universal set of 353 nuclear genes proved more powerful in recovering six monophyletic species, including D. longifolia. In phylogenetic trees based on both plastome and single-copy nuclear genes, D. wallichiana formed a clade with Coussapoa contorta and Maoutia setosa with maximum support (Figs. S2–S4) and was hence separated from all other Debregeasia species. Therefore, Debregeasia, as currently circumscribed, is polyphyletic; this is a new finding, as previous studies did not include D. wallichiana (Wu et al., 2013, 2018; Wang et al., 2020a). The taxonomic treatment of D. wallichiana has been unstable: first described as Missiessya wallichiana (Weddell, 1854), it was later placed in both Debregeasia and Morocarpus, the latter of which was later synonymized into Debregeasia (Wilmot-Dear, 1988). Debregeasia wallichiana has very similar morphology to other Debregeasia species, such as the cordate leaf base and the greenish-gray color of the abaxial leaf blade, as noted in the Flora of China (Chen et al., 2003). However, it differs in leaf shape, leaf margin, and leaf tip traits, which are key diagnostic characters in Debregeasia. Although we did not sample the full D. wallichiana species range due to sampling challenges around the southern Himalaya, the seven individuals examined did cover a wide geographical range (Fig. 2). Therefore, we strongly recommend that D. wallichiana be removed from Debregeasia, and restored to Morocarpus as M. wallichiana.

Further, ASR revealed key synapomorphies for four monophyletic species, including three traits in Debregeasiaelliptica (rounded leaf base, finely denticulate leaf margin, and elliptic leaf shape), two each in D. saeneb (thickly snow-white tomentose and sparsely spreading hirtellous), and in D. squamata (cordate leaf base and fleshy spreading protuberances mixed with appressed hair) and one trait in D. wallichiana (shortly acuminate leaf tip). In contrast, the complex of D. hekouensis, and D. libera showed rampant overlapping morphological traits with similar leaf base, tip, margin, and abaxial blade color. These two species were already known to share many similarities in general growth, habit, and reproductive structures (Chen et al., 2003). D.ebregeasia libera and D. hekouensis have similar leaf shape (broadly lanceolate) (Fig. 1) and can only be distinguished based on inflorescence length, with D. hekouensis being longer. Therefore, discrimination of these two species is morphologically complex, and they might not be truly distinct, in which case no amount of barcode data will separate them. Case in point, plastome, nrDNA sequences, and the universal set of 353 nuclear genes were unable to distinguish between these species. To fully test their status might require the additional power offered by sequencing more of the nuclear genome and samples, or perhaps all of it e.g. (Wu et al., 2024), to better reflect their evolutionary history. For instance, taxon-specific datasets, such as the whole genome (Gu et al., 2023), have been helpful in resolving critical species complexes) and could be tried here. Phenotypic clustering and ancestral state reconstruction indicate that no single character of Debregeasia's main traits is a useful tool for classifying this species complex. Hence, the relationship and boundaries of these two species should be one of the priorities for future research in Debregeasia.

Our data indicated that Debregeasia began diversifying in the mid-Miocene (12.82 Ma) (Fig. 7), with present distributions arising from at least six dispersal events. Southern China is identified as a key source of dispersal events into Southeast Asia, the Indian subcontinent, parts of Taiwan (China) and Japan, and Africa. Our dates for Debregeasia are generally older than those of previous studies that did not include D. australis (Wu et al., 2018). Insufficient sampling within a group can result in an underestimation of the crown age of the plant group (Lu et al., 2020). The convergence of land masses after the mid-Miocene is likely to have facilitated the dispersal of Debregeasia throughout Southeast Asia, with the Wallacean islands potentially serving as stepping-stones among Southeast Asian islands (Wu et al., 2019; Low et al., 2021). Additionally, the geological stability of Southeast Asia during this period suggests that the dispersal may have been facilitated by the intensification of the East Asian summer monsoon between 15 and 10 Ma (An et al., 2001). Asian monsoons have played a significant role in the diversification patterns of several plant groups, including Primulina (Gesneriaceae) (Kong et al., 2017), Melocanninae (Zhou et al., 2022), Oreocnide (Urticaceae) (Wu et al., 2022), and Cymbidium (Orchidaceae) (Chen et al., 2024). Further investigation into the mechanisms behind the long-distance dispersal of Debregeasia from Asia through Southeast Asia to Australia is necessary to understand the potential asymmetry in this exchange (Peng et al., 2021). It is, therefore, highly likely that long-distance dispersal is the main factor influencing the present distribution patterns of Debregeasia, with seeds possibly being dispersed across oceans over considerable distances (Wu et al., 2018). Also, analysis of reticulate evolution showed that hybridization might have been pervasive in Debregeasia, influenced perhaps by speciation patterns and overlaps between species ranges, contributing to the complex evolutionary history of the genus.

4.3. Phylogenetic resolution and reticulate evolution in Debregeasia

Phylogenomic analyses based on single-copy nuclear genes have provided valuable insights into the phylogenetic relationships of many plant groups where plastid genomes alone failed (Zhang et al., 2021; Gu et al., 2024; Yan et al., 2024). In this study, we utilized 352 single-copy nuclear genes to resolve phylogenetic relationships in Debregeasia. Our results recovered six monophyletic clades with stronger support than plastome and nrDNA sequences. Notably, neither the plastome nor the nrDNA sequences resolved any of D. hekouensis, D. libera, or D. longifolia as monophyletic, which was partially resolved in the universal set of 353 nuclear genes. These instead formed a mixed clade together but with differing relationships among individual accessions (Fig. 4). Such cyto-nuclear discordance can be associated with hybridization, introgression, and incomplete lineage sorting (Xiang et al., 2015; Folk et al., 2017; Loiseau et al., 2021). Coalescent simulation analysis revealed that almost none of the simulated trees supported the plastid topology Fig. 6B. This implies that ILS alone cannot account for the cyto-nuclear discordance observed in Debregeasia, which, therefore, must be driven by gene flow and plastid capture.

The overlapping geographical distributions and similar flowering and fruiting times of Debregeasia species suggest that species within the genus have adequate opportunity for interspecific gene flow. For example, D. longifolia, D. hekouensis, and D. libera are sympatric in southern Yunnan and show high conflict and gene tree discordance, as noted above, which might reflect introgression. Hybridization is common among angiosperms (Mallet, 2007) and can blur species boundaries (Petit and Excoffier, 2009; Hollingsworth et al., 2016). Introgression is plausible because D. hekouensis, D. libera, and D. longifolia occur sympatrically (Fig. 2) with overlapping flowering time (August–December), and where this occurs, no single barcode may be effective at species delimitation. The skmer analysis resolved D. libera and D. hekouensis as a single intermixed clade, sister to a clade of D. longifolia material from southern Yunnan. Sister to these was a clade comprising all other D. longifolia material, which was all from northwestern Yunnan. These two regions differ substantially in climate, with the southern region much warmer and wetter, which might explain the evolutionary patterns behind the two clades. However, clade membership varies according to analysis and region (Fig. 4), so further comprehensive analysis with wide sampling and new genomic data is needed to investigate plastome capture and reticulation patterns among these three species.

5. Conclusions

We examined plastome structural variations and investigated the efficacy of plastome, nrDNA, and single-copy nuclear genes for resolving species boundaries and evolutionary relationships based on within and between species variation in Debregeasia. Our study identified varying genetic and morphology-based species boundaries. We revealed multiple cases of overlapping diagnostic characters in Debregeasia and uncovered key synapomorphic traits for D. elliptica, D. wallichiana, D. saeneb, and D. squamata. Compared to the standard barcodes, the complete plastome proved to be an effective plant ultra-barcode for distinguishing closely related plant species. Full discrimination was not achieved, almost certainly because three species were genuinely non-monophyletic, with accessions intermixed for both plastome and nrDNA. However, single-copy nuclear genes resolved this complex and suggested six monophyletic clades in Debregeasia, excluding D. wallichiana, which our data indicates should be removed from Debregeasia. Our results indicate that hybridization might have been pervasive in Debregeasia, influenced perhaps by speciation patterns and overlaps between species ranges, providing an exciting opportunity for future studies to investigate further reticulation events in this group. In addition, estimated divergence times of Debregeasia suggest that long-distance dispersal has played an important role in forming current intercontinental distributions, with seeds dispersed across oceans. This study serves as an exemplary case for the utility of recently proposed deep genome skimming methods to retrieve nuclear and organelle genes, enabling robust phylogenetic and biogeographic reconstructions. Additionally, our study demonstrates that extensive exploration of nuclear genes can circumvent the challenges faced by ultra-barcodes and can pave the way for clarifying species relationships in evolutionarily intractable taxa.

Acknowledgements

We extend our gratitude to Mr. Tao Liu and Mr. Xuewen Liu for their effort in field sampling. We are sincerely thankful to Prof. Chia-Jui Chen, Prof. De-Zhu Li, Prof. Ran Nathan, Associate Prof. Ya-Huang Luo and Prof. Debabrata Maity for their valuable advice and insightful comments during this study. Special thanks to the Royal Botanic Garden Edinburgh (E) and the Royal Botanic Gardens, Kew (K), for providing DNA material. This research was funded by the National Natural Science Foundation of China (42171071, 42211540718), Key Research Program of Frontier Sciences, CAS (ZDBS-LY-7001), Yunnan Fundamental Research Projects (202401AT070190, 202201BC070001), Top-notch Young Talents Project of Yunnan Provincial "Ten Thousand Talents Program" (YNWR-QNBJ-2018-146, YNWR-QNBJ-2020-293), and CAS "Light of West China" Program. Jie Liu and Zeng-Yuan Wu were also supported by the China Scholarship Council (202304910135 and 202304910138) for a one-year study at the University of Toronto, Canada. Molecular experiments were conducted at the Laboratory of Molecular Biology and data analysis was partially facilitated by the iFlora HPC Center (iFlora High-Performance Computing Center), Germplasm Bank of Wild Species, Kunming Institute of Botany, Chinese Academy of Sciences.

CRediT authorship contribution statement

Amos Kipkoech: Writing – review & editing, Writing – original draft, Methodology, Formal analysis. Ke Li: Writing – review & editing, Methodology, Formal analysis. Richard I. Milne: Writing – review & editing. Oyetola Olusegun Oyebanji: Writing – review & editing, Formal analysis. Moses C. Wambulwa: Writing – review & editing. Xiao-Gang Fu: Writing – review & editing. Dennis A. Wakhungu: Writing – review & editing. Zeng-Yuan Wu: Writing – review & editing, Supervision, Funding acquisition, Data curation, Conceptualization. Jie Liu: Writing – review & editing, Supervision, Funding acquisition, Formal analysis, Data curation, Conceptualization.

Data availability

All DNA sequences generated in this study have been submitted to China National GeneBank for accession numbers presented in Tables S2 and S6.

Authors statement

All co-authors have reviewed and approved the final version of the manuscript, and those entitled to authorship have been named accordingly. This paper contains original content, and no other studies employing this database have been submitted elsewhere. The document includes all of the article's contents, raw data information, and protocols.

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