Comprehensive Transcriptome Analysis of Responses during Cold Stress in Wheat (Triticum aestivum L.)

Wheat production is often impacted by pre-winter freezing damage and cold spells in later spring. To study the influences of cold stress on wheat seedlings, unstressed Jing 841 was sampled once at the seedling stage, followed by 4 °C stress treatment for 30 days and once every 10 days. A total of 12,926 differentially expressed genes (DEGs) were identified from the transcriptome. K-means cluster analysis found a group of genes related to the glutamate metabolism pathway, and many genes belonging to the bHLH, MYB, NAC, WRKY, and ERF transcription factor families were highly expressed. Starch and sucrose metabolism, glutathione metabolism, and plant hormone signal transduction pathways were found. Weighted Gene Co-Expression Network Analysis (WGCNA) identified several key genes involved in the development of seedlings under cold stress. The cluster tree diagram showed seven different modules marked with different colors. The blue module had the highest correlation coefficient for the samples treated with cold stress for 30 days, and most genes in this module were rich in glutathione metabolism (ko00480). A total of eight DEGs were validated using quantitative real-time PCR. Overall, this study provides new insights into the physiological metabolic pathways and gene changes in a cold stress transcriptome, and it has a potential significance for improving freezing tolerance in wheat.


Introduction
Wheat (Triticum aestivum L.) is one of the world's most important food crops and possesses enormous economic significance [1]. Under climate change, cold events are likely to continue or even become more serious [2]. To survive, grow, and reproduce in the fluctuating environment yearly, plants have evolved a number of complex mechanisms to respond to abiotic stresses. The transcription of stress genes regulates the adaptation process, including developmental, physiological, and biochemical changes [3][4][5]. The stress genes are currently known to encode two main categories of proteins: one category is mainly involved in plant protection, while the other category is involved in regulating the expression of downstream genes. This second category contains transcription factors (TFs) that can bind to the promoter regions of stress-responsive genes, thereby improving stress tolerance in plants [6][7][8]. Through the analysis of stress response promoters, many cis-acting and trans-acting elements involved in the transcription of stress response genes have been identified [4,9].
Temperature change is a major seasonal change in a temperate climate. For plants, it is very important to have a mechanism to survive the freezing temperatures that occur in winter. Abiotic stress factors, such as cold and frost, severely affect the growth and development of plants, including wheat, ultimately affecting crop yields [10]. The freezing resistance of plants refers to their ability to survive and continue to function after being mechanism of cold tolerance in winter wheat, and to provide important target genes for molecular design and breeding.

Plant Materials and Cold Stress Treatment
The wheat Jing 841 used in this study is a winter variety [37]. The surface of seeds were disinfected with 75% ethanol for 1 min, then they were disinfected with 10% sodium hypochlorite for 15 min, and they were finally washed with sterile water 5 times. The seeds were then planted into Petri dishes with wet filter paper and placed in an incubator at 22 • C under 16 h light and 8 h darkness cycle. After 10 days, the seedlings were moved to another condition. The temperature of the cold stress treatment was kept at 4 • C, and the control condition was 22 • C.

Trait Measurements and Morphology Observations
The samples were observed while the length of the seedlings was about 8-10 cm, and the pictures of the seedlings were photographed using a digital camera (Nikon Coolpix 4500, Nikon Corporation, Tokyo, Japan).

Library Preparation, RNA Extraction, and Sequencing
The total RNA of the four samples was prepared using RNeasy Plant Mini Kit (QIA-GEN, Hilden, Germany) according to the manufacturer's protocol, including the control group (CK: CK0-1, CK0-2, and CK0-3), seedlings after cold stress treated for 10 days (D10: D10-1, D10-2, and D10-3), seedlings after cold stress treated for 20 days (D20: D20-1, D20-2, and D20-3), and seedlings after cold stress treated for 30 days (D30: D30-1, D30-2, and D30-3). Each sample had three biological replicates, and each sample contained more than 20 individuals. All the samples were immediately frozen in liquid nitrogen for the subsequent experiments. The RNA concentration and integrity were checked using the NanoDrop 2000 and the RNA Nano 6000 Assay Kit (NanoDrop Technologies, Wilmington, DE, USA). Qualified samples were used for library construction and Illumina sequencing.
The twelve cDNA libraries were sequenced using Illumina HiSeq X Ten from Biomarker Biotechnology Corporation (Beijing, China) [38].

