Next Article in Journal
Atypical NF1 Microdeletions: Challenges and Opportunities for Genotype/Phenotype Correlations in Patients with Large NF1 Deletions
Next Article in Special Issue
In Silico and Transcription Analysis of Trehalose-6-phosphate Phosphatase Gene Family of Wheat: Trehalose Synthesis Genes Contribute to Salinity, Drought Stress and Leaf Senescence
Previous Article in Journal
Generalizability of GWA-Identified Genetic Risk Variants for Metabolic Traits to Populations from the Arabian Peninsula
Previous Article in Special Issue
HD-ZIP Gene Family: Potential Roles in Improving Plant Growth and Regulating Stress-Responsive Mechanisms in Plants
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Transcriptome Profiling of Maize (Zea mays L.) Leaves Reveals Key Cold-Responsive Genes, Transcription Factors, and Metabolic Pathways Regulating Cold Stress Tolerance at the Seedling Stage

1
Biotechnology Research Institute, Chinese Academy of Agricultural Sciences, Beijing 100081, China
2
Maize Research Institute, Heilongjiang Academy of Agricultural Sciences, Harbin 150086, China
3
National Key Facility for Crop Resources and Genetic Improvement, Institute of Crop Science, Chinese Academy of Agricultural Sciences, Beijing 100081, China
4
National Agricultural Science and Technology Center, Chengdu 610213, China
*
Author to whom correspondence should be addressed.
Genes 2021, 12(10), 1638; https://doi.org/10.3390/genes12101638
Submission received: 13 August 2021 / Revised: 27 September 2021 / Accepted: 11 October 2021 / Published: 18 October 2021
(This article belongs to the Special Issue Genetics and Evolution of Abiotic Stress Tolerance in Plants)

Abstract

:
Cold tolerance is a complex trait that requires a critical perspective to understand its underpinning mechanism. To unravel the molecular framework underlying maize (Zea mays L.) cold stress tolerance, we conducted a comparative transcriptome profiling of 24 cold-tolerant and 22 cold-sensitive inbred lines affected by cold stress at the seedling stage. Using the RNA-seq method, we identified 2237 differentially expressed genes (DEGs), namely 1656 and 581 annotated and unannotated DEGs, respectively. Further analysis of the 1656 annotated DEGs mined out two critical sets of cold-responsive DEGs, namely 779 and 877 DEGs, which were significantly enhanced in the tolerant and sensitive lines, respectively. Functional analysis of the 1656 DEGs highlighted the enrichment of signaling, carotenoid, lipid metabolism, transcription factors (TFs), peroxisome, and amino acid metabolism. A total of 147 TFs belonging to 32 families, including MYB, ERF, NAC, WRKY, bHLH, MIKC MADS, and C2H2, were strongly altered by cold stress. Moreover, the tolerant lines’ 779 enhanced DEGs were predominantly associated with carotenoid, ABC transporter, glutathione, lipid metabolism, and amino acid metabolism. In comparison, the cold-sensitive lines’ 877 enhanced DEGs were significantly enriched for MAPK signaling, peroxisome, ribosome, and carbon metabolism pathways. The biggest proportion of the unannotated DEGs was implicated in the roles of long non-coding RNAs (lncRNAs). Taken together, this study provides valuable insights that offer a deeper understanding of the molecular mechanisms underlying maize response to cold stress at the seedling stage, thus opening up possibilities for a breeding program of maize tolerance to cold stress.

1. Introduction

Maize (Zea mays L.) is the world’s most commonly grown cereal crop, with an estimated global annual production of about 1186.86 million metric tons in 2020/2021 [1]. The high dependence on maize for human, animal, and industrial consumption makes it one of the most critical food crops. However, maize growth and yield are highly dependent on sufficient environmental factors [2]. Thus, the current and expected scarcity of water sources and arable land due to the increasing world population and the recurrent extreme weather caused by global warming is projected to increase the incidence of abiotic stresses, such as drought, cold, and freezing during the planting, flowering, and grain-filling stages, in many corn-growing areas [3]. These abiotic stresses typically serve as crucial impediments to maize production and geographical distribution [4] and restrict agricultural yields worldwide.
Since maize has a tropical origin, cold stress is a significant risk factor among the several abiotic stresses in the development of maize. A previous report has shown that cold stress adversely affects maize growth from germination to harvest, resulting in significant yield losses due to low and slow germination and poor grain filling [5]. Corn production losses can surpass 20% in the most prolonged cold temperatures [6]. Therefore, the development of high-yielding cultivars tolerant to cold stress may help in augmenting maize production in vulnerable regions and act as an essential maize-breeding target.
The optimum maize growth temperatures range from 21 to 27 °C, while sub-optimal temperatures of about 10–20 °C decrease biomass production, thereby leading to growth retardation [7]. Cold stress induces multiple abnormalities in physiological, molecular, and biochemical processes, which harm plant growth and yield. Cell membranes may become disorganized, proteins may be denatured, oxidative defense and osmotic stress may be altered, photosynthesis possibly restricted, and metabolism may become dysfunctional, all of which subsequently disrupt growth and development, decrease fertility, and cause premature senescence and even plant death [8,9,10,11]. All of these cold stress-associated changes occur through accurate gene expression regulation and are therefore genetically regulated. Thus, screening cold stress-related candidate genes may help identify essential regulators and pathways as potential targets for breeding resistant varieties adaptable to environments with fluctuating temperatures.
As a result of their sessile nature, plants have developed complex cold acclimation mechanisms, which entail the interaction of multiple biochemical pathways in an organ-, genetic-, and environmental-specific manner [12]. They sense cold stress through changes in membrane fluidity and the accumulation of calcium signatures, leading to downstream activation of cold signaling pathways [13]. The enhanced cytosolic Ca2+ levels then induce C-repeat binding factors (CBFs), which act as core regulators for expressing cold-response genes [14,15]. The stress receptors, in conjunction with the cell membrane transporters, facilitate the perception of stress signals and their transmission to target genes. Multiple protein kinases, including CaMs, CMLs, CBLs, CDPKs, and MAPKs, phosphorylate other kinases and/or various TFs, resulting in activation of the cold-responsive genes [16]. Moreover, the transcriptional factor (TF) families bHLH, CAMTA, MADS, WRKY, NAC, TRAF, C3H, and AP2/ERF are critical in cold response mechanisms, while phytohormones, such as abscisic acid (ABA), regulate specific pathways that lead to cold tolerance [17,18]. Cellular redox homeostasis is protected by synthesizing defense enzymes and other antioxidant systems, while soluble sugars serve as a stabilizer of cellular components and plasma membrane [19]. Secondary metabolites, such as lignin, anthocyanin, terpenoids, chaperones, and late embryogenesis abundant (LEA), provide cold tolerance by protecting cellular components from cold-induced cellular damage [20,21]. Transporters, such as the ATP-binding cassette (ABC) transporter, play integral roles in plant growth and development, homeostasis of phytohormones, and resistance to abiotic stress [22]. Tremendous progress has been made in elucidating the mechanisms underlying cold tolerance in plants. However, the complex molecular mechanisms of cold tolerance in maize seedlings are still elusive and require comprehensive research.
Moreover, considering the genetic diversity of maize inbred lines, it will be interesting to identify the stress-responsive genes that have a consistent function in a variety of inbred lines in spite of their genetic background. Differential transcriptome analysis using RNA-sequencing (RNA-seq) approaches has recently emerged as robust, reliable, and responsive to broader levels of gene expression [23]. This effective technology makes it easier to rapidly classify stress-responsive genes and decode metabolic pathways associated with biotic and abiotic stresses [24]. Currently, little is known about the transcriptomic responses of maize seedlings to cold stress. In this study, we used RNA-seq analysis to decipher the expression profiles of the differentially expressed genes (DEGs) responsible for the contrasting cold response of 46 maize inbred lines (24 tolerant and 22 sensitive) at the seedling stage. The common cold-responsive genes were then characterized by their patterns of expression and evaluated for their functional significance. The current study provides valuable clues for the in-depth characterization of molecular responses of maize seedlings to cold stress, which could lead to effective strategies for breeding and developing cold-tolerant maize varieties.

2. Materials and Methods

2.1. Plant Materials and Treatments

The maize inbred lines were derived from the hybrid 19NL, which was provided by Heilongjiang Academy of Agricultural Sciences. The maize hybrid 19NL is a highly suitable cultivar for the spring season in northeastern China. Field experiments were conducted at open field stations during the spring maize growing season (March–June 2019) in Heilongjiang (Harbin, China). In April 2019, approximately 6000 inbred lines of 19NL were planted in rows, with 40 cm between the rows and each row containing 20 seedlings. However, in the late spring of May 2019, due to climate change, the Heilongjiang region of Harbin was affected by a cold spell of below 10 °C for more than three days. We observed that the cold spell impacted the seedlings differently, with survival rates varying from one inbred line to another. This indicated that their response to varying degrees of cold stress was different, regardless of being from the same germplasm. When seedlings had six fully expanded leaves (40 days old), all inbred lines were sampled, and 46 inbred lines with contrasting cold tolerance were selected and classified into cold-tolerant (24 lines) and cold-sensitive (22 lines) lines based on the survival rates of the seedlings, as well as a visual observation of the phenotypic changes of the leaves. The top fully expanding leaves of the 24 cold-resistant and 22 cold-sensitive inbred lines were harvested, frozen into liquid nitrogen, and later stored at −80 °C for further use.
To validate the expression of our cold-responsive DEGs, we planted B73 and CIMBL116 maize inbred lines, which have previously been reported as being cold-sensitive and -tolerant lines, respectively [25,26]. The seeds from these two maize inbred lines were provided by the Chinese Academy of Agricultural Sciences’ crop science institute). Ten seeds from each inbred line (B73 and CIMBL116) were surface sterilized with 75% (v/v) ethanol for three minutes before being rinsed three times with distilled water. Seeds were then placed between two layers of damp paper at 25 °C and left to germinate in the dark for 3 days. Uniformly germinated seeds with 2–3 cm coleoptiles were selected and sown in pots filled with peat, vermiculite, and perlite (10:1:1 by vol.). The seedlings were then grown in a growth chamber with a controlled temperature of 25/20 °C (day/night), 450 L mol m−2 s−1 light density, and a 12/12 h (light/dark) photoperiod until the third leaves were fully developed. The seedlings from the two inbred lines were divided into two groups, with the first receiving a cold stress treatment of 4 °C for 2 h followed by 25 °C for 2 days. The other group, which served as a control, was kept in the growth chamber under the same conditions as described above. The control and cold treatment samples were harvested at the same time after 2 days and immediately frozen in liquid nitrogen before being stored at −80 °C for total RNA isolation.

2.2. RNA Extraction, Library Construction, and Illumina Sequencing

TRIzol reagent (Invitrogen, CA, USA) was used to isolate total RNA from 24 tolerant and 22 sensitive leaf samples as per the manufacturer’s standard. The samples were treated with RNase-free DNaseI (Takara, Kusatsu, Japan) to remove the genomic DNA. The NanoDrop 1000 (NanoDrop Technologies, Wilmington, DE, USA) and Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) were used to assess RNA concentration and integrity, respectively. The cDNA libraries were constructed and sequenced using the Illumina HiSeq™ 2500 platform to generate 150 bp paired-end reads. Moreover, the above procedure was also employed to extract RNA from B73 and CIMBL116 at control and treatment levels for the purpose of qRT-PCR.

2.3. Reads Processing, Mapping, and Gene Expression Quantification

