1. Introduction

The mitochondrial cytochrome oxidase 1 (CO1) gene serves as a highly effective DNA barcode for species delimitation in insects (Foottit et al., 2008; Adeniran et al., 2019), fishes (Ward et al., 2005; Decru et al., 2016), birds (Kerr et al., 2009; Tizard et al., 2019), as well as other animal groups (Govender et al., 2022; Velo-Antón et al., 2022). In plants, however, the CO1 gene exhibits significantly lower species delimitation efficacy, partially attributable to the low substitution rates characteristic of plant mitochondrial genomes (Cho et al., 2004; Fazekas et al., 2009; Lewin et al., 2022). Consequently, multiple plastid DNA regions (trnL-F, rpoC1, rbcL, matK, trnH-psbA) and the nuclear ribosomal internal transcribed spacer (nrITS), and their combinatorial applications have been prosed as candidate DNA barcodes for plants (Chase et al., 2007; Kress and Erickson, 2007; Fazekas et al., 2008; Bruni et al., 2010; Wang et al., 2010; Burgess et al., 2011; Hollingsworth et al., 2011; Li et al., 2011; Zhang et al., 2012). While numerous studies demonstrate relatively high species delimitation success rates for traditional standard DNA barcodes with in some investigated plant genera and families (Wang et al., 2010; Li et al., 2011; Liu et al., 2011, 2017; Pang et al., 2011; Yu et al., 2011; Zuo et al., 2011; Chen et al., 2015; Xu et al., 2015; Yuan et al., 2015; Gogoi et al., 2020; Le et al., 2020), the delimitation rates markedly decline in plant lineages characterized by recent rapid radiation or suspected hybridization events (Roy et al., 2010; Yu et al., 2015, 2022; Fu et al., 2022; Zhang et al., 2023a). Notably, recent rapid radiation compromises species delimitation efficacy primarily through incomplete lineage sorting and insufficient genetic divergence due to short speciation intervals, compounded by morphological convergence under accelerated adaptive evolution (Robins et al., 2014; Xu et al., 2018; Yao et al., 2024). Moreover, hybridization between species is prevalent in plants, with at least 25% of plant species (Rieseberg, 1997; Mallet, 2005) – particularly among recently diverged lineages – demonstrating hybridization and introgression patterns with congeners. Therefore, it is challenging to discriminate all species using genetic differentiations.

The accelerated advancement of next-generation sequencing technologies has catalyzed the development of innovative strategies which are now available to overcome the limitations of traditional standard DNA barcodes (Nock et al., 2011; Li et al., 2015; Hollingsworth et al., 2016). The use of the complete plastid genome (plastome-scale barcodes) – also referred to as ultra-barcoding (Kane et al., 2012), organelle-scale barcodes (Yang et al., 2013), next-generation DNA barcoding (Fu et al., 2022) or super-barcodes (Li et al., 2015; Yu et al., 2022) – has been shown to improve species delimitation efficacy compared to traditional standard barcodes in plants. Implementation of the genome skimming protocols, complemented by optimized bioinformatic workflows and pipeline standardization, has enabled cost-effective deployment of these "super barcodes" approaches (Dierckxsens et al., 2017; Sarmashghi et al., 2019; Jin et al., 2020). Therefore, plastome-scale barcodes present significant potential to enhance species discrimination capacity for hyperdiverse and taxonomically complex genera. Nevertheless, accumulating empirical findings challenge the conventional paradigm of plastome-based phylogenomic efficacy in resolving species boundaries, with limitations attributable to pervasive biological phenomena including interspecific plastome introgression via recurrent hybridization, polyploid speciation, and ancestral incomplete lineage sorting maintained through rapid radiations (Fazekas et al., 2009; Hollingsworth et al., 2011, 2016; Yu et al., 2022; Lv et al., 2023; Su et al., 2023; Zhang et al., 2024a, b; Zhu et al., 2024).

