Arabidopsis thaliana (ecotype Columbia), wild type (WT), and adh1 mutant seedlings were cultured in a growth chamber (16 h light and 8 h darkness at 23°C) after a 2–4 day vernalization period. For growth under sterile conditions, seeds were surface-sterilized (15 min in 20% bleach, 10 min in 70% ethanol, and rinsed thrice with sterile distilled water) and sown on half-strength Murashige and Skoog (MS) medium supplemented with 1% sucrose and 0.8% agar (pH of 5.7), as previously described (Song et al., 2010). Seedlings were cultivated in Petri dishes for 5 days, were subsequently transferred into potting soil and cultured in a growth house at 23°C (16 h light and 8 h darkness) for 3 weeks. Non-acclimated plants were kept in this optimal environment, while cold acclimation was achieved by subjecting 3-week-old plants to a temperature of 4°C for 7 days.

Chilling and freezing treatments were conducted in low-temp incubators (SANYO, MIR-254, Japan). For chilling stress, 3-week-old WT seedlings were subjected to 4°C for 0, 0.5, 3, 6, 12, 24, 48, and 72 h (for the ADH expression assay). Three plants were randomly selected and used for RNA extraction at indicated time points. For the freezing treatment, 3-week-old WT plants and adh1 mutants with both cold acclimation and non-acclimation were simultaneously treated in a −6°C incubator for 2 h, then all plants were transferred into the growth house at 23°C (16 h light and 8 h darkness) for a 24 h recovery period. Control plants were kept in the growth house and sampled with treatments at the indicated time.

All seedlings were rapidly harvested, flash-frozen in liquid nitrogen, and stored at −80°C until RNA and metabolite extraction.

Total RNA was extracted using the RNA purification kit (Tiangen, Beijing, China) following the manufacturer's instructions. A total of 2 μg of RNA was used for the first-strand cDNA synthesis using reverse transcriptase (Thermo Scientific, Germany) after incubation at 42°C for 1 h 20 min. Semi-quantitative RT-PCR was performed using CbADH1 and CbADH2 specific primers. PCR was performed with 30 cycles for the expression of ADH assays in the mutant seedlings. The reference gene EF4A (TAIR: AT5G60390) was amplified under identical conditions as loading control. The PCR products were separated on 1% agarose gels and stained with ethidium bromide. The utilized primer pairs were ADH1 (F: 5′-ATGTCTACCACCGGACAGATT -3′ and R: 5′- CATGGTGATGATGCAACGAA -3′), and EF4a (F: 5′-TTGGCGGCACCCTTAGCTGGATCA-3′ and R: 5′-ATGCCCCAGGACATCGTGATTTCAT -3′). The real-time qPCR analysis was performed with Platinum SYBR green I qPCR reagents (Takara Bio Inc., Tokyo, Japan). Three technical repeats per sample were used to perform qPCR analysis. UBQ5 (TAIR: AT3G62250) was used as a reference gene to normalize the expression level. The PCR conditions were 35 s at 95°C followed by 40 cycles of the following: 10 s at 95°C, 15 s at 56°C, and 20 s at 72°C. The utilized primer pairs were ADH1 (F: 5′-TATTCGATGCAAAGCTGCTGTG-3′ and R: 5′- CGAACTTCGTGTTTCTGCGGT-3′) for the expression of ADH1 assay in mutant seedlings (Licausi et al., 2010), ADH1 (F: 5′-GTCTTGGTGCTGTTGGTTTAGGC-3′ and R: 5′- GCTTGTCATGGTCTTTCGGGTT-3′) for the expression of ADH1 assay under low-temperature stress, ADH2 (F: 5′-CTTTGAGTGCATCGGGAATGTCT -3′ and R: 5′-GGGTTCGACTCTTGAAACCACCA -3′), and UBQ5 (F: 5′-AAGGTTCAGCGTTTGAGGAAGG -3′ and R: 5′-TCTTTCTGGTAAACGTAGGTGAGTC -3′) (Table S1).

Electrolyte leakage tests were performed as previously described (Rohde et al., 2004). The third to fifth leaves from 18-day old plants were added to capped test tubes (15 ml) and placed in a bath-type programmable freezer (XT5201-D31-R40C, XuTemp, China). The petiole was placed facing downward into the bottom of the tube (without damaging the leaves). Water was added to the bottom of the tube until the petiole was submerged. Samples were kept on ice for 30 min. The plants were exposed to freezing temperatures ranging from 0 to −20°C (0, −4, −6, −8, −10, −12, −14, −16, −18, and −20°C). Tubes were placed on ice after removal from the bath until removal of all tubes was completed. Leaves were immersed in 10 ml distilled water and placed on a shaker for 2 h at 4°C. Electrolyte leakage was determined as the ratio of before and after sample boiling via a conductivity meter (DDSJ-308, Leici, China). The freezing temperature causing a 50% electrolyte leakage (TEL50) was calculated from plotted data of relative electrolyte leakage as the log EC50 value of sigmoidal curve fitting to the leakage values using SPSS PASW Statistics 18.0 (IBM SPSS software, USA).