We used FastQC (V0.11.3) to evaluate the quality of raw reads, while Trimmomatic (V0.32) was utilized to eliminate low-quality and adapter-containing reads [27]. The Phred quality scores, including Q20 (99% base call accuracy), Q30 (99.9% base call accuracy), and the GC content of the clean data, were calculated. Consequently, high-quality clean data were used in all the subsequent analyses. The maize reference genome (B73_v4) was downloaded from the maize database (http://www.maizegdb.org/genome/genome_assembly/Zm-B73-REFERENCE-GRAMENE-4.0, accessed on 15 November 2019). All the clean reads obtained from the 46 samples were aligned to the maize B73_v4 reference genome using HISAT2 (V2.0.5) [28] with default parameters. The aligned reads were assembled into transcripts, and the transcripts from all samples were merged using Cufflinks [29]. The assembled transcripts were compared to the reference annotation by Cuffcompare. The HTSeq tool [30] was used to count the number of fragments mapped to each gene, and the transcripts per million (TPM) for each unigene as the expression level was calculated. Principal component analysis (PCA) of gene expression levels was performed to calculate the distance between samples using the clustering method. Differential expression analysis was performed using the DESeq2 R package to identify DEGs in general for the whole transcriptomes of the tolerant and sensitive samples. To estimate expression level, the DESeq2 program was used to normalize the number of counts of each sample gene using the base means. The difference was calculated, and the statistical significance was determined using the negative binomial distribution test [29,31]. Genes with p-value ≤ 0.05 and an absolute value of log2 fold change ≥1 or ≤−1 between tolerant and sensitive samples were considered differentially expressed.

2.4. Functional Annotation of the DEGs

For functional annotation, the 2237 transcripts that qualified to be our DEGs were annotated against the maize genome (AGPv4, B73 RefGen_v4) (http://ensembl.gramene.org/Zea_mays/Info/Index, accessed on 17 November 2019). In total, 1656 (74%) DEGs were annotated, and 581 DEGs were unannotated. To elucidate the function of the 581 unannotated genes, we applied several previously published procedures to identify high confidence lncRNAs [32]. Briefly, (i) unannotated DEG lengths were confirmed to be longer than 200 nucleotides for further analysis; (ii) DEGs that encode open reading frames (ORFs) of 120 or fewer amino acids were retained as lncRNA candidates; (iii) DEGs with similarity to known proteins based on BlastX against the SWISS-PROT database were filtered out; (iv) all the 581 unannotated DEGs were further evaluated using Coding Potential Calculator (CPC) (http://cpc.gao-lab.org/, accessed on 20 March 2020) [33], which assesses the coding probability of transcripts; (v) a total of 337 high confidence drought-responsive lncRNAs were obtained by comparing the output of the two procedures.

2.5. Gene Ontology (GO) Enrichment and KEGG Pathway Enrichment Analyses

The GO enrichment analysis of DEGs was conducted by agriGOv2 (http://systemsbiology.cau.edu.cn/agriGOv2/, accessed on 12 December 2019) [34]. Significant enriched GO terms were determined by the p-value ≤ 0.05 with the Fisher’s exact test and the Bonferroni multi-test adjustment. Redundant GO terms were removed using the web tool Revigo [35]. Significantly enriched GO terms were assigned to the GO categories of biological process (BP), molecular function (MF), and cellular component (CC). The KEGG (http://www.Genome.jp/kegg/, accessed on 15 December 2019) database [36] was used to analyze the functional involvement of DEGs in various metabolic pathways. Furthermore, the statistical enrichment of DEGs in KEGG pathways was tested using the KOBAS 3.0 webserver (http://kobas.cbi.pku.edu.cn/waitkobas.php, accessed on 15 December 2019) [37], while the criteria for substantially enriched KEGG pathways was a p-value ≤ 0.05. A co-expression network was constructed using the R package based on a weighted gene co-expression network analysis (WGCNA) to identify significant hub genes associated with cold tolerance in maize.

2.6. Validation of Cold-Responsive DEGs by Quantitative Real-Time PCR (qRT-PCR)

To validate the reliability and repeatability of the RNA-Seq data, six DEGs were randomly selected for verification by qRT-PCR. The gene-specific primers (Table S1) were designed using Primer Premier 5.0 software (Premier Biosoft International, Palo Alto, CA, USA). The qRT-PCR was conducted in triplicate using 2× ChamQ SYBR qPCR Master Mix Kit (Low ROX Premixed) (Vazyme Biotechnology Co., Nanjing, China) on an Applied Biosystems QuantStudio® 6 Flex (Thermo Lifetech, Waltham, MA, USA). The reaction system utilized was as previously described by Zhang et al. [38]. The internal reference β-actin was utilized to normalize the expression data. The relative expression levels of the six DEGs were calculated according to the 2−ΔΔCT (cycle threshold) method [39].

3. Results

3.1. Phenotypic Analysis of Maize Population under Cold Stress Conditions

Forty-six inbred lines cultivated in Heilongjiang Province, China, were selected and evaluated for cold tolerance based on their seedling survival rate and physiological response (visual observations of the leaves). From the results, 24 inbred line seedlings displayed little phenotypic changes, maintained fully expanded green leaves, had intact plant architecture, strong vigor, and high survival rates (average of 0.8) (Table S2). These inbred lines were therefore classified as cold tolerant and were labeled with an extension of −1 (Figure 1 and Figure S1A). The remaining 22 inbred line seedlings showed visible phenotypic damage, including shriveled, curled, and yellowish spots on the leaves, and low survival rates (average of 0.3) (Table S2). These seedlings were classified as cold sensitive and were labeled with an extension of −2 (Figure 1 and Figure S1B). Collectively, the susceptible lines were more severely damaged by cold stress than the tolerant lines, as evidenced by the shriveled, curled, and yellowish patches on their leaves, as well as low seedling survival rates of about 0.3 on average.

3.2. RNA-Seq Analysis and Alignment of Unique Reads to the Maize Reference Genome

The RNA for the RNA-seq analysis was extracted from the top fully expanding leaves of six-leaf-stage maize seedlings of the 24 tolerant and 22 sensitive maize inbred lines mentioned in Section 2.1 above. The cDNA libraries developed from the RNA described above were constructed and used for Illumina Genome Analyzer (HiSeq™ 2500) deep sequencing. The resulting raw data were deposited into the Genome Sequence Archive under accession number CRA003678 and are publicly accessible at https://bigd.big.ac.cn/gsa, accessed on 5 January 2021). After the filtration of low-quality sequences and adaptor sequences, the 24 tolerant and the 22 sensitive libraries produced 1.13 and 1.12 billion paired-end reads, respectively (Table S3). A total of 2.247 billion paired-end reads with a length of 2 × 150 base pairs (bp) and an average of 48.84 million clean reads per library were obtained (Table S3). The Q20 percentages (sequencing error rates lower than 1%) were more than 95.8%, while the Q30 base percentage, which is an indicator of the overall reproducibility and quality of the assay, was greater than 90%. Moreover, the GC content of each library was 52.6% on average (Table S3). About 1.8 billion clean reads (82%) were mapped to the maize B73_v4 reference genome using HISAT2. Multiple mapped clean reads in each library were excluded from further analysis, and only uniquely mapped clean reads were used for subsequent analysis.

3.3. Identification, Annotation, and Differential Analysis of DEGs

The transcription level was calculated via HTSeq-count as transcript counts, and a total of 53,037 transcripts were obtained and normalized with DESeq2. The expression patterns of these transcripts were investigated in both 24 tolerant and 22 sensitive samples. A definite expression pattern of a single transcript was established after a group comparison analysis (24 tolerant versus 22 sensitive samples), and its base mean, log2 fold change, p-value, and padj values were acquired. From this information, a gene was considered differentially expressed if the p-value was ≤0.05 and the log2 fold change value was ≥1 or ≤−1 between the 24 tolerant and the 22 sensitive samples. A total of 2237 DEGs were detected in all 46 samples. A principal component analysis (PCA) plot was generated based on the gene expression levels of the 46 samples, and the PC1 and PC2 explained 22% of the total variance (Figure 2). The PCA revealed that the correlation of the 46 samples was based on their response to cold stress. Both the tolerant and sensitive samples clustered together, implying a differential response to cold stress. A heatmap of the 2237 DEGs revealed two distinctive clusters, with 1064 DEGs possessing a positive log2 fold change enriched in the tolerant samples, and 1173 DEGs with a negative log2 fold change were enriched in the sensitive samples (Figure 3). Collectively, cold stress upregulated 1064 DEGs in the tolerant lines, while it downregulated them in the sensitive lines. Similarly, cold stress upregulated 1173 DEGs in the sensitive lines, while downregulating them in the tolerant lines.
The annotation of the 2237 DEGs with the maize reference genome B73 RefGen_v4 model resulted in 1656 (74.0%) annotated and 581 (26%) unannotated DEGs. The differential analysis of the 1656 annotated DEGs was carried out based on the two categories of cold-tolerant and cold-sensitive lines. Resultantly, 779 and 877 DEGs were significantly enhanced in the tolerant and sensitive lines, respectively. Further analysis of the 1656, 779, and 877 DEGs was carried out to establish the pathways involved in the cold stress response among the 46 samples.
For the 581 unannotated DEGs, their sequences (≥200 bp) were uploaded to the CPC website for classification as protein-coding or non-coding RNAs. A total of 271 DEGs were classified as protein-coding, while 310 were classified as long non-coding RNAs (lncRNAs) (Figure 4). Moreover, the 581 unannotated DEGs were scanned for the open reading frame (ORF). A total of 402 DEGs with an ORF greater than 120 amino acids (aa) were discarded. The remaining 179 DEGs (ORF length ≤ 120 aa) were aligned to the SWISS-PROT database for the identification of homologous proteins. A total of 56 DEGs were discarded after being homologous to known proteins (E-value ≤ 0.001), while the remaining 123 DEGs were classified as lncRNAs (Figure 4). In total, 337 putative cold-responsive lncRNAs (Figure 4) were identified from the 581 unannotated DEGs, implying the role of lncRNAs in the cold stress response. We analyzed the top-most 20 DEGs regulated by cold stress in both tolerant and sensitive lines and half were unannotated (Table 1). However, most of the annotated DEGs remain uncharacterized, which implies more research is required to unravel the molecular mechanism of cold tolerance.

3.4. Gene Ontology (GO) Analysis of DEGs

To identify the DEGs’ significantly enriched GO terms, the functions of the 1656, 779, and 877 (46 inbred lines, 24 cold-tolerant lines, and 22 cold-sensitive lines, respectively) DEGs were analyzed using AgrigoV2 software. All DEGs were classified into three main GO categories: cellular components, molecular functions, and biological processes. The GO terms related to the response to cold (GO: 0009409), homeostatic process (GO: 0042592), response to temperature stimulus (GO: 0009266), regulation of biological quality (GO: 0065008), response to abiotic stimulus (GO: 0009628), multicellular organismal process (GO: 0032501), response to stress (GO: 0006950), G-protein-coupled receptor protein signaling pathway (GO: 0007186), transcription (GO: 0006350), and cell surface receptor-linked signaling pathway (GO: 0007166) were among the most common significantly enriched terms in the biological process (BP) category of all the three groups of DEGs named above (Figure 5). Within the molecular function (MF) category, catalytic activity (GO: 0003824), protein tyrosine kinase activity (GO: 0004713), G-protein-coupled receptor activity (GO: 0004930), water binding (GO: 0050824), and transcription regulator activity (GO: 0030528) were the most significantly enriched in all the three groups of DEGs (Figure 5). In the cellular component category (CC), integral to membrane (GO: 0016021) and intrinsic to membrane (GO: 0031224) were the common significantly enriched categories (Figure 5). The GO enrichment analysis of the tolerant and sensitive lines (Figure 6) showed that the GO term of metal ion transport (GO: 0030001) was significantly enriched in the tolerant lines (Figure 6A), while the GO terms of apoptosis (GO: 0006915), proteolysis (GO: 0006508), cell death (GO: 0008219), death (GO: 0016265), and programmed cell death (GO: 0012501) were significantly enriched in the sensitive lines (Figure 6B).

3.5. KEGG Pathway Analysis of DEGs

To further explore the biological pathways of the 1656, 779, and 877 DEGs involved in maize response to cold stress, we assessed the number of DEGs in each KEGG pathway. For the 1656 DEGs, the KEGG pathways of lipid metabolism (linoleic acid and arachidonic acid), biosynthesis of other secondary metabolites (isoquinoline alkaloid, betalain, and phenylpropanoid biosynthesis), signal transduction (MAPK signaling pathway-plant), amino acid metabolism (alanine, aspartate, glutamate, glycine, serine, threonine, and phenylalanine metabolism), membrane transport (ABC transporters), terpenoids and polyketide metabolism (carotenoid biosynthesis), and metabolism of other amino acids (glutathione metabolism) were significantly enriched (Figure 7). For the 779 DEGs obtained from the cold-tolerant lines, lipid metabolism (linoleic acid, α-linolenic, ether lipid, arachidonic acid, and glycerophospholipid metabolism), replication and repair (base excision repair), biosynthesis of other secondary metabolites (monobactam and phenylpropanoid biosynthesis), membrane transport (ABC transporters), metabolism of terpenoids and polyketides (carotenoid biosynthesis), and metabolism of other amino acids (glutathione metabolism) were the most significantly enriched pathways (Figure 8A). Within the 887 DEGs from the cold-sensitive lines, the KEGG pathways of biosynthesis of other secondary metabolites (isoquinoline alkaloid biosynthesis), signal transduction (MAPK signaling pathway-plant), transport and catabolism (peroxisome), translation (ribosome), and carbon metabolism were significantly expressed (Figure 8B).

3.6. Dynamic Expression of Signaling and Transcription Factors Genes in Response to Cold Stress

Our GO analysis highlighted a significant number of GO terms related to signaling, such as G-protein-linked receptor protein signaling pathway, cell surface receptor-linked signaling pathway, and protein amino acid phosphorylation (Figure 5 and Figure 6). However, KEGG pathway analysis highlighted the MAPK signaling pathway-plant as one of the most significant pathways (Figure 8B). Stress sensing and signal transduction form crucial adaptive mechanisms in the tolerance of abiotic stresses. Cold stress causes a change in membrane fluidity and cytoskeleton rearrangement, thereby producing signals perceived in the cell membrane by either G-protein-coupled receptors (GPCRs) or osmotic sensors. In this study, 27 DEGs (11 enhanced and 16 suppressed) encoding GPCR were regulated by cold stress (Table S4). Activation of the sensors leads to the generation of reactive oxygen species (ROS), plant hormonal signaling, and cell wall integrity sensing (CWI) [40]. A substantial number of protein kinases, including 21 (11 enhanced and 10 suppressed) leucine-rich repeat protein kinase family proteins, 12 (6 enhanced and 6 suppressed) protein kinase superfamily proteins, 9 (2 enhanced and 7 suppressed) mitogen-activated protein kinases, 7 (2 enhanced and 5 suppressed) S-locus lectin protein kinase family proteins, and 5 wall-associated kinases (WAKs), were regulated by cold stress (Table S4). The high expression of protein kinase suggests that cold stress is primarily regulated at the protein level.
Transcription factors play a vital role in regulating gene expression in response to abiotic stress conditions, such as cold stress. These TFs are DNA-binding proteins that interact with cis-acting elements of genes to activate or inhibit gene transcription, hence regulating plant growth and development and response to the external environment. In the present study, we analyzed for putative TFs in the 1656 common cold-responsive genes based on the 3308 maize TFs and 56 families available in the PlantTFDB 4.0 [41]. Resultantly, 147 TFs, which fell into 32 families, were enriched from our 1656 DEGs (Table 2). ERF (18.37%), WRKY (14.29%), MYB (12.93%), NAC (8.84%), bHLH (8.84%), C2H2 (5.44%) and GRAS (2.72%) were the most abundant TF families (Table 2; Figure 9). However, the expression levels of the TF families MYB, ERF, NAC, WRKY, bHLH, and C2H2 were higher in the sensitive lines than those in the tolerant lines (Figure 9). On the contrary, the TF families MIKC_MADS, CO-like, DBB, Dof, E2F/DP, GATA, GRF, LFY, SBP, Nin-like, TALE, and ZF-HD were only induced in the tolerant lines.
Additionally, GO terms, such as regulation of transcription, transcription, and transcription factor activity, which are related to the regulation of gene expression at the transcription level, were well represented in our GO analysis (Figure 5 and Figure 6). The TF families bHLH, ERF, bZIP, and WRKY were highly regulated by cold stress in the tolerant and sensitive lines (Table S5). As observed above, DEGs encoding the TF families LFY, MIKC_MADS, GRAS, GATA, GRF, E2F/DP, TALE, and Dof were significantly induced by cold stress only in the tolerant lines (Table S5). In-depth research is required to uncover the roles of these TFs in the response of maize to cold stress at the seedling stage. Contrary, DEGs encoding the TF families AP2, B3, HSF, and RAV were all repressed by cold stress in the tolerant lines (Table S5). The tolerant lines were, to a lesser extent, affected by cold stress, which might be why fewer TFs were expressed. Otherwise, our results highlight the critical roles of TFs in the regulation of cold stress response in maize seedlings.

3.7. DEGs Related to “Response to Cold”

The GO terms of “response to cold” and “response to temperature stimulus” were well represented in our GO analysis (Figure 5 and Figure 6). These GO terms contain crucial cold-related genes, which play significant roles in the cold regulation mechanism and may contribute to the cold tolerance of maize seedlings. The expression patterns of DEGs in these GO terms (upregulated or downregulated) may shed light on the difference between the cold response of tolerant and sensitive inbred lines. Signal transduction-related proteins, such as transmembrane proteins, serine/threonine protein kinases, leucine-rich repeat receptor-like proteins, and receptor-like kinases, were regulated by cold stress (Table S6). This highlights the critical roles of the signal transduction network in the activation of various cold-responsive genes. The involvement of TFs in regulating cold stress was highlighted by the significant expression of C2H2, ERF, HD-ZIP, MYB, ERF, MYB, and NAC (Table S6). Antioxidant-related enzymes, such as glutathione peroxidase, thioredoxin, and peroxidase, were regulated by cold stress, suggesting the involvement of detoxification proteins in maize response to cold stress. Carotenoid related genes, such as β-carotene isomerase and β-carotene 3-hydroxylase, further highlighted the role of ABA as a key regulator of cold stress (Table S6). Otherwise, transporter and cell surface proteoglycan-related genes were also observed in this study (Table S6).

3.8. Expression Analysis of Genes Involved in Metabolism, Transport and Functional Impacts of Co-Expressed Gene Hubs

Our KEGG analysis highlighted significant regulation of multiple metabolism pathways, including lipid, carotenoid, ABC transport, and amino acid pathways, during cold stress conditions (Figure 7 and Figure 8). Cold stress activated the phenylpropanoid biosynthetic pathway, an essential way to accumulate various phenolic compounds during cold stress conditions. In total, 14 DEGs encoding this pathway were regulated by cold stress, and 11 of which showed high expression levels in the tolerant lines (Table S7). Moreover, there were six DEGs encoding the metabolism of alanine, aspartate, and glutamate, including β-alanine aminotransferase and glutamate decarboxylase, which are vital genes in osmotic adjustment (Table S7).
Lipid metabolism is a dynamic and complicated process involving lipid biosynthesis, transport, accumulation, turnover, and excretion, regulating the growth and tolerance of plants to different environmental stresses. In this study, 26 lipid metabolism-encoding genes were regulated by cold stress, with all of them except one being enhanced in the tolerant lines (Table S7). Among them were allene oxide synthase (AOS), an essential gene in the biosynthesis of jasmonic acid (JA), secretory phospholipase A2 (PLA2), α/β-Hydrolases, and phospholipase C (PLC), which are critical components of the signaling cascade, and 3-ketoacyl-CoA synthase (KCS) genes, which are related to the biosynthesis of cuticular wax (Table S7). Moreover, six and eight DEGs encoding ABC transporters and carotenoid biosynthesis, respectively, were regulated by cold stress. ABC transporters are associated with phytohormone homeostasis, while carotenoid genes are vital for the biosynthesis of ABA, an essential hormone in cold tolerance. Otherwise, cold stress also regulated glutathione- and peroxisome-related genes, which mediate the harmful effects of ROS. Detailed information about the metabolism genes can be found in Table S7.
WGCNA has been used to dissect the abiotic stress response in plants, thereby highlighting the power of the co-expression networks to provide deep insights into these complex processes. In this study, WGCNA identified multiple significant functional gene hubs related to the cold stress response. A total of 116 critical DEGs that were enriched by GO and KEGG into signaling, TFs, response to cold, and metabolism were defined as hub genes, highlighting the importance of their regulatory impacts on the cold stress response (Figure S2, Table S8). Thus, WGCNA elucidated the higher order relationships between genes based on their co-expression relationships and permitted a robust view of transcriptome organization in the response of cold stress. Collectively, the combination of transcriptome analysis with WGCNA represents an opportunity to achieve a higher resolution analysis that can better predict the most important functional genes that might provide a more robust bio-signature for cold tolerance in maize, thus providing more suitable biomarker candidates for future studies.

3.9. Validation of DEGs by Quantitative Real-Time PCR (qRT-PCR)

To confirm the reliability and validity of the RNA-seq results in maize seedlings, six genes were randomly selected to perform qRT-PCR. Of these, XLOC_018851, Zm00001d002748, and Zm00001d027330 all had a positive log2 fold change between tolerant and sensitive lines, suggesting that their expression was enhanced by cold stress treatment in tolerant lines and declined in sensitive lines. In contrast, XLOC_058556, Zm00001d024324, and Zm00001d024522 portrayed a negative log2 fold change between tolerant and sensitive lines, indicating that cold stress treatment increased their expression in sensitive lines but decreased their expression in tolerant lines. As a result, the fold change ratio used in this paper emphasizes the expression pattern in cold-tolerant lines, while the inverse of that ratio reflects the expression pattern in sensitive lines. We planted B73 and CIMBL116 maize inbred lines, which have previously been reported as cold-sensitive and -tolerant lines, respectively, to validate the expression pattern of these DEGs (Section 2.1 above). These two inbred lines were planted in a cold environment as well as without cold treatment (control), allowing us to compare the expression patterns of our six DEGs before (control) and after cold treatment. Resultantly, cold stress increased the expression of XLOC 018851, Zm00001d027330, and Zm00001d002748 in the tolerant line (CIMBL116), while it decreased the expression of XLOC 058556, Zm00001d024324, and Zm00001d024522 (Figure S3). A reverse expression pattern was observed in the sensitive line (B73) (Figure S3). Thus, the expression trend of our DEGs shown by our RNA-Seq results was in agreement with that shown by the qRT-PCR analyses.
Furthermore, the fold change ratio of the six DEGs between the control and the cold-treated samples was calculated and compared to the fold change obtained from RNA-Seq. As a result, the cold-sensitive line’s RNA-Seq and qRT-PCR results showed an inverse expression pattern (Figure 10A) because the fold change ratio highlighted in this study emphasizes the expression pattern in cold-tolerant lines. However, if you take the inverse of that fold change ratio, which will highlight the expression pattern in cold-sensitive lines, then the RNA-Seq and qRT-PCR of the sensitive line will have a similar expression pattern to that of the tolerant line (Figure 10B). For example, cold stress enhanced the fold change of XLOC_018851, Zm00001d002748, and Zm00001d027330 and declined that of XLOC_058556, Zm00001d024324, and Zm00001d024522 in the tolerant lines (Figure 10B). However, cold stress inversely regulated the fold change of the above-named DEGs in the sensitive lines (Figure 10A). These results validate the authenticity of the DEGs obtained in this study, as the relative fold change in qRT-PCR matched the RNA-Seq results, implying that transcript identification and abundance estimation were remarkably precise. Furthermore, the DEGs identified in this study are universal in maize seedlings of various genetic backgrounds in response to cold stress, according to the qRT-PCR results.

4. Discussion

Maize has surpassed rice and wheat as the world’s most significant cereal crop. However, cold stress affects maize at any stage of development, including the germination, vegetative, and reproductive stages. Screening for cold-tolerant maize cultivars, such as rice, is difficult due to the lack of linkage between cold resistance and developmental periods [42]. Furthermore, cold tolerance is a quantitative trait influenced by the interactions of numerous genes as well as the environment. Nevertheless, RNA-seq research, which evaluates the main genes and regulatory pathways at the transcriptome level, has been widely used to investigate the molecular basis of maize response to abiotic stress [43,44]. To gain a deeper understanding of maize’s cold stress tolerance mechanisms and to develop cold-tolerant maize cultivars, we conducted a comparative transcriptome analysis of 24 cold-tolerant and 22 cold-sensitive maize inbred lines to uncover the common cold-responsive DEGs and pathways. To the best of our knowledge, this is among the first studies to profile a broad set of maize inbred lines at the seedling stage with varying levels of resistance to cold stress in the field. Previous research has usually profiled two samples with differing tolerances to specific environmental stress [45]. As a result, our study provides valuable insight into the in-depth characterization of the molecular responses of maize seedlings to cold stress, because, to date, there remains a scarcity of data on the expression patterns of critical genes and pathways across a gradient of various genotypes with varied responses to specific stress conditions.
Plants perceive abiotic stress via cell wall receptors, which activate internal signaling components through several mechanisms. The cyclic nucleotide–gated channel (CNGC) and glutamate receptors (GLRs) are the primary cell membrane cold receptors reported in plants [46]. These receptors mediate membrane Ca2+ fluxes and produce several endogenous signals responsible for cold tolerance [46]. Multiple G-protein-coupled receptors were shown to be strongly regulated by cold stress in our study (Table S4). The regulation of these receptors may have caused Ca2+ fluxes across membranes and generated multiple signals important for maize’s cold stress response. A G-protein subunit γ gene (Zm00001d032072) and a trihelix GT-2 (Zm00001d027335) were both upregulated in the tolerant lines (Table S4), indicating their role in cold tolerance. In a previous study, transgenic cucumber plants overexpressing CsGG3.2 had increased CBF gene expression and were more tolerant to chilling stress [47]. COLD1 provided chilling tolerance in rice by encoding a signal regulator for guanine nucleotide-binding proteins (G-proteins) on the plasma membrane [48]. A GT2 family (AtGT2L) gene in Arabidopsis has been reported to interact with calcium/calmodulin, allowing plants to withstand cold and salt stress [49]. Two soybean GT2 genes provided abiotic stress resistance to transgenic Arabidopsis plants [50]. Thus, in this study, increased GPCR protein expression may have increased cytosolic calcium levels in the tolerant lines, activating a quick and diverse signaling mechanism responsible for cold acclimation.
Cold stress perceived by membrane sensors triggers an influx of Ca2+ into the cytoplasm, which then generates Ca2+ signatures that induce the activation of downstream genes, such as CBF/COR genes, in the cold signaling pathway [15]. Proteins with an EF-hand domain, such as CaMs, CMLs, CBLs, and CDPKs, act as Ca2+ sensors in response to cold stress [51,52,53]. Following the binding of Ca2+, these proteins interact with other proteins, thereby regulating downstream activities of multiple genes, which provides cold tolerance. In this study, CaMs, CMLs, CBLs, CDPKs, 21 leucine-rich repeat receptor-like kinase proteins (LRR-RLKs), 8 lectin receptor-like kinases (LecRLKs), 12 protein kinases (PKs), and 5 WAKs were regulated by cold stress (Table S4). In Camellia japonica, protein phosphorylation by CDPKs and CIPKs improve cold tolerance [54], whereas in plants, WAKs play an important function in abiotic stress tolerance [55]. Interestingly, TMK1 (Zm00001d033777, Zm00001d007313), CBL10 (Zm00001d023353, Zm00001d010459), RKL1 (Zm00001d048054), NIK3 (Zm00001d018635), and PSKR (Zm00001d018635) displayed a high expression pattern in the tolerant lines. In a previous study on Arabidopsis, CBL10 mediated salt tolerance [56], while RKL1 regulated cold and salicylic acid stresses [57]. Recent research has highlighted the key roles of CBL10 in plant abiotic stress tolerance through the regulation of Na+ and Ca2+ homeostasis [58], whereas under cold stress, TMK1 significantly regulated plant development [59]. MAPKs, including MAPKK and MAPKKK, are key players in cold tolerance [48], where they phosphorylate other kinases and/or various TFs. In this study, 14 DEGs encoding MAPKs were regulated by cold stress (Table S4). Multiple MAPKs have been implicated in improving cold tolerance in rice and Chinese jujube, according to previous research [60,61]. In this study, RBOH (Zm00001d009248), MAP3 (Zm00001d001978), and MKK3 (Zm00001d013510 and Zm00001d028026) had higher expression levels in the tolerant lines (Table S4). In a previous study, two strawberry RBOHs were reported to enhance cold stress tolerance and defense responses [62]. MKK3 was substantially expressed in tolerant lines during a comparative transcriptome investigation of two cotton cultivars with differing responses to cold [63]. The transgenic tobacco over-expressing MAP3K gene demonstrated improved tolerance to a variety of environmental stresses, including cold stress [64]. PP2C adversely controls stress-induced MAPK and SnRK2 protein kinases [4,65]. A previous study on maize has shown that stress-induced proline accumulation and tolerance to hyperosmotic stress were negatively controlled by maize PP2C [66]. This might explain why five PP2C-encoding genes were suppressed by cold stress in the tolerant lines (Table S4). Heptahelical protein 2 (HHP2) (Zm00001d046852), a FASCICLIN-like arabinogalactan (FLA) (Zm00001d016059, Zm00001d052819, and Zm00001d006009), and a membrane-associated kinase (Zm00001d044176) were enhanced in tolerant lines (Table S6). A previous study on Arabidopsis reported that the HHP2-MYB module is involved in integrating cold and abscisic acid signaling to activate the CBF–COR pathway [67]. Cold stress activated a membrane-associated kinase in rice, and FLAs were found to improve banana resistance to low temperatures by activating a cold signal pathway [68,69]. Collectively, the increased expression of Ca2+ signaling protein transcripts in maize tolerant lines activated a complex signaling cascade that regulated various downstream cold tolerance responsive genes. These genes will have significant implications for future research into maize cold tolerance at the seedling stage.
Transcription factors are key regulators of cold stress as they control multiple downstream stress-responsive genes [70]. In this study, 147 TFs belonging to 32 TF families, including AP2/ERF (27), MYB (19), bHLH (13), WRKY (21), C2H2 (8), and NAC (13), were largely regulated by cold stress (Table 2, Figure 5). Similar to our findings, in previous studies, comparative transcriptome analysis of rice, Chinese jujube, and peanut under cold stress conditions identified the above TF families to be the most regulated gene members [60,61,71,72]. These TF genes in their respective families are divided into diverse subgroups based on their specific motif structures, showing that they may perform their specific biological activities under cold stress. Moreover, individual TFs from these families have previously been reported to play a crucial function in controlling plant cold tolerance. In this study, four ERFs, namely ERF38 (Zm00001d002748), ERF022 (Zm00001d048991), DREB26 (Zm00001d018191), and ERF (Zm00001d029669), had higher expression in the tolerant lines (Table S5). They all encode for DREB elements, which enhance cold tolerance by activating CORs. ERF38, ERF022, and DREB26 effectively regulate COR genes and sugar and proline accumulation in Arabidopsis, resulting in abiotic stress tolerance [73,74]. Thus, these genes may have impacted maize cold tolerance by modulating COR genes and osmotic regulators. WRKYs are yet another important class of plant TFs with diverse roles in plant response to cold stress. In this study, 18 WRKY genes were regulated by cold stress (Table S5), with WRKY70 (Zm00001d023332) and WRKY53 (Zm00001d023336) genes showing higher expression in the tolerant lines. In previous studies, the expression of WRKY70 and TcWRKY53 was induced by cold stress in wheat and Thlaspi caerulescens, respectively [75,76]. Moreover, peanut cold stress tolerance was regulated by WRKY70 and WRKY53 via the plant–pathogen interaction pathway [72]. Moreover, MYB4 (Zm00001d041853), bHLH57 (Zm00001d027419), PIF4 (Zm00001d013130), PTF1 (Zm00001d045046), AGL22 (Zm00001d018142), GRF6 (Zm00001d000238), ATHB4 (Zm00001d002754), Scl7 (Zm00001d033834), BTB/POZ (Zm00001d023313, Zm00001d030864), and C2H2 (Zm00001d024883) were all enhanced in the tolerant lines in this study (Tables S5 and S6). SIPIF4 and ZmPTF1 have been attributed to cold and drought stress tolerance in tomatoes and maize, respectively, via the modulation of ABA synthesis and signaling pathways [77,78]. In Arabidopsis thaliana, rice MYB4 was induced by cold stress, which in return transactivated the expression of COR genes, such as RD29A, COR15a, and PAL2 [79]. The over-expression of finger millet bHLH57 caused salinity and drought stress in tobacco [80], while SlGRF6 was significantly regulated by cold stress in Solanum Lycopersicum [81]. C2H2 zinc finger proteins targeted C-repeat/DRE-binding factor genes (CBFs) to provide cold resistance in plants [82], while BTB/POZ significantly accumulated in resistant cotton cultivars during chilling stress conditions [83]. Thus, all these differentially expressed TFs identified in tolerant lines in response to cold stress could represent a useful genetic resource for breeding cold-tolerant crops. Nevertheless, many cold stress-regulating TFs are yet to be identified along with known TFs whose functions are not yet known (Table S5). Understanding the role of the above-mentioned cold-regulating TFs at the molecular level will be pivotal in improving maize performance under cold stress conditions.
Abscisic acid is an essential plant hormone that regulates cold stress via interactions between ABA-dependent and ABA-independent pathways [84]. Moreover, exogenous ABA treatment at normal temperature improves freezing tolerance [85]. In this study, carotenoid biosynthesis genes, such as ZEP (Zm00001d025968), NCED (Zm00001d042076 and Zm00001d018819), β-carotene isomerase (Zm00001d007549 and Zm00001d007560), β-carotene 3-hydroxylase (Zm00001d048469), and ABA 8′-hydroxylase (Zm00001d051554 and Zm00001d050021) were all enhanced by cold stress in the tolerant lines except NCED (Tables S6 and S7). In Arabidopsis, ABA regulates cold tolerance by improving the levels of ABA 8′-hydroxylase [86]. However, Alfalfa-related ZEP is regulated in response to drought, cold, and heat [87]. Moreover, β-carotene hydroxylase regulates the biosynthesis of a carotenoid precursor of abscisic acid called zeaxanthin. A previous report highlighted that a β-carotene hydroxylase gene caused drought and oxidative stress in rice by elevating the synthesis of ABA and xanthophylls [88]. Higher ABA levels induced cold tolerance in herbaceous plants [89]. The elevated expression of ABA-related genes was correlated to cold adaptation in a comparative transcriptome investigation of tea and tobacco plants [90,91]. Therefore, elevated carotenoid biosynthesis genes triggered ABA accumulation, and the transcriptional regulation of ABA-related gene expression is one factor that contributed to cold stress tolerance in maize. Otherwise, the suppression of NCED genes in the tolerant lines reflects the complexity of cold tolerance in plants.
Cold stress triggers rapid and intermittent ROS production that can damage plant cellular components and structures, but ROS also act as signaling molecules for abiotic stress tolerance [92]. Nevertheless, plants deploy a cascade of antioxidant machinery consisting of enzymatic and non-enzymatic defense systems to diminish the deleterious effects of ROS on plant cells [93]. The antioxidant enzymes include SOD, CAT, APX, GPX, GST, and GPX, which can trap and scavenge free radicals [94]. In this study, antioxidant genes, such as SOD (Zm00001d014632), GST (Zm00001d029699, Zm00001d043787, and Zm00001d018809), GPX (Zm00001d029089), PRX (Zm00001d008266, Zm00001d028348, Zm00001d031635, Zm00001d032406 and Zm00001d041827), and APX (Zm00001d024253), were all enhanced by cold stress in the tolerant lines (Table S7). The activities of SOD and GST were reported to reduce cold injury in cold acclimatized wheat [95], while the over-expression of GST in transgenic rice enhanced growth and development at a low temperature [96]. In cassava, chilling and oxidative stress was correlated with increased levels of SOD and APX genes [97]. PRXs play a role in phytoalexin-mediated plant defense and ROS metabolism [98]. Thioredoxin (TRX), however, functions as a redox transmitter [99]. Thus, the enhanced expression of TRX (Zm00001d011352 and Zm00001d007800) genes in the tolerant lines might be crucial in cold acclimation through redox regulation. A previous study reported that soybean TRX genes (Sb03g004670 and Sb06g029490) were significantly regulated in the cold acclimation of different accessions [100]. These findings confirm that ROS-mediated signaling could activate antioxidant enzymes, which might be responsible for imparting cold stress tolerance in maize seedlings.
The phenylpropanoid pathway and its branches of secondary metabolites are activated under cold stress, leading to the accumulation of various phenolic compounds for protection and mechanical support [101]. In this study, 14 phenylpropanoid genes, including CCR (Zm00001d019669, Zm00001d008435), PAL (Zm00001d033286), trans-cinnamate 4-monooxygenase (Zm00001d016471 and Zm00001d032468), and β-glucosidase (Zm00001d028199), and 5 PRXs were significantly enhanced in the tolerant lines (Table S7). Elevated PAL expression stimulates the biosynthesis of phenolic compounds, such as suberin and lignin, which reinforce the cell wall and prevent cell collapse during cold stress. Similarly, enhanced expression of the CCR gene has previously been correlated with lignin biosynthesis under abiotic stress [102]. Trans-cinnamate 4-monooxygenase is positioned at the turning point of phenylalanine, lignin biosynthesis, and flavonoid metabolism, making it one of the key enzymes in the synthesis of lignin and flavonoids [103]. Moreover, PRXs in the presence of H2O2 catalyze the oxidative polymerization of phenols, such as lignin precursors, which improve cell wall rigidity by boosting the cross-linking of cell wall components [104]. This increased suberin or lignin biosynthesis increases the thickness of the cell wall, preventing chilling injury and cell collapse during cold stress [105,106]. The magnitude of lignification in plants is significantly associated with their potential for cold tolerance. Previous studies have shown that β-glucosidase activates several processes, including lignin precursors [107], the release of phytohormones from inactive glycosides, and the activation of several defense compounds essential for abiotic stress tolerance [108]. However, Nicotiana benthamiana plants over-expressing CsBGlu12 displayed abiotic stress tolerance via the accumulation of antioxidant flavanols that played a crucial role in scavenging ROS [109]. Collectively, the upregulation of various transcripts encoding phenylpropanoid pathway genes in the present study indicates enhanced lignification-mediated cold acclimatization in maize seedlings.
A high accumulation of osmoprotectants, such as amino acids, polyamines, quaternary ammonium compounds, and sugars, mediates diverse functions in plant defense mechanisms under varying environmental conditions [110]. Herein, genes encoding β-alanine aminotransferase (Zm00001d038453 and Zm00001d038460) and glutamate decarboxylase (GAD) (Zm00001d031749) were enhanced by cold stress in the tolerant lines (Table S7). The enzyme β-alanine aminotransferase catalyzes the biosynthesis of pyruvate and β-alanine, with the latter product being converted to an essential osmoprotective compound (β-alanine betaine) involved in plant abiotic stress tolerance [111,112]. However, GAD catalyzes the decarboxylation of L-glutamate to form γ-aminobutyric acid (GABA), which accumulates at high concentrations under abiotic stress [113]. Glutamate decarboxylation and GABA metabolism have been reported to play a crucial role in the cold acclimation of wheat and barley [114]. Thus, GABA played a vital role in the cold acclimation of the tolerant lines. Simultaneously, the enhanced expression of β-alanine aminotransferase facilitated a β-alanine-based osmoprotectant in maize during cold stress.
During cold stress, plants adjust their lipid content to retain membrane stability and integrity. Cold tolerance in peanuts was previously found to be associated with changes in membrane modifications, such as lipid metabolism and lipid signaling [115]. In the present study, 26 lipid metabolism genes were regulated by cold stress (Table S7). Among them, PLC (Zm00001d040205), PLA2s (Zm00001d013461, Zm00001d029136), SAD (Zm00001d024273), AOS (Zm00001d028282), α/β-hydrolase (Zm00001d010840 and Zm00001d012147), KCS (Zm00001d046444 and Zm00001d032728), nsLTP (Zm00001d027332), and seven GDSL-like lipases were all enhanced by cold stress in the tolerant lines (Tables S6 and S7). AOS is a critical gene in the synthesis of jasmonic acid (JA), which affects the expression of cold-responsive genes and governs plant defense responses to various abiotic stressors [116]. In Arabidopsis, JA was found to provide cold acclimation [117]. PLC participates in signaling pathways that lead to the activation of the cold response through the CBF pathway [118]. The cold acclimation of spinach (Spinacia oleracea) leaves was found to be associated with the positive roles of PLA2 [119]. Higher expression of a SAD gene is linked to the total amount of unsaturated fatty acids (UFAs), which has been correlated to cold tolerance in tobacco plants [120]. On component change and permeability, the KCS gene catalyzes the biosynthesis of cuticular wax, which acts as a protective barrier against abiotic stresses [121]. GDSL lipases regulate plants’ development and stress response. A previous study by Kong et al. [122] highlighted an essential role of a pepper GDSL lipase gene in regulating abiotic stress tolerance. During cold stress, nsLTPs reduce lipid fluidity and membrane solute permeability, thereby reducing solute diffusion rates across the membrane and preventing osmotic membrane rupture upon thawing [123]. A previous maize study revealed that ZmLTPs have a role in response to cold stress [124]. Overall, our findings highlight the putative association of multiple lipid metabolic components and nsLTP proteins in maize cold tolerance.
Membrane transport systems help maintain cellular homeostasis in environmental stressful situations by redistributing different molecules, such as phytohormones, carbohydrates, and amino acids [125]. These unique roles of plant membrane transport systems may be leveraged to enhance productivity under unfavorable stress conditions as their impact on total plant physiology [126]. The increased expression of numerous transporters and channel protein genes has been reported in the Arabidopsis thaliana response to various abiotic stresses [127] and rice under water stress [128]. In the present study, ABCB1 (Zm00001d024600, Zm00001d025703, Zm00001d026041, Zm00001d045279, and Zm00001d049565), MATE efflux (Zm00001d031730 and Zm00001d032971), and polyol transporter (Zm00001d048774, Zm00001d029645) genes were enhanced by cold stress in the tolerant lines (Tables S6 and S7). Plant ABCB transporters transport molecules, such as plant hormones, lipids, metabolites, contaminants, and defense molecules, which play key roles in abiotic stress tolerance. Various environmental stresses were reported to enhance the expression of distinct ABCB transporters in maize [129]. Polyols play a crucial function in the symplastic and apoplastic transfer of carbon and energy in plants’ response to salt and drought [130]. Transgenic Arabidopsis overexpressing the cotton MATE gene enhanced antioxidant enzyme production and abscisic acid translocation in response to cold, drought, and salt stress [131]. Therefore, the upregulation of various transporters might be associated with cold stress tolerance, the transport of plant secondary metabolites, hormones, and general growth and development in maize.
Network analysis reveals the regulatory impacts of a group of genes on target genes, revealing unique regulatory linkages that add to our understanding of abiotic stress response. The network analysis in this study had 116 nodes and 1907 connections, with 724 activation (positive) and 1183 repression (negative) connections (Figure S2, Table S8). Some of the highly connected positive regulators were found among the 116 nodes, including TFs (bHLH, MYB4, MYB8, GATA4, TALE, and WRKY53) and signaling (respiratory burst oxidase, GPCR, BAM2, RKL1, NIK3, SRF7, SRF8, and PK), antioxidant (peroxiredoxin 6, peroxidase, thioredoxin), and metabolism/biosynthesis regulators (Table S8). The cold induction of TFs regulates a set of other downstream genes. The upregulation of MYB4 (Zm00001d041853) in the cold-tolerant line upregulated 11 DEGs related to signaling, amino acid phosphorylation, TFs, and metabolism (Table S8). In a previous study on Arabidopsis thaliana, the over-expression of OsMYB4 increased cold and chilling tolerance by increasing the expression of COR genes, such as RD29A, COR15a, and PAL2 [79]. In this study, the upregulation of the WRKY53 (Zm00001d023336) gene in the tolerant line increased the expression of 10 additional DEGs involved in signaling, amino acid phosphorylation, and transcription factors (Table S8). In a previous study, WRKY53 was highly increased in Arabidopsis thaliana under cold stress, where it interacted with hub genes, such as mitogen-activated protein kinase 3 (MPK3), WRKY33, and WRKY40, all of which are implicated in plant defense [72]. As a result, WRKY53 could have influenced maize cold tolerance via the plant–pathogen interaction route. Differential expression levels of the 116 DEGs that make up the nodes in our cold studies show that various genes respond to cold stress in different ways and have varied biological functions. These DEGs could be intriguing candidates to investigate during maize seedling cold stress responses. Further research in this regard can look into the molecular specifics of any potential role of these DEGs in the adaptation of maize seedlings to cold stress.
Non-coding RNAs (ncRNAs), such as lncRNAs, have been discovered to regulate plant response to abiotic and biotic stress by controlling the expression of functional genes [132]. A previous study on cassava reported 318 lncRNAs responsive to cold and drought stress [133], while the expression of 2088 lncRNAs in grapevine (Vitis vinifera L.) was induced by cold stress [134]. In this study, the expression of 337 putative lncRNAs was regulated by cold stress (Figure 4). Thus, these lncRNAs might have modulated multiple biological processes involved in cold acclimation in maize by influencing gene expression at the transcriptional, post-transcriptional, and epigenetic levels. Otherwise, there is considerable interest in lncRNAs among molecular biologists, plant breeders, and geneticists, and our research may have identified crucial candidates that can aid in the development of cold-tolerant cultivars. However, more research is needed to fully understand the link between these lncRNAs and cold stress.
We developed a molecular model for cold stress tolerance in maize seedlings, as shown in Figure 11, based on our main findings of the critical cold-responsive DEGs and their associated pathways, as well as the numerous published citations in the present study.

5. Conclusions

In this study, we comprehensively compared the leaf transcriptome and phenotypic response of the maize population (24 cold-tolerant and 22 cold-sensitive lines) in response to cold stress at the seedling stage. Resultantly, the tolerant lines maintained a strong vigor with higher survival rates, while the majority of the sensitive line seedlings died and had yellow spots on the leaves. Using the RNA-seq-based approach, 2237 (1656 annotated and 581 unannotated) DEGs were identified between the tolerant and sensitive samples. Moreover, cold stress significantly enhanced 779 and 887 DEGs in the tolerant and sensitive lines, respectively. Functional annotation was carried out on the three categories (1656, 779, and 887) of DEGs. In the tolerant lines, genes associated with GPCR, Ca2+ signaling, protein kinases, and ROS may have played a significant role in rapid sensing and signaling, whereas genes associated with hormones, such as ABA and JA, may have played a role in signaling and cross-talk between diverse stimuli. The activation of TFs and their binding to promoter sites of certain genes results in activation of stress-responsive genes. The upregulation of several antioxidants, transport, and osmoprotectants suggested protection of the cellular machinery, whereas genes associated with the phenylpropanoid biosynthesis pathway might be involved in providing mechanical support and protection against cold stress. Moreover, genes involved in lipid metabolism may play a critical role in cold stress resistance via membrane modification. Thus, the networks involved in the function of the genes and regulators of the above-named pathways are critical in the cold acclimation of maize at the seedling stage. Moreover, genes related to ribosome, proteolysis, peroxisome, and carbon metabolism were significantly enriched in the sensitive lines. The unannotated DEGs were more inclined in the functions of long non-coding RNAs. Our findings indicate the involvement of plant signaling, transcription factors, and protective mechanisms in the molecular mechanisms underlying cold acclimation in maize at the seedling stage. Otherwise, the essential genes and metabolic pathways identified in this study may serve as valuable genetic resources or selection targets for the genetic engineering of cold-tolerant maize cultivars.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/genes12101638/s1, Figure S1: Phenotypic analysis of maize population under cold stress conditions, Figure S2: Co-expression network analysis of our DEGs, Figure S3: Validation of expression pattern of six genes identified by RNA-Seq by Ct values obtained from qRT-PCR, Table S1: The primers of six differentially expressed genes used for qRT-PCR verification, Table S2: Survival rates of the maize seedlings after cold stress treatment, Table S3: Summary of RNA sequencing results for the forty-six maize seedling leaf samples, Table S4: Expression of signaling related genes, Table S5: Expression of transcription factor-related genes, Table S6: Expression of response to cold-related genes, Table S7: Expression of metabolism/biosynthesis pathway-related genes, Table S8: The 116 critical DEGs identified by WGCNA.

Author Contributions

Conceptualization, J.L., H.W. and J.K.W.; experiments Y.S. (Ying Sun), Y.S. (Yingly Sun), C.L. and J.K.W. Analysis of the data and writing of the original draft J.K.W.; revising and editing the manuscript H.W., J.L., C.Z., Y.S. (Ying Sun), Y.S. (Yingly Sun), Q.C., C.L. and J.K.W. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Key Research and Development Program of China, grant numbers 2018YFD1000702/2018YFD1000700 and 2016YFD0101202; National Natural Science Foundation of China, grant number 31900452; key research and development program of Xinjiang province, China, grant number 2018B01006-4; and the Agricultural Science and Technology Innovation Program of CAAS.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

This manuscript includes the essential data either as figures or as Supplementary Materials. The raw sequence reads have been deposited to the Genome Sequence Archive (GSA) under accession numbers CRA003678 (https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA004123) (accessed on 5 January 2021).

Acknowledgments

We thank Junjie Fu (Institute of Crop Science, Chinese Academy of Agricultural Sciences) for assisting us with B73 and CIMBL116 seeds.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

bHLHBasic helix-loop-helix
DEGsDifferentially expressed genes
MYBMyeloblastosis
KEGGKyoto Encyclopedia of Genes and Genomes
GOGene Ontology
PCAprincipal component analysis
GPCRG-protein coupled receptor
JAjasmonic acid
PPsprotein phosphatase
WAKwall-associated kinase
PALphenylalanine ammonia lyase
KCSketoacyl-CoA synthase
GADglutamate decarboxylase
SODSuperoxide dismutase
APXAscorbate peroxidase
GLRGlutamate receptor
ROSReactive oxygen species
ABAAbscisic acid
PKProtein kinase
PRXPeroxidases
CDPKCalcium dependent protein kinase
CaMsCalmodulin
CBLsCalcineurin B-like
CMLsCaM-related proteins
MAPKMitogen-activated protein kinase

References

  1. United States Department of Agriculture (USDA). World Agricultural Supply and Demand Estimates. 2020. Available online: https://downloads.usda.library.cornell.edu/usda-esmis/files/5q47rn72z/z890sd52n/hd76sk58s/production.pdf (accessed on 12 May 2020).
  2. Gong, F.L.; Yang, F.; Tai, X.; Wang, W. “Omics” of maize stress response for sustainable food production: Opportunities and challenges. Omics A J. Integr. Biol. 2014, 18, 714–732. [Google Scholar] [CrossRef] [Green Version]
  3. Stocker, T.F.; Qin, D.; Plattner, G.K.; Tignor, M.M.; Allen, S.K.; Boschung, J.; Nauels, A.; Xia, Y.; Bex, V.; Midgley, P.M. Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of IPCC the Intergovernmental Panel on Climate Change; Cambridge University Press: Cambridge, UK, 2014. [Google Scholar]
  4. Krasensky, J.; Jonak, C. Drought, salt, and temperature stress-induced metabolic rearrangements and regulatory networks. J. Exp. Bot. 2012, 63, 1593–1608. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Peleg, Z.; Blumwald, E. Hormone balance and abiotic stress tolerance in crop plants. Curr. Opin. Plant Biol. 2011, 14, 290–295. [Google Scholar] [CrossRef]
  6. Ma, S.Q.; Xi, Z.X.; Wang, Q. Risk evaluation of cold damage to corn in Northeast China. J. Nat. Disasters 2003, 12, 137–141. [Google Scholar]
  7. Greaves, J.A. Improving suboptimal temperature tolerance in maize-the search for variation. J. Exp. Bot. 1996, 47, 307–323. [Google Scholar] [CrossRef]
  8. Marocco, A.; Lorenzoni, C.; Fracheboud, Y. Chilling stress in maize. Maydica 2005, 50, 571–580. [Google Scholar]
  9. Sobkowiak, A.; Jończyk, M.; Adamczyk, J.; Szczepanik, J.; Solecka, D.; Kuciara, I.; Hetmańczyk, K.; Trzcinska-Danielewicz, T.; Grzybowski, M.; Skoneczny, M. Molecular foundations of chilling-tolerance of modern maize. BMC Genom. 2016, 17, 125. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  10. Leipner, J.; Stamp, P. Chilling Stress in Maize Seedlings. In Handbook of Maize: Its Biology; Springer: New York, NY, USA, 2009; pp. 291–310. [Google Scholar]
  11. Sowiński, P.; Rudzińska-Langwald, A.; Dalbiak, A.; Sowińska, A. Assimilate export from leaves of chilling-treated seedlings of maize. The path to vein. Plant Physiol. Biochem. 2001, 39, 881–889. [Google Scholar] [CrossRef]
  12. Zhang, B.; Yang, L.; Li, Y. Comparison of physiological and biochemical characteristics related to cold resistance in sugarcane under field conditions. Acta Agron. Sin. 2011, 37, 496–505. [Google Scholar] [CrossRef]
  13. Virdi, A.S.; Singh, S.; Singh, P. Abiotic stress responses in plants: Roles of calmodulin-regulated proteins. Front. Plant Sci. 2015, 6, 809. [Google Scholar] [CrossRef] [Green Version]
  14. Chinnusamy, V.; Zhu, J.; Zhu, J.K. Cold stress regulation of gene expression in plants. Trends Plant Sci. 2007, 12, 444–451. [Google Scholar] [CrossRef] [PubMed]
  15. Chinnusamy, V.; Zhu, J.K.; Sunkar, R. Gene Regulation During Cold Stress Acclimation in Plants. In Plant Stress Tolerance; Springer: Totowa, NJ, USA, 2010; Volume 639, pp. 39–55. [Google Scholar]
  16. Zhao, C.; Wang, P.; Si, T.; Hsu, C.C.; Wang, L.; Zayed, O.; Yu, Z.; Zhu, Y.; Dong, J.; Tao, W.A. MAP kinase cascades regulate the cold response by modulating ICE1 protein stability. Dev. Cell 2017, 43, 618–629. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  17. Shi, Y.; Ding, Y.; Yang, S. Cold signal transduction and its interplay with phytohormones during cold acclimation. Plant Cell Physiol. 2015, 56, 7–15. [Google Scholar] [CrossRef]
  18. Sah, S.K.; Reddy, K.R.; Li, J. Abscisic acid and abiotic stress tolerance in crop plants. Front. Plant Sci. 2016, 7, 571. [Google Scholar] [CrossRef] [Green Version]
  19. Tarkowski, Ł.P.; Van den Ende, W. Cold tolerance triggered by soluble sugars: A multifaceted countermeasure. Front. Plant Sci. 2015, 6, 203. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  20. Arnholdt-Schmitt, B.; Costa, J.H.; de Melo, D.F. AOX—A functional marker for efficient cell reprogramming under stress? Trends Plant Sci. 2006, 11, 281–287. [Google Scholar] [CrossRef] [PubMed]
  21. Al-Whaibi, M.H. Plant heat-shock proteins: A mini review. J. King Saud Univ. Sci. 2011, 23, 139–150. [Google Scholar] [CrossRef] [Green Version]
  22. Alejandro, S.; Lee, Y.; Tohge, T.; Sudre, D.; Osorio, S.; Park, J.; Bovet, L.; Lee, Y.; Geldner, N.; Fernie, A.R. AtABCG29 is a monolignol transporter involved in lignin biosynthesis. Curr. Biol. 2012, 22, 1207–1212. [Google Scholar] [CrossRef] [Green Version]
  23. Wang, B.; Guo, G.; Wang, C.; Lin, Y.; Wang, X.; Zhao, M.; Guo, Y.; He, M.; Zhang, Y.; Pan, L. Survey of the transcriptome of Aspergillus oryzae via massively parallel mRNA sequencing. Nucleic Acids Res. 2010, 38, 5075–5087. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  24. Zhang, Z.; Huang, R. Enhanced tolerance to freezing in tobacco and tomato overexpressing transcription factor TERF2/LeERF2 is modulated by ethylene biosynthesis. Plant Mol. Biol. 2010, 73, 241–249. [Google Scholar] [CrossRef]
  25. Jończyk, M.; Sobkowiak, A.; Trzcinska-Danielewicz, J.; Sowiński, P. Chromatin-Level Differences Elucidate Potential Determinants of Contrasting Levels of Cold Sensitivity in Maize Lines. Plant Mol. Biol. Report. 2021, 39, 335–350. [Google Scholar] [CrossRef]
  26. Grzybowski, M.; Adamczyk, J.; Jończyk, M.; Sobkowiak, A.; Szczepanik, J.; Frankiewicz, K.; Sowiński, P. Increased photosensitivity at early growth as a possible mechanism of maize adaptation to cold springs. J. Exp. Bot. 2019, 70, 2887–2904. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  27. Bolger, A.M.; Lohse, M.; Usadel, B. Trimmomatic: A flexible trimmer for Illumina sequence data. Bioinformatics 2014, 30, 2114–2120. [Google Scholar] [CrossRef] [Green Version]
  28. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Roberts, A.; Pimentel, H.; Trapnell, C.; Pachter, L. Identification of novel transcripts in annotated genomes using RNA-Seq. Bioinformatics 2011, 27, 2325–2329. [Google Scholar] [CrossRef]
  30. Anders, S.; Pyl, P.T.; Huber, W. HTSeq—A Python framework to work with high-throughput sequencing data. Bioinformatics 2015, 31, 166–169. [Google Scholar] [CrossRef] [PubMed]
  31. Anders, S.; Huber, W. Differential expression analysis for sequence count data. Nat. Preced. 2010. [Google Scholar] [CrossRef]
  32. Li, L.; Eichten, S.R.; Shimizu, R.; Petsch, K.; Yeh, C.T.; Wu, W.; Chettoor, A.M.; Givan, S.A.; Cole, R.A.; Fowler, J.E. Genome-wide discovery and characterization of maize long non-coding RNAs. Genome Biol. 2014, 15, R40. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  33. Kong, L.; Zhang, Y.; Ye, Z.Q.; Liu, X.Q.; Zhao, S.Q.; Wei, L.; Gao, G. CPC: Assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007, 35, W345–W349. [Google Scholar] [CrossRef]
  34. Tian, T.; Liu, Y.; Yan, H.; You, Q.; Yi, X.; Du, Z.; Xu, W.; Su, Z. agriGO v2. 0: A GO analysis toolkit for the agricultural community, 2017 update. Nucleic Acids Res. 2017, 45, W122–W129. [Google Scholar] [CrossRef]
  35. Supek, F.; Bošnjak, M.; Škunca, N.; Šmuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 2011, 6, e21800. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  36. Kanehisa, M.; Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000, 28, 27–30. [Google Scholar] [CrossRef] [PubMed]
  37. Xie, C.; Mao, X.; Huang, J.; Ding, Y.; Wu, J.; Dong, S.; Kong, L.; Gao, G.; Li, C.Y.; Wei, L. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011, 39, W316–W322. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  38. Zhang, Y.; Zhu, L.; Xue, J.; Yang, J.; Hu, H.; Cui, J.; Xu, J. Selection and Verification of Appropriate Reference Genes for Expression Normalization in Cryptomeria fortunei Under Abiotic Stress and Hormone Treatments. Genes 2021, 12, 791. [Google Scholar] [CrossRef] [PubMed]
  39. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef] [PubMed]
  40. Novaković, L.; Guo, T.; Bacic, A.; Sampathkumar, A.; Johnson, K.L. Hitting the wall—Sensing and signaling pathways involved in plant cell wall remodeling in response to abiotic stress. Plants 2018, 7, 89. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  41. Jin, J.; Tian, F.; Yang, D.C.; Meng, Y.Q.; Kong, L.; Luo, J.; Gao, G. PlantTFDB 4.0: Toward a central hub for transcription factors and regulatory interactions in plants. Nucleic Acids Res. 2016, 45, gkw982. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Cruz, R.P.D.; Sperotto, R.A.; Cargnelutti, D.; Adamski, J.M.; de FreitasTerra, T.; Fett, J.P. Avoiding damage and achieving cold tolerance in rice plants. Food Energy Secur. 2013, 2, 96–119. [Google Scholar] [CrossRef]
  43. Li, Y.; Wang, X.; Li, Y.; Zhang, Y.; Gou, Z.; Qi, X.; Zhang, J. Transcriptomic Analysis Revealed the Common and Divergent Responses of Maize Seedling Leaves to Cold and Heat Stresses. Genes 2020, 11, 881. [Google Scholar] [CrossRef]
  44. Ji, C.Y.; Kim, H.S.; Lee, C.J.; Kim, S.E.; Lee, H.U.; Nam, S.S.; Kwak, S.S. Comparative transcriptome profiling of tuberous roots of two sweet potato lines with contrasting low temperature tolerance during storage. Gene 2020, 727, 144244. [Google Scholar] [CrossRef]
  45. Yu, T.; Zhang, J.; Cao, J.; Cai, Q.; Li, X.; Sun, Y.; Duan, Y. Leaf transcriptomic response mediated by cold stress in two maize inbred lines with contrasting tolerance levels. Genomics 2021, 113, 782–794. [Google Scholar] [CrossRef] [PubMed]
  46. Zhu, J.K. Abiotic stress signaling and responses in plants. Cell 2016, 167, 313–324. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  47. Bai, L.; Liu, Y.; Mu, Y.; Anwar, A.; He, C.; Yan, Y.; Li, Y.; Yu, X. Heterotrimeric G-protein γ subunit CsGG3. 2 positively regulates the expression of CBF genes and chilling tolerance in cucumber. Front. Plant Sci. 2018, 9, 488. [Google Scholar] [CrossRef] [Green Version]
  48. Guo, X.; Liu, D.; Chong, K. Cold signaling in plants: Insights into mechanisms and regulation. J. Integr. Plant Biol. 2018, 60, 745–756. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  49. Xi, J.; Qiu, Y.; Du, L.; Poovaiah, B.W. Plant-specific trihelix transcription factor AtGT2L interacts with calcium/calmodulin and responds to cold and salt stresses. Plant Sci. 2012, 185, 274–280. [Google Scholar] [CrossRef]
  50. Xie, Z.M.; Zou, H.F.; Lei, G.; Wei, W.; Zhou, Q.Y.; Niu, C.F.; Chen, S.Y. Soybean Trihelix transcription factors GmGT-2A and GmGT-2B improve plant tolerance to abiotic stresses in transgenic Arabidopsis. PLoS ONE 2009, 4, e6898. [Google Scholar] [CrossRef] [Green Version]
  51. Sanders, D.; Pelloux, J.; Brownlee, C.; Harper, J.F. Calcium at the crossroads of signaling. Plant Cell 2002, 14, S401–S417. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Luan, S.; Kudla, J.; Rodriguez-Concepcion, M.; Yalovsky, S.; Gruissem, W. Calmodulins and calcineurin B–like proteins: Calcium sensors for specific signal response coupling in plants. Plant Cell 2002, 14, S389–S400. [Google Scholar] [CrossRef] [Green Version]
  53. Pareek, A.; Khurana, A.; Sharma, K.; Kumar, R. An overview of signaling regulons during cold stress tolerance in plants. Curr. Genom. 2017, 18, 498–511. [Google Scholar] [CrossRef]
  54. Li, Q.; Lei, S.; Du, K.; Li, L.; Pang, X.; Wang, Z.; Xu, L. RNA-seq based transcriptomic analysis uncovers α-linolenic acid and jasmonic acid biosynthesis pathways respond to cold acclimation in Camellia japonica. Sci. Rep. 2016, 6, 36463. [Google Scholar] [CrossRef]
  55. Wu, X.; Bacic, A.; Johnson, K.L.; Humphries, J. The role of brachypodium distachyon wall-associated kinases (WAKs) in cell expansion and stress responses. Cells 2020, 9, 2478. [Google Scholar] [CrossRef] [PubMed]
  56. Kim, B.G.; Waadt, R.; Cheong, Y.H.; Pandey, G.K.; Dominguez-Solis, J.R.; Schültke, S.; Lee, S.C.; Kudla, J.; Luan, S. The calcium sensor CBL10 mediates salt tolerance by regulating ion homeostasis in Arabidopsis. Plant J. 2007, 52, 473–484. [Google Scholar] [CrossRef]
  57. Lee, B.H.; Henderson, D.A.; Zhu, J.K. The Arabidopsis cold-responsive transcriptome and its regulation by ICE1. Plant Cell 2005, 17, 3155–3175. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  58. Martínez, F.A.P.; Estrada, Y.; Flores, F.B.; Ortíz-Atienza, A.; Lozano, R.; Egea, I. The Ca2+ Sensor Calcineurin B-like Protein 10 in Plants: Emerging New Crucial Roles for Plant Abiotic Stress Tolerance. Front. Plant Sci. 2020, 11, 2155. [Google Scholar]
  59. Kamal, M.M.; Ishikawa, S.; Takahashi, F.; Suzuki, K.; Kamo, M.; Umezawa, T.; Uemura, M. Large-Scale Phosphoproteomic Study of Arabidopsis Membrane Proteins Reveals Early Signaling Events in Response to Cold. Int. J. Mol. Sci. 2020, 21, 8631. [Google Scholar] [CrossRef]
  60. Kumar, M.; Gho, Y.S.; Jung, K.H.; Kim, S.R. Genome-wide identification and analysis of genes, conserved between japonica and indica rice cultivars that respond to low-temperature stress at the vegetative growth stage. Front. Plant Sci. 2017, 8, 1120. [Google Scholar] [CrossRef] [Green Version]
  61. Zhou, H.; He, Y.; Zhu, Y.; Li, M.; Song, S.; Bo, W.; Pang, X. Comparative transcriptome profiling reveals cold stress responsiveness in two contrasting Chinese jujube cultivars. BMC Plant Biol. 2020, 20, 240. [Google Scholar] [CrossRef] [PubMed]
  62. Zhang, Y.; Li, Y.; He, Y.; Hu, W.; Zhang, Y.; Wang, X.; Tang, H. Identification of NADPH oxidase family members associated with cold stress in strawberry. FEBS Open Bio 2018, 8, 593–605. [Google Scholar] [CrossRef] [PubMed]
  63. Cheng, G.; Zhang, L.; Wang, H.; Lu, J.; Wei, H.; Yu, S. Transcriptomic Profiling of Young Cotyledons Response to Chilling Stress in Two Contrasting Cotton (Gossypium hirsutum L.) Genotypes at the Seedling Stage. Int. J. Mol. Sci. 2020, 21, 5095. [Google Scholar] [CrossRef] [PubMed]
  64. Kovtun, Y.; Chiu, W.L.; Tena, G.; Sheen, J. Functional analysis of oxidative stress-activated mitogen-activated protein kinase cascade in plants. Proc. Natl. Acad. Sci. USA 2000, 97, 2940–2945. [Google Scholar] [CrossRef] [Green Version]
  65. Schweighofer, A.; Hirt, H.; Meskiene, I. Plant PP2C phosphatases: Emerging functions in stress signaling. Trends Plant Sci. 2004, 9, 236–243. [Google Scholar] [CrossRef] [PubMed]
  66. Liu, L.; Hu, X.; Song, J.; Zong, X.; Li, D.; Li, D. Over-expression of a Zea mays L. protein phosphatase 2C gene (ZmPP2C) in Arabidopsis thaliana decreases tolerance to salt and drought. J. Plant Physiol. 2009, 166, 531–542. [Google Scholar] [CrossRef] [PubMed]
  67. Lee, H.G.; Seo, P.J. The MYB 96–HHP module integrates cold and abscisic acid signaling to activate the CBF–COR pathway in Arabidopsis. Plant J. 2015, 82, 962–977. [Google Scholar] [CrossRef] [PubMed]
  68. Meng, J.; Hu, B.; Yi, G.; Li, X.; Chen, H.; Wang, Y.; Xu, C. Genome-wide analyses of banana fasciclin-like AGP genes and their differential expression under low-temperature stress in chilling sensitive and tolerant cultivars. Plant Cell Rep. 2020, 39, 693–708. [Google Scholar] [CrossRef] [PubMed]
  69. Martín, M.L.; Busconi, L. A rice membrane-bound calcium-dependent protein kinase is activated in response to low temperature. Plant Physiol. 2001, 125, 1442–1449. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  70. Agarwal, M.; Hao, Y.; Kapoor, A.; Dong, C.H.; Fujii, H.; Zheng, X.; Zhu, J.K. A R2R3 type MYB transcription factor is involved in the cold regulation of CBF genes and in acquired freezing tolerance. J. Biol. Chem. 2006, 281, 37636–37645. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  71. Pradhan, S.K.; Pandit, E.; Nayak, D.K.; Behera, L.; Mohapatra, T. Genes, pathways and transcription factors involved in seedling stage chilling stress tolerance in indica rice through RNA-Seq analysis. BMC Plant Biol. 2019, 19, 352. [Google Scholar] [CrossRef] [Green Version]
  72. Jiang, C.; Zhang, H.; Ren, J.; Dong, J.; Zhao, X.; Wang, X.; Yu, H. Comparative Transcriptome-Based Mining and Expression Profiling of Transcription Factors Related to Cold Tolerance in Peanut. Int. J. Mol. Sci. 2020, 21, 1921. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  73. Övernäs, E.; Sundås Larsson, A.; Söderman, E. Two AP2 Transcription Factors, AtERF38 and AtERF39, Have Similar Function but are Differently Regulated in ABA and Stress Responses; Uppsala University Publications: Upsala, Sweden, 2010. [Google Scholar]
  74. Krishnaswamy, S.; Verma, S.; Rahman, M.H.; Kav, N.N. Functional characterization of four APETALA2-family genes (RAP2. 6, RAP2. 6L, DREB19 and DREB26) in Arabidopsis. Plant Mol. Biol. 2011, 75, 107–127. [Google Scholar] [CrossRef]
  75. Wang, J.; Tao, F.; An, F.; Zou, Y.; Tian, W.; Chen, X.; Xu, X.; Hu, X. Wheat transcription factor TaWRKY70 is positively involved in high-temperature seedling plant resistance to Puccinia striiformis f. sp. tritici. Mol. Plant Pathol. 2017, 18, 649–661. [Google Scholar] [CrossRef]
  76. Wei, W.; Zhang, Y.; Han, L.; Guan, Z.; Chai, T. A novel WRKY transcriptional factor from Thlaspi caerulescens negatively regulates the osmotic stress tolerance of transgenic tobacco. Plant Cell Rep. 2008, 27, 795–803. [Google Scholar] [CrossRef] [PubMed]
  77. Wang, F.; Chen, X.; Dong, S.; Jiang, X.; Wang, L.; Yu, J.; Zhou, Y. Crosstalk of PIF4 and DELLA modulates CBF transcript and hormone homeostasis in cold response in tomato. Plant Biotechnol. J. 2020, 18, 1041–1055. [Google Scholar] [CrossRef]
  78. Li, Z.; Liu, C.; Zhang, Y.; Wang, B.; Ran, Q.; Zhang, J. The bHLH family member ZmPTF1 regulates drought tolerance in maize by promoting root development and abscisic acid synthesis. J. Exp. Bot. 2019, 70, 5471–5486. [Google Scholar] [CrossRef] [PubMed]
  79. Vannini, C.; Locatelli, F.; Bracale, M.; Magnani, E.; Marsoni, M.; Osnato, M.; Mattana, M.; Baldoni, E.; Coraggio, I. Overexpression of the rice Osmyb4 gene increases chilling and freezing tolerance of Arabidopsis thaliana plants. Plant J. 2004, 37, 115–127. [Google Scholar] [CrossRef] [PubMed]
  80. Babitha, K.; Vemanna, R.S.; Nataraja, K.N.; Udayakumar, M. Overexpression of EcbHLH57 transcription factor from Eleusine coracana L. in tobacco confers tolerance to salt, oxidative and drought stress. PLoS ONE 2015, 10, e0137098. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  81. Khatun, K.; Robin, A.H.K.; Park, J.I.; Nath, U.K.; Kim, C.K.; Lim, K.B.; Nou, I.S.; Chung, M.Y. Molecular characterization and expression profiling of tomato GRF transcription factor family genes in response to abiotic stresses and phytohormones. Int. J. Mol. Sci. 2017, 18, 1056. [Google Scholar] [CrossRef] [Green Version]
  82. Vogel, J.T.; Zarka, D.G.; Van Buskirk, H.A.; Fowler, S.G.; Thomashow, M.F. Roles of the CBF2 and ZAT12 transcription factors in configuring the low temperature transcriptome of Arabidopsis. Plant J. 2005, 41, 195–211. [Google Scholar] [CrossRef]
  83. Tang, S.; Xian, Y.; Wang, F.; Luo, C.; Song, W.; Xie, S.; Liu, H. Comparative transcriptome analysis of leaves during early stages of chilling stress in two different chilling-tolerant brown-fiber cotton cultivars. PLoS ONE 2021, 16, e0246801. [Google Scholar] [CrossRef]
  84. Agarwal, P.K.; Jha, B. Transcription factors in plants and ABA dependent and independent abiotic stress signaling. Biol. Plant 2010, 54, 201–212. [Google Scholar] [CrossRef]
  85. Chen, J.; Tian, Q.; Pang, T.; Jiang, L.; Wu, R.; Xia, X.; Yin, W. Deep-sequencing transcriptome analysis of low temperature perception in a desert tree, Populus euphratica. BMC Genom. 2014, 15, 326. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  86. Baron, K.N.; Schroeder, D.F.; Stasolla, C. Transcriptional response of abscisic acid (ABA) metabolism and transport to cold and heat stress applied at the reproductive stage of development in Arabidopsis thaliana. Plant Sci. 2012, 188, 48–59. [Google Scholar] [CrossRef] [PubMed]
  87. Zhang, Z.; Wang, Y.; Chang, L.; Zhang, T.; An, J.; Liu, Y.; Cao, Y.; Zhao, X.; Sha, X.; Hu, T. MsZEP, a novel zeaxanthin epoxidase gene from alfalfa (Medicago sativa), confers drought and salt tolerance in transgenic tobacco. Plant Cell Rep. 2016, 35, 439–453. [Google Scholar] [CrossRef] [PubMed]
  88. Du, H.; Wang, N.; Cui, F.; Li, X.; Xiao, J.; Xiong, L. Characterization of the β-carotene hydroxylase gene DSM2 conferring drought and oxidative stress resistance by increasing xanthophylls and abscisic acid synthesis in rice. Plant Physiol. 2010, 154, 1304–1318. [Google Scholar] [CrossRef] [Green Version]
  89. Nayyar, H.; Chander, S. Protective effects of polyamines against oxidative stress induced by water and cold stress in chickpea. J. Agron. Crop Sci. 2004, 190, 355–365. [Google Scholar] [CrossRef]
  90. Li, Y.; Wang, X.; Ban, Q.; Zhu, X.; Jiang, C.; Wei, C.; Bennetzen, J.L. Comparative transcriptomic analysis reveals gene expression associated with cold adaptation in the tea plant Camellia sinensis. BMC Genom. 2019, 20, 624. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  91. Liu, T.; Li, C.X.; Zhong, J.; Shu, D.; Luo, D.; Li, Z.M.; Ma, X.R. Exogenous 1′, 4′-trans-Diol-ABA Induces Stress Tolerance by Affecting the Level of Gene Expression in Tobacco (Nicotiana tabacum L.). Int. J. Mol. Sci. 2021, 22, 2555. [Google Scholar] [CrossRef] [PubMed]
  92. Filiz, E.; Ozyigit, I.I.; Saracoglu, I.A.; Uras, M.E.; Sen, U.; Yalcin, B. Abiotic stress-induced regulation of antioxidant genes in different Arabidopsis ecotypes: Microarray data evaluation. Biotechnol. Biotechnol. Equip. 2019, 33, 128–143. [Google Scholar] [CrossRef] [Green Version]
  93. Baek, K.H.; Skinner, D.Z. Alteration of antioxidant enzyme gene expression during cold acclimation of near-isogenic wheat lines. Plant Sci. 2003, 165, 1221–1227. [Google Scholar] [CrossRef]
  94. Rezaie, R.; Mandoulakani, B.A.; Fattahi, M. Cold stress changes antioxidant defense system, phenylpropanoid contents and expression of genes involved in their biosynthesis in Ocimum basilicum L. Sci. Rep. 2020, 10, 1–10. [Google Scholar] [CrossRef] [PubMed]
  95. Wang, R.; Ma, J.; Zhang, Q.; Wu, C.; Zhao, H.; Wu, Y.; He, G. Genome-wide identification and expression profiling of glutathione transferase gene family under multiple stresses and hormone treatments in wheat (Triticum aestivum L.). BMC Genom. 2019, 20, 986. [Google Scholar] [CrossRef] [Green Version]
  96. Takesawa, T.; Ito, M.; Kanzaki, H.; Kameya, N.; Nakamura, I. Overexpression of ζ glutathione S-transferase in transgenic rice enhances germination and growth at low temperature. Mol. Breed. 2002, 9, 93–101. [Google Scholar] [CrossRef]
  97. Xu, J.; Yang, J.; Duan, X.; Jiang, Y.; Zhang, P. Increased expression of native cytosolic Cu/Zn superoxide dismutase and ascorbate peroxidase improves tolerance to oxidative and chilling stresses in cassava (Manihot esculenta Crantz). BMC Plant Biol. 2014, 14, 208. [Google Scholar] [CrossRef] [PubMed]
  98. Almagro, L.; Gómez Ros, L.V.; Belchi-Navarro, S.; Bru, R.; Ros Barceló, A.; Pedreno, M.A. Class III peroxidases in plant defence reactions. J. Exp. Bot. 2009, 60, 377–390. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  99. Meyer, Y.; Siala, W.; Bashandy, T.; Riondet, C.; Vignols, F.; Reichheld, J.P. Glutaredoxins and thioredoxins in plants. Biochim. Biophys. Acta Mol. Cell Res. 2008, 1783, 589–600. [Google Scholar] [CrossRef]
  100. Ortiz, D.; Hu, J.; Salas Fernandez, M.G. Genetic architecture of photosynthesis in Sorghum bicolor under non-stress and cold stress conditions. J. Exp. Bot. 2017, 68, 4545–4557. [Google Scholar] [CrossRef] [PubMed]
  101. Khaledian, Y.; Maali-Amiri, R.; Talei, A. Phenylpropanoid and antioxidant changes in chickpea plants during cold stress. Russ. J. Plant Physiol. 2015, 62, 772–778. [Google Scholar] [CrossRef]
  102. Srivastava, S.; Vishwakarma, R.K.; Arafat, Y.A.; Gupta, S.K.; Khan, B.M. Abiotic stress induces change in Cinnamoyl CoA Reductase (CCR) protein abundance and lignin deposition in developing seedlings of Leucaena leucocephala. Physiol. Mol. Biol. Plants 2015, 21, 197–205. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  103. Chun-yan, F. Research progress on 4-coumarate: Coenzyme A ligase (4CL) of plants. Mod. Agric. Sci. Technol. 2010, 8, 39–40. [Google Scholar]
  104. Pandey, V.; Awasthi, M.; Singh, S.; Tiwari, S.; Dwivedi, U. A comprehensive review on function and application of plant peroxidases. Biochem. Anal. Biochem. 2017, 6, 308. [Google Scholar] [CrossRef]
  105. Naikoo, M.I.; Dar, M.I.; Raghib, M.F.; Jaleel, H.; Ahmad, B.; Raina, A.; Khan, F.A.; Naushin, F. Role and regulation of plants phenolics in abiotic stress tolerance: An overview. Plant Signal. Mol. 2019, 157–168. [Google Scholar] [CrossRef]
  106. Vogt, T. Phenylpropanoid biosynthesis. Mol. Plant 2010, 3, 2–20. [Google Scholar] [CrossRef] [Green Version]
  107. Dharmawardhana, D.P.; Ellis, B.E.; Carlson, J.E. A [β]-Glucosidase from lodge pole pine xylem specific for the lignin precursor coniferin. Plant Physiol. 1995, 107, 331–339. [Google Scholar] [CrossRef] [PubMed]
  108. Kleczkowski, K.; Schell, J.; Bandur, R. Phytohormones conjugates: Nature and function. Crit. Rev. Plant. Sci. 1995, 14, 283–298. [Google Scholar] [CrossRef]
  109. Baba, S.A.; Vishwakarma, R.A.; Ashraf, N. Functional characterization of CsBGlu12, a β-glucosidase from Crocus sativus, provides insights into its role in abiotic stress through accumulation of antioxidant flavonols. J. Biol. Chem. 2017, 292, 4700–4713. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  110. Slama, I.; Abdelly, C.; Bouchereau, A.; Flowers, T.; Savoure, A. Diversity, distribution and roles of osmoprotective compounds accumulated in halophytes under abiotic stress. Ann. Bot. 2015, 115, 433–447. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  111. Hanson, A.D.; Rathinasabapathi, B.; Chamberlin, B.; Gage, D.A. Comparative physiological evidence that β-alanine betaine and choline-O-sulfate act as compatible osmolytes in halophytic Limonium species. Plant Physiol. 1991, 97, 1199–1205. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  112. Rocha, M.; Licausi, F.; Araujo, W.L.; Nunes-Nesi, A.; Sodek, L.; Fernie, A.R.; van Dongen, J.T. Glycolysis and the tricarboxylic acid cycle are linked by alanine aminotransferase during hypoxia induced by waterlogging of Lotus japonicus. Plant Physiol. 2010, 152, 1501–1513. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  113. Renault, H.; Roussel, V.; El Amrani, A.; Arzel, M.; Renault, D.; Bouchereau, A.; Deleu, C. The Arabidopsis pop2-1 mutant reveals the involvement of GABA transaminase in salt stress tolerance. BMC Plant Biol. 2010, 10, 20. [Google Scholar] [CrossRef] [Green Version]
  114. Mazzucotelli, E.; Tartari, A.; Cattivelli, L.; Forlani, G. Metabolism of γ-aminobutyric acid during cold acclimation and freezing and its relationship to frost tolerance in barley and wheat. J. Exp. Bot. 2006, 57, 3755–3766. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  115. Zhang, H.; Dong, J.I.; Wang, J. Research progress in membrane lipid metabolism and molecular mechanism in peanut cold tolerance. Front. Plant Sci. 2019, 10, 838. [Google Scholar] [CrossRef]
  116. Hu, Y.; Jiang, Y.; Han, X.; Wang, H.; Pan, J.; Yu, D. Jasmonate regulates leaf senescence and tolerance to cold stress: Crosstalk with other phytohormones. J. Exp. Bot. 2017, 68, 1361–1369. [Google Scholar] [CrossRef] [PubMed]
  117. Hu, Y.; Jiang, L.; Wang, F.; Yu, D. Jasmonate regulates the inducer of CBF expression–c-repeat binding factor/DRE binding factor1 cascade and freezing tolerance in Arabidopsis. Plant Cell 2013, 25, 2907–2924. [Google Scholar] [CrossRef] [Green Version]
  118. Vergnolle, C.; Vaultier, M.N.; Taconnat, L.; Renou, J.P.; Kader, J.C.; Zachowski, A.; Ruelland, E. The cold-induced early activation of phospholipase C and D pathways determines the response of two distinct clusters of genes in Arabidopsis cell suspensions. Plant Physiol. 2005, 139, 1217–1233. [Google Scholar] [CrossRef] [Green Version]
  119. Gustavsson, M.H.; Sommarin, M. Characterisation of a plasma membrane-associated phospholipase A2 activity increased in response to cold acclimation. J. Plant Physiol. 2002, 159, 1219–1227. [Google Scholar] [CrossRef]
  120. Craig, W.; Lenzi, P.; Scotti, N.; De Palma, M.; Saggese, P.; Carbone, V.; Curran, N.M.; Magee, A.M.; Medgyesy, P.; Kavanagh, T.A. Transplastomic tobacco plants expressing a fatty acid desaturase gene exhibit altered fatty acid profiles and improved cold tolerance. Transgenic Res. 2008, 17, 769–782. [Google Scholar] [CrossRef] [PubMed]
  121. Bernard, A.; Joubès, J. Arabidopsis Cuticular waxes: Advances in synthesis, export and regulation. Prog. Lipid Res. 2013, 52, 110–129. [Google Scholar] [CrossRef]
  122. Hong, J.K.; Choi, H.W.; Hwang, I.S.; Kim, D.S.; Kim, N.H.; Choi, D.S.; Kim, Y.J.; Hwang, B.K. Function of a novel GDSL-type pepper lipase gene, CaGLIP1, in disease susceptibility and abiotic stress tolerance. Planta 2008, 227, 539–558. [Google Scholar] [CrossRef]
  123. Hincha, D.K.; Neukamm, B.; Sror, H.A.; Sieg, F.; Weckwarth, W.; Rückels, M.; Schmitt, J.M. Cabbage cryoprotectin is a member of the nonspecific plant lipid transfer protein gene family. Plant Physiol. 2001, 125, 835–846. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  124. Wei, K.; Zhong, X. Non-specific lipid transfer proteins in maize. BMC Plant Biol. 2014, 14, 1–18. [Google Scholar] [CrossRef] [Green Version]
  125. Osakabe, Y.; Yamaguchi-Shinozaki, K.; Shinozaki, K.; Tran, L.S.P. ABA control of plant macroelement membrane transport systems in response to water deficit and high salinity. New Phytol. 2014, 202, 35–49. [Google Scholar] [CrossRef]
  126. Schroeder, J.I.; Delhaize, E.; Frommer, W.B.; Guerinot, M.L.; Harrison, M.J.; Herrera-Estrella, L.; Sanders, D. Using membrane transporters to improve crops for sustainable food production. Nature 2013, 497, 60–66. [Google Scholar] [CrossRef] [Green Version]
  127. Seki, M.; Narusaka, M.; Ishida, J.; Nanjo, T.; Fujita, M.; Oono, Y.; Shinozaki, K. Monitoring the expression profiles of 7000 Arabidopsis genes under drought, cold and high-salinity stresses using a full-length cDNA microarray. Plant J. 2002, 31, 279–292. [Google Scholar] [CrossRef]
  128. Chai, C.; Subudhi, P.K. Comprehensive analysis and expression profiling of the OsLAX and OsABCB auxin transporter gene families in rice (Oryza sativa) under phytohormone stimuli and abiotic stresses. Front. Plant Sci. 2016, 7, 593. [Google Scholar] [CrossRef] [Green Version]
  129. Yue, R.; Tie, S.; Sun, T.; Zhang, L.; Yang, Y.; Qi, J.; Yan, S.; Han, X.; Wang, H.; Shen, C. Genome-wide identification and expression profiling analysis of ZmPIN, ZmPILS, ZmLAX and ZmABCB auxin transporter gene families in maize (Zea mays L.) under various abiotic stresses. PLoS ONE 2015, 10, e0118751. [Google Scholar] [CrossRef] [PubMed]
  130. Noiraud, N.; Maurousset, L.; Lemoine, R. Transport of polyols in higher plants. Plant Physiol. Biochem. 2001, 39, 717–728. [Google Scholar] [CrossRef]
  131. Lu, P.; Magwanga, R.O.; Kirungu, J.N.; Hu, Y.; Dong, Q.; Cai, X.; Liu, F. Overexpression of cotton a DTX/MATE gene enhances drought, salt, and cold stress tolerance in transgenic Arabidopsis. Front. Plant Sci. 2019, 10, 299. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  132. Di, C.; Yuan, J.; Wu, Y.; Li, J.; Lin, H.; Hu, L.; Zhang, T.; Qi, Y.; Gerstein, M.B.; Guo, Y. Characterization of stress-responsive lncRNAs in Arabidopsis thaliana by integrating expression, epigenetic and structural features. Plant J. 2014, 80, 848–861. [Google Scholar] [CrossRef] [PubMed]
  133. Li, S.; Yu, X.; Lei, N.; Cheng, Z.; Zhao, P.; He, Y.; Wang, W.; Peng, M. Genome-wide identification and functional prediction of cold and/or drought-responsive lncRNAs in cassava. Sci. Rep. 2017, 7, 45981. [Google Scholar] [CrossRef]
  134. Wang, P.; Dai, L.; Ai, J.; Wang, Y.; Ren, F. Identification and functional prediction of cold-related long non-coding RNA (lncRNA) in grapevine. Sci. Rep. 2019, 9, 6638. [Google Scholar] [CrossRef]
Figure 1. Maize seedling performance under cold stress. The performances were based on seedling survival rates and visual observations of the seedling leaves. The white and black dotted lines represent the cold-tolerant and cold-sensitive lines, respectively. The majority of the tolerant lines maintained a strong vigor with higher survival rates, while most of the sensitive line died, and their leaves had yellowish spots.
Figure 1. Maize seedling performance under cold stress. The performances were based on seedling survival rates and visual observations of the seedling leaves. The white and black dotted lines represent the cold-tolerant and cold-sensitive lines, respectively. The majority of the tolerant lines maintained a strong vigor with higher survival rates, while most of the sensitive line died, and their leaves had yellowish spots.
Genes 12 01638 g001
Figure 2. Principal component analysis (PCA) of pairwise genetic distance. The grouping of the 46 maize inbred line populations is indicated using blue (tolerant) and orange (sensitive). The proportion of variance captured is given as a percentage for both the first and second principal components (PC1 and PC2).
Figure 2. Principal component analysis (PCA) of pairwise genetic distance. The grouping of the 46 maize inbred line populations is indicated using blue (tolerant) and orange (sensitive). The proportion of variance captured is given as a percentage for both the first and second principal components (PC1 and PC2).
Genes 12 01638 g002
Figure 3. Heatmap showing the clustering analysis of 2237 common cold-responsive genes. The x-axis represents different maize samples. The purple color denotes the cold-tolerant lines while the pink color represents the cold-sensitive lines. The red and blue color scale represents high and low expressions, respectively.
Figure 3. Heatmap showing the clustering analysis of 2237 common cold-responsive genes. The x-axis represents different maize samples. The purple color denotes the cold-tolerant lines while the pink color represents the cold-sensitive lines. The red and blue color scale represents high and low expressions, respectively.
Genes 12 01638 g003
Figure 4. Analysis of 581 unannotated cold-responsive DEGs. CPC generated 271 coding and 310 long non-coding RNAs. The ORF method generated 123 long non-coding RNAs. In total, 337 putative long non-coding RNAs were generated by the two methods.
Figure 4. Analysis of 581 unannotated cold-responsive DEGs. CPC generated 271 coding and 310 long non-coding RNAs. The ORF method generated 123 long non-coding RNAs. In total, 337 putative long non-coding RNAs were generated by the two methods.
Genes 12 01638 g004
Figure 5. Gene Ontology (GO) enrichment analysis of the 1656 common cold-responsive genes. The GO terms shown here are the top-most biological process (BP), molecular function (MF), and cellular component (CC) categories.
Figure 5. Gene Ontology (GO) enrichment analysis of the 1656 common cold-responsive genes. The GO terms shown here are the top-most biological process (BP), molecular function (MF), and cellular component (CC) categories.
Genes 12 01638 g005
Figure 6. Gene Ontology (GO) enrichment analysis. (A) The 779 DEGs highly enriched in the tolerant lines. (B) The 877 DEGs highly enriched in the sensitive lines. The GO terms shown here are the top-most biological process (BP), molecular function (MF), and cellular component (CC) categories.
Figure 6. Gene Ontology (GO) enrichment analysis. (A) The 779 DEGs highly enriched in the tolerant lines. (B) The 877 DEGs highly enriched in the sensitive lines. The GO terms shown here are the top-most biological process (BP), molecular function (MF), and cellular component (CC) categories.
Genes 12 01638 g006
Figure 7. KEGG pathway enrichment analysis of the 1656 common cold-responsive genes. The experimental comparisons were based on the hypergeometric test, while the significance of the enrichment of the KEGG pathway was based on q value, q < 0.05. The color gradient represents the size of the q value; the color ranges from green to red, and the closer to green, the smaller the q value and the higher the significant degree of enrichment of the corresponding KEGG pathway. The “rich factor” represents the percentage of DEGs to total genes in a given pathway.
Figure 7. KEGG pathway enrichment analysis of the 1656 common cold-responsive genes. The experimental comparisons were based on the hypergeometric test, while the significance of the enrichment of the KEGG pathway was based on q value, q < 0.05. The color gradient represents the size of the q value; the color ranges from green to red, and the closer to green, the smaller the q value and the higher the significant degree of enrichment of the corresponding KEGG pathway. The “rich factor” represents the percentage of DEGs to total genes in a given pathway.
Genes 12 01638 g007
Figure 8. KEGG pathway enrichment analysis. (A) The 779 DEGs highly enriched in the tolerant lines. (B) The 877 DEGs highly enriched in the sensitive lines. The experimental comparisons were based on the hypergeometric test, while the significance of the enrichment of the KEGG pathway was based on q value, q < 0.05. The color gradient represents the size of the q value; the color ranges from green to red, and the closer to green, the smaller the q value and the higher the significant degree of enrichment of the corresponding KEGG pathway. The “rich factor” represents the percentage of DEGs to total genes in a given pathway.
Figure 8. KEGG pathway enrichment analysis. (A) The 779 DEGs highly enriched in the tolerant lines. (B) The 877 DEGs highly enriched in the sensitive lines. The experimental comparisons were based on the hypergeometric test, while the significance of the enrichment of the KEGG pathway was based on q value, q < 0.05. The color gradient represents the size of the q value; the color ranges from green to red, and the closer to green, the smaller the q value and the higher the significant degree of enrichment of the corresponding KEGG pathway. The “rich factor” represents the percentage of DEGs to total genes in a given pathway.
Genes 12 01638 g008
Figure 9. Area map of the various transcription factor families regulated by cold stress. The x-axis represents the TF families, while the number of genes per family is represented by the y-axis. The blue and orange colors indicate the TFs regulated in the tolerant and sensitive lines, respectively.
Figure 9. Area map of the various transcription factor families regulated by cold stress. The x-axis represents the TF families, while the number of genes per family is represented by the y-axis. The blue and orange colors indicate the TFs regulated in the tolerant and sensitive lines, respectively.
Genes 12 01638 g009
Figure 10. Validation of RNA-Seq results by qRT-PCR. Each log2 fold change calculated from qRT-PCR was compared with the log2 fold change of the RNA-Seq data. (A) The inverse expression patterns of the RNA-Seq and qRT-PCR results from the cold-sensitive line (B73) is because the fold change ratio highlighted in this study emphasizes the expression pattern in cold-tolerant lines. However, if you take the inverse of that fold change ratio, which will highlight the expression pattern in cold-sensitive lines, then the RNA-Seq and qRT-PCR expression trends would be identical. (B) The cold-tolerant line’s (CIMBL116) RNA-Seq and qRT-PCR results show a similar expression trend. Orange and blue bars represent RNA-Seq and qRT-PCR data, respectively.
Figure 10. Validation of RNA-Seq results by qRT-PCR. Each log2 fold change calculated from qRT-PCR was compared with the log2 fold change of the RNA-Seq data. (A) The inverse expression patterns of the RNA-Seq and qRT-PCR results from the cold-sensitive line (B73) is because the fold change ratio highlighted in this study emphasizes the expression pattern in cold-tolerant lines. However, if you take the inverse of that fold change ratio, which will highlight the expression pattern in cold-sensitive lines, then the RNA-Seq and qRT-PCR expression trends would be identical. (B) The cold-tolerant line’s (CIMBL116) RNA-Seq and qRT-PCR results show a similar expression trend. Orange and blue bars represent RNA-Seq and qRT-PCR data, respectively.
Genes 12 01638 g010
Figure 11. The schematic molecular model describing the signaling pathways involved in the acquisition of cold tolerance in maize seedlings. The model was constructed based on the main cold response components identified in this report, as well as plant abiotic stress pathway schemes previously described. The downward pointing arrows represent the sequence of events in cold tolerance in maize, from stress signal perception to acclimation mechanisms. Abbreviation key: GPCR, G-protein-coupled receptor; GLR, glutamate receptor; ROS, reactive oxygen species; ABA, abscisic acid; JA, jasmonic acid; CDPKs, calcium-dependent protein kinases; MAPK, mitogen-activated protein kinase; PPs, protein phosphatase; WAKs, wall-associated kinases; PKs, protein kinases; PRX, peroxidases; PAL, phenylalanine ammonia lyase; KCS, ketoacyl-CoA synthase; GAD, glutamate decarboxylase; SOD, superoxide dismutase; APX, ascorbate peroxidase; TRX, thioredoxin; GST, glutathione transferase; GPX, glutathione peroxidase; ZEP, zeaxanthin epoxidase.
Figure 11. The schematic molecular model describing the signaling pathways involved in the acquisition of cold tolerance in maize seedlings. The model was constructed based on the main cold response components identified in this report, as well as plant abiotic stress pathway schemes previously described. The downward pointing arrows represent the sequence of events in cold tolerance in maize, from stress signal perception to acclimation mechanisms. Abbreviation key: GPCR, G-protein-coupled receptor; GLR, glutamate receptor; ROS, reactive oxygen species; ABA, abscisic acid; JA, jasmonic acid; CDPKs, calcium-dependent protein kinases; MAPK, mitogen-activated protein kinase; PPs, protein phosphatase; WAKs, wall-associated kinases; PKs, protein kinases; PRX, peroxidases; PAL, phenylalanine ammonia lyase; KCS, ketoacyl-CoA synthase; GAD, glutamate decarboxylase; SOD, superoxide dismutase; APX, ascorbate peroxidase; TRX, thioredoxin; GST, glutathione transferase; GPX, glutathione peroxidase; ZEP, zeaxanthin epoxidase.
Genes 12 01638 g011
Table 1. Transcription factor gene families identified from 2237 DEGs in maize under cold stress.
Table 1. Transcription factor gene families identified from 2237 DEGs in maize under cold stress.
TF FamilyDEGs in the Tolerant LineDEGs in the Sensitive LineTotal
B3022
GRF101
ERF32528
DBB101
Dof101
HSF033
LBD213
LFY101
MYB71219
NAC21113
RAV011
SBP101
WOX011
TCP011
E2F/DP101
GATA101
GRAS224
bZIP123
C2H2178
bHLH8513
Nin-like101
NF-YC011
ZF-HD101
G2-like112
CO-like101
HD-ZIP224
WRKY21921
TALE101
Trihelix112
MIKC_MADS303
MYB_related134
Total47100147
Table 2. List of top 20 most regulated DEGs by cold stress.
Table 2. List of top 20 most regulated DEGs by cold stress.
Locus IDGene IDlog2 Fold Change (T/S)p_ValueChrStartEndAnnotation
XLOC_000983Zm00001d02760623.734735126.81 × 10−31Chr18,915,7198,918,057transmembrane protein
XLOC_056725-−23.73433019.38 × 10−29Chr7170,162,017170,164,579-
XLOC_046285Zm00001d017622−2.6700165399.74 × 10−25Chr5201,997,628201,998,691OSJNBa0088A01
XLOC_018249-24.010631951.05 × 10−24Chr2137,103,173137,104,433-
XLOC_058213Zm00001d021006−1.5549840723.70 × 10−22Chr7140,181,294140,183,018MTD1
XLOC_034912-23.02398335.49 × 10−21Chr4176,582,912176,584,805-
XLOC_056365Zm00001d02139422.791566517.11 × 10−21Chr7150,891,109150,895,395hypothetical protein
XLOC_067464-20.058734116.49 × 10−20Chr920,669,45520,671,879-
XLOC_049980-24.720083892.60 × 10−17Chr6114,599,197114,601,220-
XLOC_018351-24.47397645.35 × 10−17Chr2150,920,460150,923,282-
XLOC_055724-22.966325459.06 × 10−17Chr795,897,27695,897,694-
XLOC_018159Zm00001d00462022.67408871.68 × 10−16Chr2122,442,335122,447,585uncharacterized protein
XLOC_004771-−23.949512632.29 × 10−16Chr11,487,9831,491,882-
XLOC_007853Zm00001d033411−23.881069342.77 × 10−16Chr1262,279,726262,300,105hypothetical protein
XLOC_005508Zm00001d02867323.879427452.96 × 10−16Chr142,581,89842,586,482small nuclear protein G
XLOC_036989-−23.842167273.11 × 10−16Chr4108,230,941108,236,238-
XLOC_014865Zm00001d02596823.60897416.37 × 10−16Chr1134,994,735134,998,645hypothetical protein
XLOC_018011Zm00001d00433823.583045466.85 × 10−16Chr2104,459,685104,460,441hypothetical protein
XLOC_046112Zm00001d017287−4.3174815182.92 × 10−13Chr5191,750,113191,750,625uncharacterized protein
XLOC_027574Zm00001d043525−1.3502597534.09 × 10−13Chr3202,741,852202,743,157oxidative stress 3
Comparison of DEGs between tolerant (T) and sensitive (S) inbred lines after cold stress. The p-value is less than 0.05.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Waititu, J.K.; Cai, Q.; Sun, Y.; Sun, Y.; Li, C.; Zhang, C.; Liu, J.; Wang, H. Transcriptome Profiling of Maize (Zea mays L.) Leaves Reveals Key Cold-Responsive Genes, Transcription Factors, and Metabolic Pathways Regulating Cold Stress Tolerance at the Seedling Stage. Genes 2021, 12, 1638. https://doi.org/10.3390/genes12101638

AMA Style

Waititu JK, Cai Q, Sun Y, Sun Y, Li C, Zhang C, Liu J, Wang H. Transcriptome Profiling of Maize (Zea mays L.) Leaves Reveals Key Cold-Responsive Genes, Transcription Factors, and Metabolic Pathways Regulating Cold Stress Tolerance at the Seedling Stage. Genes. 2021; 12(10):1638. https://doi.org/10.3390/genes12101638

Chicago/Turabian Style

Waititu, Joram Kiriga, Quan Cai, Ying Sun, Yinglu Sun, Congcong Li, Chunyi Zhang, Jun Liu, and Huan Wang. 2021. "Transcriptome Profiling of Maize (Zea mays L.) Leaves Reveals Key Cold-Responsive Genes, Transcription Factors, and Metabolic Pathways Regulating Cold Stress Tolerance at the Seedling Stage" Genes 12, no. 10: 1638. https://doi.org/10.3390/genes12101638

APA Style

Waititu, J. K., Cai, Q., Sun, Y., Sun, Y., Li, C., Zhang, C., Liu, J., & Wang, H. (2021). Transcriptome Profiling of Maize (Zea mays L.) Leaves Reveals Key Cold-Responsive Genes, Transcription Factors, and Metabolic Pathways Regulating Cold Stress Tolerance at the Seedling Stage. Genes, 12(10), 1638. https://doi.org/10.3390/genes12101638

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop