1. Introduction

Salvia L., which includes about 1000 species, is the largest genus of Lamiaceae (Hu et al., 2018). Salvia is widely distributed in the world, mainly in Mesoamerica, South America, East Asia, the Mediterranean region, and south-western Asia (Walker and Sytsma, 2006; Wei et al., 2015). Most species of Salvia have important medicinal, ornamental, and economic values (e.g., Salvia miltiorrhiza Bung., Salvia splendens Ker Gawl., and Salvia hispanica L.) (Ge et al., 2014; Ali et al., 2012; Wang, 2010). Salvia is well known for its highly specialized flower structure and highly diverse staminal morphology which forms the staminal lever mechanism (Claßen-Bockhoff et al., 2008).

Recent molecular phylogenetic studies have suggested that traditionally defined Salvia is not monophyletic (Drew et al., 2017; Fragoso-Martínez et al., 2018; Hu et al., 2018; Jenks et al., 2012; Kriebel et al., 2019; Li et al., 2013; Walker et al., 2004; Will and Claßen-Bockhoff 2014; Will and Classen-Bockhoff, 2017). Some researchers (Drew et al., 2017; Hu et al., 2018) have proposed that Salvia should be expanded to include five additional genera (Dorystaechas Boiss. & Heldr. ex Benth., Meriandra Benth., Perovskia Karel., Rosmarinus L. and Zhumeria Rech. & Wendelbo). Will and Classen-Bockhoff (2017) have recommended dividing the genus of Salvia into six genera. Here, we adopt a broad definition of Salvia following the justification of Drew et al. (2017).

Recently, sequencing technology and analytical methods have developed rapidly, and some plastid DNA fragments, such as psbA–trnH, matK, rbcL, trnL–trnF, ycf1, and nuclear ribosomal Internal Transcribed Spacers (nrITS), have been utilized to study the phylogeny of Salvia (Chen et al., 2009). Drew et al. (2017) utilized two low-copy nuclear gene regions (PPR–AT3G09060 and GBSSI) and four plastid markers (psbA–trnH, trnL–trnF, ycf1, and ycf1–rps15) to infer the phylogenetic relationships of Salvia. Their results primarily confirmed that Salvia in the traditional sense was non-monophyletic, embedding five additional genera (Dorystaechas, Meriandra, Perovskia, Rosmarinus and Zhumeria) based on plastid and nuclear DNA data. Using nrITS and plastid marker rpl32–trnL, Will and Classen-Bockhoff (2017) found that Salvia species can be divided into four lineages, with the above-mentioned five genera clustered within them. Meanwhile, Fragoso-Martínez et al. (2017) used an anchored hybrid enrichment method to resolve relationships among subgenus Calosphace. In addition, Hu et al. (2018) used two nuclear markers (nrITS and ETS) and four plastid DNA markers (psbA–trnH, ycf1–rps15, trnL–trnF and rbcL) to reconstruct the phylogenetic analysis of East Asian Salvia. Using the same approach, Kriebel et al. (2019) provided a well-resolved backbone phylogeny of Salvia. However, previous analyses have all been based on DNA markers, whereas few genomic data sets have been used to construct a phylogeny of Salvia.

Altogether, Salvia is divided into 11 subgenera based on molecular and morphological data (Drew et al., 2017; Hu et al., 2018; Kriebel et al., 2019; Will and Classen-Bockhoff, 2017). East Asia is one of the distribution centers of Salvia, with approx. 100 species in total, more than 80% of which are in China (Hu et al., 2018). Based on phylogenetic analysis focused on Salvia in East Asia, Hu et al. (2018) established a subgenus, subg. Glutinaria, to accommodate all species of Salvia native to East Asia. There are also many important medical plants in this subgenus, such as S. miltiorrhiza, S. cavaleriei Levl and S. yunnanensis C. H. Wright (Qian et al., 2002; Wang, 2010; Zhang and Li, 1994). Few early molecular phylogenetic studies have focused on species currently placed in Salvia subg. Glutinaria. Takano and Okada (2011) used three DNA markers (rbcL, trnL–trnF and nrITS) to report that 11 Japanese Salvia species were monophyletic. Using four DNA markers (psbA–trnH, rbcL, matK and nrITS), Li et al. (2013) demonstrated that 41 species of Salvia from China and Japan formed a clade. Hu et al. (2018) used 93 taxa to reconstruct the phylogenetic relationships of Salvia subg. Glutinaria and recognized eight sections within this subgenus. However, interspecific relationships of Salvia subg. Glutinaria remain unresolved.

