BACKGROUND
low temperatures are one of the most important environmental factors that temperate plants have to cope with during their life cycles.
some plants can enhance their freezing tolerance after exposure to low but non-freezing temperatures for a period of time, a process known as cold acclimation   <cit> .
ca is a complex process that involves cellular, physiological, metabolic and molecular modifications.
when plants sense the cold temperature, a series of protective mechanisms are triggered  <cit> .
these include resetting the cellular framework; alternating the composition, structure and function of the plasma membrane; synthesizing cryoprotectant molecules such as soluble sugars, sugar alcohols and low-molecular-weight nitrogenous compounds; decreasing the ratio of free water content; improving the scavenging activity of reactive oxygen species ; and introducing antifreeze proteins.
these alterations help plants maintain a metabolic balance of substance and energy in cold environments.
a group of cold-related genes has been reported to regulate these aforementioned changes  <cit> .
moreover, changes in gene expression have been demonstrated to occur during ca in a wide range of plant species, and hundreds of cold inducible genes have been identified  <cit> .
tea is the most popular non-alcoholic health beverage in the world, and the tea plant  o. kuntze) is one of the most important economic crops in china, india, sri lanka, kenya, among others  <cit> .
as an evergreen woody plant, the tea plant can be grown in tropical to subtropical climates.
due to the local climate changes, tea plants have to cope with low temperatures during the wintertime.
low temperatures are one of the most critical environmental factors that limit its growth, survival and geographical distribution  <cit> .
thus, finding ways to improve tea plants’ resistance to low temperatures is of great importance.
like other perennial evergreen woody crops, during the ca process, the cold tolerance of tea plants enhances with the decrease in temperature and reduces with the increase in temperature.
a previous study showed that when the average air temperature decreases to around 7°c, tea plants undergo the ca process, and after the average air temperature increases to over 9°c, tea plants start the de-acclimation process  <cit> .
there are few studies that have focused on the cellular, physiological and metabolic changes during ca in tea plants.
when tea plants undergo the ca process, the thickness of palisade tissue is increased and the stability of plasma membrane is enhanced.
in addition, the concentration of the cytochylema and ratio of bound water in the cytoplasm, the amount of unsaturated fatty acids and total proteins in the plasma membrane, and the content of soluble proteins in the leaf are also increased.
meanwhile, the activities of some detoxification enzymes, such as catalase , superoxide dismutase , peroxidase  and esterase  are increased, whereas the metabolic activity is decreased  <cit> .
some cold-induced genes have been cloned in tea plants  <cit> .
as a complex biological phenomenon, the ability of tea plants to resist the cold is regulated by a series of genes involved in a complex regulatory network  <cit> .
using an ‘omics’ research strategy to understand the mechanism of ca in tea plants is the key to improving tea productivity and geographical distribution.
rna-seq is a recently developed approach using a massively parallel sequencing strategy to generate transcriptome profiles.
it has emerged as a cost-effective approach for high-throughput sequence determination and has unprecedentedly improved the efficiency and speed of gene discovery  <cit> .
digital gene expression  is a tag-based sequencing approach according to which short tags are generated by endonuclease.
the expression level of genes in the sample is measured by counting the number of tags generated from each transcript  <cit> .
this study demonstrates the first attempt to use a combination of rna-seq and dge to study the transcriptome profiles in tea plants and thereby gain a deeper insight into the molecular mechanism of ca.
the resulting transcriptome profiles from tea plants not only contributes to the in-depth knowledge of the genes involved in ca but also improves our understanding of plant-environment interactions.
RESULTS
cold tolerance changes in tea plant during the ca process
cold tolerance in tea plants varies under different temperatures and can be monitored by the relative electrical conductivity using an electrolyte leakage assay.
figure  <dig> shows a complete course of the ca process for a natural temperature change period from december  <dig> to march  <dig>  before december  <dig>  the average outdoor temperature  was above 10°c, and the relative electrical conductivity of tea-plant leaves  was at ~100%, indicating that the tea plant has a low level of cold tolerance.
after the tea plant underwent a period of time at relatively low temperatures , its relative electrical conductivity decreased, and the cold tolerance of the tea plant is enhanced.
when temperatures decreased to their lowest point, the relative electrical conductivity also reached its lowest level with the cold tolerance being at the highest level .
afterwards, the temperature rose and when the average temperature reached above 10°c, the relative electrical conductivity increased to over 80% and then maintained at a high level.
the tea plant was subsequently de-acclimated , and its cold tolerance was weak .
to obtain the transcriptomic response to the cold environment during the ca process, we selected tea-plant leaves from three stages, non-acclimated , fully acclimated  and de-acclimated  for rna-seq and digital gene expression  studies.
rna-seq and de novo assembly
we performed rna-seq analyses for ca <dig>  ca <dig> and ck using the illumina hiseq <dig> genome analyzer.
totally,  <dig>  million paired-end reads with a read length of  <dig> bp were generated .
since low-quality nucleotides from the ends of reads may lead to incorrect assembly outputs  <cit> , we trimmed the low-quality or ambiguous nucleotides at both ends of the reads.
de novo assembly was performed with the trimmed reads using trinity  <cit> .
trinity was specially developed for de novo assembly from short-read rna-seq data, which has been shown to be the best single k-mer assembler  <cit> .
in total,  <dig>  transcripts were reconstructed.
after removing the redundant transcripts caused by small variations as described in the previous study  <cit> , a final set of  <dig>  transcripts were obtained.
the average transcript size is  <dig> bp, and the n <dig> is  <dig> bp .
the transcriptome of c. sinensis was reported in a previous study by shi et al.
<cit> .
they produced rna-seq data from the mixed tissues of c. sinensis using illumina ga iix .
a combination of dataset  <dig> and dataset  <dig> was also generated, which we called dataset  <dig>  representing all available rna-seq data for c. sinensis.
short reads of dataset  <dig> and dataset  <dig> were pre-processed by the procedure described above, and then used separately for de novo assembly.
the assembly outcome from dataset  <dig> attains the longest average read length and n <dig>  while that from dataset  <dig> yields the most number of transcripts and total base pairs .
in order to evaluate the efficiency of short-read usage during the de novo assembly, we mapped our rna-seq reads back to three sets of reconstructed transcripts, respectively.
transcripts produced from dataset  <dig> achieved the best performance, with the highest mapping ratio for our short reads .
more than 10% of the short reads failed to be aligned if only dataset  <dig> was used for the de novo assembly, indicating that previous transcriptome sequences of c. sinensis are far from saturated.
although more transcriptome sequences could be produced from de novo assembly using dataset  <dig> than dataset  <dig>  the mapping ratio could not be improved , indicating that the additional transcripts from dataset  <dig> are most likely transcripts that are expressed in tissues other than the leaves of tea plants.
thus these additional transcripts are unable to contribute to this study.
based on this scenario, we chose the transcripts from dataset  <dig> to carry out the downstream analysis.
functional annotation of c. sinensis transcriptome
to predict and analyze the function of the assembled transcripts, non-redundant sequences were submitted to a blastx  search against the following databases: the ncbi’s nr database, uniref <dig>  <cit> , the arabidopsis information resource , kyoto encyclopedia of genes and genomes   <cit>  and clusters of orthologous groups from  <dig> eukaryotic complete genomes   <cit> .
we found that about one third of all non-redundant transcripts had significant homology with genes in either the nr or uniref <dig> databases .
arabidopsis thaliana is one of the most well-studied dicot plants, with a complete reference genome and comprehensively annotated gene sequences.
a blast search against genes from arabidopsis produced more definitive annotations and helped us to evaluate the quality and coverage of our assembled transcripts.
it is notable that  <dig>  arabidopsis genes located uniformly on five chromosomes were covered by  <dig>  transcripts .
a blast analysis of the assembled transcripts against the kegg database showed that  <dig>  transcripts were annotated with corresponding enzyme commission  numbers and assigned to the reference canonical kegg pathways .
a search against the kog database reported that  <dig>  transcripts had the best hits when the e-value was less than or equal to 10- <dig>  since some transcripts could be assigned multiple kog functions, altogether  <dig>  functional annotations were produced and all hit transcripts were grouped in  <dig> categories .
in total,  <dig>  transcripts got the best hits with known proteins in at least one of the five databases and  <dig>  transcripts had similarity to proteins in all of the five databases .
to functionally categorize the assembled transcripts, gene ontology   <cit>  terms were assigned to each transcript based on the best blastx hit from the nr database using blast2go  <cit> .
out of  <dig>  transcripts with nr annotation,  <dig>  transcripts were assigned  <dig>  go term annotations in three main go categories including biological process, cellular component and molecular function .
if a gene contained some conserved domains, the domain information would be useful for interpreting the gene’s function.
to annotate the potential domains inside the reconstructed sequences, the open reading frame  was predicted for each transcript , and then all transcripts with predicted orf were used to search against the pfam database based on profile hidden markov model methods.
in total,  <dig>  transcripts were assigned pfam domain information and were categorized into  <dig>  domains/families.
most domains/families were found to contain a small number of transcripts .
according to the frequency of the occurrence of c. sinensis transcripts contained in each pfam domain, pfam domains/families were ranked and the top  <dig> abundant domains/families are listed in figure 3b, with hit results similar to the previous study  <cit> .
among these domains/families, “protein kinase domain” and its subclass “protein tyrosine kinase” are known to regulate the majority of cellular pathways.
proteins with “leucine-rich repeats” domain are known to be frequently involved in the formation of protein–protein interactions, and “ppr repeat” has been reported to be a large protein family in plants with versatile functions  <cit> .
moreover, the “nb-arc” protein family, comprised of resistance proteins, was highly represented.
other protein families, such as “reverse transcriptase” and “rna recognition motif”, which have some basic functions in plants, were also found in the top ten of the list.
trinity produced all potential alternative spliced isoforms during the de novo assembly, and isoforms originated from the same gene locus were assumed to share the same “chrysalis component”, “butterfly sub-component” and some of the paths in the de bruijn graph.
we chose the longest transcript in each locus to get the unigene set, resulting in  <dig>  unigenes .
potential isoforms in each locus reported by trinity would be useful in array/primer design for quantitative gene expression and future alternative splicing analyses.
identification of genes involved in cold acclimation
the abundance estimation for reconstructed transcripts was determined by rsem software package that was shown to have the ability to effectively use ambiguously-mapping reads and to accurately estimate isoform-level abundance for de novo assembled transcripts without reference  <cit> .
the deseq package  <cit>  and the winflat program were then applied to identify differentially expressed genes.
ca-related genes were identified based on the fold change on the abundance of each gene and the corresponding false discovery rate, which resulted in  <dig>  differentially expressed genes.
of these,  <dig>  were up-regulated and  <dig> were down-regulated , indicating that more genes were activated than repressed during the ca process.
dozens of cold-regulated or cold-related genes were found in this differential expression list, including cold sensor or signal transduction genes, cold-responsive transcription factor genes, plasma membrane stabilization related genes, osmosensing-responsive genes and detoxification enzymes genes.
cold sensor or signal transduction genes
the signal transduction pathway plays a pivotal role in the response to the stress of low temperatures  <cit> .
it is well known that ca2+ acts as a key messenger in regulating growth and developmental processes and plays a crucial role in stress signaling, i.e.
cold stress  <cit> .
cold stress could activate ca2+ channels to increase the cytosolic ca2+ level, and then trigger phospholipase c and d, producing inositol triphosphate and phosphatidic acid, respectively.
inositol triphosphate could further amplify ca2+ signatures, and phosphatidic acid is proposed as a membrane-based secondary messenger molecule  <cit> .
subsequently, many signaling pathways are triggered, such as ca2+-dependent protein kinases , mitogen-activated protein kinases , calcineurin b-like protein , calmodulin, etc.
doherty et al.
<cit>  and yang et al.
<cit>  found that the calmodulin-binding transcription activator  and a novel calcium/calmodulin-regulated receptor kinase  were crucial for cold tolerance in plants.
ca2+ influx into the cell was considered to occur upstream of the expression of cbfs and cor genes in the cold signaling pathway  <cit> .
in this study,  <dig> genes, which were annotated as cdpks, cbl, calmodulin, camta, mapk and phospholipase, were identified as being involved in signal transduction upon low temperature stress.
among these genes,  <dig>  were up-regulated in ca <dig>  whereas  <dig>  were down-regulated.
plant protein kinases belong to a large superfamily, some of which have been known to play a central role in cellular signaling, for example cdpks and mapks.
in addition, a growing body of evidence has shown that receptor-like kinases  are involved in the perception of environmental signals  <cit> .
histidine kinases , being localized to the cellular membranes and endoplasmic reticulum, are the major signaling molecules and are involved in the two-component signaling pathways that mediate plant-sensed environmental signals and regulate the downstream environmental stress response  <cit> .
in this study,  <dig> rlks genes and  <dig> hks genes were differentially expressed and all of these were up-regulated in ca <dig> samples, which indicates that protein kinases play an important role in the ca process in tea plants.
cold-responsive transcription factor genes
transcription factors  play important functions in plant development and stress tolerance  <cit> .
fifty-eight genes encoding putative tfs in c. sinensis were identified.
these tfs could be divided into  <dig> groups  based on the classification of their arabidopsis homologs, and most of them have been reported to be linked to cold stress resistance in plants  <cit> .
among these tfs,  <dig> genes were up-regulated and  <dig> genes were down-regulated in our ca <dig> sample.
of the  <dig> groups of tfs, zinc finger was the most enriched tf family, containing  <dig> genes of the  <dig> cold-responsive tfs, with  <dig> genes being up-regulated and  <dig> being down-regulated.
there were  <dig> genes in the bhlh family ,  <dig> genes in the myb family ,  <dig> genes in the wrky family  and  <dig> genes in the nac family .
in addition,  <dig> genes in the bzip family,  <dig> genes in the gars family and  <dig> genes encoding heat shock proteins were all up-regulated, while  <dig> genes in the ap2/erf family were down-regulated in the ca <dig> sample.
it is interesting to find down-regulated genes in the ap2/erf family, as these suggest that the interaction of light and temperature is of special importance for plants during the ca process.
catalá et al.
<cit>  and jurczyk et al.
<cit>  have also reported that light is required for full ca in arabidopsis and festuca pratensis.
genes related to the stabilization of the plasma membrane and osmosensing-responsiveness
the plasma membrane is believed to be a primary site of injury from freezing in plants.
the process of ca can stabilize the membrane structure and prevent it from damage  <cit> .
under freezing temperatures, membranes must be kept fluid in order to sustain the functional activity of membrane proteins and membranes themselves  <cit> .
alterations occur in the composition of proteins and lipids  in the plasma membrane in response to ca, and these are associated with an increase in freezing tolerance  <cit> .
in our study, we identified  <dig> lipid-transfer protein  genes and  <dig> fatty acid desaturase  gene.
among these,  <dig> ltp genes and a fad gene were up-regulated and  <dig> ltp gene was down-regulated.
these genes were known to regulate the level of unsaturated fatty acids, and then to further mediate the regulation of membrane fluidity  <cit> .
moreover, in order to maintain the structural stabilization of the plasma membrane during the ca process, some proteins function as inhibitors to regulate the activity of ice nucleators.
these proteins are so-called anti-freezing proteins , such as β- <dig>  3-glucanase-like proteins , chitinase-like proteins , thaumatin-like proteins , polygalacturonase inhibitor proteins  and late-embryogenesis-abundant proteins .
in the ca <dig> sample, more genes encoding these proteins were up-regulated compared with genes in non-acclimated samples  <cit> .
in our study, we found  <dig> afp-related genes, including  <dig> clps,  <dig> tlp,  <dig> pgip and  <dig> lea that were up-regulated in the ca <dig> sample, indicating that during the ca process, tea plants became able to tolerate freezing temperatures through the enhancement of membrane stability.
the stabilization of the plasma membranes is also related to the osmotic equilibrium.
in order to maintain osmotic balance, plants accumulate a range of compatible solutes, including soluble sugars , sugar alcohols , and low-molecular-weight compounds  as cryoprotectant molecules in response to cold stress  <cit> .
accordingly, the expression of these metabolism related genes also changes during ca  <cit> .
we identified  <dig> genes related to the carbohydrate metabolic pathway from  <dig>  differentially expressed genes, including  <dig> galactosidases ,  <dig> amylases ,  <dig> galactinol synthase ,  <dig> raffinose synthase  and  <dig> trehalose-6-phosphate synthases .
these genes are key genes of the carbohydrate metabolic pathway, and are closely involved with the ca process  <cit> .
three monosaccharide transporter genes  were identified as well.
monosaccharide transporters play an important role in sugar transport and distribution in plants.
the expression of monosaccharide transporter genes is also regulated by cold stress  <cit> .
these results suggested that the carbohydrate metabolic pathway plays a critical role in tea plants during the ca process.
validation of rna-seq results by dge and qrt-pcr
digital gene expression  library sequencing was performed to validate the cold-regulated transcripts identified by rna-seq.
in our study, three dge libraries were sequenced: ca <dig>  ca <dig> and ck, for which  <dig> ,  <dig>  and  <dig>  million raw tags were generated, respectively .
after removing low-quality tags, the total number of clean tags per library ranged from  <dig>  to  <dig>  million.
clean tags from three dge libraries were mapped onto our assembled transcriptome sequences.
up to  <dig> % of transcripts were detected by dge tags .
of the  <dig>  differentially expressed transcripts identified by rna-seq,  <dig>  were detected by dge sequencing, but  <dig> were mapped by uncertain tags  and another  <dig> transcripts did not have enough tags  counts for all three samples) to differentiate expressions among ca <dig>  ca <dig> and ck samples.
this result illustrates that dge sequencing was limited to identify differential expression across the full scale of transcriptome profiles, especially for genes with paralogs or multiple isoforms that shared the same tags .
of the remaining  <dig> transcripts, the majority of them  showed consistent expression patterns between dge and rna-seq, with the corresponding pearson’s r being  <dig>  and  <dig>  for ca1/ca <dig> and ca1/ck, respectively , demonstrating the degree of consistency between dge and rna-seq platforms.
it is worth noting that some transcripts, though not many, showed different expression patterns in the profiling results from rna-seq and dge .
determining which method is more robust and why the two approaches yield different results would be useful for identifying the correct outcomes in this study and for other researchers to choose the appropriate approach in their future studies.
to address this,  <dig> of these transcripts that showed inconsistent results from rna-seq and dge platforms were randomly selected to assess their relative expression patterns among ck, ca <dig> and ca <dig> using quantitative rt-pcr approach .
for most of these , similar expression patterns were observed compared with those from rna-seq results, while in the other  <dig> transcripts there were only partial consistencies with either rna-seq or dge results .
in general, rna-seq outperforms dge based on the results from these  <dig> cases.
the less accurate estimation of the gene expression level by dge approach could be due to some unknown reason or to the fact that the same tags may exist in other transcripts that were partially reconstructed after de novo transcriptome assembly and lack the complete tag sequences.
since the dge approach counts all tags to the transcript with the exactly matched tag sequences, this may result in the incorrect estimation of the expression level for some transcripts.
in the remaining two genes, inconsistent expression patterns were observed among the results from the three approaches.
these genes expressed at relatively variable levels may be affected by factors other than a cold environment and these kinds of false positives could be largely avoided if more biological replicates were included.
the dge method is widely applied for studying the transcriptome.
however, it has limitations in presenting a global view of transcriptome profiles.
it is powerless to detect the abundance of transcripts when 1) there is no catg site in transcripts or 2) several transcripts share the same tag, situations involving two unrelated genes, paralogs, or alternatively spliced isoforms.
both closely related paralogs and alternatively spliced isoforms might exhibit various spatial and temporal expression patterns, or even have different functions.
thus, the ability to correctly estimate isoform expression levels will be necessary for understanding complicated biological mechanisms.
to further test the reliability of outcomes produced from the next generation sequencing platform, quantitative rt-pcr  analysis was performed for  <dig> of the  <dig>  differentially expressed transcripts.
these  <dig> transcripts were manually selected as representatives for their potential roles in cold tolerance according to their annotations.
the expression patterns of  <dig> genes detected by qrt-pcr fit well with those from rna-seq results, with one annotated as a basic helix-loop-helix dna-binding superfamily protein being inconsistent .
pathways involved during the ca process in c. sinensis
the  <dig>  transcripts were used to search the kegg pathway to determine whether the genes involved in ca were from specific pathways.
in total,  <dig> pathways were identified,  <dig> of which were significantly enriched during ca .
of these significantly enriched, metabolism was the largest category , including carbohydrate metabolism , glycan biosynthesis and metabolism , energy metabolism , amino acid metabolism , metabolism of terpenoids and polyketides , enzyme families , xenobiotics biodegradation and metabolism  and lipid metabolism .
moreover, calcium signaling pathway  and membrane transport pathway  were enriched as well.
numerous studies reported that carbohydrate metabolism plays an important role during the ca process  <cit> .
in this study, metabolic pathways for carbohydrates stood out from the enrichment analysis, including pathways for  <dig> differentially expressed transcripts , indicating that the regulation of carbohydrates is crucial for tea plants during ca.
previous studies have shown that calcium acts as a pivotal mediator in the signal transduction pathway during the ca process  <cit> .
calcium is a secondary messenger in plant signaling processes, and calcium/calmodulin-mediated signaling is believed to play an important role in plants during the cold stress response  <cit> .
in this study, the calcium signaling pathway was also enriched, and most of the genes in this pathway were up-regulated in the ca <dig> sample, proving the importance of the calcium signaling pathway for the tea plants’ response to cold stress.
CONCLUSIONS
in this study, we present a global survey for transcriptome profiles in tea plants during the ca process using rna-seq and dge.
a large number of genes from tea plants involved in diverse biological or molecular pathways were identified during the ca process, such as genes involved in cold signal sensors or transduction, genes related to the stabilization of plasma membranes, osmosensing-responsive genes, and stress-responsive transcription factor genes.
a diagram is shown to illustrate tea plants’ responses to low temperatures during the ca process .
the results showed that a series of complex regulatory networks were triggered in tea plants during ca.
our study provides insights into the molecular mechanisms of tea plants during the ca process.
it could also serve as a valuable resource for relevant research on cold-tolerance and help to explore the cold-related genes in improving the understanding of low-temperature tolerance and plant-environment interactions.

 low temperatures are one of the most important environmental factors that temperate plants have to cope with during their life cycles. some plants can enhance their freezing tolerance after exposure to low but non-freezing temperatures for a period of time, a process known as cold acclimation. ca is a complex process that involves cellular, physiological, metabolic and molecular modifications. when plants sense the cold temperature, a series of protective mechanisms are triggered. these include resetting the cellular framework; alternating the composition, structure and function of the plasma membrane; synthesizing cryoprotectant molecules such as soluble sugars, sugar alcohols and low-molecular-weight nitrogenous compounds; decreasing the ratio of free water content; improving the scavenging activity of reactive oxygen species; and introducing antifreeze proteins. these alterations help plants maintain a metabolic balance of substance and energy in cold environments. a group of cold-related genes has been reported to regulate these aforementioned changes. moreover, changes in gene expression have been demonstrated to occur during ca in a wide range of plant species, and hundreds of cold inducible genes have been identified. tea is the most popular non-alcoholic health beverage in the world, and the tea plant o. kuntze) is one of the most important economic crops in china, india, sri lanka, kenya, among others. as an evergreen woody plant, the tea plant can be grown in tropical to subtropical climates. due to the local climate changes, tea plants have to cope with low temperatures during the wintertime. low temperatures are one of the most critical environmental factors that limit its growth, survival and geographical distribution. thus, finding ways to improve tea plants’ resistance to low temperatures is of great importance. like other perennial evergreen woody crops, during the ca process, the cold tolerance of tea plants enhances with the decrease in temperature and reduces with the increase in temperature. a previous study showed that when the average air temperature decreases to around 7°c, tea plants undergo the ca process, and after the average air temperature increases to over 9°c, tea plants start the de-acclimation process. there are few studies that have focused on the cellular, physiological and metabolic changes during ca in tea plants. when tea plants undergo the ca process, the thickness of palisade tissue is increased and the stability of plasma membrane is enhanced. in addition, the concentration of the cytochylema and ratio of bound water in the cytoplasm, the amount of unsaturated fatty acids and total proteins in the plasma membrane, and the content of soluble proteins in the leaf are also increased. meanwhile, the activities of some detoxification enzymes, such as catalase, superoxide dismutase, peroxidase and esterase are increased, whereas the metabolic activity is decreased. some cold-induced genes have been cloned in tea plants. as a complex biological phenomenon, the ability of tea plants to resist the cold is regulated by a series of genes involved in a complex regulatory network. using an ‘omics’ research strategy to understand the mechanism of ca in tea plants is the key to improving tea productivity and geographical distribution. rna-seq is a recently developed approach using a massively parallel sequencing strategy to generate transcriptome profiles. it has emerged as a cost-effective approach for high-throughput sequence determination and has unprecedentedly improved the efficiency and speed of gene discovery. digital gene expression is a tag-based sequencing approach according to which short tags are generated by endonuclease. the expression level of genes in the sample is measured by counting the number of tags generated from each transcript. this study demonstrates the first attempt to use a combination of rna-seq and dge to study the transcriptome profiles in tea plants and thereby gain a deeper insight into the molecular mechanism of ca. the resulting transcriptome profiles from tea plants not only contributes to the in-depth knowledge of the genes involved in ca but also improves our understanding of plant-environment interactions.  cold tolerance changes in tea plant during the ca process cold tolerance in tea plants varies under different temperatures and can be monitored by the relative electrical conductivity using an electrolyte leakage assay. figure  shows a complete course of the ca process for a natural temperature change period from december to march before december the average outdoor temperature was above 10°c, and the relative electrical conductivity of tea-plant leaves was at ~100%, indicating that the tea plant has a low level of cold tolerance. after the tea plant underwent a period of time at relatively low temperatures, its relative electrical conductivity decreased, and the cold tolerance of the tea plant is enhanced. when temperatures decreased to their lowest point, the relative electrical conductivity also reached its lowest level with the cold tolerance being at the highest level. afterwards, the temperature rose and when the average temperature reached above 10°c, the relative electrical conductivity increased to over 80% and then maintained at a high level. the tea plant was subsequently de-acclimated, and its cold tolerance was weak. to obtain the transcriptomic response to the cold environment during the ca process, we selected tea-plant leaves from three stages, non-acclimated, fully acclimated and de-acclimated for rna-seq and digital gene expression studies. rna-seq and de novo assembly we performed rna-seq analyses for ca ca and ck using the illumina hiseq genome analyzer. totally, million paired-end reads with a read length of bp were generated. since low-quality nucleotides from the ends of reads may lead to incorrect assembly outputs, we trimmed the low-quality or ambiguous nucleotides at both ends of the reads. de novo assembly was performed with the trimmed reads using trinity. trinity was specially developed for de novo assembly from short-read rna-seq data, which has been shown to be the best single k-mer assembler. in total, transcripts were reconstructed. after removing the redundant transcripts caused by small variations as described in the previous study, a final set of transcripts were obtained. the average transcript size is bp, and the n is bp. the transcriptome of c. sinensis was reported in a previous study by shi et al.. they produced rna-seq data from the mixed tissues of c. sinensis using illumina ga iix. a combination of dataset and dataset was also generated, which we called dataset representing all available rna-seq data for c. sinensis. short reads of dataset and dataset were pre-processed by the procedure described above, and then used separately for de novo assembly. the assembly outcome from dataset attains the longest average read length and n while that from dataset yields the most number of transcripts and total base pairs. in order to evaluate the efficiency of short-read usage during the de novo assembly, we mapped our rna-seq reads back to three sets of reconstructed transcripts, respectively. transcripts produced from dataset achieved the best performance, with the highest mapping ratio for our short reads. more than 10% of the short reads failed to be aligned if only dataset was used for the de novo assembly, indicating that previous transcriptome sequences of c. sinensis are far from saturated. although more transcriptome sequences could be produced from de novo assembly using dataset than dataset the mapping ratio could not be improved, indicating that the additional transcripts from dataset are most likely transcripts that are expressed in tissues other than the leaves of tea plants. thus these additional transcripts are unable to contribute to this study. based on this scenario, we chose the transcripts from dataset to carry out the downstream analysis. functional annotation of c. sinensis transcriptome to predict and analyze the function of the assembled transcripts, non-redundant sequences were submitted to a blastx search against the following databases: the ncbi’s nr database, uniref, the arabidopsis information resource, kyoto encyclopedia of genes and genomes and clusters of orthologous groups from eukaryotic complete genomes. we found that about one third of all non-redundant transcripts had significant homology with genes in either the nr or uniref databases. arabidopsis thaliana is one of the most well-studied dicot plants, with a complete reference genome and comprehensively annotated gene sequences. a blast search against genes from arabidopsis produced more definitive annotations and helped us to evaluate the quality and coverage of our assembled transcripts. it is notable that arabidopsis genes located uniformly on five chromosomes were covered by transcripts. a blast analysis of the assembled transcripts against the kegg database showed that transcripts were annotated with corresponding enzyme commission numbers and assigned to the reference canonical kegg pathways. a search against the kog database reported that transcripts had the best hits when the e-value was less than or equal to 10- since some transcripts could be assigned multiple kog functions, altogether functional annotations were produced and all hit transcripts were grouped in categories. in total, transcripts got the best hits with known proteins in at least one of the five databases and transcripts had similarity to proteins in all of the five databases. to functionally categorize the assembled transcripts, gene ontology terms were assigned to each transcript based on the best blastx hit from the nr database using blast2go. out of transcripts with nr annotation, transcripts were assigned go term annotations in three main go categories including biological process, cellular component and molecular function. if a gene contained some conserved domains, the domain information would be useful for interpreting the gene’s function. to annotate the potential domains inside the reconstructed sequences, the open reading frame was predicted for each transcript, and then all transcripts with predicted orf were used to search against the pfam database based on profile hidden markov model methods. in total, transcripts were assigned pfam domain information and were categorized into domains/families. most domains/families were found to contain a small number of transcripts. according to the frequency of the occurrence of c. sinensis transcripts contained in each pfam domain, pfam domains/families were ranked and the top abundant domains/families are listed in figure 3b, with hit results similar to the previous study. among these domains/families, “protein kinase domain” and its subclass “protein tyrosine kinase” are known to regulate the majority of cellular pathways. proteins with “leucine-rich repeats” domain are known to be frequently involved in the formation of protein–protein interactions, and “ppr repeat” has been reported to be a large protein family in plants with versatile functions. moreover, the “nb-arc” protein family, comprised of resistance proteins, was highly represented. other protein families, such as “reverse transcriptase” and “rna recognition motif”, which have some basic functions in plants, were also found in the top ten of the list. trinity produced all potential alternative spliced isoforms during the de novo assembly, and isoforms originated from the same gene locus were assumed to share the same “chrysalis component”, “butterfly sub-component” and some of the paths in the de bruijn graph. we chose the longest transcript in each locus to get the unigene set, resulting in unigenes. potential isoforms in each locus reported by trinity would be useful in array/primer design for quantitative gene expression and future alternative splicing analyses. identification of genes involved in cold acclimation the abundance estimation for reconstructed transcripts was determined by rsem software package that was shown to have the ability to effectively use ambiguously-mapping reads and to accurately estimate isoform-level abundance for de novo assembled transcripts without reference. the deseq package and the winflat program were then applied to identify differentially expressed genes. ca-related genes were identified based on the fold change on the abundance of each gene and the corresponding false discovery rate, which resulted in differentially expressed genes. of these, were up-regulated and were down-regulated, indicating that more genes were activated than repressed during the ca process. dozens of cold-regulated or cold-related genes were found in this differential expression list, including cold sensor or signal transduction genes, cold-responsive transcription factor genes, plasma membrane stabilization related genes, osmosensing-responsive genes and detoxification enzymes genes. cold sensor or signal transduction genes the signal transduction pathway plays a pivotal role in the response to the stress of low temperatures. it is well known that ca2+ acts as a key messenger in regulating growth and developmental processes and plays a crucial role in stress signaling, i.e. cold stress. cold stress could activate ca2+ channels to increase the cytosolic ca2+ level, and then trigger phospholipase c and d, producing inositol triphosphate and phosphatidic acid, respectively. inositol triphosphate could further amplify ca2+ signatures, and phosphatidic acid is proposed as a membrane-based secondary messenger molecule. subsequently, many signaling pathways are triggered, such as ca2+-dependent protein kinases, mitogen-activated protein kinases, calcineurin b-like protein, calmodulin, etc. doherty et al. and yang et al. found that the calmodulin-binding transcription activator and a novel calcium/calmodulin-regulated receptor kinase were crucial for cold tolerance in plants. ca2+ influx into the cell was considered to occur upstream of the expression of cbfs and cor genes in the cold signaling pathway. in this study, genes, which were annotated as cdpks, cbl, calmodulin, camta, mapk and phospholipase, were identified as being involved in signal transduction upon low temperature stress. among these genes, were up-regulated in ca whereas were down-regulated. plant protein kinases belong to a large superfamily, some of which have been known to play a central role in cellular signaling, for example cdpks and mapks. in addition, a growing body of evidence has shown that receptor-like kinases are involved in the perception of environmental signals. histidine kinases, being localized to the cellular membranes and endoplasmic reticulum, are the major signaling molecules and are involved in the two-component signaling pathways that mediate plant-sensed environmental signals and regulate the downstream environmental stress response. in this study, rlks genes and hks genes were differentially expressed and all of these were up-regulated in ca samples, which indicates that protein kinases play an important role in the ca process in tea plants. cold-responsive transcription factor genes transcription factors play important functions in plant development and stress tolerance. fifty-eight genes encoding putative tfs in c. sinensis were identified. these tfs could be divided into groups based on the classification of their arabidopsis homologs, and most of them have been reported to be linked to cold stress resistance in plants. among these tfs, genes were up-regulated and genes were down-regulated in our ca sample. of the groups of tfs, zinc finger was the most enriched tf family, containing genes of the cold-responsive tfs, with genes being up-regulated and being down-regulated. there were genes in the bhlh family, genes in the myb family, genes in the wrky family and genes in the nac family. in addition, genes in the bzip family, genes in the gars family and genes encoding heat shock proteins were all up-regulated, while genes in the ap2/erf family were down-regulated in the ca sample. it is interesting to find down-regulated genes in the ap2/erf family, as these suggest that the interaction of light and temperature is of special importance for plants during the ca process. catalá et al. and jurczyk et al. have also reported that light is required for full ca in arabidopsis and festuca pratensis. genes related to the stabilization of the plasma membrane and osmosensing-responsiveness the plasma membrane is believed to be a primary site of injury from freezing in plants. the process of ca can stabilize the membrane structure and prevent it from damage. under freezing temperatures, membranes must be kept fluid in order to sustain the functional activity of membrane proteins and membranes themselves. alterations occur in the composition of proteins and lipids in the plasma membrane in response to ca, and these are associated with an increase in freezing tolerance. in our study, we identified lipid-transfer protein genes and fatty acid desaturase gene. among these, ltp genes and a fad gene were up-regulated and ltp gene was down-regulated. these genes were known to regulate the level of unsaturated fatty acids, and then to further mediate the regulation of membrane fluidity. moreover, in order to maintain the structural stabilization of the plasma membrane during the ca process, some proteins function as inhibitors to regulate the activity of ice nucleators. these proteins are so-called anti-freezing proteins, such as β- 3-glucanase-like proteins, chitinase-like proteins, thaumatin-like proteins, polygalacturonase inhibitor proteins and late-embryogenesis-abundant proteins. in the ca sample, more genes encoding these proteins were up-regulated compared with genes in non-acclimated samples. in our study, we found afp-related genes, including clps, tlp, pgip and lea that were up-regulated in the ca sample, indicating that during the ca process, tea plants became able to tolerate freezing temperatures through the enhancement of membrane stability. the stabilization of the plasma membranes is also related to the osmotic equilibrium. in order to maintain osmotic balance, plants accumulate a range of compatible solutes, including soluble sugars, sugar alcohols, and low-molecular-weight compounds as cryoprotectant molecules in response to cold stress. accordingly, the expression of these metabolism related genes also changes during ca. we identified genes related to the carbohydrate metabolic pathway from differentially expressed genes, including galactosidases, amylases, galactinol synthase, raffinose synthase and trehalose-6-phosphate synthases. these genes are key genes of the carbohydrate metabolic pathway, and are closely involved with the ca process. three monosaccharide transporter genes were identified as well. monosaccharide transporters play an important role in sugar transport and distribution in plants. the expression of monosaccharide transporter genes is also regulated by cold stress. these results suggested that the carbohydrate metabolic pathway plays a critical role in tea plants during the ca process. validation of rna-seq results by dge and qrt-pcr digital gene expression library sequencing was performed to validate the cold-regulated transcripts identified by rna-seq. in our study, three dge libraries were sequenced: ca ca and ck, for which, and million raw tags were generated, respectively. after removing low-quality tags, the total number of clean tags per library ranged from to million. clean tags from three dge libraries were mapped onto our assembled transcriptome sequences. up to % of transcripts were detected by dge tags. of the differentially expressed transcripts identified by rna-seq, were detected by dge sequencing, but were mapped by uncertain tags and another transcripts did not have enough tags counts for all three samples) to differentiate expressions among ca ca and ck samples. this result illustrates that dge sequencing was limited to identify differential expression across the full scale of transcriptome profiles, especially for genes with paralogs or multiple isoforms that shared the same tags. of the remaining transcripts, the majority of them showed consistent expression patterns between dge and rna-seq, with the corresponding pearson’s r being and for ca1/ca and ca1/ck, respectively, demonstrating the degree of consistency between dge and rna-seq platforms. it is worth noting that some transcripts, though not many, showed different expression patterns in the profiling results from rna-seq and dge. determining which method is more robust and why the two approaches yield different results would be useful for identifying the correct outcomes in this study and for other researchers to choose the appropriate approach in their future studies. to address this, of these transcripts that showed inconsistent results from rna-seq and dge platforms were randomly selected to assess their relative expression patterns among ck, ca and ca using quantitative rt-pcr approach. for most of these, similar expression patterns were observed compared with those from rna-seq results, while in the other transcripts there were only partial consistencies with either rna-seq or dge results. in general, rna-seq outperforms dge based on the results from these cases. the less accurate estimation of the gene expression level by dge approach could be due to some unknown reason or to the fact that the same tags may exist in other transcripts that were partially reconstructed after de novo transcriptome assembly and lack the complete tag sequences. since the dge approach counts all tags to the transcript with the exactly matched tag sequences, this may result in the incorrect estimation of the expression level for some transcripts. in the remaining two genes, inconsistent expression patterns were observed among the results from the three approaches. these genes expressed at relatively variable levels may be affected by factors other than a cold environment and these kinds of false positives could be largely avoided if more biological replicates were included. the dge method is widely applied for studying the transcriptome. however, it has limitations in presenting a global view of transcriptome profiles. it is powerless to detect the abundance of transcripts when 1) there is no catg site in transcripts or 2) several transcripts share the same tag, situations involving two unrelated genes, paralogs, or alternatively spliced isoforms. both closely related paralogs and alternatively spliced isoforms might exhibit various spatial and temporal expression patterns, or even have different functions. thus, the ability to correctly estimate isoform expression levels will be necessary for understanding complicated biological mechanisms. to further test the reliability of outcomes produced from the next generation sequencing platform, quantitative rt-pcr analysis was performed for of the differentially expressed transcripts. these transcripts were manually selected as representatives for their potential roles in cold tolerance according to their annotations. the expression patterns of genes detected by qrt-pcr fit well with those from rna-seq results, with one annotated as a basic helix-loop-helix dna-binding superfamily protein being inconsistent. pathways involved during the ca process in c. sinensis the transcripts were used to search the kegg pathway to determine whether the genes involved in ca were from specific pathways. in total, pathways were identified, of which were significantly enriched during ca. of these significantly enriched, metabolism was the largest category, including carbohydrate metabolism, glycan biosynthesis and metabolism, energy metabolism, amino acid metabolism, metabolism of terpenoids and polyketides, enzyme families, xenobiotics biodegradation and metabolism and lipid metabolism. moreover, calcium signaling pathway and membrane transport pathway were enriched as well. numerous studies reported that carbohydrate metabolism plays an important role during the ca process. in this study, metabolic pathways for carbohydrates stood out from the enrichment analysis, including pathways for differentially expressed transcripts, indicating that the regulation of carbohydrates is crucial for tea plants during ca. previous studies have shown that calcium acts as a pivotal mediator in the signal transduction pathway during the ca process. calcium is a secondary messenger in plant signaling processes, and calcium/calmodulin-mediated signaling is believed to play an important role in plants during the cold stress response. in this study, the calcium signaling pathway was also enriched, and most of the genes in this pathway were up-regulated in the ca sample, proving the importance of the calcium signaling pathway for the tea plants’ response to cold stress.  in this study, we present a global survey for transcriptome profiles in tea plants during the ca process using rna-seq and dge. a large number of genes from tea plants involved in diverse biological or molecular pathways were identified during the ca process, such as genes involved in cold signal sensors or transduction, genes related to the stabilization of plasma membranes, osmosensing-responsive genes, and stress-responsive transcription factor genes. a diagram is shown to illustrate tea plants’ responses to low temperatures during the ca process. the results showed that a series of complex regulatory networks were triggered in tea plants during ca. our study provides insights into the molecular mechanisms of tea plants during the ca process. it could also serve as a valuable resource for relevant research on cold-tolerance and help to explore the cold-related genes in improving the understanding of low-temperature tolerance and plant-environment interactions. 
it could also serve as a valuable resource for relevant research on cold-tolerance and help to explore the cold-related genes in improving the understanding of low-temperature tolerance and plant-environment interactions. dge and quantitative rt-pcr analysis further confirmed the results from rna-seq analysis. these include a group of cold sensor or signal transduction genes, cold-responsive transcription factor genes, plasma membrane stabilization related genes, osmosensing-responsive genes, and detoxification enzyme genes. to elucidate the molecular mechanisms of cold acclimation, we employed rna-seq and digital gene expression technologies to the study of genome-wide expression profiles during cold acclimation in tea plants. our study presents a global survey of transcriptome profiles of tea plants in response to low, non-freezing temperatures and yields insights into the molecular mechanisms of tea plants during the cold acclimation process. the tea plant o. kuntze) needs to undergo a cold acclimation process to enhance its freezing tolerance in winter. camellia sinensiscold acclimationrna-seqdgegenome-wide expression profilestea plants pathway analysis indicated that the “carbohydrate metabolism pathway” and the “calcium signaling pathway” might play a vital role in tea plants’ responses to cold stress. changes that occur at the molecular level in response to low temperatures are poorly understood in tea plants. tea is the most popular non-alcoholic health beverage in the world.
 tea is the most popular non-alcoholic health beverage in the world. the tea plant o. kuntze) needs to undergo a cold acclimation process to enhance its freezing tolerance in winter. changes that occur at the molecular level in response to low temperatures are poorly understood in tea plants. to elucidate the molecular mechanisms of cold acclimation, we employed rna-seq and digital gene expression technologies to the study of genome-wide expression profiles during cold acclimation in tea plants.  using the illumina sequencing platform, we obtained approximately million rna-seq reads. these reads were assembled into transcripts, with an average length of bp and an n of bp. in total, differentially expressed transcripts were identified, of which were up-regulated and down-regulated. these include a group of cold sensor or signal transduction genes, cold-responsive transcription factor genes, plasma membrane stabilization related genes, osmosensing-responsive genes, and detoxification enzyme genes. dge and quantitative rt-pcr analysis further confirmed the results from rna-seq analysis. pathway analysis indicated that the “carbohydrate metabolism pathway” and the “calcium signaling pathway” might play a vital role in tea plants’ responses to cold stress.  our study presents a global survey of transcriptome profiles of tea plants in response to low, non-freezing temperatures and yields insights into the molecular mechanisms of tea plants during the cold acclimation process. it could also serve as a valuable resource for relevant research on cold-tolerance and help to explore the cold-related genes in improving the understanding of low-temperature tolerance and plant-environment interactions. camellia sinensiscold acclimationrna-seqdgegenome-wide expression profilestea plants 