Transcriptomic profiling of wheat stem during meiosis in response to freezing stress

Low temperature injury in spring has seriously destabilized the production and grain quality of common wheat. However, the molecular mechanisms underlying spring frost tolerance remain elusive. In this study, we investigated the response of a frost-tolerant wheat variety Zhongmai8444 to freezing stress at the meiotic stage. Transcriptome profiles over a time course were subsequently generated by high-throughput sequencing. Our results revealed that the prolonged freezing temperature led to the significant reductions in plant height and seed setting rate. Cell wall thickening in the vascular tissue was also observed in the stems. RNA-seq analyses demonstrated the identification of 1010 up-regulated and 230 down-regulated genes shared by all time points of freezing treatment. Enrichment analysis revealed that gene activity related to hormone signal transduction and cell wall biosynthesis was significantly modulated under freezing. In addition, among the identified differentially expressed genes, 111 transcription factors belonging to multiple gene families exhibited dynamic expression pattern. This study provided valuable gene resources beneficial for the breeding of wheat varieties with improved spring frost tolerance.


Introduction
Wheat (Triticum aestivum L.) is one of the most widely grown staple crops worldwide. The spring frost injury has become the limiting factor of wheat production because it has dramatically influenced the yield, quality, and the geographical distribution of wheat (Chen et al., 2020). The Yellow and Huai River Valley Winter Wheat Zone (YHVWWZ) is the major producing area of China's winter wheat. Frost damage in spring generally occurs from February to April, when the plants are at the jointing to the booting stage of wheat growth. At this time, the wheat tissues are tender and contain a lot of water, so they are easy to get damaged (Frederiks et al., 2015). Exposure to low temperature may cause the damage of leaf, stem and young spike, which could result in the decrease in the grain number and quality. From 2001 to 2021, the frequency of spring frost injury in the YHVWWZ was as high as 40% (Liu et al., 2021). In 2018, it has caused the yield losses of 3,000 kg per hectare in wheat varieties with weak freezing tolerance when compared to freezing-resistant varieties (Ou and Wang, 2019). Therefore, breeding and promoting freezing-tolerant wheat varieties is an important approach to deal with the disaster of spring frost injury in YHVWWZ. However, at present, knowledge of the physiological mechanisms underlying the freezing tolerance is still elusive, limiting the genetic improvement of wheat varieties.
Low temperature stress includes both cold/chilling stress (>0°C) and freezing stress (<0°C). When exposed to low temperature, plants can perceive the signal and induce a series of complex changes in their morphological structure, internal substances and photosynthesis (Ding et al., 2019). The accumulation of soluble sugar and amino acid was observed in many plant species (Nakabayashi and Saito, 2015). A series of mechanisms are initiated to protect and maintain the normal metabolism and growth of plants, so as to avoid or minimize the damage of low temperature. During this complex biological process, a large number of genes are activated and play an important role for plants to survive. For instance, transcription factors (TFs) are thought to control many cold responsive genes through direct binding to cis-acting elements in the promoter regions. Regulation of cold-regulated (COR) genes by CRT/ DRE-binding factors (CBFs) constitutes the predominant cold signaling pathway and was found to be conserved in both dicots and monocots (Park et al., 2015). In diploid Triticum monococcum, three CBF genes (CBF12, CBF14, CBF15) have been identified to be highly associated with frost tolerance (Knox et al., 2008). Overexpression of CBF14 and CBF15 from winter wheat resulted in the enhanced frost tolerance in spring barley (Solteśz et al., 2013). The CBF transcription factors belong to the APETALA2/Ethylene response element binding protein (AP2/ EREBP) family. Other transcription factors including MYB, NAC, bZIP, and WRKY TFs have also been identified to be linked to the tolerance of low temperature stress (Su et al., 2010;Mao et al., 2012).
In addition, phytohormone always acts as a crucial signal to modulate multiple plant processes and has been reported to play an essential role for plant to cope with low temperature stress. A key hormone is abscisic acid (ABA), which not only controls the internal water deficiency by regulating stomatal conductance but also regulates the transcriptome of downstream stress-related genes. It has been shown that CBF transcription factors induced the biosynthesis of ABA, which subsequently stimulated the expression of COR genes (Sharabi-Schwager et al., 2010;Wang et al., 2016). The exogenous application of ABA significantly improved the low-temperature tolerance of tomato (Kim et al., 2002), pepper (Guo et al., 2012), magnolia liliiflora  and wheat (Yu et al., 2020). Moreover, accumulating evidence indicated that jasmonic acid (JA) could mediate the responses to low temperature stress with other hormones, such as ABA, auxin, ethylene, and gibberellin (Hu et al., 2017). Wang et al. (2016) reported that JA could activate the ABA-dependent CBF signaling pathway and positively modulate the reactions to cold stress in tomato. Blocking of JA biosynthesis and signaling caused the hypersensitivity to freezing stress (Hu et al., 2013). In Arabidopsis, it has been reported that ethylene and gibberellin (GA) negatively regulated plant cold tolerance (Shi et al., 2012;Richter et al., 2013). For instance, a GA-deficient mutant showed increased expression level of cold-responsive genes and enhanced resistance to cold stress (Richter et al., 2013). These studies indicate the complex nature of plant's responses to low temperature.
It has been suggested that properties of cell wall strongly influence the frost resistance of plants. When subjected to freezing, ice crystallization takes place in the intercellular space which could lead to the dehydration and mechanical damage of cells. Plasma membrane and cell wall work together to perceive the external signal of freezing stress and in turn cause the corresponding changes in cell wall content, composition and structure. This phenomenon has been observed in multiple plant species such as Arabidopsis, winter wheat, maize and tobacco (Sugiyama and Shimazaki, 2007;Bilska-Kos et al., 2016;Willick et al., 2017;Parrotta et al., 2019;Takahashi et al., 2019). There is increasing evidence showing that cell wall genes were involved in the regulation of plant low temperature tolerance. A transcriptomic study performed by Le et al. (2015) showed that several cell wall genes were dramatically up regulated in Arabidopsis accessing to sub-zero temperatures. Mutant lacking XTH19 gene, which encode a xyloglucan endotransglucosylase/ hydrolase, showed hypersensitivity to low-temperature stress compared to wild-type plants (Takahashi et al., 2021). In addition, a novel frost tolerant gene, SENSITIVE-TO-FREEZING8 (SFR8), has been cloned and demonstrated to influence pectin fucosylation (Panter et al., 2019).
Large-scale analysis of gene expression pattern by transcriptome techniques provides a cost-effective strategy for elucidating the molecular mechanisms underlying physiological processes and it can substantially increase the efficiency of identifying genes of interest. A range of plant species such as Arabidopsis, rice, wheat, maize, and cotton have been extensively subjected to transcriptional programming in response to different stresses (Jian et al., 2014;Huang et al., 2019;Sga et al., 2019;Wang et al., 2019). Zhao et al. (2019) integrated transcriptome and metabolome analyses for understanding how wheat plants respond to low temperature. 29,066 DEGs were identified, and pathways of phytohormone signaling and proline biosynthesis were shown to be significantly modulated under freezing treatment.
The objective of this study is to investigate the phenotypic changes as well as transcriptome dynamics of wheat cultivar zhongmai8444 in response to freezing temperature at the meiosis stage. Our results will provide valuable insights into understanding molecular mechanisms underlying the frost tolerance and help identify new gene resource that is beneficial for genetic improvement of spring frost tolerance in wheat.

Materials and methods
Plant materials, freezing treatment, and seed setting measurement Zhongmai8444 is a spring wheat variety breed in our lab. It showed great freezing tolerance since it could survive in the winter of Beijing (-5°C to 2°C) by field trials for five consecutive years. Zhongmai8444 was grown in a growth chamber maintaining at 22°C during the day, 15°C during the night, 16-h light, 8-h dark, 80% humidity and a PPFD of 450 µmol m -2 s -1 . Three plants were grown in each pot and their growth stages were recorded according to the Zadoks scale (Zadoks et al., 1974). Developmental stages of pollen mother cells from immature anthers were examined following Draeger and Moore (2017). When the distance from flag leaf to the adjacent leaf was around 2-3 cm (Supplemental Figure S1), Zhongmai8444 was verified to be at meiotic stage and subjected to freezing treatment shown in Supplemental Figure S2. Wheat plants grown at 22°C (normal temperature growth) prior to freezing were set as the control group (time point T1). Temperature in growth chamber was gradually decreased from 22°C to -2.5°C using 6 h (time point T2) and then maintained at -2.5°C for 48 h (time point T4). Plants were also harvested after frozen at −2.5°C for 24 h (time point T3). Finally, temperature gradually increased from -2.5°C to 22°C (time point T5) and plants were collected as recovery group. Samples were named as J1 to J5 according to the time points (T1 to T5). Stem tissues from six individual plants were pooled (with three replications) and immediately frozen in liquid nitrogen for subsequent RNA extraction. In addition, another set of plants were treated with freezing (-2.5°C for 48 h) during meiosis as described above and transferred to normal growth condition for continuing growth. Plants that were maintained at normal growth condition were used as the control. The seed setting rate was analyzed using the equation: seed setting rate=grain number each spike/flower number each spike × 100%.