The chloroplast genome, or plastome, has long been used in plant phylogenetics. Several recent studies have found that using plastomes can powerfully resolve the phylogenetic relationship of plants at different taxonomic levels (Lin et al., 2012; Wang et al., 2016; Yang et al., 2013a, 2013b). Plastomes are relatively stable, easy to obtain, and contain a large amount of information with a suitable mutation rate (Tian and Li, 2002). The plastome contains four extremely evolutionarily conserved parts: a pair of inverted repeat regions (IRa and IRb), a large single-copy region (LSC), and a small single-copy region (SSC) in most land plants. Structural variation and the presence/absence of genes in these genomes can also be used for phylogenetic inference.

In total, 53 Lamiaceae plastomes have been sequenced and are available for public use (NCBI GenBank database, accessed on 31 Jan 2020). The structure and organization of Lamiaceae plastomes is conserved and stable, and consists of the quadripartite plastome structure of most angiosperms. Lamiaceae plastomes possess between 80 and 81 unique protein-coding genes, range in length from 149, 736 to 155, 293 bp, and have a GC content between 37.8 and 38.7% (Lindqvist and Albert, 2002). However, there are only seven plastomes from Salvia publicly available. Hence, the plastome characterization of Salvia is still far from adequate and few phylogenetic trees of Salvia have been reconstructed using plastomes in previous phylogenetic relationship of Salvia, particularly in the East Asian subg. Glutinaria.

Here, we analyzed Salvia phylogeny using plastome sequence data. The main objectives of this research are to (1) investigate structural and compositional variation of plastomes in East Asian Salvia species; (2) estimate the feasibility of using plastome data for phylogenetic reconstruction of the Salvia subg. Glutinaria; and (3) find the optimal DNA regions to better construct the phylogenetic tree of Salvia at large.

2. Materials and methods 2.1. Plant materials

We sampled a total of twenty-one Salvia accessions from nineteen different species. Fresh green leaves of fourteen Salvia species and two Melissa L. were acquired from the field and sequenced. Two species of Melissa were selected as outgroup based on previous phylogenetic studies (Hu et al., 2018; Kriebel et al., 2019). Voucher specimens were deposited at the herbarium of the Kunming Institute of the Botany, Chinese Academy of Sciences (KUN) (for details see Table 1). Five Salvia plastomes (S. miltiorrhiza, S. rosmarinus Spenn., S. officinalis L., S. chanryoenica Nakai and S. bulleyana Diels) were downloaded from GenBank for analyses. Previously reported data for two Salvia species were not used because these species may have been misidentified (Salvia japonica Thunb (He et al., 2017).) or were re-sequenced in this study (Salvia przewalskii Maxim.). This data set represents all four subgenera of Salvia and each of the four sections of subg. Glutinaria.

2.2. DNA extraction, sequencing, and assembly

Chloroplast DNA was extracted following the CTAB method (Doyle, 1987). After purification, we fragmented 5 mg DNA to build short-insert libraries, following the manufacturer's instructions (Illumina). We used tags to index DNA from different individuals and pooled samples before sequencing in Illumina's Genomic Analyzer. Low-quality reads were discarded using the NGS QC Toolkit (Patel and Jain, 2012). Sequence reads were filtered from the raw sequence reads through comparison to published sequences of S. miltiorrhiza (NC_020431) (Qian et al., 2013). The remaining sequences were assembled into complete plastomes de novo using GetOrganelle (Jin et al., 2018). Complete plastomes were visualized in Bandage v.8.1. (Wick et al., 2015). Finally, the order of the complete plastomes was determined by comparison to the reference genome S. miltiorrhiza.

2.3. Chloroplast genome annotation

Plastomes were annotated using Plastid Genome Annotator (PGA) genomes and manually revised (Qu et al., 2019). The tRNA boundaries were then identified with tRNAscan-SE (v.1.3.1) (Lowe and Chan, 2016). Salvia plastome gene maps were drawn with the OrganellarGenomeDRAW tool (OGDRAW) (Wyman et al., 2004).

2.4. Comparative chloroplast genomic analyses

To investigate interspecies variation, whole plastome alignment of the 21 Salvia species and two species of Melissa was performed using the mVISTA program in Shuffle-LAGAN mode (A Frazer et al., 2004), with S. miltiorrhiza as the reference. IR expansion/contraction of these plastomes was also checked using online IRscope (Amiryousefi et al., 2018). The nucleotide diversity (Pi) of plastomes was estimated by DnaSP v.5.10 (Librado and Rozas, 2009) with step size set at 200 bp and window length set at 600 bp.

