Phylogenomic analysis of Bupleurum in Western Sichuan, China, including an overlooked new species

A comparative analysis of chloroplast (cp) genomes and 45s nuclear ribosomal DNA (nrDNA), and a phylogenomic study of six closely related species (including an overlooked new species) of genus Bupleurum from the western part of Sichuan Province in southwestern China were performed. The six species are similar morphologically and it is difficult to identify them; moreover, their genetic relationships remain unclear. It was found that the cp genomes of the six Bupleurum species were extremely similar, and they were highly homogeneous in terms of cp genome structure, genes and its arrangement. Intergenic spacer rpl32-trnL, petA-psbJ, trnK-rps16, and the coding gene ycf1 were considered highly variable. In phylogenetic trees constructed based on the complete cp genome, protein-coding sequences, nrDNA and ITS sequences, Chinese Bupleurum species all formed two major clades; among these trees, nrDNA tree had the best species resolution; the highly variable regions showed no advantage over other molecular markers. Among the six Bupleurum species, B. malconense, B. sichuanense were close relatives to B. chinense and B. yinchowense, B. chaishoui may also be a consanguinity, while B. microcephalum, B. wenchuanense, and the new species B. pseudochaishoui were closely related. At the end, the new species B. pseudochaishoui Z. Chao sp. nov. was described and illustrated, and a key to the six species was tabulated.


Introduction
Bupleurum L., comprising about 210 species accepted nowadays (POWO, 2023), is one of the largest genera of the family Apiaceae.It has long been recognized as a natural group distinctively characterized by its members having simple leaves, mostly parallel veins, and conspicuous bracts/bracteoles.Infrageneric classification and phylogeny of the genus was of interest to plant taxonomists from the early days.The first use of section names in Bupleurum can be dated back to 1848 (Grenier and Godron, 1848).Wolff published a comprehensive revision of the genus that classified it into 5 sections (Wolff, 1910), and his classification is the most widely used for about a century.Shu et al. carried out a numerical taxonomic study in 1998 and divided the Chinese Bupleurum species into subgenus Longifolia and subgenus Eubupleura (including section Ranunculoidea and section Falcata) (Shu et al., 1998).Neves & Watson suggested dividing Bupleurum into two subgenera (subg.Penninervia, subg.Bupleurum) based on phylogenetic analyses of nuclear ribosomal internal transcribed spacer (nrITS) sequences of 35 taxa (Neves and Watson, 2004).Although the previous researches have deepened our understanding of the genus, it is still hard to say being sufficient.
China harbors abundant diversity of Bupleurum.In Flora of China published in 2005, 42 species, 16 varieties and 6 forms were recorded (She and Watson, 2005).After revisions of He et al. and Pimenov, etc. in recent years, it was acknowledged there are 47 species, 1 subspecies, 6 varieties and 5 forms of Bupleurum in China (Wang et al., 2011;Ma et al., 2013;Wang et al., 2013;Ma and He, 2015;Pimenov, 2017), most of which have reputed medicinal values in treating fever and liver ails (Liang et al., 2012;Chao et al., 2014), as observed in many wild species of Apiaceae (Moazzami Farida et al., 2018;Sonigra and Meena, 2021;Perrino et al, 2023).
In the western part of Sichuan Province in Southwestern China, i.e., some areas of Aba Zang and Qiang Autonomous Prefecture near Chengdu, including Wenchuan, Mao County, Songpan, and so on, the steep mountains and plateaus lead to the differentiation and formation of Bupleurum species (Wang et al., 2008a;Yang et al., 2021).Closely related species discovered and reported in this region include  (Pan and Hsu, 1992).Most of them are quite morphologically similar and the differences between each other are subtle, thus their identification is notoriously difficult.And also, their genetic relationship remains unclear.
Phylogenomics, especially that based on the whole chloroplast (cp) genome sequences, has been proved highly efficient in resolving complicated phylogenetic problems of taxonomically complexed plant groups (Liu et al., 2022;Zhang et al., 2023), and in differentiating and identifying closely related species (Pchelin et al., 2019;Li et al., 2021).We've made some efforts in elucidating phylogeny of Chinese Bupleurum through phylogenomic means.We sequenced cp genomes of two woody Bupleurum species (B.gibraltaricum Lam. and B. fruticosum L.) endemic to Mediterranean and made comparative analysis with some Chinese species, confirming the division of 2 subgenera (subg.Penninervia, subg.Bupleurum) and revealing their divergence time (Huang et al., 2021a); we also analyzed the sequences of whole cp genomes of 4 Southwest China endemic alpine Bupleurum species and suggested phylogenetic affinity in the genus might be closely related to geographical distribution (Huang et al., 2021b); and we proposed the phylogenetic position of B. sikangense C.B. Wang, X.G.Ma & X.J.He as well (Xie et al., 2021).As a continuation of the work, we made a phylogenomic analysis of the above-mentioned species distributed in western Sichuan, based on both cp genome and 45s nuclear ribosomal DNA (nrDNA, including 16s rRNA-ITS1-5.8srRNA-ITS2-26s rRNA), in order to untangle the complicated relationship among them, and to seek suitable molecular markers for species identification.
Additionally, in the ongoing study on Bupleurum investigation, the authors encountered a population in which the individuals obviously belonged to this genus in Wenchuan county years ago.For all the time, the authors had been regarding the plants in this population as B. chaishoui.Until recently, the authors obtained some specimens of B. chaishoui in its type locality and had a deeper understanding of the species.The authors realized that the plants met before is unique in its leathery dimorphic leaves and represent an overlooked new species, and then named it as B. pseudochaishoui.It was thus described here.