Pedicularis L. (Lamiales: Orobanchaceae), one of the largest hemiparasitic angiosperm genera, consists of approximately 700 recognized species predominantly distributed across temperate and alpine ecosystems in the northern hemisphere (POWO, 2025; Wang and Yu, 2024). Notably, over 60% of Pedicularis species diversity is concentrated in the Himalaya-Hengduan Mountains region (Tsoong, 1963; Wang and Yu, 2024), which serves as the global species richness and diversification for this genus. Previous DNA barcoding study of 88 Pedicularis species from China documented that nrITS alone achieved 78.4% of discrimination rates across the genus, while combining nrITS with other plastid barcodes (matK, rbcL or trnH-psbA) increased discrimination rates to 81.8% (Yu et al., 2011). The highly diverse floral traits of Pedicularis can promote reproductive isolation among species (Eaton et al., 2012; Huang and Shi, 2013; Liang et al., 2018; Zhang et al., 2024a, b), and the uplift of the Himalaya-Hengduan Mountains creates geographical isolation, which accelerates species divergence and maintains high species richness and co-existence (Xing and Ree, 2017). However, pervasive biological phenomena such as incomplete lineage sorting, hybridization, and chloroplast genome capture events have been well-documented in the genus (e.g., Pedicularis Sect. Cyathophora in Yu et al., 2013; Pedicularis siphonantha complex in Liu et al., 2022), posing a challenge to species delineation by potentially generating discordance between morphological and molecular data. Therefore, the species discrimination rate using traditional standard DNA barcodes in Pedicularis reaches to 81.8%, which is deemed a high threshold in plants by excluding hybridization, introgression, incomplete lineage sorting, and other reticulated evolutionary processes (Mallet, 2005; Mallet et al., 2016). Given this, Pedicularis serves as an exemplary system for addressing the following three key questions: (1) Can plastome-scale barcodes enhance species discrimination capacity in Pedicularis compared to traditional standard DNA barcodes? (2) If not, what are key factors (sequence characteristics or evolutionary selection) influencing species discrimination capacity across the plastome sequences in this genus? (3) What biological factors limit the barcode discrimination rate when species cannot be discriminated?

2. Materials and methods 2.1. Taxon samplings

In this study, we analyzed 292 individuals representing 96 species of Pedicularis, with at least two samples per species. These samples cover 11 Groups and 48 Series under Tsoong's (1963) classification system (Table S1) and represent 12 out of 13 clades within the phylogenetic framework (Yu et al., 2015). DNA materials were from multiple sources, including silica gel dried fresh leave mainly collected from the Himalaya-Hengduan Mountains region in Southwest China, and supplemented by herbarium vouchers curated at the Herbarium of Kunming Institute of Botany (KUN), Chinese Academy of Sciences (CAS), and Institute of Botany (PE), CAS. Of the 292 samples used in this study, 213 were newly sequenced and assembled. The remaining 79 samples were obtained from concurrent studies, including 37 samples from Liu et al. (2024a) and 42 samples from W.-J. Wang et al. (unpublished data). For the outgroups, we downloaded the complete chloroplast genomes of five species from NCBI, i.e. Lindenbergia philippensis (Cham.) Benth. (HG530133.1), Phtheirospermum japonicum (Thunb.) Kanitz (MN075943.1), Rehmannia piasezkii Maxim. (KX636161.1), Triaenophora shennongjiaensis Xiao D. Li, Y.Y. Zan & J.Q. Li (MH071405.1), Wightia speciosissima (D. Don) Merr. (MK381318.1). Additionally, to validate our findings in the accuracy of species discrimination, we retrieved all available complete chloroplast genomes of Pedicularis species from NCBI, resulting in a dataset of 44 accessions representing 37 species (Table S5).

2.2. DNA extraction, sequencing, assembly and annotation

For the 213 new samples, the total genomic DNA of silica gel-dried leaves was extracted using a modified CTAB method (Doyle and Doyle, 1987), and the genomic DNA of herbarium materials were extracted using the DNAsecure Plant Kit (TIANGEN). The purified DNA was fragmented into about 350–500 bp for library construction using the NEBNext Ultra Ⅱ DNA Library Prep Kit for Illumina. Using the Illumina Hi-Seq 2500 platform, 150 bp of pair-end reads with more than 2 Gb raw data was generated for each sample.

Both plastome and nrITS sequences were assembled from the clean raw data using the GetOrganelle toolkit (Jin et al., 2020) with the following parameters "-w 0.65, –R 30 or 50 -k 21, 35, 55, 65, 85, 105". The complete plastome sequences were checked in Geneious v.8.1 (Kearse et al., 2012). Fragmented contigs were concatenated into complete plastome sequences in Geneious using the plastome sequence of Pedicularis tongolensis Franch. (MZ264887) as the reference. All plastomes were automatically annotated using PGA (Qu et al., 2019) and then manually adjusted in Geneious using the P. tongolensis plastome as the reference. Detailed information on plastome, samplings, vouchers, and classification for all samples is listed in Table S1.

2.3. DNA alignment and dataset preparation

Firstly, whole plastome sequences of all Pedicularis samples with an inverted repeat (IR) removed were aligned using MAFFT v.7.0 (Katoh and Standley, 2013). Then, we used DnaSP v.6.12 (Rozas et al., 2017) to assess highly variable regions and calculate nucleotide diversity (π) of the whole plastome matrix (WPM). The step and window sizes were set to 200 bp and 600 bp, respectively. Highly variable regions (HVR) were selected for evaluating barcodes. According to the plastome characteristics and previous studies of Pedicularis, seven datasets, including standard barcodes, super-barcodes and other potential combinations, were used to evaluate species discrimination: A: WPM (whole plastome matrix); B: LSC (the large single copy region of plastome); C: CDS (supermatrix of plastome coding DNA sequences); D: matK + rbcL + trnH-psbA (combination of three frequently-used plastid DNA barcode); E: HVR (matK + trnE-psbD + rpl20-clpP + ycf1, see details in Fig. S1); F: nrITS (internal transcribed spacer of nuclear ribosome DNA, consisting of ITS1-5.8S-ITS2); G: nrITS + matK + rbcL + trnH-psbA. In addition, we constructed a comprehensive plastome matrix H, which included the 292 samples used in this study and an additional 43 samples downloaded from NCBI. One of the 44 initially downloaded NCBI sequences was excluded following BLAST analysis, as it was determined to belong to the genus Dracocephalum L. (Lamiaceae). Each DNA region was aligned using MAFFT v.7.0 (Katoh and Standley, 2013). We then determined the aligned length (bp), number of variable sites (% divergence) and the number of parsimony-informative sites (% divergence) using MEGA X (Kumar et al., 2018) (Table 1).

2.4. Species delineation using genetic distance and phylogenetic tree approaches

We applied two widely used methods to evaluate species delimitation rates: distance-based and tree-based. For the distance-based method, we conducted analyses using the "Best match/Best close match" model in TaxonDNA (Meier et al., 2006). The best match model finds the closest match for each sample based on genetic distance. Here, the identification will be considered successful if the two samples have the same species name. The Best close match model finds a threshold value based on all intraspecific distances. Samples without a match below the threshold value will be considered unidentified. For the tree-based method, we constructed phylogenetic comparisons using the maximum likelihood (ML) approach in RAxML v.8.2.12 (Stamatakis, 2014). We use the GTRCAT model with 1000 bootstrap iterations. A species is considered a successful identification when all the samples of morphology-based species identification are resolved as a monophyletic group with a bootstrap support value ≥ 50%. We further recorded all such monophyletic species that achieved ≥ 70% bootstrap support.

2.5. Correlation analyses between sequence characteristics and species discrimination

To evaluate correlations between sequence characteristics and species discrimination, we conducted phylogenies of 66 plastid CDS genes, two pseudogenes, and three intergenic regions using the neighbor-joining algorithm in RapidNJ v.2.3.2 (Simonsen et al., 2008) and counted the number of identified species (NIS). We also assessed the following sequence characteristics, including alignment length (AL), number of variable sites (NVS), percentage of variable sites (PVS), number of parsimony-informative sites (NPIS) and percentage of parsimony-informative sites (PPIS) for the 71 genomic loci. Using the "cor" function in the R package stats, we conducted correlation analyses and visualized the correlation coefficient using the R package corrplot. Then, taking NIS as the response variable and the remaining as predictor variables, we conducted an All-Subsets Regression (ASR) analysis to get the best model using the "regsubsets" function in R package leaps.

2.6. Selective pressure analysis

To explore whether genetic factors may influence species discrimination, we conducted selective pressure analysis for 66 CDS genes. We used MACSE v.2 to align each CDS gene with Phtheirospermum japonicum (Thunb.) Kanitz as the reference. Then, the "Simple Ka/Ks Calculator" function of TBtools (Chen et al., 2023) was used to calculate the ratio of non-synonymous substitutions (Ka) to synonymous substitutions (Ks). Ka/Ks > 1 indicates positive selection, Ka/Ks < 1 indicates purifying selection and Ka/Ks = 1 indicates neutral selection. If Ka > 0 and Ks = 0, the value of Ka/Ks will be represented by NA.

