Transcriptional profiling of Petunia seedlings reveals candidate regulators of the cold stress response

Petunias are important ornamentals with the capacity for cold acclimation. So far, there is limited information concerning gene regulation and signaling pathways associated with the cold stress response in petunias. A custom-designed petunia microarray representing 24816 genes was used to perform transcriptome profiling in petunia seedlings subjected to cold at 2°C for 0.5 h, 2 h, 24 h, and 5 d. A total of 2071 transcripts displayed differential expression patterns under cold stress, of which 1149 were up-regulated and 922 were down-regulated. Gene ontology enrichment analysis demarcated related biological processes, suggesting a possible link between flavonoid metabolism and plant adaptation to low temperatures. Many novel stress-responsive regulators were revealed, suggesting that diverse regulatory pathways may exist in petunias in addition to the well-characterized CBF pathway. The expression changes of selected genes under cold and other abiotic stress conditions were confirmed by real-time RT-PCR. Furthermore, weighted gene co-expression network analysis divided the petunia genes on the array into 65 modules that showed high co-expression and identified stress-specific hub genes with high connectivity. Our identification of these transcriptional responses and groups of differentially expressed regulators will facilitate the functional dissection of the molecular mechanism in petunias responding to environment stresses and extend our ability to improve cold tolerance in plants.


INTRODUCTION
The petunia (Petunia hybrida) is an important ornamental that is widely used in landscaping around the world. It is also an excellent model system for studying the evolutionary development of floral organs and the biosynthesis of flavonoids. Nevertheless, because the plant originated from South America, petunias are native to warm habitats. Late-spring freezes can cause substantial damage to petunias as they are less tolerant of cold conditions and lack rapid acclimation response mechanisms as compared with species tolerant of ice in their tissues, such as the spinach and Brassicas (Pennycooke et al., 2003). Thus, frozen conditions and low temperatures are significant limiting factors for its productivity, growing period and geographical location, and there is a demand for novel petunia cultivars with improved cold tolerance.
Cold stress, including chilling (<20 • C) and/or freezing (<0 • C) temperatures, adversely affects plant growth and development and causes a wide range of physiological and biochemical responses in plants. These responses include transient increases in hormone levels (e.g., ABA) (Chandler and Robertson, 1994), changes in the compositions of membrane lipid (Lynch and Steponkus, 1987), the accumulation of compatible osmolytes, such as soluble sugars, betaine, and praline (Koster and Lynch, 1992;Nomura et al., 1995;Doerffling et al., 1997), and increases in the level of antioxidants (Pennycooke et al., 2005). However, tolerance to freezing can be enhanced when plants are prior exposed to a period of low, but non-freezing temperatures, which is called cold acclimation (CA) (Weiser, 1970). CA activates the expression of a diverse array of plant genes (Thomashow, 2001). In the past decades, numerous cold-responsive genes have been identified and their functions have been characterized using gene transfer, resulting in increased plant stress tolerance (Sanghera et al., 2011). To date, the mechanisms of CA have been extensively investigated in the model plant Arabidopsis thaliana (Gilmour et al., 1998), as well as in other crucial crop species such as barley and maize (Stanca et al., 1996;Fernandes et al., 2008). The most deeply studied mechanism involves a class of ethylene response factors (ERFs), that is, CBFs (C-repeat binding factors, also called DREB1s or dehydration-responsive elementbinding protein 1s), which can bind to cis-elements in the promoters of their downstream target genes to execute a highly coordinated transcriptional response to low temperature signals (Chinnusamy et al., 2007). In contrast, less information is available for petunias.
With the advancement of molecular biology technologies and "omics" tools, microarrays have been devised as a powerful strategy for the analysis of transcript abundance in plants. In recent years, transcriptional profiling of diverse plant species under defined stimuli such as hormone or abiotic stresses using microarray analysis has been particularly informative. For example, the expression patterns of genes were identified in rice and Arabidopsis under abscisic acid or cold, drought and high-salinity stresses (Seki et al., 2002;Rabbani et al., 2003). The transcriptomic identification of candidate genes involved in response to salt stress in cotton has also been reported (Rodriguez-Uribe et al., 2011). Additionally, transcriptome profiling of cassava treated with low temperatures highlighted the dynamic responses of tropical plants to cold stress . To date, several studies have investigated petunias using genomic tools to develop novel technologies and knowledge. For example, a cDNA library was constructed from petunias for the isolation of genes responsible for benzylbenzoate formation (Boatright et al., 2004). The capacity of a targeted transcriptomics approach to disambiguate regulation of a particular biosynthetic process was illustrated by the identification of the MYB transcription factor ODORANT1, a key regulator of floral scent biosynthesis, using a dedicated petunia microarray (Verdonk et al., 2005). Analysis based on a petunia small RNA library identified unusual tissue-specific expression patterns of conserved miRNAs and a 24-mer RNA (Tedder et al., 2009). Moreover, the petunia transcriptome has recently been employed to gain insight into potential shifts in symbiotic signaling, metabolic activity, and the defense status of plants exposed to high P levels (Breuillin et al., 2010). Nevertheless, genomewide expression profiling data for petunias under cold stress is still deficient, thereby limiting our capacity to decipher the molecular mechanisms associated with different stress responses. To address this need for a popular ornamental, comprehensive transcriptomic data describing the response of petunias to cold stress treatment were obtained and analyzed using a custom Roche Nimblegen petunia microarray. The results revealed a set of candidate regulators. Our study provides a first glimpse into the understanding of the gene regulation pathways involved in the adaptation of this plant species to cold stress, as well as the potential applications of these data in systems biology and molecular physiology.