Transcriptome Analyses
The clean reads were compared with the wheat reference genome IWGSC RefSeq v1.0 (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/, 1 January 2023) using HISAT2 [39]. The unigenes were annotated based on Nr (NCBI non-redundant protein sequences); Nt (NCBI non-redundant nucleotide sequences). The mapped reads were pieced using StringTie software and compared with the original genome annotation information, in search of the original unannotated transcription regions. KEGG orthology-based annotation of the sequence's genes was analyzed using KOBAS2.0 [40], and the amino acid sequences of the new genes were predicted. After predicting the amino acid sequences of the new genes, HMMER [41] software was used to compare with the Pfam database to complement and refine the original genome annotation information. The gene expression levels were quantified as fragments per kilobase of the transcript per million fragments mapped (FPKM) [42]. To identify differentially expressed genes (DEGs), the DEseq method was used with a false discovery rate (FDR) threshold of <0.01 and a fold change value (FC) threshold of ≥8 of significance for differential expression. Biological repetition correlation was evaluated using the Pearson correlation coefficient (PCC) and the principal components analysis (PCA) [43,44]. DEGs were screened with a threshold parameter of Log 2 FC ≥ 3 or ≤−3, and gene co-expression analysis was performed on the expression values using the K-means clustering algorithm [45]. A gene was considered valid gene if its average FPKM value among the twelve libraries was more than one. The statistical analyses of genes in different samples were performed using valid genes. All the analyses were conducted on BMKCloud (https://www.biocloud.net, accessed on 1 January 2023). The accession number of the transcriptome data in NCBI is PRJNA903114.

K-Means Clustering of DEGs
In this study, the expression values of three replicates of each sample were used for gene co-expression analysis using the K-means clustering algorithm [45]. A K-means clustering method was used to cluster the samples in Section 2.4, and the number of clusters was 8. We used log 2 (FPKM+1) and expression values normalized by the Z-score to create a heat map, and we used heat map in TBtools for the visualization.

Weighted Gene Co-Expression Network Analysis (WGCNA)
WGCNA was used to construct gene co-expression networks by analyzing pairwise correlations among genes with similar expression trends. Using the WGCNA tool, gene co-expression networks were constructed based on the pairwise correlations of gene coexpression trends in all the sampled tissues. The R package WGCNA converts adjacency into topological overlap, which measures the network connectivity of a gene by summing its adjacency with other genes. The hierarchical clustering was used to classify genes with similar expression profiles into modules based on TOM dissimilarity. At least 30 modules are used to create a gene dendrogram. Tangents were selected for merging modules by calculating the dissimilarity of the feature genes in each module. Hub genes, which are highly interconnected with nodes in a module, were considered functionally significant. Key genes were selected from the DEGs in each module, and the DEGs in the key modules were screened for further analysis [47,48].

qRT-PCR Validation
Four samples of Jing 841 with four treatments (CK, D10, D20, and D30) were prepared for qRT-PCR. Primers were designed using Primer Premier 5.0 software (www.premierbiosoft.com/primerdesign/index.html, accessed on 1 January 2023). The primer sequence is listed in Table S1. The TranScript ®® integrated First-Strand cDNA Synthesis SuperMix for qPCR was used for reverse transcription (Transgenic Biotechnology, Beijing, China). The TransStart ®® Top Green qPCR SuperMix (2×) was employed for qRT-PCR in accordance with the manufacturer's instructions on the CFX ConnectTM Real-Time System (Bio-Rad, Hercules, CA, USA). The configuration of qRT-PCR is 20 µL volumes [38]. The 2 −∆∆Ct method was utilized to calculate the gene expression level, using the wheat actin gene as an internal control [49].

