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 FIGURE 1 | Gene map of A. sinensis chloroplast (cp) genome sequence. 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. 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 singlecopy (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

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  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 FIGURE 3 | Codon distribution of all merged protein-coding genes for all considered species.  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.

DNA Extraction and Sequencing
Aquilaria sinensis fresh leaves were collected from a 2-yearold 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).

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 .
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).

Genome Sequencing and Assembly
A total of 2.48 × 10 7 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.

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 proteincoding 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.

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).

Codon Usage
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.
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 FIGURE 6 | Genome comparison of five CP genomes of Malvales to A.sinensis. 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. not frequently used. These codons were denoted in blue, which indicated RSCU value of <1 and weak codon bias. Almost twothirds 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.

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 proteincoding 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.

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).

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.   Table S1 | List of CP genomes used as reference to map to contigs.