Open Access

Chilling acclimation provides immunity to stress by altering regulatory networks and inducing genes with protective functions in Cassava

  • Changying Zeng1,
  • Zheng Chen2, 3,
  • Jing Xia2, 3,
  • Kevin Zhang2, 3,
  • Xin Chen1,
  • Yufei Zhou1,
  • Weiping Bo1,
  • Shun Song1,
  • Deli Deng1,
  • Xin Guo1,
  • Bin Wang1,
  • Junfei Zhou2,
  • Hai Peng2,
  • Wenquan Wang1,
  • Ming Peng1Email author and
  • Weixiong Zhang2, 3, 4Email author
BMC Plant Biology201414:207

https://doi.org/10.1186/s12870-014-0207-5

Received: 10 February 2014

Accepted: 22 July 2014

Published: 5 August 2014

Abstract

Background

Stress acclimation is an effective mechanism that plants acquired for adaption to dynamic environment. Even though generally considered to be sensitive to low temperature, Cassava, a major tropical crop, can be tolerant to much lower temperature after chilling acclimation. Improvement to chilling resistance could be beneficial to breeding. However, the underlying mechanism and the effects of chilling acclimation on chilling tolerance remain largely unexplored.

Results

In order to understand the mechanism of chilling acclimation, we profiled and analyzed the transcriptome and microRNAome of Cassava, using high-throughput deep sequencing, across the normal condition, a moderate chilling stress (14°C), a harsh stress (4°C) after chilling acclimation (14°C), and a chilling shock from 24°C to 4°C. The results revealed that moderate stress and chilling shock triggered comparable degrees of transcriptional perturbation, and more importantly, about two thirds of differentially expressed genes reversed their expression from up-regulation to down-regulation or vice versa in response to hash stress after experiencing moderate stress. In addition, microRNAs played important roles in the process of this massive genetic circuitry rewiring. Furthermore, function analysis revealed that chilling acclimation helped the plant develop immunity to further harsh stress by exclusively inducing genes with function for nutrient reservation therefore providing protection, whereas chilling shock induced genes with function for viral reproduction therefore causing damage.

Conclusions

Our study revealed, for the first time, the molecular basis of chilling acclimation, and showed potential regulation role of microRNA in chilling response and acclimation in Euphorbia.

Keywords

Chilling acclimation Chilling shock Gene regulation Cassava Castor bean

Background

Plants have developed complex defense systems to cope with and combat against harsh environmental stress. After environmental stress, e.g., chilly temperature, most plants gain or increase stress tolerance, resulting in stress acclimation [1],[2]. Chilling acclimation is a favorable trait that is critical for plant growth, reproduction and survival. It helps defeat various chilling-related catastrophic outcomes, such as growth retardation, chlorosis, necrosis, and yield reduction [3]. Chilling resistance can be further classified into above-zero temperature tolerance for tropical plants and sub-zero temperature tolerance for temperates [4].

Cold stress can be classified as chilling (<20°C) and freezing (<0°C) stress [5]. Tropical plants can be injured by above-zero chilling temperature; chilling-injured leaves may become purple or reddish and in some cases wilt [6]. Chilling responses of tropical plants has been documented: an oxidative signaling regulatory network triggers an early response to chilling stress in Japonica rice [7], and exogenous ABA can induce freezing tolerance in chilling sensitive rice seedlings [8]. Further, global expression profiling of chilling induced genes in the chilling tolerant Japonica rice suggests a role of ABA signaling in chilling tolerance [9]. Fatty acid desaturation of organelle membrane lipids related ACYL-LIPID DESATURASE2 is required for chilling and freezing tolerance in Arabidopsis [10]. A chloroplast-targeted protein complex stabilization related DnaJ protein contributes to maintenance of photosystem II under chilling stress in tomato [11].

Chilling response and chilling acclimation may alter gene regulatory circuitry [12]. However, the scale and mechanism of reprogramming of gene regulation as well as the molecular components involved remain to be investigated in important tropical crops, such as Cassava. Moreover, although small noncoding RNAs (sncRNAs), particularly microRNAs (miRNAs), have been recognized as essential post-transcriptional gene regulators in plant development and stress responses [12]–[14], their functions in chilling acclimation have not yet been well documented.

Most previous studies on chilling stress focus on model plants, typically Arabidopsis and rice, whereas little has been done on Euphorbia, a genus of tropical plants. Many Euphorbiaceous plants, e.g., Cassava and castor bean, are agri-economically important. Cassava (Manihot esculenta) is a major source of carbohydrates for over 500 million people in the developing countries in the tropics and sub-tropics [15]. It is also a major source of industrial material, for biofuel production for example [15],[16]. Cassava is remarkably tolerant to drought and low-fertility soils. However, it is sensitive to low temperature, and chilling injury often occurs in spring planting and autumn harvest seasons.

In this genome-wide study we analyzed transcriptome variations of Cassava plants in response to chilling and during chilling acclimation, aiming at elucidating gene regulatory networks underlying chilling acclimation. Specifically we compared the gene expression variations in responses to dramatic temperature decreases and during chilling acclimation in reference to the normal growth condition. We profiled the expression of protein-coding genes and sncRNA species using Next Generation (NextGen) sequencing. By analyzing more than 35.3 million sequencing reads from 4 mRNA libraries and 25.6 million reads from 4 sncRNA libraries, we identified differentially expressed mRNA and miRNA genes, from which we further identified and analyzed mRNA and miRNA genes that are critical to chilling acclimation.

Results

Exploring chilling response and acclimation in Cassava

We profiled the transcriptome and microRNAome of SC124, a widely planted Cassava cultivar in China. SC124 is sensitive to chilling and can be exploited to study chilling response and acclimation. The profiling experiments were carried out under three chilling stress treatments (detailed in Methods and illustrated in Additional file 1: Figure S1): 1) gradual chilling acclimation (CA) where plants grown in the normal condition of 24°C were stressed to 14°C; 2) chilling stress after chilling acclimation (CCA) where plants after 5 days of the CA treatment were transferred further from 14°C to 4°C; and 3) chilling shock (CS) where plants were experienced a dramatic temperature drop from 24°C to 4°C (see Methods for detail). For comparison, plants grown continuously under 24°C were used as the normal control (NC). Total RNA was extracted from three organs/tissues of the plants at the 6 h, 24 h and 5d of the corresponding stress treatments and the normal control in order to account for initial and secondary responses as well as functional adaption to chilling stresses.

Distinct symptoms of chilling stress of CS and CCA were observed at the end of these stress experiments (Additional file 1: Figure S1). The CCA treated plants were more chilling resistant than the CS treated plants at 4°C: fewer leaves wilted and more leaves stayed upright. Four physiological traits were measured to further evaluate the impact of chilling stress (Additional file 1: Figure S2). While there was no statistically significant changes in leaf falling (Additional file 1: Figure S2A), chlorophyll content only decreased in CS condition (Additional file 1: Figure S2B), malondialdehyde content (Additional file 1: Figure S2C) and leaf proline content (Additional file 1: Figure S2D) increased after one of the three chilling treatments. Malondialdehyde content was progressively elevated with the severity of stress, from CA, CCA to CS. Proline content exhibited different variation patterns. It increased the most under CS and had the least variation under CCA (Additional file 1: Figure S2D).

Chilling stresses trigger significant transcriptome and microRNAome variations

In order to appreciate the impact of chilling stress, we profiled the expression of mRNA genes and small-noncoding RNA (sncRNA) species of SC124 after the CA, CCA and CS treatments as well as under the normal condition (NC) using Illumina GAIIx (see Methods, sequencing data in NCBI/GEO, accession # GSE52178). Overall ~33 million (>90%) of raw sequencing reads from mRNA genes, or RNA-seq reads, were high-quality (or qualified) reads, among which more than 80% could be mapped to the Cassava reference genome (http://phytozome.net) allowing one mismatch (Additional file 1: Table S1). With the criterion of at least 10 reads per million (RPM, see Methods), the mapped reads attributed to 12,689 (37.16% of the 34,151 annotated Cassava mRNA genes), 16,023 (46.92%), 15,144 (44.34%) and 17,026 (49.85%) mRNA genes expressed under the NC, CA, CCA, and CS conditions, respectively (Figure 1A and Additional file 1: Table S1). Comparing the numbers of these expressed genes showed that at least 19% more genes were expressed in any of the three chilling treatments than in the normal condition.
Figure 1

Results of transcriptome and microRNAome profiling and their variations across the normal condition (NC) and the three chilling stress treatments – chilling shock (CS), chilling acclimation (CA) and chilling after chilling acclimation (CCA). (A) The numbers of expressed and differentially expressed (DE) genes and miRNAs across the four conditions. The two numbers in an oval are the numbers of expressed genes (the first line) and miRNAs (the second line). The two numbers on an edge are the numbers of DE genes and DE miRNAs between two conditions. (B) Relationship between the up- and down-regulated genes and that of miRNAs of CS and CA with respect to NC. The two numbers within a region are the numbers of DE genes (the first line) and miRNAs (the second line). The figure shows a substantial overlapping between the up-regulated genes (and miRNAs) and overlapping between the down-regulated genes (and miRNAs) of CS and CA. (C) Similar to (B); relationship between the up- and down-regulated genes and that of miRNAs when going from NC to AC and then from AC to ACC. The figure shows a little overlapping between the up-regulated genes (and miRNAs) and overlapping between the down-regulated genes (and miRNAs) in the comparison.

Correspondingly, four small-RNA libraries for the NC, CA, CCA and CS conditions contributed more than 25.6 million raw small-RNA sequencing reads (see Methods, sequencing data in NCBI/GEO, accession # GSE52178), among which 23,468,606 (>91% of the total) were adapter-trimmed, high-quality reads (qualified reads, Additional file 1: Table S2A). Among the qualified reads, which had lengths peaked at 21- and 24-nt, 53.40% and 73.18% could be mapped to the Cassava reference genome allowing zero and one mismatch (Additional file 1: Figure S3 and Tables S2B and S2C), respectively. Based on a set of stringent criteria (see Methods), 61 novel miRNAs from 46 miRNA families were identified. In the total of 154 (93 known and 61 novel and newly annotated) miRNAs that were detected, 145, 149, 143 and 144 were expressed in the NC, CA, CCA, and CS conditions, respectively (Figure 1A).

Severe and moderate stresses had comparable impact on transcriptome and microRNAome

Among the three stress treatments, chilling shock (CS) perturbed the transcriptome the most as it triggered the largest number of genes to express and had the largest number of differentially expressed (DE) genes, detected by two stringent criteria (see Methods), with respect to the normal condition (NC) (Additional file 1: Table S1 and Figure 1A). The mild chilling stress, chilling acclimation (CA), came in the second.

Although CA and CS reached very different temperature – the former at 14°C whereas the latter at 4°C – the amount and extent of transcriptome variations that they had induced seemed to be comparable in terms of both expressed and differentially expressed genes. Furthermore, there was only a strikingly small difference between the perturbed transcriptomes they induced. Among the 2,855 and 3,297 DE genes of CA and CS in reference to NC, respectively, 2,175 (1,134 plus 1,041 in Figure 1B, 76.2% in CA and 66.0% in CS) were common. More specifically, among these DE genes, 62.72% (1,134 out of 1,808) were down-regulated and 47.99% (1,041 out of 2,169) were up-regulated along the same direction in CA vs NC and CS vs NC. Nonetheless, none of these DE genes were regulated in the opposite directions under the mild and severe chilling treatments (Figure 1B). In concordance, only 325 genes were DE between the CA and CS treatments (Figure 1A), further reflected by the similar expression patterns of the DE genes in CA and CS (Additional file 1: Figure S4). Furthermore, the microRNAome variations caused by CA and CS, reflected by the DE miRNAs, were in concordance with that of transcriptome variations. A large portion of DE miRNAs of CA and CS were in common; in particular, 34.0% (17 out of 50) up-regulated miRNAs were shared by CA and CS, and 31.5% (18 out of 57) down-regulated miRNAs were common to CA and CS. It is important to note that neither DE mRNA genes nor DE miRNAs reversed their expression from up-regulation to down-regulation or vice versa between the CS and CA treatments (Figure 1B).

These results on transcriptome variation were confirmed by a functional enrichment analysis; the two chilling treatments perturbed biological processes of similar functions (Additional file 2: Table S3A). Most of the enriched biological processes of the DE genes of CS and CA with respect to NC were in common, which included translation (GO:0006412, FDR < 1.16x10−06), photosynthesis (GO:0015979, FDR < 0.0162), L-serine metabolic process (GO:0006563, FDR < 0.0517) and various other metabolic processes. It is not surprising to observe that the two sets of DE genes share common enriched biological processes, because more than half of these genes were in common. Indeed, translation (FDR < 1.97x10−08), photosynthesis (FDR < 0.0244), and L-serine metabolic process (FDR < 0.0244) are also enriched in these 2,175 common DE genes (Additional file 2: Table S3B). Nevertheless, the stress-specific DE genes, i.e., DE genes that were specific to CA vs NC (680 out of 2,855) and that specific to CS vs NC (1,122 out of 3,297), were also enriched with the same biological processes as the commonly shared DE genes. For example, translation (FDR < 3.84x10−05) was enriched in the DE genes specific to CA vs NC and photosynthesis (FDR < 3.72x10−09) was enriched in the DE genes specific to CS vs NC (Additional file 2: Table S3B). In short, these results suggested that CA and CS perturbed similar biological processes.

Chilling after chilling acclimation reversed the expression of a large portion of DE genes

In stark contrast, chilling stress after chilling acclimation (CCA) altered the transcriptome the least among the three stress treatments despite that it ultimately reached 4°C as CS did. Surprisingly, the number of DE genes of CCA with respect to NC was approximately one third of that of CS or CA (Figure 1A and B). Further, 2,792 genes were DE between CCA and CS and 3,044 genes were DE between CCA and CA, which are eight times more than the 325 DE genes between CS and CA (Figure 1A).

In order to appreciate the role that chilling acclimation plays in stress response, we compared the DE genes and DE miRNAs of CA and CCA. A total of 2,855 and 1,082 genes were DE in CA and CCA with respect to NC, respectively. This more than 2.6 fold difference between the two sets of DE genes alluded to a genome-wide transcriptome alteration after stress acclimation in that harsh stress after chilling acclimation reversed the transcriptome changes triggered by the moderate chilling stress. Indeed, the expressions of about two thirds of the up- or down-regulated genes in CA with respect to NC were, respectively, reversed to down- or up-regulated in CCA with respect to CA (Figure 1C). Among the 1,507 up-regulated genes after the initial temperature decrease from the normal condition to CA at 14°C, 1,160 (77.0%) changed to down-regulated when going from CA to CCA at 4°C; likewise, among the 1,348 down-regulated genes after going from NC to CA, 641 (47.6%) genes reversed to up-regulated after going from CA to CCA (Figure 1C). In contrast, even though the temperature kept decreasing from CA to CCA, only 1 down-regulated gene was further down-regulated and 3 up-regulated genes were further up-regulated (Figure 1C).

Interestingly, miRNAs might be responsible for the reversion of some of the DE genes to their original expression levels. Specifically, 14 and 16 (87.5%) miRNAs that were down- and up-regulated from NC to CA reversed, respectively, to up- and down-regulation going from CA to CCA; furthermore, 37.0% (30 out of 81) of the DE miRNAs reversed their expression directions (Figure 1C). Importantly, these 30 DE miRNAs targeted 1,198 mRNA genes, among which 48 were DE and reversed their expression directions going from NC to CA and to CCA (Figures 1C and Additional file 1: Figure S5 and Additional file 3: Table S4). This observation suggested that 48 mRNAs with the reversed expression levels might be negatively regulated by the 30 DE miRNAs, which also have reversed expression patterns. The 1,801 (641 plus 1,160) reversely regulated DE genes in the treatments from NC to CA and then from CA to CCA have enriched biological processes such as translation (FDR < 3.69x10-15), superoxide metabolic process (FDR < 0.0052), and mismatch repair (FDR < 0.061) (Additional file 2: Table S3C). We experimentally tested 4 of the reversely regulated genes related to translation (ribosomal protein L11, ubiquitin 6, ribosomal protein L31e, and zinc-binding ribosomal protein) using qRT-PCR (see Methods). The reversed expression patterns of these 4 translation-related genes are consistent with the RNA-seq data and can be confirmed by qRT-PCR assay both under CA/NC and under CCA/CA in the corresponding regulation direction regardless of the magnitude (Figure 2A and B).
Figure 2

Experiment validation of reversely expressed genes from CA/NC to CCA/CA and differentially expressed genes specific to CCA or CS. The reversed expression patterns of four translation-related genes are firstly shown in (A) RNA-Seq data and are then confirmed by (B) qRT-PCR. The expression patterns of genes, which are associated with nutrient reservoir specific to CCA and which are associated with viral reproduction specific to CS are shown in (C) RNA-Seq data and are then confirmed by (D) qRT-PCR. We remove the prefix of gene names for simplicity; e.g. “018150 m” is short for “Cassava4.1_018150m”. The expression level of each gene is averaged among three replicates and is normalized by Cassava actin gene in the qRT-PCR assay. The y-axis indicates the relative expression level of a gene in a given condition with respect to the normal control (NC) in log scale. (Gene list, 018150 m: ubiquitin 6, 017802 m: ribosomal protein L11 family protein; 019383 m: ribosomal protein L31e family protein; 020116 m: zinc-binding ribosomal protein family protein; 015731 m and 029709 m: RmlC-like cupins superfamily protein; 000174 m: unknown protein; 003690 m: ROP interactive partner 3).

Chilling acclimation prepared the plant to fend off adverse effects of further stress

The significant transcriptome and microRNAome changes, which reversed most DE genes and a substantial number of miRNAs in chilling after chilling acclimation, were in concordance with the mild symptoms of chilling stress of CCA in comparison with the symptom of CS (Additional file 1: Figure S1). This suggested that as the temperature further decreased to 4°C after the initial moderate stress, the plant was able to better adapt to further harsh stress and effectively recover some of the perturbation to biological processes or pathways that had been altered. In other words, stress acclimation (i.e., CA) helped the plant develop a kind of immunity against adverse impact of chilliness at 4°C. Moreover, miRNAs played a role in this process by regulating the expression of some mRNA genes.

A direct comparison between the transcriptomes of and the biological processes affected by CCA and CS further confirmed our observation. Even though both the CCA and CS treatments reached the same temperature of 4°C, plants that had experienced chilling acclimation at CCA exhibited drastically different transcriptome from plants that had subjected to CS. Similar to the difference between CA and CCA, there was also a huge, more than 3.0 fold, disparity between the number of DE genes of CCA (i.e.,1,082 genes) and that of CS (i.e., 3,297 genes) with respect to NC (Figure 1A). Only 446 of these DE genes were common in both chilling treatments, while 57.8% (625 out of 1,082) of them were specifically DE in CCA vs NC and 86.1% (2,840 out of 3,297) of them were specifically DE in CS vs NC (Figure 3A). This indicated that the degree of transcriptome perturbation due to CCA and CS were substantially different. The functions of the two sets of DE genes for CCA and CS helped reveal the different biological functions of the common and condition-specific genes (Figure 3B). As expected, many biological processes, such as metabolic process and responsive to stimulus, were perturbed under both chilling conditions since 446 DE genes were common to both CCA and CS (Additional file 2: Table S3A and Figure 3A). However, 8 biological processes (i.e. anatomical structure formation, cellular component biogenesis and organization, death, developmental process, multi-organismal process, multicellular organismal process and reproductive process) were not present in the common genes but exclusively appeared in the CCA/CS specific genes. These genes function in 4 different cellular components (envelope, extracellular region, extracellular region part, and membrane-enclosed lumen). The 8 biological processes, which were either associated with the anabolic, growth-promoting state or the catabolic, growth-suppressing state, play a key role in leading Cassava to distinct responses under the CCA and CS treatments. The cellular components of these CCA/CS specific DE genes are closely related to membrane rigidification caused by chilling stress.
Figure 3

Common and distinct functions of DE genes in CCA and CS with respect to NC. (A) Relationship between the up- and down-regulated genes and that of the DE mRNAs in CS and CCA with respect to NC. The two numbers within a region are the numbers of DE genes (the first line) and miRNAs (the second line). (B) Biological processes, molecular functions and cellular components that were affected by the DE genes of CCA and CS. The function of nutrient reservoir function is exclusively associated with some of the DE genes of CCA, whereas the process of viral reproduction is exclusively associated with some of the DE genes of CS.

In addition, one molecular function was specific to DE genes in CCA and one biological process was specific to CS (Figure 3B). The term nutrient reservoir was only present in DE genes in CCA condition, suggesting their possible protective activities in plants undergone chilling acclimation. Two RmlC-like cupins superfamily protein-coding genes (015731 m and 029709 m) were associated with nutrient reservoir. They were expressed at normal expression levels under the CS condition, but in contrast, they were overexpressed 4 folds under the CCA condition with respect to NC (Figure 2C and D). In contrary, the CS treatment triggered the process of viral reproduction, which did not appear in the CCA treatment. From the RNA-seq data, eight DE genes were associated with viral reproduction, and six out of these eight DE genes were overexpressed in both CA and CS, but expressed normally in CCA. The other two genes were down-regulated by 5–6 folds in CA and CS, but were again expressed normally in CCA. One up-regulated unannotated protein (000174 m) and one down-regulated ROP interactive partner 3 (003690 m) were validated by qRT-PCR methods (Figure 2C and D).

MicroRNAs contributed to chilling response and stress acclimation

As a major post-transcriptional gene regulator, miRNAs mediate the expression of their target genes in adaptation to environmental stress. To appreciate miRNA functions and gain insight into the complex regulatory networks in chilling response in Cassava, we exploited and combined the large collection of mRNA and small-RNA profiling data for an integrated transcriptome and microRNAome analysis. Based on the small-RNA profiling data, we identified 61 novel miRNAs in Cassava and 121 DE miRNAs in one of the six comparisons we considered (Figure 1A).

We identified anti-correlated pairs of DE miRNAs and DE mRNA target genes in 6 comparisons across the stress conditions (Additional file 4: Table S5A). The DE miRNAs played a role in regulating the transcriptome responses to chilling stresses. More than 30 miRNAs regulated at least two potential mRNA targets in each comparison. The regulatory effect of these DE miRNAs was most profound under CS with respect to NC and CCA conditions (Figure 4). We further investigated the miRNA-regulated DE genes that were associated with the significantly enriched pathways in the pairwise comparisons among six conditions (Additional file 4: Table S5B). These DE genes were regulated by miRNAs under at least one of the chilling stresses with respect to NC, while few of them were DE across any of two stress conditions (Additional file 4: Table S5B). Biosynthetic process (29.0%), cellular protein modification process (18.2%), response to stress (11.1%) and metabolic process (10.4%) were the top 4 biological processes that had the largest ratios of genes being regulated by miRNAs to the total number of genes on a given enriched pathway (Additional file 4: Table S5B). As expected, one miRNA can regulate one or several genes, and one mRNA gene may be targeted by multiple miRNAs. For example, 4 genes (005409 m, 006360 m, 006048 m, 005437 m and 005421 m) targeted by miR399 were enriched in the same biosynthetic process in CA vs. NC, at the same time, one gene (012052 m) enriched in oxidation-reduction process and 3 genes (033858 m, 014142 m and 000730 m) enriched in metabolic process under CCA vs. NC were all potentially targeted by miR396a/b/c/d, and novel-3 could target one gene (013577 m) enriched in metabolic process under both CCA vs. NC and CS vs. NC (Additional file 4: Table S5B). Surprisingly, one translation-related gene regulated by miR172 was differentially expressed in four comparisons: CA vs. NC, CS vs. NC, CCA vs. CA and CS vs. CCA, indicating that the interaction between miR172 and 018488 m played regulatory roles in CA and CS, but not in NC and CCA. In short, miRNAs are one of the key regulating factors during the processes of low temperature adaptation in Cassava.
Figure 4

The DE genes targeted by DE miRNAs among six comparisons. The number of miRNAs, e.g., listed under “CA/NC”, refers to the number of DE miRNAs in the CA condition with respect to NC. The number of miRNA:mRNA pairs, e.g., listed under “CA/NC”, is the number of mRNA targets that are not only DE but also anti-correlated with their targeting miRNAs in the CA condition with respect to NC.

The effect of miRNAs on their targets is reflected by the anti-correlated expression patterns between miRNAs and their mRNA targets because the major function of plant miRNAs is mRNA cleavage. We further experimentally tested 17 pairs of anti-correlated miRNAs and target mRNAs, as initially detected by the sequencing data (Figure 5). We first examined cleavage cites of miRNAs on their target genes by 5′RACE (see Methods). The cleavage sites of 13 (76.5%) of the 17 pairs were validated and the cleavage sites were within the regions of miRNA binding sites (Table 1). As expected, many of these anti-correlated miRNAs and mRNAs were related to stress responses, such as novel16-POS, where POS is associated with scavenging hydrogen peroxide, and miR398-EC, where EC is supposed to maintain the membrane potential via electron carrier. Therefore, these identified anti-correlated pairs of miRNA and mRNAs from the profiling data should be valuable candidates and subject to further investigation. NF-YA family, targeted by miR169, has been recently found as an adaptive response to adverse environmental conditions [17]. Note that the cleavage site of NFYA10 by miR169 was indeed detected in our 5’RACE assay, while the anti-correlation between miR169 and NF-YA10 was not detected due to low reads number of miR169 in the RNA-seq data (data not published).
Figure 5

Regulatory networks showing the relationship between differentially expressed miRNAs and their anti-correlated, differentially expressed target mRNAs. The yellow diamonds represent miRNAs and the blue circles denote target mRNAs. An edge between a pair of miRNA and mRNA indicates their anti-correlation relationship across two conditions, which is color coded. (A) the upper panel is referred to one of three chilling treatment (CA, CCA and CS) with respect to NC, which assigned to red, green and blue line, respectively; (B) the bottom panel is referred to in-between comparison of three chilling treatments, including CCA vs. CA (red), CS vs. CA (green) and CS vs. CCA (blue).

Table 1

The cleavage site of anti-correlated miRNA:mRNA pairs detected by 5′RACE assay

miR-Target

ID

Cleavage sites*

Anti-correlation§

Annotation

NC

CA + CCA

CS

CA/NC

CCA/NC

CS/NC

CCA/CA

CS/CA

CS/CCA

novel7a-DHDSP2

010211 m

4,±

  

+

     

dihydrodipicolinate synthase

Novel53-RPL27

016749 m

 

21,±

±

+

 

-

   

ribosomal protein large subunit 27

miR171hi-HAM3(1)

002034 m

 

10(5),7(2),-

10,7,3,-

+

+

    

GRAS family transcription factor

miR395-ATRX

000210 m

 

9(11)

9,-

+

 

+

   

P-loop containing nucleoside triphosphate

hydrolases superfamily protein

novel16-POS

033835 m

 

±

5,0,±

+

     

Peroxidase superfamily protein

Novel42-LMDH

007889 m

 

-

11,±

   

-

  

lactate/malate dehydrogenase family protein

miR394-FY

002500 m

 

4,3,1(2),-

4,-

     

+

Transducin/WD40 repeat-like superfamily protein

miR171hi-HAM3(2)

034057 m

 

10(3),7,-

10(5),7(3),-

     

+

GRAS family transcription factor

miR319-TCP4

004517 m

  

8,-

  

+

  

+

TCP family transcription factor 4

novel30-CP

014450 m

12,10

  

+

+

+

 

-

 

Cyclase family protein

miR164-NAC

029889 m

 

10(7)

10(9)

 

+

+

 

+

 

NAC domain containing protein 1

Novel53-SWEET17

014996 m

 

22,17(2),-

±

 

+

 

+

-

 

Nodulin MtN3 family protein

miR398-EC

011902 m

 

±

10(2),±

 

-

   

+

copper ion binding/electron carrier

miR319-GSTU8

015766 m

     

-

-

-

-

glutathione S-transferase TAU 8

miR172-TUA

007502 m

-

  

-

     

tubulin alpha-2 chain

miR172-LMP

002489 m

 

+

-

-

-

    

LMBR1-like membrane protein

miR394-IIL1

005847 m

  

-

-

  

-

-

 

isopropyl malate isomerase large subunit 1

The cleavage sites of 13 (out of 17, 76.5%) miRNA:mRNA pairs were detected within miRNA binding regions. Five (the first five from the top), three (the subsequent three) and five pairs (the next subsequent five) of them were anti-correlated in the three chilling stress conditions with respect to the NC condition, in the other three comparisons and both two comparisons according to RNA-seq data, respectively. The cleavage sites of the remaining four pairs (the last four) failed to be detected in the basepairing regions of corresponding target genes. The “*” column lists number of miRNAs that are complementary with target cleavage sites in binding regions; the “+” and “-” symbols refer to cleavage sites located in upstream and downstream of miRNA:target binding regions, respectively. In the “§” column, “+” means that there is an anti-correlation between a miRNA and its target gene based on sequencing data, and the cleavage site of miRNA:target is located in the basepairing region in 5′RACE assay, and “-” means that there is an anti-correlation in sequencing data, whereas the cleavage site of miRNA:target is not within the basepairing region but either upstream or downstream to the basepairing region.

Discussion

This is the first genome-wide, systematic study of chilling stress and acclimation in Euphorbia; we integrated genome-scale, sequencing-based profiling data on protein-coding genes and miRNAs on a chilling sensitive cultivar of Cassava, the most important crop in the tropical regions. The study identified critical genes and miRNAs are responsive to chilling stresses and gained insights into the biological processes underlying chilling acclimation. We expect our findings to have a potential impact on plant biology and agriculture. In particular, exposing tropical or subtropical crops to sub-optimal temperatures for stress acclimation can enhance their chilling resistance.

Comparable effects of severe and moderate stress

Based on a comparative transcriptome analyses of potato (S. tuberosum) and A. thaliana, Carvallo et al. [18] reported conserved gene expression patterns between the two species during chilling acclimation despite their independent evolution over 100 millions of years. In a large scale evolutionary studies across various species, Preston and Sandve [19] further documented that despite independent development separated by evolution over hundred millions of years, some pathways and mechanisms involved in chilling acclimation are similar between bryophytes, monocots, and eudicots. As for different cultivars of the same species with contrasting responses to chilling stress, their comprehensive comparative studies reveal relatively similar initial transcriptional responses and more diverse downstream molecular changes. While functional categories for genes responding to late phase chilling stress are diverse, ranging from functional adaptation to continuous stress, the early responses have been found to be related to transcription regulation and signal transduction [20]. Similarly, Usadel et al. [21] reported that Arabidopsis has qualitatively similar responses to 17, 14, 12, 10 and 8°C, regardless the degrees of temperature change. We also confirmed this observation in our CA and CS experiments with similar disruptions to the transcriptome. Moreover, we reported here that the chilling-sensitive Cassava could trigger similar transcriptional responses when the Cassava plant was subject to chilling stresses close to freezing temperature (4°C).

Chilling acclimation helps plants better withstand severe chilling stress

Our results revealed that the largest numbers of expressed and DE genes appeared in the chilling shock (CS) stressed plants. This suggested that a more severe transcriptome perturbation in Cassava plants that had experienced dramatic temperature decreasing (CS) compared to plants that were treated by mild chilling stress (CA) and were stressed at same temperature after acclimation with mild chilling stress (CCA). Surprisingly, more than 65% of the up- or down-regulated genes in CA vs. NC were respectively swapped to down- or up-regulation in CCA vs. CA, indicating a remarkable reversing effect on gene expression and a dramatic rewiring of gene regulatory networks in plants after experiencing chilling acclimation. Such antagonistic regulation patterns and regulatory network changes reflected plants’ profound adaptation to stress through the process of chilling acclimation.

In chilling sensitive plants, a major adjustment is the physical transition of cell membrane from a flexible liquid-crystalline to a solid gel phase, resulting in increased permeability that leads to cellular leakiness and ion imbalance. In Cassava, significant accumulation of hydrogen peroxide (H2O2) have been observed after 4 h chilling exposure; at the same time, the enzyme activities and up-regulation of four ROS-scavenging transcripts, including catalase (CAT), superoxide dismutase (SOD), and glutathione transferase (GST), have been observed and verified [22]. The enzymatic antioxidants, jointly with low molecular-mass antioxidants that were embedded in membrane and polyunsaturated fatty acids, respond to the oxidative stress. The balance of oxidative capacities and scavenging activities of antioxidants is broken as a consequence of abnormal metabolism and injured cells accumulating toxic metabolites and active oxygen species. In our study, oxidoreductase activity was significantly enriched in all of the chilling stresses compared to the control condition, indicating serious oxidative stresses. Furthermore, the DE genes that were specific to CS or CCA significantly varied in the anatomical structure formation and cellular component. In addition to the fact that the chilling injury index (Malondialdehyde content in Additional file 1: Figure S2C) increased most significantly in the CS condition, it is viable to infer that transmembrane damage and membrane permeability have a significant variation between the CCA and CS conditions. In addition, the processes of cell death and development were also significantly distinct between CCA and CS (Figure 3B), indicating that membrane and metabolism adjustment played an important part in the chilling acclimation in Cassava.

We further experimentally validated some of the genes that reversed their expression patterns under further chilling stress after stress acclimation. Ribosomal protein L11 and ubiquitin 6 were up-regulated from the NC to CA condition and subsequently down-regulated from the CA to CCA condition (Figure 2A and B). Both genes related to ubiquitination can reduce protein synthesis for normal growth and provide diversity of small protein nutrient availability for the later chilling adaption [23]. In contrary, ribosomal protein L31e superfamily and zinc-binding ribosomal proteins, which are related to the growth of meristem tissue, were first down-regulated and then up-regulated from the NC to CA and subsequently to the CCA condition. The expression pattern of these genes can help determine whether a cell is in an anabolic, growth-promoting state or a catabolic, growth-suppressing state [24]. Two previous studies also found that low temperature can induce a large number of genes involved in translation, protein synthesis or nucleosome assembly when transferring from normal condition to sub-optimal temperature in Arabidopsis[21],[25],[26]. Ribosome biogenesis is a key process for fundamental translation processes. Perturbations to the dosage of the ribosomal protein subunits regulate overall protein synthesis related biological process [24]. The hydrophilic residues of ribosomal protein L11 can interplay with p53-MDM2 function complex [27],[28] to stabilize and activate ribosomal protein-Mdm2-p53 signaling pathway to response DNA damage and ribosomal stress [23],[29]. Those genes altering the accomplishment of normal translation initiation [30],[31], elongation [32], termination [33], and probably ribosome-recycling [34], are supposed to reduce the rate of protein synthesis. Subsequently, Carroll [35] and Gerashchenko [36] reported that translational proteins have exquisitely sensitive and responsive to environmental fluctuations. Furthermore, Ferreyra [37] argued that this widespread and unequivalent translational component reprogramming can have a “turnover” effect on related mRNAs while inducing the translation of specific mRNAs in adaptation to environmental stress.

In addition to translation reprograming, plant resistance to low temperature depends on responsive speed and scope of transcripts and metabolites involved in cryoprotection and stress responses. We identified two RmlC-like cupins superfamily genes specifically DE in the CCA condition with respect to NC. The RmlC-like cupins superfamily proteins have known functions in biology processes of small nutrient molecule, such as sulfur amino acid, polysaccharide metabolic process, amine and carbohydrate biosynthetic process. These metabolites have known oxidoreductase and disomerase activities [38] and cation and hormone binding capacities [39]. They potentially help stressed plants to maintain normal osmotic potential and to defend water loss.

ROP interactive partner 3 (003690 m) was shown to significantly express in the CS condition, ROP functions as molecular switches in auxin signaling at the cell membrane to regulate many cellular and developmental processes and environmental responses [40]. The balance between the activation and inactivation of ROP helps maintain a non-toxic level of H2O2 to avoid lethal effects [41]. Yun et al. [7] report that early response triggered by oxidative signals is critical for prolonged survival under sub-optimal temperature and that many of the oxidative-mediated changes in gene expression occur within 24 hours of stress. Their transcriptome variations have striking similarities to the transcriptome variations observed in plants’ responses to disease, pathogen infection and wounding, where DE genes are involved in redox regulation. Significant down-regulation of ROPs may contribute to the imbalance of redox regulation, resulting in effects similar to viral invasion, which may be critical for enduring transient exposure to mild chilling stress. Since RmlC-like cupins can recruit genes with functions of nutrient reservoir for adapting to a lower temperature under the CCA condition and ROP genes lead plants to a catabolic, growth-suppressing state under the CS condition, the differences in the regulation of these positive and negative affecting genes may potentially cause the difference in chilling tolerance between the CCA and CS conditions.

Above-zero chilling resistance after chilling acclimation

By a conventional standard, chilling acclimation is the process by which plants acquire freezing tolerance through exposure to suboptimal, low, nonfreezing temperatures. Plant species vary widely in their abilities of chilling acclimation. The chilling-sensitive species, such as tomato and rice, may be unable to be acclimated to acquire freezing tolerance. However, chilling resistance should not be limited to sub-zero temperatures, especially for tropical and sub-tropical crops. Although being intolerant to freezing, they should be able to tolerate a much lower chilling temperature after chilling acclimation. This increased chilling resistance can allow tropical crops to withstand more severe chilling stresses than before. Unfortunately, little has been done in the past to address the impact of above-zero chilling resistance.

CBF (C-repeat/dehydration-responsive element binding transcription factor) genes are known to be ubiquitous and have conserved structures and functions in a wide range of plant species, including tropical and temperate plants [4],[42],[43] as well as other plants [18],[44]. All four CBF genes in Arabidopsis, in which CBF1, 2, 3 are cold inducible and CBF4 is drought inducible, are arranged in tandem in the Arabidopsis reference genome [45],[46]. In Cassava, 70 CBF homologues can be identified through a sequence search with an e-value less than e−20. It is perplexing that despite having such a large number of CBF-like genes, cassava is unable to withstand freezing temperatures. Nevertheless, Cassava can endure chilling temperatures to some extent, particularly after chilling acclimation as shown in our study. It is interesting to observe that none of the 154 miRNAs (93 know and 61 novel), detected by our deep sequencing data, targeted any of CBF genes. Furthermore, these CBF genes, as well as the homologues of CBF constitutive activator ICE [47],[48], had insignificant expression variation in chilling treatments partly due to their low expression levels. Those results are consistent with the results from an early Cassava chilling study [22], suggesting that the CBF regulon in Cassava may be different from its counterparts in Arabidopsis. On the other hand, not being chilling responsive at the transcriptional level did not mean that the CBFs were not responsive to chilling stress because CBFs may also be subjected to ubiquitination and sumoylation regulation at the translational and posttranslational levels, as shown in Arabidopsis[49]–[52].

The overexpression of SsDREB1 in Arabidopsis could result in chilling-regulated gene expression and freezing tolerance [53]. Conversely, constitutive overexpression of either LeCBF1 or AtCBF3 in transgenic tomato plants does not provide a detectable freezing tolerance [54]. The regulatory functions of the CBF regulon may be similar, but the whole CBF chilling-response pathway could be different between the tropical and temperate plants and the other plants. One clue to understand this difference may come from the upstream regulation components. There are probably different upstream motifs of the CBF regulon and some other transacting components in the regulatory elements of the CBF genes. Another possible reason can be attributed to the genes downstream of the CBF regulon, which lack functional CRT/DRE elements in the promoters of genes that are required in order to affect chilling acclimation response or freezing tolerance. As Zhang et al. [54] suggest that differences among plant species with respect to the DREB/CBF regulon could be due to limited DRE/CRT-like elements in the genomes of chilling sensitive species, which may be the case for rice. The upstream regulation elements, such as cis-regulatory motifs and trans-factors, should be further examined in the DE gene identified in our study.

McKhann et al. [55] analyze the CBF genes and their promoters diversity in the Arabidopsis versailles core collection containing 48 accessions that maximize the naturally occurring genetic diversity, and found polymorphisms in the CBF genes along with differences in CBF and COR gene expression. Although there may be more CBF and COR gene expressed in tolerant accessions, there is no simple correlation between CBF regulon and chilling tolerance. By selecting another closely related Solanum species that differs in its chilling acclimation ability, Thomashow and colleagues [2] find that both tested species altered gene expression in response to low temperature to similar degrees with similar kinetics and that both plants have CBF regulons composed of hundreds of genes [18]. These two studies, along with our finding that there was few miRNA or mRNA related to the CBFs’ response to the three designed chilling stress treatments, suggest that a complex network of genes are involved in the chilling acclimation [56]–[59], and that the CBF genes alone cannot explain all phenotypic differences. It suggests that more investigation is required in addition to the CBF regulon, e.g., metabolite networks that are related to low temperature tolerance.

Conclusions

The results of the transcriptome and microRNAome profiling and the integrated analysis provided three important observations. First, substantial chilling stresses, regardless whether they were severe (4°C) or moderate (14°C), could lead to comparable transcriptome and microRNAome variations and cause perturbation to biological processes of the same or similar functions. Second, dramatic temperature decreases and further temperature decreases after chilling acclimation could result in drastically distinct transcriptome and microRNAome variations even though both stresses reached at the same temperature of 4°C. Third and most importantly, chilling-sensitive Cassava could be accustomed to, through stress acclimation, severe chilling stresses close to freezing temperature (4°C) at least for a short period of time (5 days in the current study). This might be achieved by reversing the expressions of a large number of affected genes, which may potentially alter the regulatory networks, specifically by inducing or overexpressing genes that help preserve nutrients that are critical for stress tolerance.

Methods

Plant materials and stress treatments

Stem segments with three nodes of Cassava (Manihot esculenta Crantz) cultivars were extracted from 8-month-old plants, and inclined in 3-L pots filled with barren red soil: vermiculite (1:1, v/v), fertilized with Hoagland’s solution [60], to propagate and generate well-balanced seedlings. The solution was renewed with 300 ml quarter-strength solution once a week. After 2 months of planting, the uniform seedlings were subjected to chilling stress treatments. All plants were field grown in Haikou, Hainan, China, during April and June of natural conditions (11 h light, 13 h dark and 25°C during the day and 18°C at night).

Plants of Cassava cultivar SC124 were transferred to normal 24°C illumination incubator for 2 days to set a homogenous starting point, and were treated with three types of chilling stress. 1) Gradual chilling acclimation (CA in Additional file 1: Figure S1): temperature was decreased from 24°C to 14°C at the rate of −2°C/h to exert a moderate chilling stress. Temperature was then held constant at 14°C for five days to accommodate chilling acclimation. 2) Chilling stress after chilling acclimation (CCA): after 5 days of chilling acclimation and growth under 14°C, plants were watered once with Hoagland’s solution, transferred further to 4°C at the gradient of −2°C/h, and cultivated at constant 4°C for another 5 days. 3) Chilling shock (CS): normal plants under 24°C grew under 24°C for 5 days in parallel to the CA treatment, and were then subjected to dramatic temperature decline to 4°C at the rate of −4°C/h to ensure the temperature to reach 4°C at the same time as the CCA treatment finished so that these plants had the same age. The effects of the chilling stresses are shown in Additional file 1: Figure S1. In the CA treatment, total RNA was collected at 6 h, 24 h and 5d after temperature reaching 14°C; in both the CCA and CS treatments, total RNA was collected at 6 h, 24 h and 5d after temperature reaching 4°C. These time points were chosen to account for initial response, secondary response, and functional adaption to chilling stresses. In parallel, plants grown under the normal condition (NC) of 24°C were watered once with Hoagland’s solution every 5 days, and total RNA was first collected at 6 h, 24 h and 5d along with total RNA for the CA plants being collected and then at 6 h, 24 h and 5d when total RNA for the CS and CCA plants were collected so as to match samples of chilling treatments.

