Plastid genome sequence of a wild woody oil species, Prinsepia utilis, provides insights into evolutionary and mutational patterns of Rosaceae chloroplast genomes.

BACKGROUND: Prinsepiautilis Royle is a wild woody oil species of Rosaceae that yields edible oil which has been proved to possess particular benefits for human health and medical therapy. However, the lack of bred varieties has largely impeded exploiting immense potentials for high quality of its seed oil. It is urgently needed to enlarge the knowledge of genetic basis of the species and develop genetic markers to enhance modern breeding programs.
RESULTS: Here we reported the complete chloroplast (cp) genome of 156,328 bp. Comparative cp sequence analyses of P. utilis along with other four Rosaceae species resulted in similar genome structures, gene orders, and gene contents. Contraction/expansion of inverted repeat regions (IRs) explained part of the length variation in the Rosaceae cp genomes. Genome sequence alignments revealed that nucleotide diversity was associated with AT content, and large single copy regions (LSC) and small single copy regions (SSC) harbored higher sequence variations in both coding and non-coding regions than IRs. Simple sequence repeats (SSRs) were detected in the P. utilis and compared with those of the other four Rosaceae cp genomes. Almost all the SSR loci were composed of A or T, therefore it might contribute to the A-T richness of cp genomes and be associated with AT biased sequence variation. Among all the protein-coding genes, ycf1 showed the highest sequence divergence, indicating that it could accomplish the discrimination of species within Rosaceae as well as within angiosperms better than other genes.
CONCLUSIONS: With the addition of this new sequenced cp genome, high nucleotide substitution rate and abundant deletions/insertions were observed, suggesting a greater genomic dynamics than previously explored in Rosaceae. The availability of the complete cp genome of P. utilis will provide chloroplast markers and genetic information to better enhance the conservation and utilization of this woody oil plant.

Introduction

Rosaceae, containing about 3,000 species in over 100 genera, is the third economically important plant family in temperate regions [1]. A vast majority of common fruits, including Malus (apples), Pyrus (pears), Prunus (plums, peaches, cherries, almonds, and apricots), Rubus (raspberries), and (strawberries), as well as plenty sources of ornamentals, e.g., Rosa (roses), (cinquefoils), and Sorbus (mountain ashes) are produced by members of this family. From the taxonomical perspective, Rosaceae has traditionally been divided into four subfamilies: Rosoideae (fruit, achene), Prunoideae (fruit, drupe), Spiraeoideae (fruit, follicle or capsule) and Maloideae (fruit, pome) on the basis of fruit type [1,2]. While it was classified into three to twelve subfamilies historically based on molecular phylogenetic studies [1,3]. The economic importance of rosaceous crops has drawn several genome sequencing projects to these species. To date, five nucleic genomes, (strawberry) [4], Malus×domestica (apple) [5], (pear) [6], (peach), and [7], were completely sequenced. Compared with the availability of so many nucleic genome sequences, our understanding of Rosaceae’s cp genomes is left behind, as only five cp genomes in this family (strawberry [4], apple [5], pear [8], peach [9], and ) have been sequenced so far. Thus Rosaceae stands out as a worldwide distributed and economically important but cp genome-wide limited and taxonomically ambiguous family of angiosperms.

, one of the major woody oil plants in Rosaceae, is mainly distributed in mountainous regions with high elevations of 1000-3200 m in southwestern China, Pakistan, India, Nepal and Bhutan[10,11]. . In addition to , there are the other three related species, (Oliv.) Oliv. ex Bean, Batal, and Hayata, in the genus . Seed kernels of averagely yield 37.2% of semi-drying and pale yellow fatty oil which has widely been consumed for cooking [10]. Up to date, the chemical components of its seed oil (e.g., oleic acid 32.6%; linoleic acid 43.6%; palmitic acid 15.2%) have been proved to possess particular benefits or potentials for human health and medical therapy [12,13], which can be undoubtedly comparable with Mediterranean olive oil. It is well known that the majority of woody oil plants grow in tropical and subtropical regions. may represent as the best candidate oil plant which is suitable for extensive cultivation with cold weather in higher altitudes of temperate regions. The lack of bred varieties in , however, has largely impeded exploiting immense potential for high quality of seed oil. Nowadays, the collection and transplanting of wild plants have seriously destroyed natural populations, making continuous loss of numerous useful germ plasm resources. It is urgently needed to enlarge the knowledge of genetic basis of the species and develop genetic markers that will be fairly necessary to enhance modern breeding programs for the purpose of the selection and breeding of new varieties.