Cell wall thickness measurement
To investigate the cell wall morphology of stem tissue after freezing treatment, xylem cells from three individual plants were selected for further measurement. Thin cross-sections were cut from the base of wheat stem using a sterile razor blade and were subsequently incubated in 0.02% toluidine blue solution for 2 min at room temperature. After thoroughly washing, respective stem cross-sections were mounted in water prior to imaging on an OLYMPIC CX41RF microscope equipped with an EOS 450D camera (Canon). Cell wall morphology and thickness of the xylem cells were subsequently analyzed using the Image J software (National Institutes of Health; http://rsb. info.nih.gov/ij/). The cell wall thickness was analyzed from at least 20 xylem cells collected from three individual plants. For the observation by a transmission electron microscope (TEM), samples were prepared following Vigani et al. (2018) with some modifications. Briefly, stem tissues were fixed overnight in 3% (v/v) glutaraldehyde at room temperature and subsequently post fixed with 1% (w/v) osmium tetroxide for 2 h at 4°C. After dehydration in a graded ethanol series, samples were embedded in resin for ultrathin sectioning and visualization with a H-7500 TEM (HITACH, Japan) at 80 kV.
Gene expression levels were quantified as Fragments Per Kilobase of transcript sequence per million mapped fragments (FPKM) using the featureCounts v1.5.0-p3 program (Yang et al., 2014). Principal component analysis (PCA) was performed to evaluate the reproducibility of different samples under the same treatment.

Differential expression analysis
Differential gene expression analysis was performed using the DESeq2 (Love et al., 2014) package in R (1.20.0). The Benjamini and Hochberg method was used to adjust P-values to control the false discovery rate (FDR). Differentially expressed genes were identified according to the threshold of |log2FoldChange| > 1 and FDR < 0.05. Gene Ontology (GO) and Kyoto Encyclopedia of Gene and Genomes (KEGG) pathways enrichment analysis of the DEGs were performed using clusterProfiler R package (Yu et al., 2012).

Real-time quantitative PCR
Real-time quantitative PCR was performed on the stem tissues to examine expression patterns of selected genes. Samples were run in triplicate with SuperReal SYBR Green PreMix Plus (Tiangen) on a Bio-Rad CFX96 Touch Real-Time PCR Detection System. The PCR reactions were performed under the following conditions: 95°C for 15 min, followed by 40 cycles of 95°C for 10 s, 58°C for 20 s, and 72°C for 30 s. Relative quantities were calculated and normalized to the actin gene, which was used as an internal control, using the 2 -DDCt algorithm. The data are shown as means ± standard errors (SE) and primers used for qRT-PCR analysis are listed in Supplementary Table S4.

Phenotypic analysis of Zhongmai8444 after freezing treatment
In previous study, we have demonstrated that Zhongmai8444 was extremely sensitive to freezing stress at the meiotic stage/phase by evaluating the damage of freezing on plant stems, leaves, spikes and other organs during the reproductive development . Consistently, the injury of stems (Supplementary Figure S3) and cell wall thickening in the vascular tissue were observed after Zhongmai8444 was treated with -2.5°C for 48 h (Figures 1A-C). Cell walls were found to be severely damaged when plants were recovered to normal room temperature (Figures 1D,E). In addition, plant height, grain number per spike and seed setting rate were quantified after the plants were treated with freezing stress or grown in normal growth conditions. Results indicated that Zhongmai8444 exhibited significantly reduced height, grain number and seed setting rate after -2.5°C freezing treatment for 48 h at the meiotic stage (Figures 1F, G; Table 1). To further understand the molecular mechanism underlying wheat frost tolerance, we performed RNA sequencing of the stem to identify differentially expressed genes responsive to freezing stress.

Mapping of RNA-Seq data to wheat genome
Utilizing the paired-end Illumina sequencing technology, the extracted RNA from stem tissue was sequenced to compare gene expression at different time points under freezing stress, samples were named as J1 to J5 accordingly. After quality control, sequencing data yield approximately 70 million highquality clean reads per library shown in Supplemental Table S1. The clean reads shared a GC content of approximately 55% and the alignment rate uniquely mapped to the Chinese Spring iwgsc_refseqv1.0 reference genome (https://urgi.versailles.inra. fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/) was ranged from 88.96% to 90.76%. In addition, Pearson correlation analysis (PCA) and principal component analysis were performed to determine the quality of biological replications of each sample (Supplemental Figure S4, Figure  S5). Results showed that the sample J3-2 was relatively dispersed from the other two biological replicates and removed for further analysis. All the other samples had good repeatability with ePearson coefficient higher than 0.87 (Supplemental Figure S4).

Identification of differentially expressed genes under freezing
To systematically investigate the dynamics of transcriptome profiles of wheat stem, we identified genes that were differentially expressed by at least twofold (|log2FoldChange| > 1 and padj < 0.05) with -2.5°C freezing treatment. A clear difference was observed in the number of up-regulated and repressed DEGs at the indicated time points compared to control plants ( Figure 2A). A total of 2016 up-regulated and 994 down-regulated genes were detected in stem tissue after 6-h-cooling (J2 vs J1). An increase in the number of significantly up-regulated (5058 to 7000) and down-regulated (6985 to 11778) DEGs was observed with the extended duration offreezing stress (from T2 to T4), which suggested the ongoing adaptation of plants to low temperature conditions. However, there was a small change in the number of DEGs when plants were recovered at 22°C after freezing stress (J5 vs J1). We subsequently constructed the Venn analysis to investigate the number of up-regulated and downregulated genes that were commonly shared at different time points (Figures 2B,C respectively). Results showed that respectively 1010 up-regulated and 230 down-regulated genes were detected common at all time points.

Gene ontology clustering of DEGs
To functionally characterize the identified DEGs, a GO analysis was performed to classify them into the three categories: molecular function (MF), cellular component (CC) and biological process (BP). Significantly enriched GO terms identified at different time points with freezing treatment (J2 vs J1, J3 vs J1 and J4 vs J1 comparisons) were determined. The top 4 GO terms from each category were respectively presented in Figure 3 and highly overlap among samples collected from different time points. For instance, in the biological process category, the GO term "response to auxin", "response to endogenous stimulus", "response to hormone" were significantly enriched in all samples, suggesting an involvement of genes participated in the signal perception and transduction pathway. Similarly, under the cellular component and molecular function category, the most enriched GO terms "extracellular region", "apoplast", "cell wall" and "protein heterodimerization activity" were also detected in all samples. However, the "cell periphery" and "calcium ion binding" were only shared in the samples treated with 24 h or 48 h freezing stress. In addition, in the 48 h freezing treated samples, the other term "xyloglucan:xyloglucosyl transferase activity" and "glucosyltransferase activity" were also related to the cell wall metabolism and consistent with the detected phenotype of cell wall thickening. Therefore, we hypothesized that hormone signal transduction and cell wall biosynthesis process may be modulated to promote plant resistance to freezing.

Pathways enrichment analysis of DEGs
The KEGG analysis of DEGs in stem was also performed to identify pathways displaying significant changes (padj ≤ 0.05) in response to freezing (Table 2). "Plant hormone signal  transduction" pathway was significantly enriched at every time point. In addition to this pathway, "circadian rhythm-plant" was also enriched after 6-h-cooling (J2 vs J1). In 24 h freezing treated samples (J3 vs J1), four pathways significantly changed including "plant-pathogen interaction", "plant hormone signal transduction", "MAPK signaling pathway-plant", and "phenylpropanoid biosynthesis". In addition, analysis revealed that DEGs in 48 h freezing treated samples (J4 vs J1) were associated with "plant-pathogen interaction", "plant hormone signal transduction", "phenylpropanoid biosynthesis", "circadian rhythm-plant", "photosynthesis-antenna proteins", "MAPK signaling pathway-plant", and "starch and sucrose metabolism" pathways. After samples were recovered after freezing treatment (J5 vs J1), the "photosynthesis-antenna proteins", "photosynthesis", "plant hormone signal transduction", and "MAPK signaling pathway-plant" pathways were significantly enriched among DEGs. Consistent with the GO analysis, our results emphasized the involvement of plant hormone in the response of wheat to freezing.

Representative mostly up-regulated and down-regulated DEGs
Venn analysis identified 1010 up-regulated and 230 downregulated genes shared by all time points. We listed the expression of top 20 (ranked by the padj value of J4 vs J1 comparison) up-regulated (Table 3) and down-regulated DEGs (Table 4). We observed several low temperature related genes that have been significantly up regulated, including the abscisic acid (ABA) responsive protein (TraesCS5B01G502200, TraesCS5A01G488500), dehydrin (TraesCS4B01G228000, TraesCS6B01G383200) and late embryogenesis abundant (LEA) proteins (TraesCS1A01G423800, TraesCS3D01G158600, TraesCS3A01G175400). In addition, expression levels of genes involved in diverse biosynthesis process of cell wall, carbohydrates, amino acids, cytoskeleton, and other metabolites were also significantly affected. Few genes encoding the transcription factor were also identified, including those containing GRAM domain, NAC domain, and Myb-like DNA-binding domain. These results indicated that widespread genes showed substantial changes in expression level following freezing temperatures.

Identification of genes involved in hormone and signal transduction
The GO and KEGG enrichment analyses of the DEGs both highlighted that the plant hormone signal transduction pathway (dosa04075) was particularly affected by the freezing treatment ( Figure 3; Table 2). Among the 1240 DEGs shared by all time points, many genes encoding a protein kinase and probably involved in the signal transduction were detected (Supplemental Table S2). More than 40 genes predicted to encode a protein kinase were differentially expressed under freezing treatment. Most of the genes were up regulated while only 6 genes were down regulated by Histogram showing the most significantly enriched GO term from samples collected at each time point during freezing treatment. Color bars represent different functional classifications, including biological process (orange bars), cellular components (grey bars) and molecular function (blue bars). freezing stress. Genes encoding a receptor-like protein kinase constituted the largest group and may play key roles in mediating stress-associated signaling pathway, such as the leucine-rich repeat receptor-like protein kinase (Yang et al., 2022). KEGG analysis of the DEGs suggested that genes encoding protein phosphatase 2C (TraesCS3A01G362200 and TraesCS5A01G183600), which was associated in the carotenoid and ABA biosynthesis, were induced at all time points under freezing condition. The expression of  TraesCS3A01G362200 was rapidly up-regulated after 6-h cooling and progressively up-regulated by 5.2 and 10 folds after maintaining at -2.5°C for 24 h and 48 h, respectively, compared to that in control group. The other gene TraesCS5A01G183600 showed similar expression pattern and 8 folds up-regulated expression at time point T4. In addition, 66 DEGs associated with the biosynthesis or signal transduction of the phytohormones like ethylene, auxin, ABA or GA were listed in Supplemental Table S2. All the ethylene related genes were induced upon freezing treatment and determined to encode AP2/ERF domain-containing TFs. ERF genes are critical for multiple plant species to deal with low temperature stress. In Arabidopsis, ERF105 has been suggested to operate in conjunction with the well-known CBF regulon (Bolt et al., 2017).

Identification of transcription factors in response to freezing stress
Numerous families of TFs are known to play an important role when plants are challenged with various abiotic stresses. In current study, 111 putative TFs were identified to be differentially expressed (|log2FoldChange| > 1 and padj < 0.05) at all time points with freezing treatment (Table 5; Supplemental Table S3). Genes encode AP2/ERF domain-containing proteins constituted the largest group, followed by zinc finger proteins, NAC domain-containing proteins, MYB TFs, heat shock TFs, WRKY TFs, GRAS TFs, Growthregulating factors, Homeobox domain TFs, bZIP proteins and ethylene insensitive 3 TF (Table 5). The identified TFs were classified into five subgroups based on their expression patterns (Figure 4). The expression level of TFs in Subgroup I declined immediately after cooling and became more dramatically reduced after continuous freezing treatment. In contrast, the transcript abundance of TFs in Subgroup II was significantly induced at the early time point and then decreased upon the duration of freezing treatment, which suggests that they play roles in early transduction of the low-temperature signal. Subgroup III TFs were active after freezing but showed a peak of expression when plants were recovered in normal growth condition. Subgroup IV and V, the two largest subgroups, were characterized by genes showing the highest expression after 24-h and 48-h freezing treatment, respectively. In addition, we listed the 10 most highly induced TF genes in stem tissue commonly shared in all time points (Table 6). Four genes encoding AP2 domain containing proteins and three NAC domain containing TF genes were identified. The transcript level of two NAC TF encoded genes (TraesCS5D01G481200 and TraesCS5B01G480900) was most significant enhanced subjected to 24-h-freezing which showed the log2 value of fold changes as 8.52 and 7.48, respectively (Table 6). These two genes were confirmed to be the homologue of TaNAC2-5A which has been reported to play a role in plant freezing tolerance (Mao et al., 2012). Further experiments could be carried out to verify the functional role of other candidate TFs listed in Table 6.

Identification of genes associated with cell wall biosynthesis and modification
When subjected to freezing, plants could form a thick cell wall to resist low temperature stress. The GO analysis in our study also indicated that the expression of cell wall related genes was significantly modulated by freezing ( Figure 3). Therefore, we have selected genes with a putative role in cell wall remodeling shared by J3 vs J1 and J4 vs J1 comparisons (Table 7). These genes were participated in the biosynthesis or modification of major cell wall components, such as cellulose, hemicellulose, and pectin. For instance, several genes encoding cellulose synthetic enzymes such as cellulose synthase and cellulose synthase-like proteins, endoglucanases, chitinase-like proteins and cobra-like proteins were identified differentially expressed (Table 7). In addition, genes encoding hemicellulose biosynthetic enzyme xyloglucan endo-transglycosylase constituted a large group and have been reported to play an important role in frost tolerance in Arabidopsis (Takahashi D, 2021). Our study indicated that the increased amounts of enzymes involved in cell wall biosynthesis and modification is in accordance with the cell wall thickening observed in wheat stem upon freezing.

Validation of gene expression by real time quantitative PCR analysis
We also performed qRT-PCR analysis to verify the accuracy and reproducibility of transcriptome profiles. According to the literature, six genes potentially involved in low temperature tolerance were selected for the analysis, such as zinc finger proteins (TraesCS6B01G246500 and TraesCS3A01G107200), protein kinase (TraesCS5B01G381500), actin depolymerizing factor (TraesCS5A01G478500), beta-amylase (TraesCS1A01G215500) and ABA-responsive protein (TraesCS5D01G503400) (Peng et al., 2015;Byun et al., 2021;Zhang et al., 2022). Our data showed that the expression pattern of selected genes obtained by qRT-PCR was highly consistent with RNA-Seq data, with the correlation coefficients (R 2 ) > 0.8 ( Figure 5). These results support the reliability of the RNA-Seq analysis.

Discussion
Freezing injury during the reproductive phase of wheat frequently happens in spring, which has become the limiting factor of wheat production in China. To investigate the effect of freezing stress on wheat growth and reproduction, we performed the freezing treatment on Zhongmai8444 at the meiotic stage in a controlled growth condition. Transcriptomic analysis was subsequently carried out to elucidate the molecular mechanism underlying the frost tolerance of wheat plants.
Consistent with our previous study, Zhongmai8444 showed hypersensitivity to freezing stress at the meiotic stage . At this time, the developing pollen mother cells were at the meiosis phase and freezing temperature may cause the increase frequency of cellular abnormality of pollen cells. Therefore, wheat plants showed a decrease in seed setting after treated with freezing ( Figure 1F, Table 1). Consistently, it has been reported that Australian wheat cultivars exhibited significant reduction in fertility when subjected to prolonged chilling temperature during meiosis (Barton et al., 2014). The decrease of grain number per spike and 1000-grain weights was observed when winter wheats were exposed to low temperature stress at booting stage Zhang et al., 2021). Results showed that altered sucrose metabolism and inhibited starch biosynthesis in young ear may lead to the reduction of yield. In addition, low temperature stress during floral development is known to cause yield losses in different crop plants such as rice (Li et al., 2022), sorghum (Wood et al., 2006), maize (Tranel et al., 2009), bell pepper (Mercado et al., 1997) and garlic (Shemesh Mayer et al., 2013). Stem injury and reduced plant height were also obvious in freezing-treated plants (Supplemental Figure S3, Figure 1G, Table 1). Vascular tissues in stem play an important role in transporting water, minerals and product of photosynthesis throughout the plants. Injury of stem and vascular tissue would inhibit the transportation of nutrients Ethylene insensitive 3 transcription factor 1 required for plant growth and development. As a result, wheat plants exhibited reduced height and production. Cell wall provides the structural integrity as well as the flexibility to a plant cell (Lampugnani et al., 2018). Our findings provide evidence supporting the remodeling of wheat cell wall under freezing stress. Under low temperature conditions, cell wall could help the dehydrated cells maintain the morphology by enhancing its mechanical strength and reduce the rate of cell dehydration by adjusting its structure (Johnson et al., 2018). Previous study demonstrated that chilling stress led to an increase of total cell wall amount as well as structural and compositional changes in cell wall during both cold acclimation and sub-zero acclimation (Takahashi et al., 2019). We identified a range of cell wall associated genes to be dramatically up regulated upon freezing stress by transcriptome analysis (Table 7). These genes were implicated in the metabolic events related to major cell wall component biosynthesis and may be critical for plant to resist freezing stress. For instance, Takahashi et al. (2021) suggested a potential role of Arabidopsis xyloglucan endotransglucosylase/ hydrolase XTH19 in cell wall remodeling which influenced the freezing tolerance after low temperature acclimation.
To understand the molecular mechanism underlying freezing tolerance, we performed GO and KEGG enrichment analyses on identified DEGs. Our results showed that the plant hormone signal transduction was particularly affected by -2.5°C treatment, suggesting a critical role of phytohormone in conferring response of wheat to the freezing stress. Phytohormones are known to play a significant role in a wide range of adaptive responses to abiotic stresses (Peleg and Blumwald, 2011). When subjected to stresses, accumulation of A heat map Showing the expression pattern of differentially expressed TF genes under freezing stress at different time points. Columns and rows in the heat map respectively represent the five time points of freezing treatment and common differentially expressed TFs. Color bars represent scaled FPKM of each TF gene.
ABA could promote stomatal closure, enhance water balance, and induce antioxidant defense systems to alleviate oxidative injury (Lee and Luan, 2012). It has been reported that ABA regulated several ICE-CBF pathway-dependent genes, promoting plant's resistance to freezing (Knight et al., 2004). Furthermore, the involvement of ethylene signaling in lowtemperature-stress response has been implicated in several studies (Zhang and Huang, 2010;Shi et al., 2012;Cataláet al., 2014), though its effect is debated. Other hormone such as GA also plays a role in plant responses to low-temperature stress and multiple phytohormones may crosstalk in modulating the expression of cold-responsive genes (Richter et al., 2013).
As the key regulator of transcription, TFs always play a critical role in mediating abiotic stress responses. In current study, we demonstrated that most differentially expressed TF genes belong to AP2/ERF, Zinc finger and NAC family groups (Table 5 , Table 6). The CBF transcription factors belong to the AP2/ERF superfamily and regulate a spectrum of cold-regulated (COR) genes. CBFs mediated signaling has been proposed to constitute the predominant cold signaling pathway in many plant species though it has not been well characterized (Park et al., 2015). In wheat, CBF14 and CBF15 TFs were identified to be linked to the frost-tolerance locus Fr-A2 and  positively regulate low-temperature-stress responses (Knox et al., 2008;Solteśz et al., 2013). Zinc finger and NAC TFs play essential roles in response to various abiotic stresses and have been implicated in the resistance offreezing stress (Hao et al., 2011;Zhang et al., 2016). It has been reported that a C2H2-type zinc finger protein PhZFP1 could regulate cold stress tolerance in Petunia hybrida by modulating galactinol synthesis activity . Overexpression of wheat TaNAC2 gene in Arabidopsis resulted in the enhanced tolerance to multiple abiotic stress (Mao et al., 2012). Furthermore, based on the listed mostly up-regulated and downregulated DEGs (Table 3; Table 4), we proposed several other genes were also participated in freezing responsive process, such as genes encoding dehydrins, actin depolymerizing factors, LEA proteins and beta-amylase. Dehydrins are a distinct group of LEA proteins that exhibited high hydrophilicity and were proved to be an important factor to regulated plant freezing tolerance (Ndong et al., 2002). Overexpression of a dehydrin gene ShDHN resulted in the improved freezing tolerance in tomato . As the best characterized LEA protein, COR15A resided at the membrane surface during dehydration and stabilized cell membranes under freezing stress (Bremer et al., 2017a;Bremer et al., 2017b). In addition, it has been suggested that the actin depolymerizing factor protein may be required for the cytoskeletal arrangement during cold acclimation and important for plants to tolerate freezing (Byun et al., 2021). Our results provide the potential candidates associated with several metabolic pathways involved in freezing responses. Further research is needed to verify the functional roles of these predicted genes and employ them in future breeding program.
In this study, we evaluated the effect of freezing stress on wheat growth and development during meiosis phase. A model for wheat plants respond to freezing is proposed ( Figure 6). Once subjecting to sub-zero temperature, a series of stress-related genes become active to minimize the damage of freezing and maintain plant's normal growth, including those encoding TFs, hormone signaling elements, and enzymes associated with cell wall remodeling. The functional validation of these candidate genes will help us better understand the mechanism underlying freezing tolerance and hold great potential for molecular breeding in spring freezing-tolerant wheat varieties.

Data availability statement
The datasets presented in this study can be found in online repositories. All the raw data have been deposited in the NCBI Schematic diagram showing the freezing induced changes of gene expression in wheat stems during meiosis. DEGs were mainly participated in pathways of cell wall modification, hormone signal transduction and transcriptional regulation to promote the freezing tolerance of plants.
Sequence Read Archive (SRA) database under BioProject accession number PRJNA901148.

Author contributions
Conceptualization and supervision, GS, JM, XLiu and YF; methodology, JW and WP; investigation, JW, WP, BZ, XLW, XNW, XYW; bioinformatic analyses, DY and JW; data curation, DY; manuscript preparation, DY, JW and WP; writing-review and editing, GS and DY. All authors contributed to the article and approved the submitted version.