Uncovering the Gene Regulatory Network of Maize Hybrid ZD309 under Heat Stress by Transcriptomic and Metabolomic Analysis

Maize is an important cereal crop but is sensitive to heat stress, which significantly restricts its grain yield. To explore the molecular mechanism of maize heat tolerance, a heat-tolerant hybrid ZD309 and its parental lines (H39_1 and M189) were subjected to heat stress, followed by transcriptomic and metabolomic analyses. After six-day-heat treatment, the growth of ZD309 and its parental lines were suppressed, showing dwarf stature and rolled leaf compared with the control plants. ZD309 exhibited vigorous growth; however, M189 displayed superior heat tolerance. By transcriptomic and metabolomic analysis, hundreds to thousands of differentially expressed genes (DEGs) and metabolites (DEMs) were identified. Notably, the female parent H39 shares more DEGs and DEMs with the hybrid ZD309, indicating more genetic gain derived from the female instead of the male. A total of 299 heat shock genes detected among three genotypes were greatly aggregated in sugar transmembrane transporter activity, plasma membrane, photosynthesis, protein processing in the endoplasmic reticulum, cysteine, and methionine metabolism. A total of 150 heat-responsive metabolites detected among three genotypes were highly accumulated, including jasmonic acid, amino acids, sugar, flavonoids, coumarin, and organic acids. Integrating transcriptomic and metabolomic assays revealed that plant hormone signal transduction, cysteine, and methionine metabolism, and α-linolenic acid metabolism play crucial roles in heat tolerance in maize. Our research will be facilitated to identify essential heat tolerance genes in maize, thereby contributing to breeding heat resistance maize varieties.


Introduction
High temperature as a consequence of global warming has become a common concern worldwide that threatens to crop production [1][2][3]. Maize (Zea mays L.) is a well-known important crop for food, fuel, and feed [4] but is sensitive to heat stress, decreasing grain yield [2,5]. Prior data have shown that with an average temperature increase of 1 • C, maize grain yield decreased by an average of 7.4% [6]. Heat stress causes additional negative effects on maize growth and development, such as dwarf plant stature, reduced pollen viability, and poor seed set [2,3,5,7]. Adaptation to temperature gradients is the primary restriction of maize to growth around the world, regardless of latitude and altitude [4,8,9].

Estimation of the Effects of Heat Stress on Grain Yield
In the Huang-Huai-Hai River Basin area of China, the growth of maize is suffered from severe climate conditions, particularly extreme high temperatures. To estimate the heat stress effect on ZD309 grain yield, we selected 12 commercial hybrids to be subjected to heat stress (HS) treatment in the greenhouse ( Table 1). The heat-tolerant coefficient was calculated by HS grain yield/CK grain yield. Under heat stress, kernel development was repressed in other hybrids, with almost no apparent impact on the ZD309. Within the 12 tested hybrids, ZD309, DK653, Zhongkeyu505, and Zhengdan319 exhibited superior heat tolerance. However, heat stress caused notable grain yield loss in Zhengdan1002, Zhengdan326, Pionner335, and Zhengdan958.

Evaluation of the Morphological and Physiological Effects of Heat Stress on Seeding Development
To identify the effects of heat stress on ZD309 seedlings, five-week-old potted planted seedlings were subjected to heat stress treatment. After six-day-heat stress, the growth of ZD309 and its parental lines were greatly suppressed ( Figure 1A-C), exhibiting dwarf stature and rolled leaf in comparison to the corresponding controls. The ZD309 exhibited the most negative effects on stature, along with the least effects on M189 among three tested materials. ZD309 showed enhanced growth vigor compared with its parental lines. Inconsistent with the morphological changes, heat stress also mediated physiological changes ( Figure 1D-M). Among the 10 tested physiological parameters, the proline content in ZD309 and its parental lines were increased in heat stress treatment for three days (D3HS) but decreased in heat stress treatment for six days (D6HS). The malondialdehyde (MDA, an indicator of lipid peroxidation), H 2 O 2 , and soluble sugar contents in the three tested materials were increased in D3HS treatment and sustainedly increased in D6HS treatment. The chlorophyll a contents in the three tested materials were decreased in D3HS treatment and continued to be decreased in D6HS treatment. The soluble protein contents, SOD, POD, and CAT enzyme activities of ZD309 and its parental lines were increased in both heat treatments and climbed a higher level in the D3HS treatment. The chlorophyll b levels of ZD309 and its parental lines were increased in both heat stress treatments; however, no significant differences between D3HS and D6HS treatments. In both treatments, all parameters were the highest in ZD309 but the lowest in M189. The contents of MDA and H 2 O 2 in M189 in both treatments were the highest but the lowest in ZD309. The chlorophyll b contents showed no significant difference among the three tested materials. Overall, ZD309 and its parental lines displayed superior heat tolerance but varied from the morphological and physiological parameters.

DEGs of ZD309 and Its Parental Lines in Response to Heat Stress
Using the Illumina paired-end RNA-seq approach, 1228 million reads were detected, and 184.21 gigabases (Gb) data were generated. After removing the low-quality data, a total of 181.03 Gb of clean data were produced. The sequencing data have been submitted to the NCBI Short Read Archive with accession number (PRJNA791560).
The transcriptome data showed time-dependent gene expression patterns under heat treatments, and the heat stress response of ZD309 was different from those of its parental lines. The number of up-regulated DEGs is more than down-regulated DEGs at each time point in all three materials. This indicated that more genes were switched on than switched off under heat stress. In D3HS treatment, 5440 (3797 up-and 1643 down-regulated), 3200 (1667 up-and 1533 down-regulated), and 4348 (2861 up-and 1487 down-regulated) DEGs were identified within the heat-stressed H39_1, M189, ZD309, and their controls, respectively ( Figure 2A). In D6HS treatment, 3921 (2815 up-and 1106 down-regulated), 2424 (1486 up-and 938 down-regulated), and 2391 (1266 up-and 1125 down-regulated) DEGs were identified within the heat-stressed H39_1, M189, ZD309, and their controls, respectively ( Figure 2B). In D3HS treatment, 1259 common DEGs were found among the three tested materials (Figure 2A). Moreover, 2805 genes were overlapped between the inbred line H39_1 and the hybrid ZD309, 1772 genes were overlapped between M189 and ZD309, as well as 1667 genes were overlapped between H39_1 and M189. In D6HS treatment, 614 common DEGs were found ( Figure 2B). A total of 1129, 1029, and 1028 genes were shared between H39_1 and ZD309, M189, and Z309, as well as H39_1 and M189, respectively. Compared with the D6HS treatment, more DEGs were detected in ZD309 and its parental lines in the D3HS treatment. This suggested that ZD309 and its parental lines start to adapt to the stress as the treatment goes on. Overlapped DEGs between ZD309 and H39_1 were more than those of M189, which implied H39_1 contributed larger genetic gains to ZD309 with regard to heat tolerance. , as well as SOD, POD, and CAT activity. For the content of proline, soluble sugar content, BCA protein, chlorophyll a/b, the measurement unit, mg g −1 , represents the respective content in 1 g of sample's fresh weight; For the content of malondialdehyde and H 2 O 2 , the measurement unit, µmol mg −1 , represents the respective content in 1 mg of sample's fresh weight; For the activity of SOD and POD, the measurement unit, U mg −1 , represents the respective activity units detected in 1 mg of sample's fresh weight; For CAT activity, the measurement unit, U min −1 mg −1 , 1 activity unit (U) represents the A240 absorbance decreased 0.1 OD in 1 min and 1 mg of sample's fresh weight. * and ** represent that the corresponding physical character in heat-stressed maize plants are significantly and very significantly different from the control at p < 0.05 and p < 0.01 levels, respectively. To interpret heat tolerance-associated molecular pathways in ZD309, all DEGs were subjected to GO and KEGG enrichment analyses ( Figures 2C,D and 3). For GO analysis, in both D3HS and D6HS treatments, the common DEGs were strongly enriched in photosynthesis, photosystem I, photosystem II, plasma membrane, chlorophyll-binding, chloroplast thylakoid membrane, and protein-chromophore linkage terms. The unique DEGs in D3HS treatment were in protein binding, mitochondrion, zinc ion binding, RNA modification, and response to light stimulus. The unique DEGs in D6HS treatment were in carbohydrate metabolic process, plastid, thylakoid, response to heat, carbohydrate transport. For KEGG analysis, the DEGs for D3HS treatment were mainly involved in photosynthesis-antenna proteins, cysteine and methionine metabolism, ABC transporters, starch and sucrose metabolism, and protein processing in endoplasmic reticulum pathways. The DEGs for D6HS treatment were mainly involved in photosynthesis, benzoxazinoid biosynthesis, carbon fixation in photosynthetic organisms, glycolysis/gluconeogenesis, and protein processing in endoplasmic reticulum pathways. These results illustrated that heat stress has strong effects on photosynthesis, protein processing, and starch and sucrose metabolism. The overlapped DEGs among ZD309 and its parental lines in both D3HS and D6HS treatments were centered for further analysis ( Figure 4A). A total of 299 DEGs (198 upand 93 down-regulated) were suggested to be important for heat tolerance. These DEGs were enriched in photosynthesis-antenna proteins, protein processing in the endoplasmic reticulum, cyanoamino acid metabolism, and phagosome pathways ( Figure 4B). Based on the gene expression patterns, these genes were classified into four groups, continued up-regulated genes, continued down-regulated genes, prolonged heat stress suppressed up-regulated genes, and prolonged heat stress mitigated up-regulated genes ( Figure 4C). To correlate the RNA-seq results, the expression level of 10 randomly selected DEGs were analyzed by qRT-PCR ( Figure 5). Among the 198 up-regulated DEGs, 102 genes exhibited over-dominant effects (a higher expression level in ZD309 than in both parental lines), which is possible to explain the heat tolerance heterosis in the hybrid. Of these over-dominant effect genes, genes encoding for chaperone protein ClpB1, NADP-dependent oxidoreductase encoding gene, and Hsp90 were enhanced in long time heat stress. Genes encoding for NAC transcription factor and zinc finger family proteins, Zm00001d002905 (encode a protein kinase) and Zm00001d048710 (encode benzoxazinone synthesis 2) were also enhanced in the hybrid.
In addition, genes encoding for lysine histidine transporter, triose phosphate isomerase, protein kinases, and abscisic acid-inducible protein were activated by long-term heat stress. These highly up-regulated genes are involved in multiple biological processes, indicating that heat tolerance heterosis is associated with sophisticated regulatory networks.  Figure 6A). Of these DEGs, 1982 DEGs were solely expressed in the three tested materials, including 1171 up-regulated and 684 down-regulated genes ( Figure 6B). The unique DEGs were subjected to KEGG pathway enrichment ( Figure 6C). The upregulated DEGs in ZD309 and its parental lines were strongly enriched in alanine and aspartate and glutamate metabolism, benzoxazinoid biosynthesis, cysteine and methionine metabolism, and ribosome biogenesis in eukaryotes. The up-regulated DEGs in H39_1 were preferably enriched in alpha-linolenic acid metabolism, amino sugar and nucleotide sugar metabolism, valine and leucine and isoleucine biosynthesis, amino sugar and nucleotide sugar metabolism, flavonoid biosynthesis, phenylalanine and tyrosine and tryptophan biosynthesis, and pentose phosphate. The up-regulated DEGs in M189 were preferably enriched in ABC transporters and starch and sucrose metabolism. Indeed, the up-regulated DEGs in ZD309 were preferably enriched in alpha-linolenic acid metabolism, carbon fixation in photosynthetic organisms, glycolysis/gluconeogenesis, linoleic acid metabolism, purine metabolism, phenylalanine and tyrosine and typtophan biosynthesis, fructose and mannose metabolism, and tyrosine metabolism. On the other hand, the down-regulated DEGs in ZD309 and its parental lines were greatly enriched in either photosynthesis or photosynthesis-antenna proteins. In comparison to the parental lines, the down-regulated DEGs in ZD309 were mainly enriched in brassinosteroid biosynthesis, limonene, and pinene degradation. Overall, genes associated with stress response (such as response to biotic stimulus, response to water deprivation, and response to hydrogen peroxide) were preferably enriched in up-regulated DEGs in ZD309. DEGs in inbred line H39_1 share more common pathways with ZD309 than inbred line M189.

The Expression of HSFs and HSPs Genes in the DEGs
Of the detected DEGs, 15 hsf-tfs (HSFs encoding genes) were differentially expressed under heat treatments ( Table 2). In D3HS treatment, most hsf-tfs were down-regulated or insignificant-expressed. However, hsf-tf18 was enhanced in both ZD309 and H39_1; hsf-tf6 was elevated in M189; hsf-tf17 and hsf-tf9 were up-regulated specifically in H39_1; and hsf-tf7 did solely in ZD309. Interestingly, hsf-tf13 was down-regulated in the three tested materials. Compared with the D3HS treatment, more hsf-tfs were enhanced in D6HS treatment, suggesting that prolonged heat stress triggers more hsf-tfs to be differently expressed. hsf-tf18 and hsf-tf7 were up-regulated in ZD309 and its parental lines. hsf-tf20, hsf-tf24 and hsf-tf28 were up-regulated in ZD309. Most hsf-tfs were up-regulated in ZD309, followed by H39_1, then the least for M189. Note: The expression levels of the HSF and HSP/heat-response proteins encoding genes are marked in different colors. Red represents the significant up-regulated gene expression, log 2 (fold changes) > 1. Blue represents the significant down-regulated gene expression, log 2 (fold changes) < −1.

Metabolic Analysis of in ZD309 and Its Parental Lines Response to Heat Stress
To examine the metabolic response to heat stress in ZD309 and its parental lines, leaf samples were collected in D3HS and D6HS treatments for untargeted metabolite analysis based on ultra-performance liquid chromatography/mass spectrometry (UPLC/MS) platform. PLS-DA (partial least squares discriminant analysis) was performed to identify the heat-response metabolites through differential abundance. The differentially expressed metabolites (DEMs) were defined as a criterion of p-value ≤ 0.05, |log2FC (fold changes)| ≥ 1, and VIP (variable importance in projection) score ≥ 1. The metabolite data have been deposited into the CNGB Sequence Archive (CNSA) of China National GeneBank DataBase (CNGBdb) with accession number CNP0002688. In D3HS treatment, 1542, 1150, and 1401 DEMs were identified in the heat-stressed H39_1, M189, ZD309 versus the corresponding controls, respectively ( Figure 7A). Of these DEMs, 754 DEMs were shared between H39_1 and ZD309, 396 DEMs shared between M189 and ZD309, as well as 357 DEMs between H39_1 and M189. A total of 201 overlapped DEMs were found among the three tested materials. In D6HS treatment, 1782, 1801, and 2374 DEMs were identified in the heat-stressed H39_1, M189CK, ZD309CK versus the corresponding controls, respectively ( Figure 7B). Of these DEMs, 1295 DEMs were shared between H39_1 and ZD309, 1102 DEMs shared between M189 and ZD309, as well as 946 DEMs between H39_1 and M189. A total of 747 overlapped DEMs were found among the three tested materials ( Figure 7B). Of these detected DEMs, 350, 311, and 423 DEMs were up-regulated in H39_1, M189, and ZD309, respectively ( Figure 7C). A total of 150 DEMs were found across the tested materials, which might be essential for heat tolerance. Under heat stress, some metabolites were enhanced, including sugars and glycosides (e.g., mannose 6-phosphate, isopentenyladenine-9-N-glucoside, and cytarabine), methionine and organic acids (e.g., jasmonic acid (M211T234), dihydrojasmonic acid, and 9-Oxo-10(E), 12(E)-octadecadienoic acid), nucleosides and nucleotides (e.g., cytosine and 5-Methylcytosine) ( Figure 7D, Table 3). In contrast, other metabolites, amino acids (e.g., lysine, isoleucine, and aspartic acid), and alkaloids were reduced ( Figure 7D, Table 3). Furthermore, 41, 71, and 90 metabolites were found to be associated with heat stress specifically in H39_1, M189, and ZD309, respectively. After heat treatment, JA abundance was significantly increased in three tested materials versus the corresponding controls. JA abundance was also significantly increased in D6HS treatment compared with that in D3HS treatment.

Integrative Analysis of Gene Expression and Metabolic Changes under Heat Stress Conditions
To decipher the transcriptional regulation of metabolic networks, a co-enriched transcriptome and metabolome pathway analysis was conducted. In D3HS treatment, 38, 26, 26 co-enriched pathways were found in H39_1, M189, and ZD309, including starch and sucrose metabolism, arginine and proline metabolism, isoquinoline alkaloid biosynthesis, and glutathione metabolism. In D6HS treatment, 36, 22, and 37 co-enriched were found in H39_1, M189, and ZD309, including carbon fixation in photosynthetic organisms, phenylpropanoid biosynthesis, alpha-linolenic acid metabolism, and glycerophospholipid metabolism. Importantly, six co-enriched pathways were detected in both heat stress treatments, including plant hormone signal transduction, cysteine and methionine metabolism, alanine and aspartate and glutamate metabolism, ABC transporters, zeatin biosynthesis, and cyanoamino acid metabolism. These pathways might play critical roles in conferring maize heat tolerance.
Most DEGs were also identified in these co-enriched pathways. The expressions of these DEGs were highly consistent with the contents of corresponding metabolites. For instance, mannose-6-phosphate, an intermediate of fructose and mannose metabolism, was elevated among the three tested materials. Most fructose and mannose metabolism-related genes were activated under heat stress, with a remarkable level in the inbred line H39_1 (Figure 8).

The Expression of Genes Involved in JA Biosynthesis under Heat Stress
Combined transcriptome and metabolome analyses have demonstrated that JA is a key regulator to confer heat tolerance. Upon metabolomics, JA raised significantly under heat stress. Compared with M189, the JA contents (M211T234) were increased in H39_1 and ZD309 ( Figure 9A). The JA levels were enhanced along with the heat stress. The JA contents (M211T234) in D3HS treated H39_1, M189, ZD309 were increased 4.8, 2.5, and 4.6 folds versus the corresponding controls, respectively. The JA (M211T234) contents in D6HS treated H39_1, M189, ZD309 were increased 21.7, 27.9, and 73.3 folds versus the corresponding controls, respectively. The expression pattern of JA biosynthesis (αlinolenic acid metabolism) genes under heat stress was further analyzed ( Figure 9B). Overall, most genes were up-regulated among three materials, particularly in H39_1 and ZD309. Zm00001d053675 (encoding LOX) were up-regulated among the three tested materials, with the highest level in D3HS-treated ZD309 plants. However, six maize AOS encoding genes showed a higher level in H39_1 and ZD309 materials of D3HS treatment. Two AOC encoding genes, Zm00001d029594 and Zm00001d047340, were up-regulated among the three tested materials under heat stress, with a relatively lower level of Zm00001d029594. Two OPCL1 encoding genes, Zm00001d027519 and Zm00001d043376, were also up-regulated among the tested materials. In contrast, Zm00001d043376 showed a higher level in D6HS than D3HS treatments. Moreover, Zm00001d035854, encoding α-dioxygenase, was the least level in ZD309 in D6HS treatment.

ZD309 Is a Heat-Tolerant Maize Hybrid
As a sessile organism, maize suffers from several climate stresses, such as salinity, drought, heat, lodging, and flooding. Heat stress particularly hinders maize grain yield in the Huang-Huai-Hai River Basin area of China. In plants, heat stress generates excess reactive oxygen species (ROS) in the cellular compartments, such as chloroplasts, mitochondria, plasma membranes, peroxisomes, apoplast, and ER [2]. Excessive ROS alters a series of physiological responses, including photosynthesis, respiration, transpiration, membrane thermostability, and osmotic regulation [3]. Finally, the development and productivity of plants are largely inhibited [3]. In the present study, ZD309 was screened for less yield loss under heat stress and demonstrated superior heat tolerance. Further heat treatments on the ZD309 and its parental lines seedlings resulted in multiple morphological and physiological alterations. Consistent with the previous studies [1,3,[23][24][25][26][27], heat stress interrupted plant growth and yielded increased levels of proline, MDA, H2O2, soluble sugar, soluble protein, SOD, POD, and CAT, as well as decreased level of chlorophyll a contents. Although ZD309 and its parental lines indicated superior heat tolerance, they had different mechanisms to minimize heat stress. The plant height of ZD309 and H39_1 were more affected by heat stress, but they had relatively stronger growth vigor. The growth of M189 was less affected, but it produced a relatively higher content of ROS content and low levels of SOD, POD, and CAT activities. These results indicated that M189 has a stronger heat tolerance compared with ZD309 and H39_1. ZD309 and H39_1 likely minimized the effects of heat stress through growth adaptation.

Important Metabolites in Response to Heat Stress
Abiotic stresses can also cause metabolic alterations in plants, such as maintaining essential metabolism and synthesizing protective metabolites to protect themselves from stress damage. Heat stress is a major abiotic stress in crop production. By untargeted metabolomics analysis, metabolic reprogramming was identified in several crops, such as tomato, maize, barley, wheat, and soybean [16,21,[39][40][41]. Hundreds of metabolites were indeed identified that are closely associated with heat stress, such as malate, valine, isoleucine, glucose, starch, sucrose, proline, glycine, serine, drummondol, anthranilate, dimethylmaleate, galactoglycerol, guanine, and glycerone. Heat stress also caused endogenous hormone changes, including zeatin, gibberellin acid3 (GA3), gibberellin acid4 (GA4), abscisic acid (ABA), salicylic acid (SA), and jasmonic acid (JA) [15]. In the present study, 150 DEMs were identified in the tested three materials. Of these DEMs, mannose 6-phosphate, isopentenyladenine-9-N-glucoside, cytarabine, JA, octadecadienoic acid, cytosine, and 5-methylcytosine were highly accumulated under heat stress. Conversely, lysine, isoleucine, aspartic acid, alkaloids were significantly reduced. These enhanced metabolites might contribute to the heat tolerance of ZD309 and its parental lines. Most DEMs were correlated with DEGs, which mainly co-enriched in several GO terms, including plant hormone signal transduction, cysteine, and methionine metabolism, alanine and aspartate and glutamate metabolism, ABC transporters, zeatin biosynthesis, and cyanoamino acid metabolism. The key metabolites can be used to enhance maize heat tolerance through exogenous metabolites feeding.

Potential Applications of HSFs, HSPs, and Metabolites in Heat Stress Tolerance Improvement
With an increasing number of essential heat-response genes characterized, the underlying regulatory network has been partly resolved in model plants. However, there is still a gap for heat resilient variety breeding in crop plants. In the present study, some major HSFs, HSPs, and metabolites were identified through combined transcriptomic and metabolomic analysis. Of them, these identified HSFs, HSPs genes can be used for screening ideal haplotypes, thereby selecting heat-tolerant inbred lines and creating heat-tolerant hybrid varieties. In previous studies, the SA187 [42] and biochar [43] application enhanced heat tolerance in Arabidopsis and rice, respectively. Similarly, the identified metabolites can be used for improving maize tolerance by exogenous application. Taken together, our results will not only be a benefit for uncovering the mechanism of the heat-response network but also can provide the information for heat tolerance improvement.

Plant Materials and Heat Treatment
A total of 12 commercial maize hybrids, including ZD309 (Henan Academy of Agricultural Sciences, China, dent grain, summer maize), were potted planted at Yuanyang Farm of Henan Academy of Agricultural Sciences. The experimental soil belongs to yellow-cinnamon soil, which contains total N for 0.74 g kg −1 , available P for 42.6 mg kg −1 , 163.0 mg kg −1 , and organic matter for 33.0 g kg −1 (pH 8.83). To estimate heat stress effects on different maize hybrids, 50 six-week-old plants from each hybrid were selected under consistent growth conditions. The selected plants were divided into two groups, with 25 plants each. One group was moved to a sunlight greenhouse for heat stress treatment, and the other group was placed in the field as the control. In the sunlight greenhouse, the daytime setting temperature was 35 • C (actual temperature 35.0-43.6 • C, outdoor air temperature 32.5-38.7 • C), the nighttime setting temperature was 28 • C (actual temperature 28.0-33.7 • C, outdoor air temperature 24.5-30.4 • C). On average, the temperature in the greenhouse was 3-5 • C higher than the field. All plants were well watered every two days to ensure no water stress. At the kernel maturity stage, the ears from the two groups were harvested for estimating the grain yield. The potted experiment was repeated three times.
The hybrid ZD309 and its parental lines, H39_1 and M189, were germinated and grown in the greenhouse at 28/24 • C (12/12 h dark/light cycle) with 60% relative humidity and 1000 µmol m −2 s −1 light density. The uniform of five-week-old maize seedlings were selected and divided into two groups and three replications, with 24 seedlings for each group and 8 seedlings for each replication. One group was moved into a growth chamber with a 28/24 • C (day/night) cycle as the control (CK). The other group was moved into a growth chamber with a 40/35 • C(day/night) cycle as the heat treatment (HT) stress. The plants were well watered every two days to ensure no water stress. The uppermost whole leaves of maize seedlings were sampled after heat treatment for 3 days (D3HS) and 6 days (D6HS), frozen in liquid nitrogen, and stored at −80 • C for the following experiments. Each sample had three biological replicates.

Heat Stress-Related Physiological Parameters Measurement
The leaf samples of ZD309 and its parental lines were subjected to physiological parameter determination. The contents of activity of proline, malondialdehyde (MDA), H 2 O 2 , soluble sugar, soluble protein, chlorophyll, superoxide dismutase (SOD), peroxidase (POD), and catalase (CAT) were measured using the corresponding testing kit from Suzhou Comin Biotechnology Co., Ltd. (Suzhou, China. Available online: http://www.cominbio. com/a/shijihe/, accessed on 1 January 2022) by following the manufacturer's instructions.

mRNA Library Construction and Illumina Sequencing
Total RNA was extracted from maize leaves of the three materials using Trizol reagent (Invitrogen, CA, USA) according to the manufacturer's instructions. The total RNA quality was evaluated by Bioanalyzer 2100 and RNA 1000 Nano LabChip Kit (Agilent, CA, USA) with a RIN number > 7.0. Poly(A) RNA was first purified from total RNA (5 ug) with poly-T-oligo-attached magnetic beads for two rounds of purification. The full-length mRNA was next fragmented into small pieces using divalent cations under elevated temperature. The cleaved RNA fragments were then reverse transcribed into cDNA for library construction by using the mRNA-Seq sample preparation kit (Illumina, San Diego, CA, USA). The average insert size for the paired-end libraries was 300 bp (±50 bp). Paired-end sequencing was performed on the Illumina Hiseq 4000 platform (LC Sciences, Houston, TX, USA).

Quantitative Real-Time PCR
To verify the reproducibility of DEGs in the RNA-seq data, we used the same samples for evaluating the gene expression levels by qRT-PCR. Ten genes from the DEGs were selected randomly, and corresponding primers were designed by NCBI (Primer-Blast). The gene name and primer sequence of the DEGs are listed in Table S1. Total RNA of the tested samples in RNA-seq reverse transcribed using the PrimeScript™ RT reagent kit with gDNA Eraser (Perfect Real Time) (Takara, Dalian, China). Then, qRT-PCR was performed by using the Gene Applied Biosystems@7500 Fast and TransStart Top Green qPCR SuperMix (TransGen Biotech, Beijing, China). Each sample was analyzed with three technical replicates. The PCR cycle was configured as: 5 min at 95 • C followed by 40 cycles of amplification at 95 • C for 15 s, then 58 • C for 20 s, and 72 • C for 30 s. The relative expression level of the DEGs was calculated by the 2 −∆∆Ct [49]. The GAPDH gene was used as the internal control [50].

Metabolites and Metabolome Analysis
The samples were placed on ice, and the metabolites were extracted with 50% methanol buffer [51]. Briefly, 100 mg of sample was ground into powder and added 120 µL of precooled 50% methanol, vortexed for 1 min, and incubated at room temperature for 10 min; the extraction mixture was then stored overnight at −20 • C. After centrifugation at 4000× g for 20 min, the supernatants were transferred into new 96-well plates and stored at −80 • C refrigerator for further UPLC-MS analysis. In addition, pooled QC samples were also prepared by combining 10 µL of each extraction mixture. A quality control (QC) sample was included to check the repeatability within an analytical batch.
The MS data were analyzed using XCMS and CAMERA software (available online: https://bioconductor.org/, accessed on 1 January 2022) [52]. To obtain high-quality data, robust Loess signal correction (R-LSC) was used to normalize peaks from the tested samples with regard to the QC samples, and the relative standard deviation (RSD) was calculated. Those features that were detected in less than 50% of QC samples or 80% of biological samples were removed, then the remaining peaks with missing values were imputed with the k-nearest neighbor algorithm to further improve the data quality. The m/z with RSD value between 0% and 30% was subjected to principal component analysis (PCA), partial least squares discriminant analysis (PLS-DA), and variable importance in the projection (VIP) analysis. The m/z with VIP > 1, fold change > 1.5 (or <0.67) and Q value < 0.05 in any sample was retained for further metabolite identification (at identification level 2) and pathway annotation by using the online HMDB database (available online: http: //www.hmdb.ca/, accessed on 1 January 2022), METLIN (available online: http://metlin. scripps.edu/, accessed on 1 January 2022) and KEGG (www.genome.jp/kegg/, accessed on 1 January 2022) [53,54].

Statistical Analysis
All data collected from phenotypic, physiological parameters and qRT-PCR analyses were subjected to one-way variance analysis (ANOVA) and Student's t-test using software SPSS 22.0 (IBM, New York, NY, USA). p < 0.05 indicates the statistical differences to reach the significant different level, p < 0.01 for very significant different levels.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11050677/s1. Table S1. Primers used for qRT-PCR analysis in this study; Table S2. GO enrichment for DEGs identified in transcriptome sequencing; Table S3. Significant DEMs identified in metabolomic analysis.