The affinity of has been argued by many botanists since the species was reported in 1839 [14]. Based on morphological characteristics, Royle (1839) and Lindley (1853) assigned this newly described species to the tribe Chrysobalanoideae [15,16], Bentham and Hooker (1862-1867) placed it in the tribe Prunoideae [17], while Sterling (1963) assigned it to a new subfamilial group in Prinsepioideae, Rosaceae [14]. Its subfamilial characters are carpel solitary, style lateral or sub-basal; 2 ovules erect and campylotropous, the nucelli horizontal, micropyles directed toward the ventral obturator; fruit with leathery pit [14]. In recent years, it has been placed into Osmaronieae (supertribe kerriodae and subfamily Spiraeoideae) with molecular phylogenetic support [1].

Chloroplasts were derived from endosymbiosis between independent living cyanobacteria and a non-photosynthetic host [18]. This intracellular organelle encodes a number of chloroplast-specific components and has high levels of conservation [19]. Chloroplast DNA (cpDNA) of green plants provides sufficient information for genome-wide evolutionary studies. It has shown great potentials in addressing phylogenetic questions at both high and low taxonomic levels, and sometimes it is necessary to use complete cp genome sequences for resolving complex evolutionary relationships [20-22]. However, acquiring large coverage of cp genomes has typically been limited by conventional DNA sequencing. As next-generation sequencing techniques have revolutionized DNA sequencing via high-throughput capabilities and relatively low costs [23], it is now more convenient to obtain cp genome sequences and extend gene-based phylogenetics to phylogenomics.

In this study, we sequenced the complete cp genome sequence of a wild woody oil plant, , which boasted great economic and taxonomic values in Rosaceae. Comparative analyses of the five Rosaceae plastomes, including the plastome sequenced in this study, were further performed to provide insights into overall and especially evolutionary dynamics of Rosaceae cp genomes by detecting genome-wide sequence variations.

Results and Discussion

Genome organization of

The cp genome of exhibited a circular DNA molecule of 156,328 bp with a typical quadripartite structure resembling to the majority of land plant cp genomes. It consisted of two inverted repeat regions (IRa and IRb) of 26,302 bp separated by large (LSC) and small (SSC) single-copy regions of 85,239 and 18,485 bp, respectively (Figure 1). The cp genome encoded 112 unique genes, 19 of which were duplicated in the IR regions, giving a total of 131 genes (Figure 1, Table 1). Among these unique genes, 18 included one or two introns. All of these coding regions accounted for 51.3% of the whole genome. The gene contents and gene order of were identical to those of the four previously known Rosaceae plastomes (Table 2).

Figure 1

Gene map of the chloroplast genome.

Genes lying outside of the outer circle are transcribed in the clockwise direction whereas genes inside are transcribed in the counterclockwise direction. Genes belonging to different functional groups are color coded. Area dashed darker gray in the inner circle indicates GC content while the lighter gray corresponds to AT content of the genome.

Table 1

List of genes encoded by chloroplast genome.

Category for genes	Group of gene	Name of gene
Self replication	Large subunit of ribosome	rpl2*, rpl14, rpl16*, rpl20, rpl22, rpl23, rpl32, rpl33, rpl36
	Small subunit of ribosome	rps2, rps3, rps4, rps7, rps8, rps11, rps12, rps14, rps15, rps16*, rps18, rps19
	DNA dependent RNA polymerase	rpoA, rpoB, rpoC1*, rpoC2
	rRNA genes	rrn4.5, rrn5, rrn16, rrn23
	tRNA genes	trnA-UGC*, trnC-GCA, trnD-GUC, trnE-UUC, trnF-GAA, trnG-UCC, trnG-GCC*, trnH-GUG, trnI-CAU, trnI-GAU*, trnK-UUU*, trnL-UAG, trnL-CAA, trnL-UAA*, trnM-CAU, trnfM-CAU, trnN-GUU, trnP-UGG, trnQ-UUG, trnR-ACG, trnR-UCU, trnS-GGA, trnS-GCU, trnS-UGA, trnT-GGU, trnT-UGU, trnV-UAC*, trnV-GAC, trnW-CCA, trnY-GUA
Photosynthesis	Subunits of photosystemⅠ	psaA, psaB, psaC, psaI, psaJ, ycf3**, ycf4
	Subunits of photosystemⅡ	psbA, psbB, psbC, psbD, psbE, psbF, psbH, psbI, psbJ, psbK, psbL, psbM, psbN, psbT, psbZ
	Subunits of NADH-dehydrogenase	ndhA*, ndhB*, ndhC, ndhD, ndhE, ndhF, ndhG, ndhH, ndhI, ndhJ, ndhK
	Subunits of cytochrome b/f complex	petA, petB*, petD*, petG, petL, petN
	Subunits of ATP synthase	atpA, atpB, atpE, atpF*, atpH, atpI,
	Large subunit of rubisco	rbcL
Other genes	Translational initiation factor	infA
	Maturase	matK
	Protease	clpP**
	Envelop membrane protein	cemA
	Subunit of Acetyl-CoA-carboxylase	accD
	c-type cytochrom synthesis gene	ccsA
Genes of unknown function	Conserved Open Reading Frames (ORF, ycf)	ycf1, ycf2

One and two asterisks indicate one- and two-intron containing genes respectively

Table 2

Summary of the Rosaceae chloroplast genome features.

	P . utilis	P. persica	P . pyfifolia	P . rupicola	F. vesca
Genome size (bp)	156,328	157,790	159,922	156,612	155,691
LSC length (bp)	85,239	85,968	87,901	84,970	85,606
SSC length (bp)	18,485	19,060	19,237	18,942	18,175
IR length (bp)	26,302	26,381	26,392	26,350	25,555
Number of genes	131	131	131	131	131
Percent of coding regions	51.3%	50.8%	50.1%	51.2%	51.5%

From the aspect of genome size, , next to (155,691 bp), is the second smallest one among the five sequenced Rosaceae cp genomes. It is about 3.6 kb, 1.5 kb, and 0.3 kb smaller than , , and , respectively. The genome size variation can be explained mainly by differences in the length of SSC and IR regions (Table 2).

Contraction/expansion of IRs among Rosaceae cp genomes

Generally, the lengths of IR (IRa and IRb) regions differ among various plant species. Contraction or expansion of the IR regions often results in length variation of cp genomes [24]. In this study, the detailed IR–SSC and IR-LSC borders, together with the adjacent genes, were compared across the five Rosaceae cp genomes (Figure 2). Of all the species except , the junction of LSC and IRa was located in the coding region of rps19. The end of IRa of expanded 178 bp into rps19, creating a truncated rps19 pseudogene of the same length at the IRb-LSC border. A similar feature was also observed in the other three cp genomes that the end of IRa expanded 152 bp, 95 bp, and 21 bp in , and , respectively. In , LSC region possessed an intact rps19 gene and expanded to incorporate 10 bp non-coding regions apart from the IRa-LSC border. On the boundary of IRa and SSC, IRa of all the species expanded to ycf1 and resulted in a truncated ycf1 gene with the same length as its homologous part in IRb. In , , and , a part of the ndhF gene was located in IRa; in and , the entire ndhF gene was located in SSC region with 22 bp and 93 bp apart from the IRa–SSC border, respectively. A similar pattern of expansion was found in IRb. In all five species, the initiation of IRb partially contained the ycf1 gene, which was terminal with the rps19 pseudogene. In , the nonfunctional rps19 gene seemed to be lost. Overall, a similar pattern of contraction/expansion was observed when these cp genomes approximately had a comparable genome size, for example, between and , and among the other three cp genomes.

Figure 2

Detailed view of the IR junctions among the five Rosaceae chloroplast genomes.

Annotated genes or portions of genes are indicated by yellow boxes above or below the genome.

Genome-wide sequence variations among Rosaceae cp genomes

We first aligned them using MAFFT [25] to identify genome-wide sequence variation among these five Rosaceae cp genomes (as shown in Figure 3). As previously studied [26-28], most of the sequence variations were located in LSC and SSC regions, while IR regions exhibited comparatively fewer sequence variations. Interestingly, a positive association between A/T content and sequence divergence was observed when base content was introduced into sequence alignment. Since IR region of cp genomes harbors lower A/T content, this might have led to higher sequence conservation of this region. Furthermore, the low A/T content of rRNA genes (above 50%) was also consistent with its high sequence conservation (data not shown). It has been well known that cp genomes of plants are quite A/T rich (above 60%) [28], and thus base content may be a major reason for determining cp genome variations.

Figure 3

