Introduction

Poikilothermic organisms function optimally in a fairly narrow temperature range, different for different organisms. How the evolutionary history has shaped the temperature requirements of individual species remains poorly understood. This is also true for plants, which often show strikingly divergent tolerance ranges even among closely related taxa. For numerous economically important plant species their cultivation range and productivity are limited by their low-temperature sensitivity. One such species is maize (Zea mays L.), a staple crop worldwide. Most maize cultivars will tolerate spells of moderate cold, but severe cold (0–8 °C) can cause major disturbances of physiology and metabolism or even tissue damage when prolonged (Greaves 1996). At the biochemical and physiological levels, several processes have been identified as key contributors to this sensitivity. They include damage to the photosynthetic apparatus affecting redox balance and enhancing production of reactive oxygen species (Foyer et al. 2002); reduced activity of enzymes of the dark phase of photosynthesis (Leipner and Stamp 2009) and compromised phloem loading (Sowiński et al. 2001), together affecting the carbohydrate status of the leaf (Sowiński et al. 1999; Marocco et al. 2005); alteration of secondary metabolism (Christie et al. 1994); modification of the cell wall (Sobkowiak et al. 2016; Bilska-Kos et al. 2017); and water deficit (Janowiak and Markowski 1994).

A major factor affecting all aspects of plant physiology is the photoperiod, but only recently has it been fully recognized that studies of plant stress responses should not neglect rhythmic gene expression driven by the circadian clock (Bieniawska et al. 2008). Some studies on the Arabidopsis thaliana response to cold have accommodated this fact already (Bieniawska et al. 2008; Espinoza et al. 2010). Notably, expression of C-REPEAT BINDING FACTORS (CBFs) involved in the cold signaling in A. thaliana (Huang et al. 2012) has been shown to be regulated by core clock components CIRCADIAN CLOCK ASSOCIATED 1 (CCA1) and LATE ELONGATED HYPOCOTYL (LHY) through the CCA1-Binding Site (CBS) and EVENING ELEMENT (EE) regulatory elements in their promoters (Fowler 2005; Dong et al. 2011). However, the interplay between cold stress and circadian rhythm has not yet been addressed comprehensively for the evolutionarily distant thermophilic grasses, including maize. One should not expect an exact correspondence between the cold-response mechanisms of A. thaliana, capable of frost hardening upon exposure to severe cold (ca. 4 °C), and that of maize, for which such temperatures can be lethal.

That the transcription of numerous genes is cyclic has been shown for several plant species, including A. thaliana (Harmer 2000), poplar and rice (Filichkin et al. 2011), and maize (Hayes et al. 2010; Khan et al. 2010; Jończyk et al. 2011). Kocsy et al. (1997) showed that chilling (12 °C) changes the rhythmic expression of genes encoding enzymes involved in protection against oxidative stress in maize.

As orthologs of the A. thaliana clock genes have been identified in such diverse plants as Lemna gibba (Miwa 2006) and Castanea sp. (Ibañez et al. 2008) it is tempting to assume that the clock architecture and functioning are conserved in plants. Several clock components have also been found in maize (Hayes et al. 2010; Khan et al. 2010), however, due to its recent whole genome duplication, maize possibly bears two copies of each clock component, as exemplified by the GIGANTEA genes gigz1a and gigz1b (Miller et al. 2008) as well as TOC1 homologs (Jończyk et al. 2011). This makes it possible that the maize clock could exhibit certain idiosyncrasies, although some of the ortholog pairs (Hayes et al. 2010; Khan et al. 2010), but not all (Jończyk et al. 2011), show congruent expression patterns. Indeed, not all plant clocks function exactly in the same manner as in the best-studied model A. thaliana, e.g., in the CAM plants Opuntia ficus-indica and Mesembryanthemum crystallinum (Boxall 2005; Mallona et al. 2011). This could be especially true for monocots showing different diel patterns of leaf growth than dicots due to postulated different responses to the central oscillator (Ruts et al. 2012). In maize, the relations between cold-stress response and the circadian clock could be even more complicated by the fact that it has been bred to be photoperiod-insensitive to allow flowering outside its natural latitude. As flowering is inherently associated with the circadian rhythm the maize clock could be altered by the artificial selection (Colasanti and Muszynski 2009).

In this study we comprehensively investigated the maize response to severe cold in the context of diurnal gene expression by determining transcriptome changes in the leaf. Clustering of transcript profiles and subsequent analysis of enrichment of Gene Ontology categories let us identify structures and processes potentially affected by cold stress. These predictions were extended by in silico analysis of inferred changes in the metabolism. By comparing the present data with previously published transcriptomic studies on maize responses to other stresses we pinpointed the changes unique to the response to severe cold. Transcription of genes orthologous to the A. thaliana clock genes was studied in detail to find out how the low temperature affects the maize circadian clock. Finally, potential regulatory relations between the circadian clock and gene transcription changes in severe cold were identified.

Materials and methods

Growth conditions