RNA isolation, RNA library preparation and NextGen deep sequencing

Three organs/tissues – folded leaf, fully expanded leaf and roots of Cassava cultivar SC124 –were harvested at 6 h, 24 h and 5d for the three chilling treatments of CA, CCA and CS as well as the normal condition as described above. Total RNA was isolated separately for each of the three organs/tissues. The three RNA samples of equal amount for each time point of a treatment were then pooled together for expression profiling. As a result, four mRNA libraries and four small-RNA libraries, corresponding to the conditions of CA, CCA, CS and NC, were constructed.

The four mRNA libraries were sequenced by RNA-seq by Illumina GAII following Illumina RNA-seq protocol. Briefly, total RNAs were isolated, purified and reversely transcribed, the resulting cDNA products were subsequently digested with NlaIII and the 3′-cDNA fragments captured with the oligo(dT) beads, and then ligated to the Illumina GEX NlaIII Adapter 1. The junction of Illumina adapter 1 and CATG site contained the recognition site of MmeI, cutting 17 bp downstream of the recognition site (CATG) to produce tags. After removing 3′fragments with magnetic beads precipitation and MmeI digestion, an Illumina GEX adapter 2 was introduced at the end of tags. The resulting adapter-ligated cDNA tags were subjected to 15 cycles of linear PCR amplification, and then purified and sequenced with the method of sequencing by synthesis (SBS) using Illumina GAII.

