Transcriptome profiling of low temperature-treated cassava apical shoots showed dynamic responses of tropical plant to cold stress

Background Cassava is an important tropical root crop adapted to a wide range of environmental stimuli such as drought and acid soils. Nevertheless, it is an extremely cold-sensitive tropical species. Thus far, there is limited information about gene regulation and signalling pathways related to the cold stress response in cassava. The development of microarray technology has accelerated the study of global transcription profiling under certain conditions. Results A 60-mer oligonucleotide microarray representing 20,840 genes was used to perform transcriptome profiling in apical shoots of cassava subjected to cold at 7°C for 0, 4 and 9 h. A total of 508 transcripts were identified as early cold-responsive genes in which 319 sequences had functional descriptions when aligned with Arabidopsis proteins. Gene ontology annotation analysis identified many cold-relevant categories, including 'Response to abiotic and biotic stimulus', 'Response to stress', 'Transcription factor activity', and 'Chloroplast'. Various stress-associated genes with a wide range of biological functions were found, such as signal transduction components (e.g., MAP kinase 4), transcription factors (TFs, e.g., RAP2.11), and reactive oxygen species (ROS) scavenging enzymes (e.g., catalase 2), as well as photosynthesis-related genes (e.g., PsaL). Seventeen major TF families including many well-studied members (e.g., AP2-EREBP) were also involved in the early response to cold stress. Meanwhile, KEGG pathway analysis uncovered many important pathways, such as 'Plant hormone signal transduction' and 'Starch and sucrose metabolism'. Furthermore, the expression changes of 32 genes under cold and other abiotic stress conditions were validated by real-time RT-PCR. Importantly, most of the tested stress-responsive genes were primarily expressed in mature leaves, stem cambia, and fibrous roots rather than apical buds and young leaves. As a response to cold stress in cassava, an increase in transcripts and enzyme activities of ROS scavenging genes and the accumulation of total soluble sugars (including sucrose and glucose) were also detected. Conclusions The dynamic expression changes reflect the integrative controlling and transcriptome regulation of the networks in the cold stress response of cassava. The biological processes involved in the signal perception and physiological response might shed light on the molecular mechanisms related to cold tolerance in tropical plants and provide useful candidate genes for genetic improvement.


Background
Cassava (Manihot esculenta Crantz) is widely cultivated for its starchy storage roots and is a staple food and animal feed in tropical and sub-tropical areas [1]. It is also considered to be an important source of modified starches and bioethanol in China and other Southeast Asian countries [2,3]. Nevertheless, as a tropical root crop, cassava is native to a warm habitat and is categorized as a cold-sensitive species [4]. Thus, low temperatures and frozen conditions are the most important limiting factors for its geographical location and productivity. In the subtropics, where unpredictable cold weather occurs occasionally, it is important to protect the storage roots and propagation stems from chilling stress. For example, the unprecedented freezing disaster occurred in Southern China in January 2008 caused great damage to cassava stem seeds and led to yield reduction in Guangxi, Guangdong and other provinces, resulting in a loss of a billion Chinese Yuan [5]. Furthermore, to ensure a prolonged growth period (i.e., early planting and late harvesting) in the high latitude regions, novel cassava cultivars with improved cold tolerance are in demand.
Under low temperature below 10°C, many species of tropical or subtropical origin are typically injured or killed and show various symptoms of chilling injury due to the inability to adapt to non-freezing low temperatures [6]. For example, cassava exhibits obvious symptoms of damage at these temperatures, including delayed sprouting of the stem cutting, yield decrease, reduced leaf expansion, chlorosis and even necrosis in its leaves [7]. Low temperatures have also been recognized as an important facilitator of decreases in nutrient absorption rates (e.g., Boron), reductions in the leaf photosynthetic rate, and the inhibition of plant growth [4]. In addition, the physiological status of cold-stressed plants is also altered, such as transient increases in hormone levels (e.g., ABA) [8] and changes in membrane lipid composition [9]. Furthermore, the accumulation of compatible osmolytes, such as soluble sugars, betaine, and proline [10][11][12], and increases in the level of antioxidants [13] are also occurred. In contrast, temperate plants can withstand freezing temperatures following a period of low, but non-freezing, temperatures, a process called cold acclimation [14]. The mechanisms of cold acclimation have been extensively investigated in the model plant Arabidopsis thaliana [15], as well as in other important crop species such as maize and barley [16,17]. The most extensively studied mechanism involves a class of ethylene response factors (ERFs) known as the dehydration-responsive element-binding proteins/C-repeat-binding factors (DREBs/CBFs). They interact with the dehydration-responsive element/Crepeat element (DRE/CRT) in the promoters of their downstream target genes to execute a highly coordinated transcriptional response to low temperature signals. Overall, less information is available for tropical plants (e.g., cassava). Therefore, it would be important to increase our understanding of the physiological, biochemical, and molecular characteristics that are associated with cold stress in tropical cassava, especially during the early stages. Moreover, the function of homologous genes in Arabidopsis could be used to predict the function of cassava genes which may share the same biological functions or common regulatory scenarios among different plant species [18].
With the development of molecular technologies and 'omics' tools, microarray studies have become a useful strategy for the global analysis of plant gene expression. Using cDNA microarrays or whole genome arrays, the abiotic stress responses of Arabidopsis and other plants have been widely analyzed. For examples, the expression patterns of genes were identified in Arabidopsis and rice under conditions of drought, cold, high-salinity or abscisic acid treatment [19,20]. The transcriptomic identification of candidate genes involved in sunflower responses to chilling and salt stresses has also been reported [21]. The long-oligonucleotide microarray, which has many advantages compared to other microarray technologies, reliably detects transcript ratios even at one copy per cell in complicated biological samples, and could distinguish different members of gene families [22,23]. Thus, it is appropriate to explore the early cold stress responsive genes of cassava using a 60-mer long-oligonucleotide microarray approach.
Thus far, several studies have been reported in cassava using genomic tools for developing new knowledge and technologies. For instance, EST and cDNA libraries have been constructed in cassava for the identification of starch biosynthesis and biotic/abiotic-responsive genes or genes corresponding to the different developmental stages of various tissues [24][25][26][27][28]. Studies using cassava cDNA microarrays to analyze gene expression in cassava plants subjected to post-harvest physiological deterioration (PPD) or Xanthomonas axonopodis infection have also been reported [29,30]. Moreover, gene expression profiling related to the different growth stages of cassava storage roots was recently conducted using a long-oligonucleotide microarray [31]. Fortunately, a draft genome sequence for cassava has been released and recently updated (http://www.phytozome.net/cassava), which greatly facilitates cassava research worldwide. Nevertheless, comprehensive genome-wide expression profiling data for cassava under cold stress treatment are still lacking, which limits our capacity to decipher the molecular mechanisms related to various stress responses. Additionally, to stabilize cassava yield even under unfavorable growth conditions, it is necessary to study the changes of pan-genome genes or pathways in order to identify key candidates for improving cold tolerance in cassava.
In this study, a gene expression profiling of apical shoots from cassava subjected to low temperature was conducted using a custom-designed, 60-mer oligonucleotide microarray covering 20,840 cassava transcripts. In total, 508 differentially expressed genes (DEGs) were identified, and their biological functions were characterized through gene ontology (GO) annotation and KEGG pathway analysis. Our study could increase the understanding of cassava gene regulation under the condition of cold stress and reveal novel approaches to improve cold tolerance through genetic engineering.

Results and discussion
Phenotypic and physiological changes of cassava in response to cold Cassava is a typical tropical crop adapted to warm climates. Under normal conditions, 3-month-old cassava (TMS60444) plants have vigorous apical buds and fully expanded mature leaves with upward leaf fingers ( Figure  1A, Left panel). When a low-temperature treatment is applied to cassava at 7°C under weak light for 4 h, the plants displayed visible morphological changes, including weak dehydration and wilting of the apical buds and leaves, especially in the buds, as well as downward leaf fingers ( Figure 1A, Middle panel). Under prolonged exposure to the stress (9 h), the whole plants exhibited obvious phenotypic damages, including softening and downward bending of the petioles, loss of strength in the immature stems and more severe wilting leaves (Figure 1A, Right panel). Generally, the treated plants could be partially recovered and resumed growth when transferred into normal conditions (data not shown), indicating that the cellular changes caused by a 9 h cold treatment were reversible. Since apical shoots show extremely sensitive response to low temperature (Additional file 1), the apical shoots (about 4-6 cm) of 3month-old plants, which were treated for 0, 4 and 9 h, were used as material for microarray hybridization.
In addition to morphological changes, a clear impact on ultra-structure of chloroplasts was also observed in leaves of cold-stressed cassava plants. Compared to normal chloroplasts in unstressed leaves ( Figure 1B), the stacking numbers of thylakoids were dramatically reduced and less organized, and the starch granules had almost disappeared ( Figure 1C). The ultra-structural changes in cellular organelles indicate that cold stress might exert an adverse impact on plant growth, causing a reduced photosynthetic rate under the stress.
Cold stress also leads to obviously physiological changes in cassava. Two stress responsive metabolites, namely malondialdehyde (MDA) and proline, were monitored. MDA is considered to be the final product of lipid peroxidation in the plant cell membrane and is an important indicator of membrane system injuries and cellular metabolism deterioration [32]. It has also been reported that MDA contents change during cold stress treatment in plants [33]. In our experiment, compared to the control level, the MDA concentrations decreased rapidly to 50% after 4 h but increased almost 25% after 9 h; the observed levels had experienced a nearly twofold change after 24 h of stress treatment ( Figure 1D). Decrease in MDA content at 4 h might be due to a transient stress response of cassava to the low temperature shock. But, the prolonged treatment (9 h) finally led to cell damages and MDA accumulation. The accumulation of proline is frequently associated with whole plant tolerance to chilling and other stresses [34]. Unlike that of MDA, the concentration of proline increased rapidly and achieved nearly a three-fold change after 4 h; it then decreased slightly but remained at a high level until 24 h ( Figure 1E). These results were consistent with previous studies that proline accumulated in leaves exposure to cold, salt, and other stresses in Arabidopsis and other plants [35]. The accumulation of these metabolites is a good indication that the stress-treated plants were actively mounting a stress response during the periods when they were subjected to cold stress.

Microarray for transcriptomic analysis of low-temperature treated cassava
To explore the transcriptomic changes of cassava in response to cold stress, a custom long-oligonucleotide (60-mer) microarray generated by the Agilent SurePrint ink-jet technology was applied in the current study. The detailed protocol for designing the cassava microarray has been described by Yang et al [31] To validate the microarray quality, the signal-to-noise ratio (SNR) for each spot was calculated as described by Leiske et al [23] The average percentages of acceptable spots (SNR > 2.6) and highquality spots (SNR > 10) were 94.53 ± 1.84% and 79.24 ± 1.76% in all 9 arrays, which demonstrated the overall reproducibility and high quality of the array. The hierarchical clustering results indicated the clear separation of the samples from different time points (Figure 2A), suggesting that the entire experiment from sample collection to data extraction was reproducible and reliable. Five house-keeping genes (beta-actin, c15, EF1a, RuBisCO small chain precursor, and TUB6) were used as the internal controls and their expressions were consistent in the samples at different time points (Additional file 2), confirming the reliability and accuracy of these microarrays.

Transcriptomic responses to cold treatment
Before normalization, the signal intensities of each feature were filtered against negative controls on the array, (D) MDA contents in cassava apical leaves exposed to 7°C for 0, 4, 9, and 24 h. (E) Proline accumulation in cassava apical leaves exposed to 7°C for 0, 4, 9, and 24 h. The double asterisks indicate a statistically significant difference (p < 0.01) for the data of the stress-treated samples compared to those of the unstressed samples. The mean values are calculated from three biological replicates; the error bars represent the standard error of the mean (SEM). Bar = 500 nm. The red numbers are the total numbers of differentially expressed genes (DEGs); the percentages in parentheses were calculated as the ratio of regulated genes to the total number of cold-regulated genes (508). It should be noted that one gene was up-regulated at 4 h/0 h but down-regulated at 9 h/4 h on our array. and 69.87 ± 0.72%, 68.14 ± 2.59%, and 71.50 ± 5.75% of probes on the array were expressed at 0, 4 and 9 h, respectively. Two filtering criteria were used to define differentially expressed genes (DEGs) in our data analysis: a two-fold change in transcript levels among every two time points and a P value ≤ 0.05. To analyze the similarities and differences among the cold-responsive transcriptomes, a hierarchical clustering was prepared to represent the transcripts of all the differentially expressed probes at the 3 replicates of 0, 4 and 9 h. These results indicated significant differences in the gene expression profiles between the treatments of 9 h and 0 h or 9 h and 4 h, in contrast to the relatively high similarity of the expression profiles between 0 h and 4 h ( Figure 2A). Among the DEGs, 58, 252, and 92 genes were up-regulated and 13, 210, and 36 were down-regulated by analyzing 4 h/0 h, 9 h/0 h, and 9 h/4 h, respectively ( Figure 2B). Notably, compared to only 71 and 128 genes showed differentially expressed at 4 h/0 h and 9 h/4 h, 474 DEGs were found at 9 h/0 h, suggesting a prolonged stress treatment (9 h) could trigger more stress-related gene expression. Moreover, about 80% of 4 h/0 h and 75% of 9 h/4 h DEGs overlapped with those of 9 h/0 h, which indicated a strong linkage among the 3 stressed comparison points and a progressive biological process.
In total, 508 DEGs, including 292 cold-inducible and 215 cold-repressed genes, were identified on our array; only one gene (CUST_14887) was up-regulated at 4 h/0 h but down-regulated at 9 h/4 h (Additional file 3). These represent about 3.5% of the total expressed cassava transcripts on the array. This result indicated that even for tropical crops (e.g., cassava), which have adapted to warm climates, possess a large amount of cold-responsive genes as other temperate plants (e.g., Arabidopsis). Therefore, the differences in the ability of tropical species and temperate species to tolerate cold might not totally be due to the amount of genes responsive to cold stress. Instead, we speculate that plants require to timely and coordinately mobilize the biological functions and regulatory networks of stress-responsive genes against stress. Consistent with our observations, a recent study confirmed such scenario even between two close plant species [36].

Gene ontology clustering of cold-regulated genes
To explore the biological functions of cold-responsive genes, a total of 20,840 sequences were used as a query to perform an alignment with Arabidopsis proteins. Of these, 20,450 sequences had hits with 11,995 Arabidopsis proteins, and 14,813 (about 71.8%) sequences had hits with 9,158 Arabidopsis proteins at an E value ≤ 1e-5. After searching the descriptions of the 508 DEGs, 319 queries had the highest homologies with annotated proteins according to The Arabidopsis Information Resource (TAIR), which accounted for 62.8% of the responsive genes (Additional file 4). To determine the detailed function of the 319 hits with Arabidopsis gene locus identifiers, a GO annotation was performed with the GO terms of TAIR (http://www.Arabidopsis.org/ tools/bulk/go/index.jsp), in which 370 GO IDs were hit by 291 gene locus identifiers (Additional file 5). Therefore, functional clusters that were classified according to the three components ('Biological Process', 'Molecular Function', and 'Cellular Component') are presented and discussed (Additional files 6, 7, and 8, respectively). Based on the TAIR percent analysis, the categories of 'Response to abiotic and biotic stimulus' and 'Response to stress', 'Transcription factor activity', and 'Chloroplast' were the major cold-related GO Slims that strongly attracted our attention.
Cold-responsive genes related to 'Response to abiotic and biotic stimulus' and 'Response to stress' In the GO terms of 'Biological Process', the categories of 'Response to abiotic and biotic stimulus' and 'Response to stress', both accounted for 13.36% of the total GOs, were suggested to be the most relevant to cold stress ( Figure 3A). The largest proportions of 'Biological Process' terms were annotated as 'Other cellular processes' (48.09%), 'Other metabolic processes' (43.51%), 'Unknown biological process' (33.59%), and 'Other biological processes' (14.12%), indicating comprehensive changes in cassava gene expression ( Figure 3A). In addition, 'Protein metabolism' was another noticeable major category (15.27%) supposed to be also involved in the cold stress response [37].
There were 34 and 33 genes assigned to the GO terms of 'Response to abiotic and biotic stimulus' and 'Response to stress', respectively (Additional file 9). Because 27 genes were common in both categories, we combined the two categories for analysis and discussion (Table 1). Among the overlapping genes, 20 up-regulated genes and 7 down-regulated genes were included. Their encoding proteins cover a wide range of biological functions, including transcription factor (TFs, e.g., AP2 domain transcription factor), protein kinases (e.g., proline extension-like receptor kinase 1), small heat shock proteins (e.g., heat shock protein 18.2), ROS scavenging enzymes [e.g., catalase 2 (CAT2)], fatty acid desaturases [e.g., fatty acid biosynthesis 2 (SSI2)], signaling molecular proteins (e.g., MAP kinase 4), programmed cell death (PCD)-related gene (e.g., chitinase), and genes involved in the gibberellin acid (GA) and jasmonic acid (JA) metabolism pathway [e.g., gibberellic acid insensitive (GAI)]. Among the genes specifically related to the 'Response to abiotic and biotic stimulus', protein degradation genes (e.g., ubiquitin-protein ligase) and early flowering 4 were up-regulated. Meanwhile, genes involved in ABA signal transduction [e.g., ABA-responsive element binding protein 3 (AREB3)] and disease resistance protein (CC-NBS-LRR), which were downregulated by cold, were uniquely identified in GO Slim as 'Response to stress' ( Table 1). Most of these coldresponsive genes were either up-regulated or downregulated at 9 h/0 h, except 4 genes encoding the calcium ion binding protein (TCH2), ethylene response transcription activator 2 (ERF2), MAP kinase kinase 9, and pseudo-response regulator 5. It is notable that 27 out of these 40 stress-associated genes were up-regulated, indicating that the up-regulation of cold-responsive transcripts by low temperature may play a crucial role in the stress tolerance of plants [38]. Of their homologous genes in Arabidopsis, the MAP kinase 4, beta-ketoacyl-CoA synthase, CAT2, and TCH2 were previously reported to be involved in the cold stress response process. For example, MAP kinase 4, which is required for cytokinesis, was implicated in cold and salt stress tolerance [39]; beta-ketoacyl-CoA synthase, which is involved in the biosynthesis of VLCFA (very long chain fatty acids), was also responsive to cold and other osmotic stresses [40]. CAT2, a target of plastid, could increase its transcript abundance in response to cold stress [41]. TCH2 has been implicated in calcium signaling in response to diverse stimuli, such as cold, light, pathogens, and touch [42]. Gibberellin 2beta-dioxygenase, which degrades active GAs, is involved in cold response [43] and up-regulation of its transcript implied that GA level may be reduced under cold stress. Additionally, calcineurin B-like (CBL)-interacting protein kinase, which has protein serine/threonine kinase activity, has been reported to be involved in the defense response to fungus [44]. These results revealed that there are considerable conserved and common components in cold stress response mechanisms across plant species, including temperate and tropical plants.
Besides, a large number of genes with known function were firstly identified in the cold response. The AP2 domain transcription factor (CUST_2332), a member of the ERF (ethylene response factor) subfamily B-3 of ERF/AP2 transcription factor family, has been reported involved in response to fungus and chitin in Arabidopsis. In this study, the transcript abundance of the AP2 domain transcription factor was up-regulated 16-fold by cold, indicating that it may play a significant role in the regulation of novel signaling pathways to enhance cassava cold tolerance. The expression of a stearoyl-ACP desaturase gene SSI2, which involves in fatty acid desaturation, was up-regulated almost 10-fold under cold stress. This change suggested that the gene might be participated in altering membrane lipid composition to enhance membrane fluidity. Maternal effect embryo arrest 14 (MEE14), which required for embryonic development ending in seed dormancy, was also highly induced in this study. Interestingly, three transcripts (early flowering 4, gibberellin 2-beta-dioxygenase, and pseudo-response regulator 5) that respond to red or farred light in Arabidopsis were cold-inducible, indicating a crosstalk between light signaling and cold stress response.
In addition, a large number of genes with unknown function were also firstly recorded as cold responsive genes. For instance, four members of MLP-like protein (CUST_6563, CUST_15946, CUST_11303, and CUST_10998) were all up-regulated by 9 h cold treatment. Alpha-crystallin domain 31.2, Bet v I allergen family protein, and homology of yeast autophagy 18 were highly induced more than 5-fold by cold. Two members of RING-H2 finger A2A zinc ion binding proteins were also found to be down-regulated. Taken together, cassava might evolve not only conserved but also specific molecular mechanisms related to stress signaling and response.

Transcription factors response to cold stress
For the 'Molecular Function' GO terms, 'Transferase activity', 'Hydrolase activity', and 'Transcription factor activity' were the three major categories. They accounted for 13.28%, 11.72%, and 11.33% of the total GOs, respectively ( Figure 3B). Among them, genes associated with 'Transcription factor activity' may be central regulators involved in early cold signal transduction that trigger a cascade of downstream gene expression.
Transcription factors (TFs) play a significant role in plant development and stress tolerance [38]. To identify the transcription factors involved in the cold stress response, we surveyed the biological functions of putative TFs that were differentially expressed in cassava under the cold treatment. A total of 32 genes were identified as TFs when the 319 identifiers were annotated in TAIR. The number of up-regulated TFs (15) was nearly equal to the number of down-regulated ones (17), suggesting that both transcriptional activation and repression are involved. The Arabidopsis genome has more than 2,657 predicted TFs belonging to 81 gene families in the Plant Transcription Factor Database (PlnTFDB, http://plntfdb.bio.uni-potsdam.de/v3.0/index.php?sp_i-d=ATH) [45]. According to this classification, 31 coldresponsive transcription factors fell into 17 families, except 1 (CUST_19860) with non-classification (Table  2), and they were annotated with the detailed GO Slims (Additional file 10). In Arabidopsis, at least five TF families have been reported to be involved in the cold stress response process, including AP2-EREBP (APE-TALA2/ET-Responsive Element Binding Protein), MYB (Myeloblastosis), NAC (NAM, ATAF1/2, CUC2), bHLH (basic Helix-Loop-Helix), and WRKY (named after the WRKY amino acid motif) [46].
In our study, AP2-EREBP, MYB, and GRAS were the three major TF families, containing six, five, and five genes of the 32 cold-responsive TFs, respectively. The AP2-EREBP family plays a major role in the early stages of the cold response, as evidenced by many well-characterized DREBs/CBFs cold regulatory pathway. CBF proteins, belonging to A-1 subfamily of ERF/AP2 TF family, are major regulators that function in activating coldregulated effectors in Arabidopsis and other plants [15]. In our study, all members of AP2-EREBP family are grouped to different subfamilies of ERF/AP2 TF family. For instance, RAP2.11, ERF2, AP2 domain transcription factor (CUST_20765), and RAP2.4 encode members of B-6, B-3, A-5, and A-6 subfamilies of ERF/AP2 TF family, respectively. Some members have been reported to show response to cold stress in Arabidopsis, such as RAP2.11 and RAP2.4 [47,48]. However, others were firstly identified to be involved in cold response, such as AP2 domain TF (CUST_2332) and ERF2. Importantly, most members of this family were highly induced at 9 h during cold treatment, which was confirmed by realtime RT-PCR ( Figure 4A), suggesting a relatively delay in triggering transcriptional cascades in the cold response of cassava. Only the ERF9 transcription repressor was down-regulated after 4 h under cold stress. An earlier study had reported that its homologous gene (AT5G44210) of ERF9 also experienced a significant decrease in its expression after exposure to heat and freezing in Arabidopsis [49]. One of CBF homologous gene (CUST_6694, hereafter CBF-like) was found to be up-regulated following 9 h cold treatment. Due to a high P value (0.99), it was excluded from the analysis under our standard (P value ≤ 0.05). Such scenario indicates that the timely response of DREB/CBF pathway to cold might be variable among different species. Indeed, its cold-induced expression at 9 h following cold treatment, together with most of the AP2-EREBP family factors, was verified by real-time RT-PCR ( Figure 4A, B). Therefore, based on our analysis, we speculated that cassava might mount a relatively slow mobilization of stress-related regulators and their downstream genes after exposure to cold, rendering cassava vulnerable to cold stress. Additionally, the transcript of an Inducer of CBF Expression 1 (ICE1)-like gene of cassava was found to be unchanged ( Figure 4B). In Arabidopsis and other plants, ICEs have been reported to be constitutively expressed [50].
A comprehensive expression analysis of MYB TFs demonstrated that almost all MYB TFs are responsive to stresses or hormones [51]. In our study, five MYB or MYB-related family members showed differential expression, with 3 cold-inducible and 2 cold-repressed. Two transcripts, encoding PCL1 with a single myb DNA-binding domain, were required for circadian rhythms [52], and the other three MYB/MYB-related transcription factors were responsive to a diverse of hormone treatments in Arabidopsis [53]. Differential expression of MYB TFs implies that other environmental or hormonal pathways may be involved in cold response in cassava. Furthermore, the up-regulation of the CUST_3593 transcript (PCL1) was also validated by real-time RT-PCR. Members of the GRAS gene family encode transcriptional regulators that have diverse functions in plant growth and development, such as gibberellin signal transduction, root radial patterning, and axillary meristem formation [54]. Although the Arabidopsis genome encodes at least 33 GRAS protein family members, few GRAS proteins have been characterized thus far [55]. It is worth noting that all GRAS family members were down-regulated on the array ( Table 2), suggesting that cold stress inhibit plant growth and conserve energy to adapt to the adverse environment [56]. Similarly, HSF TFs were reported to be involved in the response to heat and other abiotic stresses [38,57]. The microarray analysis also showed the inducible expression of heat shock transcription factor A8 and C1 following cold treatment (Table 2). Interestingly, A8 was also found to be induced in response to low temperatures in potato and Arabidopsis [36]. The induction of heat shock TFs and heat shock proteins (Table 1) revealed that cold and heat stresses may share common responsive elements. The majority of WRKY family TFs is known to be responsive to biotic and/or abiotic stress, despite the fact that most research has focused on the role of these genes in plant-pathogen interactions [58]. WRKY7, known to encode a novel CaM-binding transcription factor of the WRKY group IId in Arabidopsis [59], was also found to have differential expression in our study. This might indicate a crosstalk between cold stress response and plant-pathogen interaction in cassava.
In addition to the above-mentioned TFs, we found two AREBs that were cold-repressed in our study. These factors belonged to group A bZIP transcription factors, which play a role in plant pathogen responses, light signaling, and ABA and abiotic stress signaling, and were induced by salt and ABA in Arabidopsis [60]. The zinc finger family has been demonstrated to be involved in the cold response in Arabidopsis [38]. In our study, we found 3 probes that fell into 2 subfamilies of the zinc finger family (C2H2 and C3H), including 2 up-regulated and 1 down-regulated genes. Zinc finger proteins involve in ROS and abiotic stress signaling in Arabidopsis [61]. The C2H2-type zinc finger protein gene ZAT12-like was induced not only by cold stress, but also by PEG and salt stress, as illustrated by our realtime RT-PCR results ( Figure 4B). This indicates that cassava needs to trigger the expression of ROS-scavenging genes to adapt oxidative stress. The plant-specific NAC transcription factor family has been implicated in plant development processes, such as shoot apical meristem (SAM) maintenance and organ differentiation [62], as well as biotic and abiotic stress responses [63]. We found only one transcript encoding for the NAC domain transcription factor, which was down-regulated on the array. In addition, five novel transcription factor families (ABI3-VP1, JUMONJI, LIM, Sigma70-like, and TCP) were also identified. Their homologous genes in other plant species have not yet been reported in response to cold stress, suggesting that these genes might be specific to cassava and are attractive targets for further functional characterization.

Cold-responsive genes related to 'Chloroplast'
The category of 'Chloroplast' was an abundant and important GO Slim. It accounted for 21.06% of total GOs in 'Cellular Component' (Figure 3C, Additional file 11), indicating that a great number of chloroplast-associated genes changed their expression to adapt to cold. This result suggests that the chloroplast might be one of the major organelles affected, as supported by the observation of damaged chloroplasts in cold-stressed leaves ( Figure 1C). Although the differential expression of chloroplast-related genes might also be affected by other environmental changes on cassava plants during the treatment, i.e. light change, low temperature was the dominant factor in this study.
Chloroplast is an important organelle unique to plant cells and is the site of photosynthesis. In previous reports, cold stress could reduce the rate of photosynthesis by interfering with the function of many photosynthesisrelated proteins [64]. In our study, 43 out of 319 coldresponsive transcripts were assigned to the 'Chloroplast' class, including 22 up-regulated and 21 down-regulated genes ( Table 3). Among them, several genes encoding protein localized to thylakoid were identified. For examples, CAB1 (chlorophyll a/b binding protein 1), a lipidassociated family protein, a thylakoid leumenal 20 kDa protein, and a hypothetical protein (CUST_7076) were all up-regulated under cold stress. In contrast, LHCA4 (photosystem I light harvesting complex gene 4), PsaL (photosystem I subunit L), and the rubredoxin family protein were dramatically down-regulated. Two of them (LHCA4 and PsaL) have also been proposed to be involved in photosynthesis in Arabidopsis. These results suggested that differential expression of chloroplast genes, especially those encoding for thylakoid associated proteins, by cold stress might lead to the processes of thylakoid distortion, chloroplast malfunction, and photosynthesis inhibition.

KEGG pathway analysis of cold-responsive genes
To determine whether the cold stress responsive genes engaged in specific pathways, 291 Arabidopsis AGI loci representing the 319 DEGs were used as objects to search against KEGG pathway maps in Arabidopsis thaliana. Finally, 44 related pathways were identified ( Figure 5, Additional file 12). Several interesting and important pathways, including 'Plant hormone signal transduction' (ath04075), 'Plant-pathogen interaction' (ath04626), 'Phenylalanine metabolism' (ath00360), 'Fatty acid biosynthesis' (ath00061), and 'Starch and sucrose metabolism' (ath00500), were involved and function in the early response to cold stress.
'Plant hormone signal transduction' (ath04075) comprised 8 genes of the 44 cold-regulated and pathway-hit genes on our array. In this pathway, the transcripts of GAI, SAUR (auxin-responsive protein), EBF1 (EIN3binding F box protein 1), EIN4 (ethylene insensitive 4), and AREB3 were determined to be down-regulated, and the transcript of ERF2 was identified as being up-regulated. Their homologous genes have been reported to be implicated in GA, auxin, ethylene (ETH), and ABAmediate hormone signal transduction, respectively [65,66]. In plants, these hormones play crucial roles in a diverse set of developmental processes, as well as biotic and abiotic stresses [67]. For example, ABA accumulates in response to abiotic stresses, such as cold, salt, and drought [68]. Among these hormone-related genes in Arabidopsis, GAI has been also reported to decrease its transcript levels in response to cold [38], and ERF2 has been identified in the response to chitin, a plant-defense elicitor [69]. Previous studies have suggested that the differential regulation of genes involved in hormone signal transduction might play key roles in the early cold stress response [38].
Two other important pathways, including 'Phenylalanine metabolism' (ath00360) and 'Plant-pathogen interaction' (ath04626), were also found in our study to be regulated. PDS1 (phytoene desaturation 1), ASP3 (Aspartate aminotransferase 3), and peroxidase 72, associated with 'phenylalanine metabolism', were all highly accumulated in response to cold. As noted in previous studies, PDS1 has been reported to be engaged in the carotenoid and plastoquinone biosynthetic process [70]. ASP3, also participating in proline metabolism, has been shown to be associated with leaf senescence and to be induced in response to darkness and ethylene-induced artificial senescence [71]. In Arabidopsis and other plants, proline levels are mainly determined by balance of biosynthetic and catabolic pathways, controlled by P5CS (△ 1 -pyrroline-5-carboxylate synthetase) and ProDH (proline dehydrogenase) genes, respectively [72]. Up-regulation of cassava homologous genes P5CS and ASP3, and down-regulation of ProDH, which were confirmed by real-time RT-PCR ( Figure 4C), might correlate with the high accumulation of proline ( Figure 1E). Peroxidase 72, a member of the Class III peroxidases, is involved in the response to oxidation stress and potassium resupply [73]. Additionally, transcripts from MAP Kinase 4, UNE14 (unfertilized embryo sac 14), and TCH2, members of the 'Plant-pathogen interaction' pathway, were also induced by cold (Table 4). MAP Kinase 4 and TCH2 have been shown to be involved in response to several stresses, including cold, salinity, and heat in Arabidopsis [39]; whereas UNE14 also encodes a calcium ion binding protein. These results indicate that calcium ion binding proteins and kinase-mediated signal transduction may play pivotal roles in cold stress response in cassava. In addition, 'Fatty acid biosynthesis' (ath00061) and 'Starch and sucrose metabolism' (ath00500) were also identified through the pathway analysis. FAB1 (fatty acid biosynthesis 1) and SSI2, which were also assigned to 'Chloroplast' category, were involved in 'Fatty acid biosynthesis'. Earlier studies in Arabidopsis indicated that FAB1 is correlated with chilling stress [74], whereas SSI2 has been shown to respond to biotic attacks, including viruses, insects, and bacteria [75]. Both were strongly induced during the entire cold treatment process, indicating that membrane modification may occur during the early cold stress response. In Arabidopsis, TPS7 [alpha, alpha-trehalose-phosphate synthase (UDPforming)/transferase] has been reported to participate in the trehalose biosynthetic process. AMY1 (alpha-amylase-like 1) has also been shown to take part in the degradation of starch to sugar and to be induced by biotic and abiotic stress [76]. In our study, TPS7 and AMY1 were identified as being involved in 'Starch and sucrose metabolism'. Their cold-induced expression was confirmed by real-time RT-PCR analysis ( Figure 4D). In cold-stressed leaves, starch granules were found to have disappeared ( Figure 1C). Even though their total sugar content remained unchanged after 4 h and was lower after 9 h of treatment, an obvious increase was detected after 24 h ( Figure 6A). Interestingly, sucrose was significantly accumulated and reached the highest level in 9 h; whereas glucose levels were only slightly increased under cold stress ( Figure 6B). These results suggested that soluble sugars might act as both osmolytes and signal molecules in the cold response of cassava, similar to other plant species [77]. It also provided insight into the hypothesis that starch degradation provides conserve energy for plants under adverse conditions. Such transition from synthesis metabolism to degradation metabolism would also explain the reduced productivity observed in plants under stress conditions.
In short, a variety of interesting and important pathways are involved and function in the cellular response to cold stress. Despite many pathways, e.g., 'Plant hormone signal transduction' (ath04075), 'Fatty acid biosynthesis' (ath00061), and 'Starch and sucrose metabolism' (ath00500), have been well-demonstrated in other plant species, regulation of different pathways or gene has been observed, illustrating that a complex and specific network is involved in the early cold response in cassava.

Validation of microarray results by real-time RT-PCR
To validate our microarray data, real-time RT-PCR analysis was performed on 18 selected differentially expressed transcripts, including 13 up-regulated and 5 down-regulated genes (Table 4). These genes belong to divergent functional categories or pathways. For instance, four genes [gibberellin 2-beta-dioxygenase, unknown protein (CUST_16591), ubiquitin-protein ligase, and GSTU19] were involved in 'Response to abiotic and biotic stimulus' and 'Response to stress'. One gene (protein phosphatase 2C) was implicated in 'Protein metabolism,' and other three genes [AP2 domain transcription factor (CUST_20765), C2H2-type zinc finger (ZAT12-like), and PCL1 (CUST_3593)] were grouped as transcription factors. In addition, 3 genes (ubiquitin-protein ligase, TCH2, and S-adenosylmethionine decarboxylase) were assigned to 'Circadian rhythmplant' (ath04712), 'Plant-pathogen interaction', and 'Arginine and proline metabolism' (ath00330), respectively. Seven genes that are unclassified in the GO Slims or pathways but displayed significant changes in transcript abundance were also tested. For these genes, the fold change (-△△Ct) measured by real-time RT-PCR and by microarray was highly consistent (72.2% and 94.4% for 4 h/0 h and 9 h/0 h, respectively). Their expression kinetics from the real-time RT-PCR results was similar to those of the microarray analysis. These results re-confirmed the accuracy of our microarray data.

Genes responsive to different stresses and their tissuespecific expressions
The cold stress-signaling pathway may interact with other signaling systems of, for example, ABA, salt, and drought [19]. The expression patterns of the 13 coldresponsive genes were analyzed by real-time RT-PCR using cassava in vitro shoot cultures treated with 100 μM ABA, 25% PEG, and 250 mM NaCl for 6 h. Genes with fold changes larger than 4 were defined as strongly responsive to the stresses. Several cold-responsive genes were also strongly induced by ABA, PEG, or salt stress (e.g., CUST_16591, novel protein with unknown function; CUST_12461, C2H2-type zinc finger, ZAT12-like; CUST_4383, gibberellin 2-beta-dioxygenase) ( Table 5). These results indicated that dynamic crosstalk exists among signaling pathways related to cold, drought, and salt in cassava because all of these stresses cause cellular  dehydration [20]. Genes uniquely responsive to cold stress drew more attention in our study, such as TCH2 (CUST_17985), AP2 domain transcription factor (CUST_20765), and JAZ7 (jasmonate-ZIM-domain protein 7, CUST_9007) ( Table 5). Both TCH2 and JAZ7 were involved in the 'Plant-pathogen interaction' in Arabidopsis (Additional file 12) and rice, suggesting that chilling-induced processes share some common features with the defense mechanism against pathogens in plants [58]. In addition, protein transmembrane transporter (CUST_13673), which was strongly induced by ABA and salt stress, was also believed to function in early signal reception.
To further analyze the expression patterns of these stress-regulated genes in different tissues, real-time RT-PCR analysis was conducted. mRNAs were extracted from the apical buds (AP), fibrous roots (FR), mature leaves (ML), stem cambia (SC), young leaves (YL), and young stems (YS) of 3-month-old cassava plants grown in a greenhouse. Interestingly, 7 tested genes were primarily expressed in FR, ML, SC, and YS rather than in AP and YL (Table 6). This tissue-specific expression of cold-responsive genes might explain why the apical shoots, including AP and YL, were more susceptible to cold stress than ML, but this difference remains to be further determined.

Changes in H 2 O 2 content and ROS scavenging enzyme activities
As previous reported, chilling conditions may lead to an accumulation of ROS such as hydrogen peroxide (H 2 O 2 ) and superoxide radical (O 2 -) [78]. To analyze and visualize H 2 O 2 produced in the cassava leaves subjected to cold stress (7°C), diaminobenzidine (DAB) was used to  stain the first or second fully expanded leaves of in vitro specimens and greenhouse-grown cassava plants. In unstressed leaves, few DAB polymerizations could be observed; in contrast, there was a significant increase in DAB polymerization detected as early as 4 h following cold treatment that continued up to 24 h ( Figure 7A). This result was consistent with the quantification of H 2 O 2 content ( Figure 7B), indicating that the oxidative stress may exert a toxic effect on cassava to adapt or survive under cold stress conditions. Plants, as well as other organisms, have evolved antioxidant systems to protect themselves against toxic species of oxygen. ROS scavenging enzymes, including catalase (CAT), superoxide dismutase (SOD), and glutathione transferase (GST), have been demonstrated to play key roles in the removal of ROS. In our DEGs, several genes encoding ROS scavenging enzymes were recorded to be target to chloroplast, including CAT2 and GSTU19 (Table 3). In addition, peroxidase 72, and other three GSTs (GSTU7, GSTU1, and microsomal GST) were also identified (Additional file 4). The upregulation of four ROS-scavenging transcripts (CAT2, peroxidase 72, GST1 and microsomal GST) was verified by real-time RT-PCR ( Figure 7C). In Arabidopsis, CAT2 was reported to be induced by cold stress and GSTU19 was responsive to many biotic and abiotic stresses, such as water deprivation and oxidative stress [79]. To further validate the function of ROS scavenging enzymes after exposure to cold stress, the activities of CAT and SOD in cassava leaves subjected to 7°C were measured. Consistent with the microarray data, an increase in CAT activity was observed after 9 h of cold treatment ( Figure  7D). Similarly, a significant elevation in SOD activity was also observed after 4 h of cold stress, and such increases continued for 24 h after the initiation of stress treatment ( Figure 7E). These results including induction of ZAT12-like gene ( Figure 4B) suggested that ROS signaling pathways, including ROS scavenging enzymes, were involved in the ROS detoxification induced by the cold stress response in cassava. Similarly, rice and maize, as chilling-sensitive tropical crops, could only withstand transient and milder cold stress. In japonica rice, an oxidative-mediated network has been proposed to play a key role in the early response to chilling stress and short-term defenses [80], and insufficient antioxidant defenses are thought to cause maize chilling sensitivity [81]. Generally, the ultra-structural changes of the chloroplast ( Figure 1C) and oxidative burst were typically characteristic of PCD in plants, indicating that PCD may be a part of an adaptive mechanism to survive stress [82].

Conclusions
Our study presented a genome-wide gene expression profiling of cassava subjected to cold stress using the microarray technology. Overall, the transcriptomic responses of cassava to cold stress were basically consistent with the changes seen in other plants under abiotic stresses [83]. A considerable amount of specific cassava genes related to different biological functions were identified. For examples, many new stress-responsive genes or TFs, such as early flower 4, ERF2, and AP2 domain transcription factor (CUST_2332), were found, suggesting that various regulatory pathways may exist in cassava together with the well-characterized CBF pathway. Based on the comprehensive transcriptomic, real-time RT-PCR and physiological analyses, our study supports the fact that crosstalks among different signaling systems play an important role in regulating the cold stress response in tropical plants. Thus, a hypothetical model for depicting the components involved in cold response networks in cassava could be established (Figure 8). Plants may perceive low temperature through cell membrane receptors, most likely transmembrane proteins, and by membrane modification through fatty acid synthesis (e.g., FAB1 and SSI2). Then, the cold signal transduction might induce cellular metabolism changes, such as the calcium signaling cascade, hormone signal transduction, and ROS signaling. The calcium signaling pathway may stimulate the calcium ion binding protein,  leading to the activation of cascade kinase activities (e.g., MAP kinase 4), which could switch on various cold stress-responsive genes and transcription factor family proteins (e.g., AP2-EREBP, HSF, and GRAS). Hormone (auxin, ABA, GA, ETH, and JA) signal transduction could also be activated to alter plant growth status in order to adapt to stress condition through the differential regulation of their downstream genes (e.g., SAUR, AREB3, GAI, ERF2, and JAZ7, respectively). Moreover, other metabolisms could also take part in the process, Figure 8 Molecular model of the early cold response in cassava. Two biological processes, namely stress perception and the physiological response, are illustrated. After signal reception, stress-activated Ca 2+ signaling, ROS signaling, and hormone signaling modulate the expression of stress-responsive genes, which include metabolic enzymes, transcription factors, kinases, and ion transporters. The physiological changes that manifested as membrane modification, chloroplast malfunction, and starch and sucrose metabolism, as well as amino acid metabolism, cause cassava to either have increased cold tolerance or to enter into accelerated programmed cell death. including starch and sucrose metabolism (e.g., TPS7 and AMY1), amino acid metabolism (e.g., ASP3 and PDS1), and the ROS scavenging system (e.g., catalase 2, SOD, and GSTs). Apparently, the balance between cell damage and tolerance might decide the fate for coldstressed cassava plants (Figure 8). On the other hand, a large number of genes (about 37.2% of total 508 DEGs) encode proteins of unknown functions. Studying these genes may reveal novel mechanisms that are fundamental to the ability of cassava to cope with cold stress.
In summary, our array study will provide the fundamental knowledge related to the biological and physiological changes of cassava under cold stress. It will also be served as a very useful genetic resource for relevant research community globally. Further studies are needed to verify the functions of candidate genes for improving cassava tolerance ability through genetic engineering.

Plant materials used for microarray hybridization
In vitro plants of a cassava cultivar (TMS60444) were planted in plastic pots at 28°C under a 16 h light photoperiod (110-150 μmol·m -2 ·s -1 ) for 3 months in the greenhouse. Plants with a uniform growth status were transferred to a chamber for cold treatment at 7°C under weak light (cool-white fluorescent light at approximately 35 μmol·m -2 ·s -1 ). The apical shoots (about 4-6 cm in size, Additional file 1), including apical buds, stems, and the first and second expanded leaves, were harvested after 4 h and 9 h of treatment, then frozen in liquid nitrogen and held at -80°C prior to transporting to Shanghai Biochip Corporation (Shanghai, China) for RNA extraction. Untreated apical shoots were harvested as controls (0 h). More than 3 plants were harvested and pooled for each time point, and the collection was repeated 3 times as biological replicates.

Physiological analyses of cold-treated cassava plants
To analyze the physiological changes of cassava under cold treatment prior to the microarray study, malondialdehyde (MDA) and proline, indicators of the cold response in plants, were measured. The MDA content was determined by the thiobarbituric acid (TBA) reaction with minor modifications of the method of Dhindsa et al [84] Proline concentrations were measured by the sulfosalicylic acid-acid ninhydrin method according to Bates et al [85] with slight modifications.

Determination of hydrogen peroxide (H 2 O 2 ) content and diaminobenzidine (DAB) staining
Hydrogen peroxide (H 2 O 2 ) levels were determined according to Velikova et al [86] The first and second expanded leaves (1 g) were homogenized in an ice bath with 10 mL of 0.1% (w: v) trichloroacetic acid (TCA).
The homogenate was centrifuged at 12,000 × g for 15 min, and 1 mL of the supernatant was added to 1 mL of 10 mM potassium phosphate buffer (pH7.0) and 2 mL of 1 M KI. The reaction solution containing the supernatant was read at 390 nm using a UV-vis spectrophotometer. The content of H 2 O 2 was given on a standard curve.
H 2 O 2 was visualized by staining with DAB, which undergoes a polymerization reaction to yield a darkbrown color once it encounters H 2 O 2 [87]. The first or second fully developed leaves harvesting from the in vitro and greenhouse-grown cassava plants, which were treated with cold stress (7°C) for 0, 4, 9, and 24 h, were excised. Their leaf petiole cuttings were immersed in 1 mg/mL of DAB solution (pH3.8) for 8 h under continuous light. The leaves were immersed in 96% (w/v) boiling ethanol for 10 min to decolorize the chloroplasts. After cooling, the leaves were kept in ethanol and photographed.

Assay for catalase (CAT) and superoxide dismutase (SOD) activities
For the measurement of enzyme activities, the first and second fully expanded leaves (1 g) from 3-month-old greenhouse plants under cold stress for 0, 4, 9, and 24 h were homogenized in 5 mL of 10 mM potassium phosphate buffer (pH7.0) containing 4% (w: v) polyvinylpyrrolidon (Mr 25 000). The homogenate was centrifuged at 10,000 rpm for 15 min, and the supernatant obtained was used as the enzyme extract. All steps in the preparation of the enzyme extract were carried out at 0-4°C . CAT activity was determined by directly measuring the decomposition of H 2 O 2 at 240 nm as described by Cheng et al [88] The reaction mixture contained 50 mM potassium phosphate buffer (pH7.0), 10 mM H 2 O 2 , and 200 μL of enzyme extract in a 2 mL volume. SOD activity assay was based on the method described by Beauchamp and Beaucham et al [89], which measures the inhibition of the photochemical reduction of nitro blue tetrozulium (NBT) at 560 nm. Three milliliters of reaction mixture contained 50 mM phosphate buffer (pH7.8), 0.1 mM EDTA, 13 mM methionine, 75 μM NBT, 16.7 μM riboflavin, and 300 μL of enzyme extract.

Measurement of total soluble sugars, sucrose and glucose
Soluble sugars were extracted from the leaf tissues by a hot ethanol extraction as follows. Leaves were sampled from 3-month-old cassava, which was treated by cold stress (7°C) at 0, 4, 9, and 24 h and quickly dried at 110°C for 15 min and then 70°C overnight. Exactly 100 mg of homogenized dry leaf powder for each sample was extracted with 8 mL of 80% ethanol (v/v) at 85°C for 40 min. The extracts were then centrifuged at 12,000 × g for 10 min. The ethanol extraction step was repeated two times. The 3 resulting supernatants were combined and diluted with 80% ethanol to a volume of 50 mL; they were then treated with 20 mg of activated charcoal at 80°C for 30 minutes. The soluble sugar analysis was conducted according to Ebell [90] with minor modifications. After removing the activated charcoal with a 0.2 μm filter, an aliquot of extraction buffer was reacted with anthrone-sulfuric acid at 95°C for 15 min and then cooled to room temperature. The absorbance of the reaction solution was read at 620 nm. The total soluble sugar content was calculated by a standard curve. Sucrose and glucose content were measured by glucose oxidase method as described by Johnson et al [91] with slightly modification.

Oligonucleotide microarray preparation, hybridization, and data extraction
The custom-designed 60-mer microarray [31] was constructed based on public sequence information from a large collection of cassava ESTs from NCBI (71,520 ESTs, released 28 March 2008) and TIGR (Manihot_es-culenta_release_5, containing 5,189 assemblies and 10,214 singletons, released 1 June 2007), as well as a 35,400 full-length cDNA RIKEN library [27]. Total RNA was extracted from both control and cold-treated tissue samples using an RNeasy Mini Kit (Qiagen, Valencia, CA, USA). The RNA quality was checked on a 1.2% agarose gel using an RNase-free electrophoresis system. RNA labeling and hybridization were conducted by the Shanghai Biochip Corporation (Shanghai, China) following the manufacturer's instructions. The arrays were incubated at 65°C for 17 h in Agilent hybridization chambers (G2545A) and then washed according to the protocol at room temperature. Hybridized microarray slides were scanned at 5 μm resolution with an Agilent Technologies Scanner (G2505B), and the images were saved in JPG format. Both the 10% and 100% Photomultiplier tube (PMT) settings were selected, and the combined images were exported. The signal intensities of all spots on each image were quantified with the Feature Extraction software (Agilent Technologies, Santa Clara, California, USA), and the data were saved as .txt files for further analysis.

Data normalization and identification of differentially expressed genes (DEGs)
The signal intensity of each gene was globally normalized within the GeneSpring GX Software (Agilent Technologies, Santa Clara, California, USA) following the workflow guide [92]. The signal-to-noise ratio (SNR) was calculated as follows: the difference of the median signal minus the background median signal, divided by the background standard deviation [22]. Pair comparison was used to analyze the normalized and averaged data from the three types of samples (0, 4, and 9 h). The P values were calculated using a t test, and the fold changes between each comparison for each gene were compared. The genes that were induced or suppressed at levels equal to or greater than a twofold ratio (twofold change cutoff, FCC) were taken to be differentially expressed when SNR > 2.6 and P value ≤ 0.05.

Microarray data analysis
All non-redundant sequences that were considered to be unique cassava genes were locally blasted in the TAIR protein database (27,217 Arabidopsis protein sequences), which was download from ftp://ftp.arabidopsis.org/ home/TAIR, using the blastx program in the blastall package (version 2.2.9). The top hits were used for gene annotations, and the corresponding Arabidopsis gene locus identifiers were mapped to the GO function annotation (http://www.arabiopsis.org/) and to the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The TAIR percent of cold-responsive genes was calculated as follows: The TAIR percent (%) = the number of genes annotated to terms in the GO Slim category divided by N × 100, where N represents the total number of genes from the input list annotated to any term in this ontology.
Utilization of the microarray MIAME information about the cassava transcriptome microarray used here has been deposited in the Gene Expression Omnibus (GEO) of NCBI. The accession numbers are: Platform, GPL11271; Series, GSE31073; Samples, GSM769563-GSM769571.

RNA preparation and real-time RT-PCR
To verify the expression patterns of DEGs identified from the microarray analysis, real-time RT-PCR was conducted using plant materials from either in vitro shoot cultures or greenhouse-grown plants. Four-weekold in vitro seedlings were transferred into plastic jars containing 150 mL of MS solution [1 × Murashige and Skoog basal salts (MS, Duchefa), 1% sucrose, pH5.7] in a chamber at 25°C, 70% relative humidity, and a light intensity of 125 μmol·m -2 · s -1 on a 16 h light/8 h dark cycle for 4 days. The solution was then replaced with fresh MS medium (pH5.7) supplemented with 100 μM ABA, 25% PEG, or 200 mM NaCl. Meanwhile, the seedlings were directly transferred to a 7°C chamber for cold stress treatment. After 6 h, whole seedlings were harvested and frozen in liquid nitrogen; they were then stored at -80°C for abiotic stress analysis. Plants in the same hydroponic system with no stress treatment were harvested and used as controls. The apical buds (AP), fibrous roots (FR), mature leaves (ML), stem cambia