Fresh leaves of Glycyrrhiza uralensis Fisch. ex DC., Glycyrrhiza glabra L., Glycyrrhiza inflata Batalin, and Glycyrrhiza pallidiflora Maxim. were obtained from the Beijing Medicinal Plant Garden, and Glycyrrhiza squamulosa Franch. and Glycyrrhiza yunnanensis S.H.Cheng & L.K.Tai ex P.C.Li were collected from Hengshui City in Hebei Province and Lijiang City in Yunnan Province, respectively. All samples were identified by Prof. Yulin Lin from the Institute of Medicinal Plant Development (IMPLAD). Voucher specimens were deposited in the herbarium at IMPLAD. Total DNA of the six species was extracted by a DNeasy Plant Mini Kit (Qiagen Co., Germany), and the DNA concentration and quality were assessed through Nanodrop 2000C spectrophotometry and electrophoresis in 1% (w/v) agarose gel, respectively.

The DNA of six Glycyrrhiza species was used to generate libraries with average insert size of 500 bp and sequenced using Illumina Hiseq X in accordance with the standard protocol. Low-quality reads resulting from all samples were trimmed by Trimmomatic software (Bolger et al., 2014). GetOrganelle (Jin et al., 2020) and NOVOPlasty (Dierckxsens et al., 2017) were used to assemble the chloroplast genomes. GeSeq (Tillich et al., 2017), CPGAVAS2 (Shi et al., 2019), and PGA (Qu et al., 2019) software were used to annotate the sequences initially and correct them manually. The final complete chloroplast genomes of G. uralensis, G. glabra, G. inflata, G. pallidiflora, G. squamulose, and G. yunnanensis were submitted to GenBank and obtained the accession numbers of PP119344, PP119342, PP119340, PP119341, PP119343, and PP119345, respectively.

Chloroplast genome maps were generated using Chloroplot (Zheng et al., 2020) and then manually corrected. CodonW software (Sharp and Li, 1987) was adopted to analyze the usage of codon. Simple sequence repeats (SSRs) were detected by MISA (Beier et al., 2017) with the definition of ≥10 repeat units for mononucleotide; ≥5 repeat units for dinucleotide; and ≥4 repeat units for trinucleotide; and ≥3 repeat units for tetranucleotide, pentanucleotide, and hexanucleotide SSRs. The long repeated sequences were detected by REPuter (Kurtz et al., 2001) with length of ≥30 bp and over 90% similarity between two copies.

