Physiological, Epigenetic, and Transcriptome Analyses Provide Insights into the Responses of Wheat Seedling Leaves to Different Water Depths under Flooding Conditions

Flooding stress, including waterlogging and submergence, is one of the major abiotic stresses that seriously affects the growth and development of plants. In the present study, physiological, epigenetic, and transcriptomic analyses were performed in wheat seedling leaves under waterlogging (WL), half submergence (HS), and full submergence (FS) treatments. The results demonstrate that FS increased the leaves’ hydrogen peroxide (H2O2) and malondialdehyde (MDA) contents and reduced their chlorophyll contents (SPAD), photosynthetic efficiency (Fv/Fm), and shoot dry weight more than HS and WL. In addition, FS increased catalase (CAT) and peroxidase (POD) activities more than HS and WL. However, there were no significant differences in the contents of H2O2, MDA, SPAD, and Fv/Fm, and the activities of superoxide dismutase (SOD) and POD between the HS and WL treatments. The changes in DNA methylation were related to stress types, increasing under the WL and HS treatments and decreasing under the FS treatment. Additionally, a total of 9996, 10,619, and 24,949 genes were differentially expressed under the WL, HS, and FS treatments, respectively, among which the ‘photosynthesis’, ‘phenylpropanoid biosynthesis’, and ‘plant hormone signal transduction’ pathways were extensively enriched under the three flooding treatments. The genes involved in these pathways showed flooding-type-specific expression. Moreover, flooding-type-specific responses were observed in the three conditions, including the enrichment of specific TFs and response pathways. These results will contribute to a better understanding of the molecular mechanisms underlying the responses of wheat seedling leaves to flooding stress and provide valuable genetic and epigenetic information for breeding flood-tolerant varieties of wheat.


Introduction
Flooding stress is a major ecological threat that restricts crop growth and yield in high-rainfall zones across the globe [1].Based on the depth of the water table, flooding can be classified as waterlogging, when it is superficial and covers only the roots, or as submergence, when water completely covers the aerial plant tissues [2].Flooding can cause the accumulation of reactive oxygen species (ROS), such as superoxide (O 2 •− ) and hydrogen peroxide (H 2 O 2 ), increasing lipid peroxidation and resulting in plant cell damage [3,4].Flooding also results in a reduction in stomatal conductance, the CO 2 assimilation rate, photosynthesis rate, and nutritional imbalance [5].With the onset of climate change, the average and extreme precipitation intensities and the frequency of extremely heavy precipitation have significantly increased [6].Therefore, understanding the response mechanisms of plants to flooding stress has important implications for developing flood-tolerant varieties and promoting their effective adaptation to climate change.
In recent decades, the physiological and molecular mechanisms of adaptation to flooding stress have been reported [7][8][9][10][11][12].For example, increases in the contents of superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), and other antioxidant enzymes, have been suggested as a key strategy through which plants can effectively resist waterlogging stress [8].For example, γ-aminobutyric acid was found to enhance waterlogging tolerance in maize by increasing the activities of several antioxidant enzymes (SOD, POD, CAT, etc.) [13].A transcriptomic analysis showed that Chinese wingnut plants enhanced their waterlogging tolerance by increasing their synthesis of alpha-linolenic acids and flavonoids and activating plant hormone signaling pathways [8].In cucumber, waterlogging-induced differentially expressed genes were especially related to enhanced glycolysis, adventitious root development, and amino acid metabolism [11].In Phalaris arundinacea, waterlogging-stress-induced differentially expressed genes were mainly involved in carbohydrate metabolism, hormone signaling regulation, and the scavenging of reactive oxygen species [12].Moreover, TFs, including MYB, bHLH, NAC, WRKY, ERF, and bZIP, are also reported to be the major regulators involved in waterlogging tolerance [12].For instance, the ERFs of group VII are well known for facilitating ethylene signal transduction and enhancing waterlogging tolerance in maize [14] and wheat [15].
Epigenetic alterations have been found to be associated with abiotic stresses in crop species [16][17][18][19].DNA methylation is one of the critical epigenetic mechanisms for the regulation of gene expression under abiotic stresses [17,18,20,21].In wheat, DNA demethylation significantly increases waterlogging-related gene expression in tolerant genotypes under hypoxic stress, such as the expression of ERF1, ACC1, and CKX2.3 [18].In rice, DNA methylation in response to drought stress regulates the expression of unique genes responsible for drought stress tolerance [16,19].Furthermore, the response of DNA methylation to stress is also related to the type of stress.For example, total DNA methylation increases under drought stress in sesame but decreases under waterlogging stress [17].These studies show that epigenetic regulation is important in the abiotic stress responses in plants.However, the potential mechanisms underlying epigenetic regulation under flooding stress remain largely unclear.
In recent years, water depth has been recognized as one of the important factors limiting the growth and development of flooded plants [22][23][24].For example, the leaves of mulberry seedlings under shallow submergence were healthy, while the leaves of mulberry seedlings treated with half submergence and full submergence showed waterlogging symptoms, to varying degrees, in their middle [22].Partially submerged Paspalum dilatatum plants showed a reduction in the starch concentration in their leaves, but their biomass was unaffected, whereas completely submerged plants did not survive [23].All Melilotus siculus accessions were able to reorient petioles towards the vertical axis under both partial and full submergence [24].However, petiole extension rates were maintained under the partial submergence treatment and decreased under the full submergence treatment.Moreover, the sugar contents of Melilotus siculus leaflets rose during partial submergence but were depleted during full submergence [24].These studies implied that different performance and response mechanisms exist for the response to flooding stresses at different water depths, but the underlying adaptation mechanisms remain largely unknown.
Wheat (Triticum aestivum L.) is an important food crop, and major source of calories, that is grown worldwide [25].In the case of wheat, the seedling stage is one of the important stages at which waterlogging is most detrimental to wheat yields, after germination [9].Approximately 15-20% of the total wheat cropping area is under threat of flooding, which could reduce the global grain yield by up to 43% [26,27].Thus far, physiological, transcriptional, and epigenetic regulatory responses to waterlogging/hypoxic stress in wheat have been investigated [9,15,18].However, whether different flooding depths affect the physiology, transcription, and epigenetics of wheat is a question that is yet to be fully answered, and one that is related to the screening and evaluation of plant flooding tolerance and the breeding of resistant varieties.Thus, in this study, to gain insights into the molecular responses of wheat to flooding, three types of flooding stress treatments, based on water depth, were designed: soil surface submergence (the water was superficial and covered only the roots and soil surface), half submergence (the water covered half of the plant), and full submergence (the water completely covered the plant tissues).We systematically investigated the changes in the physiology, epigenetic regulation (DNA methylation), and genomic transcription of wheat seedling leaves under controlled conditions and three types of flooding conditions.This study aims to provide a comprehensive research framework for understanding the physiological changes, epigenetic regulation, and global transcription in wheat under flooding stress, which is of great significance for the genetic improvement and resistance evaluation of flood-resistant wheat varieties.

The Distinct Effects of FS and HS/WL on Biomass and Photosynthesis
After a 7-day treatment, the Efumai 1 plants under FS treatment showed wilting and yellowing, whereas the plants under the WL and HS treatments had slightly wilted leaves (Figure 1A).The shoot dry weight significantly decreased by 58.26% under the FS treatment compared to the control plants, followed by the HS (40.10%) and WL (38.70%) treatments (Figure 1B).Their SPAD values significantly decreased by 26.33%, 33.41%, and 64.41%, and their Fv/Fm values significantly decreased by 20.46%, 23.85%, and 34.87%, under the WL, HS, and FS treatments, respectively (Figure 1C,D).Interestingly, no differences were observed between WL and HS.
of flooding, which could reduce the global grain yield by up to 43% [26,27].Thus far, physiological, transcriptional, and epigenetic regulatory responses to waterlogging/hypoxic stress in wheat have been investigated [9,15,18].However, whether different flooding depths affect the physiology, transcription, and epigenetics of wheat is a question that is yet to be fully answered, and one that is related to the screening and evaluation of plant flooding tolerance and the breeding of resistant varieties.Thus, in this study, to gain insights into the molecular responses of wheat to flooding, three types of flooding stress treatments, based on water depth, were designed: soil surface submergence (the water was superficial and covered only the roots and soil surface), half submergence (the water covered half of the plant), and full submergence (the water completely covered the plant tissues).We systematically investigated the changes in the physiology, epigenetic regulation (DNA methylation), and genomic transcription of wheat seedling leaves under controlled conditions and three types of flooding conditions.This study aims to provide a comprehensive research framework for understanding the physiological changes, epigenetic regulation, and global transcription in wheat under flooding stress, which is of great significance for the genetic improvement and resistance evaluation of flood-resistant wheat varieties.

The Distinct Effects of FS and HS/WL on Biomass and Photosynthesis
After a 7-day treatment, the Efumai 1 plants under FS treatment showed wilting and yellowing, whereas the plants under the WL and HS treatments had slightly wilted leaves (Figure 1A).The shoot dry weight significantly decreased by 58.26% under the FS treatment compared to the control plants, followed by the HS (40.10%) and WL (38.70%) treatments (Figure 1B).Their SPAD values significantly decreased by 26.33%, 33.41%, and 64.41%, and their Fv/Fm values significantly decreased by 20.46%, 23.85%, and 34.87%, under the WL, HS, and FS treatments, respectively (Figure 1C,D).Interestingly, no differences were observed between WL and HS.

Physiochemical Changes under the Three Flooding Conditions
The H 2 O 2 content of each sample was measured, and the results show that the FS treatment resulted in a significant increase in H 2 O 2 content (by 44.85%), followed by the WL (19.72%) and HS (19.60%) treatments, compared with the control (CK) (Figure 2A).Similarly, the MDA content significantly increased under the three flooding conditions, and the highest MDA content was found in the FS treatment (344.62% increase), followed by the HS (159.07%) and WL (90.10%) treatments (Figure 2B).Compared with CK, the highest CAT (72.18% increase) and POD (303.13%)activities were also found in the FS treatment (Figure 2D,E), while the highest SOD activity was identified in the WL (55.44%) and HS (62.21%) treatments (Figure 2C).Notably, except for CAT, the H 2 O 2 and MDA contents and SOD activities did not differ between the HS and WL treatments (Figure 2).

Physiochemical Changes under the Three Flooding Conditions
The H2O2 content of each sample was measured, and the results show that the FS treatment resulted in a significant increase in H2O2 content (by 44.85%), followed by the WL (19.72%) and HS (19.60%) treatments, compared with the control (CK) (Figure 2A).Similarly, the MDA content significantly increased under the three flooding conditions, and the highest MDA content was found in the FS treatment (344.62% increase), followed by the HS (159.07%) and WL (90.10%) treatments (Figure 2B).Compared with CK, the highest CAT (72.18% increase) and POD (303.13%)activities were also found in the FS treatment (Figure 2D,E), while the highest SOD activity was identified in the WL (55.44%) and HS (62.21%) treatments (Figure 2C).Notably, except for CAT, the H2O2 and MDA contents and SOD activities did not differ between the HS and WL treatments (Figure 2).