Whole chloroplast genome alignment of the five Rosaceae species with content plotted along the sequences.

The alignment was wrapped corresponding to sequence similarity and GC/AT content was calculated in every 100 bp sequence window. The green lines indicate AT content, while the blue lines refer to GC content.

To look sequence divergence into details, we further plotted sequence divergence along the whole genomes with genes annotated (Figure S1). Consistent with other angiosperms, sequence divergence in intergenic regions was higher than that in genic regions of these five genomes. Intergenic regions with higher levels of divergence were matK–rps16, rps16–psbK, psbI–atpA, atpH–atpI, petN–psbM, ndhC–atpE, and rpl32–ccsA. Genic regions were ycf1, rpl22, and matK. The large sequence divergence in both coding and non-coding regions was mainly due to indels, as the sites where large indels (≥100 bp) occurred were mostly located in LSC or SSC regions. In addition, the IR–SSC and IR-LSC borders also harbored more sequence variations, which coincided with quite often contraction or expansion of the IR regions.

Simple sequence repeats (SSRs) of Rosaceae cp genomes

Cp SSRs have been proven to be highly valuable in plant population genetic [29,30] and phylogenetic studies [31,32]. In this study, we calculated the number of SSRs of together with that of the previously sequenced four Rosaceae cp genomes by using thresholds of mononucleotide repeats ≥ 8 nt, dinucleotide repeats ≥ 5 nt, and trinucleotide repeats ≥ 5 nt. We also calculated the number of longer repeats with thresholds of mononucleotide repeats ≥ 12 nt, dinucleotide repeats ≥ 10 nt, and trinucleotide repeats ≥ 10 nt [33]. There is no particular criterion about which nucleotide length or repeat unit may have a significant variation [34]. However, it has been suggested that SSRs longer than 8 nt tend to cause slip-strand mispairing during DNA replication [35], which is generally acknowledged as the major mechanism by which length variants are generated [36,37]. Based on the first criteria, we found 182 (168+14+0), 176 (159+17+0), 199 (180+19+0), 174 (166+8+0), 161 (147+13+1) 8, 5, 5 SSRs in , , , , and , respectively. While according to the second criteria, only 21, 21, 35, 19 and 9 mononucleotide repeats were observed in these five species, respectively (Table 3).

Table 3

Number of SSRs present in Rosaceae chloroplast genomes.

	8, 5, 5, SSRs (number)						12, 10, 10 SSRs (number)
Taxon	mono	di	tri	total	A or T mono. repeat	coding region	mono	di	tri	total	A or T mono. repeat	coding region
P . utilis	168	14	0	182	177	55	21	0	0	21	21	4
	(92%)	(8%)			(97%)	(31%)	(100%)				(100%)	(19%)
P. persica	159	17	0	176	171	47	21	0	0	21	21	2
	(90%)	(10%)			(97%)	(36%)	(100%)				(100%)	(10%)
P . pyfifolia	180	19	0	199	195	42	35	0	0	35	35	1
	(90%)	(10%)			(98%)	(21%)	(100%)				(100%)	(3%)
P . rupicola	166	8	0	174	160	44	19	0	0	19	19	1
	(95%)	(5%)			(96%)	(26%)	(100%)				(100%)	(5%)
F. vesca	147	13	1	161	152	35	7	2	0	9	9	0
	(91%)	(8%)	(1%)		(94%)	(21%)	(78%)	(22%)			(100%)

Among these cpSSR nucleotide units, the longest repeat stretch was polyA and (or) polyT with 16–17 bp, and most of the repeat units (> 95%) were also composed of polyA or polyT in these cp genomes; when the threshold was increased, all the repeat units were composed of polyA or polyT (Table 3). This was consistent with the observation in previous studies that cpSSRs were dominated by A or T mononucleotide repeats [38]. Furthermore, cpSSRs were unevenly distributed along the whole genome with about 70% in noncoding regions (Table 3). The composition and distribution of cpSSRs were quite similar within these Rosaceae cp genomes. Among these identified cpSSRs, about 80% were shared among these five Rosaceae cp genomes (similar repeat units located in similar genomic regions) (data not shown), suggesting the potential of these cpSSRs to be applied in further studies of more Rosaceae species.

