Complete Chloroplast Genome Sequence of Aquilaria sinensis (Lour.) Gilg and Evolution Analysis within the Malvales Order.

Aquilaria sinensis (Lour.) Gilg is an important medicinal woody plant producing agarwood, which is widely used in traditional Chinese medicine. High-throughput sequencing of chloroplast (cp) genomes enhanced the understanding about evolutionary relationships within plant families. In this study, we determined the complete cp genome sequences for A. sinensis. The size of the A. sinensis cp genome was 159,565 bp. This genome included a large single-copy region of 87,482 bp, a small single-copy region of 19,857 bp, and a pair of inverted repeats (IRa and IRb) of 26,113 bp each. The GC content of the genome was 37.11%. The A. sinensis cp genome encoded 113 functional genes, including 82 protein-coding genes, 27 tRNA genes, and 4 rRNA genes. Seven genes were duplicated in the protein-coding genes, whereas 11 genes were duplicated in the RNA genes. A total of 45 polymorphic simple-sequence repeat loci and 60 pairs of large repeats were identified. Most simple-sequence repeats were located in the noncoding sections of the large single-copy/small single-copy region and exhibited high A/T content. Moreover, 33 pairs of large repeat sequences were located in the protein-coding genes, whereas 27 pairs were located in the intergenic regions. Aquilaria sinensis cp genome bias ended with A/T on the basis of codon usage. The distribution of codon usage in A. sinensis cp genome was most similar to that in the Gonystylus bancanus cp genome. Comparative results of 82 protein-coding genes from 29 species of cp genomes demonstrated that A. sinensis was a sister species to G. bancanus within the Malvales order. Aquilaria sinensis cp genome presented the highest sequence similarity of >90% with the G. bancanus cp genome by using CGView Comparison Tool. This finding strongly supports the placement of A. sinensis as a sister to G. bancanus within the Malvales order. The complete A. sinensis cp genome information will be highly beneficial for further studies on this traditional medicinal plant. Moreover, the results will enhance our understanding about the evolution of cp genomes of the Malvales order, particularly with regard to the role of A. sinensis in plant systematics and evolution.

Introduction

Plant chloroplasts (cps) are key organelles for photosynthesis and carbon fixation (Neuhaus and Emes, 2000). The cp genome sequence contains useful information in plant systematics because of its maternal inheritance in most angiosperms (Corriveau and Coleman, 1988; Zhang et al., 2003). Substitution rates in plant cp genomes are much lower than those in nuclear genomes (Wolfe et al., 1987). Plant cp genomes are valuable sources of genetic markers for phylogenetic analyses because of their very low level of recombination (Provan et al., 2001; Ravi et al., 2008). The cp DNA sequence was initially discovered during physical mapping of the Zea mays cp, which was produced by digestion with multiple restriction enzymes (Bedbrook and Bogorad, 1976). Subsequently, the first complete nucleotide sequence of Nicotiana tabacum was determined by the clone sequencing of plasmid and cosmid libraries (Kumano, 1986). Over 600 plant cp genomes have been sequenced and deposited in the NCBI Organelle Genome Resources. The cp genome evolution in land plants may be elucidated using these database resources. The cp in angiosperms exhibits a conserved quadripartite structure ranging from 115 to 165 kb in length and consists of one large single-copy (LSC) region, one small single-copy (SSC) region, and two copies of inverted repeat (IR; Palmer, 1991; Raubeson and Jansen, 2005). The Arabidopsis thaliana cp genome contains a circular DNA composed of 154,478 bp with 87 potential protein-coding genes (Sato et al., 1999). The cp DNA from maize (Z. mays) consists of 140,387 bp with a total of 104 genes (Bedbrook and Bogorad, 1976). The complete cp DNA of Cedrus deodara is circular molecule of 119,298 bp with 114 genes (Ching et al., 2010). However, some parasitic plants, such as Conopholis americana, which demonstrate unique life cycles, are exceptions because the sizes of their cp genomes are beyond 115–165 kb, with the smallest plastome of 45 kb in land plants (Wicke et al., 2013). The development of DNA sequencing technology has resulted in the extensive use of cp genomes for molecular marker and molecular phylogenetic studies (Tangphatsornruang et al., 2009; Takano and Okada, 2011; Awasthi et al., 2012; Jheng et al., 2012; Chen and Melis, 2013; Turner et al., 2013; Gaudeul et al., 2014).