The six small-RNA libraries (four for Cassava and two for castor bean) were subjected to small-RNA deep sequencing by Illumina Genome Analyzer (GAIIx). Briefly, total RNA was isolated using RNA plant Reagent kit (TIANGEN, Beijing, China). Small RNAs were enriched by poly-ethylene glycol precipitation, separated on 15% denaturing PAGE, and visualized by SYBR-gold staining. Small RNAs of 16- to 28-nt were gel-purified. Purified small RNAs were ligated to a 5′ adaptor and a 3′ adaptor sequentially, reverse transcription amplified, and sequenced.

All mRNA and small-RNA sequencing data have been deposited into NCBI/GEO under accession # GSE52178.

mRNA differential expression assay

The RNA plant Reagent kit (TIANGEN, Beijing, China) was used for total RNA isolation. The quantity and quality of extracted total RNAs were detected by 1% agarose gel electrophoresis and spectrophotometer. To assess the expression of an mRNA, the first-strand cDNA products were prepared using Fermentas reverse transcriptase kit (K1611). Real-time PCR was performed following the standard SYBR Premix Ex Taq™ kit (TaKaRa) protocol. The reactions were incubated in 0.1 ml tubes of Rotor-gene 6000 machine. The procedure ended by melt-curve. A negative control (no template) was included for each primer combination.

Gene actin served as a reference for mRNA in each sample, genes were amplified in parallel for 3 replicates. The relative concentration was calculated as 2 powered -CT, where CT = (CT sample - CT control), CT = CT (target)-CT (reference) in each sample. If the CT value was greater than that of no template control (NTC), the mRNA was considered not expressed. The gene specific primers used in the target validation are provided in Additional file 1: Table S6.

Experimental miRNA target validation

RNA ligase-mediated rapid amplification of 5′cDNA ends (RLM-RACE) GeneRacer Kit (Invitrogen, USA) was used to validate miRNA-guided mRNA cleavage, which differed with traditional 5′RACE of full-length cDNA by omitting the 5′ phosphates of truncated mRNA removal and the 5′ cap structure of full-length mRNA removal treatments. Briefly, total RNA was extracted with RNAplant regent (TIANGEN, DP407-01), and PolyA RNA was isolated using polyAtract mRNA isolation system III (Promega, USA) to eliminate contaminated non-mRNA. Ligation with a 5′ RNA adapter and a reverse transcription were performed. The resulting cDNA was used as a template for PCR amplification. Two ~100 bp spaced gene specific reverse primers (GSP1 and GSP2) for each target, designed based on the downstream sequence of the miRNA:target binding site at the target gene sequence. Combining with two GeneRacer 5′ forward primers (included in GeneRacer kit) to specifically nest amplify the 3′ cleavage product of the target mRNA. The amplified PCR products were gel purified, cloned and sequenced (Sangon, China). Gene specific primers that we used are provided in Additional file 5: Table S7.

