Background & Summary

Tea is the most popular non-alcoholic caffeine-containing and the oldest beverage in the world since 3000 B. C.1,2. The production of tea made from the young leaves of Camellia sinensis var. sinensis and C. sinensis var. assamica, together with ornamentally well-known camellias (e.g., C. japonica, C. reticulata and C. sasanqua) and worldwide renowned wooden oil crop C. oleifera3 has made the genus Camellia possess huge economic values in Theaceae. Besides its industrial, cultural and medicinal values, botanists and evolutionary biologists have increasingly paid attention to this genus. As a result of frequent hybridization and polyploidization, Camellia is almost commonly regarded as one of the most taxonomically and phylogenetically difficult taxa in flowering plants4. Thus, it has long been problematic for the taxonomic classification of the Camellia species based on the morphological characteristics5. The chloroplast (cp) genomes are able to provide valuable information for taxonomic classification, tracing source populations6,7 and the reconstruction of phylogeny to resolve complex evolutionary relationships8,9,10 due to the conservation of genomic structure, maternal inheritance and a fairly low recombination rate. Genetically speaking, cp genomes are comparatively conserved than plant mitochondria (mt) genomes which are more heterogeneous in nature. However, the presence of NUPT (nuclear plastid DNA) into cp genomes argues that cp genomes assembled from WGS data may include the heterogeneity due to the nuclear cp DNA transferred to the nucleus, resulting in erroneous phylogenetic inferences11. It has long been acknowledged that mtDNA has the propensity to integrate DNA from various sources through intracellular and horizontal transfer12,13,14. Partially due to these reasons, the mt genomes vary from ~200 Kbp to ~11.3 Mbp in some living organisms15,16,17. The dynamic nature of mt genome structure has been recognized, and plant mt genomes can have a variety of different genomic configurations due to the recombination and differences in repeat content18,19. These characteristics make the plant mt genome a fascinating genetic system to investigate questions related to evolutionary biology. The first effort has been made to sequence the 13 representative Camellia chloroplast genomes using next-generation Illumina genome sequencing platform, which obtained novel insights into global patterns of structural variation across the Camellia cp genomes4. The reconstruction of phylogenetic relationships among these representative species of Camellia suggests that cp genomic resources are able to provide useful data to help to understand their evolutionary relationships and classify the ‘difficult taxa’. Increasing interest in the Camellia plants have made up to thirty-eight of cp genomes be sequenced up to date20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37. Recently, we decoded the first nuclear genome of C. sinensis var. assamica cv. Yunkang10, providing novel insights into genomic basis of tea flavors38. Besides the lack of the C. sinensis var. assamica cp genome among thirty-eight cp genomes that were sequenced in this genus4,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37, up to data, none of mt genome has been determined in the genus Camellia.

In this study, we filtered cpDNA and mtDNA reads from the WGS genome sequence project38 and de novo assembled the mt genome and cp genome of C. sinensis var. assamica. The information of both cp and mt genomes will help to obtain a comprehensive understanding of the taxonomy and evolution of the genus Camellia. These genome sequences will also facilitate the genetic modification of these economically important plants, for example, through chloroplast genetic engineering technologies.

Methods

Plant materials, DNA extraction and genome sequencing

Young and healthy leaves of an individual plant of cultivar Yunkang10 of C. sinensis var. assamica were collected for genome sequencing in April, 2009, from Menghai County, Yunnan Province, China. Fresh leaves were harvested and immediately frozen in liquid nitrogen after collection, followed by the preservation at −80 °C in the laboratory prior to DNA extraction. High-quality genomic DNA was extracted from leaves using a modified CTAB method39. RNase A and proteinase K were separately used to remove RNA and protein contamination. The quality and quantity of the isolated DNA were separately checked by electrophoresis on a 0.8% agarose gel and a NanoDrop D-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE). A total of eleven paired-end libraries, including four types of small-insert libraries (180 bp, 260 bp, 300 bp, 500 bp) and seven large-insert libraries (2 Kb, 3 Kb, 4 Kb, 5 Kb, 6 Kb, 8 Kb, 20 Kb), were prepared following the Illumina’s instructions, and sequenced using Illumina HiSeq. 2000 platform by following the standard Illumina protocols (Illumina, San Diego, CA). We totally generated ~707.88 Gb (~229.31×) of raw sequencing data38. Further reads quality control filtering processes yielded a total of ~492.15 Gb (~159.43×) high-quality data retained and used for subsequent genome assembly.

De novo chloroplast and mitochondria genome assemblies

