Fresh leaves of S. cusia from multiple individuals were collected from Guangxi Medicinal Plant Garden, Nanning, Guangxi, China and stored at 4°C for cp DNA isolation. S. cusia is not an endangered or protected species, and no specific permissions were required for collection. A voucher specimen was deposited at the Institute of Medicinal Plant Development, Beijing, China. Total DNA was extracted using a plant genomic DNA kit (Tiangen Biotech, Beijing, Co., Ltd.). DNA purity was evaluated using 1.0% agarose gel, and DNA concentration was measured using a Nanodrop spectrophotometer 2000 (Thermo Fisher Scientific, United States).

Further steps of library preparation were performed using a TruSeq DNA Sample Prep Kit (Illumina, Inc., United States) according to the manufacturer’s instructions. The DNA was sheared to yield approximately 500 bp-long fragments for paired-end library construction. The library was sequenced on Illumina HiSeq 3000 (Illumina Inc.). In total, 9,912,889 paired-end reads (2 × 150 bp) were obtained.

We first downloaded 1,006 plastid genomes from GenBank in February 2016. These plastid genome sequences were used to search against Illumina paired-end reads using BLASTN with an E-value cutoff of 1e-5. The genome sequence of A. paniculata (Accession number: NC_022451) had the highest overall sequence similarity to the reads and was used as a reference for the downstream genome assembly.

AbySS (v1.5.2) (Simpson et al., 2009) and CLC Genomics Workbench (v7) was used for the de novo genome assembly. Using Gepard (Krumsiek et al., 2007), we identified 7 contigs from the assembly of contigs of AbySS and CLC Genomics Workbench, respectively, that nearly spanned the entire cp genome. All the 14 contigs were assembled using Seqman module of DNASTAR (v11.0). Then, we obtained only one sequence corresponding to the large single-copy (LSC), IR, and small single-copy (SSC) regions of S. cusia. Regions corresponding to the IR/SSC and IR/LSC boundaries were confirmed via direct PCR amplifications.

PCR amplifications were performed using the sequence-specific primers (Supplementary Table S1) under the following conditions: pre-denaturation at 94°C for 2 min, 35 cycles of amplification at 94°C for 30 s, 55°C for 30 s, and 72°C for 30 s, followed by a final extension at 72°C for 2 min. The PCR reaction mixture contained 25 μL of Taq MasterMix (2 × ), 2 μL of forward primer (10 μM), 2 μL of reverse primer (10 μM), and purified cp DNA ( < 1 μg). RNase-free water was added to a final reaction volume of 50 μL.

The CpGAVAS web service (Liu et al., 2012) was used to annotate the S. cusia cp genome. Cutoffs for the E-values of BLASTN and BLASTX were 1e-10. The number of top hits to be included in the reference gene sets for annotation after the pre-filtering step was 10. Meanwhile, tRNA genes were identified using tRNAscan-SE (Lowe and Eddy, 1997) and ARAGORN (Laslett and Canback, 2004). Manual corrections on the positions of the start and stop codons, and for the intron/exon boundaries were performed based on the entries in the cp genome database (Cui et al., 2006) using the Apollo program (Lee et al., 2009). Moreover, the circular cp genome map of S. cusia was drawn using OrganellarGenomeDRAW (Lohse et al., 2013). Furthermore, codon usage and GC content (that is, the percentage of Guanines and Cytosines) were analyzed using the Cusp and Compseq programs provided by EMBOSS (Rice et al., 2000). Final genome assembly and genome annotation results were deposited in the GenBank (accession number: MG874806).

Simple sequence repeats (SSRs) were detected using MISA Perl Script available at http://pgrc.ipk-gatersleben.de/misa/ with the following thresholds: 8 repeat units for mononucleotide SSRs, 4 repeat units for di- and trinucleotide repeat SSRs, and 3 repeat units for tetra-, penta-, and hexanucleotide repeat SSRs. Tandem repeats were analyzed using Tandem Repeats Finder (Benson, 1999) with parameter settings of two for matches and seven for mismatches and indels. The minimum alignment score and maximum period size were set at 50 and 500, respectively. All the identified repeats were manually verified and nested, or redundant results were removed. REPuter (Kurtz et al., 2001) was employed to identify the IRs in S. cusia by forward versus reverse complement (palindromic) alignment. The minimal repeat size was set at 30 bp, and the cutoff for similarities among the repeat units was set at 90%.