Here, we found that most of the cpSSRs were A or T mononucleotide repeats, and they were mainly located in intergenic regions which harbored the vast majority of genomic variation. Previous studies also suggested that most repeats and indels were located in intergenic regions which were A + T-rich in cp genome [39]. Therefore, the observation that the genome divergence was associated with AT content in this study may indicate a common mutational pattern in the cp genomes of angiosperms.

Phylogenetic analysis based on the complete cp genome sequences

To explore phylogenetic relationship of with other major angiosperms clades, we employed protein-coding sequences to resolve relationships among these divergent angiosperm clades. Since our interest focused on phylogenetic relationship between Rosaceae and its close sister clades, with an especial emphasis on its relationship with other species of rosids, we included most of the angiosperms clades previously studied instead of all the sequenced cp genomes [20,40]. A dataset including 78 protein-coding genes comprised of 65,713 aligned nucleotide characters was finally used for phylogenetic analyses of 78 taxa.

Maximum likelihood (ML) analysis resulted in a single tree with -lnL = 833985.25 (Figure 4). Bootstrap analyses indicated that most of the nodes (67/74) were supported by values of 95% or greater and 57 of these had a bootstrap value of 100%. Maximum parsimony (MP) analysis also generated a single, fully resolved tree with a length of 74,348, a consistency index of 0.354 (excluding uninformative characters) and a retention index of 0.614 (Figure S2). Bootstrap analyses indicated that 58 of the 78 nodes were supported by values of 95% or greater, and 56 of these had a bootstrap value of 100%. Both the ML and MP trees produced a similar topology (Figure 4 and Figure S2) and one of major differences concerned the position of as described in earlier papers [20,40]. The major clades including basal angiosperms, monocots, eudicots were strongly supported in this study.

Figure 4

ML phylogram of the angiosperms using whole chloroplast genome sequences.

The tree has two outgroup taxa, and . Numbers above each node indicate the ML bootstrap support values.

In Rosaceae, both ML and MP analyses strongly supported that was sister to the rest of the four species. Besides, ML tree divided the rest four species into two sister clades with 100% bootstrap support. The first clade included and , and the second clade was consisted of the other two species. MP tree made a little difference that and were grouped together but and did not form one sister clade. By comparing with previous studies, ML tree might represent a more correct topology in constructing Rosaceae phylogeny by using these genes [1,2]. Thus, only the ML tree was included in our latter analyses.

To identify indel variations along the Rosaceae phylogeny, a total of 95 indels in the gene regions of which the largest gap size was over 110 bp were found and mapped to the cp genome-based phylogenetic tree using as an outgroup (Figure 5). These 95 indels were distributed in 23 different genes, three (rpoC2, ycf2, and ycf1) of which contained as many as 12, 14, and 27 indels, respectively. Among these 7 branches, the branch leading to contained most of the indels (21 of 95). The second most were contained by the common branch leading to two sister clades (19 of 95). The branch leading to the clade of and only harbored three indels. Up to now, there is no clear consensus about whether indels are useful or not in phylogenetic analyses [41,42], and the doubtfulness is essentially based upon studies in which only one or several DNA fragments were utilized. In this study, 33 of all the indels observed were unambiguously mapped to monophyletic groups and they might have intensive benefit for nucleotide-based phylogenetic analyses. The remained indels might be homoplasies and possibly had negative effects on the reconstruction of phylogenetic tree. Previous studies also suggest that indels should be cautiously used in phylogenetic studies especially when few DNA fragments were analyzed.

Figure 5

Phylogenetic distribution of exon coding indels in the five Rosaceae chloroplast genomes.

The phylogenetic tree was a subtree of the one in Figure 4. Gene names are given above boxes, sizes of indels (bp) and polarity (‘+’ = insertion, ‘-’ = deletion) are given below boxes, and numbers following * in parentheses indicate indels which appeared no less than twice in one gene. Polarity of mutations was determined by comparison to the outgroup, F. vesca.

Identification of chloroplast molecular markers