Processing of RNA-seq and small RNA-seq data

Raw RNA-Seq reads with low quality were firstly discarded. Sequencing adaptors were trimmed using an in-house method that recursively searches for the longest substring of the adaptor appearing within a sequence read. If a raw sequence read did not have a substring of the adaptor longer than 6 nt, it was considered as having no adaptor. The adaptor-trimmed sequences with no ambiguous reads, which were referred to as qualified reads, were then mapped to the cDNA sequences and reference genome of Cassava using Bowtie 0.12.7 [61] allowing no more than one mismatch. For small RNA sequencing data, the raw sequence reads with no 3′ sequencing adaptor, of low quality, or shorter than 17 nt were removed. The adaptor trimming for the small RNA-seq reads was done in the same way as for the mRNA reads. The qualified small RNA reads were mapped to the Cassava genome using Bowtie allowing no more than one mismatch.

Identification of expressed and differentially expressed mRNA genes

The number of reads mapped onto each mRNA transcript was recorded as a raw read count. A gene was considered as expressed if its CPM (Count Per Million mapped reads) was not less than 10. For genes that had CPM less than 10, we considered the aligned reads as noise and set their expression level to 0. For expressed genes, we normalized the raw read counts using the upper-quartile normalization method [62]. Given two conditions to be compared, a gene was considered differentially expressed if any of the following two criteria was satisfied: (1) the change of normalized counts was no less than 4 folds if the gene was expressed in both conditions, or (2) the gene was not expressed (CPM < 10) under one condition and over-expressed (CPM > 40) under the other condition.

Identification of miRNA targets and anti-correlated miRNA-mRNA pairs

The TargetFinder program [63] was used for computational prediction of miRNA targets. TargetFinder uses a scoring scheme that charges a penalty of 1 for a mismatch or insertion/deletion and a penalty of 0.5 for a wobble base pairing, and it doubles these penalties for base-pairs in the seed region (the 2-7 nt from the 5′-end) of a miRNA. To reduce false positive targets, we set the threshold for target alignment score to 4 so that only target genes that have alignment score no greater than 4 were retained as putative targets. A pair of miRNA and its putative target was considered anti-correlated if both the miRNA and the target were DE and the target was up or down-regulated in the opposite direction with respect to the miRNA expression change.

GO function enrichment analysis

The Gene Ontology annotations of Cassava mRNAs were retrieved from Phytozome website v4.1 [64]. The statistical significance of GO annotation enrichment was measured by Fisher’s exact test. For each ontology category, the p-value was calculated as the probability under which we would observe at least k genes to have a given GO term, t, if we randomly selected m genes from the given M genes in the genome, where n genes were associated with term t. Given a term and a list of genes, k follows the hyper-geometric distribution. We also reported the False Discovery Rate (FDR) based on Benjamini-Hochberg multiple testing correction [65].

Supporting data

