Deep-sequencing transcriptome analysis of low temperature perception in a desert tree, Populus euphratica
© Chen et al.; licensee BioMed Central Ltd. 2014
Received: 17 July 2013
Accepted: 23 April 2014
Published: 1 May 2014
Compared with other Populus species, Populus euphratica Oliv. exhibits better tolerance to abiotic stress, especially those involving extreme temperatures. However, little is known about gene regulation and signaling pathways involved in low temperature stress responses in this species. Recent development of Illumina/Solexa-based deep-sequencing technologies has accelerated the study of global transcription profiling under specific conditions. To understand the gene network controlling low temperature perception in P. euphratica, we performed transcriptome sequencing using Solexa sequence analysis to generate a leaf transcriptome at a depth of 10 gigabases for each sample.
Using the Trinity method, 52,081,238 high-quality trimmed reads were assembled into a non-redundant set and 108,502 unigenes with an average length of 1,047 bp were generated. After performing functional annotations by aligning all-unigenes with public protein databases, 85,584 unigenes were annotated. Differentially expressed genes were investigated using the FPKM method by applying the Benjamini and Hochberg corrections. Overall, 2,858 transcripts were identified as differentially expressed unigenes in at least two samples and 131 were assigned as unigenes expressed differently in all three samples. In 4°C-treated sample and -4°C-treated sample, 1,661 and 866 differently expressed unigenes were detected at an estimated absolute log2-fold change of > 1, respectively. Among them, the respective number of up-regulated unigenes in C4 and F4 sample was 1,113 and 630, while the respective number of down-regulated ungenes is 548 and 236. To increase our understanding of these differentially expressed genes, we performed gene ontology enrichment and metabolic pathway enrichment analyses. A large number of early cold (below or above freezing temperature)-responsive genes were identified, suggesting that a multitude of transcriptional cascades function in cold perception. Analyses of multiple cold-responsive genes, transcription factors, and some key transduction components involved in ABA and calcium signaling revealed their potential function in low temperature responses in P. euphratica.
Our results provide a global transcriptome picture of P. euphratica under low temperature stress. The potential cold stress related transcripts identified in this study provide valuable information for further understanding the molecular mechanisms of low temperature perception in P. euphratica.
KeywordsIllumina/Solexa Low temperature Populus euphratica Transcriptome
Low temperature has strong influence on the geographical distribution, growing season, quality and yield of plants. Previous reports have shown that plants may develop acquired freezing tolerance after exposure to non-freezing low temperatures [1, 2]. Plant cells cope with cold stress by regulating the expression of transcription factors and effectors during non-freezing low temperatures . However, the transcriptome-level changes that underlie perception of temperatures below zero, which may be related to the ability to survive under extremely low temperatures , is poorly understood comparing that with cold sensing above freezing temperature.
A variety of cold-responsive (COR) genes that are under the control of some key transcription factors (TFs) are thought to be involved in cold sensing . For example, the well characterized TF DREB1/CBF can regulate a subset of COR genes by binding the DRE/CRT cis-elements in promoter regions of COR genes [6–8]. By studying the DREB1/CBF network pathway, the roles of cellular or environmental factors, e.g. calcium, light, and circadian rhythm, are revealed in cold acclimation [1, 9]. The DREB1/CBF pathway in chilling response seemed well characterized in some plants, and its regulon has been identified in Brassica napus, rice, and poplar. However, the molecular mechanisms of cold-acclimation response are not well understood at the whole genome or transcriptome level since only 12% of cold responsive genes are members of the CBF regulon [10, 11]. It remains to be answered whether low temperature perception occurs below freezing temperature, and if so, whether it occurs by a similar molecular mechanism as above freezing temperature.
Populus euphratica Oliv. is naturally distributed in semiarid areas and plays an important role in maintaining local arid ecosystems . They distinguish themselves considerably from other species by growing in deserts with extremely salty soil and wide environmental temperature ranges. Thus, P. euphratica has been widely considered as a model species for elucidating abiotic resistance mechanisms in trees [13–17]. Screening for cold responsive genes in P. euphratica can be a useful approach to understand the responses of woody plants to low temperatures and can also help elucidate the difference in cold perception between below- and above- freezing temperatures.
Recently, the development of Illumina/Solexa-based deep-sequencing technologies has made it possible to capture an unbiased view of the RNA transcript profile of a species under a given condition at the whole genome level . Using this method, ESTs and numerous novel transcripts have been discovered in a tissue-specific manner [19, 20]. In the current study, we sought to identify genes linked to low temperature (below or above zero) perception and to explore the regulatory and sensory mechanisms involved in low temperature response processes by performing de novo assembly of the P. euphratica transcriptome using Solexa data. Two-year-old plants were subjected to temperatures of 4°C and a further drop to -4°C to conduct comprehensive analysis of transcriptional responses. The acquired information may facilitate attempts to elucidate response mechanisms of this species to chilling stress and will help in the development of strategies for improving of freezing tolerance in trees.
Results and discussion
Reads assembly and poplar databses alignment
Overview of the sequencing and assembly
Total Raw Reads
Total Clean Reads
Total Clean Nucleotides (nt)
Total Number of contig
Total Length of contig (nt)
Mean Length of contig (nt)
N50 of contig
Total Number of unigne
Total Length of unigne (nt)
Mean Length of unigne (nt)
N50 of unigene
To estimate the representation of unigenes in the collection, all unique sequences generated from different assemblages were subjected to a BLAST comparison to compare EST collections from a variety of Populus species. The results indicated that our assemblages covered most P. euphratica transcripts (Additional file 2). By performing BLASTx against the Populus trichocarpa v3 dataset with an E-value of 1.0E-5 as threshold, 83,618 ESTs were assigned with an identity score ≥ 75%, covering 77.07% of assembled unigenes (Additional file 3). Of these unigenes, 79, 389 (97.3%) members shared >90% identity with their homologs from P. trichocarpa. Meanwhile, 3,6559 homologs (>80%) of P. trichocarpa v3 gene models have been sequenced. All these results indicate that our RNA-sequencing data has high contiguity, coverage, and could be used for further analyses.
Functional annotation and classification of the unigenes
To obtain a better understanding of the biological functions of the unigenes, a KEGG pathway-based analysis was also performed. Based on a comparison against the KEGG database using BLASTx with an E-value cutoff of <1.0E-5, 39,313 (36.23%) of the 108,502 unigenes had significant matches in the database and were assigned to 127 KEGG pathways. Of the 8,220 metabolism pathway unigenes, 2,726 were involved in plant hormone signal transduction pathways, including tryptophan metabolism, zeatin biosynthesis, diterpenoid biosynthesis, carotenoid biosynthesis, cysteine and methionine metabolism, brassinosteroid biosynthesis, α-Linolenic acid metabolism, and phenylalanine metabolism.
The three samples had 68 members in common when the 100 most abundant transcripts were compared (Additional file 5). The 23 unique members highly expressed in the control were involved in auxin signaling, cell division, and biogenesis. In contrast, the 19 unique members highly expressed in the C4 sample were stress (e.g., arginine decarboxylase, and dehydration) -induced genes. The 28 unique members highly expressed in the F4 sample were also stress-related genes, e.g., the glucanase, zinc finger protein, and E3 ubiquitin-protein ligase genes. These results indicate that our data are reliable.
Protein coding sequence prediction
Differentially expressed gene among three samples
We measured gene expression levels based on fragments per kilobase of exon model per million mapped reads (FPKM). After applying the chi-square test and Benjamini-Hochberg multiple testing corrections using R program among three samples simultaneously, we identified 2,858 genes as reliable DEGs in at least two samples (assigned as either DEGs) regardless of fold change (Additional file 6). Of these DEGs, 131 were expressed differentially in all three samples (assigned as all-DEGs, Additional file 7). Given a standard at an estimated absolute log2-fold change of >1, the respective DEGs of CK vs. C4, CK vs. F4 and C4 vs. F4 were 1,661, 866, and 1,161 (Additional files 8, 9, 10). The number of up-regulated unigenes in C4 and F4 samples was 1,113 and 630, respectively.
To accurately identify DEGs, we selected the 50 most significantly up-regulated transcripts that could be well-annotated by poplar database or NR database. As a result, those coding for the chlorophyll a/b binding protein (e.g., Unigene50811, CL12828.Contig3, Unigene50363, and Unigene55266), rubisco activase (CL4046.Contig4, Unigene50527, and Unigene55538), AP2/ERF transcription factors (Unigene26311,Unigene22719, CL9386.Contig2, Unigene18453, and CL9876.Contig3), and some other transcription factors (CL1721.Contig8, and Unigene27837) were the most up-regulated interpretable transcripts in C4 sample (Additional file 11). As for the top 50 up-regulated transcripts in the F4 sample (Additional file 12), the annotated transcripts focused on transcription factors (DREB1 transcription factors e.g. unigene26567 and unigene26311; WRKY transcription factor Unigene18620) and xyloglucan endotransglycosylases (Unigene19292, Unigene14078 and CL29.Contig1).
Gene ontology and pathway enrichment analyses of differentially expressed unigenes
Over-representative GO terms of DEGs in low temperature stressed P. euphratica GO ID
DEGs in 4°C treatment
response to stress
response to organic substance
response to stimulus
response to carbohydrate stimulus
response to endogenous stimulus
response to hormone stimulus
organic substance metabolic process
starch metabolic process
response to chemical stimulus
glucan metabolic process
regulation of metabolic process
regulation of macromolecule metabolic process
transferase activity, transferring hexosyl groups
structural molecule activity
DEGs in -4°C treatment
response to stress
response to stimulus
response to carbohydrate stimulus
structural molecule activity
enzyme inhibitor activity
external encapsulating structure
regulation of gene expression
regulation of macromolecule metabolic process
regulation of metabolic process
response to abiotic stimulus
cellular macromolecule biosynthetic process
macromolecule biosynthetic process
regulation of biological process
By performing the KEGG pathway analyses, we identified twelve pathways that changed significantly (q ≤ 0.05) under 4°C treatment, including the members involved in carbohydrate, energy, vitamin, hormone, and nitrogen metabolism (Additional file 13). “Plant pathogen interaction”, “hormone signal transduction”, and “biosynthesis of unsaturated fatty” pathways had the top three most differentially expressed unigene numbers and thus seem to play important roles in low temperature perception above freezing point. As for -4°C treatment, only 3 pathways changed significantly (Additional file 13). “Plant pathogen interaction” was assigned as a major pathway that changed significantly in both treated groups, indicating that low temperature stress response signal network may overlap with plant-pathogen interactions signals in P. euphratica. This is a notable finding considering that little is known about the overlap in signal transduction between abiotic and biotic stresses. Additionally, the transcripts of all of unsaturated fatty acid pathway genes increased significantly in the C4 sample. This result is in accordance with previous reports that plants undergoing low temperature stress preferentially accumulate poly-unsaturated and unsaturated fatty acids, which enhance low temperature tolerance under chilling conditions [24, 25].
Transcription factors responding to low temperature stress
Transcription factors play crucial roles in the regulation of target gene expression via specific binding to cis-acting elements in their promoters . Many of the COR genes contain cis-elements, such as dehydration-responsive elements/C-repeat elements (DRE/CRT, A/GCCGAC) and myeloblastosis (MYB, C/TAACNA/G) [27, 28] in their promoters that can be regulated by DREB and MYB transcription factors. Analysis of these transcription factors could provide useful information on the complex regulatory networks involved in P. euphratica cold stress responses.
Distribution of differentially expressed transcription factors
Previous studies have shown that not all cold-regulated gene expression is under the direct control of the CBF/DREB family [11, 29]. Besides the AP2/ERF family, it is likely that the WRKY and NAC transcription factors also play important roles in the transcriptional regulation of genes in early cold response in P. euphratica because they were overrepresented in the up-regulated gene list (log2-fold change > 1). In the WRKY family, 20 and 12 members were up-regulated in the C4 and F4 groups, respectively (Table 3). In comparison, none was found down-regulated in the respective groups. In the NAC transcription factor family, the transcripts of 9 and 5 members were up-regulated, while none was found down-regulated in both treated samples. Evidence that the WRKY and NAC transcription factor gene families may play important roles in the regulation of transcriptional reprogramming associated with cold stress responses is incremental [30–33]. The role of several WRKY genes, such as AtWRKY25, AtWRKY33, HvWRKY38, and GmWRKY21 [34–36], and several NAC genes, such as TaNAC2, SsNAC23, and OsNAC5 [30, 37, 38], in plant cold response has been established. The overrepresentation (log2-fold change > 1, FDR ≤ 0.05) of these two transcription factor families in the up-regulated gene list may indicate their importance in P. euphratica cold perception.
It is notable that transcription factors occupied about 20% members (9 out of 50 in C4, and 11 out of 50 in F4) in top 50 up-regulated gene lists (Additional files 11, 12). However, our analysis also showed that some transcription factors were not always up-regulated. For example, HD-ZIP transcription factors from both samples were found down-regulated. Some down-regulated members were also identified in C3H and bHLH families. According to previous reports, the bHLH family is involved in long term responses to low temperatures . Some ICE1-like bHLH transcription factors may be involved in CBF regulation . Here, several members of this transcription factor family have altered their expression in C4 sample but not in F4 sample. It seems likely that some bHLH transcription factors may play important roles in response of low temperature stress above freezing point. These results imply that the mechanism by which these genes mediate low temperature perception is complex.
Photosystem response to low temperature
The eventual effect of an abiotic stress on plant growth not only depends on the extent of the damage, but also on the capacity for recovery after the damage occurs . Cold-tolerant plants usually show a lower decrease in the rates of photosynthetic electron transport and photosynthetic carbon metabolism [42, 43]. These changes facilitate better recovery from chilling stress in these plants compared with cold-sensitive genotypes .
GO enrichment results showed that “carbon fixation” was one of the main biological processes after 4°C treatment. It was consistent with the over-representation of numbers of Ribulose bisphosphate carboxylase/oxygenase (Rubisco) genes  listed in the top 50 up-regulated unigenes of C4 sample. However, none was found in the top 50 up-regulated unigenes of F4 sample. Thus, we conclude that the expression elevation of Rubisco genes under low temperature below freezing point may be part of an adaptive mechanism that promotes P. euphratica survival.
Chilling-induced photoinhibition and subsequent recovery has been studied in maize. No changes or even increases in the photosystem I (PSI) activity have been observed . As for the photosystem II (PSII) system, an increase in proportion of its inactive centers and a lower capacity for repair of its damaged centers were observed and these changes were assumed to have led to persistent depression of photosynthetic efficiency . According to our data, a number of genes encoding photosystem I or II protein complexes exhibited up-regulation under 4°C treatment, while only a few of them up-regulated after -4°C treatment. Differential expression patterns of photosystem genes in the C4 and F4 samples suggest that photosystem could percept low temperature above but not below freezing point. This result is consistent with the findings that the photosystm of cold-acclimated plants was less affected than that of the nonacclimated plants .
ABA signal transduction components in P. euphratica cold response
Current evidence suggests that multiple mechanisms are involved in activating the cold-acclimation response. Both ABA-dependent and -independent pathways could regulate cold-responsive genes. ABA is an important stress hormone that mediates abiotic stress responses in plants . Plants grown under cold stress may experience altered hormone homeostasis and/or signal transduction . Chen et al.  found that ABA levels increased in cold-treated plants. The exogenous application of ABA at normal temperature has been found to enhance the freezing tolerance and thus aid in cold-acclimation in several plant species. However, whether ABA has a fundamental role in cold sensing and cold responsive network is unresolved.
To investigate whether ABA is involved in cold perception, we first investigated expression patterns of ABA biosynthesis genes annotated by NR or poplar transcript annotation. Sharp increase in fold-change of two unigenes (CL5772.Contig2 and CL5772.Contig1) encoding the key rate-limiting enzymes of ABA synthesis (namely, 9-cis epoxycarotenoid dioxygenase or NCED1) were observed in both samples . These results are consistent with those obtained in Arabidopsis in previous studies, suggesting that ABA biosynthesis may be involved in cold response [51, 52]. Recently, clade A protein phosphatases type 2C (PP2Cs) function as key negative regulators of ABA signaling by interacting with ABA receptors has been established [53, 54]. Here, all candidate PP2C genes identified in C4 and F4 samples showed an up-regulated expression pattern (Additional files 8, 9). Taken together, our results suggest ABA may function as an important mediator of low temperature perception and promoter of chilling and freezing tolerance in P. euphratica.
Calcium signal transduction components in P. euphratica cold response
As a second universal messenger involved in abiotic stress responses in plants, there is increasing evidence that calcium functions as an important messenger in a low temperature signal transduction pathway [1, 55]. In Arabidopsis, cytoplasmic calcium levels increase rapidly in response to low temperatures, largely due to an influx of calcium from extracellular stores . Calcium signals are perceived by calcium sensors that relay the signals into downstream targets such as phosphorylation cascades and gene expression regulation . Three major families of Ca2+ sensors function in abiotic stress responses in higher plants: calmodulins (CaMs) and CaM-like proteins, calcineurin B-like (CBL) proteins, and calcium-dependent protein kinases (CDPK) [58, 59]. CaMs and CBLs are small proteins that transmit the Ca2+ signal by interacting with target proteins and regulating their activity. CDPKs are monomeric proteins containing a CaM-like domain with four EF-hand motifs .
To investigate whether these Ca2+ sensor are involved in cold perception of P. euphratica, we searched for calcium related genes in differentially expressed gene list. As a result, we have found 41 members including CDPK homologs, calmodulin (CaM)-binding protein genes, CBL-interacting protein kinases (CIPKs), and Ca2+-ATPase up-regulated (Additional file 8). Only 3 of them showed down-regulated after 4°C treatment. As for -4°C treatment, 19 out of 21 calcium related genes showed up-regulated while only 2 were down-regulated (Additional file 9). These results are consistent with previous studies on Arabidopsis and rice showing that low temperature enhances the gene expression of Ca2+ sensors such as Ca2+-ATPase, CaMs, and CDPKs [61, 62]. Surprisingly, none of Ca2+ sensors from the CBL family showed up-regulated in both samples although most of its target kinases (e.g., CIPK10 homolog CL6537.contig1, CIPK6 homolog unigene17612, CIPK11 homolog unigene35076, and CIPK5 homolog unigene29689) were up-regulated. Given that the overexpression of CBL1 could result in attenuation of CBF/DREB induction , we speculate that the negative regulatory mechanism of CBL in cold perception may exist in P. euphratica. The complex cold responsive network composed of CBL calcium sensors and their target kinases require detailed analysis at the molecular level. Our results indicate that Ca2+ binding proteins are one of the signaling components induced at early stages of low temperature stress.
P. euphratica is a perennial desert tree with a high capacity for low temperature resistance compared with other poplar species . Here, we profiled the P. euphratica transcriptome under low temperature treatments in order to understand cold perception in this species. We obtained 108,502 assembled unigenes using the Trinity de novo assembly method and identified numerous potential cold sensing transcription factor genes and various key signal transduction components at the transcriptome level. These data provide useful resources for gene mining to improve low temperature stress tolerance in plants. Our results indicate not only that the CBF orthologs play key roles, but also that unsaturated fatty acids, ABA and calcium signaling components underlie a rapid and flexible cold response mechanism in P. euphratica. Compared with transient chilling stress response, fewer genes were found to have altered expression under temperature below freezing point. Particularly, only 3 pathways (q ≤ 0.05) clustered significantly with DGEs under freezing treatment. Taken together, our results provide valuable information on the molecular adaptation of P. euphratica to low temperature stress.
Stress treatment and RNA isolation
Two year-old P. euphratica seedlings, obtained from the Xinjiang Autonomous Region of China, were planted in individual pots (15 L) containing loam soil and placed in a greenhouse at Beijing Forestry University. Each pot contained four individuals. The temperature in the greenhouse was 22°C – 25°C with a 16 hour photoperiod (110–150 μmol.m-2.s-1). Potted plants were irrigated at a three-day interval based on evaporation demand for two months before treatment. In the chilling stress treatment group, P. euphratica plants were subjected to 4°C temperature under weak light. Simultaneously, some P. euphratica plants were subjected to a further drop of 8°C (-4°C) temperature under weak light conditions. The samples were harvested 6 h later and stored at -80°C until RNA extraction. Untreated P. euphratica plants harvested at 22°C under weak light was used as controls.
Total RNA was extracted from leaves using the CTAB method , treated with RNAase free DNAase. The A260/A280 ratios of the RNA samples were examined by NanoDrop 2000 and those with values ranging from 1.9 to 2.1 were chosen. The integrity of the RNA samples was examined with an Agilent 2100 Bioanalyzer. All RNA integrity number (RIN) values were over 8.0 and no signs of degradation were observed. Approximately 25 μg RNA sample with a concentration of ≥ 750 ng/μl from the CK, C4, and F4 samples were used for cDNA library construction.
Illumina cDNA library preparation and sequencing
An Illumina Hiseq 2000 library was constructed for Solexa sequencing. Poly (A) mRNA was first enriched with oligo (dT), and then was fragmented into small pieces of 200–700 bp using divalent cations at an elevated temperature. Based on these cleaved RNA fragments, we used random hexamer-primer and reverse transcriptase (Invitrogen) to synthesize cDNA. Three paired-end cDNA libraries with an insert size of 200 bp were constructed, and then sequenced with the Illumina HiSeq™ 2000 to generate an average read length of 90 bp as raw data.
De novo assembly and assessment
Raw data generated from Solexa sequencing were preprocessed to remove nonsense sequences including (1) adapters that were added for reverse transcription and sequencing, (2) sequences containing too many unknown bases (>5%), and (3) low-quality bases (>50% of the bases with a quality score ≤5). The preprocessed sequences were then assembled by Trinity  to construct contigs with default or optimal parameter. The contigs were then realigned to construct unigenes by Trinity. To fill the intra-scaffold gaps, we then used the paired-end information to retrieve read pairs that had one read well-aligned on the contigs and another read located in the gap region, and then locally assembled the collected reads. After gap closure, we constructed a non-redundant Unigene set from all three assembled datasets using the EST assembly program TGICL . All clean reads were mapped back to the assembled transcripts using standalone Bowtie v0.12.8 and the percentage of reads that used in assembling was calculated. To decide the sequential orientation for each unigene, a set of sequential BLASTx searches against GenBank’s non-redundant database, the Swiss-Prot protein database (http://www.expasy.ch/sprot), the KEGG pathways database , and the COG database (http://www.ncbi.nlm.nih. gov/COG/) were carried out. The sequential orientation orders of the unigenes not found in any of the four databases were deciphered by ESTScan software .
To assess the quality of the de novo assembly, a similarity search against P. trichocarpa transcripts (version 3) was further conducted using BLASTN algorithm with E-value less than 1.0E-5. The unigenes with identities higher than 75% were listed. The P. trichocarpa transcripts (version 3) and annotated gene set were downloaded from the Department of Energy Joint Genome Institute (DOE-JGI; http:://www.phytozome.net/poplar). All unigenes were also compared with data in publicly available databases including 13,991 unassembled P. euphratica ESTs from GenBank and assembled ESTs from other poplar species in the TIGR Plant Transcript Assemblies database.
Unigene annotation and function classification
To find the most descriptive annotation for each sequence, all assembled unique sequences were searched against NR, Swiss-Prot, KEGG, and COG using the BLASTx and against Nt using the BLASTn algorithms with a threshold of E < 1.0E-5. The protein with the highest sequence similarity was retrieved. For NR annotation, Blast2GO software  was first used to perform gene annotation, and then the WEGO and AGRIGO software  were used to conduct GO functional classification to interpret the distribution of gene functions defined by molecular function, cellular component, and biological process ontologies. Unigene sequences were also aligned to the COG database to predict and classify gene functions. Pathway assignments were carried out according to the KEGG pathway database using BLASTx with an E-value threshold of 1.0E-5.
Protein-coding region prediction and transcription factor analysis
For protein coding sequence prediction, unigenes were searched against the NR, Swiss-Prot, KEGG, and COG protein databases in that order using a BLASTx algorithm (E-value < 1.0E-5). Unigenes that aligned to a high priority database were not aligned to databases of lower priority. The correct reading frame of the nucleotide sequences (5’-3’direction) of unigenes was defined by the highest rank in the BLAST results, and the corresponding protein sequences were obtained from the standard codon table.
Transcription factors were predicted according to protein sequences obtained from CDS prediction. We used hmmsearch to search for domains of the plant transcription factors (http://plntfdb.bio.uni-potsdam.de/v3.0/) and classified unigenes according to the gene family information. Genes that were believed to be associated with cold stress were selected for further investigation.
Gene expression analysis
where C is the number of reads that uniquely aligned to one unigene; N is the total number of reads that uniquely aligned to all unigenes; L is the base number in the CDS of one unigene. The goal of this transformation is to normalize the counts in regard of the differing library sizes and the length of the transcripts .
We identified DEGs from different samples (CK, C4 and F4) using R program. The Pearson’s chi-square test was applied to assess the lane effect. For each gene, the P-value was computed. After that, Benjamini–Hochberg false discovery rate (FDR) was applied to correct the results for p value. FDR method is widely used in deep-sequencing studies because of its power in finding over-representative unigenes [70–73]. The transcripts that were induced or suppressed at an estimated absolute log2-fold change of >1 and FDR adjusted p-value ≤ 0.05 were considered to be differentially expressed .
GO and KEGG analyses for differentially expressed unigenes
In order to find the significantly enriched GO terms in DEGs against a genome background, the DEGs were annotated to GO database (http://www.geneontology.org/) using hypergeometric test for statistical analysis . For p value correction, we used the rigorous Bonferroni correction method. The cutoff p value after correction was 0.05. GO terms fulfilling this condition were defined as being significantly enriched. The KEGG pathway enrichment analysis of DEGs was also performed with the whole genome background as a reference to find the main biochemical pathways and signal transduction pathways in which DEGs involved. After multiple testing corrections, we defined pathways with q-value ≤0.05 as being those significantly enriched in DEGs.
Quantitative PCR analysis
Quantitative PCR (qPCR) was performed to determine the expression level of selected unigenes. The qPCR was conducted using a power SYBR Green PCR Kit (ABI) in a MicroAmp™ 96-well plate with a StepOnePlus™ Real-Time PCR System (ABI). The relative expression value was calculated by the 2-ΔΔCt method using PeActin (GenBank accession number EF148840) as an internal control . Gene-specific primers used in the qPCR analysis are listed in Additional file 14. RNA pools used in the qPCR analyses were extracted from three independent samples which were different from those used for RNA-seq. Three technical replicates were used for each sample.
CRT binding transcription factor
Differentially expressed gene
False discovery rate
Fragments per kilobase of exon region in a given gene per million mapped fragments
Kyoto Encyclopedia of Genes and Genomes
Quantitative real-time PCR
We thank the technical assistance of Chu-Yu Ye, Chao Liu and Sha Tang from College of Biological Sciences and technology in Beijing Forestry University. We also appreciate the assistance of Hongqing Xie and Tong Xu from Beijing Genomics Institute in data assessment. This work was supported by Beijing Higher Education Young Elite Teacher Project (YETP0754), National Natural Science Foundation of China (31100492 and 31370597), Fundamental Research funds for the Central Universities (DT2012-01 and YX2011-22), National R&D Key Project of China (2011BAD38B01), and Research Fund for the Doctoral Program of Higher Education of China (20110014120002).
- Thomashow MF: Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999, 50: 571-599. 10.1146/annurev.arplant.50.1.571.PubMedView Article
- Conesa A, Gotz S, Garcia-Gomez JM, Terol J, Talon M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinfo. 2005, 21 (18): 3674-3676. 10.1093/bioinformatics/bti610.View Article
- Korn M, Peterek S, Mock HP, Heyer AG, Hincha DK: Heterosis in the freezing tolerance, and sugar and flavonoid contents of crosses between Arabidopsis thaliana accessions of widely varying freezing tolerance. Plant Cell Environ. 2008, 31 (6): 813-827. 10.1111/j.1365-3040.2008.01800.x.PubMed CentralPubMedView Article
- Hannah MA, Wiese D, Freund S, Fiehn O, Heyer AG, Hincha DK: Natural genetic variation of freezing tolerance in Arabidopsis. Plant Physiol. 2006, 142 (1): 98-112. 10.1104/pp.106.081141.PubMed CentralPubMedView Article
- Viswanathan C, Zhu JK: Molecular genetic analysis of cold-regulated gene transcription. Philos Trans R Soc Lond B Biol Sci. 2002, 357 (1423): 877-886. 10.1098/rstb.2002.1076.PubMed CentralPubMedView Article
- Carvallo MA, Pino MT, Jeknic Z, Zou C, Doherty CJ, Shiu SH, Chen TH, Thomashow MF: A comparison of the low temperature transcriptomes and CBF regulons of three plant species that differ in freezing tolerance: Solanum commersonii, Solanum tuberosum, and Arabidopsis thaliana. J Exp Bot. 2011, 62 (11): 3807-3819. 10.1093/jxb/err066.PubMed CentralPubMedView Article
- Jaglo-Ottosen KR, Gilmour SJ, Zarka DG, Schabenberger O, Thomashow MF: Arabidopsis CBF1 overexpression induces COR genes and enhances freezing tolerance. Science. 1998, 280 (5360): 104-106. 10.1126/science.280.5360.104.PubMedView Article
- Hsieh TH, Lee JT, Yang PT, Chiu LH, Charng YY, Wang YC, Chan MT: Heterology expression of the Arabidopsis C-repeat/dehydration response element binding factor 1 gene confers elevated tolerance to chilling and oxidative stresses in transgenic tomato. Plant Physiol. 2002, 129 (3): 1086-1094. 10.1104/pp.003442.PubMed CentralPubMedView Article
- Lee BH, Henderson DA, Zhu JK: The Arabidopsis cold-responsive transcriptome and its regulation by ICE1. Plant Cell. 2005, 17 (11): 3155-3175. 10.1105/tpc.105.035568.PubMed CentralPubMedView Article
- Chinnusamy V, Zhu JK, Sunkar R: Gene regulation during cold stress acclimation in plants. Methods Mol Biol. 2010, 639: 39-55. 10.1007/978-1-60761-702-0_3.PubMed CentralPubMedView Article
- Fowler S, Thomashow MF: Arabidopsis transcriptome profiling indicates that multiple regulatory pathways are activated during cold acclimation in addition to the CBF cold response pathway. Plant Cell. 2002, 14 (8): 1675-1690. 10.1105/tpc.003483.PubMed CentralPubMedView Article
- Browicz K: Chorology of populus euphratica Olivier. Arboretum Kornickie. 1977, 22: 5-27.
- Sun J, Wang MJ, Ding MQ, Deng SR, Liu MQ, Lu CF, Zhou XY, Shen X, Zheng XJ, Zhang ZK, Song J, Hu ZM, Xu Y, Chen SL: H2O2 and cytosolic Ca2+ signals triggered by the PM H+-coupled transport system mediate K+/Na+ homeostasis in NaCl-stressed Populus euphratica cells. Plant Cell Environ. 2010, 33 (6): 943-958. 10.1111/j.1365-3040.2010.02118.x.PubMedView Article
- Ye CY, Zhang HC, Chen JH, Xia XL, Yin WL: Molecular characterization of putative vacuolar NHX-type Na+/H+ exchanger genes from the salt-resistant tree Populus euphratica. Physiol Plant. 2009, 137 (2): 166-174. 10.1111/j.1399-3054.2009.01269.x.PubMedView Article
- Wu Y, Ding N, Zhao X, Zhao M, Chang Z, Liu J, Zhang L: Molecular characterization of PeSOS1: the putative Na+/H+ antiporter of Populus euphratica. Plant Mol Biol. 2007, 65 (1–2): 1-11.PubMedView Article
- Ottow EA, Brinker M, Teichmann T, Fritz E, Kaiser W, Brosche M, Kangasjarvi J, Jiang X, Polle A: Populus euphratica displays apoplastic sodium accumulation, osmotic adjustment by decreases in calcium and soluble carbohydrates, and develops leaf succulence under salt stress. Plant Physiol. 2005, 139 (4): 1762-1772. 10.1104/pp.105.069971.PubMed CentralPubMedView Article
- Gu R, Fonseca S, Puskas LG, Hackler L, Zvara A, Dudits D, Pais MS: Transcript identification and profiling during salt stress and recovery of Populus euphratica. Tree Physiol. 2004, 24 (3): 265-276. 10.1093/treephys/24.3.265.PubMedView Article
- Wilhelm BT, Marguerat S, Watt S, Schubert F, Wood V, Goodhead I, Penkett CJ, Rogers J, Bahler J: Dynamic repertoire of a eukaryotic transcriptome surveyed at single-nucleotide resolution. Nature. 2008, 453 (7199): 1239-1243. 10.1038/nature07002.PubMedView Article
- Wu T, Qin Z, Zhou X, Feng Z, Du Y: Transcriptome profile analysis of floral sex determination in cucumber. J Plant Physiol. 2010, 167 (11): 905-913. 10.1016/j.jplph.2010.02.004.PubMedView Article
- Wang XW, Luan JB, Li JM, Bao YY, Zhang CX, Liu SS: De novo characterization of a whitefly transcriptome and analysis of its gene expression during development. BMC Genomics. 2010, 11: 400-10.1186/1471-2164-11-400.PubMed CentralPubMedView Article
- Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen ZH, Mauceli E, Hacohen N, Gnirke A, Rhind N, Palma FD, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A: Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol. 2011, 29 (7): 644-652. 10.1038/nbt.1883.PubMed CentralPubMedView Article
- Iseli C, Jongeneel CV, Bucher P: ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proc Int Conf Intell Syst Mol Biol. 1999, 7: 138-148.
- Zhao Z, Tan L, Dang C, Zhang H, Wu Q, An L: Deep-sequencing transcriptome analysis of chilling tolerance mechanisms of a subnival alpine plant, Chorispora bungeana. BMC Plant Biol. 2012, 12: 222-10.1186/1471-2229-12-222.PubMed CentralPubMedView Article
- Kodama H, Hamada T, Horiguchi G, Nishimura M, Iba K: Genetic enhancement of cold tolerance by expression of a gene for chloroplast [omega]-3 fatty acid desaturase in transgenic tobacco. Plant Physiol. 1994, 105 (2): 601-605.PubMed CentralPubMed
- Williams JP, Khan MU, Mitchell K, Johnson G: The effect of temperature on the level and biosynthesis of unsaturated fatty acids in diacylglycerols of Brassica napus leaves. Plant Physiol. 1988, 87 (4): 904-910. 10.1104/pp.87.4.904.PubMed CentralPubMedView Article
- Yamaguchi-Shinozaki K, Shinozaki K: Transcriptional regulatory networks in cellular responses and tolerance to dehydration and cold stresses. Annu Rev Plant Biol. 2006, 57: 781-803. 10.1146/annurev.arplant.57.032905.105444.PubMedView Article
- Stockinger EJ, Gilmour SJ, Thomashow MF: Arabidopsis thaliana CBF1 encodes an AP2 domain-containing transcriptional activator that binds to the C-repeat/DRE, a cis-acting DNA regulatory element that stimulates transcription in response to low temperature and water deficit. Proc Natl Acad Sci U S A. 1997, 94 (3): 1035-1040. 10.1073/pnas.94.3.1035.PubMed CentralPubMedView Article
- Yamaguchi-Shinozaki K, Shinozaki K: Organization of cis-acting regulatory elements in osmotic- and cold-stress-responsive promoters. Trends Plant Sci. 2005, 10 (2): 88-94. 10.1016/j.tplants.2004.12.012.PubMedView Article
- Kreps JA, Wu Y, Chang HS, Zhu T, Wang X, Harper JF: Transcriptome changes for Arabidopsis in response to salt, osmotic, and cold stress. Plant Physiol. 2002, 130 (4): 2129-2141. 10.1104/pp.008532.PubMed CentralPubMedView Article
- Hu H, You J, Fang Y, Zhu X, Qi Z, Xiong L: Characterization of transcription factor gene SNAC2 conferring cold and salt tolerance in rice. Plant Mol Biol. 2008, 67 (1–2): 169-181.PubMedView Article
- Şahin-Çevik M: A WRKY transcription factor gene isolated from Poncirus trifoliata shows differential responses to cold and drought stresses. Plant Omics J. 2012, 5 (438): 45-
- Wang JY, Wang JP, He Y: A Populus euphratica NAC protein regulating Na(+)/K(+) homeostasis improves salt tolerance in Arabidopsis thaliana. Gene. 2013, 521 (2): 265-273. 10.1016/j.gene.2013.03.068.PubMedView Article
- Aslam M, Grover A, Sinha VB, Fakher B, Pande V, Yadav PV, Gupta SM, Anandhan S, Ahmed Z: Isolation and characterization of cold responsive NAC gene from Lepidium latifolium. Mol Biol Rep. 2012, 39 (10): 9629-9638. 10.1007/s11033-012-1828-0.PubMedView Article
- Marè C, Mazzucotelli E, Crosatti C, Francia E, Cattivelli L: HvWRKY38: a new transcription factor involved in cold-and drought-response in barley. Plant Mol Biol. 2004, 55 (3): 399-416. 10.1007/s11103-004-0906-7.PubMedView Article
- Zhou QY, Tian AG, Zou HF, Xie ZM, Lei G, Huang J, Wang CM, Wang HW, Zhang JS, Chen SY: Soybean WRKY-type transcription factor genes, GmWRKY13, GmWRKY21, and GmWRKY54, confer differential tolerance to abiotic stresses in transgenic Arabidopsis plants. Plant Biotechnol J. 2008, 6 (5): 486-503. 10.1111/j.1467-7652.2008.00336.x.PubMedView Article
- Jiang Y, Deyholos MK: Functional characterization of Arabidopsis NaCl-inducible WRKY25 and WRKY33 transcription factors in abiotic stresses. Plant Mol Biol. 2009, 69 (1–2): 91-105.PubMedView Article
- Nogueira FT, Schlögl PS, Camargo SR, Fernandez JH: SsNAC23, a member of the NAC domain protein family, is associated with cold, herbivory and water stress in sugarcane. Plant Sci. 2005, 169 (1): 93-106. 10.1016/j.plantsci.2005.03.008.View Article
- Song SY, Chen Y, Chen J, Dai XY, Zhang WH: Physiological mechanisms underlying OsNAC5-dependent tolerance of rice plants to abiotic stress. Planta. 2011, 234 (2): 331-345. 10.1007/s00425-011-1403-2.PubMedView Article
- Hannah MA, Heyer AG, Hincha DK: A global survey of gene regulation during cold acclimation in Arabidopsis thaliana. PLoS Gen. 2005, 1 (2): e26-10.1371/journal.pgen.0010026.View Article
- Chinnusamy V, Ohta M, Kanrar S, Lee B-h, Hong X, Agarwal M, Zhu J-K: ICE1: a regulator of cold-induced transcriptome and freezing tolerance in Arabidopsis. Genes Dev. 2003, 17 (8): 1043-1054. 10.1101/gad.1077503.PubMed CentralPubMedView Article
- Zhang S, Scheller HV: Photoinhibition of photosystem I at chilling temperature and subsequent recovery in Arabidopsis thaliana. Plant Cell Physiol. 2004, 45 (11): 1595-1602. 10.1093/pcp/pch180.PubMedView Article
- Ribas-Carbo M, Aroca R, Gonzalez-Meler MA, Irigoyen JJ, Sanchez-Diaz M: The electron partitioning between the cytochrome and alternative respiratory pathways during chilling recovery in two cultivars of maize differing in chilling sensitivity. Plant Physiol. 2000, 122 (1): 199-204. 10.1104/pp.122.1.199.PubMed CentralPubMedView Article
- Pietrini F, Iannelli M, Battistelli A, Moscatello S, Loreto F, Massacci A: Effects on photosynthesis, carbohydrate accumulation and regrowth induced by temperature increase in maize genotypes with different sensitivity to low temperature. Funct Plant Biol. 1999, 26 (4): 367-373.
- Aroca R, Irigoyen JJ, Sánchez-Dı́az M: Photosynthetic characteristics and protective mechanisms against oxidative stress during chilling and subsequent recovery in two maize varieties differing in chilling sensitivity. Plant Sci. 2001, 161 (4): 719-726. 10.1016/S0168-9452(01)00460-5.View Article
- Ellis RJ: The most abundant protein in the world. Trends Biochem Sci. 1979, 4: 241-244. 10.1016/0968-0004(79)90212-3.View Article
- Hola D, Kocova M, Rothova O, Wilhelmova N, Benesova M: Recovery of maize (Zea mays L.) inbreds and hybrids from chilling stress of various duration: photosynthesis and antioxidant enzymes. J Plant Physiol. 2007, 164 (7): 868-877. 10.1016/j.jplph.2006.04.016.PubMedView Article
- Fryer MJ, Oxborough K, Martin B, Ort DR, Baker NR: Factors associated with depression of photosynthetic quantum efficiency in Maize at low growth temperature. Plant Physiol. 1995, 108 (2): 761-767.PubMed CentralPubMed
- Kuk YI, Lee JH, Kim HY, Chung SJ, Chung GC, Guh JO, Lee HJ, Burgos NR: Relationships of cold acclimation and antioxidative enzymes with chilling tolerance in cucumber (Cucumis sativus L.). J Am Soc Hortic Sci. 2003, 128 (5): 661-666.
- Chen HH, Li PH, Brenner ML: Involvement of abscisic Acid in potato cold acclimation. Plant Physiol. 1983, 71 (2): 362-365. 10.1104/pp.71.2.362.PubMed CentralPubMedView Article
- Nitsch LM, Oplaat C, Feron R, Ma Q, Wolters-Arts M, Hedden P, Mariani C, Vriezen WH: Abscisic acid levels in tomato ovaries are regulated by LeNCED1 and SlCYP707A1. Planta. 2009, 229 (6): 1335-1346. 10.1007/s00425-009-0913-7.PubMedView Article
- Park HY, Seok HY, Park BK, Kim SH, Goh CH, Lee BH, Lee CH, Moon YH: Overexpression of Arabidopsis ZEP enhances tolerance to osmotic stress. Biochem Biophys Res Commun. 2008, 375 (1): 80-85. 10.1016/j.bbrc.2008.07.128.PubMedView Article
- Lang V, Mantyla E, Welin B, Sundberg B, Palva ET: Alterations in water status, endogenous abscisic acid content, and expression of rab18 Gene during the development of freezing tin Arabidopsis thaliana. Plant Physiol. 1994, 104 (4): 1341-1349.PubMed CentralPubMed
- Pandey S, Nelson DC, Assmann SM: Two novel GPCR-type G proteins are abscisic acid receptors in Arabidopsis. Cell. 2009, 136 (1): 136-148. 10.1016/j.cell.2008.12.026.PubMedView Article
- Ma Y, Szostkiewicz I, Korte A, Moes D, Yang Y, Christmann A, Grill E: Regulators of PP2C phosphatase activity function as abscisic acid sensors. Science. 2009, 324 (5930): 1064-1068.PubMed
- Ye CY, Yang X, Xia X, Yin W: Comparative analysis of cation/proton antiporter superfamily in plants. Gene. 2013, 521 (2): 245-251. 10.1016/j.gene.2013.03.104.PubMedView Article
- Monroy AF, Dhindsa RS: Low-temperature signal transduction: induction of cold acclimation-specific genes of alfalfa by calcium at 25 degrees C. Plant Cell. 1995, 7 (3): 321-331.PubMed CentralPubMed
- Gong D, Guo Y, Schumaker KS, Zhu J-K: The SOS3 family of calcium sensors and SOS2 family of protein kinases in Arabidopsis. Plant Physiol. 2004, 134 (3): 919-926. 10.1104/pp.103.037440.PubMed CentralPubMedView Article
- Luan S: The CBL-CIPK network in plant calcium signaling. Trends Plant Sci. 2009, 14 (1): 37-42. 10.1016/j.tplants.2008.10.005.PubMedView Article
- Sheen J: Ca2+-dependent protein kinases and stress signal transduction in plants. Science. 1996, 274 (5294): 1900-1902. 10.1126/science.274.5294.1900.PubMedView Article
- Hrabak EM, Chan CW, Gribskov M, Harper JF, Choi JH, Halford N, Kudla J, Luan S, Nimmo HG, Sussman MR, Thomas M, Walker-Simmons K, Zhu JK, Harmon AC: The Arabidopsis CDPK-SnRK superfamily of protein kinases. Plant Physiol. 2003, 132 (2): 666-680. 10.1104/pp.102.011999.PubMed CentralPubMedView Article
- Komatsu S, Yang G, Khan M, Onodera H, Toki S, Yamaguchi M: Over-expression of calcium-dependent protein kinase 13 and calreticulin interacting protein 1 confers cold tolerance on rice plants. Mol Genet Genomics. 2007, 277 (6): 713-723. 10.1007/s00438-007-0220-6.PubMedView Article
- Abbasi F, Onodera H, Toki S, Tanaka H, Komatsu S: OsCDPK13, a calcium-dependent protein kinase gene from rice, is induced by cold and gibberellin in rice leaf sheath. Plant Mol Biol. 2004, 55 (4): 541-552. 10.1007/s11103-004-1178-y.PubMedView Article
- Chen J, Xue B, Xia X, Yin W: A novel calcium-dependent protein kinase gene from Populus euphratica. Biochem Biophys Res Commun. 2013, 441 (3): 630-636. 10.1016/j.bbrc.2013.10.103.PubMedView Article
- Chang S, Puryear J, Cairney J: A simple and efficient method for isolating RNA from pine trees. Plant Mol Biol Rep. 1993, 11 (2): 113-116. 10.1007/BF02670468.View Article
- Pertea G, Huang X, Liang F, Antonescu V, Sultana R, Karamycheva S, Lee Y, White J, Cheung F, Parvizi B, Tsai J, Quackenbush J: TIGR Gene Indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinfo. 2003, 19 (5): 651-652. 10.1093/bioinformatics/btg034.View Article
- Kanehisa M, Araki M, Goto S, Hattori M, Hirakawa M, Itoh M, Katayama T, Kawashima S, Okuda S, Tokimatsu T, Yamanishi Y: KEGG for linking genomes to life and the environment. Nucleic Acids Res. 2008, 36 (Database issue): D480-D484.PubMed CentralPubMed
- Ye J, Fang L, Zheng H, Zhang Y, Chen J, Zhang Z, Wang J, Li S, Li R, Bolund L: WEGO: a web tool for plotting GO annotations. Nucleic Acids Res. 2006, 34 (Web Server issue): W293-W297.PubMed CentralPubMedView Article
- Mortazavi A, Williams BA, McCue K, Schaeffer L, Wold B: Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008, 5 (7): 621-628. 10.1038/nmeth.1226.PubMedView Article
- Soneson C, Delorenzi M: A comparison of methods for differential expression analysis of RNA-seq data. BMC Bioinfo. 2013, 14: 91-10.1186/1471-2105-14-91.View Article
- Junttila S, Laiho A, Gyenesei A, Rudd S: Whole transcriptome characterization of the effects of dehydration and rehydration on Cladonia rangiferina, the grey reindeer lichen. BMC Genomics. 2013, 14: 870-10.1186/1471-2164-14-870.PubMed CentralPubMedView Article
- Reiner A, Yekutieli D, Benjamini Y: Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinfo. 2003, 19 (3): 368-375. 10.1093/bioinformatics/btf877.View Article
- Pang T, Ye C-Y, Xia X, Yin W: De novo sequencing and transcriptome analysis of the desert shrub, Ammopiptanthus mongolicus, during cold acclimation using Illumina/Solexa. BMC Genomics. 2013, 14 (1): 488-10.1186/1471-2164-14-488.PubMed CentralPubMedView Article
- Lu T, Lu G, Fan D, Zhu C, Li W, Zhao Q, Feng Q, Zhao Y, Guo Y, Li W: Function annotation of the rice transcriptome at single-nucleotide resolution by RNA-seq. Genome Res. 2010, 20 (9): 1238-1249. 10.1101/gr.106120.110.PubMed CentralPubMedView Article
- Wu J, Zhang Y, Zhang H, Huang H, Folta KM, Lu J: Whole genome wide expression profiles of Vitis amurensis grape responding to downy mildew by using Solexa sequencing technology. BMC Plant Biol. 2010, 10: 234-10.1186/1471-2229-10-234.PubMed CentralPubMedView Article
- Chen J, Xia X, Yin W: A poplar DRE-binding protein gene, PeDREB2L, is involved in regulation of defense response against abiotic stress. Gene. 2011, 483 (1–2): 36-42.PubMedView Article
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited.