Comparison of Gene Expression Changes in Three Wheat Varieties with Different Susceptibilities to Heat Stress Using RNA-Seq Analysis

Wheat is highly susceptible to heat stress, which significantly reduces grain yield. In this study, we used RNA-seq technology to analyze the transcript expression at three different time-points after heat treatment in three cultivars differing in their susceptibility to heat stress: Jopum, Keumkang, and Olgeuru. A total of 11,751, 8850, and 14,711; 10,959, 7946, and 14,205; and 22,895, 13,060, and 19,408 differentially-expressed genes (log2 fold-change > 1 and FDR (padj) < 0.05) were identified in Jopum, Keumkang, and Olgeuru in the control vs. 6-h, in the control vs. 12-h, and in the 6-h vs. 12-h heat treatment, respectively. Functional enrichment analysis showed that the biological processes for DEGs, such as the cellular response to heat and oxidative stress—and including the removal of superoxide radicals and the positive regulation of superoxide dismutase activity—were significantly enriched among the three comparisons in all three cultivars. Furthermore, we investigated the differential expression patterns of reactive oxygen species (ROS)-scavenging enzymes, heat shock proteins, and heat-stress transcription factors using qRT-PCR to confirm the differences in gene expression among the three varieties under heat stress. This study contributes to a better understanding of the wheat heat-stress response at the early growth stage and the varietal differences in heat tolerance.


Introduction
Global climate change is expected to have a direct effect on crop production around the world. In particular, wheat is highly susceptible to heat stress, and it is estimated that global wheat production will decrease by 4% to 8.5%/ • C increase in temperature [1-3]. Among abiotic stress factors, heat is one of the most widely studied in plants. Thus, it is well known that heat causes morphological, physiological, and biochemical changes that lead to significant reductions in wheat yield [4]. The effect of heat on plants depends on the duration of the exposure to heat and the specific plant growth-stage at which the stress is experienced [5][6][7][8][9][10]. The optimum temperature range for wheat flowering and grain filling is 12-22 • C [11]. However, an increase in temperature of 1-2 • C shortens the grain filling period, resulting in lower grain weight [12]. In addition, grain yield can be significantly reduced by short-term exposure to high (>35 • C) temperatures [13]. Overall, heat stress disrupts productivity at all developmental stages, including embryonic cells [14], the early stages of meiosis [15], microspore and pollen cell development at the time of floral initiation [16], grain filling [17], root growth [18], and the survival of productive tillers [19], all of which result in the reduction of grain yield. Additionally, heat causes the interruption of protein folding and synthesis [11], leading to the production and accumulation of stressing agents that immediately disrupt major cellular metabolic processes, as well as DNA replication, transcription, and mRNA transport and translation until the cells recuperate [20]. Furthermore, heat shock causes changes in membrane potential (depolarization), lipid (susceptibility) and 'Sakai 143' (susceptibility); 'Keumkang' was obtained by crossing varieties of 'Sakai 75' (resistance) and 'Geuru' (susceptibility); and 'Jopum' was inherited from 'Sakai 75' (resistance) [37]. In this study, we used these varieties and an RNA-seq expression analysis to compare their gene expression profiles in young seedlings growing under heat stress. First, to confirm the heat response of young seedlings of the three varieties, we observed the extent of leaf drying after repeated heat treatment for 3 h at 45 • C and their subsequent recovery three days after return to normal temperature conditions. Differences in the green index among the three varieties were not significant after the first or second heat treatments; however, after the third heat treatment, the green index of OL decreased compared to those of JP and KK, and after the fourth heat treatment, JP showed only a few dry leaves, KK looked normal, and OL showed many dry leaves ( Figure 1A); concomitantly, green leaves of untreated plants from the three cultivars became more abundant as plants continued to grow ( Figure 1B). Seedlings of OL showed a significant decrease in green index compared with those of JP ( Figure 1C). Consistent with the abovementioned classification of the tested varieties in terms of their level of heat tolerance, our results indicate that, with respect to young seedlings-as might be expected-the tested cultivars ranked in the following order in extent of leaf senescence observed in response to heat stress: JP < KK < OL.