Agarwood is widely used as a sedative, digestive, and antiemetic traditional drug. Agarwood sculpturing is valuable for interior decoration and is also popularly used as incense and perfume in Asia. The stems, branches, or roots of Aquilaria and Gyrinops trees are wounded and infected by fungi to produce agarwood (the wounds can be caused by wind, lightning strikes, gnawing of ants or insects, or microorganism invasion). Aquilaria sinensis is the only certified source for producing agarwood listed in China Pharmacopoeia (China pharmacopoeia Committee, 2010). All A. sinensis species are endangered because of the high demand for agarwood products; hence, these species are regulated under the Convention on International Trade in Endangered Species of Wild Fauna and Flora. However, the genomic resources for A. sinensis are limited, and little is known about the composition and organization of its cp genomes and their evolution. In this study, we report the complete cp genome sequence of A. sinensis (GenBank accession number: KT148967) in accordance with the Illumina Hiseq2500 standard protocol. Overall, the results provide basic genetic information on A. sinensis cp and the role of A. sinensis in plant systematics and evolution.

Materials and methods

DNA extraction and sequencing

Aquilaria sinensis fresh leaves were collected from a 2-year-old tree at the Experimental Farm of the Chinese Academy of Tropical Agriculture Sciences, Hainan, PR China. The leaves were cleansed, frozen in liquid nitrogen, and ground using a tissue lyser. DNA was extracted using a Plant Genomic DNA Kit (Foregene Biotech, China). DNA was used to generate 500 bp (insert size) paired-end library in accordance with the Illumina Hiseq2500 standard protocol. Approximately 3.1 Gb of raw data were generated with pair-end 125 bp read length.

De novo CP genome assembly

The obtained nucleotide sequencing reads were qualitatively assessed and assembled to contigs by using SOAPdenovo2 (Luo et al., 2012) with kmer length of 83. The assembled contigs included a mixture of sequences from organellar and nuclear genomes. The average coverage of cp genomes is usually much higher than that of nuclear genomes because many cps are found in a single cell (Steele et al., 2012; Straub et al., 2012). Thus, a complete de novo assembly of the cp genomes was performed using the assembly quality-filtered reads that exhibit high coverage for the cp genomes. We sorted the assembled contigs by contig-read depth analysis of assemblies by using the high correlation between sequencing depth and number of copies in the genome. The quality-filtered reads were remapped to the assembled contigs to calculate the sequencing depth with BWA (Li and Durbin, 2009). Thus, the cp contigs with high coverage (more than 500 ×) were isolated from the nuclear contigs by using the difference of read depths between contigs (Figure S1). All published cp genome sequences of dicotyledons were used as references to map the contigs with BLAST (Table S1) and thus confirm the cp genome contigs. Finally, all isolated cp contigs were combined, and reads were recaptured to isolate more cp DNA reads. Contigs were reassembled and extended to obtain a complete cp genome sequence.

Genomic annotation and analysis

Preliminarily gene annotation was performed using the online program Dual Organellar Genome Annotator (OGDRAW v1.2; Wyman et al., 2004) and cp Genome Annotation, Visualization, Analysis, and GenBank Submission Tool (Cheng et al., 2013) with plastid/bacterial genetic code and default conditions. Putative gene and protein sequences were BLAST-searched in non-redundant nucleotide database and non-redundant protein database to verify the exact gene and exon boundaries. All tRNA genes were further confirmed through online Trnascan-SE and tRNADB-CE search server (Griffiths-Jones et al., 2003; Schattner et al., 2005; Abe et al., 2011). The graphical map of the circular plastome was drawn using Organellar Genome DRAW (Lohse et al., 2007).

Identification of simple-sequence repeats (SSRs)