The raw sequencing and processed data from this publication have been deposited into GEO database (http://www.ncbi.nlm.nih.gov/geo/) with accession number GSE52178.

Additional files

Abbreviations

CA: 

Chilling acclimation

CCA: 

Chilling stress after chilling acclimation

CS: 

Chilling shock

NC: 

Normal condition

Declarations

Acknowledgements

The work was supported in part by Chinese National Basic Research and Development Program (973) (2010CB126602), Natural Science Foundation of China (31300999 and 31101193), National Nonprofit Institute Research Grant of CATAS-ITBB (ITBB140101), the municipal government of Wuhan, Hubei, China (2014070504020241 and the Talent Development Program), the Hainan Province Natural Science Foundation (313075), the Major Technology Project of Hainan (ZDZX2013023-1), and an internal research grant of Jianhan University, Wuhan, China, as well as by United States National Institutes of Health (R01GM100364) and United States National Science Foundation (DBI-0743797).

Authors’ Affiliations

(1)
The Institute of Tropical Bioscience and Biotechnology, Chinese Academy of Tropical Agricultural Sciences
(2)
Institute for Systems Biology, Jianghan University
(3)
Department of Computer Science and Engineering, Washington University in St. Louis
(4)
Department of Genetics, Washington University School of Medicine

References

  1. Levitt J: Responses of Plants to Environmental Stress. Chilling, Freezing and High Temperature Stresses. Academic, New York, 1980:497-607.Google Scholar
  2. Thomashow MF: Plant cold acclimation: freezing tolerance genes and regulatory mechanisms. Annu Rev Plant Physiol Plant Mol Biol. 1999, 50: 571-599.View ArticlePubMedGoogle Scholar
  3. Sanghera GS, Wani SH, Hussain W, Singh NB: Engineering cold stress tolerance in crop plants. Curr Genomics. 2011, 12: 30-43. 10.2174/138920211794520178.PubMed CentralView ArticlePubMedGoogle Scholar
  4. Le MQ, Engelsberger WR, Hincha DK: Natural genetic variation in acclimation capacity at sub-zero temperatures after cold acclimation at 4 degrees C in different Arabidopsis thaliana accessions. Cryobiology. 2008, 57: 104-112. 10.1016/j.cryobiol.2008.06.004.View ArticlePubMedGoogle Scholar
  5. Chinnusamy V, Zhu J-K, 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 CentralView ArticlePubMedGoogle Scholar
  6. Costello LR, Perry EJ, Matheny NP, Henry JM, Geisel PM: Abiotic Disorders of Landscape Plants-a Diagnostic Guide. University of California Division of Agriculture and Natural Resources; 2003:133–138. 3420.Google Scholar
  7. Yun K-Y, Park MR, Mohanty B, Herath V, Xu F, Mauleon R, Wijaya E, Bajic VB, Bruskiewich R, de Los Reyes BG: Transcriptional regulatory network triggered by oxidative signals configures the early response mechanisms of japonica rice to chilling stress. BMC Plant Biol. 2010, 10: 16-10.1186/1471-2229-10-16.PubMed CentralView ArticlePubMedGoogle Scholar
  8. Shinkawa R, Morishita A, Amikura K, Machida R, Murakawa H, Kuchitsu K, Ishikawa M: Abscisic acid induced freezing tolerance in chilling-sensitive suspension cultures and seedlings of rice. BMC Res Notes. 2013, 6: 351-10.1186/1756-0500-6-351.PubMed CentralView ArticlePubMedGoogle Scholar
  9. Chawade A, Lindlöf A, Olsson B, Olsson O: Global expression profiling of low temperature induced genes in the chilling tolerant japonica rice Jumli Marshi. PLoS One. 2013, 8: e81729-10.1371/journal.pone.0081729.PubMed CentralView ArticlePubMedGoogle Scholar
  10. Chen M, Thelen JJ: Acyl-lipid desaturase2 is required for chilling and freezing tolerance in Arabidopsis. Plant Cell. 2013, 25: 1430-1444. 10.1105/tpc.113.111179.PubMed CentralView ArticlePubMedGoogle Scholar
  11. Kong F, Deng Y, Zhou B, Wang G, Wang Y, Meng Q: A chloroplast-targeted DnaJ protein contributes to maintenance of photosystem II under chilling stress. J Exp Bot. 2014, 65: 143-158. 10.1093/jxb/ert357.PubMed CentralView ArticlePubMedGoogle Scholar
  12. Sunkar R, Chinnusamy V, Zhu J, Zhu J-K: Small RNAs as big players in plant abiotic stress responses and nutrient deprivation. Trends Plant Sci. 2007, 12: 301-309. 10.1016/j.tplants.2007.05.001.View ArticlePubMedGoogle Scholar
  13. Sunkar R, Zhu J-K: Novel and stress-regulated microRNAs and other small RNAs from Arabidopsis. Plant Cell. 2004, 16: 2001-2019. 10.1105/tpc.104.022830.PubMed CentralView ArticlePubMedGoogle Scholar
  14. Zhu Q-H, Spriggs A, Matthew L, Fan L, Kennedy G, Gubler F, Helliwell C: A diverse set of microRNAs and microRNA-like small RNAs in developing rice grains. Genome Res. 2008, 18: 1456-1465. 10.1101/gr.075572.107.PubMed CentralView ArticlePubMedGoogle Scholar
  15. El-Sharkawy MA: Cassava biology and physiology. Plant Mol Biol. 2004, 56: 481-501. 10.1007/s11103-005-2270-7.View ArticlePubMedGoogle Scholar
  16. El-Sharkawy MA: International research on Cassava photosynthesis, productivity, eco-physiology, and responses to environmental stresses in the tropics. Photosynthetica. 2006, 44: 481-512. 10.1007/s11099-006-0063-0.View ArticleGoogle Scholar
  17. Leyva-González MA, Ibarra-Laclette E, Cruz-Ramírez A, Herrera-Estrella L: Functional and transcriptome analysis reveals an acclimatization strategy for abiotic stress tolerance mediated by Arabidopsis NF-YA family members. PLoS One. 2012, 7: e48138-10.1371/journal.pone.0048138.PubMed CentralView ArticlePubMedGoogle Scholar
  18. Carvallo MA, Pino M-T, Jeknic Z, Zou C, Doherty CJ, Shiu S-H, Chen THH, 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: 3807-3819. 10.1093/jxb/err066.PubMed CentralView ArticlePubMedGoogle Scholar
  19. Preston JC, Sandve SR: Adaptation to seasonality and the winter freeze. Front Plant Sci. 2013, 4: 167-PubMed CentralPubMedGoogle Scholar
  20. Zhang T, Zhao X, Wang W, Pan Y, Huang L, Liu X, Zong Y, Zhu L, Yang D, Fu B: Comparative transcriptome profiling of chilling stress responsiveness in two contrasting rice genotypes. PLoS One. 2012, 7: e43274-10.1371/journal.pone.0043274.PubMed CentralView ArticlePubMedGoogle Scholar
  21. Usadel B, Bläsing OE, Gibon Y, Poree F, Höhne M, Günter M, Trethewey R, Kamlage B, Poorter H, Stitt M: Multilevel genomic analysis of the response of transcripts, enzyme activities and metabolites in Arabidopsis rosettes to a progressive decrease of temperature in the non-freezing range. Plant Cell Environ. 2008, 31: 518-547. 10.1111/j.1365-3040.2007.01763.x.View ArticlePubMedGoogle Scholar
  22. An D, Yang J, Zhang P: Transcriptome profiling of low temperature-treated Cassava apical shoots showed dynamic responses of tropical plant to cold stress. BMC Genomics. 2012, 13: 64-10.1186/1471-2164-13-64.PubMed CentralView ArticlePubMedGoogle Scholar
  23. Bursać S, Brdovčak MC, Pfannkuchen M, Orsolić I, Golomb L, Zhu Y, Katz C, Daftuar L, Grabušić K, Vukelić I, Filić V, Oren M, Prives C, Volarevic S: Mutual protection of ribosomal proteins L5 and L11 from degradation is essential for p53 activation upon ribosomal biogenesis stress. Proc Natl Acad Sci U S A. 2012, 109: 20467-20472. 10.1073/pnas.1218535109.PubMed CentralView ArticlePubMedGoogle Scholar
  24. Deisenroth C, Zhang Y: The ribosomal protein-Mdm2-p53 pathway and energy metabolism: bridging the Gap between feast and famine. Genes Cancer. 2011, 2: 392-403. 10.1177/1947601911409737.PubMed CentralView ArticlePubMedGoogle Scholar
  25. Provart N: Correlation networks visualization. Front Plant Sci. 2012, 3: 240-10.3389/fpls.2012.00240.PubMed CentralView ArticlePubMedGoogle Scholar
  26. Kumar SV, Wigge PA: H2A.Z-containing nucleosomes mediate the thermosensory response in Arabidopsis. Cell. 2010, 140: 136-147. 10.1016/j.cell.2009.11.006.View ArticlePubMedGoogle Scholar
  27. Zhang Q, Xiao H, Chai SC, Hoang QQ, Lu H: Hydrophilic residues are crucial for ribosomal protein L11 (RPL11) interaction with zinc finger domain of MDM2 and p53 protein activation. J Biol Chem. 2011, 286: 38264-38274. 10.1074/jbc.M111.277012.PubMed CentralView ArticlePubMedGoogle Scholar
  28. Bhat KP, Itahana K, Jin A, Zhang Y: Essential role of ribosomal protein L11 in mediating growth inhibition-induced p53 activation. EMBO J. 2004, 23: 2402-2412. 10.1038/sj.emboj.7600247.PubMed CentralView ArticlePubMedGoogle Scholar
  29. Dai M-S, Shi D, Jin Y, Sun X-X, Zhang Y, Grossman SR, Lu H: Regulation of the MDM2-p53 pathway by ribosomal protein L11 involves a post-ubiquitination mechanism. J Biol Chem. 2006, 281: 24304-24313. 10.1074/jbc.M602596200.PubMed CentralView ArticlePubMedGoogle Scholar
  30. Harding HP, Novoa I, Zhang Y, Zeng H, Wek R, Schapira M, Ron D: Regulated translation initiation controls stress-induced gene expression in mammalian cells. Mol Cell. 2000, 6: 1099-1108. 10.1016/S1097-2765(00)00108-8.View ArticlePubMedGoogle Scholar
  31. Smirnova JB, Selley JN, Sanchez-Cabo F, Carroll K, Eddy AA, McCarthy JEG, Hubbard SJ, Pavitt GD, Grant CM, Ashe MP: Global gene expression profiling reveals widespread yet distinctive translational responses to different eukaryotic translation initiation factor 2B-targeting stress pathways. Mol Cell Biol. 2005, 25: 9340-9349. 10.1128/MCB.25.21.9340-9349.2005.PubMed CentralView ArticlePubMedGoogle Scholar
  32. Casati P, Walbot V: Gene expression profiling in response to ultraviolet radiation in maize genotypes with varying flavonoid content. Plant Physiol. 2003, 132: 1739-1754. 10.1104/pp.103.022871.PubMed CentralView ArticlePubMedGoogle Scholar
  33. Wurtmann EJ, Wolin SL: RNA under attack: cellular handling of RNA damage. Crit Rev Biochem Mol Biol. 2009, 44: 34-49. 10.1080/10409230802594043.PubMed CentralView ArticlePubMedGoogle Scholar
  34. Grousl T, Ivanov P, Malcova I, Pompach P, Frydlova I, Slaba R, Senohrabkova L, Novakova L, Hasek J: Heat shock-induced accumulation of translation elongation and termination factors precedes assembly of stress granules in S. cerevisiae. PLoS One. 2013, 8: e57083-10.1371/journal.pone.0057083.PubMed CentralView ArticlePubMedGoogle Scholar
  35. Carroll AJ: The Arabidopsis cytosolic ribosomal proteome: from form to function. Front Plant Sci. 2013, 4: 32-PubMed CentralPubMedGoogle Scholar
  36. Gerashchenko MV, Lobanov AV, Gladyshev VN: Genome-wide ribosome profiling reveals complex translational regulation in response to oxidative stress. Proc Natl Acad Sci U S A. 2012, 109: 17394-17399. 10.1073/pnas.1120799109.PubMed CentralView ArticlePubMedGoogle Scholar
  37. Falcone Ferreyra ML, Pezza A, Biarc J, Burlingame AL, Casati P: Plant L10 ribosomal proteins have different roles during development and translation under ultraviolet-B stress. Plant Physiol. 2010, 153: 1878-1894. 10.1104/pp.110.157057.View ArticlePubMedGoogle Scholar
  38. Iyer LM, Abhiman S, de Souza RF, Aravind L: Origin and evolution of peptide-modifying dioxygenases and identification of the wybutosine hydroxylase/hydroperoxidase. Nucleic Acids Res. 2010, 38: 5261-5279. 10.1093/nar/gkq265.PubMed CentralView ArticlePubMedGoogle Scholar
  39. Uberto R, Moomaw EW: Protein similarity networks reveal relationships among sequence, structure, and function within the cupin superfamily. PLoS One. 2013, 8: e74477-10.1371/journal.pone.0074477.PubMed CentralView ArticlePubMedGoogle Scholar
  40. Choi Y, Lee Y, Kim SY, Lee Y, Hwang J-U: Arabidopsis ROP-interactive CRIB motif-containing protein 1 (RIC1) positively regulates auxin signalling and negatively regulates abscisic acid (ABA) signalling during root development. Plant Cell Environ. 2013, 36: 945-955. 10.1111/pce.12028.View ArticlePubMedGoogle Scholar
  41. Yang Z, Fu Y: ROP/RAC GTPase signaling. Curr Opin Plant Biol. 2007, 10: 490-494. 10.1016/j.pbi.2007.07.005.PubMed CentralView ArticlePubMedGoogle Scholar
  42. Li C, Rudi H, Stockinger EJ, Cheng H, Cao M, Fox SE, Mockler TC, Westereng B, Fjellheim S, Rognli OA, Sandve SR: Comparative analyses reveal potential uses of Brachypodium distachyon as a model for cold stress responses in temperate grasses. BMC Plant Biol. 2012, 12: 65-10.1186/1471-2229-12-65.PubMed CentralView ArticlePubMedGoogle Scholar
  43. Koehler G, Wilson RC, Goodpaster JV, Sønsteby A, Lai X, Witzmann FA, You J-S, Rohloff J, Randall SK, Alsheikh M: Proteomic study of low-temperature responses in strawberry cultivars (Fragaria x ananassa) that differ in cold tolerance. Plant Physiol. 2012, 159: 1787-1805. 10.1104/pp.112.198267.PubMed CentralView ArticlePubMedGoogle Scholar
  44. Yamori W, Noguchi K, Hikosaka K, Terashima I: Phenotypic plasticity in photosynthetic temperature acclimation among crop species with different cold tolerances. Plant Physiol. 2010, 152: 388-399. 10.1104/pp.109.145862.PubMed CentralView ArticlePubMedGoogle Scholar
  45. Medina J, Bargues M, Terol J, Pérez-Alonso M, Salinas J: The Arabidopsis CBF gene family is composed of three genes encoding AP2 domain-containing proteins whose expression Is regulated by low temperature but not by abscisic acid or dehydration. Plant Physiol. 1999, 119: 463-470. 10.1104/pp.119.2.463.PubMed CentralView ArticlePubMedGoogle Scholar
  46. Haake V, Cook D, Riechmann JL, Pineda O, Thomashow MF, Zhang JZ: Transcription factor CBF4 is a regulator of drought adaptation in Arabidopsis. Plant Physiol. 2002, 130: 639-648. 10.1104/pp.006478.PubMed CentralView ArticlePubMedGoogle Scholar
  47. 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: 1043-1054. 10.1101/gad.1077503.PubMed CentralView ArticlePubMedGoogle Scholar
  48. Fursova OV, Pogorelko GV, Tarasov VA: Identification of ICE2, a gene involved in cold acclimation which determines freezing tolerance in Arabidopsis thaliana. Gene. 2009, 429: 98-103. 10.1016/j.gene.2008.10.016.View ArticlePubMedGoogle Scholar
  49. Guo Y, Xiong L, Ishitani M, Zhu J-K: An Arabidopsis mutation in translation elongation factor 2 causes superinduction of CBF/DREB1 transcription factor genes but blocks the induction of their downstream targets under low temperatures. Proc Natl Acad Sci U S A. 2002, 99: 7786-7791. 10.1073/pnas.112040099.PubMed CentralView ArticlePubMedGoogle Scholar
  50. Ishitani M, Xiong L, Lee H, Stevenson B, Zhu JK: HOS1, a genetic locus involved in cold-responsive gene expression in arabidopsis. Plant Cell. 1998, 10: 1151-1161. 10.1105/tpc.10.7.1151.PubMed CentralView ArticlePubMedGoogle Scholar
  51. Miura K, Jin JB, Lee J, Yoo CY, Stirm V, Miura T, Ashworth EN, Bressan RA, Yun D-J, Hasegawa PM: SIZ1-mediated sumoylation of ICE1 controls CBF3/DREB1A expression and freezing tolerance in Arabidopsis. Plant Cell. 2007, 19: 1403-1414. 10.1105/tpc.106.048397.PubMed CentralView ArticlePubMedGoogle Scholar
  52. Dong C-H, Agarwal M, Zhang Y, Xie Q, Zhu J-K: The negative regulator of plant cold responses, HOS1, is a RING E3 ligase that mediates the ubiquitination and degradation of ICE1. Proc Natl Acad Sci U S A. 2006, 103: 8281-8286. 10.1073/pnas.0602874103.PubMed CentralView ArticlePubMedGoogle Scholar
  53. Zhang L, Li Z, Li J, Wang A: Ectopic overexpression of SsCBF1, a CRT/DRE-binding factor from the nightshade plant Solanum lycopersicoides, confers freezing and salt tolerance in transgenic Arabidopsis. PLoS One. 2013, 8: e61810-10.1371/journal.pone.0061810.PubMed CentralView ArticlePubMedGoogle Scholar
  54. Zhang X, Fowler SG, Cheng H, Lou Y, Rhee SY, Stockinger EJ, Thomashow MF: Freezing-sensitive tomato has a functional CBF cold response pathway, but a CBF regulon that differs from that of freezing-tolerant Arabidopsis. Plant J. 2004, 39: 905-919. 10.1111/j.1365-313X.2004.02176.x.View ArticlePubMedGoogle Scholar
  55. McKhann HI, Gery C, Bérard A, Lévêque S, Zuther E, Hincha DK, De Mita S, Brunel D, Téoulé E: Natural variation in CBF gene sequence, gene expression and freezing tolerance in the Versailles core collection of Arabidopsis thaliana. BMC Plant Biol. 2008, 8: 105-10.1186/1471-2229-8-105.PubMed CentralView ArticlePubMedGoogle Scholar
  56. Doherty CJ, Van Buskirk HA, Myers SJ, Thomashow MF: Roles for Arabidopsis CAMTA transcription factors in cold-regulated gene expression and freezing tolerance. Plant Cell. 2009, 21: 972-984. 10.1105/tpc.108.063958.PubMed CentralView ArticlePubMedGoogle Scholar
  57. Dong C-H, Hu X, Tang W, Zheng X, Kim YS, Lee B, Zhu J-K: A putative Arabidopsis nucleoporin, AtNUP160, is critical for RNA export and required for plant tolerance to cold stress. Mol Cell Biol. 2006, 26: 9533-9543. 10.1128/MCB.01063-06.PubMed CentralView ArticlePubMedGoogle Scholar
  58. Gong Z, Dong C-H, Lee H, Zhu J, Xiong L, Gong D, Stevenson B, Zhu J-K: A DEAD box RNA helicase is essential for mRNA export and important for development and stress responses in Arabidopsis. Plant Cell. 2005, 17: 256-267. 10.1105/tpc.104.027557.PubMed CentralView ArticlePubMedGoogle Scholar
  59. Kidokoro S, Maruyama K, Nakashima K, Imura Y, Narusaka Y, Shinwari ZK, Osakabe Y, Fujita Y, Mizoi J, Shinozaki K, Yamaguchi-Shinozaki K: The phytochrome-interacting factor PIF7 negatively regulates DREB1 expression under circadian control in Arabidopsis. Plant Physiol. 2009, 151: 2046-2057. 10.1104/pp.109.147033.PubMed CentralView ArticlePubMedGoogle Scholar
  60. Hoagland DR, Arnon DI: The Water-Culture Method for Growing Plants without Soil. 1950, Department of Agriculture experimental station, CaliforniaGoogle Scholar
  61. Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009, 10: R25-10.1186/gb-2009-10-3-r25.PubMed CentralView ArticlePubMedGoogle Scholar
  62. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L: Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010, 28: 511-515. 10.1038/nbt.1621.PubMed CentralView ArticlePubMedGoogle Scholar
  63. Fahlgren N, Howell MD, Kasschau KD, Chapman EJ, Sullivan CM, Cumbie JS, Givan SA, Law TF, Grant SR, Dangl JL, Carrington JC: High-throughput sequencing of Arabidopsis microRNAs: evidence for frequent birth and death of MIRNA genes. PLoS One. 2007, 2: e219-10.1371/journal.pone.0000219.PubMed CentralView ArticlePubMedGoogle Scholar
  64. Goodstein DM, Shu S, Howson R, Neupane R, Hayes RD, Fazo J, Mitros T, Dirks W, Hellsten U, Putnam N, Rokhsar DS: Phytozome: a comparative platform for green plant genomics. Nucleic Acids Res. 2012, 40 (Database issue): D1178-D1186. 10.1093/nar/gkr944.PubMed CentralView ArticlePubMedGoogle Scholar
  65. Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B. 1995, 57: 289-300.Google Scholar

Copyright

© Zeng et al. 2014

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/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0) applies to the data made available in this article, unless otherwise stated.