Global Changes in DNA Methylation under the CK, WL, HS, and FS Conditions
A total of 453 clear and reproducible bands were successfully obtained under the CK, WL, HS, and FS conditions (Figure S1).Most of the CCGG sites were shown to be largely methylated, with values ranging between 86.75% and 93.38% (Table 1), and there was a greater number of fully methylated bands (average 65.45%) (types III and IV) than hemi-methylated bands (type II) (average 23.95%).Compared with the CK (88.30%), the total DNA methylation level was slightly increased to 89.18% and 93.38% under the WL and HS conditions, respectively.This ratio decreased to 86.75% under the FS treatment (Table 1).To further investigate the difference in wheat DNA methylation in response to the WL, HS, and FS treatments, a total of 16 alternative band patterns between the CK and the flooding stress treatments were identified, which could be classified into three groups: no change, hypomethylation, and hypermethylation (Table 2).Compared with the CK, 69.50% of the detected loci showed an altered DNA methylation status under the HS conditions, followed by 51.11% under FS, and 43.79% under WL.The percentages of bands hypermethylated were 27.09% and 52.52% under the WL and HS treatments, respectively, higher than that under FS (25.33%).More hypomethylation events were detected in the FS (25.78%) than in the WL (16.70%) and HS conditions (16.97%).These results suggest the divergent reprogramming of the methylation pattern under the different types of flooding stresses.After filtering, an average of 81.90 million clean reads were obtained from each sample.The average GC value was 57.00%.The Q30 percentage of the 12 libraries varied from 93.38% to 93.75%.More than 84% of the clean reads were uniquely mapped to the wheat reference genome sequence (Table S1).A principal component analysis (PCA) of all the genes in the 12 libraries revealed a clear separation between the control and flooding treatments.However, the WL-and HS-treated samples were closest to each other, indicating the highest degree of similarity in their transcription patterns (Figure S2A).Furthermore, a Pearson correlation analysis showed a high correlation among the biological replicates (the Pearson correlation coefficient ranged from 0.97 to 0.99) (Figure S2B).These results indicate the good reproducibility and quality of the RNA-seq data.
Differential expression analysis between the flooding-treated and control samples revealed a total of 9996, 10,619, and 24,949 DEGs in the WL, HS, and FS conditions, respectively (Figure 3A and Table S2).A Venn diagram showed that a total of 32,708 DEGs (17,252 upregulated and 15,456 downregulated) were identified across the three flooding conditions (Figure 3B,C).A cross-comparison between the different gene sets showed that the highest proportion of DEGs (52.66%) was unique to the FS treatment, while the WL (9.88%) and HS (7.82%) treatments exhibited fewer uniquely expressed genes.Furthermore, 2342 upregulated and 821 downregulated DEGs were detected under all three flooding conditions, which suggested that these genes might play an important role in the functioning of wheat under flooding stress.To validate the reliability of the expression profiles obtained using RNA-Seq, sixty DEGs were randomly selected to perform qRT-PCR on.Pearson's correlation coefficients showed that the qRT-PCR and RNA sequencing data for these genes were highly correlated (r = 0.91) (Figure 3D).

KEGG Enrichment Analysis of the WL-, HS-and FS-Induced DEGs
KEGG enrichment analysis showed that a total of 36, 37, and 66 pathways were significantly enriched (p value < 0.05) under the WL, HS, and FS conditions, respectively (Tables S3-S5).The top 25 significantly enriched KEGG pathways are presented in Figure 4E-G, in which 24 of the enriched pathways-including the 'biosynthesis of secondary metabolites', 'metabolic pathways', 'phenylpropanoid biosynthesis', 'plant hormone signal transduction', 'starch and sucrose metabolism', 'flavonoid biosynthesis', and 'photosynthesis'-were common to all three flooding treatments.Moreover, a number of DEGs were enriched in four pathways unique to the WL treatment-'brassinosteroid biosynthesis', 'biotin metabolism', 'ether lipid metabolism', and 'glycosaminoglycan degradation' (Figure 4E and Table S3)-whereas three pathways were specific to the HS treatment-'taurine and hypotaurine metabolism', 'circadian rhythm-plant', and 'selenocompound metabolism' (Figure 4F and Table S4).Moreover, 32 specific pathways were significantly enriched in response to the FS treatment, including 'carbon metabolism', 'glyoxylate and dicarboxylate metabolism', 'carbon fixation in photosynthetic organisms', 'porphyrin and chlorophyll metabolism', and 'peroxisome' (Figure 4G and Table S5).

GO Enrichment Analysis of the WL-, HS-and FS-Induced DEGs
GO classification analysis showed that there were similar processes enriched in the wheat leaves exposed to SS, HS, and FS, including metabolic process and cellular process in the biological process category, binding and catalytic activity in the molecular function category, and membrane and membrane parts in the cellular component category.However, the genes enriched in these GO terms were mainly upregulated under SS and HS, while they were mainly downregulated under FS (Figure S3).Furthermore, GO enrichment analysis showed that 116 and 64 GO terms were significantly enriched in the upregulated and downregulated genes under WL conditions, respectively (Tables S12 and S13).Of these, 'catalytic activity', 'extracellular region', 'cell wall', 'external encapsulating structure', and 'hydrolase activity, acting on glycosyl bonds' were the top five significantly enriched GO terms for the genes upregulated under WL conditions (Table S12), and 'glucosidase activity', 'sucrose alpha-glucosidase activity', 'pyruvate kinase activity', 'potassium ion binding', and 'alkali metal ion binding' were the top five significantly enriched GO terms for the genes downregulated under WL conditions (Table S13).A total 83 and 143 GO terms were significantly enriched in the HS-induced upregulated and downregulated genes, respectively (Tables S14 and S15).Of these, 'transferase activity, transferring hexosyl groups', 'solute:cation symporter activity', 'solute: proton symporter activity', 'symporter activity', and 'transferase activity, transferring glycosyl groups' were the top five significantly enriched GO terms in the upregulated genes (Table S14).Meanwhile the 'lipid biosynthetic process', 'catalytic activity', 'transferase activity, transferring acyl groups other than aminoacyl groups', 'transferase activity, transferring acyl groups', and 'fatty acid biosynthetic process' were the top five significantly enriched GO terms for the downregulated genes (Table S15).Under the FS condition, a total of 153 and 194 GO terms were significantly enriched in the upregulated and downregulated genes, respectively (Tables S16 and S17).Of these, 'protein kinase activity', 'catalytic activity', 'kinase activity', 'phosphotransferase activity, alcohol group as acceptor', and 'transferase activity' were the top five significantly enriched GO terms for the upregulated genes (Table S16), and 'photosynthesis', 'thylakoid plastid', 'plastid part', and 'thylakoid membrane' were the top five significantly enriched GO terms for the upregulated genes (Table S17).

Analysis of the DEGs Related to the Photosynthesis Pathway in Response to WL, HS, and FS Stresses
Two photosynthesis pathways, 'photosynthesis-antenna proteins' and 'photosynthesis', were significantly enriched in all flooding conditions in the KEGG analysis (Tables S3-S5), while 'photosynthesis-antenna proteins' were only significantly enriched under HS and FS (Tables S3-S5).In total, 155 DEGs were related to KEGG photosynthesis pathways (Figure 4A and Table S6), including 35 genes related to photosystem I (PSI) (Figure 4C), 68 genes related to photosystem II (PSII) (Figure 4D), 8 genes related to cytochrome b6/f complex (Figure 4E), 31 genes related to photosynthetic electron transport (Figure 4F), and 13 genes related to F-type ATPase (Figure 4G).A total of 76 DEGs were involved in the 'photosynthesis-antenna proteins' pathway (Figure 4B and Table S7), including 17 genes related to light-harvesting chlorophyll protein complex I (LHCI) (Figure 4H) and 59 genes related to LHCII (Figure 4I).Most of the DEGs involved in the both the 'photosynthesis-antenna proteins' and 'photosynthesis' pathways were significantly downregulated, especially under the FS treatment (Figure 4).
significantly downregulated, especially under the FS treatment (Figure 4).