Conserved sequences were identified between the cp genomes of Astragalus membranaceus and those of A. paniculata (NC_022451.2), R. breedlovei (KP300014.1), Tanaecium tetragonolobum (NC_027955.1), Dorcoceras hygrometricum (NC_016468.1), Salvia miltiorrhiza (NC_020431.1), Olea europaea (NC_013707.2), Sesamum indicum (NC_016433.2), and Scrophularia takesimensis (NC_026202.1) by using BLASTN with an E-value cutoff of 1e-10. The homologous regions and gene annotations were visualized using a web-based genome synteny viewer GSV (Revanna et al., 2011).

A total of 31 species were examined for expansion or contraction of the IR, representing 29 species from Lamiales and 2 species from Actinidiaceae. The annotations of these cp genomes were downloaded from RefSeq database.

A total of 28 complete cp DNA sequences belonging to the Lamiales order were obtained from RefSeq database. For the phylogenetic analysis, 65 protein sequences were shared among all these 31 species, and S. cusia was aligned using the CLUSTALW2 (v2.0.12) program. The 65 proteins included ATPA, ATPB, ATPE, ATPF, ATPH, ATPI, CCSA, CEMA, MATK, NDHA, NDHB, NDHC, NDHE, NDHF, NDHG, NDHH, NDHI, NDHJ, NDHK, PETA, PETD, PETG, PETL, PETN, PSAA, PSAB, PSAC, PSAI, PSAJ, PSBA, PSBC, PSBD, PSBE, PSBF, PSBH, PSBJ, PSBK, PSBL, PSBM, PSBN, PSBT, PSBZ, RBCL, RPL14, RPL2, RPL20, RPL22, RPL23, RPL32, RPL33, RPL36, RPOB, RPOC1, RPOC2, RPS11, RPS14, RPS15, RPS18, RPS2, RPS3, RPS7, RPS8, YCF2, YCF3, and YCF4 (Supplementary File 1). The alignment was manually examined and adjusted. Then, the evolutionary history was inferred using the Maximum Likelihood method implemented in RaxML (v8.2.4) (Stamatakis, 2015). The detailed parameters were “raxmlHPC-PTHREADS-SSE3 -f a -N 1000 -m PROTGAMMACPREV -x 551314260 -p 551314260 -o A_thaliana, N_tabacum -T 20”. The tree with the highest log likelihood (-126826.029496) was shown. The significance level for the phylogenetic tree was assessed by bootstrap testing with 1000 replications. The bootstrap values exceeding 50% were shown next to the corresponding nodes.

The complete plastid genome of S. cusia is a circular molecule of 144,133 bp in length. The genome possesses an LSC region of 91,666 bp, separated from the 17,811 bp SSC region by two IRs, each at 17,328 bp (Figure 1). The genome was deposited in GenBank under accession number: MG874806. A total of 113 genes are contained within the S. cusia cp genome, including 79 protein-coding genes, 30 tRNA genes, and 4 RNA genes (Table 1). Five protein-coding and seven tRNA-coding genes are duplicated on the IR regions. The LSC region contains 64 protein-coding genes and 22 tRNA genes, whereas the SSC region contains 11 protein-coding genes and one tRNA gene. Fifteen and three genes in the S. cusia cp genome possess one and two introns, respectively. Compared with the genes of R. breedlovei, the two genes matk and ycf2 are shorter at 453 and 663 bp, respectively, at the 5′-end sequence. The sequences of ycf2 and matk are confirmed via direct PCR amplifications and Sanger sequencing using the primers listed in Supplementary Table S1.