Since cp genomes of higher plants are quite conserved, the identification of molecular markers with higher sequence variation could be fairly valuable for population genetic and phylogenetic studies. Our genome-wide comparative analyses found that some genes boasted superior sequence divergence which could be definitely suitable for phylogeny studies (Figure S1). The values of overall mean distance (standing for genetic divergence) of 76 protein-coding genes in the five Rosaceae cp genomes, were further calculated (Figure S3). Among them, ycf1, one of the largest genes in the cp genomes (about 5.5 kb), harbored the highest sequence divergence. To explore its phylogenetic utility, Maximum likelihood (ML) analysis for ycf1 was conducted in 65 species resulting in a single tree with -lnL = 179121.15 and high bootstrap values (Figure S4). 37 of the 62 nodes possessed bootstrap values of 95% or greater and 31 of which had 100%. Despite the limited number of characters, overall topology of the phylogeny was quite similar to the one constructed with 78 protein-coding genes (Figure 4), indicating that ycf1 showed the potential for phylogenetic analysis in angiosperms. By comparing with another 9 potential molecular markers of which sequence divergence was no lower than 0.65 (rps18, rpl32, ccsA, ndhF, matK, rpl22; Figure S3) or indels were no less than 6 (matK, rpoC1, rpoC2, ycf2; Figure 5), ycf1 still performed the best to differentiate these five Rosaceae species (Figure S5). While other markers presented diverse topologies in reconstructing phylogeny of Rosaceae and with lower node supports, ycf1 showed the same topology as the subtree constructed by 78 protein-coding genes with comparatively higher node supports (Figure S5). Thus, ycf1 may serve as the best gene marker in distinguishing phylogenetic relationships within both Rosaceae and angiosperms.

Finding variable and easily aligned DNA sequences has been a challenging task in molecular genetics. Until now, most of DNA sequences used in plant phylogenetic analysis have come from the plastid genes, e.g. matK, rbcL, and trnL-F which have been extensively used. Although utility of these DNA regions has often obtained well-resolved phylogenies, the need for more data to address questions at multiple phylogenetic levels still persists. Previous studies indicated that ycf1 may be the most variable coding region in the cp genome, which is consistent with our analysis in the Rosaceae [43-45]. The usage of ycf1 has been confined to only orchids [43], Pinus [44], and magnoliid [45] where it has been demonstrated to be phylogenetically useful in analyses at the species, genus, and subfamily levels. Our study cans broad its uses as a new marker at least in Rosaceae.

Conclusion

In this study, we reported the complete chloroplast genome sequence of , a wide woody oil plant in Rosaceae. The plastome was characterized by comparing with other sequenced plastomes in Rosaceae and the evolutionary position of in Rosaceae was determined by phylogenetic analysis of concatenated alignments of gene sequences. While comparisons of plastomes in Rosaceae have identified many regions with higher sequence divergence, ycf1 was proved to boast the potential to solve phylogenetic issues ranging from lower taxonomic ranks to the relationships among genera and species. The availability of the complete cp genome of will provide new resources for the utilization of cp genome sequence to further evolutionary studies and genetic breeding programs in Rosaceae.

Materials and Methods

Ethics statement

The species is naturally grown in Kunming Botanical Garden of Kunming Institute of Botany, the Chinese Academy of Sciences. The voucher specimens were deposited at the Herbarium of Kunming Institute of Botany (KUN). Permit for collection of the plant leaves was obtained from Wei-bang Sun, Kunming Botanical Garden, the Chinese Academy of Sciences.

Chloroplast DNA Extraction, Genome Sequencing and Assembly

About 20g fresh mature leaves were sampled from a single plant. cpDNA was extracted by use of a modified high salt method reported by our group [46].

After the cpDNA isolation, approximately 5-10µg of DNA was sheared, followed by adapter ligation and library amplification, and then subjected to Illumina Sample Preparation Instructions. The fragmented cpDNAs were sequenced at both single-read using the Illumina Genome Analyzer IIx platform at the in-house facility at The Germplasm Bank of Wild Species in Southwestern China. Files containing the raw read sequences are available from the National Center for Biotechnology Information (NCBI) Short Read Archive with the accession number: SRX277333. The obtained paired-end reads (2×100 bp read lengths) were assembled using SOAP de novo [47] to the cp genome sequence of , which was downloaded from GenBank (NC_014697) and served as a reference genome. Regions with ambiguous alignment (conflicted reads mapped to the same genomic region) were trimmed off manually and considered as gaps. Polymerase chain reaction (PCR) amplified fragments yielded by primers (Table S1) derived from the terminal ends of contigs, and the fragments were then sequenced to flank the gap regions. The PCR amplification reactions were template denaturation at 80°C for 5 min followed by 30 cycles of denaturation at 95°C for 30 sec, primer annealing at 55°C for 30 sec, and primer extension at 65°C for 1 min; followed by a final extension step of 5 min at 65°C [48]. PCR products were separated by electrophoresis in 1.5% agarose gel and sequenced on an Applied Biosystems (ABI) 3730 automated sequencer. Finally, the average coverage depth reached approximately 855×, which is greatly higher than the proposed minimum for cpDNAs (30×) [49].