We also used the HyPhy program (Pond et al., 2005) with the hypothesis testing framework RELAX (Wertheim et al., 2015) on Datamonkey (https://www.datamonkey.org/) to test selection strength (relaxation or intensification). We set five species (Lindenbergia philippensis, Phtheirospermum japonicum, Rehmannia piasezkii, Triaenophora shennongjiaensis, Wightia speciosissima) as reference branches and all samples of Pedicularis as test branches. K < 1 indicates relaxation of selection, K = 1 indicates neutral selection, and K > 1 indicates intensification of selection. Moreover, we examined potential variations in discriminatory power among genes experiencing differential selection intensities (relaxation vs. intensification).

3. Results 3.1. Dataset characteristics

The plastid genomes of the 213 samples produced in this study allowed for the generation of a complete circular genome structure. Of the 292 plastomes of Pedicularis, the genome size ranged from 140, 174bp to 162, 218bp, including 74–97 protein-coding genes, eight rRNA genes and 36–39 tRNA genes (Table S1). The LSC region size ranged from 78, 698bp to 83, 843bp, where the IR region varied in length, and the SSC region was absent in a few samples (e.g., P. chenocephala W2021138 and P. lutescens LR202070803). The GC content of Pedicularis species ranged from 37.8% to 38.6% (Table S1).

Of the seven datasets, the WPM dataset (A) had the greatest number of variable sites and parsimony-informative (PI) sites (53, 074 bp (31.62%) variable sites and 41, 762 bp (24.88%) PI sites). Still, the nrITS dataset (F) had the highest percentage of variable sites and PI sites (385 bp (58.69%) variable sites and 342 bp (52.13%) PI sites) (Table 1). For the combination of plastid genes and regions, the HVR dataset (E) had a higher percentage of variable sites and PI sites (5277 bp (41.98%) variable sites and 4191 bp (32.72%) PI sites), followed by the CDS dataset (C) had variable sites and PI sites (15, 411 bp (26.44%) variable sites, and 12, 265 bp (21.04%) PI sites). The combination of three traditional plastid barcodes (dataset D: matK + rbcL + trnH-psbA) (1571 bp (38.06%) variable sits, and 1259 bp (30.50%) PI sites), which is more than the LSC dataset (B) in the percentage of PI sites (36, 739 bp (32.26%) variable sites, and 28, 936 bp (25.41%) PI sites) (Table 1).

In addition, as the largest hemiparasitic plant genus, 14 genes have undergone varying degrees of gene loss and pseudogenization in Pedicularis (Orobanchaceae), including eleven NDH (NA(D)H dehydrogenase) genes associated with photosynthesis, ccsA, accD and ycf15. Based on the degree of pseudogenization and loss, the state of each gene is categorized into four stages: (0) full length with function, (1) full length with a premature stop codon, (2) partial loss, and (3) complete loss (Fig. 2).

3.2. Species delineation 3.2.1. Tree-based method

Ninety-six species could be attributed to 12 clades using seven datasets. Of them, six datasets (A, B, C, D, E, and G) (Figs. S2–S6 and S8) which is relatively consistent with the topology of Yu et al. (2015). The remaining dataset (F) produced different topologies with low support values among different clades (Fig. S7). Moreover, all datasets have high support values for species monophyly.

Based on the ML method, the plastome-scale barcodes (datasets A, B and C) showed comparable or slightly higher discrimination rates compared to different combinations of standard barcodes (datasets D and G) and HVR datasets. Among the seven datasets, the dataset G exhibited the highest discriminatory rate, distinguishing 78 species (81.25% of total species, hereafter), followed by the dataset B (77 species, or 80.21%), the datasets A, C and D (76 species or 79.17%). The datasets D discriminated 74 species (77.08%), while the dataset F discriminated 60 species (62.50%) (Table 2). Applying a ≥ 70% bootstrap threshold, discrimination failed for just three species in the traditional barcode datasets (D, F and G), a single species in the LSC and CDS datasets, and remained unchanged in the WPM and HVR datasets (Table 2 and Fig. 3).

For the dataset H, which includes 43 samples downloaded from NCBI, 23 samples represent 20 species, which have been sampled in our study and show consistency with our sampled samples. Of these 23 samples, three samples (marked in red in Fig. S9) should be misidentified or mislabeled (error rate: 13.04%), and four samples (marked in blue in Fig. S9) had incorrect phylogenetic placements in accordance with the phylogenetic study of Yu et al. (2015). Overall, eight of the 44 initially downloaded NCBI sequences (error rate: 18.18%) showed clear identification errors. This minimum observed error rate of 18% underscores the necessity for rigorous species identification within the genus Pedicularis and highlights the urgent need to establish a comprehensive, verified reference dataset.

3.2.2. Distance-based method

In the distance-based method, all datasets showed a high discrimination rate with relatively small differences. For the best match model, the result showed that dataset G had the highest discrimination rate (89.38%), followed by dataset C (88.01%), and dataset E had the lowest discrimination rate (79.10%) (Table 2). The best close match model exhibited a trend similar to the best match model, with the sole distinction being that the dataset G achieved the highest discrimination rate (88.69%). The coding sequences dataset C ranked second with a discrimination rate of 88.01%, followed by LSC (85.95%), WPM (84.24%), nrITS (81.50%), matK + rbcL + trnH-psbA (80.82%) and HVR (79.10%), respectively (Table 2).

3.3. Correlation and selective pressure analyses

We assessed species discrimination efficacy (NIS) by analyzing 71 genomic loci (including 66 CDS genes, two pseudogenes, and three intergenic regions). The matK and ycf1 genes exhibited the highest discriminatory power, distinguishing 75 species (78.13%), whereas the petN gene showed limited utility, discriminating only one species (1.04%) (Tables S2 and S3). Correlation analyses revealed strong positive associations between AL and both NVS (r = 0.88) and NPIS (r = 0.87), but weak correlations with proportional metrics (PVS: r = 0.14; PPIS: r = 0.13). Species discrimination capacity (NIS) showed moderate correlations with all five variables analyzed (r = 0.47–0.67). Notably, PVS and PPIS showed strong linear relationship (r = 0.98), while NVS and NPIS were perfect correlated (r = 1.00). To further disentangle these interdependencies, the ASR analysis selected the composite model integrating AL and PPIS as optimal, indicating their synergistic influence on species discrimination efficacy (NIS) (Fig. 4).

The selective pressure analyses revealed that all genes were subject to purifying selection in nearly all examined species, with the exception of atpF, infA, rps18, petL, and rpl23 genes in some species (Fig. 5 and Table S6). Variable amplitude of Ka/Ks values increased with a concomitant decline in species discrimination (Fig. 5). The RELAX analyses revealed that 36 genes underwent intensified selection pressure, while 30 genes experienced relaxed selection (Fig. 6). Genes under intensified selection showed lower efficacy of species discrimination with an average of 29.17 identified species (median = 27.50). In contrast, genes subject to relaxed selection exhibited higher efficacy of species discrimination, with an average of 40.43 identified species (median = 46.50) (Fig. 6). A statistically significant difference was observed between these two gene groups (Mann–Whitney U test, P = 0.041).

4. Discussion 4.1. Parsimony-informative sites determining the power of species discrimination

The phylogenetic tree-based method for species discrimination in Pedicularis revealed that the supermatrix combining nrITS and three traditional standard barcodes (dataset G) achieved the highest discrimination rates at 81.25%, outperforming the WPM dataset (dataset A) and the LSC dataset (dataset B), by discriminating two and one additional species, respectively. Furthermore, the matK dataset and the concatenated matrix of three standard plastid barcodes (dataset D: matK + rbcL + trnH-psbA) successfully discriminated 75 (78.13%) and 74 species (77.08%), respectively. Compared to standard DNA barcodes, previous studies have shown that the plastome-scale barcodes generally exhibit higher discrimination rates. For example, the plastome barcodes have improved species discrimination from 33% to 55% in Rhododendron L. (Ericaceae) (Fu et al., 2022); 58%–68% in Cymbidium Sw. (Orchidaceae) (Zhang et al., 2023a); 0%–27% in Schima Reinw. (Theaceae) (Yu et al., 2022), and 56%–78% in Calligonum L. (Polygonaceae) (Song et al., 2020). However, our findings present a notable exception that there are minimal differences in species discrimination rates between traditional DNA barcodes and plastome-scale barcodes. This discrepancy appears to be linked to the specifically biological habits of Pedicularis and may reflect unique patterns of their plastid genome.

Pedicularis, the largest hemiparasitic and herbaceous genus within the family Orobanchaceae (Lamiales) (Li, 1951), has undergone extensive pseudogenization and gene loss (Zhang et al., 2020; Li et al., 2021). These genomic changes have led to instability in plastid genome structure, particularly evident in the expansion and contraction of the inverted repeat (IR) and the small single copy (SSC) regions (Zhang et al., 2020; Li et al., 2021). Therefore, the evolutionary stability of the large single copy (LSC) region (dataset B) has more PVS and PPIS, so that it can discriminate one more species (77 of 96 species) than the complete plastid genome (dataset A: 76 of 96 species) in Pedicularis. However, what is particularly noteworthy is that the supermatrix datasets combining standard DNA barcodes (i.e., datasets D and G) exhibited discriminatory power comparable to the plastome-scale barcodes (i.e., datasets A, B, and C), and the single matK gene or the ycf1 gene also showed high discriminatory power (75 of 96 species). These findings indicate that while simply increasing the number of variable sites does not necessarily enhance species discrimination in Pedicularis. Instead, the barcodes with more PIS and PPIS serve as critical determinants for successful species discrimination. Four standard barcodes were originally proposed as universal plant barcodes due to their sufficient variation rates and PCR amplification and sequencing feasibility, and high performance in species discrimination (Kress et al., 2005; CBOL Plant Working Group, 2009; Li et al., 2011). However, these genes exhibit dramatically elevated evolutionary rates in parasitic plants compared to their non-parasitic relatives (Nickrent and Starr, 1994; Bromham et al., 2013; Wicke et al., 2016), particularly in the plastid genes, leading to the accumulation of more variable sites. The parasitic lifestyle of Pedicularis likely explains the high discrimination capacity of standard DNA barcodes derived from plastid genes. Our findings further revealed that the integration of AL and PPIS for predicting discriminatory power, combined with the unexpected strength of conserved genes under relaxed selection, provides a transferable framework beyond Pedicularis. To apply this to other groups: (1) Pre-screen candidate loci to calculate AL and PPIS; prioritize loci within the target zone. (2) Identify key genes: Perform selection pressure analysis (dN/dS) on pre-screened loci, focusing on traditionally conserved genes exhibiting relaxed selection. This framework efficiently directs resources to the most promising markers.

Phylogenetic relationships within Pedicularis have been significantly clarified through plastome phylogenies in comparison with previous phylogenies using several DNA regions (Yang et al., 2003; Ree, 2005; Yang and Wang, 2007; Tkach et al., 2014; Robart et al., 2015; Yu et al., 2015). The plastome phylogeny revealed shorter internodes along main branches, contrasted by extended terminal branch lengths, which suggests an early rapid radiation event in this genus. Unlike other species-rich groups that undergo extensive interspecific hybridization or introgression (Fu et al., 2022; Zhang et al., 2023b; Liu et al., 2024b; Wang et al., 2024) and show a relatively chaotic state of phylogeny, the monophyly of most Pedicularis species have received relatively strong support. However, given the widespread occurrence of natural hybridization across plant lineages and the fact that approximately 25% of angiosperm species are involved in hybridization and introgression (Mallet, 2005), the species discrimination accuracy up to 80.21% achieved in this study (based on plastid dataset) likely approaches the theoretical maximum resolution attainable through genetic markers. The residual unresolved taxa (19.79%) may reflect both inherent limitations of plastid data and potential inconsistencies in current taxonomic delineations.

4.2. Key constraints on species discrimination using DNA barcodes: interspecific chloroplast genome capture and ancestral incomplete lineage sorting

Pedicularis is a highly diversified group (Fig. 1) that has recently undergone rapid diversification in the Himalaya-Hengduan Mountains region, driven by orogenic uplifts and monsoon dynamics (Xing and Ree, 2017; Ding et al., 2020). Although reproductive interference could be effectively reduced by divergent floral trait assembly among sympatrically distributed Pedicularis species (Eaton et al., 2012), this may lead to low level of interspecific gene flow and hybridization among co-flowering and sympatric species. Field experiments in pollination biology demonstrated that mechanical isolation mediated by floral trait divergence are imperfect in Pedicularis (Huang and Shi, 2013; Liang et al., 2018), while post-pollination isolation can play a critical role among coexisting Pedicularis species with pollen exchange (Liang et al., 2018). Nevertheless, ancient hybridization/introgression events have been well documented in Pedicularis sect. Cythophyalla (Yu et al., 2013; Li et al., 2021), and the P. siphonantha lineage (Yu et al., 2018; Liu et al., 2022). In the section Cyathophylla, cytoplasmic-nuclear genome discordance manifests through chloroplast capture events - Sichuan populations of P. cyathophylloides retain nrITS monophyly while bearing captured plastomes from an ancestral P. cyathophylla lineage (Yu et al., 2013). Such chloroplast genome capture events likely extend to some unidentified species in this research, such as the P. giraldiana clade, the P. rhynchodonta clade, and the P. verticillata clade (Yu et al., 2015). The P. rex complex exemplifies even greater phylogenetic complexity: subsp. rex displays plastid paraphyly due to incomplete lineage sorting (ILS) during population expansions, whereas var. rockii exhibits polyphyly through ancient hybridization (Yu et al., 2013). These dual mechanisms (incomplete lineage sorting vs. chloroplast genome capture) likely explain the mosaic plastid phylogeny of P. rex observed in our analyses.

Notably, our data revealed a distinct geographical signature in hybridization/introgression patterns, where population clustering reflects spatial proximity rather than taxonomic affinities. This phenomenon is particularly evident in the phylogenetic discordance observed across multiple species (Mallet et al., 2016). For example, both plastid phylogenies and nrITS phylogeny separate five samples of P. anas into two groups, three northern Sichuan samples and two southern Sichuan samples (Figs. S2–S8). Three northern Sichuan samples clustered with P. cheilanthifolia and P. giraldiana, and two southern Sichuan samples clustered with P. rupicola. Similar patterns are observed in four samples of P. roylei (two from Sichuan, and two from Xizang), and four samples of P. cheilanthifolia (two samples from Qinghai, and two samples from Xinjiang) (Figs. S2–S8), demonstrating consistent biogeographic clusters with potential hybridization/introgression events.

4.3. DNA barcoding challenges cryptic and recently diverged species

The taxonomic delimitation of recently diverged species is often complicated by their highly similar morphological traits, which can lead to significant challenges in accurate identification. Taxonomically, diagnostic traits of Pedicularis are mainly based on vegetative and reproductive characters, including root morphology, leaf architecture, phyllotaxy, inflorescence structure, floral tube dimensions, galea curvature, and calyx features. However, some key diagnostic traits (e.g., the degree of galea curvature) show minimal variation among closely related species, which might be difficult to quantify using specific terms and to distinguish from close relatives. Consequently, species discrimination based on plastid DNA markers or other neutral DNA regions tend to be failed or unresolved (Ross, 2014; Yu et al., 2022), especially considering the presence of cryptic species in Pedicularis (Yu et al., 2013, 2018; Wang et al., 2015; Liu et al., 2022). For example, a new species from the Jiaozi Mountain should be separated from P. miliana through integrative morphological, molecular, and biogeographic evidence in the P. siphonantha complex (Liu et al., 2022). In particular, herbarium specimens have lost key diagnostic traits for species identification, which could cause taxonomic confusion among this lineage (Yu et al., 2018; Liu et al., 2022). Furthermore, our analyses revealed that ten samples of P. davidii form a monophyletic clade based on the plastid datasets. In contrast, these samples split into two separated clades in the nrDNA phylogeny. Field observations provided additional insights, showing that two clades correspond to high-altitude populations and low-altitude populations, respectively. This pattern highlights a striking morphological differentiation between high- and low-altitude populations, suggesting a cryptic species in the P. davidii lineage by a recently phylogenetic divergence and ecological adaptation. These findings highlight the critical need for refined taxonomic frameworks in Pedicularis to reconcile plastid-based species discrimination with evolutionary realities.

Hybridization dynamics further complicate species boundaries, as evidenced by field documentation of introgressive populations exhibiting intermediate morphologies. For example, multiple hybrid populations exist between Pedicularis milliana and P. sigmoidea in the Yulong Snow Mountain region. High-altitude hybrids backcrosse with the ancestral lineage of P. milliana by inheriting their parental morphological traits, while low-altitude hybrids are in the process of speciation and display morphological differences from P. milliana and P. sigmoidea (Liu et al., 2022). Undoubtedly, the stochastic inheritance and independent evolution of the parental traits in hybrid populations complicate species identification. Although this phenomenon is common in highly diversified groups, which may also harbor undescribed species (Struck et al., 2018), the efficacy of plastid DNA barcodes in discriminating cryptic Pedicularis species that lack visually distinguishable morphological differences appears limited. The development of nuclear markers may provide a promising avenue for disentangling complex evolutionary history within the genus.

5. Conclusions

This study highlights the efficacy of standard DNA barcodes (e.g., nrITS, matK, rbcL, and trnH-psbA) in species discrimination within the hemiparasitic genus Pedicularis, demonstrating performance comparable to plastome-scale barcodes. Key factors influencing species resolution included alignment length (AL) and the percentage of parsimony-informative sites (PPIS), with conserved genes under relaxed selection exhibiting higher discriminatory power. The study also revealed challenges of species discrimination posed by rapid diversification, hybridization, and incomplete lineage sorting in Pedicularis. Phylogenetic discordance among chloroplast and nuclear markers (e.g., nrITS) underscored the impact of historical introgression and geographic isolation, particularly in recently diverged lineages. Morphological convergence and cryptic speciation further complicated taxonomic delineation, as evidenced by divergent molecular clades within morphologically uniform species like P. davidii. These findings suggest that plastid DNA barcodes alone may approach their theoretical resolution limits (~80%) in Pedicularis, with residual uncertainties likely reflecting taxonomic inconsistencies or limitations of plastid data. Future studies should integrate nuclear markers to resolve hybridization-driven complexities and validate cryptic species. Enhanced taxonomic scrutiny, combined with genomic approaches, will be critical for refining species boundaries in this rapidly radiating genus.

Acknowledgements

This study was supported by the National Natural Science Foundation of China (32371700, 32071670 and 31870196), the Strategic Priority Research Program of the Chinese Academy of Sciences (XDB31000000), the Science and Technology Basic Resources Investigation Program of China (2021FY100200), Yunnan Revitalization Talent Support Program "Young Talent" and "Innovation Team" Projects (202405AS350019), the 14th Five-Year Plan of Xishuangbanna Tropical Botanical Garden, Chinese Academy of Science (XTBG-1450101), the Key R & D program of Yunnan Province, China (202103AC100003) and the Key Basic Research program of Yunnan Province, China (202101BC070003). We are grateful to Jie Cai, Li-Na Dong, Lian-Ming Gao, Hua-Jie He, Jun He, Wei Jiang, Rong Li, Bin Liu, Cheng Liu, En-De Liu, Jie Liu, Min-Lu Liu, Lu Lu, Yang Luo, Hui Tang, Chun-Lei Xiang, Ji-Dong Ya, Qiu-Lin Yang, Xiu-Long Yang, Ting Zhang, and Shu-Dong Zhang for their help with fieldwork and providing plant samples, and to Jing Yang and Zhi-Rong Zhang for their help and suggestions in the lab work, and to the physical support from the Germplasm Bank of Wild Species, Kunming Institute of Botany, Chinese Academy of Sciences, and the Information Center, Xishuangbanna Tropical Botanical Garden, Chinese Academy of Sciences.

Data availability statement

All DNA sequences generated in this article have been uploaded to China National Center for Bioinformation (CNCB) and the GenBank accession numbers can be found in the supplementary material of this article.

CRediT authorship contribution statement

You Wu: Investigation, Data curation, Formal analysis, Visualization, Writing original draft, Writing − review & editing. Rong Liu: Resources, Methodology, Investigation. Wei-Jia Wang: Resources, Investigation. De-Zhu Li: Resources, Methodology, Investigation. Kevin S Burgess: Writing − review & editing, Supervision. Wen-Bin Yu: Methodology, Investigation, Visualization, Writing − review & editing, Supervision, Funding acquisition. Hong Wang: Investigation, Writing − review & editing, Supervision, Funding acquisition.

Declaration of competing interest

The authors Dezhu Li and Hong Wang are Editors for Plant Diversity and were not involved in the editorial review or the decision to publish this article. The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.pld.2025.09.005.