The overall GC content of the S. cusia cp genome is 38%. The GC content of the IR regions (46%) is higher than that of the LSC and SSC regions (37% and 32%). Overall, 49.4% of the S. cusia cp genome sequence is composed of genes that encode proteins. The overall GC content of the S. cusia cp genome is 38.38%, whereas that of the protein-coding regions is 38.00%. Within the protein-coding regions, the GC contents for the first, second, and third positions of codons are 46.17%, 38.51%, and 30.47%, respectively. The codon usage and codon-anticodon recognition patterns of the S. cusia cp genome are summarized in Supplementary Table S2. The 30 tRNA genes contain codons corresponding to all 20 amino acids that are necessary for protein biosynthesis.

Simple Sequence Repeats are valuable molecular markers of high-degree variations within the same species; these markers have been used in population genetics and polymorphism investigations (Xue et al., 2012). We analyzed the occurrence, type, and distribution of SSRs in the S. cusia cp genome. In total, 160 SSRs were identified in S. cusia cp genome by using the software tool MISA (Table 2). Among these SSRs, the majority consisted of mono- and dinucleotide repeats, which were found 94 and 43 times, respectively. By contrast, tri- (10), tetra- (11), and pentanucleotide repeat sequences (2) were observed with lower frequency. Most mononucleotide repeat sequences consisted of A/T repeats (92.6%). Similarly, 55.8% of the dinucleotide repeat sequences consisted of AT/AT repeats (Table 2). These results are in agreement with the previous findings stating that the cp SSRs are generally composed of short poly-A or poly-T repeats and rarely contained tandem G or C repeats (Kuang et al., 2011).

The locations of SSRs in the cp genome are shown in Supplementary Table S3. The corresponding locations were further classified as intergenic (IGS), exonic, and intronic regions. Several SSRs that were separated with distances less than 100 bp were considered together as a compound SSR. In total, 21 compound SSRs were identified. In terms of the distribution, 64, 46, and 18 repeats were located in the IGS, exon, and intron regions, respectively. In particular, two of them spanned over the IGS and CDS regions (repeats no. 51 and 101, Supplementary Table S3).

By definition, SSRs possess a repeat unit size equal to or less than five. We analyzed repeats with unit sizes exceeding five by using software tool Tandem Repeats Finder and REPuter. The repeats were divided into two types: tandem repeats and dispersed repeats. Tandem repeats are defined as those with two or more adjacent, approximate copies of a pattern of nucleotides. In total, 20 tandem repeats with repeat unit size longer than 12 bp (included) were identified, with the overall percentage of similarities between the adjacent repeat units set to 60% (Supplementary Table S4). The repeat unit was enclosed with parentheses. Among them, 4, 14, and 2 repeats were located in the CDS, intergenic, and intron regions, respectively. In contrast to tandem repeats, dispersed repeats is defined as with two or more approximate copies of a pattern of nucleotides separated by a certain number of nucleotides. Dispersed repeats can be further divided into four types, namely, forward, reverse, complement, and palindromic repeats (Kurtz et al., 2001). In total, two forward repeats and six palindromic repeats were identified using REPuter with a total length cutoff of 30 bp (Supplementary Table S5). For the palindromic repeat, the complementary structures were shown using pairing parentheses. This large set of repeats provides a rich resource for the development of novel molecular markers and to understand the genome evolution in S. cusia.

Analysis of the genome revealed a total of 18 genes with introns (Table 3). Among the genes, the rps12 is a trans-spliced gene; its 5′ end is located on the LSC region, and the 3′ end is located on the IR region. Two other genes, clpP and ycf3, presented two introns and three exons each. Nine genes, namely, atpF, ndhA, ndhB, petB, petD, rpl16, rpl2, rpoC1, and rps16, each with one intron, were identified. Finally, six tRNA genes with one intron and two exons each, namely, trnA-UGC, trnC-ACA, trnI-GAU, trnK-UUU, trnL-UAA, and trnS-CGA, were observed. The lengths of these introns range from 476 bp to 2391 bp, with the longest intron found in the trnK-UUU gene.