Expression of the Plant Hormone Signal Transduction Pathway and Transcription Factor Genes in Response to Different Flooding Stresses
In addition to the phenylpropanoid biosynthesis pathway, KEGG enrichment results also showed that genes related to the plant hormone signal transduction pathway were significantly enriched in all the flooding treatments.A total of 173 DEGs were identified to be involved in eight plant hormone signal transduction pathways, including those for auxin (IAA), cytokinin (CTK), gibberellin (GA), abscisic acid (ABA), ethylene (ETH), brassinosteroid (BR), jasmonic acid (JA), and salicylic acid (SA) (Figure 6A and Table S20).
WL condition.Nine TGAs and twelve PR-1s were only downregulated or upregulated under the FS treatment.Furthermore, three NPR1s, one TGA, and two PR-1s were coupregulated under both the HS and FS conditions (Figure 6B and Table S20).
A total of 12, 2, 16, and 6 DEGs were identified in the CK, GA, ETH, and JA signaling pathways, respectively, and most of these were only differentially expressed under FS treatment (Figure 6B and Table S20).Most of the genes in the CTK, GA, and JA signaling pathways were downregulated, while the 16 ETH-related DEGs were upregulated under the FS treatment.
For the BR signal, four genes, encoding one BR-signaling kinase (BSK) and three xyloglucan-xyloglucosyl transferase TCH4 proteins (TCH4s), were detected (Figure 6B and Table S20).Two TCH4s (TraesCS7A02G427600 and TraesCS7B02G327700) were co-upregulated under all the flooding conditions.TraesCS1B02G192300, encoding a BSK, was only downregulated under the FS treatment (Figure 6B and Table S20).
For the SA signal, 34 DEGs, including 6 regulatory NPR1 proteins (NPR1s), 11 TGA transcription factors, and 17 pathogenesis-related protein 1s (PR-1s), were detected in this study (Figure 6B and Table S20).Most of these genes were upregulated under one or more flooding conditions.Five genes, including one NPR1 and four PR-1s, were co-upregulated in all flooding conditions, while two NPR1s were only differentially expressed under the WL condition.Nine TGAs and twelve PR-1s were only downregulated or upregulated under the FS treatment.Furthermore, three NPR1s, one TGA, and two PR-1s were coupregulated under both the HS and FS conditions (Figure 6B and Table S20).

The Submergence of Functional Leaves Was the Key Factor Determining Photosynthetic Efficiency after Submergence
Flooding usually reduces the rate of photosynthesis in plants [7].Chlorophyll fluorescence is suggested to be a sensitive indicator of the stress-induced damage to photosystem II, and a reduction in F v /F m is a good indicator of the photosynthetic impairment resulting from waterlogging stress [4].Previous studies have found that waterlogging stress significantly reduces SPAD and F v /F m values [4,10].Moreover, waterlogging stress damages the PSII reaction center of mulberry seedlings (F o , F m , and F v /F o ), reducing their ability to accept electrons, and the degree of damage is proportional to the submergence depth [22].Similar results were also observed in the present study, in which FS induced a greater reduction in the SPAD and F v /F m values of wheat leaves and led to a higher loss of biomass than HS and WL (Figure 1), which is consistent with the results of the transcriptome analysis, where the DEGs involved in 'photosynthesis-antenna proteins' and 'photosynthesis' pathways were significantly downregulated, especially under the FS treatment (Figure 4).Furthermore, more DEGs involved in the photosystem II process were downregulated than in other sub-pathways of the 'photosystem' pathway (Figure 4D).Meanwhile, no significant difference in these parameters was detected between the WL and HS treatments (Figure 1C,D), which is consistent with a previous study showing that the biomass (shoot and root) of Paspalum dilatatum plants did not change after 30 days of waterlogging and partial submergence treatments, whereas completely submerged plants did not survive [23].Additionally, a significant difference in plant weight between partial and full submergence conditions was also detected in Melilotus siculus plants, and the plant weight showed a greater reduction under the FS condition than under the partial submergence condition [24].These results may indicate that the effect of flooding on phenotypic changes may be the same when functional leaves are not submerged.

The Role of Antioxidant Mechanisms in Response to the Three Types of Flooding Conditions
Flooding conditions lead to hypoxic and anoxic conditions for plants and induce the production of ROS in plant cells, and ROS can directly attack membrane lipids, leading to the accumulation of MDA [4,7,13].MDA is commonly used as a marker of lipid peroxidation, reflecting the degree of cell membrane damage incurred in response to different environmental stresses [28].For instance, waterlogging stress significantly increased the contents of O 2 •− , • OH, and H 2 O 2 in maize leaves, leading to an accumulation of MDA [13].Wang et al. (2020) found that, compared to salt stress, saline-alkali stress more significantly upregulated MDA levels in Triarrhena sacchariflora leaves, indicating that saline-alkali stress caused greater harm to yellow horn seedlings [29].In the present study, a significant increase in H 2 O 2 content was observed under the FS treatment, leading to a higher content of MDA in wheat leaves, which clearly indicated the existence of oxidative stress and greater damage in response to FS stress than when plants were exposed to WL and HS stresses (Figures 1 and 2).
Antioxidant enzymes, including SOD, POD, CAT, APX, and GPX, play vital roles in eliminating ROS [4,7,30].The high activities of these enzymes under stress can improve stress tolerance in plants [4,13,30,31].SOD catalyzes the conversion of peroxide anions to H 2 O 2 and O 2 , whereas POD, CAT, and APX catalyze the conversion of H 2 O 2 to oxygen and water [29].In this study, SOD, POD, and CAT activities significantly increased under all the flooding stresses compared to the control.The result was similar to those of previous reports on soybean [32], peach [4], maize [13], Triarrhena sacchariflora [7], and wheat [30] under waterlogging stress.Thus, the expression of these antioxidant enzymes may be induced by flooding stress, and these enzymes may be critical for the survival of wheat leaves under flooding conditions [7].Interestingly, a discrepancy between the gene expression levels and the corresponding enzyme activities was observed.For example, the activity of SOD increases only slightly, but the corresponding gene expression dramatically increases (Figure 2C and Table S19).These results may be due to the fact that the translation of expressed mRNA into active proteins is a complicated process that may be affected by multiple factors, such as post-transcriptional modifications, the stoichiometric balance of protein biosynthesis, and protein degradation [33].The expressed mRNA is not necessarily translated into a similar amount of detected enzyme activity [34].The activities of CAT and POD under FS were significantly higher than those under HS and WL.Furthermore, in terms of the H 2 O 2 and MDA contents and SOD and POD activities, there was no significant difference between the WL and HS treatments.Similar results were also observed for the expression of ROS-scavenging-related genes (Tables S18 and S19).These results suggest that the ROS regulation mechanisms induced by WL/HS and FS are different.