Stress responses of plants were examined during the short-term cold stress of −6°C for 2 h, subsequently removing the stress and analyzing the recovery response over the following 24 h at 23°C, with cold acclimation (4°C 7 d) and de-acclimation, respectively.

The sampling, extraction, storage, derivatization, standard addition, sample injection and metabolite identification/annotation were performed as previously described (Lisec et al., 2006; Fernie et al., 2011). All shoots samples were collected and froze immediately with liquid nitrogen and kept at −80°C until further analysis. Samples (about 100 mg of fresh weight) stored at −80°C were ground using liquid nitrogen and transferred into pre-chilled 2 ml centrifuge tubes. One thousand two hundred microliters of 100% methanol (pre-cooled at −20°C) was then added and vortexed for 10 s, followed by adding 60 μl of ribitol (Sigma, 0.2 mg/ml stock in dH2O) as an internal quantitative standard and vortexed for 10 s. The tubes were placed into an ultrasound machine at 70°C for 30 min, and then centrifuged for 10 min at 11,000 g and 1 ml supernatant was transferred into 10 ml glass centrifuge tubes. After adding 750 ul chloroform (pre-cooled at −20°C) and 1500 μl deionized water (dH2O) (4°C), the tubes were vortexed for 30 s, then centrifuged for 15 min at 3000 g. One thousand microliters supernatant was transferred into a new eppendorf tube. Samples were blow-dried by moderate nitrogen. Fifty microliters of 15 mg/ml methoxyamine pyridine solution was then added, vortexed for 30 s and reacted for 90 min at room temperature. Finally, 50 μl BSTFA reagent (containing 1% TMCS) was added into the mixture, reacted for 60 min at 70°C. After the above reactions, samples were determined for contents of metabolites using Agilent 7890A GC system coupled to an Agilent 5975C inert XL EI/CI mass spectrometric detector (MSD) system (Agilent Technologies, Santa Clara, CA, USA). Furthermore, a mixed n-alkane standard solutions C8–C20 and C21–C40 (Sigma-Aldrich) was used for the determination of Retention indices (RI).

For GC-MS analysis, Gas chromatography was performed on a HP-5MS capillary column (5% phenyl/95% methylpolysiloxane 30 m × 250 μm i.d., 0.25 μm film thickness, Agilent J & W Scientific, Folsom, CA, USA) to separate the derivatives at a constant flow of 1 mL/min helium. One microliter of sample was injected in split mode in a 20:1 split ratio by the auto-sampler. Injection temperature was 280°C, the interface set to 150°C and the ion source adjusted to 230°C. The programs of temperature-rise was followed by initial temperature of 80°C for 5 min, 20°C /min rate up to 300°C and staying at 300°C for 6 min. Mass spectrometry was determined by full-scan method with range from 35 to 500 (m/z).

For quality control (QC) samples, aliquots (about 40 μl) of all prepared sample extracts were mixed as previously described (Sangster et al., 2006). Furthermore, empty vials referred to as “blank controls” were included into the measurement sequence to test for laboratory contaminations.

For metabolite annotation process, raw GC/MS data was converted into CDF format (Net CDF) files via Agilent GC/MS 5975 Data Analysis software (MSD ChemStation workstation) and subsequently processed peaks identification, peaks filtration and peaks alignment via XCMS (www.bioconductor.org). XCMS's default set with the following changes: xcmsSet (fwhm = 3, snthresh = 3, mzdiff = 0.5, step = 0.1, steps = 2, max = 300), group (bw = 2, minfrac = 0.3, max = 300). The signal integration area of each metabolite was normalized to the internal standard (ribitol) for each sample. Identification of metabolites using the Automatic Mass Spectral Deconvolution and Identification System (AMIDS) was searched against commercial available databases such as National Institute of Standards and Technology (NIST) and Wiley libraries. Metabolites were confirmed via comparison of mass spectra and retention indices to the spectra library using a cut-off value of 70% (Abu Dawud et al., 2012).