The loss of genes in cp genomes may result from the transfer of genes into the nuclear genome or the deletion of genes. To understand which genes may become indispensable during evolution, we performed hierarchical clustering of lost genes by using Ward’s method among 29 species in the Lamiales order (Figure 2). As shown in Figure 2, a total of 15 genes that have lost one gene in at least one species are clustered into five groups. The first group includes genes accD and ycf15, which are lost in three and two species, respectively. Genes psbB and ycf1, which are lost in R. breedlovei, are included in the second group. Genes psbI and rps19, which are lost in D. hygrometricum, are included in the third group. Genes rps12 and infA, which are included in the fourth group, are lost in five species, including Hesperelaea palmeri, three O. europaea subsp. maroccana, and Olea woodiana. The last group includes seven genes, namely, clpP, ndhD, petB, rpl16, rpoA, rps16, and rps4, which are all lost in T. tetragonolobum. In S. cusia, only the ycf15 gene is lost, indicating that the gene loss rate in S. cusia is lower than those in other species. The ycf15 gene belongs to the PFAM protein family PF10705, in which its function remains unknown. In certain plant species, the ycf15 gene is probably not a protein-coding gene because the proteins in these species possess premature stop codons (Steane, 2005). The loss of ycf15 gene results in further undermining of the notion that ycf15 plays important functions in the genome of S. cusia.

To identify the possible occurrence of genome rearrangement, we selected the cp genome sequences of S. cusia and eight other species belonging to Lamiales for synteny analyses. These eight species included A. paniculata, R. breedlovei, T. tetragonolobum, D. hygrometricum, S. miltiorrhiza, O. europaea, S. indicum, and S. takesimensis, which are members of the Acanthaceae, Bignoniaceae, Gesneriaceae, Lamiaceae, Oleaceae, Pedaliaceae, and Scrophulariaceae families, respectively (Figure 3). Members of the Orobanchaceae and Lentibulariaceae families were not included. Orobanchaceae are parasites; without chlorophyll, their cp genomes are considerably reduced both in size and gene content and suffer horizontal gene transfer from its host, gene loss, and gene pseudogenization (Li et al., 2013). An exclusive combination of losses and pseudogenization of the plastid NAD(P)H-dehydrogenase (ndh) gene complex was observed in the cp genome of the Lentibulariaceae family (Silva et al., 2016). Therefore, we did not discuss the species belonging to family Lentibulariaceae and Orobanchaceae in the phylogenetic analysis and IR contraction/extension analysis part. Pairwise cp genome comparison between S. cusia and the other eight cp genomes recovered a high degree of synteny (Figure 3). However, the S. cusia cp genome exhibits two unique features: first, the genes of psbA and trnH-GUG are duplicated in the S. cusia cp genome, which are located near the boundaries of the IRa-LSC, and IRb-LSC regions; second, the IR sequence of S. cusia is 17,328 bp in length, which is significantly shorter than those of A. paniculata (25,340 bp) and R. breedlovei (22,704 bp).

By far, only two cp genomes have been completely sequenced for species from Acanthaceae, a family under the order Lamiales. To determine the phylogenetic position of S. cusia in Lamiales, we obtained 29 complete cp genome sequences belonging to the Lamiales from the RefSeq database (Supplementary Table S6). The distributions of the 29 species among different families are as follows: Lamiaceae (16), Bignoniaceae (1), Acanthaceae (3), Gesneriaceae (1), Oleaceae (6), Pedaliaceae (1), and Scrophulariaceae (1), with the number shown in the enclosed parentheses representing the number of species in the corresponding clade. We extracted 65 protein sequences present among the 29 cp genomes and those of Nicotiana tabacum and Arabidopsis thaliana, which served as the outgroup taxa. Multiple sequence alignment of the amino acid sequences resulted in a total of 19,277 positions in the final dataset. The phylogenetic tree was constructed using the Maximum Likelihood method and is shown in Figure 4. The closest sister species of S. cusia are R. breedlovei and A. paniculata in Acanthaceae, (bootstrap values: 100%), and this finding was consistent with the expected relationships.