Maize (Zea mays spp. indentata) inbred line CM109 was used in experiments. This line was used before by several groups as a chilling-sensitive model line in physiological studies (Sowiński 1995; Kocsy et al. 1996; Verheul et al. 1996; Janda et al. 1998; Haldimann 1999; Bilska-Kos et al. 2017). Kernels were germinated in wet sand in darkness at 25 °C and seedlings were transferred to pots containing Knop’s solution supplemented with Hoagland’s micronutrient solution. Seedlings were grown in a climate chamber (14 h light/10 h darkness, temperature: 24 °C/22 °C, light irradiance: 250 μmol quanta m−2 s−1, relative humidity: 60–80%) to the V3 stage (fully developed third leaf). With the onset of the dark period randomly selected plants were harvested as time-zero sample. Randomly selected half of the remaining plants were transferred to a cold chamber (light/dark temperature: 8 °C/6 °C, other parameters as above) and the other half were grown as above and served as controls. Further samples were taken from chilled and control plants after 200, 400, 600, 810, 1020, 1230, and 1440 min, i.e., the experiment lasted 24 h. The samples comprised the middle part of the third leaf, pooled from three plants and frozen in liquid nitrogen for further use. Four independent biological replications of the experiment were carried out.

Evaluation of physiological response of seedlings to severe cold

Physiological measurements were carried out on plants grown as above. Maximum quantum efficiency of PSII primary photochemistry (Fv/Fm) and PSII effective quantum yield (ΦPSII) were measured with a fluorometer (PAM 200, H. Walz, Germany) in control plants and in plants subjected to the cold for 14 and 38 h (at 10:00 am). For the measurements of Fv/Fm, the saturating one-second light flash intensity was about 3500 µmol quanta m−2 s−1. Before measurements, the plants were dark-adapted for 30 min at 24 °C. ΦPSII was assessed from induction curves with quenching analysis at a 10 ms/pulse sampling rate and actinic light similar to growth conditions. Measurements were performed at room temperature. Experiments were repeated two times using 5 or 6 plants per experimental variant.

Water (Ψw) and osmotic (Ψπ) potential were measured (Ling et al. 1994) with leaf psychrometer L-51A connected to dewpoint microvoltmeter HR-33 T (Wescor, Logan, UT, USA). Measurements were performed at 10:00 am and 6:00 pm on the last day of growth at optimal temperature and on two consecutive days under severe cold. A single 10-mm leaf disc was used for each measurement. Osmotic potential was determined from the same material after freezing the tissue in liquid N2 and thawing. Experiments were repeated two times using five plants (2 or 3 discs per plant) per experimental variant.

To test significance of changes of the physiological parameters we used one-way ANOVA test and Tukey HSD multiple comparison test with significance cut-off set to 0.05. For photosynthetic parameters, Fv/Fm and ΦPSII we performed all pairwise comparisons between control and two cold variants (day 1 and day 2). Similarly for water potential and osmotic potential we compared pairwisely all six time points. Tests were run in R (R Development Core Team 2011).

Microarray description

Two-color microarrays (version 1, 46k) were designed and produced by the Maize Oligonucleotide Array Project (University of Arizona, Tucson, AZ, USA). The microarrays comprised 46,128 probes (mainly 70-mer, but also some 40- and 50-mer) including positive and negative controls (Gardiner et al. 2005).

RNA extraction, labeling, and hybridization

Leaf samples were ground in liquid N2. Further RNA extraction, purification, amplification, Cy3 and Cy5 dye coupling and hybridization to microarrays were done following procedures described by the microarray manufacturer with minor modifications (Sobkowiak et al. 2016). Samples were hybridized in the reference design, i.e., each sample was paired with the time-zero sample from a given biological replication. To cancel potential dye-bias, the labeling scheme was swapped for half of the replications (i.e., two replications: reference-Cy3 vs sample-Cy5 and two opposite ones). For each replication 14 hybridizations were done (56 in all). Microarrays were scanned with a GenePix 4000B scanner, image analysis and fluorescence reading was done using GenePix Pro 3.0 software (Axon Instruments).

Data preprocessing and statistical analysis

Median spot fluorescence minus median local background fluorescence was used in further analyses. Occasionally, when the background-corrected value was negative the log-transformed value could not be determined. In those cases, to fill-in the missing data and obtain a complete expression matrix we used the JBPCAfill program (Oba et al. 2003) which estimates missing values by the PCA method. The completed data was normalized using loess method as described earlier (Jończyk et al. 2011). Further analyses were run on expression fold-change values, i.e., relative to the time-zero sample, in JMP Genomics 6.0.3 (SAS Institute).

To identify genes with transcript level changing over time in response to cold stress, mixed model ANOVA was run with time, temperature, and the time × temperature interaction as fixed effects and replication as random effect. To limit false-positive calls we adjusted the p values with an FDR (False Discovery Rate) multiple comparison correction (Benjamini and Hochberg 1995) set at 0.05. As a result the cold-responsive genes were divided into two groups according to statistically significant effects: of temperature (“temperature”)—showing difference in mean transcript level between cold-treated and control plants, and of the time and temperature interaction (“time × temperature”)—showing major differences in the shape of transcript level curves over time between the cold and control samples.