2.5. Plastid phylogenomic analyses

In this study, 23 plastome sequences, including those of two outgroups from subtribe Melissinae, were aligned for phylogenomic analysis. To examine the performance of different plastome regions in phylogenetic resolution, 10 matrices were created for phylogenetic analyses (Table 2). The first five matrices of plastid sequence data were as follows: (1) complete plastid DNA sequences containing only one IR (Plastomes); (2) protein-coding exons (CR); (3) introns and spacers (IS); (4) the large single-copy region and the small single-copy region (SCR); and (5) the inverted repeat region (IR). We created five additional matrices by eliminating rapidly evolving sites from each of the preceding data sets using Gblocks (Castresana, 2000) with default parameters. Plastome sequences were aligned using the Multiple Sequence Alignment Program (MAFFT v.5) and manually adjusted when necessary (Katoh et al., 2005). Ambiguously aligned regions were excluded.

Maximum likelihood (ML) analyses were conducted using RAxML v.8.2.9 (Stamatakis et al., 2005), under the general time-reversible model of nucleotide substitution with a gamma distribution of rate variation among sites (GTRGAMMA) with 1000 replicates and other default parameters.

Maximum parsimony (MP) analysis were performed with PAUP∗4.0b10 (Swofford, 2002). A 1000–replicate heuristic search was applied with random addition, and tree bisection-reconnection (TBR) branch swapping with the MUL-trees was used to search for MP trees. Parameter settings remained unchanged.

Bayesian inference (BI) analysis was carried out in MrBayes v.3.2.6 (Ronquist et al., 2012) in PhyloSuite (Zhang et al., 2020), using the ModelFinder (Kalyaanamoorthy et al., 2017) to select the best model based on Akaike information criterion (AIC). GTR + I + F was found to be the best-fit model for IR Regions and IR Regions (using Gblocks), and GTR + I + G + F for the other data sets (Table 2). The Markov chain Monte Carlo (MCMC) algorithm was run for 2, 000, 000 generations. One tree was sampled every 1000 generations. The first 25% of generations were discarded as burn-in. To construct a consensus tree, we used trees that remained after reaching a stagnant state, i.e., when the average standard deviation of the splitting frequency was < 0.01. Reconstructed trees were visualized in FigTree v.1.4.2 (Rambaut, 2012).

3. Results 3.1. Features of the plastomes

All newly sequenced Salvia and Melissa plastomes possess the typical quadripartite structure of angiosperms, with two single-copy regions (LSC and SSC) and a pair of IR regions (IRa and IRb) (Fig. 1). The length of the plastomes in Salvia ranged from 150, 829 bp (S. yunnanensis) to 154, 984 bp (Salvia castanea Diels f. castanea). The LSC regions ranged from 82, 078 bp in Salvia tiliifolia Vahl to 83, 879 bp in S. sp.; the SSC regions ranged from 17, 500 bp in S. officinalis to 19, 031 bp in Salvia mekongensis Stib.; and the IR regions ranged from 24, 978 bp in S. mekongensis to 27, 265 bp in S. castanea f. castanea. The GC content ranged from 37.9 to 38.1%, with the lowest GC content observed in Salvia subpalmatinervis Stib. and S. mekongensis (37.9%) and the highest GC content observed in S. yunnanensis, S. rosmarinus 1 and S. castanea f. castanea (38.1%). The plastomes of the two Melissa species were 151, 885 bp and 151, 985 bp, and the GC content was 38.0% and 38.1%, respectively (Table 3). All plastomes contained the same 114 unique genes, including 80 protein-coding genes, 30 tRNA genes and 4 rRNA genes (Fig. 1, Table 3).

3.2. Comparative genomic analyses

In all plastomes, the coding region of the rps19 gene is located at the LSC/IRb boundary, which is stable at 42–43 bp, resulting in a pseudogene fragment (Ψrps19) at the LSC/IRa border. In all plastomes except Salvia mekongensis and S. rosmarinus 2, the ndhF gene spanned the IRb/SSC boundary and ranged from 32 to 64 bp. In all plastomes except S. rosmarinus 1, the coding region of the ycf1 gene spans the IRa/SSC region with a length of 402–1167 bp, which results in a pseudogene (Ψycf1) at the IRa/SSC boundary (Fig. 2). In all plastomes, sequences were more conserved in the IR regions than in the LSC and SSC regions; furthermore, sequence variation in non-coding regions was more divergent than in coding regions. Plastome sequences varied most in the intergenic region (e.g., rps16–trnQ–UUG, atpH–atpI, petN–psbM, psaA–ycf3, accD–psaI, ycf4–cemA, petA–psbJ, rps12–trnV–GAC, rpl32–trnL) (Fig. 3).

The average nucleotide diversity (Pi) of these plastomes was 0.009383 (Fig. 4). Pi values were mostly greater than 0.008 in the LSC and SSC regions. Moreover, two regions (ycf1 and rpl32–trnL–UAG) with relatively high variability (> 0.03) were located in SSC regions. The highest divergence level in SSC regions was 0.04112 (ycf1) (Fig. 4). In the LSC region, one gene (rpl16) and one intergenic region (trnL–UUU–rps16) had a relatively high variability. The IR regions also had one highly variable locus, rps12–trnV–GAC.

3.3. Phylogenetic analysis

Ten matrices from 23 complete plastomes were used for phylogenetic analyses (Table 2). The complete plastid DNA sequences only containing one IR (Plastomes) data set was the longest alignment with a length of 137, 226 bp and the shortest alignment was the IR region with Gblocks (IR (Gblocks)) with a length of 24, 540 bp. The constant sites from large to small were plastomes data set with 113, 108 bp to IR (Gblocks) dataset with 23, 993 bp. The parsimony-informative varied from the plastomes dataset with 3644 bp to IR/IR (Gblocks) dataset with 119 bp.

ML, MP and BI analyses resulted in a very similar phylogenetic relationship. The ML topology was therefore selected for discussion, and ML bootstrap (MLBS), MP bootstrap (MPBS), and posterior probabilities (PP) values are given above branches (Fig. 5, Fig. 6, Fig. 7). In this study, we defined well supported as MLBS ≥85%, MPBS ≥85%, and PP ≥ 0.99; moderately supported as 70% < MLBS ≤85%, 70% < MPBS ≤85%, 0.90 < PP ≤ 0.99; and weakly-supported as MLBS < 70%, MPBS < 70%, and PP < 0.90.

In all matrices, Salvia formed a monophyly with strong support (Fig. 6, Fig. 7, MLBS, MPBS, PP = 100%, 100%, 1.00, respectively). Three clades were recovered (Fig. 5). The first clade, which is strongly supported (100%, 100%, 1.00), consists of two subgenera, S. subg. Rosmarinus represented by S. rosmarinus and S. subg. Salvia represented by S. officinalis. The second clade is the Salvia subg. Calosphace represented by Salvia tillifolia. The third clade is the subgenus Glutinaria with strongly supported values (100%, 100%, 1.00), consisting of 16 species (Fig. 5, Fig. 6, Fig. 7). The topologies of ten phylogenetic trees obtained by different matrices are highly similar, with only the support values varying in different trees (Fig. 6, Fig. 7). The phylogenetic tree based on the matrix of LSC and SSC region with Gblocks (SCR(Gblocks)) had the highest support values. All relationships were well supported (100%, 96–100%, 1.00) with the exception of that between Salvia umbratica Hance and S. castanea f. castanea, which was moderately supported (Fig. 7, 81%, 71%, 0.99). The poorest resolved phylogenetic trees were constructed using IR and IR (Gblocks); some branches had weak-support, and only five relationships were strongly supported (100%, 100%, 1.00) (Fig. 6, Fig. 7).

4. Discussion 4.1. Genome organization

All Salvia plastomes examined here possess the highly conserved quadripartite structure typical of most angiosperms. Salvia plastomes shared 114 unique genes that include 80 protein-coding genes, 30 tRNAs, and 4 rRNAs (Lukas and Novak, 2013), which is consistent with the plastomes of two Melissa (Lamiaceae, Nepetoideae) species and Lamiaceae plastomes such as Lavandula angustifolia Mill. (Lamiaceae, Nepetoideae) (Ma, 2018). Salvia plastome gene order is consistent with that of species throughout Lamiaceae (He et al., 2016; Ma, 2018). Moreover, the 21 Salvia plastomes had similar GC content. Furthermore, the GC content of Saliva was to those of other genus within the Lamiaceae family, even as well as to that of Boea hygrometrica (Bunge) R. Br. in Gesneriaceae family from the same order (Lamiales) (Zhang et al., 2012).