The genomic sequence was applied to exploit potential SSRs by using MISA software (http://pgrc.ipk-gatersleben.de/misa/). Tandem repeats of 1–6 nucleotides were considered as microsatellites. The minimum numbers of repeats were set to 10, 6, 5, 5, 5, and 5 for mono-, di-, tri-, tetra-, penta-, and hexa-nucleotides, respectively.

Long repeat analysis

Web-based REPuter (http://bibiserv.techfak.uni-bielefeld.de/reputer/) was used to analyze the repeat sequences, which included forward, reverse, and tandem repeats with minimal lengths of 30 bp and edit distances of less than 3 bp.

Codon usage

Codon usage was determined for all protein-coding genes. Statistical analyses of the distributions and visualization of codon usage in the form of heatmaps of 28 species of Angiosperms and histogram were conducted using R language with relative synonymous codon usage (RSCU) value (Sharp and Li, 1987).

RSCU is a simple measure of non-uniform usage of synonymous codons in a coding sequence. The RSCU value is the number of times a particular codon is observed, relative to the number of times that the codon would be observed for a uniform synonymous codon usage (i.e., all codons for a given amino acid exhibit similar probabilities). The RSCU value in the absence of any codon usage bias is 1.00, which is the case for the CDS sequence in the following example. A codon used less frequently than expected will achieve RSCU of <1.00, whereas codons used more frequently than expected may reach RSCU of >1.00.

Phylogenetic analysis

The jModeltest 0.1.1 software was employed to analyze the general GTR+G+I model for nucleotide sequence and HIVb+I+G model for protein sequence by using optimized parameters (Posada, 2008). Phylogenetic analysis was subsequently performed using Maximum likelihood (ML) and Bayesian inference (BI) methods. ML analysis was conducted using RAxML8.1.5 with 1000 bootstrap replicates (Stamatakis, 2014). BI analysis was conducted using Phylobayes 4.1b with two chain max diff <0.01 (Lartillot et al., 2009).

CGView comparison tool (CCT) map

The A. sinensis cp genome was compared with other available cp genomes of Malvales by using CCT (Grant and Stothard, 2008). Genes were signed by Clusters of Orthologous Groups, and BLAST was used to align other genomes to A. sinensis. The results are shown as a circular map. AT distributions were measured on the basis of AT skewed using the equation: AT-skew = (A−T)/(A+T).

Results and discussion

Genome sequencing and assembly

A total of 2.48 × 107 reads with an average read length of 125 bp were obtained after low-quality bases and adapter sequences were trimmed. De novo assembly produced 691,722 contigs (2.78%). The size of the A. sinensis cp genome was 159,565 bp (Figure 1). The genome included an LSC region of 87,482 bp, an SSC region of 19,857 bp, and a pair of IRs (IRa and IRb) of 26,113 bp each (Table 1). The GC content was 37.11% (Table S2). However, the GC content was unevenly distributed in the entire cp genome, with the highest value in the IR regions (42.86%), followed by the LSC (34.95%) and SSC (31.58%) regions. The frequency of codon usage was deduced for the cp genome on the basis of the CDS sequences. Notably, the AT contents were 54.64, 62.31, and 69.34% at the first, second, and third codon positions, respectively, within the protein-coding regions (Table S2). Bias toward higher AT content at the third codon position was consistent with the enrichment of A and T, which has been widely observed in many other sequenced land plant cp genomes (Morton, 1998; Tangphatsornruang et al., 2009; Nie et al., 2012; Qian et al., 2013). The sequences of the A. sinensis cp genome were deposited in GenBank with accession number KT148967.

Figure 1

Gene map of . Organization of the cp genome of A. sinensis showing annotated genes. Genes drawn inside the circle are transcribed clockwise, and those outside are counter clockwise. Genes belonging to different functional groups are color-coded. The inner circle shows the locations of the large single-copy region, small single-copy, and the pair of inverted repeats (IRa and IRb). The darker gray in the inner circle corresponds to GC content, whereas the lighter gray corresponds to AT content.

Table 1

Gene conten.

Category of genes	Group of gene	Name of gene
Self-replication	Small subunit of ribosome	rps7* rps11 rps15 rps14 rps19 rps18 rps4 rps3 rps2 rps8 rps12* rps16
	Large subunit of ribosome	rpl23* rpl20 rpl22 rpl2* rpl36 rpl32 rpl33 rpl14 rpl16
	DNA-dependent RNA polymerase	rpoB rpoA rpoC2 rpoC1
	Ribosomal RNA genes	rrn16S* rrn5S* rrn23S* rrn4.5S*
	Transfer RNA genes	trnQ-TTG trnR-ACG* trnM-CAT trnY-GTA trnH-GTG trnA-TGC* trnP-TGG
		trnS-GCT trnN-GTT* trnL-TAG trnG-GCC trnW-CCA trnK-TTT
		trnT-TGT trnL-CAA*trnS-TGA trnI-CAT* trnfM-CAT trnF-GAA trnS-GGA
		trnV-GAC* trnC-GCA trnD-GTC trnE-TTC trnI-GAT* trnT-GGT
Genes for photosynthesis	Subunits of NADH dehydrogenase	ndhB* ndhC ndhA ndhF ndhG ndhD ndhE ndhJ ndhK ndhH ndhI
	Large subunit of Rubisco	rbcL
	Subunits of photosystem II	psbE psbD psbF psbA psbC psbB psbM psbL psbN psbI psbH psbK psbJ psbT psbZ
	Subunits of photosystem I	psaI psaJ psaB psaC ycf4 ycf3 psaA
	Subunits of ATP synthase	atpI atpH atpB atpA atpF atpE
	Subunits of cytochrome	petD petG petA petB petL petN
Other genes	Envelope membrane protein	cemA
	C-type cytochrome synthesis gene	ccsA
	Subunit of acetyl-CoA	accD
	Protease	clpP
	Maturase	matK
	Component of TIC complex	ycf1
Genes of unknown function	Conserved open reading frames	ycf2* orf42 orf56
Pseudogenes		ycf1* ycf15*

Duplicated gene.

Genomic annotation

The draft genome was drawn using OGDRAW v1.2 (Figure 1). The single collapsed IR contig was separated into two repeat regions. Assembly of the two IRs and LSC and SSC contigs covered the complete sequence without gaps. The positions of all genes identified in the cp genome and functional categorization of these genes are presented in Figure 1. The A. sinensis cp genome was 159,565 bp long with a typical quadripartite structure. A total of 113 functional genes were identified, which comprised 82 protein-coding genes, 27 tRNA genes, and 4 rRNA genes (Table 1). Comparing to the genes in other species (Figure S2), little change was found in gene structure. The very low level of recombination was also reported in the cp genome of land plant (Provan et al., 2001; Ravi et al., 2008). Among the 82 protein-coding genes, 75 were single-copy genes, and 7 were duplicates. Among the 31 RNA genes, 20 were unique, and 11 were duplicates. Among the 113 unique genes, 9 genes contained 1 intron (7 protein-coding and 2 tRNA genes), and 1 gene (ycf3) contained 2 introns (Table S3). The ycf3 gene was similar to those in Globe artichoke and Metasequoia glyptostroboides (Chen et al., 2015; Curci et al., 2015). Out of the 10 genes with introns, 3 protein-coding genes were located in the LSC, 1 in the SSC, and 6 (4 protein-coding genes and 2 tRNAs) in the IR region. The ndhA gene presented the largest intron (1148 bp). In addition, ndhB and rpl2 were identified as duplicate genes.

SSR analysis

SSRs consist of 1–6 nucleotide repeat units, which are also known as microsatellites and short tandem repeats (Chen et al., 2006). SSRs are important in plant typing (Yang et al., 2011; Xue et al., 2012) and widely used for genetic molecular markers in population genetics (Doorduin et al., 2011; He et al., 2012). A total of 45 SSR regions were identified using the microsatellite identification tool (MISA) in A. sinensis cp genome (Table 2), accounting for 499 bp of the total sequence (0.3%), and 37 SSRs were only composed of A or T bases. Two SSRs were composed of C bases, and six SSRs were composed of dinucleotide (AT/TA/TC) repeats. Therefore, SSRs in A. sinensis cp genome were rich in AT. Poly(A)/(T) had been reported to exhibit higher proportion relative to poly(G)/(C) in many plant families (Kumar et al., 2009; Melotto-Passarin et al., 2011; Nie et al., 2012; Martin et al., 2013). Among these SSRs, 36 SSRs were located in noncoding sections of the LSC/SSC region, and 9 SSRs in protein-coding genes (rpoC2, rpoB, psbF, cemA, psbN, rps19, and ycf1). No tri- or tetra-nucleotide repeats over 15 bp long were found. The SSRs identified in this study may provide a new perspective to refine the phylogeny and elucidate the origin of cultivars.

Table 2

Simple sequence repeats in .

cpSSR ID	Repeat motif	Length (bp)	Start	End	Region	Annotation
1	(A)10	10	1883	1892	LSC
2	(T)10	10	2023	2032	LSC
3	(A)12	12	4127	4138	LSC
4	(A)12	12	4803	4814	LSC
5	(C)10	10	5260	5269	LSC
6	(A)10	10	6662	6671	LSC
7	(A)15	15	7644	7658	LSC
8	(A)11	11	8033	8043	LSC
9	(T)11	11	8380	8390	LSC
10	(T)11	11	9080	9090	LSC
11	(T)10	10	9863	9872	LSC
12	(TA)6	12	10743	10754	LSC
13	(A)10	10	13927	13936	LSC
14	(T)10	10	14100	14109	LSC
15	(TC)7	14	17620	17633	LSC
16	(T)11	11	18009	18019	LSC	rpoC2
17	(T)10	10	19917	19926	LSC	rpoC2
18	(T)10	10	27649	27658	LSC	rpoB
19	(T)10	10	30717	30726	LSC
20	(A)10	10	32194	32203	LSC
21	(T)10	10	44610	44619	LSC
22	(A)10	10	45233	45242	LSC
23	(A)10	10	47349	47358	LSC
24	(A)11	11	47627	47637	LSC
25	(T)10	10	48204	48213	LSC
26	(T)12	12	51015	51026	LSC
27	(T)11	11	51316	51326	LSC
28	(T)10	10	51846	51855	LSC
29	(A)11	11	61867	61877	LSC
30	(T)11	11	62016	62026	LSC
31	(A)11	11	62578	62588	LSC
32	(TA)6	12	62964	62975	LSC
33	(T)10	10	64002	64011	LSC	psbF
34	(C)10	10	66343	66352	LSC
35	(T)12	12	67133	67144	LSC	cemA
36	(T)11	11	71791	71801	LSC
37	(TA)6	12	77233	77244	LSC	psbN
38	(T)15	15	87280	87294	LSC	rps19
39	(TA)7	14	113802	113815	SSC
40	(A)10	10	116349	116358	SSC
41	(A)10	10	116494	116503	SSC
42	(T)10	10	117050	117059	SSC
43	(AT)7	14	117644	117657	SSC
44	(T)10	10	129560	129569	SSC	ycf1*
45	(T)13	13	132293	132305	SSC	ycf1

Duplicated gene.

Large repeat analysis

Large repeat sequences showed repeats with length of ≥30 bp each. Sixty pairs of large repeat sequences with sequence identity of >90% were found in the A. sinensis cp genomes (Table 3). The repeats ranged from 30 to 600 bp in length and were repeated twice. A total of 33 large repeat sequences were located in protein-coding genes (e.g., ycf1 and ycf2), and 27 large repeat sequences were located in the intergenic regions. Numerous repeated sequences were identified in cp genomes, particularly in the intergenic spacer regions, and have been reported in several angiosperm lineages (Yang et al., 2013).

Table 3

Long repeat sequences in .

id	Repeat Start 1	Type	Size(bp)	Repeat Start 2	Mismatch(bp)	E-value	Gene	Region
1	1471	F	41	71846	−3	4.26E-10	IGS	LSC
2	1480	F	41	71855	−3	4.26E-10	IGS	LSC
3	1574	F	163	71953	−3	1.00E-81	IGS	LSC
4	1587	F	156	71966	−3	1.44E-77	IGS	LSC
5	8609	F	31	37664	−3	1.88E-04	IGS	LSC
6	41152	F	41	43376	−2	1.09E-11	psaB(CDS); psaA(CDS)	LSC
7	46415	F	36	101941	−1	1.64E-10	ycf3(intron); IGS	LSC; IRA
8	46415	F	36	124843	−3	2.92E-07	ycf3(intron); ndhA(inton)	LSC; SSC
9	46424	F	30	101950	−1	5.59E-07	ycf3(intron); IGS	LSC; IRA
10	70177	F	30	70366	−3	6.81E-04	accD(CDS)	LSC
11	70195	F	67	70273	−2	6.54E-27	accD(CDS)	LSC
12	70195	F	36	70291	−2	8.60E-09	accD(CDS)	LSC
13	70199	F	53	70259	−1	1.40E-20	accD(CDS)	LSC
14	70218	F	34	70260	−1	2.47E-09	accD(CDS)	LSC
15	70260	F	31	70296	−1	1.44E-07	accD(CDS)	LSC
16	95892	F	31	95928	−3	1.88E-04	ycf2(CDS)	IRA
17	96540	F	31	150227	−3	1.88E-04	ycf15(CDS)	IRA; IRB
18	97906	F	30	148863	−2	2.43E-05	IGS	IRA; IRB
19	101938	F	41	124840	−2	1.09E-11	IGS; ndhA(intron)	IRA; SSC
20	112678	F	31	112720	0	1.55E-09	IGS	IRA
21	113385	F	55	113454	0	5.52E-24	IGS	IRA
22	118510	F	64	118590	−1	4.04E-27	ccsA(CDS)	SSC
23	119996	F	32	120064	−3	5.20E-05	IGS	SSC
24	131641	F	30	131707	−3	6.81E-04	ycf1(CDS)	SSC
25	132792	F	31	132831	−2	6.50E-06	ycf1(CDS)	SSC
26	133153	F	39	133216	−3	5.85E-09	ycf1(CDS)	SSC
27	133168	F	38	133231	−1	1.08E-11	ycf1(CDS)	SSC
28	133282	F	62	133351	0	3.37E-28	ycf1(CDS)	SSC; IRB
29	150839	F	31	150875	−3	1.88E-04	ycf2(CDS)	IRB
30	8610	I	30	47858	−1	5.59E-07	IGS	LSC
31	37665	I	30	47858	−3	6.81E-04	IGS	LSC
32	46415	I	36	144821	−1	1.64E-10	ycf3(intron); IGS	LSC; IRB
33	46424	I	30	144818	−1	5.59E-07	ycf1(CDS)	LSC; IRB
34	89862	I	30	156882	−3	6.81E-04	ycf2(CDS)	IRA; IRB
35	89886	I	30	156906	−3	6.81E-04	ycf2(CDS)	IRA; IRB
36	94593	I	38	152149	−2	6.00E-10	ycf2(CDS)	IRA; IRB
37	94611	I	38	152167	−2	6.00E-10	ycf2(CDS)	IRA; IRB
38	95892	I	31	150839	−3	1.88E-04	ycf2(CDS)	IRA; IRB
39	95900	I	41	150839	0	1.48E-15	ycf2(CDS)	IRA; IRB
40	95918	I	41	150857	0	1.48E-15	ycf2(CDS)	IRA; IRB
41	95928	I	31	150875	−3	1.88E-04	ycf2(CDS)	IRA; IRB
42	112668	I	41	134047	−1	1.82E-13	IGS	IRA; IRB
43	112668	I	62	134047	−2	5.73E-24	IGS; ycf1(CDS),IGS	IRA; IRB
44	112678	I	31	134047	0	1.55E-09	IGS	IRA; IRB
45	112689	I	62	134068	−2	5.73E-24	IGS	IRA; IRB
46	112699	I	31	134047	−1	1.44E-07	IGS	IRA; IRB
47	112710	I	41	134089	−1	1.82E-13	IGS	IRA; IRB
48	112720	I	31	134089	0	1.55E-09	IGS	IRA; IRB
49	112720	I	31	134068	−1	1.44E-07	IGS	IRA; IRB
50	113385	I	70	133274	−1	1.08E-30	IGS; ycf1(CDS)	IRA
51	113454	I	55	133358	0	5.52E-24	IGS; ycf1(CDS)	IRA; SSC
52	113526	I	530	117403	0	0.00E+00	IGS,rpl(CDS)	IRA; SSC
53	124840	I	41	144819	−2	1.09E−11	ndhA(CDS,intron); IGS	SSC; IRB
54	59836	T	31	59866	0	1.55E-09	rps18	LSC
55	112668	T	41	112710	−1	1.82E−13	IGS	IRA
56	118496	T	78	118576	−3	1.61E−31	IGS	SSC
57	132782	T	41	132821	−3	4.26E−10	ycf1(CDS)	SSC
58	133193	T	44	133235	0	2.31E−17	ycf1(CDS)	SSC
59	133274	T	70	133343	−1	1.08E-30	ycf1(CDS)	SSC; IRB
60	134047	T	41	134089	−1	1.82E-13	IGS	IRB

F, Forward; I, Inverted; T, Tandem; IGS, intergenic space.

Most protein-coding genes in these basal eudicots employ the standard ATG as the initiator codon. However, ATA, ATC, TTG, and ATT are also used as alternatives to ATG as the start codon. Among the A. sinensis cp protein-coding genes, nine genes were used alternatively to ATG as the start codon as follows: ATA for atpF; ATT for ycf1 and petB; ATC for rpl16; GTG for rps8, psbC, and ndhD; and TTG for ndhA and rpoC1. In the N. tabacum cp genome, GTG was used as start codon for rps19, psbC, and ycf15, whereas ACG was used for psbL and ndhD (Sugiura et al., 1998). ACG and GTG were used as start codon for rpl2 and rps19, as reported in Oryza sativa (Liu and Xue, 2004).

Furthermore, the codon usage patterns of the 82 distinct cp protein-coding genes in A. sinensis were examined. All the protein-coding genes were composed of 26,160 codons. Interestingly, as synonymous codons, almost each of these codons contained half synonymous codon, which ended with A or T with high RSCU values, and the other half ended with C or G with low RSCU values (Table S4). These codon usage patterns may be driven by the composition bias of the high proportion of A/T similar to those of other reported cp genomes (Raubeson et al., 2007; Delannoy et al., 2011) and mitochondrial genomes (Barth and Berendonk, 2011).

Figure 2 shows that the RSCU value increased with the number of codons that code a particular amino acid. The high RSCU value was probably attributed to the function of the amino acid or the structure of the peptide to avoid error in transcription.

Figure 2

Codon content of 20 amino acid and stop codon of 82 coding genes of . Color of the histogram is corresponding to the color of codons.

Statistical analyses of the distributions and visualization of codon usage in the form of heatmaps of 28 species of Angiosperms (Figure 3) showed that approximately half of the codons were not frequently used. These codons were denoted in blue, which indicated RSCU value of <1 and weak codon bias. Almost two-thirds of all the codons with high RSCU values ended with purine (A/T). Thus, we hypothesized that the codon in A. sinensis cp genome bias ended with A/T. This phenomenon was also found in many other plant and algal lineages (Morton, 1998). The distribution of codon usage in A. sinensis cp genome was most similar to the codon usage in Gonystylus bancanus cp genome.

Figure 3

Codon distribution of all merged protein-coding genes for all considered species. Color key: higher red values indicate higher RSCU value and lower blue values indicate lower RSCU value; A Hierarchical clustering (Lance and Williams, 1967; average linkage method) was performed for codon patterns (x-axis).

Phylogenetic analyses

cp genomes are significant in the research of phylogenetics, evolution, and molecular systematics. Numerous analytical studies have been conducted in the past decade to address phylogenetic issues at deep nodes by comparing multiple protein-coding genes (Lemieux et al., 2000; Delas Rivas et al., 2002; Moore et al., 2010) and complete sequences in cp genomes (Goremykin et al., 2004; Moore et al., 2007). These studies enhanced the understanding about enigmatic evolutionary relationships among Angiosperms. We examined the phylogenetic position of A. sinensis and relationships within Angiosperms. We previously selected 82 protein-coding genes commonly present in cp genomes of 29 species, including the A. sinensis cp genome sequenced in the current study. ML and BI nucleic acid analyses were performed, and the results are summarized in Figures 4, 5. Similar phylogenetic topologies were found in the ML and BI nucleic acid analyses. Bootstrap values were very high, and 28 of 29 nodes with 100% bootstrap values were found using ML. Up to 25 out of the 29 nodes with bootstrap values of ≥99% were found using BI. A. sinensis and G. bancanus were grouped in Malvales, with 100% bootstrap values in the ML and BI phylogenetic trees. Similarly, in the ML protein analyses, 27 of the 29 nodes yielded bootstrap values of 100%, and 24 nodes reached ≥99% in the BI. ML and BI protein analyses showed that A. sinensis can also be grouped with G. bancanus within Malvales. The phylogenetic results strongly support the position of A.sinensis within the Malvales order. However, the results were inconsistent because A.sinensis was classified into Myrtales according to the traditional morphological classification of China.

Figure 4

Phylogenetic position of . The first number above the lines indicates the BI bootstrap value of the nucleic acid analysis for each clade, whereas the second number indicates the ML bootstrap value. The position of A.sinensisis shown in boldface.

Figure 5

Phylogenetic position of . The first number above the lines indicates the BI bootstrap value of the protein analysis for each clade, whereas the second number indicates the ML bootstrap value. The position of A.sinensis is shown in boldface.

CCT map

Five available cp genomes of Malvales were compared with the A.sinensis cp genome by using CCT (Figure 6). The sequence identity between the A.sinensis cp and other representatives of the Malvales was analyzed. The results showed that G. bancanus cp genome achieved the highest sequence similarity (>90%), which was consistent with the result of the phylogenetic analysis. The highest similar region across all five cp genomes occurred in the IR region. LSC and SSC regions were less conservative; thus, many regions with identity results were lower than 90%, such as the ycf1 gene in the SSC region and intergenic space regions. This evolutionary conserved feature was reported in the IR region in cp (Curtis and Clegg, 1984; Palmer, 1985).

Figure 6

Genome comparison of five CP genomes of Malvales to . From the outer to the inner color ring: Gonystylus bancanus, Theobroma cacao, Gossypium longicalyx, Hibiscus syriacus, and Gossypium bickii. BLAST was used to align other genomes to A. sinensis, and the results are shown with a circular map. The color codes are based on the similarity score, that is, dark red and blue depict similarity scores of 100%, above 90%, and below 90%, respectively. The four outer narrow rings are the protein-coding gene positions based on the A. sinensis cp genome. The color codes are based on Clusters of Orthologous Groups. The innermost ring is AT skew in the A. sinensis. AT skew+ indicate A>T, AT skew- indicate A < T.

Conclusions

The complete cp sequence of traditional medicinal plant A. sinensis was assembled, annotated, and analyzed. The cp genome consisted of one LSC, one SSC, and two IR regions. A total of 45 polymorphic SSR loci and 60 pairs of large repeats were identified in the A. sinensis cp genome. These repeat motifs can be selected to develop markers and conduct phylogenetic analysis. Both ML and BI phylogenetic analyses strongly supported the position of A. sinensis as a sister to G. bancanus within the Malvales order. The distribution of codon usage in A. sinensis cp genome was most similar to that in G. bancanus cp genome. CCT analytical results also indicated that the A. sinensis cp genome achieved higher sequence similarity than the G. bancanus cp genome. However, the traditional morphological classification of China classified A. sinensis into Myrtales. The data obtained in this study will be beneficial for further investigations on A. sinensis. Moreover, the results will help expand the understanding about the evolutionary history of the Malvales order, particularly regarding the role of A. sinensis in plant systematics and evolution.

Author contributions

YW designed and performed the experiment as well as drafted the manuscript. XC and SP conceived the study and revised the manuscript. DZ, XJ, WM, and HD prepared the samples and analyzed part of the data. All the authors have read and approved the final manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.