Genes Related to Phenylpropanoid Biosynthesis Were Differentially Regulated in Response to Different Flood Stresses in the Wheat Seedling Stage
In plants, the phenylpropanoid biosynthesis pathway is suggested to be vital for adaptation to environmental stresses [35].The accumulation of phenolic compounds through the activation of this pathway plays an important role in neutralizing harmful ROS and protecting the plant from oxidative damage caused by ROS [19].Under waterlogging stress, the phenylpropanoid biosynthesis pathway was frequently enriched in different species, such as Triarrhena sacchariflora [7], alfalfa [10], wheat [9], cucumber [11], and Phalaris arundinacea [12].Consistent with these studies, the phenylpropanoid biosynthesis pathway was significantly enriched under the three flooding conditions in this study (Figure 5), and 21 genes were co-upregulated by WL, HS, and FS, suggesting that the activation of the phenylpropanoid pathway may be an important mechanism underlying plant responses to flooding stress.
Three preliminary steps of the phenylpropanoid biosynthesis pathway, catalyzed by phenylalanine ammonia-lyase, 4-coumarate CoA ligase, and trans-cinnamate 4-monooxygenase, are essential for subsequent branching and regulating the production of antioxidant phenolic compounds such as flavonoids, lignins, and tannins.A previous study demonstrated that waterlogging induced the high expression of isochorismate synthase-and phenylalanine ammonia lyase-related genes in wheat and enhanced the SA content in the roots, which promoted the formation of axile roots and surface adventitious roots [36].Consistent with this, a total of 7, 13, and 23 phenylalanine ammonia-lyase genes were upregulated by WL, HS, and FS, and 9 of these were uniquely regulated by FS (Figure 5B).Additionally, two and three 4-coumarate CoA ligase genes were uniquely upregulated and downregulated under the WL and FS treatments, respectively.Three genes encoding trans-cinnamate 4-monooxygenase were only downregulated under the FS condition (Figure 5B).Lignin is the main component of the plant cell wall and is one of the most important products in the phenylpropanoid biosynthesis pathway, which plays an important role in the response of wheat to waterlogging stress [37].Caffeic acid O-methyltransferase and caffeoyl-CoA O-methyltransferase are two important enzymes that participate in lignin biosynthesis, controlling the syringyl and guaiacyl units of the lignin polymer, respectively [38].The S/G ratio is demonstrated to be a major determinant of lignin quality [39].The present investigation showed that three caffeic acid 3-O-methyltransferases were specifically upregulated under WL stress but downregulated under FS stress (Figure 5B).Two caffeoyl-CoA O-methyltransferase genes were uniquely upregulated under FS stress (Figure 5B).Shikimate O-hydroxycinnamoyltransferase is involved in the production of methoxylated monolignols, which are precursors to syringyl-and guaiacyl-unit lignins.In the present study, three DEGs encoding shikimate O-hydroxycinnamoyltransferases were only upregulated by FS (Figure 5B).Additionally, it was also observed that three 5-O-4coumaroyl-D-quinate 3 -monooxygenases, which catalyze the conversion of p-coumaryl CoA to p-coumaryl shikimate, were uniquely downregulated under the FS condition (Figure 5B).These results suggest that the synthesis of different phenolic components may be an important mechanism of the plant response to different types of flooding.However, further research is needed to confirm this hypothesis.

Specific Responses of Plant Hormone Signaling Genes to Different Depths of Flooding Stress
Plant hormones, including IAA, CTK, GA, ABA, ETH, BRs, JA, and SA, play an important role in the stress resistance of many plants [39].Many genes involved in plant hormone signaling were found to be affected in flooded plants [11,12,40].ABA is widely recognized as an important hormone and chemical signal in plants responding to waterlogging stress [41].Under waterlogging stress, stomatal closure, presumably regulated by abscisic acid (ABA), is triggered to prevent dehydration [42].The Arabidopsis abi2-1 mutant, which has a higher stomatal transpiration due to impaired ABA signaling, showed enhanced plant survival after submergence [43].Waterlogging stress can promote the accumulation of ABA, which leads to a reduction in stomatal conductance in leaves [7].In the ABA signaling pathway, PYR/PYL, the most upstream regulator in this pathway, interacts with PP2C to reduce the inhibition of SnRK2, thereby regulating downstream factors such as ABF [44].When plants are exposed to waterlogging, the DEGs regulating PYL are significantly activated and upregulated [39].The level of endogenous IAA was increased in Prunus persica leaves under waterlogging stress [45].AUX/IAA and SAUR are the two major classes of primary auxin response genes.The genes in the IAA signaling pathway are negatively regulated by waterlogging [8].In plants, SA is a common phenolic compound.Under waterlogging conditions, the content of SA was significantly increased in waterlogging-tolerant soybean cultivars, further stimulating adventitious roots, enhancing gas exchange, and ultimately conferring tolerance to waterlogging stress [46].Exogenous SA promotes the formation of axile roots and surface adventitious roots in wheat under waterlogged conditions [34].Chen et al. (2022) found that a decreased SA content caused by low-expressed PAL might impact the resistance of Styrax tonkinensis seedlings to waterlogging stress [47].Zhang et al. (2017) further suggested that JA, SA, and BR are involved in the waterlogging response, given that many differentially expressed genes associated with JA, SA, and BR were upregulated in stressed cotton plants [48].Consistent with these studies, the plant hormone transduction pathway was significantly enriched in wheat seedlings under the three flooding conditions (Figure 6B and Table S8).Most of the DEGs associated with IAA, GA, ABA, BRs, JA, and SA were upregulated under one or more flooding conditions, highlighting the importance of plant hormones in the plant responses to flooding stress.
On the other hand, a large number of plant hormone signal-related genes were flooding-type-specific.For example, GA-, JA-, ETH-, and CTK-related DEGs were mainly induced by FS, and a few genes associated with ABA, SA, and IAA were coregulated under the three flooding conditions (Figure 6B and Table S20).We speculate that wheat responds to different types of flooding stress by activating different hormone signaling pathways.

Flooding-Type-Dependent Response Pathways and TFs
Previous research has suggested a possible connection between the number of responsive genes and their association with the complexity and intensity of the imposed stress treatment [49].For example, experiments in barley exposed to combined water deficits and salt stress showed that the duration of individual or combined stresses increased the number of differentially expressed genes [49].Consistent with these studies, the present study revealed that the negative effect of FS on the growth of wheat seedlings was greater than those of HS and WL (Figure 1A).At the same time, we observed clear differences in the numbers and types of DEGs that were upregulated with differences in water depth (Figure 3B,C).Furthermore, previous studies demonstrated that the response of plants to stress is related to the type of stress [49,50].For instance, Luo et al. (2019) found that the CDPK, MAPK, CIPK, and PYL-PP2C-SnRK2 signaling pathways were involved in osmotic stress, while the SOS core pathway was activated by ionic stress [50].Consistent with our study, a KEGG enrichment analysis indicated that several pathways were specially enriched in the WL, HS, and FS treatments, such as the 'brassinosteroid biosynthesis' pathway under the WL treatment, 'taurine and hypotaurine metabolism' under the HS treatment, and 'carbon metabolism' under the FS treatment (Figure 3E-G), suggesting that type-specific response pathways exist in plant responses to flooding stress.
Transcription factors (TFs) are crucial controllers of abiotic stress, including waterlogging, and participate in the regulation of downstream stress-responsive genes [7,9].These proteins have been identified as among the most promising targets for the improvement of plant performance under waterlogging stress [7].Here, a total of 1754 differentially expressed TFs from 35 different families were identified in wheat seedlings under WL, HS, and FS conditions (Figure 6C and Table S21).The greatest difference was the significant over-representation of the bHLH, NAC, MYB, WRKY, and AP2/ERF families in the upregulated DEGs and downregulated DEGs in response to flooding.Consistent with previous studies, these TF families are involved in abiotic stress and may positively improve plant tolerance [7,9].In soybean, Glyma.01G198000,encoding a bHLH TF, was predicted to be a candidate gene involved in seed waterlogging tolerance [51].The constitutive expression of the TaERFVII.1 gene in wheat enhanced the tolerance to waterlogging in transgenic wheat without negative impacts on its development and yield [15].However, in this study, the three copies of the TaERFVII.1 gene (TraesCS5B02G315500, TraesCS5A02G314600, and TraesCS5D02G32080) were barely expressed under both control and flooded conditions (Table S2).This result suggests that wheat tolerance to flooding stress is regulated by other genes.Moreover, a large number of flooding-type specific TFs were also identified across the three flooding conditions; for example, two TF families, ARR-B and YAABY, were only identified under the FS treatment.Stress-type-specific-response TFs were also reported by Osthoff et al. (2019) and Luo et al. (2019), between the NaCl and mannitol treatments [49], and between the water deficit and salt treatments [50], respectively.These results may indicate that plants respond to stress via a complex regulatory mechanism by varying the combination and concentration of TFs according to the stress type.
3.6.Type-Dependent Alterations in DNA Methylation Levels under the WL, HS, and FS Treatments Epigenetic regulators are remarkably diverse in plants, facilitating the phenotypic plasticity of plant development, survival and reproduction in unfavorable environments [52].In plants, alterations in histone modification and DNA methylation are coordinated with changes in the expression of stress-responsive genes to adapt to environmental changes [53].DNA methylation is a well-studied epigenetic mechanism underlying plant stress responses [20,21,54].Furthermore, the response of DNA methylation to stress is related to the intensity and type of stress.For example, chromium stress increased the DNA methylation level in kenaf in a chromium-concentration-dependent manner [21].Li et al. (2022) showed that heat stress increased the DNA methylation level in a heat-sensitive rice group but decreased it in a tolerant rice group [20].The total DNA methylation in sesame increased under drought stress but decreased under waterlogging stress [17].In this study, both WL and HS increased the total DNA methylation level in wheat seedlings, but this was decreased under FS (Table 1).Thus, we deduced that flooding-type-dependent alterations in DNA methylation levels may be a tolerance strategy in wheat.In future studies, it would be valuable to compare the differences in genome-wide methylation under different flooding conditions using high-throughput methylation sequencing.

Plant Material and Flooding Treatments
A commercial spring wheat (Triticum aestivum L.) variety Efumai 1 was used in this study.Seeds of Efumai 1 were surface-sterilized with 2% sodium hypochlorite for 15 min and rinsed with deionized water three times.The seeds were then planted in growth pots (6 plants in each pot), filled with the same amount of soil, and placed in a greenhouse with a 14 h/10 h and 22 • C/18 • C day/night and light/temperature cycle.Relative humidity was maintained at 70%.Six pots were established for each replicate.After 18 days of growth (at the stage with three fully expanded leaves), uniform seedlings were kept and divided into four groups, including the CK and three flooding treatments.The three types of flooding stress treatments were designed based on the depth of the water table, described as waterlogging (WL: keeping the water level at the surface soil level ±1 cm), half submergence (HS: keeping the water level at half of the plant height), and full submergence (FS: keeping the plant fully submerged).For the control, the seedlings were grown in ideal conditions.A completely randomized design with three biological replicates was applied in this experiment.

Physiological Measurements and Growth Parameters
After a 7-day treatment, the maximum photochemical efficiency of PSII (F v /F m ) and relative chlorophyll content (SPAD) were measured in the first completely-extended top leaves using a MINI-PAM-II (Walz, Eiffeltrich, Germany) and a chlorophyll meter (Minolta SPAD-502, Tokyo, Japan), respectively, according to the manufacturer's protocols.The upper, middle, and lower parts of the wheat leaves were measured twice, and the average value was taken as one repeat.Six plants were randomly selected for each biological repeat.Shoot (both stem and leaves) dry weight was measured for each treatment (10 plants of each treatment) after drying at 105 • C in an oven for 72 h.
On the other hand, the first top leaves of Efumai 1 were randomly selected and stored at −80 • C to further detect the physiological and molecular changes that occurred under the CK and three flooding conditions.H 2 O 2 was extracted by homogenizing a 0.3 g sample with 2.7 mL of phosphate buffer (50 mM, pH 6.5) at 4 • C. The H 2 O 2 concentration was determined via a colorimetric method using a commercial kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China), according to the manufacturer's instructions.The absorbance value of each sample was calibrated to a standard concentration curve to calculate the H 2 O 2 concentration.The results are expressed as mmol g −1 protein.The MDA content was extracted using a 5% trichloroacetic acid (TCA) buffer.The activities of enzymes (SOD, POD, and CAT) were extracted using a 50 mM potassium phosphate buffer (pH 7.8, containing 1% (w/v) polyvinylpolypyrrolidone (PVP), and 0.1 mM ethylenediaminetetraacetic acid (EDTA), and 0.3% (w/v) Triton X100), respectively.The absorbance of each sample was measured using a UV-1800 spectrophotometer (Shimadzu, Kyoto, Japan) and quantified according to the method described by Ren et al. (2020) [55].MDA contents and enzyme activities are expressed as nmol g −1 FW and U g −1 FW, respectively.

Methylation-Sensitive Amplified Polymorphism (MSAP) Analysis
The total DNA of each sample was isolated using the CTAB procedure.The procedure of MSAP analysis, which uses two pairs of restriction endonucleases (EcoRI + HpaII and EcoRI + MspI) (Thermo Fisher Scientific, Waltham, MA, USA) for the restrictive digestion of the DNA of each sample, is described in previous research [20].The adaptors and primers are listed in Table S22.The PCR products were separated using a Fragment Analyzer Automated CE System (AATI, Ames, IA, USA) and the DNF-900 dsDNA Reagent Kit 35~5500 bp (AATI, USA), according to Li et al. (2022) [20].The MSAP data were exported using PROSize version 2.0 software (AATI, USA) and transformed into a binary character matrix, using "1" or "0" to indicate the presence or absence of bands.Only the consistent epilocus among the three biological repeats was used for future analysis.The three types of MSAP bands were defined as nonmethylation, hemi-methylation, and full methylation (Table 1) according to Tang et al. (2022) [21].

RNA Extraction, Library Preparation, and Sequencing
The total RNA of each sample was isolated using a TRNzol Universal RNA extraction Kit (Tiangen, Beijing, China) according to the manufacturer's instructions.The quantity and purity of each RNA sample were determined using a Nanodrop 2000 (ThermoFisher Scientific, Waltham, MA, USA) and 2100 Bioanalyzer instrument (Agilent Technologies, Santa Clara, CA, USA), respectively.cDNA Library preparation for RNA-Seq was conducted using a NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (NEB, Ipswich, MA, USA), following the manufacturer's recommendations.The final library quality was assessed using the Agilent Bioanalyzer 2100 system (Agilent Technologies, USA), and the library was sequenced on an Illumina NovaSeq platform (150 bp pair end).Three biological replicates were used for RNA-Seq experiments.