The chloroplast reads were filtered from whole genome Illumina sequencing data of C. sinensis var. assamica, we mapped all the sequencing reads to the reference genomes4 using bowtie2 (version 2.3.4.3)40. The mapped chloroplast reads were assembled into a circular contig of 157,100 bp in length with an overall GC content of 37.29% using CLC Genomics Workbench v. 3.6.1 (CLC Inc., Rarhus, Denmark) (Fig. 1). For mitochondria genome assembly, the PE and MP sequencing reads were used separately. Briefly, we first performed de novo assembly with VELVET v1.2.0841, which was previously described42,43. Scaffolds were constructed using SSPACE v.3.044. False connection was manually removed based on the coverage and distances of paired reads. Gaps between scaffolds were then filled with GapCloser (version 1.12)45,46 using all pair-end reads. We obtained the two complete circular scaffolds (701719 bp and 177329 bp) of the C. sinensis var. assamica mt genome from the de-novo assembly of the filtered mitochondrial reads (Figs 24). The two scaffolds of the mt genome had overall GC contents of 45.63% and 45.81%, respectively. The completed chloroplast and mitochondria genomes are publicly available in NCBI GenBank under accession numbers MH019307, MK574876 and MK574877 and BIG Genome Warehouse WGS000271, WGS000272.

Fig. 1
figure 1

Genome map of C. sinensis var. assamica cv. Yunkang10. 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.

Fig. 2
figure 2

The assembly and annotation pipeline of the tea tree mitochondrial genome.

Fig. 3
figure 3

Circular map of scaffold 1 in the C. sinensis var. assamica cv. Yunkang10 mitochondrial genome. Gene map showing 54 annotated genes with different functional groups that are color-coded on outer circle as transcribed clock-wise (outside) and transcribed counter clock-wise (inside). The inner circle indicates the GC content as dark grey plot.

Fig. 4
figure 4

Circular map of scaffold 2 in the C. sinensis var. assamica cv. Yunkang10 mitochondrial genome. Gene map showing 17 annotated genes with different functional groups that are color-coded on outer circle as transcribed clock-wise (outside) and transcribed counter clock-wise (inside). The inner circle indicates the GC content as dark grey plot.

Genome annotation and visualization

The complete chloroplast genome of C. sinensis var. assamica was preliminarily annotated using the online program DOGMA47 (Dual Organellar Genome Annotator) followed by manual correction. A total of 141 genes were annotated, of which 87 were protein-coding genes, 46 were tRNA genes and eight were rRNA genes (Table 1). MITOFY15 was used to characterize the complement of protein-coding and rRNA genes in the mitochondrial genome. A tRNA gene search was carried out using the tRNA scan-SE software (version 1.3.1)48. We annotated a total of 71 genes, including 44 protein-coding genes, 24 tRNAs and 3 rRNAs (Table 2). Circular genome maps were drawn with OrganellarGenomeDRAW49 (Figs 34).

Table 1 Gene annotation of the C. sinensis var. assamica cp genome.
Table 2 Gene content of the C. sinensis var. assamica mt genome.