Raw data from the microarray experiments, described in compliance with MIAME (Minimum Information About Microarray Experiment; Brazma et al. 2001) guidelines, have been deposited at the ArrayExpress database (http://www.ebi.ac.uk/arrayexpress, Parkinson et al. 2010) under experiment accession numbers E-MEXP-3212 (24 °C/22 °C samples only) and E-MTAB-1252 (24 °C/22 °C and 8 °C/6 °C samples). We have previously documented excellent correspondence of the microarray data with RNA quantitation obtained by RT-qPCR (Trzcinska-Danielewicz et al. 2009; Sobkowiak et al. 2014, 2016).

Functional annotation of probes

The microarrays used in this project were designed and characterized basing on the state of knowledge of 2005 (Gardiner et al. 2005). Publication of the B73 maize genome sequence (Maize Genome Sequencing Project; Schnable et al. 2009) and numerous other projects have improved the maize gene annotation. This has led us to carry out a comprehensive re-annotation of the microarray probes.

Probe sequences were aligned to cDNA sequences from ten maize inbred lines (Table S1, Online Resource 1) using BLAST (Altschul et al. 1990). We annotated successfully 80% (34,548 out of 43,405) unique probes (details of the annotation procedure are given in Online Resource 1). Although we matched the probes to genes, further ANOVA and clustering analyses were done at the probe level to avoid being confined to the present state of maize genome annotation and to allow future use of the data.

Clustering and functional analyses

Transcript profiles from the “time × temperature” group were clustered with the nonhierarchical Fuzzy C-Means (FCM) procedure using Mfuzz package (Kumar and Futschik 2007) for R environment (R Development Core Team 2011). The optimal number of clusters and fuzzification parameter were determined by clustering of randomized data according to the standard Mfuzz workflow.

GO enrichment analyses were done with Ontologizer using the Parent–Child Union method (Grossmann et al. 2007) which takes into account the hierarchical structure of the GO system, and the FDR multiple comparison correction (Benjamini and Hochberg 1995) set at 0.05. Only non-redundant sequences were considered, i.e., when several probes matched a single gene, it was counted only once.

Potential changes in metabolism were investigated by in silico analysis using MaizeCyc standalone version 2.0.1 (Monaco et al. 2013) in Pathway Tools visualization environment (Karp et al. 2010). Only genes showing a minimum twofold expression change were used.

Genes without an annotation in MaizeCyc were analyzed based on grouping according to putative function as described by GO, InterPro, MaizeGDB, “Gene Models” (listed on http://www.maizegdb.org; Sen et al. 2009), “Classical Maize Genes” set (Schnable and Freeling 2011), and Grassius Transcription Factor database (Yilmaz et al. 2009).

Search for transcriptional regulation

Potential promoter regions (3 kb upstream) of genes whose transcript showed a significant “time × temperature” effect were searched for putative cis-acting elements in two separate groups: those showing flattening of the transcript level curve over time in response to the cold, and those with a shift of the transcript maximum. TFBSTools (Tan and Lenhard 2016) and JASPAR2016 (Tan 2015) packages in R (R Development Core Team 2011) were used. Redundancy was removed and for each potential cis-acting element Fisher Exact test of over-representation was performed with a p value cutoff of 0.05, with corresponding sequences for all the genes represented on the microarray used as background.

Meta-analysis of transcriptomic data on maize responses to different stresses

To identify transcriptomic changes specific to severe cold we analyzed data deposited in ArrayExpress (Parkinson et al. 2010) and Gene Expression Omnibus (Barrett et al. 2010) databases. We chose projects on maize responses to stress meeting the following criteria:

  1. 1.

    Plants were at similar developmental stages.

  2. 2.

    Seedlings were subjected to treatment of similar duration in different projects investigating a given stress type.

  3. 3.

    Raw data (before normalization) was available.

  4. 4.

    Experiments using inbred lines with known increased stress-tolerance were excluded.

In all, data from ten projects were analyzed, concerning response to severe cold, moderately low temperature, and fungal infection (Table S2, Online Resource 2). Raw data from each project was quality controlled and normalized using standard procedures for each microarray type (details in Online Resource 2). To allow direct comparison the data was transformed to Probability Of Expression scale (POE, Parmigiani et al. 2002) using metaArray package (Choi et al. 2007). Further analysis was done at the gene level. Transcript levels were compared between stressed and control variants for each stress using t-test (multtest package, Pollard et al. 2004) with an FDR correction (Benjamini and Hochberg 1995) of 0.05.

Results

To verify the hypothesis that the maize response to severe cold combines a direct effect of the stress and its effect on the diurnal rhythm, we compared genome-wide transcriptomic profiles during the first 24 h of plant growth in the cold with those in optimal (control) conditions. The genes showing diurnal rhythm of transcription in control conditions have been discussed earlier (Jończyk et al. 2011). Among the cold-responsive genes a statistical analysis distinguished genes whose diurnal profile of transcription was changed in the cold from genes whose transcription showed no temporal changes but was altered by cold. These two gene groups were analyzed in respect to enrichment of Gene Ontology categories.

Physiological reaction to low temperature

The maize inbred line used in this study, CM109, has been used in several studies before and its reaction to low temperatures under both field and controlled conditions is well characterized (Sowiński 1995; Kocsy et al. 1996; Verheul et al. 1996; Janda et al. 1998; Haldimann 1999; Bilska-Kos et al. 2017). Nevertheless, to verify how the plants react to the cold stress under the exact conditions used further to study the transcriptomic response, we performed basic measurements characterizing the photosynthetic efficiency and water potential in the cold. The results (Figs. S1, S2, S3, Online Resource 3) confirmed earlier reports by showing substantial deterioration of all the parameters. Notably, this effect was visible only on the second day of cold exposure, except for osmotic potential which changed already in the first day.

General analysis of transcript profiles

The first group (hereafter referred to as “time × temperature”, reflecting a significant interaction effect) contained 4116 genes whose transcription was changed as an effect of low temperature on their diurnal pattern. Clustering of their circadian transcription profiles under optimal and low temperatures produced eight clusters (in all, 4036 profiles) with either the diurnal transcript profile flattened or the maximum or minimum of transcription delayed by cold. They are shown in Figs. 1, 2, 3 and 4 in the order of the timing of their transcript peak in the control, i.e., night, dawn, day and dusk. Significantly enriched GO terms in the cluster containing genes with transcript maximum shifted from night toward day in the cold (cluster 1, Fig. 1) were related to protein modification, phosphorus metabolism, nucleobase, nucleoside and nucleotide metabolism, and plasma membrane. A substantial number of genes showed just flattening of their transcript profiles in the cold; they lacked any enriched GO categories (cluster 2, Fig. 1). In the cluster of genes peaking at dawn in the control and with a delayed transcript maximum in the cold (cluster 3, Fig. 2) GO categories related to transmembrane transport and protein complex were enriched. The sister cluster with a flattened transcript pattern in the cold contained one enriched GO category, related to hydrolase activity. The two clusters comprised of genes with a peak of transcription during the day in the control and no cycling in the cold (clusters 5 and 6, Fig. 3) contained several enriched GO categories representing various aspects of photosynthesis and heat response. In cluster 7 containing genes with transcript minimum at dawn and a peak at dusk in the control and a flattened profile in the cold (Fig. 4) enriched were GO categories related to biological processes linked to transcription and translation, which was also reflected in Cellular Component and Molecular Function GO classes. In cluster 8 with a transcript peak at dusk in the control (Fig. 4) and progressive down-regulation in the cold no enriched GO categories were found.

Fig. 1
figure 1

Expression profiles of genes responding to cold treatment from the “time × temperature” group with peak of expression at night in control conditions. Diel transcript level profiles were subjected to fuzzy C-means clustering as described in "Materials and methods". Optimal fuzzification parameter was set at 1.20. Transcript level at each time point is shown relative to time-zero sample (log2 [sample/time-zero]) for control conditions (left) and cold treatment (right). Dark period is marked with a black bar. Number of profiles is shown for each cluster. Color scale membership value indicating profile fit to the cluster. Gene Ontology categories enriched in a cluster are listed on the right

Full size image
Fig. 2
figure 2

Expression profiles of genes responding to cold treatment from the “time × temperature” group with peak of expression at dawn in control conditions. See legend to Fig. 1 for details

Full size image
Fig. 3
figure 3

Expression profiles of genes responding to cold treatment from the “time × temperature” group with peak of expression at day in control conditions. See legend to Fig. 1 for details

Full size image
Fig. 4
figure 4

Expression profiles of genes responding to cold treatment from the “time × temperature” group with peak of expression at dusk in control conditions. See legend to Fig. 1 for details

Full size image

To find out if the observed massive and diversified modification of the diurnal gene transcription patterns could be related to disturbances of the central oscillator by the cold, we investigated the transcript profiles of orthologs of the main Arabidopsis clock genes (Fig. 5). Among them we included GRMZM2G145041 recognized as an ortholog of the RVE2 gene, although RVE2 is not considered a component of the Arabidopsis central oscillator. However, the maize gene shows high sequence similarity also to the bona fide clock components RVE4, RVE6 and RVE8, and also its transcript profile in control conditions is very similar to the profiles of RVE4, 6 and 8 (Ko et al. 2016). Notably, for all the genes inspected here their patterns in the cold deviated significantly from those in the control. However, the character of the alteration was diverse: LHY (3 probes; when a gene was represented by more than one probe, mean profile is shown), CCA1 (3 probes), RVE2 (2 probes), and LUX (1 probe) showed a delay of the transcript extremum by several hours and did not return to the time-zero reference level at dusk; TOC1 (3 probes), PRR5 (4 probes), PRR7 (4 probes), and PRR9 (2 probes) exhibited a diminished amplitude of transcript level without a major change in the timing of the maximum; and GI (5 probes) stopped cycling altogether.

Fig. 5
figure 5

Expression profiles of maize orthologs of Arabidopsis main clock components. Transcript level at each time point (min) is shown relative to time-zero sample (log2 [sample/time-zero]). Red profile—control conditions, blue—cold treatment. A–G Profiles of individual genes, as indicated, except A, where mean values for LHY1 (GRMZM2G474769) and CCA1 (GRMZM2G014902) are shown, and F, where three PRR genes are shown. Dark period is marked with a black bar

Full size image

In addition to those responding to both the time of day and the cold stress (“time × temperature”) many genes showed a significant effect of the cold treatment alone (“temperature”, 6523 probes): 3010 were up-regulated and 3513 down-regulated by the low temperature. The enriched GO categories for the up-regulated genes were related to protein refolding, DNA packaging, nucleosome organization, and response to stimulus (Fig. S4, Online Resource 4). Among the down-regulated genes, the photosynthesis and cellular homeostasis GO categories were enriched (Fig. S5, Online Resource 5).

Metabolic pathway analysis

To identify potential metabolic changes deducible from the transcriptomic response to the cold we performed an in silico pathway analysis. Here we analyzed a combined set comprising genes whose diurnal rhythm of transcription was changed by cold (“time × temperature” group) and genes whose transcription was altered by cold alone. The source data are shown in supplementary materials (Table S3, Online Resource 6). The pathways found to be affected were grouped manually into assemblies related to redox balance, metabolism of growth regulators, carbohydrates, and lipids, tRNA charging, and cell wall modifications (Tables 1, 2, 3, 4, 5, 6). Selected genes from these groups, potentially most important for the maize response to low temperature, will be highlighted below and their relevance considered in the Discussion.

Table 1 Metabolic pathways related to redox balance affected by cold inferred from changes in transcription profiles of relevant genes
Full size table
Table 2 Metabolic pathways related to metabolism of growth regulators affected by cold inferred from changes in transcription profiles of relevant genes
Full size table
Table 3 Metabolic pathways related to metabolism of carbohydrates affected by cold inferred from changes in transcription profiles of relevant genes
Full size table
Table 4 Metabolic pathways related to metabolism of lipids affected by cold inferred from changes in transcription profiles of relevant genes
Full size table
Table 5 Metabolic pathways related to tRNA charging affected by cold inferred from changes in transcription profiles of relevant genes
Full size table
Table 6 Metabolic pathways related to metabolism of cell wall components affected by cold inferred from changes in transcription profiles of relevant genes
Full size table

The redox-balance-related genes showed diurnal cycling in control plants with the peak of expression during the light-phase (dawn, day or dusk), lost the clear-cut peak in the cold and instead showed a gradual increase or decrease of transcript level towards the end of the day. We draw attention to a set of four genes encoding violaxanthin de-epoxidase and zeaxanthin epoxidase, enzymes involved in violaxanthin cycle (Table 1).

Several pathways related to plant growth regulators were affected by cold, including genes related to hormone activation and inactivation. Remarkably, ethylene biosynthesis appeared to be induced, since several genes encoding 1-aminocyclopropane-1-carboxylic acid (ACC) synthase and ACC oxidase were up-regulated (Table 2). One should note, however, that the increase of their transcript level was in both cases gradual and slow, suggesting that ethylene is not engaged in the fast response of maize to low temperature. Similar slow up-regulation of genes encoding enzymes related to the biosynthesis of brassinosteroids, jasmonic acid and abscisic acid was found. In contrast, a gene related to ABA inactivation, also up-regulated, was one of the fastest-responding ones. Genes related to IAA biosynthesis changed their diel pattern of transcription markedly. In the control a clear-cut minimum was at dawn, which was replaced by a roughly uniform slow drop throughout 24 h in the cold; as a result, during the night the transcript was more abundant in the cold than in the control. At the same time, however, strong up-regulation of genes related to IAA inactivation was found. Cytokinin metabolism was also affected in a complex way, combining strong and fast down-regulation of cytokinin dehydrogenase involved in cytokinin degradation, and slow and weak up-regulation of zeatin-O-glucosyltransferase responsible for reversible inactivation of cytokinins.

The most remarkable effect of cold treatment on carbohydrate metabolism was a simultaneous induction of genes encoding enzymes engaged in the synthesis and degradation of trehalose and stachyose (Table 3). One of those enzymes, galactinol-raffinose galactosyltransferase, is capable of catalyzing the two opposing processes of stachyose synthesis and degradation. This gene was one of the fastest-responding ones.

Several genes involved in lipid metabolism (Table 4) were affected in a complex manner. Those engaged in biosynthetic pathways showed both induction and repression by the cold, while phospholipases generally showed induction—either gradual throughout the 24 h of cold action (phospholipases C and D) or strong and rapid one followed by a gradual decrease (phospholipase A). Remarkably, in control conditions these three phospholipases showed strikingly different diurnal patterns of expression.

A distinct group comprised 12 genes involved in protein biosynthesis and encoding aminoacyl-tRNA synthetases (ligases) (Table 5). Most of these genes showed a transcript peak at dusk in control plants and their almost universal response to the cold was a loss of cycling.

The largest group of genes affected by cold in this analysis were those related to cell wall modification (Table 6). Most of these genes showed gradual up-regulation in the cold despite showing remarkably diverse diurnal transcript patterns in control conditions.

Transcriptomic response of maize to diverse stresses—a meta-analysis

When responding to unfavorable conditions, plants show a two-pronged reaction: a generalized stress response common to a number of different stresses, and a stress-specific response. In order to identify the transcriptomic response of maize specific to severe cold stress we compared published microarray data concerning severe cold, moderate chilling (ca. 14 °C) and fungal attack. The datasets analyzed shared 16,865 genes; for further analyses we selected genes showing a minimum two-fold statistically significant (p ≤ 0.05) expression change upon any stress. These criteria were met by 971 induced and 385 repressed genes. There was a substantial overlap of the gene sets responding to the cold and chilling stresses, while the vast majority of genes affected by a fungal attack were unique to this stress.

The potentially affected metabolic pathways were analyzed in silico separately for the gene sets specific to each stress and for the set comprising genes common to all three stresses (Table 7). The response specific to severe cold comprised the up-regulation of lipid and fatty acid metabolism pathways, lycopene synthesis, and the synthesis of energy carriers and substrates for other reactions (pentose phosphate pathway, non-oxidative branch; salvage pathways of pyrimidine ribonucleotides). As in the comprehensive pathway analysis presented above, degradation of ABA (phaseic acid biosynthesis) was also induced. Also the simultaneous induction of trehalose synthesis and degradation pathways (see above) was specific to severe cold stress.

Table 7 Metabolic pathways affected specifically by severe cold and those affected by all analyzed stresses (severe cold, chilling, and fungal infection) inferred from changes in transcription profiles of relevant genes
Full size table

Prediction of cis and trans regulatory factors

To identify potential cis-elements engaged in the modulation of the temporal pattern of gene transcription by cold we analyzed putative promoter regions of the cold-responsive genes from the “time × temperature” set for over-represented transcription-factor (TF)-binding sites. The genes were divided into two sets of clusters (see above) showing either amplitude flattening (clusters 2–7) or a shift of transcript maximum (clusters 1, 3, and 8) in the cold. These two sets were analyzed separately as we considered it highly likely that their respective responses to low temperature could be fundamentally different.

Three TF-binding sites were over-represented in the “amplitude flattening” group (Table 8), recognized by WRKY8, CCA1, and NAC043. Genes encoding maize orthologs of two of the cognate TFs, an ortholog of WRKY8 and an ortholog of NAC043, did not show a statistically significant change of transcript profiles in the cold, while an ortholog of LHY1 and CCA1 and thus an element of the core oscillator showed a delay of transcript maximum.

Table 8 cis-elements significantly over-represented in promoter regions of genes showing significant “time × temperature” effect
Full size table

In the “shift of expression maximum” group, five TF-binding sites were overrepresented (Table 8), recognized by FUS3, PIF4, DOF18, OJ1058_F058, and bZIP68. Among the cognate TFs, only bZIP71 (an ortholog of bZIP68) showed a significant change (drop) of transcript level during the day at low temperature. One should note that the relative abundance of all these TF-binding motifs was only moderately, albeit statistically significantly, higher than that expected by chance.

In a partly complementary approach we attempted to identify trans-acting factors possibly involved in the response to the cold by inspecting the meta-analysis data presented above for genes encoding maize TFs affected specifically by severe cold (Table 9). Thirty-four bona fide TFs showed an expression change in severe cold, representing eleven families. The EREB family stood out not only by the number of TFs responding (13), but also by the character of this response: four TFs from this family (ZmEREB92, ZmEREB115, ZmEREB126, and ZmEREB160) showed particularly strong and fast, albeit non-uniform, responses to the cold treatment. ZmEREB115 (C-repeat/DRE binding factor 2) is a maize ortholog of AT4G25470.1 (aka ATCBF2, CBF2, DREB1C or FTQ4) from Arabidopsis.

Table 9 Genes encoding transcription factors altered specifically by severe cold
Full size table

Discussion

While the physiological and biochemical responses of maize to low temperatures are fairly well understood, surprisingly little is known regarding the molecular level. To address this deficiency we followed gene expression profiles in juvenile leaves of seedlings treated with severe cold and used the obtained transcriptomic data to infer which biochemical pathways and physiological processes are likely affected and by what mechanism. This in silico analysis indicated tens of metabolic pathways potentially altered by the cold. It is worth underscoring that we studied the first 24 h of cold action, when the physiological changes due to cold stress are poorly pronounced. Thus, we were looking at molecular changes preceding the physiological reaction and not secondary gene expression changes in response to the physiological alterations. Only the metabolic pathways modified specifically under severe cold stress identified by the meta-analysis of published microarray projects concerning severe cold (<10 °C; Sobkowiak et al. 2014 and this work), moderate chilling (ca. 14 °C; Sobkowiak et al. 2016) and fungal attack (Doehlemann et al. 2008) stresses will be discussed here.

Cold stress and circadian pattern of gene expression

Based on literature data (Kocsy et al. 1997) and our earlier observations (Sobkowiak et al. 2014) we posited that an important part of the maize reaction to severe cold could consist in a disturbance of the biological clock. Here we found that the cold modifies the transcription pattern of a major portion of genes showing diurnal oscillation in control conditions and thus likely controlled by the central oscillator. The change of the expression profiles of the cycling genes was strikingly dual and involved reduction of the peak—essentially flattening the pattern—or replacement of the diel cycle with a fairly monotonous drop or increase of transcript level. Among the genes showing the loss of cycling were those related to the key biological processes of transcription (methylation, DNA packaging) and translation (ncRNA metabolic process, macromolecular complex subunit organization). Since these genes have the transcript maximum at dusk in control conditions, their cold-induced loss of cycling effectively causes relative up-regulation during the day. The presence of “DNA packaging”-related genes in this group suggests that modification of chromatin state is important for the early maize response to severe cold involving both massive up- and down-regulation of genes.

Another important group of genes losing cycling under stress were those related to photosynthesis and showing a transcript peak during the day in control conditions. In their case the flattening of expression profiles amounts to a down-regulation during the day. Also the genes related to photosynthesis that do not undergo diurnal changes of transcript level under control conditions showed repression by the cold. These changes fully support the accepted concept that disturbances of the photosynthetic apparatus are among the most important effects of cold in maize (Leipner and Stamp 2009), here confirmed by direct assessment of activity of photosynthetic apparatus. This effect is usually attributed to photoinhibition resulting from PSII damage (Murata et al. 2007). In our earlier study on the effects of severe cold on two maize inbred lines differing in cold-sensitivity we also observed a drop of photosynthesis-related transcripts, similar in both lines (Sobkowiak et al. 2014). Our results show that the inhibition of photosynthesis under severe cold, in particular during the first hours of treatment is likely not a simple biophysical response but could have a deeper basis in the repression of photosynthesis-related genes. We have demonstrated here that the phenomenon is in part due to a loss of the diurnal profile of expression of these genes.

Anti-oxidant defense is a protective mechanism that in principle could prevent the decrease of the photosynthetic apparatus efficiency under severe cold, but the transcriptomic results do not unequivocally indicate induction of such defense. A notable exception is the “trans-lycopene biosynthesis II (plants)” pathway, involved in the biosynthesis of carotenoids critical in photo-protection and specifically induced in severe cold. Progressive alterations in secondary metabolism on long-term chilling leading to changes in plant pigmentation, including carotenoids, have been described (Haldimann 1999). One may thus hypothesize that the 24 h-long cold treatment studied here could lead to a long-term acclimation response to compensate for a compromised short-term defense of the photosynthetic apparatus.

In addition to the genes which stopped cycling in the cold, a large group reacted by overall up-regulation during the day. Among such genes that could be classed to physiologically meaningful groups (i.e., more specific than those simply ascribed, e.g., a hydratase activity), those related to the plasma membrane and transmembrane transport seemed particularly interesting. Circadian regulation of transporter activity has been found earlier (Haydon et al. 2011) and transcription of some of their genes peaked around dawn (Haydon et al. 2013). After severe-cold treatment, the GO category “organic substance transport” was found to be enriched among cold-regulated genes in two other maize inbred lines (Sobkowiak et al. 2014). Membranes have been postulated as a target for cold stress in maize by a physiological study (Shabala and Shabala 2002). Thus, the present and earlier transcriptomic data seem to confirm that modulation of transport through membranes in leaves of maize seedlings is one of the maize responses to severe cold.

Severe cold stress signaling

As a first attempt at deciphering the potential regulatory mechanisms involved in the maize response to severe cold we followed transcript profiles of genes encoding transcription factors, among them circadian clock components, genes involved in the metabolism of growth regulators, and other potential elements of signal transduction pathways.

The experimental setup did not focus on the immediate response to the cold as the first sampling was done ca. three hours after the temperature lowering. Nevertheless, some genes showed strong up- or down-regulation already at that point, indicating their possible involvement in cold-stress signal transduction. This could certainly be said of the phospholipase A gene showing sharp initial up-regulation. Phospholipase A is a part of the phospholipid stress signaling pathway (Kacperska 2004) and can be involved in the induction of either cell proliferation or apoptosis/necrosis depending on the severity of stress and membrane status (Spiteller 2003). In addition to the fast up-regulation of phospholipase A, phospholipase C and D genes also showed induction, albeit more gradual, further confirming the role of membranes and phospholipids in the response to severe cold in maize.

Potential participation of carbohydrates in cold-stress signaling was suggested by dramatically modified transcript profiles of several genes related to stachyose and trehalose biosynthesis/degradation. Modification of carbohydrate metabolism in leaves is one of the most pronounced maize responses to low temperature (Hodges et al. 1997; Sowiński et al. 1999; Bilska-Kos et al. 2017). It has also been postulated that the changes of the leaf carbohydrate status could affect the water potential balance (Bilska-Kos et al. 2017). Our results confirm that carbohydrate metabolism could be one of the targets of cold in maize leaves, but the meta-analysis has shown that this modification is not cold-specific but rather is part of a general stress response in maize. Only trehalose biosynthesis/degradation and stachyose biosynthesis seem to be altered specifically by severe cold. Both these sugars are known osmoregulators and one could envisage their role in the regulation of water potential under severe cold. Indeed, our measurements showed a pronounced increase of osmotic potential in leaves already during the first day of cold treatment, and a drop of water potential on the second day of treatment. However, one should bear in mind that trehalose is hardly detectable in plants and therefore its role as an osmoprotectant is doubtful. Instead, it could play a role in stress signaling as suggested earlier (Avonce 2004). Stachyose and trehalose have not been mentioned until now in the context of cold stress in maize, but some literature data do show their participation in stress signaling and a metabolic regulatory network in drought stress in other plants (Avonce 2004; Seki et al. 2007). Since severe cold leads to water potential changes (Pérez de Juan et al. 1997), it is not unlikely that trehalose, stachyose, or both take part in the signaling or modulation of maize metabolism under cold stress. A similar role may be played by myo-inositol, an alcohol derivative of carbohydrates. Potential modification of its biosynthesis was found to be one of the specific effects of severe cold in maize as well.

Another large and diversified group of genes involved in the cold stress response comprises those related to growth regulators. Our results indicate profound changes in the metabolism of all major plant growth regulators like ethylene, brassinosteroids and jasmonate and are consistent with their proposed roles in the cold-stress response of maize (Szalai et al. 2000; Janowiak et al. 2002; Du et al. 2013; Vardhini and Anjum 2015; Li et al. 2016). The picture is much less clear regarding the metabolic pathways related to the growth stimulators auxins, cytokinins and gibberellins that appeared to be altered under severe cold in a complex and inconsistent manner. While these data confirm results of other molecular studies on the maize response to severe cold (Li et al. 2000; Sobkowiak et al. 2014), the question of their physiological role in this process awaits a direct study. The most interesting is the case of abscisic acid, since both ABA biosynthesis and degradation appeared to be up-regulated. Notably, up-regulation of the transcript encoding abscisic acid 8-hydrolase involved in ABA inactivation was also found in our earlier report (Sobkowiak et al. 2014) and was stronger in a cold-sensitive inbred line. The CM109 inbred line used in this work is also cold-sensitive, therefore one could hypothesize that inhibition of ABA signaling by cold may occur in maize materials susceptible to low temperature, which is in line with the physiological data that showed a positive correlation between the ABA content in leaves and the level of genotype cold-tolerance (Szalai et al. 2000; Janowiak et al. 2002).

Regulation of transcription under severe cold

Low temperature changes the transcript levels of tens of genes encoding transcription factors (TFs), but here we concentrate only on those affected specifically by severe cold. The EREB family (ERF/AP2 - ethylene responsive factor/apetala2), already known to be induced in maize at low temperature (Leipner and Stamp 2009; Li et al. 2016), was by far the most numerous. Most of these genes were from the APETALA subfamily, and two, ZmEREB56 (dbp4 - DRE-binding protein4) and ZmEREB115 (ereb115 - AP2-EREBP-transcription factor 115), from the DREB subfamily. Notably, the latter, an ortholog of CBF2/DREB1C, a negative regulator of the TFs CBF1/DREB1B and CBF3/DREB1A (Janská et al. 2010), was one of the most strongly and quickly up-regulated genes. These three TFs are components of the CBF/DREB pathway engaged in the A. thaliana response to cold. In maize, ZmDREB1A was found as part of an ABA-independent pathway of cold signaling (Qin et al. 2004), while CBF1 and ZmDREB2A were parts of an ABA-dependent pathway (Qin et al. 2007), indicating that also in maize the cold-stress signal transduction is CBF-dependent (Sharma et al. 2005). On the other hand, some authors (Tondelli et al. 2011) highlight the data suggesting that while CBF genes are indeed involved in the main mechanism of response to the cold in temperate grasses, this is not true for tropical grasses.

Except for the EREB transcription factors discussed above most other TFs showing transcript level changes specific to severe cold have not been ascribed a well-defined role. However, some of them seem likely to be involved in the metabolic changes discussed above. Thus, ZIM15 and ZIM16, upregulated by cold, are orthologs of A. thaliana TFs involved in jasmonate signaling. Three others, GRAS 75, ABI17 and EREB64 are orthologs of TFs involved in the response to gibberellins (in rice), auxins (Arabidopsis), and cytokinins (Arabidopsis), respectively (data from MaizeGDB). As discussed above, the metabolism of these hormones appeared to be modified in severe cold. The parallel changes of transcript levels of relevant TFs support the role of these growth regulators in the maize response to cold stress, postulated by our present and earlier results (Sobkowiak et al. 2014), and literature data (Li et al. 2000).

As one of the initial assumptions of this study, later strongly supported by the clustering data, was that severe cold affects the maize circadian clock, we specifically analyzed the transcript profiles of genes orthologous to the main components of the central oscillator in Arabidopsis. Indeed, we found disruption of the expression patterns of these clock genes in maize by severe cold. Some of them lost their pronounced temporal profile and instead were expressed at a nearly constant level (TOC1, GI, PRR5, PRR7 and PRR9), similarly as has been observed in A. thaliana at 4 °C (Bieniawska et al. 2008); other clock genes, however, responded in a markedly different manner by undergoing a gradual increase (LHY/CCA1, RVE2) or decrease (LUX) of expression, suggesting in all four cases a temporal shift of the transcript extremum by ca. 3 h. This response was not unique to the CM109 line studied here. Also in other maize inbreds studied by us earlier the transcript levels of the clock genes showed similar changes upon treatment with severe (Sobkowiak et al. 2014) or moderate cold (Sobkowiak et al. 2016; Jończyk et al. unpublished data) when analyzed at a single time point (Fig. 6).

Fig. 6
figure 6

Summary of transcript level changes of circadian clock genes in response to severe cold or moderate chilling in various maize inbred lines. Fold-change (expressed as log2) relative to control conditions at 9:30 or 10:00 am is shown. Studied lines were: CM109 (8 °C/6 °C, this report), DH7 and DL3 (8 °C/6 °C, Sobkowiak et al. 2014), CM109 (14 °C/12 °C, Jończyk et al. unpublished data), S6891, S50676, and S160 (14 °C/12 °C; Sobkowiak et al. 2016)

Full size image

Thus, the clock genes behaved in the same manner as did the other maize genes from the “time × temperature” group. This dual mode of transcript level changes of the clock components cannot, however, be easily reconciled with the expression changes of the other genes. In this context we would like to draw attention to the CCA1 gene, as its binding site (CBS, CCA1-binding site) is over-represented in promoter regions of genes from the “time × temperature” group; these are, however, genes showing flattening of the transcript pattern in the cold, while the CCA1 gene itself shows a gradual up-regulation during the 24 h of cold action. This discrepancy indicates that CCA1 does not act as an independent transcription regulator of these genes but rather as a triggering element which, when expressed too late, cannot activate its cognate genes.

Conclusions

Comprehensive transcriptome profiling supported by bioinformatical analysis has indicated remarkable complexity of the reaction to severe cold in maize, suggesting that no single master regulator is involved. The biochemical pathways and regulatory cis- and trans-elements indicated here to participate in the cold response should be investigated directly in future studies. Of particular interest is the link between the massive repression of photosynthesis-related genes during the day and photoinhibition, and the participation of trehalose and stachyose in low-temperature signaling in maize. The obtained results underline that circadian regulation of gene expression is one of the main targets of severe cold and that the transcriptome changes largely precede the physiological response. A large proportion of rhythmically expressed genes, including major components of the central oscillator itself, are affected by the cold in an ordered dual mode involving either flattening of their transcript profiles or a delay of the transcript level maximum/minimum, indicating fundamental differences from the best-studied Arabidopsis model.