Genome Annotation

Annotation of the assembled genome was performed with Dual Organellar GenoMe Annotator (DOGMA) using default parameters to predict protein-coding genes, transfer RNA (tRNA) genes, and ribosome RNA (rRNA) genes [50]. BLASTX and BLASTN searches against a custom database of previously published cp genomes identified locations of these putative genes. For genes with low sequence identity, manual annotation was performed to determinate the position of start and stop codons depending on the translated amino acid sequence using the chloroplast/bacterial genetic code. Gene map of this plastome was drawn by OGDraw v 1.2 [51]. The plastid genome sequence has been deposited in GenBank (accession number KC571835).

Whole-genome Sequence Alignment and Detection of SSRs

mVISTA program (http://genome.lbl.gov/vista/index.shtml) was utilized for pair-wise whole-genome alignment of the five Rosaceae cp genomes with serving as an outgroup. This program takes both fasta format sequences and gene annotation text as input files, but here, only gene annotation was used. To further estimate sequence divergences, the common protein-coding genes were extracted from these five species and aligned in MEGA5 [52]. Gene divergences represented by overall mean distance were calculated with Jukes-cantor model with all gaps and missing data treated as complete deletion. SSRs were detected using Msatfinder on-line v.2.0 (http://www.genomics.ceh.ac.uk/cgi-bin/msatfinder/msatfinder.cgi). We calculated the number of SSRs by using two thresholds of 1) mononucleotide repeats ≥ 8 nt, dinucleotide repeats ≥ 5 nt, and trinucleotide repeats ≥ 5 nt; 2) mononucleotide repeats ≥ 12 nt, dinucleotide repeats ≥ 10 nt, and trinucleotide repeats ≥ 10 nt.

Phylogenetic Analyses

For the purpose of phylogenetic analyses, two data sets were utilized. The first set included 78 common protein-coding genes which were extracted from 76 angiosperm cp genomes and 2 gymnosperm outgroups, representing different plant lineages of which cp genomes were sequenced (see Table S2 for a complete list of the genomes examined). For all the species, every gene was translated into amino acid sequence separately, and then aligned in MSWAT (http://mswat.ccbb.utexas.edu/). This alignment was next constrained back to the nucleotide alignment. Some genes (e.g., infA, rpl22, and ) that got lost from the majority of species were excluded from the analysis. The second data set included the ycf1 gene sequence from the same taxa as the 78-gene data set except that those taxa without ycf1 were deleted. These sequences were aligned using the program MAFFT version 5 [25] and adjusted manually where necessary.

Maximum likelihood (ML) analyses were implemented in RAxML version 7.2.6 [53]. RAxML searches relied on the general time reversible (GTR) model of nucleotide substitution with the gamma model of rate heterogeneity. Non-parametric bootstrapping as implemented in the ‘‘fast bootstrap’’ algorithm of RAxML used 1,000 replicates. The aligned data matrices are available upon request. MP analyses were performed with PAUP*4.0b10. Heuristic tree searches were conducted with 1,000 random-taxon-addition replicates and tree bisection-reconnection (TBR) branch swapping, with ‘‘multrees’’ option in effect. Non-parametric bootstrap analysis was conducted under 1,000 replicates with TBR branch swapping.

Visualization of alignments among the five VISTA-based identity plots show sequence identity among the five sequenced chloroplast genomes with as a reference. Genome regions are color-coded as coding and non-coding regions.

(TIF)

Click here for additional data file.

MP phylogram of the angiosperms using whole chloroplast genome sequences. Numbers above each node indicate the MP bootstrap support values.

(TIF)

Click here for additional data file.

Gene divergences among the five The genes are oriented according to their locations in the chloroplast genome.

(TIF)

Click here for additional data file.

ML phylogram of the angiosperms using the Numbers at the nodes are ML bootstrap support values.

(TIF)

Click here for additional data file.

ML phylogram of the five Each tree was conducted using as outgroup. Numbers above each node indicate the ML bootstrap support values.

(TIF)

Click here for additional data file.

Primers used for validating genome assembly and gaps closing.

(XLS)

Click here for additional data file.

Taxa included in the phylogenetic analyses with GenBank accession numbers.

(DOC)

Click here for additional data file.