Simple sequence repeats (SSRs) were identified and located using MISA (http://pgrc.ipk-gatersleben.de/misa/). All the annotated SSRs were classified by the size and copy number of their tandemly repeated: monomer (one nucleotide, n ≥ 8), dimer (two nucleotides, n ≥ 4), trimer (three nucleotides, n ≥ 4), tetramer (four nucleotides, n ≥ 3), pentamer (five nucleotides, n ≥ 3), hexamer (six nucleotides, n ≥ 3). A total of 214 SSRs were identified in cp genome with 74.42% of which were monomers, 19.07% of dimers, 0.47% of trimers, 4.65% of tetramers and 0.93% of hexamers (Table 3). There were no pentamers found in the cp genome. In mt genome, we obtained 665 SSRs distributed into monomers, dimers, trimers, pentamers, tetramers and hexamers with 31.53%, 45.35%, 4.95%, 15.17%, 2.70% and 0.15%, respectively (Table 3). Repeat sequences including forward and palindromic repeats, were also searched by REPuter50 with the following parameters: minimal length 50 nt; mismatch 3 nt. Long repeat sequences (repeat unit > 50 bp) of forward and palindromic repeats were further annotated, resulting in 149 bp from 4 paired repeats in the cp genome (Table 4) and 37,878 bp from 58 paired repeats in the mt genome (Online-only Tables 12). Our repeat content analyses indicate that the mt genome is more abundant in repeat sequences and more variable than the cp genome of C. sinensis var. assamica (Table 4; Online-only Tables 12).

Table 3 Statistics of SSR motifs in the C. sinensis var. assamica mt and cp genomes.
Table 4 Long repeats (repeat unit > 50 bp) in the C. sinensis var. assamica cp genome.

Prediction of RNA-editing sites

Putative RNA editing sites in protein-coding genes were predicted using the PREP-cp and PREP-mt Web-based program (http://prep.unl.edu/)51,52. To achieve a balanced trade-off between the number of false positive and false negative sites, the cutoff score (C-value) was set to 0.8 and 0.6, respectively53.

Almost all transcripts of protein encoding genes in the plant mitochondria are subject to RNA editing except the T-urf13 gene54. Our results showed that the extent of RNA editing varied by gene for both cp and mt genomes of C. sinensis var. assamica. In the C. sinensis var. assamica cp genome, we detected 54 RNA-editing sites in 21 protein-coding genes, ranging from one editing site in atpF, atpI, petB, psaI, psbE, psbF, rpoA, rps2 and rps8 to 8 editing sites in ndhB (Online-only Table 3). In the C. sinensis var. assamica mt genome, we predicted 478 RNA-editing sites in 42 protein-coding genes; they varied from two editing site in atp9 (of scaffold2), sdh3 (of scaffold1 and scaffold2, respectively) and rps14 (of scaffold2) to 35 editing sites in ccmFn (of scaffold1) (Online-only Table 45).

Phylogenetic analyses

To further determine the phylogenetic position of C. sinensis var. assamica we performed phylogenomic analysis of 20 complete cp genomes using the GTR + R + I model under the maximum likelihood (ML) inference in MEGA v.7.055. Besides C. sinensis var. assamica cv. Yunkang 10, we selected cp genomes from the eighteen Camelia species (C. oleifera, C. crapnelliana, C. szechuanensis, C. mairei, C. elongata, C. grandibracteata, C. leptophylla, C. petelotii, C. pubicosta, C. reticulata, C. azalea, C. japonica, C. cuspidata, C. danzaiensis, C. impressinervis, C. pitardii, C. yunnanensis and C. taliensis) using Apterosperm oblata as outgroup. Our results showed that C. sinensis var. assamica was grouped with C. grandibracteata with 100% bootstrap support (Fig. 5).

Fig. 5
figure 5

Phylogenetic relationships of 20 complete chloroplast genomes. Maximum likelihood phylogenetic tree of C. sinensis var. assamica cv. Yunkang 10 with 18 species in the genus Camellia based on complete chloroplast genome sequences. The chloroplast sequence of Apterosperma oblata was set as outgroup. The position of C. sinensis var. assamica cv. Yunkang 10 is shown in bold and bootstrap values are shown for each node.

The same method was used for phylogenetic analysis with mt genome. A total of thirteen conserved mt protein-coding genes among C. sinensis var. assamica and 14 other plant species were individually aligned with ClustalW56, and then concatenated to construct a contiguous sequence in the order of cob, cox1, cox2, cox3, nad1, nad2, nad3, nad4, nad4L, nad5, nad6, nad7 and nad9. The selected 14 species includes Cycas taitungensis, Ginkgo biloba, Triticum aestivum, Oryza sativa, Sorghum bicolor, Zea mays, Gossypium arboretum, G. barbadense, Carica papaya, Vitis vinifera, Hevea brasiliensis, Bupleurum falcatum, Glycine max and Salvia miltiorrhiza. The alignment file was used for the construction of Neighbor-Joining Tree at 1000 bootstrap replicates with MEGA 7.0.2655. Our results showed that C. sinensis var. assamica is clearly grouped with other dicots that were separated from monocots of the angiosperms while the two gymnosperms (Cycas taitungensis and Ginkgo biloba) were formed the basal clade (Fig. 6).

Fig. 6
figure 6

Phylogeny inferred from 13 genes common in the 15 plant mitochondrial genomes. Neighbor-joining tree of C. sinensis var. assamica cv. Yunkang 10 with other 14 species based on 13 conserved protein-coding gene sequences with bootstrap support values on each node. The mt sequence of Cycas taitungensis and Ginkgo biloba were set as outgroup.

Data Records

Raw reads from Illumina are deposited in the NCBI Sequence Read Archive (SRA)57,58,59,60,61,62 and BIG Genome Warehouse63. Assembled cp genome sequences and accompanying gene annotations of C. sinensis var. assamica are deposited in the NCBI GenBank64 and BIG Genome Warehouse65. The mt genome final assembly and accompanying gene annotations are deposited at NCBI GenBank66,67 and BIG Genome Warehouse68. The alignment and tree files of the chloroplast genome and mitochondrial genome form the Camellia genus were deposited in Figshare database69.

Technical Validation

Quality filtering of raw reads

The initially generated raw sequencing reads were evaluated in terms of the average quality score at each position, GC content distribution, quality distribution, base composition, and other metrics. Furthermore, the sequencing reads with low quality were also filtered out before the genome assembly and annotation of gene structure.

Assembly and validation

The chloroplast reads were filtered from whole genome Illumina sequencing data of C. sinensis var. assamica. We mapped all the cleaned reads to the reference chloroplast sequence4 using bowtie2 (version 2.3.4.3)40 with default parameters. The mapped chloroplast reads were de novo assembled into the complete chloroplast genome.

For mitochondria genome assembly, the PE and MP sequencing reads were used separately. Briefly, we first performed de novo assembly with VELVET v1.2.0841, which was previously described42,43. Scaffolds were constructed using SSPACE v.3.044. False connection was manually removed based on the coverage and distances of paired reads. Gaps between scaffolds were then filled with GapCloser (version 1.12)45,46 using all pair-end reads.