Evaluation of the Heat-Stress Response in the Wheat Varieties Tested
Previously, Ko et al. [37] reported the effects of high temperature on chlorophyll content, grain weight, and shoot dry weight at the early grain filling stage in 11 Korean wheat varieties. Based on the rate of reduction observed in these parameters, cultivars of Jopum (JP), Keumkang (KK), and Olgeuru (OL) were classified as heat tolerant, moderately heat tolerant, and heat susceptible, respectively [37]. 'OL' was obtained by crossing 'Geuru' (susceptibility) and 'Sakai 143' (susceptibility); 'Keumkang' was obtained by crossing varieties of 'Sakai 75' (resistance) and 'Geuru' (susceptibility); and 'Jopum' was inherited from 'Sakai 75' (resistance) [37]. In this study, we used these varieties and an RNA-seq expression analysis to compare their gene expression profiles in young seedlings growing under heat stress. First, to confirm the heat response of young seedlings of the three varieties, we observed the extent of leaf drying after repeated heat treatment for 3 h at 45 °C and their subsequent recovery three days after return to normal temperature conditions. Differences in the green index among the three varieties were not significant after the first or second heat treatments; however, after the third heat treatment, the green index of OL decreased compared to those of JP and KK, and after the fourth heat treatment, JP showed only a few dry leaves, KK looked normal, and OL showed many dry leaves ( Figure. 1A); concomitantly, green leaves of untreated plants from the three cultivars became more abundant as plants continued to grow ( Figure 1B). Seedlings of OL showed a significant decrease in green index compared with those of JP ( Figure 1C). Consistent with the abovementioned classification of the tested varieties in terms of their level of heat tolerance, our results indicate that, with respect to young seedlings-as might be expected-the tested cultivars ranked in the following order in extent of leaf senescence observed in response to heat stress: JP < KK < OL. The 1st, 2nd, 3rd, and 4th indicate repeating 1st, 2nd, 3rd, and 4th heat treatments at 45 °C for 3 h. Asterisks indicate significant differences at p < 0.01 (**) as calculated using Student's t-test. KK, Keumkang; OL, Olgeuru. The 1st, 2nd, 3rd, and 4th indicate repeating 1st, 2nd, 3rd, and 4th heat treatments at 45 • C for 3 h. Asterisks indicate significant differences at p < 0.01 (**) as calculated using Student's t-test.

Differential Expression and Cluster Analyses of Inducible Genes by RNA-Seq
To identify DEGs in 6-h (H6) and 12-h (H12) heat-treated seedlings of JP, KK, and OL relative to the untreated control seedlings (CT), we calculated the gene transcript abundance using Kallisto (v. 0.46.0) [38] (Table S1) and identified DEGs in R (v. 3.6.3) by setting a threshold of log2 fold-change > 1 and FDR (padj) < 0.05. Differentiallyexpressed genes (DEG) were identified by pairwise comparisons of the three experimental conditions in all three cultivar datasets separately. As shown in Figure 2 and Table S2

Differential Expression and Cluster Analyses of Inducible Genes by RNA-Seq
To identify DEGs in 6-h (H6) and 12-h (H12) heat-treated seedlings of JP, KK, and OL relative to the untreated control seedlings (CT), we calculated the gene transcript abundance using Kallisto (v. 0.46.0) [38] (Table S1) and identified DEGs in R (v. 3.6.3) by setting a threshold of log2 fold-change > 1 and FDR (padj) < 0.05. Differentially-expressed genes (DEG) were identified by pairwise comparisons of the three experimental conditions in all three cultivar datasets separately. As shown in Figure 2 and   To provide further insight into the genes whose expression changed after heat treatment, we used Trinity (version 2.14.0) [39] to group DEGs into five clusters according to the change in their expression pattern ( Figure S1). The results showed that genes in cluster 1 showed a decreasing expression at H6 and H12 compared to the CT group. Specifically, genes in cluster 4 showed significant downregulation at 6 h after heat treatment, followed by upregulation thereafter. In contrast, genes in cluster 2 consistently showed upregulated or constant levels of expression. Meanwhile, genes in cluster 3 were upregulated at H6, and then downregulated at H12. Most DEGs identified in OL were grouped in cluster 1 (66.14%), followed by those grouped in cluster 2 (27.18%), with very few DEGs in clusters 3-5 (Supplementary File S1). Conversely, most DEGs identified in JP were grouped in cluster 1, followed by those grouped in cluster 2 (Supplementary File S2). Interestingly, in KK, DEGs were predominantly grouped in cluster 5 (78.73%), and to a lesser proportion in clusters 1 and 2 (Supplementary File S3).

Functional Annotation of the DEGs and Shared and Specific Biological Processes (BPs)
To elucidate the biological functions related to heat stress at the molecular level, genes that were differentially expressed in Jopum, Keumkang, and Olgeuru were included in a functional enrichment analysis performed using the GO databases (Table S3). GO terms were grouped using REVIGO [40] to remove the redundancy.
The analysis of significantly-enriched GO functions with respect to ROS metabolic processes (GO:1901031, GO:0000302, GO:0019430, GO:1901671) and cellular response to heat (GO:0009408, GO:0034605) were performed for all of the comparative combinations (Supplementary File S6). These genes were significantly enriched among the three comparisons as well as the three cultivars (Supplementary File S6).
We also identified 30 HSPs that were highly responsive to heat stress at 6 h from treatment initiation, compared with their basal levels of expression in the three varieties. Interestingly, most expression levels detected by RNA-seq analysis suggested downregulation ( Figure 6B, Table S7). One HSP101 (TraesCS3D02G273600.1), HSP90 (TraesCS2B02G047400.1), four HSP70 (TraesCS4A02G066100.1, TraesCS4B02G243400.1, TraesCS4D02G243000.1, and TraesCS1D02G284000.2), and eight sHSPs (TraesCS6D02G169100.1, TraesCS7A02G202200.1, TraesCS3A02G033900.1, TraesCS4B02G225400.1, TraesCS7A02G177700.1, TraesCS3D02G114900.1, TraesCS6D02G322300.1, and TraesCS4D02G226000.1) were significantly decreased after 6 h of heat treatment. Unfortunately, however, corresponding data could not be obtained for these genes because the filtering CPM (counts per million) did not exceed 2 after 12 h of heat treatment. In contrast, four small HSPs (sHSP; TraesCS1D02G319400.1, TraesCS1B02G331900.1, TraesCS1A02G319600.1, and TraesCS2B02G247300.1) increased to a similar extent at 6 h in the three varieties but showed a pattern similar to that of the control group after 12 h of heat treatment ( Figure 6B). One HSP90, TraesCS5B02G258900.1, showed different expression patterns in JP and OL. The expression of this gene increased at 6 h in both JP and OL and continued to increase at 12 h in JP, but decreased in OL ( Figure 6B).