PLANT MATERIALS USED FOR MICROARRAY HYBRIDIZATION
More than 20 inbred lines were established by successive generations of selfing. To compare the differences of cold resistance among those inbred lines, they were planted in open field in winter under natural low temperatures as well as in artificial climate chambers. One line (coded H) with the best cold tolerance was screened out.
In vitro seedlings of line H were planted in plastic pots at 25 • C under a 14 h photoperiod with a light intensity of 2000-2500 lux provided by cool-white fluorescent lights for 1 month in the tissue culture room. Plants at 4-5 pairs of true leaf stage with a uniform growth status were transferred to a chamber for cold treatment at 2 • C under weak light (cool-white fluorescent light at approximately 500-1000 lux). Whole plants were harvested after 0.5 h, 2 h, 24 h, and 5 d of treatment, then frozen in liquid nitrogen and held at −80 • C prior to transportation to the Bioassay Laboratory of CapitalBio Corporation (Beijing, China) for RNA extraction. Untreated whole plants were harvested as controls (0 h). More than three plants were harvested and pooled for each time point, and the collection was repeated three times as biological replicates.

MICROARRAY PREPARATION, HYBRIDIZATION, AND DATA EXTRACTION
The custom-designed four-plex microarray with 72,000 features was constructed based on the EST sequences generated in Breuillin's et al. (2010) study together with all of the P. hybrida and P. axillaris sequences available from Genbank, TIGR, and the Solanaceae genomics network SGN. The compiled sequence set was ultimately composed of 24,816 non-redundant unique sequences (unigenes). Total RNA was extracted from both the control and cold-treated samples with Trizol Reagent (Invitrogen, http://www.invitrogen.com). 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 Bioassay Laboratory of CapitalBio Corporation (Beijing, China) following the manufacturer's instructions. Total RNA was used to generate Cy5/Cy3-labeled DNA using a cRNA amplification labeling kit (CapitalBio). Briefly, the 1st-strand cDNA was reverse transcribed with the T7 Oligo (dT) primer, then the 2nd-strand was synthesized and purified as the template for cRNA transcription with the T7 Enzyme Mix. Two micrograms of cRNA was reverse transcribed to cDNA with random primers, and the cDNA was used as the template for Cy5/Cy3labeled DNA synthesis with Cy5/Cy3-dCTP. The Cy5/Cy3-labeled DNA in 80 μL hybridization solution (3 × SSC, 0.2% SDS, 5 × Denhart's, and 25% formamide) was hybridized to the 4 × 72K Petunia Array (Roche Nimblegen). After hybridization at 42 • C overnight, the samples were washed with 2 × SSC (plus 0.2% SDS) and 0.2 × SSC for 5 min each. The GeneChips were used for scanning with the LuxScan 10 KA (CapitalBio). Three biological replicates were included in this study.

DATA COLLECTION AND IDENTIFICATION OF DIFFERENTIALLY EXPRESSED GENES
The probe intensity files resulting from RNA hybridization were first loaded into R. Then, background correction, quantile normalization and summarization were performed using the Robust Multi-array Average (RMA) method included in the Bioconductor Affy package (Bolstad et al., 2003;Irizarry et al., 2003;Gautier et al., 2004;Gentleman et al., 2004) based on the homemade chip description library. The correlation coefficient of paired samples was calculated to test the repeatability between samples. Then, euclidean distances were calculated, and the clustering tree was constructed using the hierarchical cluster method of "complete linkage clustering." Differentially expressed genes were calculated using the Bioconductor RankProd package function. RP. Rank product is a non-parametric method returning up-and down-regulated genes, their fold changes (FC), p-values and the percentages of false predictions (PFP). This method was shown to perform better than other methods, including significance analysis of microarrays (SAM), Fisher's Inverse χ 2 test and t-based hierarchical modeling (Hong and Breitling, 2008). A gene was considered to be up-or down-regulated if the p-value of the RankProd analysis was <0.05 and the fold change of average expression was >2 (logfold change >1 for up-regulated genes and log-fold change < −1 for down-regulated genes).

FUNCTIONAL ENRICHMENT ANALYSIS
All probe set sequences were annotated by Blast2GO by searching against the nr database using Blastx with parameter E = 10e-3. Then, all probe sets were annotated with three principal GO categories: biological process, molecular function and cellular component. We used the terms in the biological process category for GO analysis. Finally, GO enrichment analysis was performed using a weighted method in combination with Fisher's exact test with a P-value of 0.05 provided by the Bioconductor topGO package (Alexa et al., 2006).

RNA PREPARATION AND REAL-TIME RT-PCR
To validate the expression patterns of the differentially expressed genes identified from the microarray analysis, real-time RT-PCR was performed using plant materials from in vitro shoot cultures. In vitro seedlings with 4-5 pairs of true leaves were transferred into plastic jars containing 100 mL of MS solution (1 × Murashige and Skoog basal salts and 3% sucrose, pH 6.0) in a chamber at 25 • C, 70% relative humidity, and a light intensity of 2000-2500 lux on a 14 h light/10 h dark cycle for 5 days. Then, the solution was replaced with fresh MS medium (pH 6.0) supplemented with 200 mM NaCl, 400 mM mannitol, 100 μM ABA, or 100 μM MeJA. For cold stress treatment, the seedlings were directly transferred into a 2 • C chamber. For dehydration stress treatment, the seedlings were transferred and dehydrated onto dried filter paper in covered plastic dishes. All abiotic stress samplings were harvested from whole seedlings and frozen in liquid nitrogen after 3, 6, and 12 h and stored at −80 • C for abiotic stress analysis. Plants in the same tissue culture system with no stress treatment were collected and used as controls. All of the stress treatments included three biological and temporal replicates.
The expression of cold-responsive genes was validated by realtime RT-PCR using RNA samples extracted with the EASYspin Plant RNA Mini kit according to the manufacturer's protocol (Aidlab, Beijing, China). The RNA concentration and integrity were estimated with the NanoDrop 2000 UV-Vis Spectrophotometer (Thermo scientific, USA) and denaturing agarose gel electrophoresis. Reverse transcription was performed using the Prime Script® RT reagent kit with gDNA Eraser (Perfect Real Time) (TaKaRa, Dalian, China) according to the manufacturer's protocol. Gene-specific primers (Table S1) were designed with ABI primer express 3.0 software (Applied Biosystems, USA). The product size from each amplicon ranged from 105 to 215 base pairs. The real-time RT-PCR was performed with SYBR® Premix Ex Taq™ (Perfect Real Time) (TaKaRa, Dalian, China) as described in the users' instructions for the ABI 7500 Fast Real-Time PCR System (Applied Biosystems, USA). The PCR reaction mixture contained 1 μl cDNA template, 5 μl SYBR®Premix Ex Taq™ (2x), 0.2 μl forward and reverse primer (10 μmol/μl each), 0.2 μl ROX Reference Dye II (50x), and water up to 10 μl. The amplification conditions were as follows: an initial incubation at 95 • C for 30 s, followed by 40 cycles of 95 • C for 3 s, 60 • C for 30 s, and 72 • C for 20 s. A melting curve (60 • -95 • C with a heating rate of 0.05 • C·s −1 and a continuous fluorescence measurement) was performed for each qRT-PCR reaction to determine PCR performances. The glyceraldehyde-3-phosphatedehydrogenase gene (GAPDH) was used as the internal control. All of the samples were measured in triplicate, and the experiments were performed on three biological replicates. The comparative Ct method was adopted to calculate the relative gene expression level across the samples. The relative expression level of each gene in one sample ( Ct) was calculated as follows: Ct target gene -Ct GAPDH. The relative expression of each gene in two diverse samples ( Ct) was calculated as follows: Ct (sample 1) -Ct (sample 2).

CO-EXPRESSION NETWORK ANALYSIS
To identify co-expression of the petunia genes on the array, log-transformed gene expression values for each condition were used in a meta-analysis to perform Weighted Gene Co-expression Network Analysis (WGCNA) (Langfelder and Horvath, 2008). Briefly, WGCNA can be used for finding clusters (modules) of highly correlated genes, for summarizing such clusters using the module eigengene or an intramodular hub gene, for calculating module membership measures, and for relating modules to one another and to external sample traits (using eigengene network methodology). The blockwiseModules function of the WGCNA package in R was used to generate the modules with a power of six. The correlation coefficient relationship of the module eigengene was calculated, then euclidean distances were calculated and the clustering tree was constructed using the hierarchical cluster method of "complete linkage clustering" included with the R package.

RESULTS
To explore the transcriptomic changes in petunias in response to cold stress, a custom four-plex microarray with 72,000 features generated by Roche NimbleGen was employed in the current study. The detailed protocol for designing the petunia microarray has been described by Breuillin et al. (2010). Our microarray data have been deposited to Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/query/ acc.cgi?acc=GSE64963) and the accession number is GSE64963. The results of hierarchical clustering indicated the clear separation of the samples from various time points (Figure 1A), which suggested that the whole experiment from sample collection to data extraction was reliable and reproducible.

TRANSCRIPTOMIC RESPONSES OF PETUNIA SEEDLINGS TO COLD TREATMENT
Two filtering criteria were used in our data analysis to define differentially expressed genes (DEGs): a two-fold or greater change in transcript levels between any two time points and a P-value <0.05. To analyze the differences and similarities among the coldresponsive transcriptomes, a hierarchical clustering approach was utilized to depict the transcripts of all of the differentially expressed probes of the three replicates at 0 h, 0.5 h, 2 h, 24 h, and 5 d. These results demonstrated distinct differences in gene expression profiles between the treatments of 5 d and 0 h or 24 h and 0 h that were in contrast to the relatively highly similar expression profiles between 2 h and 0 h or 0.5 h and 0 h ( Figure 1A). Among the DEGs, 198, 357, 682, and 675 genes were considered to be up-regulated and 53, 17, 416, and 622 genes were considered to be down-regulated based on the analysis of  (Table S2).
Notably, compared with only 251 and 374 genes that showed differential expression at 0.5 h/0 h and 2 h/0 h, 1098, and 1297 genes were found to be differentially expressed at 24 h/0 h and 5 d/0 h, suggesting that prolonged stress treatment triggered the expression of more stress-related genes. To identify both unique and common genes showing differential expression patterns at all time points compared with 0 h, the numbers were calculated and presented using a Venn diagram ( Figure 1B). The results showed that 75 DEGs were commonly induced in all of the comparisons for each time point compared with 0 h, which demonstrated a strong linkage among the four stressed comparison points and a progressive biological process. In contrast, 42, 102, 212, and 268 DEGs exhibiting up-regulated expression were unique at 0.5 h/0 h, 2 h/0 h, 24 h/0 h, and 5 d/0 h; thus, these genes may play special roles at those time points. Only 2 suppressed DEGs were found to be common to all of the comparisons for each time point compared with 0 h, whereas approximately 81% (43/53), 29% (5/17), 58% (243/416), and 73% (455/622) of DEGs exhibiting downregulated expression were unique at 0.5 h/0 h, 2 h/0 h, 24 h/0 h, and 5 d/0 h, respectively.
In total, 2071 DEGs, including 1149 cold-inducible and 922 cold-repressed genes, were identified using our array. These genes represent approximately 8.3% (2071/24,816) of the total expressed petunia transcripts on the array. This result indicated that Petunia hybrida is a species with the capacity for cold acclimation (Yelenosky and Guy, 1989) that possess a similarly large amount of cold-responsive genes similar to that reported for other plants (e.g., Arabidopsis). Indeed, the petunia arrays employed here included probes for only a portion of the total genes. Thus, if the probes on the array covered the whole genome, many more genes responsive to cold would have been found. Because differences in the ability of diverse plant species to tolerate cold might not be totally due to the amount of genes responsive to cold stress (Carvallo et al., 2011), it has been suggested that plants require to duly and coordinately mobilize the biological functions and regulatory networks of stress-responsive genes.