For multivariate statistical analysis, the XCMS output was further processed using Microsoft Excel (Microsoft, Redmond, WA, USA) and the normalized data was imported into the Simca-P software version 11.0 (Umetrics AB, Umea, Sweden, www.umetrics.com/simca) (Eriksson et al., 2006) for multivariate statistical analyses, including PCA and PLS–DA. All data were mean-centered and unit variance (UV)-scaled before PCA and PLS–DA applied in order to guard against overfitting. A default 7-fold (Leave-1/7th samples-Out) cross validation procedure and 100 random permutations testing were carried out to guard against over-fitting of supervised PLS-DA models. These discriminating metabolites were obtained using a statistically significant threshold of variable influence on projection (VIP > 1.0). Values were obtained from the PLS–DA model and were further validated via student's t-test (P < 0.05). The metabolites with VIP values above 1.0 and P values below 0.05 (threshold) were selected as discriminating metabolites between two classes of samples.

Hierarchical cluster analysis (HCA) was performed and visualized using the pheat map package available in R 3.0.3 version (www.r-project.org). The data were centered by the average of the controls to create heatmap. Pearson's product-moment correlation (Pearson's ρ) coefficient of 0.9 was performed to calculate the correlation. Corresponding P-values and FDR (false discovery rates) (Benjamini et al., 2001) of each correlation were also calculated using the Cor. test function in R.