In angiosperms, gene loss-and-gain events and structural rearrangement of plastomes have been found in some species or genera (Wicke et al., 2011), but these events have not been observed in Salvia species. These results suggest that the plastome structure and content of Salvia are largely conservative, even in the order Lamiales.

IR contraction-expansion is a pervasive and primary cause of plastome length variation (Raubeson et al., 2007; Wang et al., 2008). The IR/LSC junctions of the 21 Salvia plastomes were highly conserved (Fig. 2). In all Salvia plastomes, the Ψrps19 and trnH genes were present at the IRs/LSC boundaries. The IRb/SSC junctions of 19 Salvia plastomes were conserved, and the ndhF gene spanned the IRb/SSC junctions. IR contraction was observed in two species (S. rosmarinus 2 and S. mekongensis). The IRa/SSC junctions of the 20 Salvia plastomes were conserved and the ycf1 gene spanned the IRa/SSC border, IR contraction was only observed found in S. rosmarinus 1. The SSC/IR boundaries of the two S. rosmarinus accessions differed, which may be due to intraspecies variation. The largest difference in Salvia plastome size was 4155 bp (S. yunnanensis versus S. castanea f. castanea); in the IR region the largest difference was 2287 bp (S. mekongensis versus S. castanea f. castanea); in the SSC region, the largest difference was 1531 bp (S. officinalis versus S. mekongensis); and in the LSC region, the largest difference was 1801 bp (S. tiliifolia versus S. sp.). These differences in Salvia plastid genome lengths can be largely explained by differences in IR length, which was the region of greatest variation. Salvia plastome sequence variation was highly consistent (Fig. 3); IR regions were mostly conserved and non-coding regions were more highly divergent than coding regions. These patterns suggest that the content and structure of Salvia plastomes are highly conserved. The Pi values of these sites (ycf1 and rpl32–trnL–UAG) exceed 0.03, with one Pi value, ycf1, higher than 0.04. We believe that these sites with large variation can be used to identify species in Salvia.

4.2. Phylogenetic implications

In this study, we used ten data matrices to generate high resolution phylogenetic trees to evaluate the effectiveness of different plastome sequences in phylogenetic inference of Salvia, and in particular the Salvia subg. Glutinaria. All phylogenetic trees have highly similar topologies and are generally well supported (Fig. 5, Fig. 6, Fig. 7). Recently, East Asian Salvia have been treated as a subgenus, S. subg. Glutinaria (Hu et al., 2018). Our results comfirm that S. subg. Glutinaria is a monophyletic group with strong support in all trees (Fig. 5, Fig. 6, Fig. 7), corroborating previous work based on two nuclear ribosomal spacers and four plastid markers (Hu et al., 2018). Previous studies have been unable to determine the interspecific relationships between species within the subg. Glutinaria (e.g., S. cavaleriei, S. miltiorrhiza and S. yunnanensis) (Hu et al., 2018). In our study, S. cavaleriei, S. miltiorrhiza and S. yunnanensis form a strongly supported clade with Salvia plectranthoides, in which S. miltiorrhiza occupies the basal position, followed by S. cavaleriei, which is sister to S. yunnanensis and S. plectranthoides in all trees. The phylogenetic trees generated in our study are more highly resolved than those generated by nuclear genes in previous studies and differ in several ways (Hu et al., 2018). 1) Salvia cyclostegia is not in sect. Eurysphace (which includes S. przewalskii, S. subpalmatinervis, S. bulleyana, S. castanea, S. mairei, and Salvia flava). 2) S. cavaleriei is embedded in sect. Drymosphace (which includes S. miltiorrhiza, S. yunnanensis, and S. plectranthoides), and 3) S. umbratica is embedded in sect. Eurysphace. Many factors may contribute to these discrepancies, including incomplete lineage sorting, chloroplast capture, and hybridization (Albaladejo et al., 2005; Deng et al., 2015; Xiang et al., 2013). S. castanea Diels f. tomentosa and S. castanea f. castanea do not form a monophyletic group, which is consistent with previous phylogenies based on nuclear ribosomal spacers (Hu et al., 2018). This implies that these two sub-species underwent cryptic speciation rather than hybridization, as the two forms are morphologically similar, but genetically distinct. Previous molecular phylogenetic studies have used both nuclear genome fragments and plastid fragments; however, clades were not highly supported and relationships were not well resolved (Drew et al., 2017; Hu et al., 2018; Will and Classen-Bockhoff, 2017). The main reason for low support of previous studies is the lack of a sufficient number of informative sites to construct phylogenetic tree. The use of plastid phylogenomics in the present study has provided better resolution of the phylogenetic tree of Salvia subg. Glutinaria.