Quantitative Real-Time PCR Analysis to Validate Gene Expression Levels Found by RNA-Seq Analysis
To validate the gene expression levels obtained by RNA-seq analysis, 21 DEGs were selected to perform quantitative real-time polymerase chain reaction (qRT-PCR) analysis. Subsequently, we determined the transcriptional levels of ROS-scavenging genes encoding APX (TraesCS2B02G457600.1, TraesCS6A02G412900.1), CAT (TraesCS6D02G048300, TraesCS7A02G549800, and TraesCS5A02G498000.1), GPX (TraesCS2B02G604800 and TraesCS 2A02G582100.1), GST (TraesCS3D02G299900.1), GR (TraesCS6A02G383800.2), and POX (TraesCS1A02G077700.2, TraesCS1D02G096400.1, and TraesCS6B02G303900.1) (Figure 7). Two CATs (CAT6D: TraesCS6D02G048300 and, CAT-7A: TraesCS7A02G549800) and one POX (POX-1D: TraeCS1D02G096400.1) showed progressive and time-dependent expression increases in JP and KK, but with different expression patterns in OL ( Figure 7C,D,K). In all three varieties, the expression of two GPXs (GPX-2A: TraesCS2A02G582100.1 and GPX-2B: TraesCS2B02G604800) was progressively reduced in a time-dependent manner by heat treatment ( Figure 7G-H). In turn, the transcription levels of the three POX genes were significantly lower in OL than in JP or KK ( Figure 7J-L). Meanwhile, the expression of one POX (POX-6B: TraeCS6B02G303900.1) decreased progressively in both JP and KK, but showed very low basal expression levels in OL, and was observed showing a slight increase at 3 h and 6 h, and then a decrease in expression level ( Figure 7L). Lastly, the expression of other genes showed slightly different patterns depending on the variety, but the differences were not significant. Specifically, the transcriptional levels of CAT-5A and POX-1D in OL were not detected after 9 h of heat treatment ( Figure 7E,K).  Figure 7L). Lastly, the expression of other genes showed slightly different patterns depending on the variety, but the differences were not significant. Specifically, the transcriptional levels of CAT-5A and POX-1D in OL were not detected after 9 h of heat treatment ( Figure 7E,K). Transcript levels of HSP and HSF genes, which are up-or downregulated upon heat treatment, were also investigated by DEG analysis. We selected six highly-responsive genes to heat stress, as well as three previously reported genes: TaHSF3 [41], TaHsfA6f [42], and TaHSP23.9 [43] (Figure 8). The transcript levels of five HSPs, HSP101c, HSP90.1-B1, HSP90.6, HSP70-8, and TaHSP23.9, were markedly increased 3 h after heat treatment, Relative gene expression level of APXs, TraesCS2B02G457600.1 (A), and TraesCS6A02G412900.1 (B); CATs TraesCS6D02G0483600 (C); TraesCS7A02G549800 (D); and TraesCS5A02G49800.1 (E); GR, TraesCS6A02G383800.2 (F); GPXs TraesCS5A02G604800 (G); and TraesCS2A02G582100.1 (H); GST TraesCS3D02G299900.1 (I); POXs TraesCS1A02G77700.2 (J); TraesCS1D02G096400.1 (K); and TraesCS6B02G303900.1 (L). Wheat leaves were heat-treated at 35 • C during 3, 6, and 9 h. Data are representative results from three independent experiments. Error bars represent SD (n = 3). Different letters on the bars indicate significant differences between treatments (ANOVA followed by Tukey's test, p < 0.05). APX, ascorbate peroxidase; CAT, catalase; GPX, glutathione peroxidase; GST, glutathione S-transferase; GR, glutathione reductase. Transcript levels of HSP and HSF genes, which are up-or downregulated upon heat treatment, were also investigated by DEG analysis. We selected six highly-responsive genes to heat stress, as well as three previously reported genes: TaHSF3 [41], TaHsfA6f [42], and TaHSP23.9 [43] (Figure 8). The transcript levels of five HSPs, HSP101c, HSP90.1-B1, HSP90.6, HSP70-8, and TaHSP23.9, were markedly increased 3 h after heat treatment, and then gradually decreased in all three cultivars ( Figure 8E-I). Similarly, the transcript levels of HSFA-2c increased in the three varieties after 3 h of heat treatment, and then decreased in KK and OL, but not in JP ( Figure 8A). Similarly, the relative expression levels of HsfB-2b and TaHSF3 significantly increased after 3 h of heat treatment in KK and OL and then decreased, but increased up to 6 h in JP and decreased thereafter; furthermore, their expression levels were lower than those in JP and OL ( Figure 8B,C). The transcript level of TaHsfA6f was higher under basal conditions and then gradually decreased with increasing duration of heat treatment ( Figure 8D). These results indicate that the expression of HSFs and HSP genes varied depending on the variety after heat treatment, with KK showing the highest level of expression for these genes, followed by OL, and JP was the lowest. increasing duration of heat treatment ( Figure 8D). These results indicate that the expression of HSFs and HSP genes varied depending on the variety after heat treatment, with KK showing the highest level of expression for these genes, followed by OL, and JP was the lowest. . Different letters on the bars indicate significant differences between treatments (ANOVA followed by Tukey's test, p < 0.05). HSP, heat shock protein; HSF, heat stress transcription factor.

