*3.5. Comparative Analysis of Differential Genes*

We performed principal component analysis (PCA) on the gene expressions of the six samples to identify outliers and to discriminate between clusters of samples with high similarity. We found significant differences between the HT and HT+BT groups, with the two groups forming separate clusters. Moreover, the spatial distribution of the three replicates within the groups was more concentrated, which indicates that the results are reproducible (Figure 4A). In the differential gene volcano plots, we can see that the 10 mM betaine treatment induced the differential expressions of a total of 528 genes under heat stress, of which 332 genes were downregulated, and 196 genes were upregulated. This implies that betaine under heat stress promotes rice seed germination through the differential expressions of these 528 genes (Figure 4B).

**Figure 4.** (**A**) HT and HT + BT treatments were the sources of gene expression variation under heat stress. We used principal component analysis to explore the impact on the variation. (**B**) Volcano plots of differentially expressed genes (DEGs) between HT and HT + BT; each point represents a gene. The upregulated genes are represented by red dots, the downregulated genes are represented by blue dots, and the genes with no significant differences are represented by gray dots. Under the corrected Q value, the DEG is considered significantly different.

#### *3.6. GO Pathway Analysis*

Functional annotation of transcriptome sequences is an important aspect of functional genomics. GO analysis is currently the main method of functional annotation. GO is a standard classification system of gene functions, which can comprehensively describe the properties of biological genes and gene products [40]. GO analysis includes the classification annotation of GO function and the significance enrichment analysis of GO function. The GO function classification annotation can obtain the gene list of the number statistics of genes with a certain GO function. The GO function significance enrichment analysis gives the GO items of genes that are significantly enriched compared with genomic background, which can obtain the biological functions of genes. The GO database is based on GO terms, which is a tree structure with redundancy. There are three root nodes in the GO database, which describe the molecular function, cellular component and biological process of genes. The GO term is a gene set, which is generally used to find functional changes caused by different genes [41,42]. For example, the molecular function of a gene may be catalytic activity. Its cellular components are localized to the cell membrane in the cell, and the biological process involved is protein trafficking. The three levels of GO are not related to each other, and there is no definition of the mutual relationship of genes. In order to visualize the data and obtain complete functional information, we performed a GO analysis to unify the gene attributes and classify differential genes into putative functional taxa. We divided the functional annotations into three broad categories: biological processes, cellular components and molecular functions. We annotated 528 significantly differentially expressed genes (DEGs) to 40 functional groups within the three broad functional categories (Figure 5A,B). A total of 16 were biological processes, 14 were cellular components and 10 were molecular functions. The majority of DEGs in biological processes were associated with cellular processes (29.06%) and metabolic processes (28.34%); the majority of DEGs in cellular components were associated with cells (34.29%), membranes (16.67%), membrane fractions (14.87%) and organelles (18.35%); the majority of DEGs in molecular functions were associated with binding (40.28%) and catalytic activity (38.89%).

**Figure 5.** GO term comparison of differential genes with and without betaine soaking treatment under heat stress: (**A**) GO classification of differential genes; (**B**) bubble diagram of differential gene GO enrichment.

#### *3.7. KEGG Pathway Analysis*

KEGG is a signaling pathway database with an extremely rich mapping of signaling pathways and of the interactions between the genes contained in the pathways [38]. KEGG analysis allows us to visualize the expressions of genes and their regulatory pathways. According to the results, 290 DEGs were enriched in 90 pathways (Figure 6A,B). Seed germination is regulated by factors such as ambient temperature and endogenous hormones. In the signal transduction pathway of environmental information, we found that many different genes were enriched in plant hormone signal transduction, the plant mitogenactivated protein kinase (MAPK) signal pathway, carotenoid, diterpene biosynthesis and other plant hormone synthesis pathways.

**Figure 6.** KEGG term comparison of differential genes with and without betaine soaking treatment under heat stress: (**A**) KEGG classification of differential genes; (**B**) bubble diagram of differential gene KEGG enrichment.

#### *3.8. Differential Genes Related to Reactive Oxygen Species in MAPK Plant Signaling Pathway*

As an important part of cell signal transduction, ROS are inevitably produced during the normal growth and development of plants and under stress. A certain concentration of ROS is necessary during plant growth, but if excessive ROS cannot be cleared in time, it will cause oxidative damage to cells [9]. H2O2, as a signal substance, can enter cells through aquaporins. In cells, H2O2 triggers reactive oxygen species (ROS) signals and responses, thus interfering with ROS stabilization [10]. MAPK signaling acts on receptor-like protein kinases (RLKs) and downstream of ROS signals to regulate ROS-related gene expression and programmed cell death (PCD) [43]. ROS signaling is closely related to the MAPK cascade pathway. H2O2 can induce and activate some MAPKs, and the MAPK cascade pathway can also regulate the dynamic balance of ROS [44]. In tobacco, the two signaling pathways, NPK1–MEK1–Ntf6 and Mek2–SIPK, can promote the accumulation of ROS by inducing the expression of the *NbRbohb* gene, thus enhancing its stress resistance [45]. In *Arabidopsis*, *AtMPK8* is also involved in regulating the dynamic balance of ROS, and the MKK1–MPK6 signaling pathway is involved in inducing *CAT1* expression and H2O2 production in the ABA signaling pathway [46,47].

Heat stress produces an oxidative stress response in seeds, which leads to the production of large amounts of reactive oxygen species, which ultimately inhibit seed germination. Therefore, reducing the accumulation of excess reactive oxygen species due to heat stress is an effective way to improve seed germination rates under stress. Our previous results show that 10 mM of betaine soaking significantly increased the activity of some antioxidant enzymes in seeds under heat stress and reduced the MDA content to promote rice seed germination. In addition to the enzymatic aspects, according to our transcriptomic results (Figure 7A), betaine infiltration significantly downregulated the expression level of *MEKK1* in the H2O2 signaling pathway in the plant MAPK signaling pathway. It further reduces the production of reactive oxygen species under heat stress and improves seed adaptation to high temperatures.

**Figure 7.** Expression profiles of DEGs in reactive oxygen species accumulation and three endogenous hormone-related pathways: (**A**) expression profiles of DEGs related to the H2O2 signal transduction pathway under HT and HT + BT treatments; (**B**) expression profiles of DEGs related to the GA pathway under HT and HT + BT treatments; (**C**) expression profiles of DEGs related to the ABA pathway under HT and HT + BT treatments; (**D**) expression profiles of DEGs related to the IAA pathway under HT and HT + BT treatments. The sample names are shown at the bottom of the figure. We indicate changes in expression levels by changes in color: green indicates a lower expression level, whereas red indicates a higher expression level. All data shown reflect the average mean of three biological replicates (*n* = 3). Means with different letters in each treatment represent a significant difference of *p* ≤ 0.05.

#### *3.9. Differential Genes Related to Hormone Signaling Pathway*

To further explore the action pathway between betaine and plant hormones, we analyzed the differentially expressed genes related to plant hormones. Among the differential genes related to GAs, one is expressed in the GA synthesis pathway, and two are expressed in the GA signal transduction pathway (Figure 7B).

According to the gene expression heat map, betaine infiltration under heat stress significantly downregulated the expression level of gibberellin (*Gibberellin 2-beta-dioxygenase* (*OsGA2ox3*)) [48] and inhibited GA inactivation. We verified this by the significantly higher GA content in the betaine treatment compared with the heat-stress treatment. The betaine treatment decreased the expression level of *DELLA* [49] in the GA signal pathway, and it enhanced the output of the GA signal. In the ABA synthesis pathway, betaine induced the decrease in the expression level of *9-cis-epoxycarotenoid dioxygenase* (*NCED*) [50], the rate-limiting enzyme of ABA synthesis, and inhibited the synthesis of ABA (Figure 7C). However, at the same time, betaine also induced the downregulation of the *(+)—abscisic acid 8 –hydroxylase* (*CYP707A*) gene [51] in the ABA metabolic pathway, slowing down the metabolic rate of ABA.

In addition, in the ABA signal transduction pathway, the expression levels of three genes of *protein phosphatase 2c* (*PP2C*) [52] were downregulated in the betaine treatment, and the ABA output signal was weakened. In the IAA signal transduction pathway, we observed the expressions of two IAA auxin early-response genes (*GH3*) [53], in which the expression level of *TLD1* [54] in the HT treatment was significantly lower than that in the HT+BT treatment, while *GH3–4* [53] was significantly higher than that in the HT+BT treatment (Figure 7D). Moreover, the expression level of *SAUR* [55] was significantly higher in the HT+BT treatment than in the HT treatment. The high expression of *TLD1* promotes IAA inactivation and reduces the level of IAA in seeds. These results are consistent with the assay results of the phytohormone levels.

#### *3.10. Validation of DEGs by qRT-PCR*

To verify the accuracy and reproducibility of the RNA-seq results, we randomly selected nine genes in the pathways related to seed germination for qRT-PCR validation (Figure 8). We calculated the expression levels of the selected genes using the 2−ΔΔCt method. According to the qRT-PCR results, the expression levels of nine DEGs were consistent with the results of the RNA-seq, confirming the reproducibility of the data.

**Figure 8.** The expressions of 9 differentially expressed genes in the H2O2 pathway and 3 endogenous hormone pathways were randomly verified by qRT-PCR. The x-axis is the name of the treatment, in which HT is the 38 ◦C non-soaking treatment, and HT + BT is the 38 ◦C betaine seed-soaking treatment. The y-axis shows the relative expression of a specific gene relative to the reference gene actin. The last panel shows the results of the RNA-seq\qRT-PCR correlation analysis (lower right corner). The x-axis is the log2-fold-change value of the qRT-PCR, and the y-axis is the log2-foldchange value of the RNA-seq. According to the Pearson correlation (*R*<sup>2</sup> = 0.60; *p* < 0.05), there is a significant positive correlation between the multiple changes in expression. Each parameter \* indicates significant statistical difference, \* (*p* < 0.05).