The chloroplast genomes of all published Glycyrrhiza species (Supplementary Table S1) were compared using the online genome comparison tool mVISTA software (Frazer et al., 2004). The nucleotide diversity values (Pi) of the chloroplast genomes were computed using DnaSP (Librado and Rozas, 2009). KaKs_Calculator v. 2.0 was used to determine the ratio of non-synonymous substitutions (Ka) and synonymous substitutions (Ks) (Wang et al., 2010). Chloroplast genome sequence homology and collinearity were analyzed using Mauve software (Darling et al., 2004). A total of 87 chloroplast genomes of Fabaceae species, including all tribes of Fabaceae with published chloroplast genomes (Supplementary Table S1), were used to construct a phylogenetic tree. MAFFT software (Katoh and Standley, 2013) was used to compare the common protein-coding genes shared by the chloroplast genomes of these Fabaceae species. Maximum likelihood (ML) analysis was carried out by IQ-TREE (Nguyen et al., 2015) based on these common protein-coding genes with a bootstrap of 1000 repetitions. ML analysis was conducted based on the GTR + F + R5 model. The MCMCTREE program in the PAML package (Yang, 2007) was used to estimate the differentiation time of Fabaceae. The topology of the legume ML tree was used to demarcate the differentiation time nodes, and the results were viewed and edited by FigTree software (http://tree.bio.ed.ac.uk/software/figtree/).

Correctly assembling complete chloroplast genomes ensures the reliability and reproducibility of research on the evolution of eukaryotic chloroplast genomes, as well as downstream studies based on chloroplast genomes. Currently, the most widely used tools for assembling organelle genomes are GetOrganelle and NOVOPlasty. This study employed GetOrganelle and NOVOPlasty for the assembly of six Glycyrrhiza species. The results revealed that the species successfully assembled using GetOrganelle were G. uralensis, G. squamulosa, and G. yunnanensis, whereas those successfully assembled using NOVOPlasty were G. glabra, G. inflata, G. pallidiflora, and G. squamulosa. The final assembled chloroplast genome sizes for the six Glycyrrhiza species ranged from 126,380 bp (G. glabra) to 129,115 bp (G. yunnanensis), with an overall GC content of approximately 34%. The coding sequences (CDS) accounted for the largest proportion, comprising around 52% of the total genome. The chloroplast genomes of the six Glycyrrhiza species showed a circular structure (Figure 1). Similar to the published chloroplast genomes of Glycyrrhiza species, these six chloroplast genomes also lacked the IRs. The success rates of chloroplast genome assembly for these six Glycyrrhiza species were similar between GetOrganelle and NOVOPlasty. However, a study based on testing with a publicly available dataset of reads from 50 plant species showed that, under slightly high computational resource consumption, the default parameter of GetOrganelle achieved a much higher complete circularity rate compared with the best parameter result of NOVOPlasty (Jin et al., 2020). In the benchmark study evaluating mainstream chloroplast genome assembly tools, GetOrganelle achieved significantly higher circularity and accuracy rates compared with other tools; thus, it was recommended as the optimal assembly tool (Freudenthal et al., 2020). Additionally, the output results from GetOrganelle included excellent visualization, allowing for an intuitive assessment of the assembly success without the need for additional verification. This feature significantly reduces the occurrence of false positives.

For chloroplast genome annotation, we employed the widely used organelle genome annotation tools, namely, GeSeq, CPGAVAS2, and PGA. The results indicated that the annotation effectiveness varied among different software tools for the six Glycyrrhiza species (Supplementary Table S2). GeSeq provides a comprehensive annotation of genes, and its real-time updated organelle genome database facilitates the selection of closely related species as references, which is crucial for the accuracy of chloroplast genome annotation. Combining these three annotation tools with manual correction, a total of 109–110 genes were annotated, including 76 protein-coding genes, 29–30 tRNAs, and 4 rRNAs, similar to the published Glycyrrhiza chloroplast genomes (Supplementary Table S3). Introns can accumulate more mutations than exons and play an important role in gene expression regulation (Kelchner, 2002). Previous studies have shown that introns can increase the expression level of foreign genes in eukaryotic genomes (Xu et al., 2003). In the six Glycyrrhiza species, there are a total of 16 genes with introns, with one gene (ycf3) having 2 introns, and the remaining 15 genes (atpF, ndhA, ndhB, petB, petD, rpl16, rpl2, rpoC1, rps12, trnK-UUU, trnV-UAC, trnL-UAA, trnG-UCC, trnI-GAU, and trnA-UGC) each having 1 intron (Supplementary Table S4). There were two special genes: rps12 gene was a trans-splicing gene; trnK-UUU gene sequence contained the matK gene, similar to that of many other plants (Wu et al., 2021).

In summary, various chloroplast genome assembly and annotation tools yield varied results for different species or datasets, even within the same genus. Therefore, it is recommended to use multiple tools for chloroplast genome assembly and annotation, and the results from different tools should be combined to ensure the accuracy of chloroplast genome assembly and annotation.

Codon usage bias refers to the phenomenon that all synonymous codons in protein transcripts are not uniformly used to encode all amino acids in the protein except methionine and tryptophan (Jia and Xue, 2009). Given the difference in the degree of variation and selection pressure, the variation in compositional constraints between different genomes is a key factor in the formation of codon usage bias (Sharp and Li, 1987; Hershberg and Petrov, 2008; Paul et al., 2018). Codon usage bias can be used to investigate the evolutionary history of organisms, predict the expression level, and understand the evolutionary processes acting on genomes at the molecular level (Jia and Xue, 2009; Leffler et al., 2012).

The relative synonymous codon usage (RSCU) of the chloroplast genomes of all published Glycyrrhiza species was calculated based on all protein-coding genes (Figure 2). The results showed that the chloroplast genomes of Glycyrrhiza species contained 64 types of codons encoding 20 amino acids. Of all amino acid codons, leucine had the highest number of codons, whereas cysteine had the lowest number of codons. Thirty codons were found with an RSCU of >1, of which 29 were A/U-ending codons; 34 codons were found with an RSCU of ≤1, of which 31 were G/C-ending codons. Thus, compared with the G/C-ending codons, the chloroplast genome exhibited a greater bias toward the A/U-ending codons. The highest RSCU value was recorded for UUA and the lowest was recorded for CUC, which both encode leucine. In addition, the codon usage of the chloroplast genomes in section Glycyrrhiza was relatively close, whereas that of section Pseudoglycyrrhiza was close.

GC reflects the strength of directional mutation pressure, whereas GC3s is closely related to codon bias and used as an important basis for analyzing codon usage pattern (Shang et al., 2011; Zhao et al., 2016). GC and GC3s in the codons of these 23 chloroplast genomes were all less than 0.5, indicating that the chloroplast genomes of Glycyrrhiza species tended to use A/U bases and A/U-ending codons. Codon adaptation index values and effective number of codon values indicated a slight bias in codon usage in the Glycyrrhiza species. The frequency of optimal codons was relatively low. In addition, the hydrophobicity of the protein and the aromatic protein had little effect on codon usage bias (Table 1).

SSRs are widely distributed in the chloroplast genome and usually used as important molecular markers for species identification (Yang et al., 2011; Jiao et al., 2012). A total of 69–96 SSRs were identified in the chloroplast genomes of the Glycyrrhiza species (Figure 3A). In addition, the base composition of the repeating motifs from mononucleotide repeats to tetranucleotide repeats, which existed in all Glycyrrhiza species, had a certain A-T base preference. In these SSRs, mononucleotide repeats were the largest in number, and they were found 31–43 times in these chloroplast genomes. A/T repeats were the most common mononucleotide repeats, whereas the majority of dinucleotide repeats, trinucleotide repeats, and tetranucleotide repeats comprised AT/AT, AAT/ATT, and AAAT/ATTT, respectively. These results were consistent with A-T enrichment in complete chloroplast genomes (Qian et al., 2013). Pentanucleotide repeats were found in G. squamulosa, G. yunnanensis, Glycyrrhiza asymmetrica Hub.-Mor., Glycyrrhiza echinata L., Glycyrrhiza gontscharovii Maslenn., Glycyrrhiza lepidota Pursh, and Glycyrrhiza macedonica Boiss. & Orph., and hexanucleotide repeats only existed in G. pallidiflora, G. asymmetrica, and G. gontscharovii.

Except for SSRs, some repetitive sequences with length ≥30 bp and sequence similarity ≥90% are called long repetitive sequences, including complement (C), forward (F), palindrome (P), and reverse (R). In Glycyrrhiza species, our results revealed 61 (Glycyrrhiza astragalina Gillies ex Hook. & Arn.)–138 (G. pallidiflora) long repeats, most of which were F (35–71) and P (21–35) repeats. No C repeat was found in the chloroplast genome of G. lepidota. The length of C and R repeats was mainly within the range of 30–39 bp (Figure 3B).

In this study, the complete chloroplast genomes of all published Glycyrrhiza species were compared with the G. uralensis genome as the reference genome (Figure 4). Overall, the comparative genomic analysis revealed that the 23 Glycyrrhiza chloroplast genomes were relatively conserved. Intergenic spacers and intron regions showed more variations than protein-coding regions. Most protein-coding regions had a very high degree of conservation, and rRNA genes were highly conserved with almost no variation.

Nucleotide diversity (Pi) of shared genes and intergenic spacers of the chloroplast genomes of the 23 Glycyrrhiza species were calculated. Figure 5 shows the intergenic spacers and genes with Pi > 0. The intergenic spacers had more polymorphisms (average Pi = 0.0127) than the gene regions (average Pi = 0.0049), consistent with mVISTA analysis. Combining mVISTA and Pi, four highly variable regions were screened out, namely, accD–psaI, ndhD–ccsA, rpl36–rps8, and rrn5–trnR-ACG, for species identification and relationship analysis.

Selective pressure refers to the external forces exerted on a species during the evolutionary process, driving the species to adapt to its natural environment. In genetics, ω = Ka/Ks represents the ratio between non-synonymous mutations (Ka) and synonymous mutations (Ks). In general, synonymous mutations are considered not subject to natural selection, whereas non-synonymous mutations are influenced by natural selection. ω > 1 indicates positive selection; ω = 1 implies neutral evolution, where no selection is acting; 0<ω < 1 is considered negative or purifying selection, with small ω values indicating high negative selection pressure (Wang et al., 2010).

The results of selective pressure analyses on Glycyrrhiza chloroplast genomes showed that almost all chloroplast genes in Glycyrrhiza species underwent negative or purifying selection, which indicated these genes were relatively conserved. Only petL in G. glabra, G. pallidiflora, G. uralensis, G. inflata, Glycyrrhiza foetida Desf., Glycyrrhiza glandulifera Waldst. & Kit., Glycyrrhiza iconica Hub.-Mor., Glycyrrhiza laxissima Vassilcz., Glycyrrhiza triphylla Fisch. & C.A.Mey., G. aspera, Glycyrrhiza alaschanica Grankina, Glycyrrhiza gobica Grankina, G. gontscharovii, Glycyrrhiza zaissanica Serg., and Glycyrrhiza korshinskyi Grig. and ycf2 in G. gobica were positively selected genes, suggesting that some advantageous mutations were positively favored by selection (Figure 6). The positively selected genes may be associated with the adaptation of plants to low-light and low-temperature environments (Yang et al., 2023).

Traditionally, the Fabaceae family is classified into three subfamilies based on the different floral morphologies: Mimosoideae, Caesalpinioideae, and Papilionoideae. Some traditional classification systems, such as the Cronquist system, elevate these three subfamilies to the rank of family (Taubert, 1894). However, molecular studies indicate that Mimosoideae and Papilionoideae are generally monophyletic, but Caesalpinioideae is a highly basal and paraphyletic group. Therefore, the practice of dividing Fabaceae into three subfamilies or three families does not adhere to the principle of monophyly (Taubert, 1894). Given the keen interest of the legume taxonomic community in the subordinal classification system, the Legume Phylogeny Working Group, composed of legume experts worldwide, published a consensual phylogenetic tree of the legume family in 2013, and three schemes for subfamily division were proposed (Group, 2013). Among these schemes, the 15-subfamily scheme retains the traditional concept of the Mimosoideae subfamily, whereas the Caesalpinioideae subfamily is split into several subfamilies (with over half of them distributed in China). The 6-subfamily scheme merges the Mimosoideae subfamily into a newly defined Caesalpinioideae, and it can avoid excessive splitting of the Caesalpinioideae subfamily. Considering that the species diversity of the traditional Mimosoideae and Caesalpinioideae in the Chinese region is not high, the 6-subfamily scheme is more convenient for Chinese scholars compared with the 15-subfamily scheme. In 2017, they officially adopted the 6-subfamily scheme, publishing the comprehensive Linnaean system for the classification of the legume family. They refined the hierarchical system below the subfamily level and above the genus level based on the research results accumulated in the academic community over the years (Group, 2017).

The chloroplast genome is an important tool to explore the phylogeny of species (Hu et al., 2016). The Fabaceae family comprises 6 subfamilies and 79 tribes. Previous studies have focused on the phylogenetic relationship within the Glycyrrhiza genus (Jo et al., 2018; Jia et al., 2019; Duan et al., 2020; Jiang et al., 2020; Duan et al., 2023). Research on phylogenomics of the Fabaceae family using chloroplast genome data is limited. To date, chloroplast genomes of species from 65 tribes have been published. In this study, chloroplast genomes from species of these 65 tribes were used to construct a phylogenetic tree, including all published species of the Glycyrrhiza genus (Figure 7). The results indicated that species from the six subfamilies formed distinct clusters, consistent with the classification scheme of the six subfamilies. In the chloroplast genomes of Fabaceae plants, the absence of the IR regions was not unique to the Glycyrrhiza genus. Species lacking the IR regions were primarily distributed in various tribes, i.e., Glycyrrhizeae tribe, Adinobotryeae tribe, Wisterieae tribe, Hedysareae tribe, Caraganeae tribe, Astragaleae tribe, and Fabeae tribe in the Papilionoideae subfamily. The results of phylogenetic analysis showed that these tribes lacking the IR regions clustered together. Among them, the Glycyrrhizeae tribe clustered with the Wisterieae tribe and then with the Adinobotryeae tribe. The Hedysareae tribe, Caraganeae tribe, Astragaleae tribe, and Fabeae tribe formed a distinct cluster. In Glycyrrhiza species, the species of the section Glycyrrhiza and section Pseudoglycyrrhiza, distributed in China, respectively clustered together. The North American species G. lepidota, although had a lower content of glycyrrhizic acid (Hayashi et al., 2005) and was in another group, indicating that the groups containing glycyrrhizic acid were not monophyletic, consistent with previous studies (Duan et al., 2020). The Flora of China treats G. alaschanica, G. korshinskyi, and G. gobica as synonyms of G. uralensis; G. glandulifera as a synonym of G. glabra; and G. laxissima and G. zaissanica as synonyms of G. aspera (Li et al., 2010). However, according to the results of this study, only G. gobica and G. uralensis clustered together, G. laxissima and G. aspera clustered together, and others should be considered as independent species.

In addition to being frequently used for phylogenetic analysis and species identification, the chloroplast genome is employed for studying species divergence times (Mo et al., 2022). This study explored the divergence times of Fabaceae species using the chloroplast genome (Figure 8). The results indicated that the subfamily Detarioideae was the first to diverge, followed by the subfamily Cercidoideae, subfamily Duparquetioideae, subfamily Dialioideae, subfamily Papilionoideae, and subfamily Caesalpinioideae. In the subfamily Papilionoideae, the IRLC population was the last to differentiate. In Glycyrrhiza species, the ones without glycyrrhizic acid differentiated earlier than those with glycyrrhizic acid. This result suggested that the common ancestor of Glycyrrhiza lacked glycyrrhizic acid, consistent with previous studies (Duan et al., 2023). Medicinal groups in the Eurasian continent share a common ancestor, and their descendants result from two rapid differentiation events within the last million years. This phenomenon has led to significant morphological over-diversification and taxonomic confusion within this group (Duan et al., 2023).

A co-linear analysis of the 23 Glycyrrhiza chloroplast genomes was conducted, and the G. uralensis genome was used as the reference genome. Results showed that the entire genome sequence was a homologous region with no big indels. The 23 chloroplast genomes connected with a line, indicating that the chloroplast genomes of Glycyrrhiza species were relatively conserved, and no rearrangement and inversion occurred in gene organization (Supplementary Figure S1). Additionally, a co-linear analysis was conducted on the chloroplast genomes of other species within the Fabaceae family (Figure 9). These species were from the six subfamilies: Afzelia quanzensis Welw. (subfamily Detarioideae), Cercis chinensis Bunge (subfamily Cercidoideae), Duparquetia orchidacea Baill. (subfamily Duparquetioideae), Dialium schlechteri Harms (subfamily Dialioideae), Biancaea sappan (L.) Tod. (subfamily Caesalpinioideae), and Lotus japonicus (Regel) K.Larsen (subfamily Papilionoideae), as well as all the tribes within the subfamily Papilionoideae that lack IR regions. Thus, the chloroplast genomes of A. quanzensis, C. chinensis, D. orchidacea, D. schlechteri, and B. sappan were relatively conservative, with minimal gene rearrangements and inversions. However, when compared with species in the subfamily Papilionoideae, especially Vicia costata Ledeb., the relative positions of some homologous blocks changed, which indicated rearrangements in the sequence. Additionally, homologous blocks below the line indicated reverse complementary alignments, also known as inversion blocks, suggesting possible inversion events. Species with gene rearrangements and inversions may be prone to the structural mutations in DNA strands and accumulate constantly due to the low GC content in their chloroplast genomes under the long evolutionary history (Ding et al., 2021).