Plant materials
In this study, 65 plant samples, representing 29 Bupleurum species were collected in the summer of the years of 2014 through 2020, from North, Northeastern, Northwestern and Southwestern China, including 20 samples of the six above-mentioned species in Western Sichuan (Table 1).The voucher specimens were deposited in the herbarium of the University (Southern Medical University herbarium, SMU, to be listed in the Index Herbarium).
Specimens of other Chinese Bupleurum species, especially those distributed in Western Sichuan, which were deposited in the major herbaria of China including China National Herbarium (PE), Herbarium of Department of Biology, Sichuan University (SZ), Herbarium of Chengdu Institute of Biology, Chinese Academy of Sciences (CDBI), Kunming Institute of Botany Herbarium (KUN) and Herbarium of College of Life Sciences, Northwest Agriculture & Forestry University (WUK), were also consulted and inspected on site or through Chinese Virtual Herbarium (CVH) website.

Morphological inspection
The morphological characters of taxonomic significance in Bupleurum (Pan et al., 2002), i.e., the number, shape and size of bracteoles, comparison of bracteoles length and umbellule size, the number and length of rays, the number and size of florets, the shape, size and texture of leaves, of B. pseudochaishoui were inspected, and were compared with the related species, B. wenchuanense, B. chaishoui, B. malconense, B. sichuanense, and B. microcephalum.The data were acquired from at least 10 individuals of living materials or herbarium specimens, and were also checked in accordance with the descriptions in Flora of China (She and Watson, 2005).For measurements of some quantitative traits, such as plant height, leaf size, the number of small flowers, etc., 50 pieces of data were collected for each item, then the variation ranges were summarized and recorded.

DNA extraction
Genomic DNA was isolated from the silica gel dried leaves of each sample using DNAsecure Plant Kit (TIANGEN Biotech Co, Ltd., Beijing, China).After purification, concentration of the extracted DNAs was determined on a NanoDrop 2000C spectrophotometry (Thermo Fisher, US).

ITS region amplification and sequencing
ITS region of some samples (those with ITS accession number marked by "*" in Table 1) were amplified and sequenced.The universal primers ITS-p4 and ITS-p5 (Cheng et al., 2015) were used.PCR amplification was performed in a 15 mL reaction system containing 6 mL of 2×Taq Plus PCR MasterMix (Tiangen Biotech Co, Ltd., Beijing, China), 2 mL of DNA template, 0.5 mL of each primer (10 mM), and 6 mL of ddH 2 O.The PCR products were purified and sequenced bidirectionally by Ruibo Biotech Co. Ltd. (Guangzhou, China).The sequences have been deposited in GenBank, and the accession numbers were listed in Table 1.We also downloaded several ITS sequences of the related species from GenBank (Table 1).

Skimmed genome sequencing
The purified genomic DNAs were subjected to genome skimming performed on an BGISEQ-500 sequencer by BGI Genomics (Shenzhen, China).For each sample, a total of 3.0 Gb sequence data was obtained.

Whole cp genome and nrDNA assembly and annotation
Whole cp genomes of the samples were assembled using GetOrganelle v1.7.4 software (Jin et al., 2020).The published cp genome of B. chinense DC. (GenBank accession number: MN337347) was used as the reference genome (Zhang et al., 2019).Bowtie2 (Langmead and Salzberg, 2012) was first called to filter out plastid-like reads, and then SPAdes (Bankevich et al., 2012) was called for assembly.The cp genome was annotated using GeSeq (Tillich et al., 2017) and Plastid Genome Annotator (PGA) software (Qu et al., 2019), and the tRNAs in the genome were identified by tRNAscan-SE v2.0.7 (Lowe and Chan, 2016) and ARAGORN v1.2.38 (Laslett and Canbäck, 2004).The annotated genomes, together with the reference genome, were input into Geneious v22.2.2 (Kearse et al., 2012) for comparison to determine start and stop codons, as well as to manually correct the possible errors.Finally, the cp genome was physically mapped into a loop using the online software OGDRAW (Lohse et al., 2007).The assembled and annotated cp genome sequences were deposited in GenBank (http://www.ncbi.nlm.nih.gov/)under accession numbers OQ460227-OQ460239 (Table 1).nrDNAs were also assembled using GerOrganelle v.1.7.3 software, and the output contigs were compared and annotated with the reference genome in Geneious v22.2.2.

Genome structure and comparative analysis
The cp genome characteristics of the six species of Bupleurum were described, including the overall genome structure, the size of each part of the genome, and the GC content of the genes, which was acquired using Geneious v22.2.2.Multiple genome alignments of 13 sequences were done using the Mauve v2.4.0 plugin (Darling, 2004) in Geneious v22.2.2, with default parameters.
To determine the expansion/contraction of IR, genes distributed within and next to the boundaries of LSC, SSC and IR regions were compared and manually mapped to visualize them.
To investigate the highly variable regions (divergence hotspots), 13 complete chloroplast sequences of six Bupleurum species, were compared using MAFFTv.1.3.7 (Katoh and Standley, 2013).Subsequently, the cp genomes of Bupleurum species were paired and visualized using the Shuffle LAGAN mode in mVISTA (Frazer et al., 2004), with the cp genome of B. chinense (GenBank access number: MN337347) as a reference.Using DnaSP v6.0 (Rozas et al., 2017), the nucleotide diversity (pi) in the cp genome sequence was calculated through sliding window analysis, with a window size of 600 bp and a step size of 200 bp.

Molecular phylogenetic analysis
In order to get insight into the phylogenetic relationship among the Bupleurum species, maximum likelihood (ML) analyses and Bayesian inference (BI) analyses based on the whole cp genome, the coding sequences (CDS), the highly variable regions, the nrDNA, and the ITS region were performed respectively.Besides the sequences obtained in this study, the previously reported cp genome and ITS sequences of Chinese Bupleurum species were also downloaded from Genbank and incorporated into phylogenetic analysis (Table 1).The best model was obtained under AIC criteria using ModelFinder (Kalyaanamoorthy et al., 2017) for molecular modeling of the dataset (Table 2).For ML analyses, IQ-tree 1.6.12(Nguyen et al., 2015) was used, and the reliability of each branch was assessed using bootstrap values (number of replicates set to 1000).As for BI analyses, MrBayes 3.2.7a(Ronquist et al., 2012) was used, and Markov chain Monte Carlo (MCMC) analysis was run for a total of 20 million generations, with samples taken every 1000 generations.In the phylogenetic tree construction, Centella asiatica were used as outgroup species, and the resulting tree files were visualized and embellished using iTOL v.6.1.2software (Letunic and Bork, 2021).
Genes and its arrangement in the cp genomes of the six Bupleurum species are almost the same.Each cp genome encodes 114 genes, including 80 protein-coding genes (PCG), 4 rRNA genes, and 30 tRNA genes.There're 60 PCGs and 22 tRNA genes in the LSC region, 12 PCGs and 1 tRNA gene in the SSC region.In the IR region, for 12 of the 13 cp genomes, there're 20 gene repeats, including 8 PCGs, 7 tRNA genes and 4 rRNA genes, except for one B. wenchuanense sample (voucher 2020080639), which has 2 extra PCGs (rpl22 and rps3) and thus a total of 22 genes in the IR region.Among these 114 genes, eighteen genes are found with introns, 15 of which containing one intron, while the other three (clpP, ycf3, and rps12) containing two introns.Incomplete copies of ycf1 and rps19 in the IR region are considered as pseudogenes.(Table 3, Table 4, Figure 1).

Comparative analysis of chloroplast genome and nrDNA structure
Results of Mauve pairwise analysis indicate that no rearrangement has occurred in coding and noncoding regions of the cp genomes of the six Bupleurum species (Figure 2), implying that the cp genomes are conservative.
The genes rpl22 and rpl2 are located on the left and right sides of the LSC/IRb junction respectively, while the gene rps19 crosses the JLB line (70 bp occurring in the IR region).The ndhF gene is at the right side of the IRb/SSC junction and 14-26 bp away from the boundary.Traversing the SSC and IRa region is the gene ycf1; a fragment of 1877-1902 bp in length of the gene is located in the IR region.The gene trnH lies 4-bp away from the JLA line (IRa/LSC boundary) to the right (Figure 3).
There's an obvious expansion in the IR region of a B. wenchuanense sample (voucher 2020080639).In this sample, as mentioned above, genes rpl22 and rps3 are found in the IR region.Generally, these two genes were located in the LSC region.The rps3 gene is crossing the junction of LSC/IRb (JLB) and LSC/IRa (JLA).(Figure 3).We compared and mapped whole chloroplast gene alignments among 13 sequences of six Bupleurum species using mVISTA, using the published cp genome of B. chinense (No.MN337347) as the reference.The results showed that the divergence in the noncoding regions was greater than that in the coding regions, and a high degree of divergence in the 13 cp genomes appeared at intergenic spacers, including trnK-rps16, atpF-atpH, atpH-atpI, trnC-petN, petN-petM, ndhC-trnV, petA-psbJ, accD-psaI, rps15-ycf1, rpl32-trnL, and trnV-rps7.The coding regions were more conserved, and trnL, rpoC, petD, ycf2, ndhA and ycf1 were the more divergent coding regions of these species.(Figure 4).
In order to assess the sequence divergence level of these Bupleurum species, the nucleotide diversity (Pi) of the cp genomes and nrDNAs were calculated.
The Pi values of the whole cp genome of the six species range from 0.00-0.06,and the two IR regions were significantly less divergent than the LSC and SSC regions.The spacer regions between genes were found with higher Pi values than the coding regions.Totally there're 7 spacer regions with Pi values greater than 0.01; the spacer with the highest Pi value was rpl32-trnL (Pi> 0.06), followed by petA-psbJ (Pi> 0.02), and consecutively, trnK-rps16, atpF-atpH, atpH-atpI, ndhC-trnV and rps15-ycf1 (Pi> 0.01).Only two coding regions (trnL-UAA and ycf1) had Pi values higher than 0.01, and the ycf1 gene has higher value.(Figure 5A).
The concatenated nrDNA sequences of the six Bupleurum species have 65 variation sites and 54 parsimony-informative sites, accounting for 1.12% and 0.93% of the total sequences, respectively.The Pi value of the whole sequence was 0-0.044, and the variant regions were mainly concentrated in the transcribed spacer region (ITS), with a mean value of 0.029 (range 0-0.044) for ITS1, 0.007 (range 00-0.036) for 5.8S rRNA and 0.023 (range 0-0.044) for ITS2.18S rRNA and 26S rRNA were more conserved, with a mean value of 0.001 (range 0-0.025) for 18S rRNA and 0.002 (range 0-0.036) for 26S rRNA.(Figure 5B).Topology of clade D in the chloroplast-based trees and nuclear gene-based trees showed difference to some extent.The number of branches, the composing species of each branch, and the position of the branches, were quite different; for example, the position of B. chaishoui, the relationship between B. microcephalum and B. wenchuanense (See below, Discussion section, and Figure 6).Nevertheless, all four phylogenetic trees showed that the new The nrDNA phylogenetic tree has the highest resolution in differentiating the six Bupleurum species, in which each species formed a monophyletic branch.

Discussion
Differences in gene and structure of chloroplast genomes among Bupleurum species We described for the first time the cp genomes of six species of Bupleurum from western Sichuan Province, China, providing important insights into the cp genome characteristics of members of this genus.Our results show that the cp genomes of Bupleurum species are extremely similar, indicating that they are highly homogeneous in terms of cp genome structure, gene content and tendency of SDRs and SSRs.The cp genomes of the six Bupleurum species differ in size by only a few hundred bp (155,012-155,712 bp), consistent with the description that the cp genomes of Bupleurum species show only minor differences in size (~1 kb) (Huang et al., 2021a).Many mutational events occurring in the cp genome, including substitutions, insertions and deletions (InDels), inversions, genomic rearrangements and translocations (Abdullah et al., 2021), did not occur in any of the six Bupleurum species.
It has been reported that IR region contraction and expansion is a common phenomenon in cp genomes, and that changes in the size
of the cp genome are usually the result of IR region expansion and that these variants can be observed in both closely and distantly related plant species and may lead to pseudogene production, gene duplication and deletion of individual copies of genes (Plunkett and Downie, 2000;Wang et al., 2008c;Saina et al., 2018).Comparative analysis of IR boundaries showed the same distribution of genes at the LSC/IR junction in the cp genomes of the six Bupleurum species, with minor differences in gene length (ycf1 and rps19) at the SSC/IR junction.The ycf1 sequence located at the IRa and SSC boundaries was identified as a pseudogene because it was truncated and was an incomplete duplicate of a normal copy.In the cp genome of a sample (B.wenchuanense voucher 2020080639), expansion of the IR region was found occurred, resulting in inconsistency between LSC/IR and SSC/IR border genes to other Bupleurum.The expansion was due to duplication of genes in IR (rpl22 and rps3).
Variation in the number of duplicated genes in angiosperm cp genomes is common and different numbers of duplicated genes can be observed in different species (Ahmed et al., 2012;Menezes et al., 2018;Abdullah et al., 2020), but IR expansion within species is somewhat uncommon.This phenomenon was once observed in the cp genome of Sinopodophyllum hexandrum (Ma et al., 2022).We performed read mapping analysis and confirmed such an expansion in the very sample of B. wenchuanense (Figure S4).However, the overall gene number and order remains consistent and reflects the significant conserved nature of chloroplasts in this species.

Highly variable regions in cp genome presumed as potential molecular markers
Inspired by Dong's research work (Dong et al., 2012), highly variable regions in cp genome have been generally thought to be potential molecular markers for phylogeny and identification of plant species.In some cases, the highly variable regions worked well in resolving phylogeny and identifying species (Dong et al., 2021;Song et al., 2022), but in some other cases they were not so effective (Wu et al., 2020); more often, their effectiveness has not been evaluated (Liu K et al., 2021;Liu X.J. et al., 2021;Wu et al., 2022).In order to assess whether the they can work in Bupleurum phylogenetic analysis and species identification, we checked the highly variable regions in cp genomes of the six Bupleurum species.
We found that the non-coding intergenic spacer regions have considerable diversity over the protein-coding regions, which was consistent with most previous studies (Daniell et al., 2016;Javaid et al., 2022).Seven spacer regions (rpl32-trnL, petA-psbJ, atpF-atpH, trnK-rps16, atpH-atpI, ndhC-trnV, and rps15-ycf1), and two coding regions (trnL-UAA and ycf1) were found to be highly variable (Pi > 0.01).However, the discovered highly variable regions were not all the same as the previously reported ones in Bupleurum (Li et al., 2020;Huang et al., 2021a;Huang et al., 2021b;Xie et al., 2021).The spacer regions petA-psbJ, trnK-rps16 and the gene ycf1 are the highly variable regions found in common in these studies on Bupleurum cp genomes.There is evidence supporting the ycf1 gene as one of the core plastid DNA barcodes in land plants (Neubig et al., 2009;Franck et al., 2012;Dong et al., 2015).The rpl32-trnL, petA-psbJ gene have also been reported to show very high nucleotide per site diversity in several genera (Dong et al.,FIGURE 4 Comparison of the six Bupleurum species cp genomes using the mVISTA alignment program, with B. chinense as a reference. 2012; Tan et al., 2020;Wang et al., 2022;Niu et al., 2023), and indeed the Pi values of the rpl32-trnL spacer regions of these six Bupleurum species were also much higher than the other regions in the present study, so rpl32-trnL together with petA-psbJ, trnK-rps16, and ycf1 were considered as having potential to be used as high-resolution DNA barcodes to identify species of Bupleurum.
We constructed phylogenetic trees of the six Bupleurum species using the four screened potential molecular markers independently and in combination, and found the trees were in similar topology to the whole cp genome and CDS trees, i.e., the six species were divided into two clades.But in each tree, there's species group that cannot be clearly separated (Figures S3A-F).The "potential" disappears in reality.Therefore, we propose that in similar studies in the future such presumed molecular markers should be verified.

Phylogenetic relationship inferred from cpDNA and nrDNA
The backbone of phylogenetic trees constructed based on the complete cp genomes, protein-coding sequences, nrDNA of the six Bupleurum species were highly silimar; and they were consistent with those acquired in the previous studies (Neves and Watson, 2004;Wang et al., 2008a;Wang et al., 2008b;Wang et al., 2011).The results confirmed Bupleurum is a monophyletic group and supported the treatment of its independent tribe status in subfamily Apioideae.The two main subclades represented the two subgenera in the genus, subgenus Penninervia and subgenus Bupleurum, respectively (Neves and Watson, 2004;Wang et al., 2011).The two clades (clade C and clade D in Figure 6) of Chinese Comparative analysis of the nucleotide diversity (Pi) value among the six Bupleurum species (A) cp genome; (B) nrDNA X-axis: nucleotide position, Y-axis: nucleotide diversity value.
Bupleurum, as well as the species included in each one, were also in concordance with the earlier researches (Wang et al., 2008b;Wang et al., 2011).

B. microcephalum
In the six Bupleurum species in Western Sichuan, B. microcephalum and B. malconense were considered to be the most closely related species (Li and Sheh, 1979;She and Watson, 2005).
The leaves of B. microcephalum are linear, extremely long, soft and thin, while the leaves of B. malconense are more than half shorter than that of B. microcephalum, and the texture is hard.But in all phylogenetic trees, B. microcephalum forms a clade with B. wenchuanense, and B. pseudochaishoui; in the tree based on nrDNA and ITS, it forms a sister branch to the latter two species.This is not consistent with the previous concepts that it's close to B. malconense.The three species in this clade all have solitary stem, which is different from species that have caespitose rootstock and numerous stems in the other clade.B. sichuanense and B. malconense sichuanense was established in 1992 by Pan and Hsu (Pan and Hsu, 1992).In the original species description, it was pointed out B. sichuanense was closely related to B. malconense, and was differ from the latter in its short and tender cauline leaves, more umbel rays (6~7), and bracteoles not longer than the pedicel (< 2 mm).In Flora of China, it was synonymized with B. malconense (She and Watson, 2005).Our molecular phylogenetic analysis proved their intimate relationships.In the phylogenetic trees constructed with the whole cp genomes, the CDS regions and the high variable regions, neither of the two formed a separate clade.In the nrDNA ITS tree, more than 10 sequences derived from these two species tangled together; but in the tree based on the nrDNA, B. sichuanense and B. malconense were clearly separated, implying they were two independent species.We checked the specimens of B. malconense and B. sichuanense deposited in Chinese Virtual Herbarium and those we collected and were of the opinion that the nrDNA correctly reflects the relationship between the two species.The results also implied the nrDNA might serve as a more suitable molecular marker for Bupleurum phylogenetic analysis and species identification.
The close relationship of B. sichuanense and B. malconense to B. chinense and B. yinchowense was disclosed almost 15 years ago (Xie et al., 2009).Their intimate relationship has once again been proven in this study.B. chinense is the officially recognized origin species of the crude drug of Chaihu (Chinese Pharmacopoeia Commission, 2020).There is a need for molecular marker to distinguish B. chinense from the other three closely related species (Chao et al., 2014).It seems nrDNA and ITS perform better than cpDNA does for such a purpose.

B. chaishoui and the nucleocytoplasmic conflict
Phylogenetic trees constructed from cp genome and nrDNA show different topologies, indicating the existence of nuclear cytoplasmic conflicts.The position of B. chaishoui is an obvious example.Phylogenetic trees constructed from both nrDNA and ITS sequences support that B. chaishoui is closely related to B. sikangense and B. commelynoideum, while in the phylogenetic trees constructed from the whole cp genome sequence and CDS sequence, B.chaishoui is closely related to B. malconense; the latter is consistent with the opinion in FRPS that stated B. chaishoui and B. malconense are similar in morphology (Li and Sheh, 1979).The two species, together with the third species in the same clade, B. sichuanense, share a common morphological feature, i.e., the caespitose rootstock.
The inconsistencies in the topology of the gene trees acquired from the two different set of chloroplast and nuclear genomic data is a widespread phenomenon caused by a variety of reasons, including gene duplication and loss, horizontal gene transfer, incomplete lineage sorting, and hybridization introgression, etc. (Steenwyk et al., 2023).The most common reason for this phenomenon was believed to be chloroplast capture caused by hybridization (Fehrer et al., 2007;Galtier and Daubin, 2008;Soltis and Soltis, 2009;Stull et al., 2020).It was also reported cp genome-based phylogenetic relationships were often found to be related to geography rather than morphology (Liu et al., 2020).The Hengduan Mountains-Western Sichuan area is one of the diversity centers of Bupleurum species, and more than half of the species in China can be found in this region.The complex terrain and landforms, and diverse climate of ancient and modern times in the region have created opportunities for natural interspecific hybridization, and the results of phylogenetic trees suggest the possibility of ancient chloroplast capture events (Wang et al., 2008b;Sun et al., 2017;Yang et al., 2019;Xie et al., 2021).Moreover, the widespread presence of hybridization in biological taxa has gained increasing recognition (Fehrer et al, 2007;Acosta and Premoli, 2010).However, compared to its complex biodiversity, there are still very few clearly described hybridization events in the genus Bupleurum.In-depth researches in this aspect are expected in the future to reveal the important value of hybridization for the formation and maintenance of biodiversity of Bupleurum.

The overlooked new species
The materials collected from Weizhou Town and Shuixi Town of Wenchuan (voucher Chao Zhi 1682502 and Chao Zhi & Huang Rong 2020080849) had been overlooked and mistakenly acknowledged as B. chaishoui for years, recently we realized it represents a new species and named it B. pseudochaishoui.These two samples formed distinct clade in all phylogenetic trees, supporting its status as an independent species.
B. pseudochaishoui shows distinguishable morphological character to its allies.Its leaves are coriaceous and some shiny.The basal leaves are long petiolated and the blades are quite large, blade broadly ovate-elliptic to lanceolate elliptic; the middle and upper cauline leaves are much smaller, long lanceolate, with no petiole (Figure 7).In Chinese Bupleurum, only B. longiradiatum Turcz.(Northeastern, Northern and Northwestern China) and B. aureum Fisch.ex Hoffm.(Xinjiang) have basal leaves with blades of such a similar shape and size to B. pseudochaishoui.
Molecular data indicated that B. wenchuanense is the most intimately related species to B. pseudochaishoui, but it differs from the new species in its numerous, rosette-caespitose oblanceolate basal leaves and subulate to squamose middle and upper leaves, and the (1-)2-3 very unequally rayed umbels, 1-4-flowered umbellules.
B. chaishoui, the species which we previously misrecognized as B. pseudochaishoui, has caespitose numerous stems, and lanceolate to elliptic cauline leaves; the shape and size of cauline leaves vary greatly, especially those born on the same node or adjacent nodes are very unequal; the upper ones often grow in a reflexed way (Figure 8).On the contrary, B. pseudochaishoui has solitary stem, and its cauline leaves often are not reflexed.
In the areas neighbor to Southwestern China, there are another two Bupleurum species having basal leaves with broadly oval lamina, namely B. lanceolatum Wall.ex DC. in Himalayas in India and Western Pakistan and B. gilesii H. Wolff in Western Pakistan and Afghanistan (Missouri Botanical Garden, 2023).The upper cauline leaves are of similar shape to the basal and lower leaves in B. lanceolatum, and short condensed main stem and numerous branches from the base in B. gilesii, make them easily to be distinguished from B. pseudochaishoui.

Conclusion
Comparative analysis of the cp genomes and nrDNA sequences of six Bupleurum species from western Sichuan Province, China revealed their conserved structure.Phylogenetic analysis based on nrDNA proved the monophyly of each species.Incongruence between nuclear and plastome gene trees was observed, which may be attributed to chloroplast capture.For overall consideration, B. sichuanense was thought to be an independent species different from B. malconense; B. chaishoui was supposed to be close relative to these two species, and they all shared the characteristic of having caespitose numerous stems.The other three species with solitary stem, B. microcephalum, B. wenchuanense and B. pseudochaishoui, were close in consanguinity.The new species, B. pseudochaishoui, was described and illustrated.
Habitat and distribution:-The new species was found only in its type locality.It inhabits shrubs on mountain slopes, at an elevation about 2000 m.Fl. and fr.July-Sept.
Etymology:-The epithet refers to it has been previously mistaken as B. chaishoui.

Figures
Figures6A-Dshows the results of maximum likelihood (ML) analysis and Bayesian inference (BI) analysis of chloroplast whole genomes, protein-coding genes, the nrDNA and ITS region respectively, for the studied Bupleurum species.The phylogenetic trees estimated by ML and BI analyses showed similar topologies with high bootstrap support values and strong posterior probabilities for most nodes.The trees constructed based on the cp genome, protein-coding genes, and nrDNA were essentially alike.Bupleurum species formed a monophyletic clade which was sister to the other three species of Apioideae and could be divided into two main subclades, one consisting of B. fruticosum and B. gibraltaricum from the Mediterranean region (clade A), and the other consisting of all the remaining Chinese species (clade B).The clade of Chinese species in subgenus Bupleurum can be further divided into clades C and D. In all trees, the constituent species of these two clades are consistent.A small one (clade C) was composed of B. hamiltonii N.P. Balakr., B. candollei Wall.ex DC., B. yunnanense Franch., B. marginatum Wall.ex DC. and its variety B. marginatum var.stenophyllum (H.Wolff) R.H. Shan & Y. Li.A large one (clade D) included more than 20 species; the six species produced in western Sichuan are all in this large clade.Topology of clade D in the chloroplast-based trees and nuclear gene-based trees showed difference to some extent.The number of branches, the composing species of each branch, and the position of the branches, were quite different; for example, the position of B. chaishoui, the relationship between B. microcephalum and B. wenchuanense (See below, Discussion section, and Figure6).Nevertheless, all four phylogenetic trees showed that the new

FIGURE 1 cp
FIGURE 1 cp genome map of the six Bupleurum species Genes shown outside of the larger circle are transcribed clockwise, while genes shown inside are transcribed counterclockwise.Thick lines of the smaller circle indicate IRs and the inner circle represents the GC variation across the genic regions.

FIGURE 3
FIGURE 3Comparison of the LSC, SSC and IR junction among the six Bupleurum species JLB, junction line between LSC and IRb; JSB, junction line between SSC and IRb; JSA, junction line between SSC and Ira; JLA, junction line between LSC and IRa.
FIGURE 6 Phylogenetic trees showing the relationship between Bupleurum species.BS>70, PP>0.95 obtained by ML/BI methods are shown above the branches.(A) Phylogenetic tree based on the whole cp genome; (B) Phylogenetic tree based on protein-coding genes; (C) Phylogenetic tree based on the nrDNA; (D) Phylogenetic tree based on the ITS region.

TABLE 1
Sources of material and GenBank accession numbers.

TABLE 2
Best models for ML and BI analyses.

TABLE 3
Summary on chloroplast features of six Bupleurum species.

TABLE 4
List of genes encoded in six Bupleurum plastomes.
FIGURE 2 MAUVE alignment of the six Bupleurum species Locally collinear blocks are represented by continuously colored regions.