The plant materials of 13 species of the Lamiaceae used in this study were collected from the Tibet Autonomous Region of China (Tibet) (Table 1, Fig. 1) and identified by Professor Guoyue Zhong of Jiangxi University of Chinese Medicine (JUTCM). The certificate specimens were all stored in the herbarium of JUTCM.

Fig. 1.

The collected 13 Lamiaceae species.

Leaves dried with silica gel were used to extract DNA using a Plant Genomic DNA Kit (Tiangen Biochemical Technology (Beijing) Co., Ltd.). The quality of genomic DNA was detected by a spectrophotometer and 1% agarose gel electrophoresis. According to the standard genomic DNA library preparation procedure provided by Illumina, intact DNA samples with concentrations greater than 20 ng$\cdot$µL${}^{-1}$ were lysed by ultrasonic treatment, and fragmented DNA was purified and end–repaired. Libraries with an insertion size of 350 bp were prepared and then sequenced using the Illumina NovaSeq 6000 high–throughput sequencing platform (Illumina, San Diego, CA, USA). Sequencing was performed by Novgene Biotech Co., Ltd.(Nanjing, Jiangsu Province, China).

Trimmomatic v0.36 software [27] (Aachen and Institute of Bio- and Geosciences: Plant Sciences, Forschungszentrum Jülich, Leo-Brandt-Straße, Jülich, Germany) was used for quality filtering of raw data in Illumina sequencing. Single bases with a quality score lower than 20 were deleted from both ends of the sequence, as well as regions of sequences with more than three consecutive uncalled bases. All reads less than 40 bp were discarded. Filtered data were mapped to available cp genomes of the closest species in the Lamiaceae using Bowtie2 v.2.2.3 software [28] (Department of Computer ScienceUniversity College London, Gower Street, London, UK), and reads from nuclear and mitochondrial origins were excluded. The chloroplast genomes were then assembled and reconstructed using GetOeller 1.7.5 software [29] (Chinese Academy of Sciences, Kunming, Yunnan, China). CPGAVAS2 software [30] (Faculty of Technology, University of Bielefeld, Germany) (http://47.96.249.172:16019/analyzer/annotate/) was used to complete automatic annotation. Geneious 11.0.5 software [31] (Oxford, England) was used for manual correction, and OGDRAW V1.1 software [32] (Max-Planck-Institut für Molekulare Pflanzenphysiologie, Am Mühlenberg, Germany) (https://chlorobox.mpimp-golm.mpg.de/OGDraw.html) was used to map the physical structure of the chloroplast genome.

MISA software [33] was used to screen microsatellite sequences, and the parameters were set as ten repeats for single–nucleotide SSRs, five repeats for dinucleotide single sequence repeats (SSRs), four repeats for trinucleotide SSRs, and three repeats for four– and five–nucleotide SSRs.

Forward, reverse, complementary and palindrome repeat sequences were identified using REPuter software [34] (Faculty of Technology, Bielefeld, Germany) (http://bibiserv.techfak.uni-bielefeld.de/reputer/).

CodonW software 1.4.2 [35] (Paul Sharp lab, Dept of Genetics, University of Nottingham, UK) was used to identify codon usage patterns and calculate codon bias (RSCU). The CODEML program in PAML V4.973 software [36] (Department of Biology, Galton Laboratory, London, UK) was used to calculate ratio of nonsynonymous nucleotide substitution rate to synonymous nucleotide substitution rate (Ka/Ks) for each gene. Pairwise comparisons were performed, dN/dS ratio of each protein–coding gene were calculated using the yn00 program in PAML and the detailed parameters were: icode = 10, weighting = 0, commonf3x4 = 0, and other parameters in the CODEML control file were left at default settings.

Nucleotide variability (Pi) of cp genomes was assessed by alignment of sequences using MAFFT V7 software [37] (Bioinformatics Center, Institute for Chemical Research, Kyoto University Uji, Kyoto 611-0011, Japan), followed by manual adjustment of sequences using BioEdit software [38] (Borland company, Scotts Valley, California, USA). Sliding window analysis was completed using DnaSP version 5.1 software [39] (Departament de Genètica, Facultat de Biologia and Institut de Recerca de la Biodiversitat, Universitat de Barcelona, Barcelona, Spain). The step size was set to 200 bp, and the window length was set to 600 bp.

Using IRSCOPE software [40] (https://irscope.shinyapps.io/irapp/), the IR (Inverted repeat) boundaries were drawn, and the gene type and location in the border area were identified. The expansion and contraction of genes were also discovered.

The differences in cp genome sequences of the 13 Lamiaceae species evaluated in this study were compared using mVISTA software [41] (National Energy Research Scientific Computing Center Genome Sciences Department, Berkeley, USA) (http://genome.lbl.gov/vista/mvista/submit.shtml).

Selecting the published cp genomes of 84 Lamiaceae species and the 13 cp genomes from the Lamiaceae species in Tibet, the MAFFT V7 program [37] was used to compare the entire cp genomes, with manual adjustment when necessary. After being estimated in Modelfinder, the general time reversible (GTR) + r model of nucleotide substitution was used for phylogeny. Using the cp genome of L. ruthenicum (NCBI accession number NC039651.1) as the outgroup, a maximum likelihood (ML) phylogenetic tree was constructed using PhyloSuite v1.2.2 [42] (Key Laboratory of Aquaculture Disease Control, Chinese Academy of Sciences, Wuhan, China) with bootstrap values of 10,000 replicates.

Consistent with most land plants, the cp genomes of the 13 Lamiaceae species showed a conserved quadripartite structure [43] (Fig. 2), with lengths ranging from 149,081 bp to 152,312 bp, indicating that the lengths of these cp genomes were conserved [23]. The LSC region (80,908–83,416 bp) and SSC region (17,177–17,692 bp) were separated by two IRs (IRa and IRb, 24,956–25,790 bp), and the GC content was 37.83–38.49%, much lower than the AT content.

Fig. 2.

The cp genome map of 13 Lamiaceae species. Note: The black thick line indicates that the large single copy area (LSC) and small single copy area (SSC) were separated by two reverse repeat areas (IRa and IRb). The dark inner circle indicates the GC content in the chloroplast genome. Genes drawn inside the circle are transcribed clockwise, and those drawn outside the circle are transcribed counter clockwise.

The cp genome sequence of each species encoded 131–133 genes, including 88 protein–coding genes, 37 tRNA genes and 8 rRNA genes (Table 2). Taking the cp genome of A. nubigena as an example, the encoded genes can be divided into various types according to their product functions. There were 49 genes associated with photosynthesis, of which four genes (ycf1, ycf15, ycf2 and ndhB) had two copies. Fifty–nine genes were involved in self–replication, 16 of which had two copies. There were also five additional genes (Table 3). In 12 of the species, except for E. densa, there were 20 multicopy genes, and 18 genes (trnK–UUU, rps16, trnG–UCC, atpF, rpoC1, ycf3, trnL–UAA, trnV–UAC, clpP, petB, petD, rpl16, rpl2, ndhB, rps12, ndhA, trnA–UGC, trnI–GAU) contained one or two introns. Among the 18 genes, three (rps12, clpP, ycf3) contained two introns. E. densa has an extra gene (trnH–GUG) that contains an intron, which is consistent with previous studies [22].

A total of 34–57 simple repeat sequences of 6 types were detected in the 13 cp genomes, among which E. fruticosa had the highest SSR number and M. complanatum and N. hemsleyana had the lowest SSR numbers (Fig. 3). Among the six types of repeats, single–nucleotide repeats accounted for the highest proportion (61.07%), followed by dinucleotide (16.42%) and tetranucleotide (16.78%) repeats. Trinucleotide, pentanucleotide and hexanucleotide repeats accounted for 3.50%, 2.02% and 0.18%, respectively. Only one hexanucleotide repeat was detected in G. pseudogentiana. The length of single–nucleotide tandem repeats accounted for 41.88% of the total sequence length, followed by the trinucleotide sequence length (30.26%), dinucleotide sequence length (15.87%) and pentanucleotide sequence length (9.96%). The tetranucleotide sequence length accounted for 2.03%. In single–nucleotide sequences, A/T accounted for 97.58%.

Fig. 3.

Type and number of SSRs in cp genome of 13 Lamiaceae species.

SSR sequences were not evenly distributed in these cp genomes. They were distributed mostly in the LSC region, accounting for approximately 65–85%, while the remaining 15–35% were distributed in the LSC and IR regions. The SSR sequences of E. fruticosa, M. complanatum and P. betonicoides were not distributed in the IR region. SSR loci showed high polymorphism in the cp genomes of the 13 Lamiaceae plants, which provided a basis for the subsequent development of SSR molecular markers [22].

In the scattered long repeats in the cp genome, a total of 477 long repeats were detected in the 13 species analysed by REPuter software. E. fruticosa had the largest number of long repeats, while E. densa had the smallest number (Fig. 4). There were four types of long repeats, including forward, palindromic, reverse and complement. Among them, E. fruticosa and E. densa had all four types. A. nubigena, E. eriostachya, G. bifida and S. sikkimensis only had forward and palindromic types; and the other seven species had forward, palindromic and reverse types.

Fig. 4.

Long repeats in the cp genomes of 13 Lamiaceae species.

The RSCU value refers to the ratio between the actual usage frequency of a certain codon and its theoretically expected usage frequency, which is often used as an important parameter to measure codon bias. An RSCU value $>$1 indicates that the codon usage bias is strong. Such codons are high–frequency codons [44]. A total of 819 codons were found in the 13 species of the Lamiaceae. In these codons, 403 codons were high–frequency codons with RSCU values $>$1, among which 61 codons ended in A/T(U), accounting for 98.13% (Supplementary Fig. 1), indicating that high–frequency codons preferred to end in A/T(U). The lowest RSCU value was 0.30, and the highest RSCU value was 3.97.

During plant evolution, genes are affected by nonsynonymous substitutions (Ka) and synonymous substitutions (Ks). The ratio Ka/Ks can be used to indicate the evolutionary rate to determine whether gene sequences have selective pressure in the process of evolution [45]. When Ka = Ks, the gene is not under selective pressure; when Ka $>$ Ks, the gene is under positive selective pressure. The Ka/Ks ratio of genes in the cp genomes of the 13 species was 0–0.58 (Fig. 5), indicating that most genes are subjected to purification selective pressure. The Ka/Ks ratios of 8 photosynthesis–related genes (petG, petL, petN, psaC, psbF, psbI, psbM, psbN), a ribosomal large gene (rpl36) and a small subunit gene (rps7) were all 0. The Ks values of 3 genes (psaI, rpl23 and ycf15) were also 0.

Fig. 5.

Analysis of gene selection pressure of 13 Lamiaceae species. Note: The x–axis represents the gene, and the y–axis represents the value of Ka/Ks.

Nucleotide diversity (Pi) values were calculated to identify divergence hotspots. With Pi values ranging from 0.000 to 0.230, a total of 754 variation sites were detected (Fig. 6). In addition, there were 8 variation sites with Pi values $>$0.12. Between 4810–5400, 7001–7600, 8001–8600, 71,801–72,400, 133,601–134,200, 135,401–136,000, 135,601–136,200, and 135,801–136,400, four genes were distributed in the LSC region, two in the SSC region and two in the IRb region. Higher levels of genetic variation were detected in the LSC and SSC regions, suggesting that rapid nucleotide substitution in the Lamiaceae species may be involved, which may provide further consideration for the identification of the Lamiaceae species and the construction of phylogeny.

Fig. 6.

Nucleotide diversity analysis of 13 Lamiaceae species.

The expansions and contractions of the IR regions often result in genome size variations among various plant lineages. Therefore, the analysis of the boundary changes of the cp genome can reflect the evolutionary process [46]. In the cp genomes of N. laevigata, rpl22/trnH and trnH/psbA were located on the two sides of the LSC/IRb and LSC/IRa boundaries, respectively (Fig. 7). In the cp genomes of other species, rps19/rpl2 and rpl2/trnH were discovered on the two sides of the two boundaries. The expansion of the IR boundary allows it to enter the coding regions of nearby genes, which will lead to the formation of pseudogenes in the other IR copy region due to the inverted repeat characteristic of IR [23]. In this study, the LSC/IRb boundary entered the coding region of the rps19 gene, and the SSC/IRa boundary entered the coding region of the ycf1 gene; therefore, both ycf1 and rps19 became pseudogenes. The IR boundaries of the cp genomes of the other 12 Lamiaceae species except N. laevigata were essentially conserved, suggesting that the evolutionary direction of N. laevigata may be different from those of the other 12 species.

Fig. 7.

Comparative analysis of cp genome boundary regions (LSC, SSC, IR) in 13 Lamiaceae species. Note: The coloured squares indicate that the large single copy (LSC) area and the small single copy (SSC) area are separated by two reverse repeat areas (IRa and IRb). JLB (LSC and IRb boundary), JSB (IRb and SSC boundary), JSA (SSC and IRa boundary), JLA (IRa and lSC boundary).

The order of genes in the 13 cp genomes were essentially identical, no rearrangement or inversion were detected among these genomes (Supplementary Fig. 2).

To investigate intergeneric differences among cp genomes, the identity percentage of these species was plotted using the mVISTA program, taking the cp genome of A. nubigena as a reference. A high degree of similarity was detected among the 13 cp genomes. More variations were detected in the LSC/SSC region than in the IR region. The variation in the non–coding region was significantly higher than that in the coding region, indicating that the 13 cp genomes were conserved (Supplementary Fig. 3).

In this study, we constructed an ML phylogenetic tree (Fig. 8) and analysed the genetic relationship. During sequence alignment, the alignment length was 18,755,965 bp. The results showed that all species comprise eight clades, which consisted of the Nepetoideae subfamily, Lamioideae subfamily, Ajugoideae subfamily, Scutellarioideae subfamily, Premnoideae subfamily, Viticoideae subfamily, Tectona genus, and Callicarpa genus (Fig. 8). The major clade made from the Nepetoideae subfamily was composed of species of 12 genera, e.g., the Salvia, Melissa and Mentha genera; the Viticoideae subfamily was composed of the Vitex genus; the Lamioideae subfamily was composed of species of 13 genera, e.g., the Stenogyne, Phyllostegia and Haplostachys genera; the Scutellarioideae subfamily was composed of the Scutellaria and Holmskioldia genera; the Ajugoideae subfamily was composed of species of four genera, e.g., the Ajuga, Teucrium and Caryopterisgenera; and the Premnoideae subfamily was composed of the Gmelina and Premna genera. These results are consistent with the findings of Li [47]. Except for three branch nodes, the bootstrap value of each node was greater than 70, indicating that the constructed ML tree was reliable. In addition, 97 species could also be divided into eight tribes (Fig. 8). Except for Lamium galeobdolon (NC036972.1) which was not clustered with the Lamium genus, the phylogenetic relationships of other species in the ML tree at the tribe level were consistent with those in Li’s study [47].

Fig. 8.

The ML phylogenetic tree of cp genomes of 97 Lamiaceae species. Note: Coloured squares represent subfamilies and genera, coloured horizontal lines represent families. Brown fonts represent sister groups and orange fonts represent the complex group.