RNA-Seq Read Mapping, Sequence Assembly, and Differential Expression
After filtering the sequencing adapters and low-quality reads, valid data were aligned to the wheat reference genome (https://urgi.versailles.inra.fr/download/iwgsc/IWGSC_RefSeq_Assemblies/v1.0/, accessed on 1 May 2023) using HISAT2, and the alignments were sorted using SAMTools v1.8 with default parameters.The uniquely mapped reads were aligned with high-confidence genome annotation files (IWGSC RefSeq v1.1, https: //wheat-urgi.versailles.inra.fr/Seq-Repository/Annotations,accessed on 5 May 2023).The reads mapped to each gene were counted using featureCounts.The expression level of each gene was estimated from the fragments per kilobase of transcript per million fragments mapped (FPKM) value.Differential expression analysis between the CK and the treatments was performed using DESeq2 software (version 1.26.0).An adjusted p value < 0.05 and |log2-fold fold change| ≥1 were used as criteria for identifying differentially expressed genes (DEGs).

Functional Analysis of DEGs
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of the DEGs were performed using the free online platform, OmicShare tools (http://www.omicshare.com/tools,accessed on 10 May 2023), where the threshold was a corrected p value < 0.05.Transcription factors (TFs) were predicted and classified into different families using PlantTFDB (http://planttfdb.gao-lab.org/,accessed on 11 May 2023).Heatmaps of gene expression were generated using TBtools (version v2.008) [56].

Quantitative Real-Time PCR (qRT-PCR) Validation
To validate the repeatability and reproducibility of the gene expression data obtained via RNA-Seq in the four conditions, a total of 20 genes were randomly selected, and gene-specific primers were designed with the online tool Primer 3 (http://primer3.ut.ee, accessed on 25 May 2023).First-strand cDNA was synthesized using the UEIris RT mix with a Dnase (All-in-One) kit (US Everbright, Jiangsu, China) according to the manufacturer's protocol.qRT-PCR was performed with the QuantStudio™ 6 Flex Real-Time PCR System (Applied Biosystems, Foster, CA, USA) using Biolife Green I (US Everbright, Jiangsu, China).Actin was used as an internal control to normalize the gene expression level.The relative expression levels of selected genes were determined using the 2 −∆∆CT method.The primers for qRT-PCR are presented in Table S23.All reactions were performed in triplicate.

Statistical Analysis
All statistical analyses were performed using SPSS 18.0.The data were analyzed using a one-way analysis of variance (ANOVA, Duncan's test).Differences with p < 0.05 were considered significant.

Conclusions
Our research provides the first characterization of the physiological, epigenetic, and transcriptomic responses of wheat seedlings to different levels of flooding stress.The results suggest that there was no difference in the effects of WL and HS on biomass and physiology without the submergence of functional leaves.However, significant differences in epigenetic and transcriptional responses were identified among the three flooding conditions.FS increased the levels of H 2 O 2 and MDA in the leaves and led to greater reductions in the chlorophyll content (SPAD), photochemical efficiency (F v /F m ), and biomass production than HS and WL.Type-dependent alterations in DNA methylation were observed between the three flooding conditions.Similarly, the common and specific response pathways of wheat to the flooding stresses of different water depths were also observed.We found that the DEGs involved in 'photosynthesis', 'phenylpropanoid biosynthesis', and 'plant hormone signal transduction' were closely related to the flooding response, and these genes were involved in specific responses to different flooding stress types.Thus, our results provide a comprehensive view of the complex molecular events involved in the responses of wheat leaves to flooding stress, which will promote research on the development of flood-resistant crops and provide new information for wheat breeders.In the future, a systemic investigation of epigenetic regulation and flooding tolerance would be valuable.Moreover, flood depth is also a factor worth considering in future flood assessments.

Figure 1 .
Figure 1.Biomass and photosynthesis analysis of Efumai 1 seedlings in the control and three flooding treatments.(A) Images of Efumai 1 plants after a 7-day treatment under CK, WL, HS, and FS conditions; (B) shoot dry weight; (C) SPAD value; and (D) Fv/Fm value.Significant differences were assessed using ANOVA with p < 0.05; The different letters above the bars indicated significant differences (p < 0.05).Three independent biological replicates for each treatment were performed and six pots were established for each replicate.

Figure 2 .
Figure 2. Physiological responses of wheat seedlings under the control and three flooding treatments.(A) H2O2 content; (B) MDA content; (C) SOD activity; (D) CAT activity; and (E) POD activity.Values are expressed as mean ± SD.The different letters above the bars indicated significant differences using ANOVA with p < 0.05.

Figure 2 .
Figure 2. Physiological responses of wheat seedlings under the control and three flooding treatments.(A) H 2 O 2 content; (B) MDA content; (C) SOD activity; (D) CAT activity; and (E) POD activity.Values are expressed as mean ± SD.The different letters above the bars indicated significant differences using ANOVA with p < 0.05.

Figure 3 .
Figure 3. Overview of the effects of flooding treatments on gene expression.(A) Number of DEGs in the three comparisons.(B) Overlapping upregulated DEGs in the three comparisons.(C) Overlapping downregulated DEGs in the three comparisons.(D) Correlation analysis of differentially expressed genes between RT-PCR analysis (2 −ΔΔCt value) and RNA-seq experiment (Log2 fold change value).The top 25 significantly enriched KEGG pathways of DEGs induced by WL (E), HS (F), and FS (G) treatments, from the outside to the inside.The first circle represents the top 25 enrichment pathways, and the number outside the circle is the coordinate ruler of the number of genes.The second circle represents the number and −log10 (p value) of background genes in this pathway.The third circle represents the number of upregulated and downregulated DEGs in this pathway.The fourth circle represents the value of the Rich Factor in each pathway.

Figure 3 .
Figure 3. Overview of the effects of flooding treatments on gene expression.(A) Number of DEGs in the three comparisons.(B) Overlapping upregulated DEGs in the three comparisons.(C) Overlapping downregulated DEGs in the three comparisons.(D) Correlation analysis of differentially expressed genes between RT-PCR analysis (2 −∆∆Ct value) and RNA-seq experiment (Log 2 fold change value).The top 25 significantly enriched KEGG pathways of DEGs induced by WL (E), HS (F), and FS (G) treatments, from the outside to the inside.The first circle represents the top 25 enrichment pathways, and the number outside the circle is the coordinate ruler of the number of genes.The second circle represents the number and −log10 (p value) of background genes in this pathway.The third circle represents the number of upregulated and downregulated DEGs in this pathway.The fourth circle represents the value of the Rich Factor in each pathway.

Figure 4 .
Figure 4. DEGs related to the photosynthesis pathways.The KEGG pathway map of 'photosynthesis' (A) and 'photosynthesis-antenna protein' (B).The DEG expression pattern was used to annotate and generate a corresponding map.The green box with gene symbols denotes downregulated expression in the chlorosis group, while the red box denotes upregulated expression.The genes that were not significantly altered are displayed with a white box.(C-G) Expression profiles (Log2 fold change) of genes related to 'photosynthesis' pathway.(H,I) Expression profiles (Log2 fold change) of genes related to 'photosynthesis-antenna proteins' pathway.

Figure 4 .
Figure 4. DEGs related to the photosynthesis pathways.The KEGG pathway map of 'photosynthesis' (A) and 'photosynthesis-antenna protein' (B).The DEG expression pattern was used to annotate and generate a corresponding map.The green box with gene symbols denotes downregulated expression in the chlorosis group, while the red box denotes upregulated expression.The genes that were not significantly altered are displayed with a white box.(C-G) Expression profiles (Log2 fold change) of genes related to 'photosynthesis' pathway.(H,I) Expression profiles (Log2 fold change) of genes related to 'photosynthesis-antenna proteins' pathway.

Figure 5 .
Figure 5. DEGs involved in the phenylpropanoid biosynthesis pathway.(A) The location of DEGs in phenylpropanoid biosynthesis pathway.(B) The expression pattern of the phenylpropanoidbiosynthesis-related genes.

Figure 5 .
Figure 5. DEGs involved in the phenylpropanoid biosynthesis pathway.(A) The location of DEGs in phenylpropanoid biosynthesis pathway.(B) The expression pattern of the phenylpropanoidbiosynthesis-related genes.

Figure 6 .
Figure 6.DEGs (A,B) and transcription factors (C) involved in the plant hormone signal transduction pathway.Figure 6. DEGs (A,B) and transcription factors (C) involved in the plant hormone signal transduction pathway.

Figure 6 .
Figure 6.DEGs (A,B) and transcription factors (C) involved in the plant hormone signal transduction pathway.Figure 6. DEGs (A,B) and transcription factors (C) involved in the plant hormone signal transduction pathway.

Author
Contributions: Y.X. and C.J. conceived the project, designed the experiments, and revised the manuscript.B.L. and W.H. performed the experiments, analyzed the data, and drafted the manuscript.S.Z., L.X. and C.Y. performed the experiments.Z.Z., Y.G. and M.Z.analyzed the data and revised the manuscript.All authors have read and agreed to the published version of the manuscript.

Table 1 .
MSAP-based cytosine methylation levels in wheat leaves under CK, WL, HS, and FS treatments.

Table 2 . Alternations in DNA methylation patterns induced by the WL, HS, and FS treatments. Description of Groups CK Flooding-Treated WL HS FS HpaII MspI HpaII MspI
Analysis of DEGs in Wheat in Response to the Three Types of Flooding Stress