Discussion
Heat stress is generally considered the most detrimental abiotic stress factor for economically important crops, especially wheat [44]. Padaria et al. [45] compared the differences in expression profiles induced by heat treatment at different developmental stages.  .1 (H), and TaHSP23.9 (I) in wheat leaves after heat treatment at 35 • C during 3, 6, and 9 h. Data are representative of results from three independent experiments. Error bars represent SD (n = 3). Different letters on the bars indicate significant differences between treatments (ANOVA followed by Tukey's test, p < 0.05). HSP, heat shock protein; HSF, heat stress transcription factor. Transcript levels of HSP and HSF genes, which are up-or downregulated upon heat treatment, were also investigated by DEG analysis. We selected six highly-responsive genes to heat stress, as well as three previously reported genes: TaHSF3 [41], TaHsfA6f [42], and TaHSP23.9 [43] (Figure 8). The transcript levels of five HSPs, HSP101c, HSP90.1-B1, HSP90.6, HSP70-8, and TaHSP23.9, were markedly increased 3 h after heat treatment, and then gradually decreased in all three cultivars ( Figure 8E-I). Similarly, the transcript levels of HSFA-2c increased in the three varieties after 3 h of heat treatment, and then decreased in KK and OL, but not in JP ( Figure 8A). Similarly, the relative expression levels of HsfB-2b and TaHSF3 significantly increased after 3 h of heat treatment in KK and OL and then decreased, but increased up to 6 h in JP and decreased thereafter; furthermore, their expression levels were lower than those in JP and OL ( Figure 8B,C). The transcript level of TaHsfA6f was higher under basal conditions and then gradually decreased with increasing duration of heat treatment ( Figure 8D). These results indicate that the expression of HSFs and HSP genes varied depending on the variety after heat treatment, with KK showing the highest level of expression for these genes, followed by OL, and JP was the lowest.
Transcript levels of HSP and HSF genes, which are up-or downregulated upon heat treatment, were also investigated by DEG analysis. We selected six highly-responsive genes to heat stress, as well as three previously reported genes: TaHSF3 [41], TaHsfA6f [42], and TaHSP23.9 [43] (Figure 8). The transcript levels of five HSPs, HSP101c, HSP90.1-B1, HSP90.6, HSP70-8, and TaHSP23.9, were markedly increased 3 h after heat treatment, and then gradually decreased in all three cultivars ( Figure 8E-I). Similarly, the transcript levels of HSFA-2c increased in the three varieties after 3 h of heat treatment, and then decreased in KK and OL, but not in JP ( Figure 8A). Similarly, the relative expression levels of HsfB-2b and TaHSF3 significantly increased after 3 h of heat treatment in KK and OL and then decreased, but increased up to 6 h in JP and decreased thereafter; furthermore, their expression levels were lower than those in JP and OL ( Figure 8B,C). The transcript level of TaHsfA6f was higher under basal conditions and then gradually decreased with increasing duration of heat treatment ( Figure 8D). These results indicate that the expression of HSFs and HSP genes varied depending on the variety after heat treatment, with KK showing the highest level of expression for these genes, followed by OL, and JP was the lowest.

Discussion
Heat stress is generally considered the most detrimental abiotic stress factor for economically important crops, especially wheat [44]. Padaria et al. [45] compared the differences in expression profiles induced by heat treatment at different developmental stages. However, related transcriptional profiles in wheat have not been widely studied among cultivars, particularly for comparison purposes. RNA sequencing (RNA-seq) has been primarily used to identify novel and conserved stress-responsive genes associated with heat stress tolerance in wheat [33,35,45,46]. In this study, we evaluated heat-stress-induced gene expression by RNA sequencing over time after heat treatment in three cultivars differing in susceptibility to heat stress. The mapping ratio of the sample to the reference genome ranged from~50% to 86%, with an average of 80% (Table S1). To further identify DEGs at 6 h (H6) and 12 h (H12) after heat-treatment initiation relative to the untreated control (CT) seedlings, we analyzed DEGs using a volcano plot (Figure 2). There were significant differences in the number of DEGs in the H6 vs. CT, H12 vs. CT, and H6 vs. H12 groups. These results indicate the differences in wheat regulatory mechanisms in response to heat under different treatment durations. The analysis of significantly-enriched GO functions with respect to ROS metabolic processes (GO:1901031, GO:0000302, GO:0019430, GO:1901671) and cellular response to heat (GO:0009408, GO:0034605) was significantly enriched among the three comparisons and three cultivars. Thus, we identified the DEGs associated with enzymatic ROS-scavenging activity, HSPs, and HSFs that respond significantly to heat stress.
A direct result of stress-induced cellular changes is the overproduction of ROS in plant tissues, which are continuously reduced/eliminated by the plant antioxidant system, thus maintaining a steady-state level of ROS under stress conditions [47]. Consistently, an increased activity of many antioxidant enzymes has been observed in many plants [48,49]. Here, RNA-seq analysis showed that nine CATs, twelve APXs, seven POX, two SODs, five GPXs, two GSTs, and one GR were identified in the three varieties under study, which were highly responsive to heat treatment. Furthermore, these DEGs were validated by qRT-PCR. Transcript levels of the identified DEGs showed various expression patterns in the three cultivars. In particular, two CAT (CAT-6D:TraeCS6D02G048300 and CAT-7A: TraeCS7A02G549800) and one POX (POX-1D: TraeCS1D02G096400.1) were upregulated by heat treatment, with their expression levels differing by cultivar (Figure 7).
Previously, Schleiff et al. [50] reported that HSFs are primarily involved in stimulating the rapid synthesis and accumulation of HSPs. Indeed, HSFs are rapidly induced within 1 h of heat stress initiation to activate downstream pathways [44]. Thus, 56 HSF genes have been identified in wheat [51] and divided into 3 classes: HsfA, HsfB, and HsfC. However, few Hsfs have been cloned and characterized because of the complexity of the wheat genome. Particularly, TaHsfA-2d [52], TaHsfA6b [53], TaHsfA2-1 [54], TaHsfA2e-5D [55], and TaHSFA6f [43] reportedly play important roles in heat stress. Using RNA-seq analysis, here, we identified three new candidate genes that might be involved in the heat-stress response of wheat young seedlings. One from the HsfA group (HsfA-6e, TraesCS1A02G375600), one from the HsfB group (HsfB-2b, TraesCS7A02G270100.1), and one from the HsfC group (HsfC-1b, TraesCS3A02G280800.1). As per qRT-PCR analysis, the relative expression of HsfB-2b and HsfC-1b was significantly increased upon heat stress, but the expression pattern differed by cultivar. The HSF genes identified and investigated in this study suggest that they play positive roles in regulating thermotolerance, especially in wheat. Plants synthesize many stress-responsive proteins, including HSPs, in response to heat stress. These HSPs exhibit tissue-specific and developmental-stage-specific expression [34]. However, in wheat, very limited information is available regarding the HSP family, although it is by and large the most important staple food crop in the world. Here, using RNA-seq analysis, we identified a new candidate gene involved in the heat-stress response in wheat plants and evaluated it by qRT-PCR analysis.
As a result of RNA-seq and qRT-PCR analyses, we showed that the expression level of HSF and HSP genes after heat treatment differed by variety. However, the ranking of varieties for gene expression level and phenotype did not match the degree of heat sensitivity estimated. Nonetheless, these results will help understand the role of HSPs in wheat under heat stress conditions. Moreover, our findings show that, although RNA-seq methods and data analysis approaches allowed the identification of a significant number of expressed genes, there were cases in which the results observed were not consistent with the results obtained by qRT-PCR analysis, suggesting that further study is required.

Evaluation of Heat Stress in Three Wheat Varieties
To investigate leaf drying after heat treatment, a 50-hold plastic tray (54 × 28 × 6 cm) was filled with a 1:1 substrate garden bed soil (Punong, Kyoungju, Korea) and agricultural soil (Seoul Bio, Korea), and the seeds of three Korean wheat varieties [56], 'Jopum' (JP, accession no. 2002-874), 'Keumkang' (KK, accession no. 1997-4435), and 'Olgeuru' (OL, accession no. 1997-4432), were sown in a row; ten seeds were shown per variety. To ensure uniform germination, the trays with the sown seeds were incubated at 4 • C for three weeks in a plant growth chamber (Dasol Scientific, Hwaseong, Korea), and for a further two days in a controlled growth room at 22 • C under a 22-h photoperiod [57]. Trays were watered from underneath and maintained under well-watered conditions. Heat treatment was performed at 45 • C for 3 h in a growth chamber (Dasol Scientific, Korea), followed by recovery for three days in the controlled growth room. The treatment was repeated for three additional cycles of heat treatment and recovery periods, respectively, with the same set of seedlings. The experiment was repeated three times. The greenness index was calculated as a percentage by dividing the area of the green color using the formula proposed by Laohaichi et al. [58] by different numbers of the plant RGB image+ and openCV of Python 3.8 [59].

RNA Preparation and Quality Control for RNA-Seq Transcriptional Profiling
For RNA extraction, seedlings of the three cultivars, JP, KK, and OL, were grown in pots for 19 days at 22 • C under a 22-h photoperiod in the controlled growth room [57]. The seedlings were subjected to heat treatment at 35 • C [60] for 6 and 12 h; seedlings without heat treatment served as the control. Subsequently, leaves from the control and heat-treated plants were collected and rapidly frozen in liquid nitrogen and then stored in a −80 • C freezer until the RNA extraction procedure. Total RNA was extracted using TRIzol TM reagent (Thermo Scientific, Waltham, MA, USA) and the RNA quality and quantity were analyzed using a NanoDrop 1000 spectrophotometer (Thermo Scientific, MA, USA). RNA with an RNA integrity number (RIN) >7.0 was determined using a 2100 Expert Bioanalyzer (Agilent Technologies Inc., Santa Clara, CA, USA) for further use in transcriptome profiling.

cDNA Library Preparation and Illumina Sequencing
A total of 27 RNA-Seq samples from three wheat cultivars were analyzed. A quantity of 1 µg of total RNA was used to construct sequencing libraries using a Truseq Stranded mRNA library prep kit (Illumina, San Diego, CA, USA). Paired-end (2 × 100 bp) reads of samples were sequenced using the Illumina NovaSeq 6000 platform. Low-quality reads, redundant reads, and adapter sequences were eliminated using Trimmomatic (v. 0.38) by following default parameters: removing a read with an average base quality (Q20) below 20 and removing a read less than 50 bp [61]. FastQ files were generated via Illumina bcl2fastq2 (version 2.20) and the quality of trimmed reads was assessed using FastQC (version 0.11.9). To further analyze RNA-Seq data, clean reads were mapped to Triticum aestivum genome sequences [62] using Kallisto (version 0.46.0) [38]. Library construction and sequencing were performed by DNA LINK, Inc. (Seoul, Korea).

Differentially-Expressed Gene (DEG) Identification
Transcript abundance was determined using Kallisto (version 0.46.0). A differential gene expression analysis among libraries was performed in R (version 3.6.3) software using the edgeR package (version 3.30) [63]. Differentially-expressed genes (DEG) were identified by pairwise comparisons of the three experimental conditions from the Jopum, Keumkang, and Olgeuru datasets separately. The trimmed means from the M-value (TMM) normalization methods were used to calculate the normalization factors. Genes were filtered by above/below 2-fold, above CPM (counts per million) 2, and FDR (falsediscovery rate) < 0.05. The Benjamini-Hochberg method was used to estimate the FDR when identifying DEGs.

Cluster Analysis, GO Enrichment Analysis of DEGs, and Heatmap
A DEG cluster analysis was used to isolate the cluster analysis expression patterns of genes at two different sampling time-points after heat treatment of the three cultivars using Trinity (version 2.14.0 [39]. Hierarchical clustering using the Euclidean method of normalized gene expression was achieved using centralized and log2 (FPKM+1) transforms and tree cutting at 60% depth. Genes were filtered by two conditions of above/below 2-fold change and FDR < 0.05, and grouped under similar expression patterns. For ontology analysis, the Databaseor Annotation, Visualization, and Integrated Discovery (DAVID) v2022q2 [64] was used to select genes with above 2-fold change and FDR (padj) < 0.05. The genes were categorized according to their assigned ontology (GO) terms (cellular components, biological processes, and molecular functions). GO terms were grouped using REVIGO [40] to remove the redundancy. The heatmap was created in MS Excel using conditional formatting and a color scale using selected genes (Tables S5-S7).