Transcriptome Data Overview
A total of 114.57 Gb clean data and about 382 million reads were obtained from the twelve libraries, averaging 8.59 Gb and 31 million reads per sample. The G + C content of all the samples varied from 56.13% to 57.97%, and the average Q30 percentage was 94.73% (Table S2). The percentage of unique mapped reads also exceeded 94.70% (Table S3). Based on the mapped results, 22,928 new genes were discovered, of which 15,242 have functional annotations (Table S4).
A good correlation was clearly shown among the transcriptome sequencing samples, as shown in the three-dimensional map of PCA ( Figure 1). The colored points on the PCA map represented the distinct samples, and the three eigenvectors (22.25%, 14.90%, 10.59%) distinctly separated them, which was consistent with the morphological pattern of the different tissues. Furthermore, PCC analysis demonstrated that each correlation coefficient exceeded 0.98 among all the replicated samples, indicating that the transcriptome sequencing reads were of a good quality for the subsequent identification of DEGs and regulatory network construction. on the mapped results, 22,928 new genes were discovered, of which 15,242 have functional annotations (Table S4).
A good correlation was clearly shown among the transcriptome sequencing samples, as shown in the three-dimensional map of PCA ( Figure 1). The colored points on the PCA map represented the distinct samples, and the three eigenvectors (22.25%, 14.90%, 10.59%) distinctly separated them, which was consistent with the morphological pattern of the different tissues. Furthermore, PCC analysis demonstrated that each correlation coefficient exceeded 0.98 among all the replicated samples, indicating that the transcriptome sequencing reads were of a good quality for the subsequent identification of DEGs and regulatory network construction.  The transcriptional network of Jing 841 was established, and the key genes regulating cold stress treatment were identified. From the twelve libraries, 133,718 genes (unigenes) were obtained, and 125,405 unigenes were BLAST-annotated in multiple databases (Table S5). In the four samples, 53,625, 51,516, 51,609, and 48,193 valid genes were identified ( Figure 2a). The number of valid genes only expressed in CK, CK-D10, CK-D20, and CK-D30 was 4542, 452, 1357, and 3044, respectively ( Figure 2b). The number of DEGs was the least between CK-D10, while the number of DEGs was the most between CK-D30.