To identify coordinated cold treatments induced metabolic perturbation, correlation-based network was used (Hochberg et al., 2013). Pearson's product–moment correlation (Pearson's ρ) was performed on each of the six matrices of data profiles obtained from the two genotypes (WT and adh1) under the three conditions [the control group (CK), cold shock (S), and cold acclimation (A)]. P value was corrected with false discovery rate (FDR). Correlation matrix was calculated using cor. and cor. test function available in R 3.0.3 version (www.r-project.org). The adjacency matrices were construct via the significant correlation (|r| > 0.9 and FDR < 0.05). Network communities and properties were computed by the igraph R package V0.7.1 (http://www.igraph.org/r/). Average node degree or average links: the average number of edges per node in the network; network density: characterizing the proportion of edges in a network in relation to the total number of all possible edges in the network (Barabási and Oltvai, 2004; Hochberg et al., 2013). The above properties were investigated. For comparative analysis, the network intersections and network symmetric difference were determined among the optimum temperature, cold shock and cold acclimation, between adh1 and WT. Between-genotype differences are visualized in the symmetric difference networks (SDN) by R. Cytoscape (Shannon et al., 2003) version 3.2.0 (http://www.cytoscape.org/).

Identified metabolites were mapped onto general biochemical pathways according to the annotation in KEGG (Kyoto Encyclopedia of Genes and Genomes) (Kanehisa and Goto, 2000) (http://www.genome.jp/kegg/). Metabolic pathway maps were constructed by incorporating the identified and annotated metabolites using R 3.0.3 version (www.r-project.org) and Power Point (Microsoft, Redmond, WA, USA).

In order to verify whether the ADH gene was induced by low temperatures, the expression levels of ADH1 and ADH2 at 4°C were tested via qPCR analysis. Conditions were optimized to show induction relative to basal levels (time 0 h). Following treatment with 4°C for 72 h, the ADH1 transcript level significantly increased, whereas ADH2 showed no significant increase (Figure 1). The 4°C stress caused ADH1 transcript induction at 3 h following stress treatment, and the expression level declined after initial rise due to 24 h treatment.

To further reveal the in vivo role of ADH1 for plant cold tolerance, we isolated T-DNA mutants with T-DNA insertion in the exon of ADH1 [adh1-1 (SALK_052699) and adh1-2 (SALK_825718)] (Figures 2A–C). Using semi-quantitative RT-PCR analysis, the transcript of ADH1 was tested in adh1-1 and adh1-2 mutants (Figure 2B). Quantitative real-time PCR analysis showed that the transcript level of ADH1 was completely inhibited in the adh1-2 mutant (Figure 2C).

To determine the cold hardiness of the adh1 mutant, the electrolyte leakage rates in WT, adh1-1, and adh1-2 were tested after freezing treatments in the panel, and their semi-lethal temperatures (LT₅₀) were calculated based on changes in electrolyte leakage rates under different low-temperature stresses using the logistic equation.

The results revealed that with decreased temperature and prolonged treatment time, the electrolyte leakage rates showed an S-shaped curve (Figure 3). The semi-lethal temperatures (LT₅₀) of WT, adh1-1 (SALK_052699), and adh1-2 (SALK_825718) were not significantly different (Figures 3A,C). However, electrolyte leakage rates of WT, adh1-1, and adh1-2 subsequent to cold-acclimation (4°C for 7 d) were −6.7°, −6.42°, and −5.95°C, respectively (Figures 3B,D). The adh1-2 mutant showed a higher percentage of ion leakage than WT after freezing treatment, suggesting that plasma membranes were more damaged in the mutant. Without cold adaptation, adh1-2 also showed lower basal freezing tolerance than WT plants.

To better evaluate the role of ADH1 in plant freezing tolerance, the metabolites profiling adh1-2 responses were analyzed via GC-MS. Treatments were −6°C for 2 h with non-acclimation and cold acclimation (cold shock and cold acclimation) followed by growth at 23°C for 24 h (recovery). In total, 263 resolved peaks were identified, among which 78 could be identified as discrete metabolites with known chemical structure (Supplementary File 1). A considerable amount of material has not been annotated as a result of low signal or missing database. The results revealed that amino acids, organic acids, and several other substances changed in the treatment groups (Figure 4).

Principal component analysis (PCA) of all samples was performed (Figure 5), revealing that the metabolite differences were not significant between WT and adh1-2. However, they were significantly different between cold shock, cold acclimation, and room temperature controls. The time series effects (freezing treatment for 2 h and recovery for 24 h in optimal temperature) clearly contributed most to total variance within the dataset (Figure 5).

In the GC/MS-based analysis of central metabolites, the first principal component (PC1) and the second principal component (PC2) explain the highest amount of variance (50.2%) across the dataset, separating the samples (Figure 5B). PC1 explained 36.6%, and Phenylalanine, Isoleucine, and 4-hydroxybenzoic acid were the main metabolites contributing to the dispersion of samples on PC1 (Figure 5B). The second component (PC2) explained 13.6% of the variance and separated the samples (Figure 4A). Changes in the abundance of glucose, maltose, and gentiobiose were predominantly responsible for sample dispersion along PC2 (Figure 5B). Detailed information is presented in the supplementary file (Supplementary File 2). The PCA results are consistent with PLS–DA analysis (Figure S1).

To better evaluate the changes in the level and metabolism of each metabolite, students' T-test was used to reveal that the levels of 78 metabolites significantly varied (p < 0.05). These included several soluble sugars, amino acids, organic acids, nucleotides, phosphates, and glycerol, as well as sugars such as sucrose, glucopyranose, melibiose, gentiobiose, maltose, fructose, cellobiose, xylose, inositol, maltotriose, and glucose (Table 1), suggesting that ADH1 affected the accumulation of soluble sugars as cold stress response. For example, sucrose (P < 0.05, fold-change = −14.723 in the control group at 23°C) was significantly reduced in the background level in WT compared to adh1 seedlings. However, when cold acclimation preceded freezing stress, this gap of sucrose content had narrowed, and sucrose was significantly increased in WT compared to adh1 mutants of seedlings (P = 0.382, fold-change = 2.751 in the freezing treatment −6°C for 2 h in the cold acclimation group; P = 0.446, fold-change = 1.959 in recovery at 23°C for 24 h). Of particular interest were the cold shock results, where the reduced state was maintained continuously (P < 0.05, fold-change = −1.558 in the cold shock group −6°C for 2 h; P = 0.241, fold-change = −4.141 in recovery at 23°C for 24 h) (Table 1 and Figure 6A). A similar situation was detected for glucopyranose and xylitol. However, gentiobiose increased in the background level in WT compared to adh1 seedlings but significantly decreased in the freezing treatment with cold acclimation of 24 h, and changed only a marginally in the cold shock treatment (Table 1 and Figure 6B). Glycerol showed an increase in cold treatments but remained basically the same between WT and adh1 in the control group at 23°C (Table 1 and Figure 6C). Detailed information is presented in Table 1 as well as in Supplementary File 1.

Correlation-Based Network Analysis (CNA) to Identify the Freezing Treatment Induced Metabolic Perturbation

To further screen the calculated correlations, seven networks were generated (Supplementary File 3, Figure 7, and Figure S2). At a false discovery rate (FDR) < 0.05 and r > 0.8, the network of the adh1 mutant and WT changed little in both the control group and the cold acclimation group (Figures 7A,B). However, cold shock caused a significant increase in the number of edges in the adh1 mutant from 334 to 1410, and in WT from 378 to 2700. Network density changed in the adh1 mutant from 0.105 to 0.309, and in WT from 0.088 to 0.669. The average node degree changed from 11.8 to 42.5, and from 11.5 to 84.4, respectively. In summary, under cold shock, WT networks were characterized by increased numbers of edges compared to adh1 mutant networks (2700: 1410), which markedly increased by 50% (Figure 7C and Supplementary File 3), suggesting a significant effect of adh1 mutation on the metabolism under cold shock.

The difference of adh1 mutant and wild type in metabolic response to cold treatment can be appreciated by comparing the graphs. To further visualize the between-sample effects on metabolic networks, symmetric difference networks (SDN) were utilized (Figure 7), in which blue edges were specific to adh1 and red edges were specific to WT. The SDN is mainly used for network connections under different conditions and emphasizes the presence of different connected samples specific nodes. E.g., under cold shock, SDN comprised 132 nodes and 1870 edges, 68 nodes were specific to adh1 and 64 nodes were specific to WT, 290 edges were specific to adh1, and 1580 were specific to WT, resulting in a ratio of 1 to 5.45 (Supplementary File 2). Notably, all 40 edges of maleic acid (MA) and all 44 edged of asparagine (Asn) were specific to adh1, while 26 edges of maltotriose (Mat) and 42 edges of gentiobiose (Geb) were specific to WT (Supplementary File 3: (S)WT-ADH).

This result revealed the differences of extensive topological networks between adh1 and WT due to different cold treatments. E.g., in both strains (adh1 and WT) adenosine played a central role in the structure of the network under cold treatment (cold shock and cold acclimation), but under control conditions (constant mild temperature), most of the connectivity was lost (Supplementary File 3, Figure 7, and Figure S2). Metabolites in adh1 altered between cold shock (−6°C for 2 h) and cold acclimation (4°C for 7 d, then treated at −6°C for 2 h). According to the connectedness of correlation-based network analysis, melibiose, fumaric acid, succinic acid, glycolic acid, and xylose lost the relations with the network of cold acclimation, but enhanced relations with the network of cold shock. These results have not been found for the WT network.

The metabolite abundance indicates changes in metabolic activity of a pathway. To assess how plant metabolism deals with cold treatments and recovers in mild temperatures (23°C for 24 h), a comparative metabolic analysis was performed for both adh1 and WT (Figure 8).

Firstly, the change fold metabolites of WT/adh1 were insignificant in this analysis except for sucrose in the control group, which showed decreasing levels (the background level in 23°C), suggesting high enrichment of sucrose in the adh1 mutant. Interestingly, under freezing treatment (with or without cold acclimation), the change-fold of metabolites was significantly increased; however, after recovery for 24 h, all metabolites returned to a level consistent with the control group. The response of increasing levels of succinic acid, glutaric acid, and maleic acid were observed in cold acclimation in WT/adh1 analysis, suggesting their enrichment was affected in adh1 mutation, as were the decreases in the levels of tyrosine, phenylalanine, valine, asparagine, lysine, fumaric acid, and pyroglutamate during the recovery phase in the acclimation group (cold acclimation 4°C for 7 d, followed by treatment at −6°C for 2 h). The most severe recovery responses were observed for maleic acid, and following the removal of cold treatment, this strongly reduced compared to the control level, causing a corresponding change of fumaric acid.

Secondly, the changes of metabolites were detected between the freezing treatment with cold acclimation (A: cold acclimation at 4°C for 7 d, then treated at −6°C for 2 h) and cold shock without cold acclimation (S: cold shock −6°C for 2 h) in adh1 and WT. The fold-change of metabolites of S/A in this analysis is presented in Figure 8. Almost all detected metabolites of S/A decreased in adh1, i.e., these metabolites accumulated more during cold acclimation than during cold shock. However, the metabolites increased only after cold treatment was removed and recovery had lasted for 24 h in adh1 mutant seedlings. Similar results were found for WT. The only difference is that asparagine did not increase, but decreased during recovery. However, for sucrose, the fold-changes of metabolites of S/A were essentially identical for both cold acclimation group and cold shock in adh1, but were reduced in WT.

Furthermore, the fold-change of metabolites of A/CK in adh1 and WT was analyzed. galactinol, melibiose, cellobiose, gentiobiose, maltose, gluconic acid, proline, and glycine levels increased both in adh1 and WT under cold acclimation, and subsequently returned to control level or maintained the response level after removal from the cold treatment and 24 h recovery. In contrast, the fold-change of metabolites of S/CK in adh1 and WT, beta-alanine, and glycine significantly changed from initial reduction to a dramatic increase in cold shock and recovery. However, they have been increasing in the cold acclimation group.

In summary, ADH1 abundantly affected some metabolites, especially numerous sugars (e.g., sucrose), and amino acids (e.g., asparagine). These metabolites are intermediates of glycolysis, tricarboxylic acid cycle, and anaerobic respiration.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.