One distinct feature of the S. cusia cp genome is its IR length of 17,328 bp, which is significantly shorter than those of most angiosperms [ca. 20–28 kb, (Chumley et al., 2006)]. The shortening of the IR region suggests the occurrence of an IR contraction event. By contrast, the IRs of S. cusia spanned to include the trnH and psbA genes, indicating the occurrence of another IR expansion event. To identify any potential evolutionary events, we analyzed the gene contents at the border areas in detail. The regions around the IR-LSC and IR-SSC borders extending from all the available cp genomes in Lamiales were analyzed, except those from the family Lentibulariaceae (the bladderwort family) and Orobanchaceae (a family of parasitic plants) (Figure 5A). S. cusia is the only species in the Lamiales with an SSC region containing two copies of the trnH genes and two partial psbA genes. To determine the uniqueness of this pattern, we searched the SSC regions of all cp genomes and identified three species, Coriandrum sativum, Actinidia deliciosa, and Actinidia chinensis, whose SSC regions contain two copies of the trnH and psbA genes (Figure 5B).

The genes at the IR/SSC regions are rather conserved (Figure 5), given that ndhF (not shown in the alignment) and ycf1 genes are found in the SSC regions near the borders in almost all examined cp genomes examined; both genes appeared on the left side of the alignment in Figures 5A,B. By contrast, various patterns of expansion and contraction of the IR regions into or from the LSC regions were observed, as shown on the right side of the alignment in Figures 5A,B.

LSC/IR regions can be divided into four different types based on the patterns for the 29 cp genomes (Figure 5A) belonging to the Lamiales. Type I, the most dominant pattern, is found in 22 species, including D. hygrometricum, T. tetragonolobum, S. takesimensis, A. paniculata, Scutellaria insignis, Stenogyne kanehoana, Stenogyne haliakalae, Stenogyne bifida, Phyllostegia velutina, Haplostachys haplostachya, Stachys chamissonis, Stachys byzantina, Stachys coccinea, Stachys sylvatica, Premna microphylla, Tectona grandis, Jasminum nudiflorum, H. palmeri, O. woodiana subsp. woodiana, O. europaea, O. europaea subsp. maroccana, and O. europaea subsp. cuspidate. For this type, the genes included rpl2, rpl23, trnI, ycf2, trnL, ndhB, rps7, 3-rps12, trnV, rrn16, trnI, trnA, rrn23, rrn4.5, rrn5, trnR, trnN, and ycf1 (pseudogene), which are duplicated in the IR regions, except for rps19. The LSC/IR border patterns of O. woodiana subsp. woodiana, O. europaea subsp. maroccana, and O. europaea subsp. cuspidata are the same as that of O. europaea. As a result, their structures of the LSC/IR regions are not shown in Figure 5A.

The Type II pattern is found in five species, including Scutellaria baicalensis, S. miltiorrhiza, Rosmarinus officinalis, Lavandula angustifolia, and S. indicum. This pattern shows an rps19 pseudogene at the IR/LSC border. Meanwhile, Type III pattern is only found in R. breedlovei, in which the rps19, rpl2, and rpl22 are not found in the IR region, and an ycf2 pseudogene is found at the IR/LSC border.

Type IV is found only in the S. cusia cp genome. On the one hand, the genes of rps19, rpl2, rpl23, and ycf2 are found in the IRs. On the other hand, the IR extends into the LSC regions to include a trnH gene and a truncated psbA pseudogene at the IR/LSC border. A similar pattern has been observed in several species other than Lamiales (Figure 5B). For example, in the cp genomes of C. sativum and A. chinensis, the IR region extends into the psbA gene and engulfs a fragment of the psbA gene at the IR/LSC border. However, these patterns slightly vary. In the cp genome of C. sativum, a serial of genes, including rps19, rpl2, rpl23, ycf2, ndhB, rps7, and 3-rps12, are excluded from the IR regions. Meanwhile, in A. chinensis, the IR-LSC border expanded to the ycf2 gene. Compared with the cp genome of other Lamiales, the loss of genes rpl2, rpl23, and ycf2 in the IR regions may be responsible for the considerably shortened length of S. cusia cp genome. The various patterns of the IR/LSC border reflect the highly dynamic nature of the IR/LSC border regions.