Phylogenetic analysis using complete plastome sequences contained over 100 times more simplified informative sites than those of previous analysis based on DNA fragments. Consequently, each node of our phylogenetic trees is highly supported. Specifically, four matrices have high support (Fig. 6, Fig. 7) (Plastomes, SCR, Plastomes (Gblocks), and SCR (Gblocks)), and only one or two nodes do not have strongly supported values (100%, 100%, 1.00). The other six matrices all have more than two nodes weakly supported. These analyses indicate that plastome sequences are conducive to resolving phylogenetic relationships of Salvia subg. Glutinaria. Saturation of sequence mutations can impact support values of phylogenetic trees (Knight and Mindell, 1993; Simon et al., 1994). From the constructed phylogenetic trees, we found that resolving phylogenetic relationships is not a matter of having the most variable sites possible. Some saturated variation will affect the resolution of systematic relationships, a phenomenon that now seems to be evident in Salvia. For example, although the plastome matrix had the greatest number of informational sites (3644 bp), two relationships are still only moderately supported. Phylogenetic trees based on the SCR matrix, which had a lower number of informational sites (3505 bp), have higher support values, and only the relationship between S. castanea f. castanea and S. umbratica was weakly-supported (Fig. 6, 61%, 50%, 0.86). Gblocks reduced the number of informative sites in the SCR (Gblocks) data set (3425 bp); however, the clade that contains S. castanea f. castanea and S. umbratica became moderately supported (Fig. 7, 81%, 71%, 0.99). Conversely, there can be too few informative sites. Most clades based on IR and IR (Gblocks) (119 bp) were weakly supported because there were insufficient numbers of informative sites (Table 2). Previous studies have suggested that the selection of suitable information sites is important to phylogenetic reconstruction. For example, in Cymbidium Sw., using introns and interval sequences of plastomes can effectively enhance the resolution of the tree (Yang et al., 2013a). Here the phylogenetic tree constructed with the SCR (Gblocks) matrix resolved relationships within the Salvia subg. Glutinaria very well (Fig. 6, Fig. 7). Therefore, our findings show that using sequences with suitable mutations and informative sites is helpful for resolving phylogenetic relationships. In fact, because saturation of sequence mutations may affect phylogenetic reconstruction, deleting saturation mutation information appropriately is more conducive for phylogenetic reconstruction.

Although all phylogenetic relationships cannot be resolved by using complete plastomes (Petersen et al., 2011; Richard et al., 2009), our results suggest that plastome-wide analysis will provide resolution for some disputed relationships. Further research is needed that combines nuclear genes and more plastomes to elucidate the phylogeny and evolution of Salvia.

5. Conclusion

The present study has shown that the plastome structure of Salvia is largely conserved. The plastome of 21 Salvia has a typical quadripartite structure with 114 coding genes. The length of plastomes in Salvia is 150, 829–154, 984 bp, GC content is 37.9–38.1%, and the variation of the non-coding regions is higher than coding regions. These characteristics are similar to those of the plastomes of other Lamiaceae species. We found that using complete plastomes can better resolve the phylogenetic relationships of Salvia subg. Glutinaria. In addition, our results show that using a plastome data set that includes the large single-copy and small single-copy regions with rapidly evolving sites removed is the most suitable approach to reconstructing the phylogeny of Salvia. Our demonstration that comparative plastomic analysis provides insight into the phylogeny of Salvia, the largest genus of Lamiaceae, indicating that using suitable mutations with informative sites are helpful in resolving phylogenetic problems.

Author contributions

Study conception: HW, H-TL, D-ZL; acquisition of molecular data: HW, H-TL; plastomic and phylogenetic analysis: HW, H-TL, P-FM; drafting of manuscript and pictures: HW, H-TL, P-FM, G-XH, D-ZL; critical revisions: HW, P-FM, D-ZL. All authors have read, commented and approved the final manuscript.

Declaration of Competing Interest

The authors declared that they have no conflict of interest.

Acknowledgements

We are grateful to the Germplasm Bank of Wild Species at the Kunming Institute of Botany (KIB) and Molecular Biology Experiment Center, Germplasm Bank of Wild Species for facilitating this study. This research was supported by the Large-scale Scientific Facilities of the Chinese Academy of Sciences (Grant No: 2017–LSFGBOWS–02). We thank Dr. Chun-Lei Xiang for reviewing the first draft for us.