COMPARISON OF SIGNIFICANT GENE ONTOLOGY TERMS OF DEG SETS AT ALL TIME POINTS
To gain insight into the putative functions of the 2071 coldresponsive genes identified by the microarray analysis, gene ontology annotation for all of the differentially expressed genes was performed using Blast2gGO against the local nr database. Accordingly, 1388 annotations were obtained based on the gene enrichment analysis that was performed for each comparison group between all time points and 0 h. All three principal GO categories ("Biological Process," "Molecular Function," and "Cellular Component") were collected for each gene, and the biological process was selected for functional enrichment analysis.
We found 198 unique GO terms enriched by DEGs when analyzed gene sets that were up-or down-regulated separately at each time point (Table S3). At 0.5 h/0 h, the most significant GO terms found in the up-regulated gene set were regulation of RNA biosynthetic process (GO:2001141), regulation of RNA metabolic process (GO:0051252) and regulation of transcription, DNA-templated (GO:0006355). The most significant GO terms in the down-regulated sets were transition metal ion transport, and especially zinc ion transmembrane transport (e.g., GO:0000041 and GO:0071577). At 2 h/0 h, only the upregulated gene sets could be enriched for GO analysis. These DEGs were mainly involved in positive regulation of nucleobasecontaining compound metabolic process (GO:0045935), positive regulation of gene expression (GO:0010628), positive regulation of RNA metabolic process (GO:0051254) and positive regulation of transcription, DNA-templated (GO:0045893). The functional enrichment analysis for the down-regulated gene sets was not available, partially due to the lack of annotation for a number of petunia genes. At 24 h/0 h, the terms, "flavone metabolic process" (GO:0051552), "pigment biosynthetic process" (GO:0046148), "flavonol metabolic process" (GO:0051554), "flavone biosynthetic process" (GO:0051553), and "flavonol biosynthetic process" (GO:0051555) were significantly enriched in the up-regulated gene sets (shown in oxygen-containing compound, alcohol and lipid (GO:1901700, GO:0097305 and GO:0033993) and "jasmonic acid biosynthetic process" (GO:0009695) were significantly enriched in the down-regulated gene sets. The top significant GO terms found in the up-regulated gene sets at 5 d/0 h were identical to those at 24 h/0 h (Table S3), whereas the terms significantly enriched in the down-regulated gene sets at 5 d/0 h were mainly represented by response to oxygen-containing compound (GO:1901700).

TRANSCRIPTION FACTOR RESPONSES TO COLD STRESS
We surveyed the putative transcription factors (TFs) that were differentially expressed in petunias at one or more time points under cold treatment based on the annotations (Table S2); a total of 61 DGEs were identified as TFs involved in the cold stress response. The number of up-regulated TFs (48) was prominently greater than the number of downregulated TFs (13), suggesting that transcriptional activation but not repression was involved. There are more than 2657  ing 15, 11, 9, 9, 5, and 5 genes of the 61 total cold-responsive TFs, respectively. Members of the AP2-EREBP family belong to different subfamilies of the ERF/AP2 TF family. For example, AtERF025, AtERF003, CaEREBP-C2, and NsERF3 encode members of the A-4, B-6, B-3, and B-1 subfamilies of the ERF/AP2 TF family. Importantly, we observed that more than half of the members of this family were diversely induced at the early stage (0.5 h or 2 h) during cold treatment, suggesting a relatively fast reaction in triggering transcriptional cascades in the cold response of petunias. We arbitrarily selected two AP2/EREBP TFs (cn9604 and GO_dr001P0015C16_F_ab1) for validation by qRT-PCR. As shown in Table 3, both microarray and qRT-PCR analyses were in agreement concerning the substantial inducibility of the transcripts of these two members following cold treatment. Moreover, the highly induced expression of one CBF homologous gene (GO_dr004P0029M08_F_ab1), named PhCBF1 hereafter, was also verified by real-time RT-PCR ( Table 3), implying that the timely response of the DREB/CBF pathway to cold might be similar at least in part between petunias and other plants.
A total of 11 probes fell into three subfamilies of the zinc finger family (C2H2, C3H and B-box). Among these 11 up-regulated zinc finger proteins, we observed that the C2H2 subfamily was particularly enriched in cold-responsive transcripts. Additionally, another C2H2-type zinc finger TF (cn1485) was conspicuously up-regulated following 2 h, 24 h, and 5 d of cold treatment. Because the P-value was higher than 0.05, this TF was excluded from this analysis. However, its cold-induced expression, together with the other four zinc finger proteins (three C2H2 and one Bbox), was verified by real-time RT-PCR (Table 3). We also identified one C3H-type zinc finger gene (GO_drs31P0004P18_R_ab1) whose expression level was evidently increased by cold stress and verified by real-time RT-PCR ( Table 3).
The majority of MADS-box TFs detected by our array were upregulated following 0.5 h or 2 h of cold treatment. Only FBP22 was found to be down-regulated after 2 h under cold stress. Interestingly, all of the MADS-box TFs were first identified to be involved in cold responses. Nine MYB or MYB-related family members showed differential expression, with five exhibiting cold-inducible and four cold-repressed expression. Several transcripts (i.e., cn2063 and cn8223) act as transcriptional regulators in phenylpropanoid metabolism. Five transcripts encoding the NAC domain transcription factor were identified, of which only one was down-regulated on the array. Similarly, four GRAS transcription factors detected on our array were up-regulated by cold, whereas one was down-regulated.
In addition to the above-mentioned TFs, we identified two HLHs (Helix-Loop-Helix) that were cold-repressed in our study. In contrast, two probes were found to encode SBP-box proteins, both of which were up-regulated. Additionally, three other types of transcription factor families, including AUX/IAA, HSF, and TCP, were also identified. Previous studies have linked some HSF TFs to cold stress (Carvallo et al., 2011).

GENES RESPONSIVE TO COLD AND DIFFERENT ABIOTIC STRESSES
To validate the results obtained through microarray analysis, real-time RT-PCR was conducted on 12 selected differentially expressed transcripts, most of which were transcription factors (e.g., AP2-EREBP, C2H2-type, and C3H-type zinc finger protein) and exhibited up-regulated expression ( Table 2). The expression of 75% of the cold-responsive genes (9 of 12) based on the microarray data was confirmed using real-time RT-PCR. The log2 ratios from both the results of the microarray and the qRT-PCR analysis for these genes are shown in Table 3. Their expression kinetics from the real-time RT-PCR results were similar to those of the microarray analysis. These results re-confirmed the veracity of our microarray data. However, the fold changes of three cold-responsive genes measured by microarray and by real-time RT-PCR were not consistent. The relative low correspondence may be partially due to the low cut-off value used for the identification of DEGs.
The expression profiles of seven cold-responsive genes were analyzed by real-time RT-PCR using petunia in vitro shoot cultures treated with various stresses, including cold (2 • C), dehydration, salinity (200 mM NaCl), osmotic (400 mM mannitol), ABA (100 μM), and JA (100 μM MeJA) for 3 h, 6 h and 12 h ( Table 4). From these results, one gene, referred hereafter as PhZFP1 (GO_drpoolB-CL3762Contig1), positively and exclusively responded to cold stress. This gene was predominantly up-regulated by cold and more or less downregulated by other stresses. Several cold-responsive genes were also strongly induced by the dehydration, salt, osmotic, ABA, or JA stresses. For instance, the transcripts of two genes (GO_dr001P0015C16_F_ab1 and cn9788) were differentially accumulated under dehydration, salinity and osmotic stresses and were remarkably induced by ABA. Two genes encoding zinc finger proteins (cn1485 and GO_drs21P0009P06_F_ab1) were strongly induced by all of the abiotic stresses tested (including both ABA and JA) in this study.

NETWORK CLASSIFICATION AND IDENTIFICATION
To identify sets of functionally related genes, we constructed a petunia gene co-expression network using our microarray samples as input. Prior to the analysis, the gene numbers of each module were transformed to the log2 base scale. Then, a clustering approach of the weighted correlation network was undertaken, resulting in a total of 65 modules of highly correlated genes. In Figure 2, the modules are denoted by different colors. Obviously, the numbers of genes (probe sets) per module varied, and more than half of the modules contained less than 500 genes (probe sets) ( Table S4). To explore the co-expression relationships between modules, a module's representative expression pattern was summarized using the first principal component of

www.frontiersin.org
March 2015 | Volume 6 | Article 118 | 7  all of the module's gene members. Using the complete linkage method, all of the module eigengenes were clustered, thereby allowing the characterization of similar structures between the modules (Figure 3). Among the 65 modules, the purple module (size 491) strongly attracted our attention. This module included several transcription factors such as PhZFP1 (GO_drpoolB-CL3762Contig1) and other zinc finger TFs (cn1485 and GO_drs31P0004P18_R_ab1) displaying substantial responsiveness to cold stress in our analysis. PhZFP1 was regarded as a particularly promising candidate gene worthy of further investigation based on its unique stressinduced expression profile compared with other cold-responsive genes tested in our study (Table 4). To quantify the association between PhZFP1 and other petunia genes in the purple module, the top 50 genes positively correlated with PhZFP1 were chosen and are shown in Table S5. Strikingly, PhCBF1 (GO_dr004P0029M08_F_ab1) was also found to be contained in the purple module and included in the top 50 genes positively correlated with PhZFP1. Because PhZFP1 and PhCBF1 were included in the same module and positively correlated with each other, they most likely have a certain relatedness in functionality.

RELATIONSHIP BETWEEN FLAVONOID METABOLISM AND COLD STRESS RESPONSE IN PETUNIAS
The results of the biological process for functional enrichment analysis showed that the categories associated with flavonoid metabolism, including "pigment biosynthetic process" and "flavonol metabolic process," were the most remarkable GO terms. This result suggests that the activation of flavonoid metabolism at the relatively late stage may be an integral part of the plant adaptation to low temperatures in petunias. Previous work has mainly indicated that flavonoids are induced in response to UV irradiation. It was considered that these UV-absorbing compounds could provide protection against cell death stemming from UV-B damage by protecting DNA from dimerization and breakage (Dixon and Paiva, 1995). Other stresses such as cold that induce phenylpropanoids have been less well studied. Increased anthocyanin levels were detected following cold stress in blood oranges (Crifo et al., 2011), but the reasons for this increase are still unclear. In this study, several genes encoding key enzymes implicated in flavonoid metabolism were found among the coldinduced genes at 5 d/0 h (Table S2), such as CHS (chalcone synthase, GI_TC432), CHI (chalcone isomerase, cn302) and F3H (flavanone 3-hydroxylase, cn2956) ( Table S2). These genes are held in common in the biosynthetic pathway of anthocyanins, flavonols, and flavan-3-ols. Our results in petunias are to a great extent consistent with previous observations concerning the increased accumulation of two flavonoid pathway proteins (CHS and F3H) in strawberries in response to cold stress (Koehler et al., 2012). Thus, we speculate that the well-defined induction of these genes in petunias subjected to cold stress would probably lead to a rise in flavonoid levels; in turn, the rising flavonoid levels might contribute to the control of cell osmotic potential, thereby effectively coping with the imposed adverse conditions (Chalker-Scott, 1999). However, it is worthwhile to note that in the present study the noticeable major category predicted to be involved in the cold stress response in petunias was associated with flavonoid metabolism, whereas other significant cold-related GO terms were rare. This may be because a large portion of DEGs share homology with genes encoding putative proteins of unknown function or share no significant homology with any database accession. Future studies will be focused on detailing the quantitative and qualitative transcriptomes of petunias under cold stress. These efforts will be greatly facilitated by the completion of petunia genome sequencing and the rapidly developing technologies for RNA sequencing.

ESTABLISHED AND NOVEL PLAYERS INVOLVED IN COLD STRESS REGULATION IN PETUNIAS
Transcription factors play important roles in plant development and stress tolerance , and genes related to  "Transcription factor activity" may be central regulators involved in upstream cold signal transduction which are capable of triggering a cascade of downstream gene expression. In this study, we found 61 transcription factor genes that were differentially regulated across different time points. The six most highly represented TF families were AP2-EREBP, zinc finger, MADS-box, MYB/MYB-related, NAC, and GRAS. These cold-responsive TFs (shown in Table 2) represent potentially important regulators of the cold-stress response in petunias. The AP2-EREBP family plays a predominant role in the early stages of the cold response, which has been evidenced by some of the best characterized CBFs (CRT/DRE binding factors) in the cold regulatory pathway. CBF transcription factors, attached to the A-1 subfamily of the ERF/AP2 TF family, are major regulators that exert their effects via activation of cold-regulated effectors in Arabidopsis and other plants (Gilmour et al., 1998). Here, one CBF homologous gene (GO_dr004P0029M08_F_ab1) PhCBF1, was highly induced, suggesting that the response of the DREB/CBF pathway to cold in petunias is at least partially similar to that observed in other plants, although it is probably not the only pathway involved. Certain member of the AP2-EREBP family has been reported to show responses to cold stress in Capsicum annuum, such as CaEREBP-C2 (http://www.ncbi.nlm. nih.gov/protein/AAX20035). Other members were first identified as being involved in the cold response, such as NtERF5, which was previously characterized as a pathogen-regulated factor (Fischer and Droge-Laser, 2004). Additionally, NsERF3, which is homologous to two AP2/EREBP family members (cn9604 and GO_dr001P0015C16_F_ab1), was shown to be involved in the regulation of gene expression by stress factors as well as components of stress signal transduction pathways. This gene most likely acts as a transcriptional repressor and may regulate other AtERFs (based on similarity), whereas the regulation mechanism remains obscure (Kitajima et al., 2000), it might be of interest to further study these two AP2/EREBP TFs.
Zinc finger family members in Arabidopsis have been demonstrated to be involved in the cold response . Recently, the genes encoding zinc finger proteins in blueberries were also detected to be cold acclimation-responsive based on transcriptome database mining (Die and Rowland, 2014). In our study, we identified 11 up-regulated zinc finger proteins, among which members of the C2H2 subfamily were especially abundant. Previous studies revealed that many C2H2 zinc finger proteins regulate responses to multiple abiotic stresses, such as ZAT10 and ZAT12 in Arabidopsis (Davletova et al., 2005;Mittler et al., 2006). Thus, C2H2 zinc finger proteins are crucial factors that link various signal transduction pathways in response to different stresses and participate in several cellular processes, whereas only a few C2H2 proteins seem to be associated with responses to one specific stress (Kielbowicz-Matuk, 2012). Consistent with this finding, several petunia cold-responsive C2H2-type zinc finger genes found in our array were induced not only by cold stress but also by other abiotic stresses, including dehydration and salinity. The exception was one C2H2-type zinc finger gene called PhZFP1 (GO_drpoolB-CL3762Contig1). PhZFP1 was exclusively and positively responsive to cold. Interestingly, our findings that PhZFP1 and PhCBF1 were contained in the same module and tightly correlated with each other suggested that this gene may be a critical regulator of the cold stress response that functions in an ABA-independent manner similar to CBFs. Abundant evidence has illustrated that the transcription network of CBF plays a significant role in cold acclimation of evolutionarily different plant species. Nevertheless, gene expression analysis revealed several transcription factors besides CBFs that are induced during cold acclimation. Transgenic analysis of cold-inducible TFs helped to validate their functions in cold tolerance and suggested that other classes of transcription factors also play a significant role in cold acclimation (Chinnusamy et al., 2010). For example, constitutive overexpression of the soybean C2H2-type zinc finger SCOF1 (soybean cold-inducible factor-1) in Arabidopsis transgenic plants elevated the expression of COR (cold-regulated) genes and conferred constitutive freezing tolerance. SCOF1 interacts with soybean G-box binding bZIP transcription factor SGBF1, which is induced by both cold and ABA. It was suggested that SCOF-1 might function as a positive regulator of COR gene expression mediated by the ABA responsive element via protein-protein interactions, which in turn improves cold tolerance of plants (Kim et al., 2001). The large proportion of C2H2-type zinc finger family members with cold-inducible transcripts identified in our study makes this class of TFs an attractive target for further functional characterization. For instance, the question of whether PhZFP1 is involved in the CBF-dependent or CBF-independent transcriptional pathway and the mechanism by which the PhZFP1 transcriptional networks operate during cold acclimation and cold stress tolerance in petunias are unknown. Compared with the relatively large C2H2 subfamily, knowledge about the biological role of the smaller C3H and B-box subfamilies is very limited. C3H-type zinc fingers are a type of plant cold shock domain protein (CSDP). Although their biological function and modes of action are not totally revealed, evidence from Arabidopsis suggests that CSDPs may confer cold tolerance by functioning as RNA chaperones (Kim et al., 2007). All of these results suggest an important role for zinc finger genes in petunia responses to cold stress.
The notable MADS-box family members are known to be key regulators of several plant development processes. The best studied plant MADS-box TFs are those involved in floral organ identity determination (Parenicova et al., 2003). However, a recent study of the Arabidopsis MADS-box gene SUPPRESSOR OF OVEREXPRESSION OF CONSTANS1 (SOC1) revealed that this representative floral activator could also function as a negative regulator of the cold response pathway through the direct repression of CBFs (Seo et al., 2009). This result hints that the responsiveness of the MADS-box transcription factors to cold stress observed in this study might not be a coincidence, and offers the possibility that some MADS-box TFs, for instance, FBP22 (GI_NP1240016) which displays high homology to SOC1 and other SOC1-like proteins are the key regulators of crosstalk between the cold response and flowering time regulation. MYB proteins are crucial factors controlling development, metabolism and responses to biotic and abiotic stresses in regulatory networks (Dubos et al., 2010). A comprehensive analysis of MYB transcription factor gene expression demonstrated that almost all MYB TFs are responsive to hormones or stresses (Chen et al., 2006). In this study, five up-regulated and four down-regulated MYB or MYBrelated family members showed diverse expression patterns. One MYB-related transcription factor (GO_dr004P0017J05_F_ab1) was reported to be a regulator of the tobacco defense-related genes (Sugimoto et al., 2000). Several transcripts, such as cn2063 and cn8223, were identified to be involved in the regulation of anthocyanin biosynthesis (Nakatsuka et al., 2008) (http://www.ncbi. nlm.nih.gov/protein/ABY40371). The differential expression of MYB TFs implies that overlapping environmental pressures or certain metabolic pathways may be integrated into petunia responses to cold stress mediated by the regulation of MYBs.
The appearance of five transcripts encoding the NAC domain TF on the list of transcription factors regulated in response to cold stress is in line with the existing knowledge concerning plantspecific NAC proteins, which are renowned for their roles in biotic and abiotic stress responses (Puranik et al., 2012). In contrast, the currently available knowledge concerning the function of the GRAS gene family is not abundant. Most characterized members of this family encode transcriptional regulators which have functions in plant growth and development, such as axillary meristem formation, root radial patterning and gibberellin signal transduction (Tian et al., 2004). It is worth noting that there have been some results supporting a function for GRAS TFs in plant responses to biotic and abiotic stresses (Kim et al., 2001).
An additional five families (HLH, SBP-box, AUX/IAA, HSF, and TCP) accounted for approximately 11.5% of the total number of cold-responsive TFs detected in our array. Although these families play various roles in plant developmental processes and environmental responses, for instance, the SBP-box gene family has been revealed to control flower and fruit development as well as other significant physiological processes (Chen et al., 2010), most have not been linked to cold stress resistance in plants. Nevertheless, it is important to note the knowledge concerning the function of some members in other aspects has started to gradually increase. For example, several grape SBP-box genes were identified and presumed to be potentially involved in the defense against different stresses, including cold treatment (Hou et al., 2013).
A great number of genes (approximately 37.8% of the total 2071 DEGs) encode proteins of unknown function. The study of these genes may reveal new mechanisms that are fundamental to petunia plants coping with low temperatures. Indeed, the microarray that was used to monitor gene expression was based on ESTs from Petunia hybrida and Petunia axillaries, which was not sufficiently reflective of the overall genome. Thus, the present analysis was significantly limited in scope. Nevertheless, our results support the fact that the large-scale dataset strategy can be useful in obtaining meaningful biological information. In this sense, it is possible that our results have identified for the first time a set of potential regulators which may play a role in the regulation of cold stress tolerance in petunias. Further investigations will be focused on the functional validation of candidate genes for improving petunia tolerance through genetic engineering.

CROSSTALK AMONG DIFFERENT SIGNALING SYSTEMS IN REGULATING THE COLD STRESS RESPONSE IN PETUNIAS
In plants, the existence of complex signaling networks and multiple defense strategies results in their enhanced defense capacity. Because the cold stress-signaling pathway may interact with other signaling systems (i.e., drought, salt, and ABA) (Seki et al., 2002), we conducted qRT-PCR to investigate the expression patterns of seven cold-responsive genes under different abiotic stresses. Most were induced to varying degrees by at least one of the five common abiotic stress treatments, including exposure to dehydration, salinity, osmotic, ABA, and JA in addition to cold, suggesting that they may participate in the regulation of a wide spectrum of responses to diverse abiotic stresses. These results indicate that dynamic crosstalk exists among signaling pathways associated with cold, dehydration, and salt stress in petunias, which may induce common alterations, such as cellular dehydration (Rabbani et al., 2003).
Induced defense responses are controlled by a network of interconnected signal transduction pathways in which hormonal signals, such as ABA or JA, could coordinately activate the transcription of various defense-related genes (Glazebrook, 2005). Plant hormones offer a vital function in signaling networks implicated in plant responses to a variety of biotic and abiotic stresses. ABA plays a critical role in the abiotic stress response by inducing stomatal closure which results in the reduction of transpiration (Pantin et al., 2013) and by regulating root growth, ion channels and gene expression (Duan et al., 2013), whereas JA plays an important role in the defense response, especially against biotic stresses, such as necrotrophic pathogens and wounding (Thaler et al., 2012). In this study, we determined that the expressions of the majority of the cold-responsive genes analyzed (5/7) were significantly modified by ABA treatment (Table 4). Thus, our result is in agreement with the known effect of ABA on cold stress signaling, which indicates their possible involvement in ABA-dependent pathways. However, recent findings have increasingly suggested that JA also plays a role in the regulation of plant abiotic stress tolerance, including stresses, such as salinity and cold (Kang et al., 2005;Khan et al., 2012). These findings may explain our observation of the substantial elevated expression level of two cold-responsive zinc finger genes (cn1485 and GO_drs21P0009P06_F_ab1) following MeJA treatment (Table 4). For example, it was demonstrated that jasmonate acts as an upstream signal of the ICE-CBF/DREB1 transcriptional pathway to positively regulate freezing stress responses in Arabidopsis (Hu et al., 2013). Because the two zinc finger protein genes mentioned above were also responsive to all of the other abiotic treatments, the result suggests that their involvement in the response to cold stress may be dependent on multiple signaling pathways including MeJA molecular signals. Taken together, our analyses indicate that crosstalk among different signaling systems plays a significant role in regulating the cold stress response in petunias.