Identification of DEGs
DEGs were identified through a pairwise comparison of the twelve libraries. When the filtering threshold remained unchanged at FDR < 0.01, the number of identified DEGs in the four tissues subjected to the four different treatments varied: 71,557, 33,231, 19,174, and 12,364. The number of DEGs decreased as the fold change increased. (Figure 2c).
To identify key genes related to cold stress development, significant DEGs (Log 2 FC ≥ 3 or ≤3) were screened by comparing samples treated with cold stress at different times with CK. In the sample pairs of CK-D10, CK-D20, and CK-D30, a total of 3108, 4197, and 11,869 DEGs (Table S6 and Figure 2c) were identified (FDR < 0.01 and FC ≥ 8). TF genes have long been an important research focus as key regulatory components in the regulation of gene expression. In this study, we annotated TFs based on specific DNA-binding domain signatures from the Plant-TFDB database (http://planttfdb.cbi.pku.edu.cn, accessed on 1 January 2023) using RNA-seq. We detected 854 unique differentially expressed TFs among the DEGs (Table S7). It showed that the transcription of TFs is important in the wheat cold stress response. Among the TFs of differential genes during cold stress treatment for 10 days, the bHLH-, C2H2-, and MYB-related families seem to play a more important role in the cold response, because there were 21, 17, and 11 in these 3 families, respectively. Each member has a differential expression. Among the TFs of DEGs during the cold stress treatment, the bHLH, MYB, and NAC families may play a relatively important role in responding to cold stress because each of these 3 families have 73, 59, and 56 members with a differential expression ( Figure 3).

Identification of DEGs
DEGs were identified through a pairwise comparison of the twelve libraries. When the filtering threshold remained unchanged at FDR < 0.01, the number of identified DEGs in the four tissues subjected to the four different treatments varied: 71,557, 33,231, 19,174, and 12,364. The number of DEGs decreased as the fold change increased. (Figure 2c).
To identify key genes related to cold stress development, significant DEGs (Log2FC ≥ 3 or ≤3) were screened by comparing samples treated with cold stress at different times with CK. In the sample pairs of CK-D10, CK-D20, and CK-D30, a total of 3108, 4197, and 11,869 DEGs (Table S6 and Figure 2c) were identified (FDR < 0.01 and FC ≥ 8). TF genes have long been an important research focus as key regulatory components in the regulation of gene expression. In this study, we annotated TFs based on specific DNA-binding domain signatures from the Plant-TFDB database (http://planttfdb.cbi.pku.edu.cn, accessed on 1 January 2023) using RNA-seq. We detected 854 unique differentially expressed TFs among the DEGs (Table S7). It showed that the transcription of TFs is important in the wheat cold stress response. Among the TFs of differential genes during cold stress treatment for 10 days, the bHLH-, C2H2-, and MYB-related families seem to play a more important role in the cold response, because there were 21, 17, and 11 in these 3 families, respectively. Each member has a differential expression. Among the TFs of DEGs during the cold stress treatment, the bHLH, MYB, and NAC families may play a relatively important role in responding to cold stress because each of these 3 families have 73, 59, and 56 members with a differential expression ( Figure 3). DEGs were classified into 45 out of 52 subcategories using the GO database ( Figure  S1), while the KEGG database classified these genes into six major metabolic pathways.
In addition, the further enrichment and analysis of DEGs was performed. We present the top enrichment categories of the DEGs in biological processes, cellular components, and molecular functions in Table 1. In the comparison between the three groups of differ- DEGs were classified into 45 out of 52 subcategories using the GO database ( Figure S1), while the KEGG database classified these genes into six major metabolic pathways.
In addition, the further enrichment and analysis of DEGs was performed. We present the top enrichment categories of the DEGs in biological processes, cellular components, and molecular functions in Table 1. In the comparison between the three groups of differently treated samples and the control group, the top enrichment terms of DEGs in the biological process, cell composition, and molecular function were the same in each group. They were nucleosome assembly (GO:0006334), nucleosome (GO:0000786), and protein heterodimerization activity (GO:0046982). The nucleosome assembly was very important to reconstruct the functional chromatin structure, which affected many biological processes, including DNA replication, repair, recombination, transcription, cell proliferation and differentiation, and individual development. These enrichment results suggested that most genes were associated with many biological processes, such as individual development. The top metabolic pathways (Table 2) were classified based on gene annotation in the KEGG database. The DEGs of CK-D10 involved the top three enriched pathways that were photosynthesis-antenna proteins (ko00196), starch and sucrose metabolism (ko00500), and phenylpropanoid biosynthesis (ko00940). The DEGs of CK-D20, which involved the top three enriched pathways, were photosynthesis-antenna proteins (ko00196), DNA replication (ko03030), and benzoxazinoid biosynthesis (ko00402). The DEGs of CK-D30 involved the top three enriched pathways that were fatty acid degradation (ko00071), porphyrin and chlorophyll metabolism (ko00860), and starch and sucrose metabolism (ko00500). During this cold stress treatment, the reason for the increase in the starch and sucrose content was to better prevent freezing damage.

DEGs Co-Expression Clusters
To identify DEGs with similar expression patterns in response to cold stress, we performed K-means clustering tests for 12,926 DEGs (Table S7) and generated 8 optimal clusters (Figure 4). Compared with the control CK, clusters K1, K5, and K3 contained 4043 genes, and their expression levels gradually increased with the accumulation of cold stress time. The expression of clusters K2 (687 genes) and K4 (2033 genes) with the accumulation of cold stress time showed a trend of first decreasing, then being gradual, and then decreasing. Cluster K6 (3825 genes) showed a gradual decrease in gene expression with the accumulation of cold stress time. The gene expression of cluster K7 (1166 genes) accumulated with the cold stress time showed a rapid increase at first, then gradually decreased in the next 20 days. Cluster 8 had 1172 genes, which decreased rapidly in the first 10 days of responding to cold stress and stabilized in the next 20 days.

DEGs Co-Expression Network Analysis with WGCNA
The WGCNA was conducted on 12,927 differential genes (Table S8). The cluster tree diagram showed seven different modules (marked with different colors), with each branch representing a module, and each leaf in the branch denoting a gene. The characteristic gene of a module is the first main component of a given module and can be regarded as a representative of the gene expression profile of the module. These seven modules were well correlated with the tissue expression profiles of their representative genes in the four samples (Figure 5b,c). The correlation coefficient between the blue module and the D30 sample was 0.99 (the error was 1 × 10 −9 ), which indicated a high degree of correlation between the module and the sample. In addition, the correlation coefficient between the red module and the D10 sample was 0.95 (the error was 2 × 10 −6 ), while the correlation coefficient between the green module and the D20 sample was 0.94 (the error was 5 × 10 −6 ). To investigate the metabolic pathways in which the clustered genes may participate, several pathways were identified using the KEGG database ( Figures S2 and 4b). Among the eight clusters, clusters K1 and K3 exhibited similar expression patterns and were found to be involved in glutathione metabolism (ko00480), aminoacyl-tRNA biosynthesis (ko00970), and tyrosine metabolism (ko00350). Meanwhile, clusters K2, K4, and K6 had the largest number of genes and similar expression patterns, and they were found to participate in starch and sucrose metabolism (ko00500) and phenylpropanoid biosynthesis (ko00940). Interestingly, the metabolism of starch and sucrose (ko00500) and phenylpropanoid biosynthesis (ko00940) were involved in six out of the eight clusters. Glutathione metabonomics (ko00480) in the eight clusters was involved in three clusters.

DEGs Co-Expression Network Analysis with WGCNA
The WGCNA was conducted on 12,927 differential genes (Table S8). The cluster tree diagram showed seven different modules (marked with different colors), with each branch representing a module, and each leaf in the branch denoting a gene. The characteristic gene of a module is the first main component of a given module and can be regarded as  (Figure 5b,c). The correlation coefficient between the blue module and the D30 sample was 0.99 (the error was 1 × 10 −9 ), which indicated a high degree of correlation between the module and the sample. In addition, the correlation coefficient between the red module and the D10 sample was 0.95 (the error was 2 × 10 −6 ), while the correlation coefficient between the green module and the D20 sample was 0.94 (the error was 5 × 10 −6 ).
Genes 2023, 14, 844 9 of 15 a representative of the gene expression profile of the module. These seven modules were well correlated with the tissue expression profiles of their representative genes in the four samples (Figure 5b,c). The correlation coefficient between the blue module and the D30 sample was 0.99 (the error was 1 × 10 −9 ), which indicated a high degree of correlation between the module and the sample. In addition, the correlation coefficient between the red module and the D10 sample was 0.95 (the error was 2 × 10 −6 ), while the correlation coefficient between the green module and the D20 sample was 0.94 (the error was 5 × 10 −6 ).  We enriched different modules in the KEGG database. The main pathways for genes in the turquoise module were DNA replication (ko03030), benzoxazinoid biosynthesis (ko00402), cyanoamino acid metabolism (ko00460), and phenylpropanoid biosynthesis (ko00940). The main ways for genes in the red module to be enriched were photosynthesisantenna proteins (ko00196), circadian rhythm-plant (ko04712), and flavonoid biosynthesis (ko00941). The main pathways for genes in the green module were photosynthesis-antenna proteins (ko00196), plant hormone signal transduction (ko04075), and monoterpenoid biosynthesis (ko00902). The main pathways for genes in the blue module were glutathione metabolism (ko00480), cutin, suberin, and wax biosynthesis (ko00073), and aminoacyl-tRNA biosynthesis (ko00970).

Screening the Key Genes of Cold Treatment Changes with WGCNA
Complex traits in crops are typically controlled by multiple transcriptional networks. To identify the co-expression networks related to the cold resistance, we utilized the R WGCNA software, which was based on the FPKM and phenotypic data from the four treatments: before cold stress treatment, and after 10, 20, and 30 days of cold stress treatment, respectively. The clustering of the samples and correlation coefficients showed that the biological replicates were highly repeatable, and no outliers needed to be eliminated ( Figure S3a). We used the automatic R package blockwiseModules network construction We enriched different modules in the KEGG database. The main pathways for genes in the turquoise module were DNA replication (ko03030), benzoxazinoid biosynthesis (ko00402), cyanoamino acid metabolism (ko00460), and phenylpropanoid biosynthesis (ko00940). The main ways for genes in the red module to be enriched were photosynthesisantenna proteins (ko00196), circadian rhythm-plant (ko04712), and flavonoid biosynthesis (ko00941). The main pathways for genes in the green module were photosynthesis-antenna proteins (ko00196), plant hormone signal transduction (ko04075), and monoterpenoid biosynthesis (ko00902). The main pathways for genes in the blue module were glutathione metabolism (ko00480), cutin, suberin, and wax biosynthesis (ko00073), and aminoacyl-tRNA biosynthesis (ko00970).

Screening the Key Genes of Cold Treatment Changes with WGCNA
Complex traits in crops are typically controlled by multiple transcriptional networks. To identify the co-expression networks related to the cold resistance, we utilized the R WGCNA software, which was based on the FPKM and phenotypic data from the four treatments: before cold stress treatment, and after 10, 20, and 30 days of cold stress treatment, respectively. The clustering of the samples and correlation coefficients showed that the biological replicates were highly repeatable, and no outliers needed to be eliminated ( Figure S3a). We used the automatic R package blockwiseModules network construction method to identify co-expression modules ( Figure S3c). This allowed the module to be visualized with a color scheme that shows highly correlated genes in the same color and weakly correlated genes in a different color ( Figure S3b). The module construction process showed that the functional color modules were divided. After merging modules with similar expression patterns, it produced seven color modules, each of which consisted of genes with similar expression patterns (Figure 5b).
To identify hub genes in our interested modules, we selected the top 150 genes according to the kME value as hub genes in each module. In the 6618 genes of the turquoise module, based on the iTAK results, a total of 333 TFs, 64 TRs, and 257 kinase genes were identified. TraesCS7B01G016800 (UDP-Glc glucosyltransferase) was selected as the hub gene of the turquoise module. In the primary TraesCS7B01G016800 network, TraesCS7B01G016800 was directly co-expressed with 149 genes, including 2 TFs (belonging to the bHLH and HB-HD-ZIP families, respectively), 2 TRs, and 1 kinase genes (Figure 6a). To identify hub genes in our interested modules, we selected the top 150 genes according to the kME value as hub genes in each module. In the 6618 genes of the turquoise module, based on the iTAK results, a total of 333 TFs, 64 TRs, and 257 kinase genes were identified. TraesCS7B01G016800 (UDP-Glc glucosyltransferase) was selected as the hub gene of the turquoise module. In the primary TraesCS7B01G016800 network, TraesCS7B01G016800 was directly co-expressed with 149 genes, including 2 TFs (belonging to the bHLH and HB-HD-ZIP families, respectively), 2 TRs, and 1 kinase genes ( Figure  6a). In the 233 genes of the red module, a total of 11 TFs, 4 TRs, and 2 kinase genes were uncovered. TraesCS6B01G383200 (cold shock protein CS66 in wheat, Gene Name = CS66) was selected as the hub gene of the red module. In the primary TraesCS6B01G383200 network, TraesCS6B01G383200 was co-expressed with 149 genes, including 7 TFs (bHLH and HSF families), 1 TR, and 2 kinase genes (Figure 6b).
In the 253 genes of the green module, 13 TFs, 1 TR, and 12 kinase genes were identified. TraesCS6B01G274200 (ERF034 in Arabidopsis thaliana) was selected as the hub gene of the green module. In the primary TraesCS6B01G274200 network, TraesCS6B01G274200 was directly co-expressed with 149 genes, including 8 TFs (bHLH, HSF, WRKY, and APETALA2/ETHYLENE RESPONSE FACTOR (AP2/ERF) families) and 8 kinase genes (Figure 6c).
The correlation coefficient of the blue module was the highest with the samples treated with cold stress for 30 days. Most genes in this module were enriched in glutathione metabolism (ko00480). In the 4076 genes of the blue module, a total of 215 TFs, 32 TRs, and 290 kinase genes were identified. TraesCS1D01G189600, TraesCS1D01G190100, TraesCS5B01G426300, TraesCS7D01G030800, and TraesCS7D01G050800 were selected as In the 233 genes of the red module, a total of 11 TFs, 4 TRs, and 2 kinase genes were uncovered. TraesCS6B01G383200 (cold shock protein CS66 in wheat, Gene Name = CS66) was selected as the hub gene of the red module. In the primary TraesCS6B01G383200 network, TraesCS6B01G383200 was co-expressed with 149 genes, including 7 TFs (bHLH and HSF families), 1 TR, and 2 kinase genes (Figure 6b).
In the 253 genes of the green module, 13 TFs, 1 TR, and 12 kinase genes were identified. TraesCS6B01G274200 (ERF034 in Arabidopsis thaliana) was selected as the hub gene of the green module. In the primary TraesCS6B01G274200 network, TraesCS6B01G274200 was directly co-expressed with 149 genes, including 8 TFs (bHLH, HSF, WRKY, and APETALA2/ETHYLENE RESPONSE FACTOR (AP2/ERF) families) and 8 kinase genes (Figure 6c).
The correlation coefficient of the blue module was the highest with the samples treated with cold stress for 30 days. Most genes in this module were enriched in glutathione metabolism (ko00480). In the 4076 genes of the blue module, a total of 215 TFs, 32 TRs, and TraesCS1D01G189600, TraesCS1D01G190100, TraesCS5B01G426300, TraesCS7D01G030800, and TraesCS7D01G050800 were selected as the hub genes of the red module. In the primary network, these hub genes belong to glutathione S-transferase. The co-expression network of the blue module comprised 150 genes, including 5 TFs (bHLH, MADS-MIKC, HB-HD-ZIP, MYB, and PLATZF families) and 1 kinase gene (Figure 6d).
The hub genes of each module were basically enriched in the key KEGG pathway of each module. The most interesting thing was that bHLH was involved in the expression of four important modules, which also showed that bHLH played important roles in the abiotic stress responses. There were also some key TFs involved in abiotic stress including WRKY and AP2/ERF-ERF TFs.

Validation of Expression Level of DEGs
To validate the results of our RNA-seq, eight DEGs were chosen for transcription validation using qRT-PCR. The relative expression levels of these eight key DEGs were analyzed using qRT-PCR across the four cold stress treatments of the seedlings to assess the accuracy of the RNA-seq results. The relative expression patterns of the tested DEGs were found to be positively correlated with the fold change variations obtained from the RNA-seq results (Figure 7). The correlation coefficients between qRT-PCR and RNA-seq for the 15 DEGs were observed to be within the range of 0.74 to 0.99. This indicated that our RNA-seq results were both reliable and accurate.
Genes 2023, 14, x FOR PEER REVIEW 12 of 17 the hub genes of the red module. In the primary network, these hub genes belong to glutathione S-transferase. The co-expression network of the blue module comprised 150 genes, including 5 TFs (bHLH, MADS-MIKC, HB-HD-ZIP, MYB, and PLATZF families) and 1 kinase gene (Figure 6d). The hub genes of each module were basically enriched in the key KEGG pathway of each module. The most interesting thing was that bHLH was involved in the expression of four important modules, which also showed that bHLH played important roles in the abiotic stress responses. There were also some key TFs involved in abiotic stress including WRKY and AP2/ERF-ERF TFs.

Validation of Expression Level of DEGs
To validate the results of our RNA-seq, eight DEGs were chosen for transcription validation using qRT-PCR. The relative expression levels of these eight key DEGs were analyzed using qRT-PCR across the four cold stress treatments of the seedlings to assess the accuracy of the RNA-seq results. The relative expression patterns of the tested DEGs were found to be positively correlated with the fold change variations obtained from the RNA-seq results (Figure 7). The correlation coefficients between qRT-PCR and RNA-seq for the 15 DEGs were observed to be within the range of 0.74 to 0.99. This indicated that our RNA-seq results were both reliable and accurate.

Discussion
The plant acclimatization response to different abiotic stresses involves complex molecular regulatory networks. For example, many isomers of plant thiol peroxidase play a central role in plant adaptation to stresses such as low temperatures, high amounts of light, or osmotic pressures [50]. In order to maintain plant life under stress conditions, antioxidant enzymes such as superoxide dismutase (SOD), catalase (CAT), peroxidase (POX), etc., play a crucial role [51]. This study reported an important economic gramineous wheat transcriptome of cold stress. High-quality transcription modules were achieved by sequencing with Super Illumina sequencing. Our present results can serve to not only enrich the existing transcriptome data resources for wheat but also offer valuable insights into the molecular mechanism of wheat cold stress resistance.

Discussion
The plant acclimatization response to different abiotic stresses involves complex molecular regulatory networks. For example, many isomers of plant thiol peroxidase play a central role in plant adaptation to stresses such as low temperatures, high amounts of light, or osmotic pressures [50]. In order to maintain plant life under stress conditions, antioxidant enzymes such as superoxide dismutase (SOD), catalase (CAT), peroxidase (POX), etc., play a crucial role [51]. This study reported an important economic gramineous wheat transcriptome of cold stress. High-quality transcription modules were achieved by sequencing with Super Illumina sequencing. Our present results can serve to not only enrich the existing transcriptome data resources for wheat but also offer valuable insights into the molecular mechanism of wheat cold stress resistance.
WGCNA showed that four expression modules had the highest correlation with cold stress. The combined analysis of WGCNA and K-means showed that five metabolic pathways were significantly associated with traits of interest, including phenylpropanoid biosynthesis (ko00940), benzoxazinoid biosynthesis (ko00402), flavonoid biosynthesis (ko00941), hormone signal transduction (ko04075), and glutathione metabolism (ko00480). This indicates that the transcription regulatory factors of related genes may be potential targets to improve the cold resistance of wheat crops.
The correlation coefficient of the blue module was the highest treated with cold stress for 30 days. The genes presented in the blue module exhibited a significant enrichment in glutathione metabolism (ko00480). Glutathione is one of the major endogenous antioxidants in plants, and plays a crucial role in plant defense mechanisms [52]. Glutathione S-transferase (GST), glutathione peroxidase (GPX), and glutathione reductase (GR) utilize glutathione to safeguard plants against abiotic stress [53]. Cold stress is known to produce reactive oxygen species (ROS) in plants and ROS accumulation may play a role in regulating cell death [54]. GST and GPX, as key enzymes of cellular detoxification systems, are essential in defending cells against ROS [55].
Another module that we were interested in was the green module. The genes present in the green module exhibited a significant enrichment in plant hormone signal transduction (ko04075) including the AP2/ERF TF family TF [3]. The ERF TF family plays a crucial role in regulating plant development, hormonal signaling, and environmental responses [56,57]. The AP2/ERF family is a large family of plant-specific TFs, which share a conserved DNA-binding domain [58]. It is worth noting that the regulation of cuticular wax accumulation is attributed to certain members of the AP2/ERF family that were originally isolated from Arabidopsis, namely WAX INDUC-ER1/SHINE1 (WIN1/SHN1) [59,60]. Three SHINE genes AtSHN1, AtSHN2, and AtSHN3 were classified into the ERF-B6 clade [57,61]. The overexpression of SHINE1 (SHN1) in Arabidopsis thaliana increased the deposition of epidermal wax, increased the permeability of stratum corneum, decreased the stomatal index, and enhanced drought resistance [60]. Studies have indicated that TaSHN1 may act as a positive regulator of drought stress tolerance in wheat, potentially by promoting an increase in alkane accumulation and a decrease in stomatal density. Therefore, TaSHN1 holds promise as a valuable candidate gene for developing new wheat cultivars with improved drought tolerance [62]. The gene WIN1 is a positive regulator of Arabidopsis epidermal wax biosynthesis. TraesCS6A01G181400 and its homologous sequence TraesCS6B01G210300 are homologous to WIN1 in Arabidopsis thaliana and are differentially expressed in this module. Thus, we inferred that the regulatory network of the WIN1 involved in the complex may have been evolutionarily conserved between Arabidopsis and wheat. Considering the differences between Arabidopsis and wheat, the conserved WIN1 gene might be the best operable target locus in the breeding of stress-resistant varieties in wheat and related crops.
The cold shock gene was up-regulated during cold stress in wheat, which enhanced the cold adaptability of wheat. A cold shock affected the gene expression along with the starch and sucrose metabolic pathways. During cold adaptation, a metabolic readjustment occurs, which slows down growth and development and activates the defense process [63].

Conclusions
In this study, we identified 12,926 DEGs from transcriptome comparisons of seedling Jing 841 accessions with cold stress at four developmental stages. By studying the gene expression patterns and comparing the whole transcriptome in different cold stress periods and seedlings without cold stress treatment, it was determined that prolonged cold stress forced the suppression of plant metabolism, especially photosynthesis and sugar metabolism, and induced stress response genes, such as ERF and bHLH TFs. However, the regulation of CBF gene expression was not the only cold signaling factor. At the same time, it was also affected by various plant hormones, which indicated that CBFs were the center of the interaction between chilling and the hormone signaling pathways, and they played a variety of roles in regulating plant growth and cold tolerance. The aim of this study was to elucidate the relationship between the developmental stages of winter and spring habit varieties in the wheat seedling period and frost resistance.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes14040844/s1, Figure S1: Functional classification of the differentially expressed genes (DEGs) referring to the Gene Ontology (GO) database; Figure S2: Differentially expressed genes (DEGs) involved top 20 enriched pathways of the eight K-means clusters referring to KEGG database; Figure S3: Screening the key genes of cold treatment changes with WGCNA; Table S1: Primers used in qRT-PCR; Table S2: Statistic data of the cDNA sequencing; Table S3: Statistics of sequence alignment; Table S4: Summary of the new unigene annotation; Table  S5: Summary of the unigene annotation; Table S6: Analysis of DEGs between samples treated with cold stress at different times; Table S7: Identification of DEGs with similar expression patterns during cold stress development; Table S8: The DEGs classified by K-means clustering (K = 8).