Comparative analysis of chloroplast genome and new insights into phylogenetic relationships of Ajuga and common adulterants

Introduction The potential contamination of herbal medicinal products poses a significant concern for consumer health. Given the limited availability of genetic information concerning Ajuga species, it becomes imperative to incorporate supplementary molecular markers to enhance and ensure accurate species identification. Methods In this study, the chloroplast (cp) genomes of seven species of the genus Ajuag were sequenced, de novo assembled and characterized. Results exhibiting lengths ranging from 150,342 bp to 150,472 bp, encompassing 86 - 88 protein-coding genes (PCGs), 35 - 37 transfer RNA, and eight ribosomal RNA. The repetitive sequences, codon uses, and cp genomes of seven species were highly conserved, and PCGs were the reliable molecular markers for investigating the phylogenetic relationship within the Ajuga genus. Moreover, four mutation hotspot regions (accD-psaI, atpH-atpI, ndhC-trnV(UAC), and ndhF-rpl23) were identified within cp genomes of Ajuga, which could help distinguish A. bracteosa and its contaminants. Based on cp genomes and PCGs, the phylogenetic tree preliminary confirmed the position of Ajuga within the Lamiaceae family. It strongly supported a sister relationship between Subsect. Genevense and Subsect. Biflorae, suggesting the merger of Subsect. Biflorae and Subsect. Genevenses into one group rather than maintaining separate categorizations. Additionally, molecular clock analysis estimated the divergence time of Ajuga to be around 7.78 million years ago. Discussion The species authentication, phylogeny, and evolution analyses of the Ajuga species may benefit from the above findings.


Introduction
The Ajuga genus is a member of the Lamiaceae family and encompasses numerous economically and medically significant plants.Many Ajuga species have been employed in Traditional Chinese Medicine for relieving cough, reducing sputum, and arresting bleeding (Dai et al., 2010).Among the Ajuga species, A. bracteosa is particularly prevalent and extensively utilized in folk medicine (Millen et al., 2001).Previous investigation has revealed that the morphologies of most Ajuga species are similar and indistinguishable, potentially leading to misclassification (Talebi et al., 2019).According to the World Health Organization, adulterating herbal products threatens consumer safety (Duan et al., 2017).In recent years, molecular identification techniques, particularly molecular markers, have made notable advancements in traditional Chinese medicine (TCM), which were used to recognize variations among individuals of distinct species or populations (Cheng et al., 2018;Liu et al., 2019).Previous studies have employed the nuclear ribosomal DNA internal transcribed spacer (ITS) 2 to identify certain Ajuga species and related taxa, such as A. ciliate, A. decumbens, and A. lupulina (Han et al., 2008).Nonetheless, some frequent adulterants were not probed, and single-locus DNA barcodes inherently possess limitations.
Additionally, elucidating the taxonomic relationships between species of the genus Ajuga is crucial for understanding and harnessing the medicinal properties of the different species.However, the presence of similar morphological features coupled with the dearth of molecular information has prevented the accurate identification of taxonomic relationships in this genus.Hence, it is essential to develop a more accurate and effective method for identifying and classifying the Ajuga taxa.
The organelle chloroplast (cp) is vital in plant photosynthesis and biochemical processes (Zoschke and Bock, 2018).Compared to conventional DNA portions, the cp genome exhibits higher conservation with slight variations, making it applicable in various research areas, including species authentification and developing DNA markers (Qiao et al., 2016).Plastid genome analysis has been widely employed to identify Paris, Isodon, and respective adulterants and clarify phylogenetic relationships (Jiang et al., 2022;Zhou et al., 2022).Although previous studies have reported plastid genomes of Ajuga (Tao et al., 2019), their focus was primarily on characterizing a single genome, and cp genomes have not been used to differentiate Ajuga taxa and their frequent adulterants.In addition, cp genomes may lead to a wrong inference of phylogenetic relationships due to DNA length disparity, gaps representing insertions/deletions (indels), and improper models of DNA evolution in merged datasets.Previous investigations have indicated that protein-coding genes (PCGs) offer an improved resolution for understanding phylogenetic relationships due to the genetic divergence in gene-encoding regions occurring at a slower rate compared to the non-coding areas (Lockhart and Penny, 2005;Cui et al., 2019).However, no reports exist on using PCGs to comprehend the interspecific relationships among the Ajuga species.
Herein, cp genomes of seven Ajuga species were sequenced, de novo assembled, and annotated.Subsequently, we compared the architecture and evolutionary connections of these genomes.This research aims (i) to contribute novel and full-length plastid genomes of Ajuga to understand more about the genome structure of relevant species, (ii) to expound the phylogenetic relationship of Ajuga by comparing genome sequences, and (iii) to develop promising DNA markers for distinguishing A. bracteosa from its contaminants.Our results significantly increase the genome information of Ajuga, facilitate evolutionary scrutiny and authentification of Ajuga, and ensure the safe utilization of A. bracteosa.

Plant materials and DNA extraction
Fresh leaves from seven species, namely A. forrestii, A. nubigena, A. campylantha, A. macrosperma, A. bracteosa, A. nipponensis, and A. ovalifolia, were gathered from the Germplasm Resource Garden of Kunming Zhifen Biotechnology Co., Ltd, Yunnan, China (102°48′58″E, 24°49′55″N).The voucher specimens were identified by Professor Baozhong Duan and preserved at the herbarium of Dali University.Information on each sample is detailed in Table S1 and Figure S1.Approximately one gram of fresh leaves of each species was collected, instantly frozen in liquid nitrogen, and stored for subsequent DNA extraction.Genomic DNA was extracted from samples with a Plant Genomic DNA kit (Tiangen, Beijing, China) following the manufacturer's instructions.The extracted DNA was checked with high-sensitivity Qubit 4.0 fluorometry (Life Technologies, Inc.).

Identification and validation of barcode for species discrimination
The intergenic spacers (IGS) were obtained from seven Ajuga species with PhyloSuite v1.2.2 (Zhang et al., 2019).Primers were designed based on the variable intergenic regions using Snapgene 6.2.1 (Snapgene, Insightful Science, available at http:// www.snapgene.com,last used in 2023).PCR amplifications were conducted in a final volume of 25 mL, consisting of 12.5 mL of 2×Taq Plus PCR Master Mix, 1 mL of each primer, 2 mL of template DNA, and 8.5 mL of ddH 2 O (Mei5 Biotechnology, Co., Ltd).All amplifications were performed using a RePure-A PCR system (Applied Biogener, Hangzhou, China) under the following conditions: an initial denaturation at 94°C for 3 min, followed by 35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 1 min, with a final extension at 72°C for 5 min.PCR products were examined by 1% agarose gel electrophoresis to confirm the amplification of the target fragments.The purified PCR products were sequenced in both directions on a 3730XL DNA Sequencer (Applied Biosystems, Waltham, USA) using the same primers at Sangon Biotech Co., Ltd.(Shanghai, China).

Phylogeny and divergence time estimation
The phylogenetic analysis involved a total of 35 taxa, consisting of 28 species obtained from NCBI (refer to Table S2) and seven species that we newly sequenced, as detailed in Table 1.The selection of these species for phylogenetic analysis is based on the classification system of Lamiaceae within the Angiosperm Phylogeny Group IV system (APG IV).Additionally, two species, Callicarpa macrophylla (GenBank NC058323.1)and C. arborea (GenBank NC058321.1)served as outgroups.The 68 common PCGs of 35 cp genomes were extracted based on the annotation files.The aligned sequences were generated using the MAFFT program and verified manually.Phylogenetic analysis was performed using Maximum likelihood (ML), Bayesian inference (BI), and Neighbor joining (NJ).For the ML tree reconstruction, IQtree was employed with default settings, 1,000 iterations, and 1,000 replications.Model selection was based on the best-fit approach (Katoh and Standley, 2013).BI analysis was performed using MrBayes v.3.2.6 (Ronquist et al., 2012).The most appropriate model for sequence substitution in plastid genomes (GTR + G + I) and PCGs (GTR + G) was determined using MEGA X v.10.2.6 (Mao et al., 2023).The parameters were set for five million generations, with sampling every 1,000 generations.The initial 25% of each run was discarded as burn-in (Jiang et al., 2023).Moreover, the NJ tree was inferred with MEGA X v.10.2.6 and subjected to the bootstrap test of 1,000 repetitions (Kumar et al., 2018).

Results and discussion
3.1 Genome structure 3.1.1Genome characteristic Around 2.51 -4.54 Gb data were obtained from each species.The A. macrosperma and A. ovalifolia were reported for the first time.As shown in Figure 1, cp genomes of seven species are circular DNAs ranging from 150,342 bp to 150,472 bp and exhibit the typical quadripartite structure commonly observed in most angiosperm cp (Palmer, 1985;Henriquez et al., 2020;Yun and Kim, 2022), consisted of two IRs (IRa and IRb) separated by large single copy (LSC) and small single copy (SSC) regions, respectively.

FIGURE 1
The cp genome map of Ajuga.(Tao et al., 2019), indicating a high degree of conservation among Ajuga species' cp genomes.Notably, the IR regions exhibited a significantly higher GC content of 43.3% in contrast to the LSC (36.4%) and SSC (32.1% -32.3%) regions.This discordance may be attributed to the fact that four ribosomal RNA (rRNA) genes (rrn23, rrn16, rrn5, rrn4.5) with high GC content were located in the IR regions, which was similar to the most plant species (Yan et al., 2019;Gui et al., 2020;Chen et al., 2022).
In addition, 129 -133 genes were detected, including 86 -88 PCGs, 35 -37 transfer RNAs (tRNAs), and eight ribosomal RNAs (Table 1).Variations in gene numbers across these species can be attributed to the expansion and contraction of IRs.These genes can be categorized into three groups: 45 associated with photosynthesis, 27 involved in self-replication, and the remaining genes serving various other functions.Notably, three genes (chIB, chIL, ycf68) were absent in the seven Ajuga species.The missing chIB and chIL may be a distinctive feature of flowering plants (Jansen et al., 2007).Furthermore, the ycf68 gene was also absent in the cp genomes of Miscanthus sinensis and M. floridulus (Sheng et al., 2021).In addition to these observations, it's worth mentioning that two tRNAs, trnF-GAA and trnfM-CAU, were absent in three Ajuga species.Specifically, A. nubigena lacked the trnF-GAA gene, while trnF-GAA and trnfM-CAU genes were missing in A. forrestii and A. campylantha.These findings revealed that the structure of cp genomes was highly conserved in Ajuga species, albeit with some alterations that have accrued during the angiosperm evolution, which was also supported by previous studies (Millen et al., 2001).

Inverted repeats regions contraction and expansion
The expansion/contraction of IR regions is frequently observed during evolution and may account for the disparity in the size of plastid genomes (Menezes et al., 2018).As depicted in Figure 2, the rpl2 gene is entirely situated within the IR regions across all species, while the trnH gene exclusively occupies the LSC region, consistent with the cp genomes of most angiosperms (Wang J. et al., 2022;Fang et al., 2023).Besides, the trnN gene was entirely located within the IRa region in all species except for A. ciliata, A. lupulina, and A. campylanthoides, where it was utterly localized in the IRb region at varying distances from the junction of SSC/IRb.
Among the seven Ajuga species, a truncated copy of ndhF genes was identified at the junction of SSC/IRb, starting from SSC and integrating into the IRb region.Conversely, in the remaining species, including A. ciliata, A. lupulina, and A. campylanthoides, ndhF was found at the junction of IRa/SSC.In addition, the ycf1 gene was observed at the IRb/SSC junction in all species.It originated from the IRb region and integrated into SSC, with sizes ranging from 5 to 4,357 bp.Notably, the ycf1 gene was also present at the IRa/SSC junction in all species except A. bracteosa, in which this gene was exclusively present in SSC.This result implied that the ycf1 could potentially serve as a marker for distinguishing A. bracteosa.A previous study also highlighted ycf1 as a powerful barcode for land plants (Dong et al., 2015).Furthermore, the rps19 gene was consistently located at the boundaries of the IRs/LSC in all Ajuga species.This pattern aligns with the cp genomes of other Lamiaceae species, such as S. mekongensis, Mentha spicata, and Dracocephalum heterophyllum (Wu H. et al., 2021;Fu et al., 2022).In summary, the sizes of the cp genomes in the 11 Ajuga species vary, and there are noteworthy variations in the junction regions.These findings provide evidence of a distinctive pattern of IR contraction/expansion within the cp genomes of Ajuga species, which can be employed to investigate species-specific gene loci.

Codon usage bias of the cp genomes
Analyzing codon usage is essential to evaluate the evolution of the cp genome (Wang N. et al., 2022).In genes of seven Ajuga species, 64 codons were identified, of which 61 encoded 24 amino acids.Leucine exhibited the highest frequency among all the amino acids encoded by cp genomes, whereas cysteine was found to be a rare amino acid (Table S4).This observation is consistent with the codon usage bias reported by previous studies (Oresǐčand Shalloway, 1998;Liu et al., 2020).Moreover, pronounced bias towards A or T at the tertiary position of codon was observed, which could be attributed to the high AT proportion in plastid genomes.Similar results were observed in other angiosperm taxa (Wu L. et al., 2021;Jiang et al., 2022).As shown in Figure 3, UUA had the highest frequency, followed by AGA, while GCA had the lowest frequency.In Ajuga, RSCU values of 30 codons were higher than 1.00, 32 codons had values below 1.00, and two had values of 1.00.An RSCU value below 1.0 suggests that the codon usage frequency is less than expected, while an RSCU value over 1.0 means that the codon usage frequency is more than expected (Ma et al., 2015;Parvathy et al., 2022).This variability in RSCU values reflects evolutionary information resulting from mutation and selection, which is essential in studying organismal evolution (Morton, 2003).
Additionally, the GC proportion of the GC3s was closely associated with codon bias and was an important parameter for evaluating the codon use pattern (Ma et al., 2015;Parvathy et al., 2022).In Ajuga, GC3s values ranged from 26.9% to 27.0%, indicating that the genus Ajuga preferred A/U ending codons (Zhang et al., 2022).Previous studies have highlighted that high AT content is the primary reason for synonymous codons ending in A/ U, potentially linked to natural selection and mutation during evolution (Zhang et al., 2018).Additionally, the Enc varied between 49.80% and 49.94%, while the CAI and the optimal frequency were lower than 0.5.These results suggested a slight bias in codon usage within the seven Ajuga taxa.

Repeat sequences
Large and complex repeat sequences in the cp genome are potential markers for revealing gene rearrangements and losses during evolution (McDonald et al., 2011;Yi et al., 2013).Analysis of oligonucleotide repeat in seven cp genomes indicated that the number and length of repeat sequences differed among genomes and were distributed randomly, with repeat sequences ranging from 30 to 82 bp and most being within 30 -46 bp (Figure 4B).Meanwhile, 274 long repeats were identified, including 149 P repeats, 124 F repeats, and 1 R repeat, with P and F being more than R and C, a pattern consistent with most plastid genomes of angiosperms (Vu et al., 2020;Luo et al., 2022).R repeats were only present in A. campylantha, while C repeats were absent in all seven Ajuga species (Figure 4A).These findings provide a molecular basis for identifying the Ajuga species.
SSRs, also known as microsatellites, are popular genetic indicators because of their significant polymorphism, repeatability, and reliability, which can be used to detect genetic diversity, population, and polymorphisms at intraspecific, distant phylogenetic relationships and cultivar levels (Wang et al., 2009;Peng et al., 2022). 34, 36, 28, 33, 34, 33, and 33 SSRs were found in A. forrestii, A. nubigena, A. campylantha, A. macrosperma, A. bracteosa, A. nipponensis, and A. ovalifolia, respectively (Figure 5A).Among these, mononucleotide (A/T/C) repeats were more than other types of repeats, accounting for almost 62%, as previously reported by Zhou et al. (Zhou et al., 2022).The second most common was tetranucleotide repeat (23.53% -29.41%), with a predominant motif of AAAG/CTTT and AAAT/ATTT, followed by dinucleotide repeats (8.33% -14.29%), with a dominant motif of AT/AT.Hexanucleotides (2.78% -3.03%) were only present in plastomes of A. nubigena (AACTAT/AGTTAT) and A. nipponensis (AAAAAT/ATTTTT), while trinucleotide or pentanucleotide repeats were absent in all seven Ajuga species.These findings indicate that mononucleotide repeats are more frequent than other types, and A/T motifs were the most abundant in the mononucleotide repeats, which is consistent with previous studies (Munyao et al., 2020;Ren et al., 2022).It was suggested that the high amount of mononucleotide repeats in the cp genome may contribute to heritable variations (Bi et al., 2018).
Additionally, the number of SSRs exhibited significant variation across distinct structural and functional regions of the cp genomes (Sun et al., 2022).Frequency analysis revealed that SSRs were more prevalent in LSC (44.55%) than in IR (28.79%) and SSC (26.66%) regions (Figure 5B).Besides, most SSRs were located in gene regions and IGS, with an average number of 16 and 9, respectively.In contrast, the regions encompassing introns and exons contained the fewest SSRs, with an average count of 3 and 2, respectively (Figure 5C).These findings are consistent with those observed in other Lamiaceae species (Fu et al., 2022).Previous research has emphasized the suitability of SSRs as a genetic marker in plant molecular studies (Khan et al., 2019), particularly in noncoding regions exhibiting high intraspecific variation (Eguiluz et al., 2017).Our investigation indicated that most SSRs were situated within non-coding regions, with a limited presence in the coding areas.Consequently, these SSRs have the potential to serve as markers for discerning various evolutionary changes, such as genetic diversity, and they may even facilitate species differentiation within Ajuga.

Comparing genomes and nucleotide diversity
The comparative analysis of cp genomes is a practical approach to investigating the genetic structure and phylogenetic kinships of plants (Daniell et al., 2016).Overall sequence variation in plastid genomes of Ajuga indicated a high level of conservation, with the protein-coding regions exhibiting more significant conservation than non-coding regions, except ndhF, ycf1, and ycf2 genes (Figure 6).This observation aligns with a previous report on the cp genome of S. miltiorrhiza within the Lamiaceae family (Qian et al., 2013).Besides, the most significant discrepancy was mainly found in IGS, e.g., trnH(GUG)-psbA, rps16-trnQ(UUG), atpH-atpI, rbcL-accD, ndhC-trnV(UAC), accD-psaI, trnF(GAA)-ndhJ, and ndhF-rpl32.Previous studies have identified IGS as hotspots with numerous nucleotide substitutions and indel mutations, making them valuable markers with high resolution for phylogenetic analyses (Drummond, 2008).For example, rps16-trnQ(UUG) is highly variable in most plants and has been utilized for DNA barcoding in phylogenetic studies across various angiosperm genera (Dong et al., 2012).Similarly, rbcL-accD and trnH(GUG)-psbA have been proposed as critical molecular markers for phylogenetic analyses of Viola species (Cao et al., 2022).Consequently, these highly variable regions are expected to offer ample genetic information for conducting studies on species delimitation and the phylogenetic evolution of Ajuga.

IGS-based species authentication
IGS is often employed as indicators in phylogeny inference at various taxonomic hierarchies, which is more variable and can offer more evolutionarily revealing characters (Fang et al., 2023).Previous molecular studies of the genus Pogostemon (Zhang et al., Number of long repetitive repeats on the cp genome of seven Ajuga species.(A) the number of repeat types (P-palindromic repeats, F, forward repeats; R, Reverse); (B) Frequency of the repeats more than 30 bp long.2020) and Clerodendrum (Chen et al., 2023) have demonstrated the high identification capabilities of cp genetic markers.In this study, eight IGSs were extracted from seven Ajuga species.ML analyses were conducted on each IGS using IQtree (Nguyen et al., 2015).As illustrated in Figure S2.1-2.8, the results showed that A. bracteosa could be differentiated from its common adulterants based on ndhF-rpl32, ndhC-trnV(UAC), accD-psaI, atpH-atpI, and trnH (GUG)-psbA, while the remaining IGSs were incapable with weak bootstrap values (<70).Previous studies have also identified the ndhF-rpl32 region as a useful molecular marker for distinguishing Magnolia polytepala and its closely related species (Sun et al., 2020) and ndhC-trnV(UAC) for distinguishing Isodon rubescens and its adulterants (Zhou et al., 2022).Moreover, accD-psaI or atpH-atpI were probable markers for identifying additional species (Alzahrani, 2021;Zheng et al., 2022).
Although general DNA barcodes (e.g., ITS) can distinguish A. ciliata and related taxa (Han et al., 2008), some common adulterants were not investigated.Our results revealed that the studied IGSs exhibited more variability than ITS.The ML tree was constructed based on five IGSs that could differentiate A. bracteosa from its common adulterants.The tree demonstrated that all A. bracteosa species formed a monophyletic clade, and A. macrosperma formed independent branches.Strong support was observed for a sister relationship between A. macrosperma and A. bracteosa (Figure S3).These findings indicate that five combining IGSs can successfully distinguish A. bracteosa and its frequent dopants.
Additionally, we designed primers for five IGS and conducted amplification and sequencing experiments (Table 2).It is worth noting that due to the lack of residual DNA in three Ajuga samples, our supplementary investigation was limited to four species with remaining DNA (A.bracteosa, A. forrestii, A. macrosperma, and A. campylantha).The results demonstrated that except for the trnH (GUG)-psbA primers, the other four fragments (accD-psaI, atpH-atpI, ndhC-trnV(UAC), and ndhF-rpl23) produced products of the expected sizes in the selected Ajuga species (refer to Figure S4).Both amplification and sequencing achieved a 100% success rate.The sequence data from the four Ajuga species align with the cp genome results, with each species exhibiting distinct base differences (refer to Figure S5.1-5.4).These findings confirm that the four DNA barcodes are ideal tools for distinguishing A. bracteosa from its adulterants.

Phylogenetic analysis
Cp genomes have a wealth of phylogenetic information and are extensively employed for reconstructing phylogenies and conducting plant population analyses (Yan et al., 2019).Here, ML and BI trees were reconstructed using cp genomes and common PCGs of 35 species, respectively.Figure 8 illustrates ML and BI trees based on 68 shared PCGs, showing similar topologies.Among the three subfamilies, Ajugoideae and Lamioideae Harley emerged as sister taxa, while Nepetoideae (Dumort.)Burnett appeared as the Phylogenetic analysis also demonstrated that all A. bracteosa species formed a monophyletic clade.A. bracteosa (GenBank OR038702) was an independent branch of the phylogeny and deeply nested within clade B, with strong support (BS = 100) for the sister relationship between A. bracteosa and A. macrosperma, suggesting that the cp genome could distinguish A. bracteosa from other species.Notably, the samples of A. forrestii (Genbank MN814855.1)did not form a monophyletic group and were Global comparison of complete genomes of Ajuga.Previous studies also confirmed that intraspecific diversity existed in Artemisia argyi (Chen et al., 2022), Isodon rubescens (Zhou et al., 2022), and Phyllanthus urinaria collected from different geographical areas (Fang et al., 2023).This phenomenon could potentially be attributed to the influence of the geographical area of origin on the variation in A. forrestii.
Additionally, a strong support value of 100 was observed for the sister relationship between A. nubigena from Subsect.Biflora and A. ovalifolia from Subsect.Genevenses, contradicting the taxonomic findings presented in Flora of China (FOC) (Li and Hedge, 1994).The BI tree also confirmed the same result with the posterior probabilities value of 1.This finding suggests a potential limitation in the current classification of Subsect.Biflora as independent entities.Therefore, based on this evidence, we propose an alternative classification scheme that amalgamates Subsect.Biflorae and Subsect.Genevenses into a single group rather than maintaining separate categorizations.
Notably, Clerodendreae is classified as part of the Verbenaceae family in the FOC (Chen and Gilbert, 1994).Our study confirms its placement within the subfamily Ajugoideae, aligning with the APG IV classification (Angiospermm, 2016).Therefore, we recommend reclassifying the genus Clerodendreae under Lamiaceae.
Furthermore, the phylogenetic trees constructed using cp genomes (Figure S6.1-6.2) and common PCGs (Figure 8 and Figure S6.3) showed a high degree of similarity.Previous studies have suggested that using cp genomes may lead to missing relationships due to length variations, gaps/index deletions, and inappropriate models of DNA evolution in concatenated datasets (Goremykin et al., 2005;Lockhart and Penny, 2005).Given that the genetic divergence in gene-encoding regions occurs more slowly than in non-coding sequences, we considered that utilization of common PCGs is more appropriate for the identification and phylogeny analysis of Ajuga.
In conclusion, our findings serve as a valuable basis and reference for utilizing plastid genomes and common PCGs in species identification, contributing to a better understanding of Ajuga's phylogeny.

Divergence time
The 35 cp genomes from Lamiaceae family plants, including 16 Ajuga species, were utilized to estimate the divergence time based on the ML tree.As shown in Figure 9, the subfamilies of Ajugoideae and Lammiodeae shared a common ancestor in the late Eocene (39.60 Mya), and the split between Ajuga, Clerodendrum, and Rotheca could occur in the Oligocene (38.65 Mya).The process of speciation of the Ajuga genus was estimated to originate at 7.78 Mya in the late Miocene.Notably, the Three main lineages (clade I: 6.38 Mya; clade II: 0.84 Mya; clade III: 0.61 Mya) within the Ajuga genus diverged in the late Miocene and continued throughout the Pleistocene.This time frame aligns with the final stages of rapid uplift in the Qinghai-Tibetan Plateau (QTP), which has been recognized as a region abundant in biodiversity and a source of various herbal resources (Wen et al., 2014;Favre et al., 2015).Biasatti et al. also suggested that QTP could have influenced the geographical environment, climate, and distribution/divergence of species (Biasatti et al., 2012).These findings imply that the interspecific divergence of the Ajuga species may have a close association with the uplift of the QTP.
In addition, the Pleistocene has been proposed as a critical period for refugial isolation and subsequent lineage formation, leading to modern species diversity (Johnson and Cicero, 2004).According to the Pleistocene speciation model, glacial cycles during this era acted as a 'species pump,' contributing significantly to the diversity of organisms inhabiting temperate regions (Pyron and Burbrink, 2009).Previous research has also indicated that the extreme climatic fluctuations of the Pleistocene played a pivotal role in driving diversification between lineages in specific taxa (Hewitt, 1999;Bryson et al., 2012).
Consequently, it is postulated that the intense uplift of the QTP, coupled with the climatic oscillations of the Pleistocene, may have influenced diversification and facilitated the radiation of Ajuga species.

Conclusions
In this study, plastid genomes of seven Ajuga species were de novo assembled based on short sequencing reads, and cp genome sequences of A. macrosperma and A. ovalifolia were reported for the first time.These plastid genomes were broadly conserved, displaying comparable gene organization and content.We demonstrated the utility of PCGs integration in phylogeny investigations of Ajuga.Phylogeny analysis of 68 common PCGs strongly supported the taxonomic placement of Ajuga within the Lamiaceae family.It explicitly supported a sister relationship between A. nubigena from Subsect.Biflora and A. ovalifolia from Subsect.Genevense.Consequently, we propose amalgamating Subsect.Biflorae and Subsect.Genevenses into a single group, advocating against their separate categorization.Four highly variable cp loci, including atpH-atpI, accD-psaI, ndhC-trnV(UAC), and ndhF-rpl23, were identified, which hold promise as markers for distinguishing A. bracteosa from its common adulterants.The divergence time of Ajuga occurred in the early Pliocene, possibly due to the intense uplift of QTP and the global cooling event.In summary, this study provides a valued reference for ensuring the efficacy and safety of clinical application while also facilitating bioprospecting and conservation efforts concerning the Ajuga species.

FIGURE 2
FIGURE 2Comparisons of the borders of LSC, SSC, and IRa/b regions among the 11 Ajuga plastid genomes.

FIGURE 3
FIGURE 3Heat map of the RSCU values among Ajuga cp genome.
FIGURE 5 Compares SSR distribution in the cp genomes of seven species.(A) the number of different SSR types; (B) Frequency of SSRs in the LSC, IR, and SSC region; (C) Number of SSRs in the intergenic regions, PCGs, and introns.

FIGURE 7
FIGURE 7Sliding window analysis of Ajuga cp genome.

FIGURE 8
FIGURE 8ML and BI phylogenetic tree based on 68 common PCGs of 35 species.

FIGURE 9
FIGURE 9Divergence times estimation based on cp genomes.

TABLE 1
Summary of cp genome features.
The length of LSC regions ranged from 82,080 bp (A.campylantha) to 82,170 bp (A.ovalifolia), SSC regions ranged from 17,183 bp (A.nubigena) to 17,153 bp (A.macrosperma), and IRa and IRb regions ranged from 25,573 bp (A.nipponensis) to 25,544 bp (A.campylantha).The overall GC content of seven Ajuga cp genomes was 38.3%, largely concordant with the prior study ofTao et al.

TABLE 2
Primer design by SnapGene.