Leaf transcriptome analysis of a subtropical evergreen broadleaf plant, wild oil-tea camellia (Camellia oleifera), revealing candidate genes for cold acclimation
© The Author(s). 2017
Received: 7 July 2016
Accepted: 9 February 2017
Published: 28 February 2017
Cold tolerance is a key determinant of the geographical distribution range of a plant species and crop production. Cold acclimation can enhance freezing-tolerance of plant species through a period of exposure to low nonfreezing temperatures. As a subtropical evergreen broadleaf plant, oil-tea camellia demonstrates a relatively strong tolerance to freezing temperatures. Moreover, wild oil-tea camellia is an essential genetic resource for the breeding of cultivated oil-tea camellia, one of the four major woody oil crops in the world. The aims of our study are to identify variations in transcriptomes of wild oil-tea camellia from different latitudes and elevations, and discover candidate genes for cold acclimation.
Leaf transcriptomes were obtained of wild oil-tea camellia from different elevations in Lu and Jinggang Mountains, China. Huge amounts of simple sequence repeats (SSRs), single-nucleotide polymorphisms (SNPs) and insertion/deletions (InDels) were identified. Based on SNPs, phylogenetic analysis was performed to detect genetic structure. Wild oil-tea camellia samples were genetically differentiated mainly between latitudes (between Lu and Jinggang Mountains) and then among elevations (within Lu or Jinggang Mountain). Gene expression patterns of wild oil-tea camellia samples were compared among different air temperatures, and differentially expressed genes (DEGs) were discovered. When air temperatures were below 10 °C, gene expression patterns changed dramatically and majority of the DEGs were up-regulated at low temperatures. More DEGs concerned with cold acclimation were detected at 2 °C than at 5 °C, and a putative C-repeat binding factor (CBF) gene was significantly up-regulated only at 2 °C, suggesting a stronger cold stress at 2 °C. We developed a new method for identifying significant functional groups of DEGs. Among the DEGs, transmembrane transporter genes were found to be predominant and many of them encoded transmembrane sugar transporters.
Our study provides one of the largest transcriptome dataset in the genus Camellia. Wild oil-tea camellia populations were genetically differentiated between latitudes. It may undergo cold acclimation when air temperatures are below 10 °C. Candidate genes for cold acclimation may be predominantly involved in transmembrane transporter activities.
Cold tolerance is a key determinant of the geographical distribution range of a plant species and crop production [1–3]. Many plant species demonstrate an increase in freezing tolerance upon a period of exposure to low nonfreezing temperatures, a phenomenon known as cold acclimation . In nature, the process of cold acclimation helps plants to prepare for the coming of winter so as to survive under seasonal freezing temperatures. On the other hand, freezing damage on crops can lead to serious loss of crop production . Therefore, intensive researches have been conducted in many plant species for understanding the molecular mechanisms of cold acclimation [1, 2]. The majority of the previous studies focused on herbaceous plant species (e.g. Arabidopsis). A review by Wisniewski et al.  indicated that the molecular mechanisms of cold acclimation in woody plants were more complex than in herbaceous plants, and it was still not clear about what made a woody plant more cold hardy than an annual, herbaceous plant.
As an evergreen broadleaf shrub or small tree, oil-tea camellia (Camellia oleifera) is one of the representative plant species in subtropical evergreen broadleaf forests  with relatively strong tolerance to cold climates. Oil-tea camellia is widely distributed in the subtropical mountain areas of the Yangtze River basin and South China, with elevation ranging from about 200 to 2000 m [5, 6]. The northern range of oil-tea camellia is located in the mountain areas of the north subtropical region in China, where the mean annual air temperature is 14–16 °C, the mean January air temperature is 0–4 °C, and the minimum air temperature is as low as −17 °C . It has been reported that oil-tea camellia could survive under −26 °C in the USA, and was used to cross with C. sasanqua and C. hiemalis (susceptible to winter injury) for producing ornamental camellia varieties with cold tolerance [7, 8]. As an evergreen broadleaf plant species, oil-tea camellia has green leaves even in cold winter. Unlike most of the flowering plant species in the world, it flowers in autumn and winter. The fatty acid contents in oil-tea camellia seeds showed significant correlations with latitudes , which may be due to natural selection on seed germination temperature . Therefore, oil-tea camellia can be used as a model to study the molecular basis of cold tolerance in evergreen broadleaf plant species. However, the candidate genes related to cold acclimation and the gene expression patterns are still unknown in oil-tea camellia.
Cultivated oil-tea camellia is regarded as one of the world’s four major woody oil crops together with oil palm, coconut and oil olive, and is the top one woody oil crops in China . The utilization of oil-tea camellia seed oil (camellia oil) as cooking oil has a history of more than 2300 years in China . Camellia oil is rich in unsaturated fatty acids (more than 80% of total oil content), containing mainly monounsaturated fatty acid (i.e. oleic acid, contributing to more than 68% of total oil content) and some polyunsaturated fatty acid (i.e. linoleic acid and linolenic acid) [9, 11]. Its fatty acid composition is similar to olive oil, and it is therefore known as “oriental olive oil” [11, 12]. Camellia oil also contains other functional components such as camellia saponin, tea polyphenol and squalene . It has been shown that the intake of camellia oil is good for health, for instance, helping to reduce blood lipid and prevent cardiovascular diseases . Currently, China has about 3 million hectare cultivated oil-tea camellia, producing about 0.26 million ton camellia oil per year . To meet the rapidly increasing demands for healthy vegetable oil, the Chinese government plans to increase the cultivation of oil-tea camellia to more than 4 million hectare by 2020, with a yearly camellia oil production up to 2.5 million ton . The key issues for the development of oil-tea camellia cultivation are how to accelerate the breeding processes of varieties suitable for various regions, increase the yield and quality of camellia oil, and improve the resistance to diseases and pests.
Crop wild relatives are valuable genetic resources for crop breeding, for instance, helping to improve disease and pest resistances, and increase yield and quality of crops . Wild oil-tea camellia (C. oleifera) is an essential genetic resource for cultivated oil-tea camellia breeding. However, the patterns of genetic differentiation along latitude and elevation gradients in wild oil-tea camellia are still unknown, which is the basis for the utilization of wild oil-tea camellia resources. Currently, the evaluation of oil-tea camellia genetic resource for selective breeding is based on phenotypic traits. Due to the complex interactions between genotype and environment, phenotypic traits may not reflex the actual level of genetic variation in a population. On the other hand, as a perennial woody plant, oil-tea camellia has a juvenile phase of about 5 years . Many valuable economic traits need to be evaluated in the adult phase, such as fruit and seed characteristics, leading to a traditionally long breeding process of cultivated oil-tea camellia. Therefore, a huge amount of molecular markers shall be developed and applied for the purposes of genetic resource evaluation and marker-assisted breeding, so as to dramatically accelerate the breeding processes. Xia et al.  published the first transcriptome sequencing dataset of oil-tea camellia, providing the major genetic information of this species in public databases. However, their study used samples from a single oil-tea camellia individual in the botanic garden, and so the genetic variations were underestimated and could not represent the genetic differentiation along latitude and elevation gradients in nature .
Our study sequenced the leaf transcriptomes of wild oil-tea camellia from different latitudes and elevations, and analyzed the variations in gene sequences and gene expressions. The objectives of our study were to: 1) obtain the leaf transcriptomes of wild oil-tea camellia for functional genomics studies; 2) detect simple sequence repeats (SSRs), single-nucleotide polymorphisms (SNPs) and insertion/deletions (InDels) suitable for analyzing genetic differentiations along latitude and elevation gradients in wild oil-tea camellia; and 3) compare the gene expression patterns among different temperatures in leaves of wild oil-tea camellia and discover the candidate genes related to cold acclimation.
Study sites and sampling
Wild oil-tea camellia samples of different latitudes and elevations for the transcriptome sequencing
In each mountain, leaf samples were collected from flowering wild oil-tea camellia at different elevations within 3 h in the afternoon. Three to five fresh leaves without obvious damage were randomly picked from each plant, covered by aluminum foil and immediately placed in a vacuum bottle with liquid nitrogen. Latitude, longitude and elevation of each sampling plant were recorded. At the same time, air temperature next to each sampling plant was measured. All samples were stored at −80 °C in the lab.
RNA extraction and transcriptome sequencing
Each leaf was mixed with liquid nitrogen and ground into fine powder. About 100 mg tissue powder of each leaf was used for RNA extraction. Total RNAs were extracted using the EASYspin Plus Plant RNA Kit (Aidlab, Beijing, China). To account for the gene expression variations among leaves of the same plant, the RNAs from two leaves of the same plant were equimolarly pooled and used as a single sample for the transcriptome sequencing. In total, eight samples were used for the transcriptome sequencing (Table 1). According to the differences in air temperature, the samples could be divided into five temperature groups: T2, T5, T10, T14 and T18 (Table 1).
The cDNA libraries were constructed from RNA samples for Illumina paired-end (PE) sequencing following the Illumina protocol. PE sequencing (2 × 100 bp) was carried out on the Illumina HiSeq 2000 platform (Illumina, San Diego CA, USA) at Novogene Bioinformatics Technology Co., Ltd (Beijing, China).
Sequence assembly and Unigene annotation
Raw reads were processed to remove reads containing adaptors, with more than 10% ambiguous bases (N), or of low quality (more than 50% bases with small Qphred ≤ 5). All the downstream analyses were based on the resulting clean reads. Clean reads were assembled using Trinity (version r2012-10-05) with min_kmer_cov = 2 . The longest assembled transcript of a gene was taken as a unigene. All the assembled unigenes were used as reference sequences for the leaf transcriptome of wild oil-tea camellia.
Functional annotations of unigenes were based on the following databases: Nr (NCBI non-redundant protein sequences), Nt (NCBI non-redundant nucleotide sequences), Pfam (Protein family: http://pfam.xfam.org/) , KOG (euKaryotic Ortholog Groups)/COG (Clusters of Orthologous Groups of proteins) (http://www.ncbi.nlm.nih.gov/COG/) , Swiss-Prot (a manually annotated and reviewed protein sequence database: http://www.ebi.ac.uk/uniprot) , KEGG (Kyoto Encyclopedia of Genes and Genomes: http://www.genome.jp/kegg/) , and GO (Gene Ontology: http://geneontology.org/). NCBI blast 2.2.28+ was used for the alignments of unigenes to Nr, Nt, Swiss-Prot, and KOG. The E-value threshold was set to 1E − 5 in the alignments to Nr, Nt, and Swiss-Prot. For the alignments to KOG, the E-value threshold was 1E − 3. The hmmscan in HMMER 3.0 was used to search Pfam . The GO annotations were performed with Blast2GO v2.5  based on the Nr and Pfam annotations. KAAS (KEGG Automatic Annotation Server: http://www.genome.jp/kegg/kaas/) was used for the KEGG annotations .
Detection of SSRs, SNPs and InDels
MISA 1.0 was used to detect SSRs in unigenes. The minimum repeat number for unit size of mono-, di-, tri-, tetra-, penta-, and hexanucleotide was set to 10, 6, 5, 5, 5, and 5, respectively. Primer3 (2.3.5) was used to design primers around SSRs with default settings.
Clean reads of each sample were aligned to the reference sequences (unigenes) using bowtie 2 (mismatch 0) . The alignments were processed with SAMtools  and Picard tools for sequence sorting and duplicate removing. SNP and InDel callings were then performed using GATK2 . Those with QUAL < 30.0 and QD < 5.0 were removed.
Genetic structure analysis
In order to illustrate the genetic structure of wild oil-tea camellia samples from different latitudes and elevations, phylogenetic analysis was carried out based on the SNP data. SNP positions were chosen with no more than 2 alleles and the number of reads per sample ≥ 6. Then, an R  script was written to genotype each sample at each SNP position using IUPAC (International Union of Pure and Applied Chemistry) nucleic acid codes. Using the SNP genotypes of different samples, the Bayesian estimation of phylogeny was performed in MrBayes 3.2.5 by sampling across the entire general time reversible (GTR) model space . The resulted consensus tree was viewed and edited in FigTree v1.4.2.
Gene expression analysis
RSEM  was used to calculate the read count of each unigene in a sample and transform it into FPKM (expected number of fragments per kilobase of transcript per million fragments mapped) . In RSEM analysis, bowtie was used with mismatch 2. The resulting FPKM values were used to represent the gene expression levels of unigenes in different samples. To examine whether sequencing depth was sufficient for gene expression analysis, varied percentages of mapped reads were randomly taken from each sample and fraction of genes with an expression level within 10% of the final expression value (according to 100% mapped reads) was calculated. A curve illustrating the relationship between fraction of genes within 10% of the final value and percentage of mapped reads was made for each sample. If a curve became flat (saturation) with the increase in percentage of mapped reads, the sequencing depth should be sufficient for gene expression analysis.
To examine the effects of air temperature on gene expression patterns of leaves, density distributions of FPKM values were compared among different temperature groups. Differential gene expression analyses between different temperature groups were performed using DESeq . The threshold adjusted p-value for significance was 0.05. Hierarchical clustering and Venn diagram were used to illustrate the differential gene expression patterns between different temperature groups. According to the functional annotations of unigenes, putative functions of the differentially expressed genes (DEGs) were inferred to discover candidate genes for cold acclimation.
A new method was developed to figure out the major functional groups of genes involving in cold acclimation by comparing the GO annotations of DEGs with the GO annotations of all expressed genes detected in our study. To account for the effects of random sampling, all expressed genes were randomly sampled for 10000 times with a size equaling to the number of DEGs. Then, the number of DEGs in a GO functional group was tested for whether it was significantly different from the number of genes in the same GO functional group resulted from the random sampling of all expressed genes. The significant level was adjusted using the Bonferroni correction (0.05 divided by number of tests). An R script was written and used for the random sampling and statistical analysis.
Quantitative real-time PCR analysis
In order to validate the DEGs identified from transcriptome sequencing, quantitative real-time PCR (qRT-PCR) analysis was performed. Independent wild oil-tea camellia leaf samples at different air temperatures were used: JG05 at 17.8 °C and JG06 at 14.7 °C from Jinggang Mountain; LS05 at 11.0 °C and LS06 at 4 °C from Lu Mountain. Total RNAs were extracted as described before for the samples used in transcriptome sequencing. Using the PrimeScript™ RT reagent qPCR Kit with gDNA Eraser (Takara, Dalian, China), genomic DNA was removed from total RNAs (300 ng RNAs of each sample) and cDNA was synthesized. The PCR mixture contained 12.5 μL SYBR® Premix Ex Taq™ II (Tli RNaseH Plus) (Takara, Dalian, China), 9.5 μL ddH2O, 1 μL of each gene-specific primer (10 μM) and 1 μL cDNA template. The qRT-PCR assays were performed in a CFX96 Touch™ RT-PCR Detection System (Bio-RAD, USA) with the following program: 94 °C for 2 min; 40 cycles of 94 °C for 20 s, 57 °C for 20 s and 72 °C for 30 s. A commonly used reference gene, glyceraldehyde-3-phosphate dehydrogenase (GAPDH) gene, was used to normalize the expression levels of target genes . The relative expression levels of target genes were calculated with the 2−∆∆Cq method .
Summary of sequences and assembly
Summary of the sequencing data from different wild oil-tea camellia samples
No. of clean reads
Clean bases (Gb)
Error rate (%)a
GC content (%)
Functional annotation of unigenes
Annotation of unigenes in different databases
No. of annotated unigenes
Percentage of annotated unigenes (%)
At least onea
Detection of SSRs, SNPs and InDels
We discovered 661280 SNPs. About 54.3% of the SNPs were non-coding SNPs. For the coding SNPs, the ratio of non-synonymous to synonymous SNPs was 0.604, indicating most SNPs were synonymous. There were 103442 SNP positions with 2 alleles and number of reads ≥ 6 per sample (Additional file 2: Table S2). Genes containing the SNPs and ratio of non-synonymous to synonymous SNPs in each gene were summarized in Additional file 3: Table S3. Such data can help to develop SNP markers for oil-tea camellia.
We discovered 47056 InDels. For the purposes of molecular marker development, a gene containing only one InDel with two alleles was chosen and there were 6534 such InDels in total (Additional file 4: Table S4).
Differential gene expression
According to the relationships between fraction of genes within 10% of the final expression value (based on 100% mapped reads) and percentage of mapped reads, the curves became saturation when FPKM >3 for all samples (Additional file 5: Figure S1). Such results indicated that the sequencing depth was sufficient for gene expression analysis.
Significant functional groups of differentially expressed genes (DEGs) at T5/T2 versus T10/T14/T18. All products of the DEGs are integral components of membranes with transmembrane transporter activities
DEGs at T5a
DEGs at T2
comp200120_c0, comp203054_c0, comp206804_c0, comp209071_c1, comp212144_c0, comp212939_c0, comp214280_c0, comp215715_c0, comp216420_c0, comp218213_c0, comp220377_c0
comp183939_c0, comp207638_c0, comp209330_c0, comp214001_c0, comp214280_c0, comp215310_c0, comp218213_c0, comp220377_c0
Amino acid transporter
comp193115_c1, comp199501_c0, comp217974_c0
Small conductance calcium-activated potassium channel
Dicarboxylic acid transmembrane transporter
ATPase coupled transporter
ATP synthase coupled hydrogen ion transporter
Quantitative real-time PCR analysis
Leaf transcriptome of wild oil-tea camellia
Our study obtained a high-quality and large transcriptome dataset of wild oil-tea camellia, representing the leaf transcriptome variation from different latitudes and elevations. The amount of data obtained (57.3 Gb), the number (177258) and the total length (91.6 Mb) of unigenes assembled are all much larger than those reported in a previous study on the transcriptome of oil-tea camellia (170 Mb data and 104842 unigenes with a total length of 38.9 Mb) by Xia et al. (2014) . Such differences may be mainly due to the fact that the Illumina Hiseq 2000 high-throughput sequencing platform used in our study can produce a much larger amount of sequencing data than the 454 GS-FLX sequencing platform applied in the previous study. In our study, 83352 unigenes were annotated (Table 3), and 1504 unigenes were found to be involved in the lipid metabolism pathways (Fig. 1), which can serve as a basis for understanding the processes of fatty-acid biosynthesis in this important oil crop. Moreover, our study used more diverse samples of wild oil-tea camellia from various latitudes and elevations, and therefore the SSRs, SNPs and InDels identified in our study may be more useful for developing molecular markers to detect the genetic structure of wild oil-tea camellia along latitude and elevation gradients.
The genus Camellia has about 120 species. Besides of oil-tea camellia, the genus Camellia has many other important economic species. Tea plant, C. sinensis, is one of the most important economic crops, generating the most popular non-alcoholic beverage in the world. Tea plant is also the best-studied species in the genus Camellia with the largest amount of genomic data available in public databases. A deep transcriptome sequencing of tea plant was published in 2011, where 2.32 Gb high-quality data were generated and assembled into 127094 unigenes with a total length of 45.1 Mb . Wang et al.  reported the leaf transcriptomes of tea plant in response to cold acclimation. They obtained 4.96 Gb high-quality sequencing data and assembled 216831 non-redundant transcript sequences with a total length of 77.4 Mb . After combining all available transcriptome data of the tea plant, they got 282395 non-redundant transcript sequences with a total length of 94.7 Mb. Therefore, the transcriptome sequencing data of oil-tea camellia obtained in our study also provide one of the largest transcriptome dataset in the genus Camellia. Such a large transcriptome dataset of oil-tea camellia can facilitate the functional genomic studies and the molecular breeding of many other economically important Camellia species, especially for those more closely related to oil-tea camellia (subgenus Camellia) than tea plant (subgenus Thea), such as the famous ornamental plants C. japonica and C. sasanqua.
Genetic structure of wild oil-tea camellia
Based on 90000 SNPs, a phylogenetic tree was constructed with the wild oil-tea camellia samples from different latitudes and elevations. The wild oil-tea camellia samples seem to be separated mainly by latitudes (between Lu and Jinggang Mountains) and then by elevations (within Lu or Jinggang Mountain) (Fig. 3). Previous studies indicated that air temperature was one of the major factors affecting the growth and development of oil-tea camellia, which may lead to genetic differentiation along air temperature gradients . Because air temperature decreases 1 °C with an increase of about 167 m in elevation or about 145 km in latitude, gene flow may be more restricted between different latitudes than between different elevations along the same air temperature gradients . Thus, genetic differentiation may become more distinct between different latitudes than between different elevations . However, the genetic differentiation of wild oil-tea camellia between Lu and Jinggang Mountains is incomplete (Fig. 3). Wild oil-tea camellia sample LS02 from Lu Mountain is grouped together with JG01 from Jinggang Mountain. Such results suggest that wild oil-tea camellia from the two mountains may be isolated recently so that the effects of genetic drift and natural selection could not completely drive them apart to form two distinct branches. Although the genetic structure was constructed on a huge amount of SNPs, only eight wild oil-tea camellia samples from two mountains were used in our study. A large number of wild oil-tea camellia samples from various latitudes and elevations should be used to exactly identify the genetic structure in future studies.
Differential gene expressions under cold acclimation
Tropical or subtropical plants may experience chilling stress when air temperatures fall below 10 °C . A previous study also showed that tea plant (C. sinensis) underwent cold acclimation when air temperatures were below 10 °C . Many genes were differentially expressed in tea plant leaves during the cold acclimation process, and the number of up-regulated genes was about twice the number of down-regulated genes . In our study, the gene expression patterns in the leaves of wild oil-tea camellia obviously changed when air temperatures were below 10 °C (Fig. 4), suggesting that oil-tea camellia may also undergo cold acclimation when air temperature was below 10 °C. Again, almost all the DEGs at T2 and T5 were significantly up-regulated, showing that most of the genes were activated instead of being repressed during the cold acclimation process as also indicated in the above study of tea plant. In addition, more genes in oil-tea camellia leaves were differentially expressed at T2 than at T5, suggesting a stronger cold stress at T2. This may be due to the fact that wild oil-tea camellia underwent a longer period at lower air temperatures in the habitats of higher elevations.
Among the DEGs in cold acclimation, transmembrane transporter genes were found to be predominant in oil-tea camellia leaves. Many of the genes encode transmembrane sugar transporters (Table 4). Researches on tea plants, another Camellia species, also found that sugar transporter genes were differentially expressed in leaves during the cold acclimation process and may lead to sugar accumulation in cells [34, 36]. During the cold acclimation process, the soluble sugar content was found to be constantly elevated in leaves of tea plants . In addition, amino acid transporters were found to be significantly up-regulated in oil-tea camellia leaves at both T2 and T5, which may lead to amino acid accumulation in cells (Table 4). The accumulation of soluble sugars and amino acids in cells is thought to be one of the most predominantly metabolic changes in many plant species during cold acclimation [4, 36–39]. The massive increase of such solutes may help to adjust osmotic pressure in cells to decrease freezing temperature of cell sap and protect membranes against freeze-induced damages . Further metabolic studies are needed to find out whether or not soluble sugars and amino acids do increase in oil-tea camellia leaves during cold acclimation.
On the other hand, five transporter genes were found to be significantly up-regulated in oil-tea camellia leaves only at T2 (Table 4). Among the products of these genes, putative sodium/calcium exchanger and small conductance calcium-activated potassium channel may be related to the calcium signaling pathway. Calcium is an important second messenger in the low-temperature signal transduction pathway regulating the cold-acclimation response . Most of the genes in the calcium signaling pathway were also found to be up-regulated in tea plants during cold acclimation . In addition, putative sodium/calcium exchanger and dicarboxylic acid transmembrane transporter were found to increase in leaf tonoplasts of Arabidopsis thaliana upon cold acclimation . Therefore, the transporter genes only up-regulated in oil-tea camellia leaves at T2 may also be owing to the effects of cold acclimation. Moreover, a putative C-repeat binding factor (CBF) gene was significantly up-regulated in oil-tea camellia leaves only at T2. The CBFs are the major transcriptome factors regulating the expression of a large number of genes in cold acclimation . The CBF gene expression can be activated in about 15 min after transferring plants to low temperature (4 °C) . The significant up-regulation of genes concerned with cold acclimation in oil-tea camellia leaves only at T2 suggests that the air temperature at T5 may not be sufficiently low to remarkably activate the expressions of those genes.
Differential gene expression can be due to differences in genotype, environmental condition and the interaction between them. According to the genetic structure of wild oil-tea camellia (Fig. 3), most of the genetic differentiation occurred between Lu and Jinggang Mountains. However, the gene expression patterns could not be well explained by such a genetic structure. Sharp changes were observed in gene expressions of wild oil-tea camellia samples within the Lu Mountain instead of between the two mountains, and the gene expression patterns of two samples (T10: LS01 and LS02) from the Lu Mountain was similar to those (T14 and T18) from the Jinggang Mountain (Figs. 4 and 5). Meanwhile, the gene expression patterns could be well explained by the differences in air temperature. In general, when air temperature decreased to 2–5 °C (T2 and T5) from around 10–18 °C (T10, T14 and T18), many genes had increased expression levels. To minimize the random effects on gene expression of environmental factors other than air temperature and the interactions with genotypes, we selected genes showing consistently differential expression patterns between low (T2 or T5) and high (T10, T14 or T18) temperatures. The distributions of these DEGs in different functional groups are quite similar between T2 and T5 (Fig. 7), demonstrating that the differential gene expression patterns are not due to random effects. The predominantly functional groups of the DEGs were determined, and most of the genes were also reported to be differentially expressed during cold acclimation in previous studies. The qRT-PCR analysis using independent wild oil-tea camellia leaf samples at different air temperatures showed that the relative expression levels of the selected sugar transporter genes dramatically increased at 4 °C (Fig. 8). Such results demonstrated that the DEGs identified were reliable. Thus, we can conclude that the differential gene expressions in wild oil-tea camellia leaves observed in our study may be mainly due to the differences in air temperature. The DEGs identified may represent the major candidate genes concerned with the cold acclimation process in oil-tea camellia.
The leaf transcriptomes of wild oil-tea camellia obtained in our study provide one of the largest transcriptome dataset in the genus Camellia. Such a dataset can facilitate the functional genomic studies and the molecular breeding of many economically important Camellia species. Large amounts of SSRs, SNPs and InDels were identified. The phylogenetic analysis based on SNPs showed genetic differentiation between latitudes. Such a result suggests that the sequence variations identified can be used to develop molecular markers for analyzing genetic differentiations along latitude and elevation gradients in wild oil-tea camellia. Wild oil-tea camellia may undergo cold acclimation when air temperatures are below 10 °C. During cold acclimation, many genes may be up-regulated and the number of DEGs increases with the decrease in air temperature. We provide a new method for identifying significant functional groups of DEGs. Based on the new method, our study clearly showed that candidate genes for cold acclimation may be predominantly involved in transmembrane transporter activities. Our study can serve as a basis for studying molecular mechanisms of cold tolerance in evergreen broadleaf woody plants.
We are grateful to colleagues in Novogene Bioinformatics Technology Co., Ltd (Beijing, China) for supports in transcriptome sequencing and data analyses. This work was supported by the National Natural Science Foundation of China (NSFC) grant no. 31460072 and “Gan-Po Talent 555” Project of Jiangxi Province, China.
This work was supported by the National Natural Science Foundation of China (NSFC) grant no. 31460072 and “Gan-Po Talent 555” Project of Jiangxi Province, China.
Availability of data and materials
The datasets generated during the current study are available in the NCBI Short Read Archive (SRA) database with the accession number: SRR2146977 (LS01), SRR2146978 (LS02), SRR2146979 (LS03), SRR2146980 (LS04), SRR2146973 (JG01), SRR2146974 (JG02), SRR2146975 (JG03), and SRR2146976 (JG04).
JC did the sampling, conducted the experiments, analyzed the data and wrote the manuscript. XY did the sampling, conducted the experiments and analyzed the data. XH did the sampling and conducted the experiments. SD helped to conceive and design the experiments and did the sampling. CL helped to conceive and design the experiments and did the sampling. JC helped to conceive and design the experiments. JR conceived and designed the experiments, did the sampling, analyzed the data and wrote the manuscript. All authors read and approved the final manuscript.
The authors declare that they have no competing interests.
Consent for publication
Ethics approval and consent to participate
Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- Knight MR, Knight H. Low-temperature perception leading to gene expression and cold tolerance in higher plants. New Phytol. 2012;195:737–51.View ArticlePubMedGoogle Scholar
- Theocharis A, Clément C, Barka EA. Physiological and molecular changes in plants grown at low temperatures. Planta. 2012;235:1091–105.View ArticlePubMedGoogle Scholar
- Wisniewski M, Nassuth A, Teulières C, Marque C, Rowland J, Cao PB, Brown A. Genomics of cold hardiness in woody plants. Crit Rev Plant Sci. 2014;33:92–124.View ArticleGoogle Scholar
- Thomashow MF. Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999;50:571–99.View ArticlePubMedGoogle Scholar
- Ming TL. Monograph of the genus Camellia. Kunming: Yunnan Science and Technology Press; 2000.Google Scholar
- Zhuang RL. Oil-tea camellia in China. 2nd ed. Beijing: China Forestry Publishing House; 2012.Google Scholar
- Ackerman WL, Egolf D. Winter's Rose', 'Snow Flurry', and 'Polar Ice' Camellias. Hortscience. 1991;26:1432–3.Google Scholar
- Ackerman WL, Egolf DR. Winter's Charm', Winter's Hope', and Winter's Star' Camellias. Hortscience. 1992;27:855–6.Google Scholar
- Yao X, Wang Y, Wang K, Ren H. Effects of geographic latitude and longitude on fat and its fatty acid composition of oil-tea camellia seeds. China Oils and Fats. 2011;36:31–4.Google Scholar
- Linder CR. Adaptive evolution of seed oils in plants: accounting for the biogeographic distribution of saturated and unsaturated fatty acids in seed oils. Am Nat. 2000;156:442–58.View ArticleGoogle Scholar
- Ma J, Ye H, Rui Y, Chen G, Zhang N. Fatty acid composition of Camellia oleifera oil. Journal für Verbraucherschutz und Lebensmittelsicherheit. 2011;6:9–12.View ArticleGoogle Scholar
- Li H, Zhou G-y, Zhang H-y, Liu J-a. Research progress on the health function of tea oil. J Med Plant Res. 2011;5:485–9.Google Scholar
- State Forestry Adminstration of the People's Republic of China. National oil-tea camellia industry development plan (2009-2020). 2009.Google Scholar
- Hajjar R, Hodgkin T. The use of wild relatives in crop improvement: a survey of developments over the last 20 years. Euphytica. 2007;156:1–13.View ArticleGoogle Scholar
- Xia E-H, Jiang J-J, Huang H, Zhang L-P, Zhang H-B, Gao L-Z. Transcriptome analysis of the oil-rich tea plant, Camellia oleifera, reveals candidate genes related to lipid metabolism. PLoS One. 2014;9:e104150.View ArticlePubMedPubMed CentralGoogle Scholar
- Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011;29:644–52.View ArticlePubMedPubMed CentralGoogle Scholar
- Finn RD, Bateman A, Clements J, Coggill P, Eberhardt RY, Eddy SR, Heger A, Hetherington K, Holm L, Mistry J, et al. Pfam: the protein families database. Nucleic Acids Res. 2013;42:D222–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Koonin EV, Fedorova ND, Jackson JD, Jacobs AR, Krylov DM, Makarova KS, Mazumder R, Mekhedov SL, Nikolskaya AN, Rao BS, et al. A comprehensive evolutionary classification of proteins encoded in complete eukaryotic genomes. Genome Biol. 2004;5:R7.View ArticlePubMedPubMed CentralGoogle Scholar
- Boeckmann B, Bairoch A, Apweiler R, Blatter M-C, Estreicher A, Gasteiger E, Martin MJ, Michoud K, O'Donovan C, Phan I, et al. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res. 2003;31:365–70.View ArticlePubMedPubMed CentralGoogle Scholar
- Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28:27–30.View ArticlePubMedPubMed CentralGoogle Scholar
- Finn RD, Clements J, Eddy SR. HMMER web server: interactive sequence similarity searching. Nucleic Acids Res. 2011;39:W29–37.View ArticlePubMedPubMed CentralGoogle Scholar
- Götz S, García-Gómez JM, Terol J, Williams TD, Nagaraj SH, Nueda MJ, Robles M, Talón M, Dopazo J, Conesa A. High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 2008;36:3420–35.View ArticlePubMedPubMed CentralGoogle Scholar
- Moriya Y, Itoh M, Okuda S, Yoshizawa AC, Kanehisa M. KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res. 2007;35:W182–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Meth. 2012;9:357–9.View ArticleGoogle Scholar
- Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078–9.View ArticlePubMedPubMed CentralGoogle Scholar
- McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20:1297–303.View ArticlePubMedPubMed CentralGoogle Scholar
- R Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2015.Google Scholar
- Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Hohna S, Larget B, Liu L, Suchard MA, Huelsenbeck JP. MrBayes 3.2: Efficient Bayesian Phylogenetic Inference and Model Choice Across a Large Model Space. Syst Biol. 2012;61:539–42.View ArticlePubMedPubMed CentralGoogle Scholar
- Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinf. 2011;12:323.View ArticleGoogle Scholar
- Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.View ArticlePubMedPubMed CentralGoogle Scholar
- Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biol. 2010;11:R106.View ArticlePubMedPubMed CentralGoogle Scholar
- Jiang X, Liu Y, Li W, Zhao L, Meng F, Wang Y, Tan H, Yang H, Wei C, Wan X, et al. Tissue-specific, development-dependent phenolic compounds accumulation profile and gene expression pattern in tea plant [Camellia sinensis]. PLoS One. 2013;8:e62315.View ArticlePubMedPubMed CentralGoogle Scholar
- Shi C-Y, Yang H, Wei C-L, Yu O, Zhang Z-Z, Jiang C-J, Sun J, Li Y-Y, Chen Q, Xia T, et al. Deep sequencing of the Camellia sinensis transcriptome revealed candidate genes for major metabolic pathways of tea-specific compounds. BMC Genomics. 2011;12:131.View ArticlePubMedPubMed CentralGoogle Scholar
- Wang X-C, Zhao Q-Y, Ma C-L, Zhang Z-H, Cao H-L, Kong Y-M, Yue C, Hao X-Y, Chen L, Ma J-Q, et al. Global transcriptome profiles of Camellia sinensis during cold acclimation. BMC Genomics. 2013;14:415.View ArticlePubMedPubMed CentralGoogle Scholar
- Jump AS, Mátyás C, Peñuelas J. The altitude-for-latitude disparity in the range retractions of woody species. Trends Ecol Evol. 2009;24:694–701.View ArticlePubMedGoogle Scholar
- Yue C, Cao H-L, Wang L, Zhou Y-H, Huang Y-T, Hao X-Y, Wang Y-C, Wang B, Yang Y-J, Wang X-C. Effects of cold acclimation on sugar metabolism and sugar-related gene expression in tea plant during the winter season. Plant Mol Biol. 2015;88:591–608.View ArticlePubMedGoogle Scholar
- Winfield MO, Lu C, Wilson ID, Coghill JA, Edwards KJ. Plant responses to cold: transcriptome analysis of wheat. Plant Biotechnol J. 2010;8:749–71.View ArticlePubMedGoogle Scholar
- An D, Yang J, Zhang P. Transcriptome profiling of low temperature-treated cassava apical shoots showed dynamic responses of tropical plant to cold stress. BMC Genomics. 2012;13:64.View ArticlePubMedPubMed CentralGoogle Scholar
- Schulze WX, Schneider T, Starck S, Martinoia E, Trentmann O. Cold acclimation induces changes in Arabidopsis tonoplast protein abundance and activity and alters phosphorylation of tonoplast monosaccharide transporters. Plant J. 2012;69:529–41.View ArticlePubMedGoogle Scholar
- Scholz FG, Bucci SJ, Arias N, Meinzer FC, Goldstein G. Osmotic and elastic adjustments in cold desert shrubs differing in rooting depth: coping with drought and subzero temperatures. Oecologia. 2012;170:885–97.View ArticlePubMedGoogle Scholar
- Thomashow MF. Molecular basis of plant cold acclimation: insights gained from studying the CBF cold response pathway. Plant Physiol. 2010;154:571–7.View ArticlePubMedPubMed CentralGoogle Scholar