Transcriptomic, Physiological, and Metabolomic Response of an Alpine Plant, Rhododendron delavayi, to Waterlogging Stress and Post-Waterlogging Recovery

Climate change has resulted in frequent heavy and prolonged rainfall events that exacerbate waterlogging stress, leading to the death of certain alpine Rhododendron trees. To shed light on the physiological and molecular mechanisms behind waterlogging stress in woody Rhododendron trees, we conducted a study of Rhododendron delavayi, a well-known alpine flower species. Specifically, we investigated the physiological and molecular changes that occurred in leaves of R. delavayi subjected to 30 days of waterlogging stress (WS30d), as well as subsequent post-waterlogging recovery period of 10 days (WS30d-R10d). Our findings reveal that waterlogging stress causes a significant reduction in CO2 assimilation rate, stomatal conductance, transpiration rate, and maximum photochemical efficiency of PSII (Fv/Fm) in the WS30d leaves, by 91.2%, 95.3%, 93.3%, and 8.4%, respectively, when compared to the control leaves. Furthermore, the chlorophyll a and total chlorophyll content in the WS30d leaves decreased by 13.5% and 16.6%, respectively. Both WS30d and WS30d-R10d leaves exhibited excessive H2O2 accumulation, with a corresponding decrease in lignin content in the WS30d-R10d leaves. At the molecular level, purine metabolism, glutathione metabolism, photosynthesis, and photosynthesis-antenna protein pathways were found to be primarily involved in WS30d leaves, whereas phenylpropanoid biosynthesis, fatty acid metabolism, fatty acid biosynthesis, fatty acid elongation, and cutin, suberin, and wax biosynthesis pathways were significantly enriched in WS30d-R10d leaves. Additionally, both WS30d and WS30d-R10d leaves displayed a build-up of sugars. Overall, our integrated transcriptomic, physiological, and metabolomic analysis demonstrated that R. delavayi is susceptible to waterlogging stress, which causes irreversible detrimental effects on both its physiological and molecular aspects, hence compromising the tree’s ability to fully recover, even under normal growth conditions.


Introduction
Climate change has had a profound impact on the global hydrological cycle, resulting in a marked increase in heavy or recurrent precipitation events in certain regions. Such rainfall frequently exceeds the soil's water-holding capacity, resulting in oversaturation of terrestrial ecosystems, especially in zones with poor drainage, clay-rich soils, and lowlying areas. This oversaturation can trigger severe waterlogging stress as a result of soil ponding [1][2][3].
Waterlogging stress impedes the gas exchange between plant tissues and soil, causing anoxic or hypoxic conditions in the waterlogged tissues. These conditions can cause an "energy crisis" due to the inhibition of ATP generation [1,4]. Plant responses to waterlogging stress involve complex signal transduction pathways that alter the synthesis and balance of plant hormones, resulting in stomatal closure and morphological changes [5,6]. Additionally, waterlogging stress induces the accumulation of plant metabolites such as flavonoids, vitamin B6, and sugars that play a crucial role in the response to waterlogging stress [7]. Previous studies have shown that hormone mediated transcription factors, including MYBs, ERFs, WRKYs, and NACs, regulate plant responses to waterlogging stress [8,9]. Moreover, waterlogging triggers anaerobic metabolism in roots, leading to the production of ethanol, acetaldehyde, and other harmful substances that can poison root cells and damage protein structure [10,11]. Ultimately, this leads to accelerated root decay, damage, and even death [12].
Currently, extensive research has been conducted on the molecular and physiological responses to waterlogging stress in model plants such as Arabidopsis and crops such as Solanum lycopersicum, Cucumis sativus, and Glycine max [2,13]. However, trees play a crucial ecological role in global climate change and are essential for urban environments, such as greening streets and parks. Therefore, it is crucial to understand the mechanisms of underlying tree response to waterlogging stress. Although some studies have reported on the physiological response mechanisms of waterlogging-sensitive species or varieties [3,[14][15][16], only a few have investigated changes in the transcriptome and metabolome of woody economic species, such as poplar, cotton, oak, kiwifruit, and Guazuma ulmifolia, exposed to waterlogging stress [17][18][19][20]. Thus, there is a need to investigate the response mechanisms of specific woody trees to waterlogging stress, as given that their tolerance mainly depends on the species and degree of waterlogging [21].
Rhododendron, the largest genus in the Ericaceae family and one of the most notable alpine flowers, comprises approximately 1143 species of evergreen and deciduous woody plants. In the northwest Guizhou province of China, wild Rhododendron trees in the Baili Azalea Nature Reserve (BANR) represent significant potential tourism resources for local communities. Although these trees typically grow in high mountain areas, some low-lying areas or clay-rich soils have experienced long-term ponding in recent years, hindering the growth of Rhododendron species and leading to the death of some trees. Despite the importance of Rhododenron trees in this region, there is currently a scarcity of research on the physiological and molecular response to waterlogging stress and post-waterlogging recovery in this genus.
Rhododendron delavayi, a crucial member of the Rhododendron genus in the BANR, has unfortunately experienced mortality in low-lying areas or clay-rich soils due to continuous rainfall. This has prompted us to urgently investigate the response mechanisms of R. delavayi to waterlogging stress. In this study, R. delavayi, an alpine plant, has been utilized as a model system to explore the physiological and molecular responses to waterlogging stress and post-waterlogging recovery. This approach aims to provide key insights into how alpine plants respond to waterlogging stress both at physiological and molecular levels. Ultimately, our findings can contribute a fundamental understanding of the death of R. delavayi caused by waterlogging stress.

Waterlogging Stress Caused the Aged Leaf Wilting in R. delavayi Seedlings
To explore the physiological and molecular responses of R. delavayi to waterlogging stress, R. delavayi seedlings were subjected to various durations of waterlogging stress (WS10d, WS20d, and WS30d) followed by a post-waterlogging recovery period (WS10d-R10d, WS20d-R10d, and WS30d-R10d). During WS10d, WS20d, WS10d-R10d, and WS20d-R10d, no stress symptoms were observed when compared to their pre-waterlogging state ( Figure 1). However, prolonged waterlogging stress (WS30d) led to the wilting of aged leaves in some seedlings (as indicated by white arrow in Figure 1) and even 10 days of post-waterlogging recovery was insufficient to fully restore the wilted leaves to their prewaterlogging state (as indicated by red arrow in Figure 1). These results suggested that waterlogging stress lasting for 30 days can severely damage R. delavayi seedlings.
WS20d-R10d, no stress symptoms were observed when compared to their pre-waterlogging state ( Figure 1). However, prolonged waterlogging stress (WS30d) led to the wilting of aged leaves in some seedlings (as indicated by white arrow in Figure 1) and even 10 days of post-waterlogging recovery was insufficient to fully restore the wilted leaves to their pre-waterlogging state (as indicated by red arrow in Figure 1). These results suggested that waterlogging stress lasting for 30 days can severely damage R. delavayi seedlings. Figure 1. The leaf phenotype of R. delavayi after waterlogging stress and post-waterlogging recovery. Control (CK), waterlogging stress for 10 days (WS10d), 20 days (WS20d), and 30 days (WS30d) were indicated. Post-waterlogging recovery for 10 days was indicated by WS10d-R10d, WS20d-R10d, WS30d-R10d. Wilting of leaves after WS30d and WS30d-R10d was represented by white and red arrow, respectively. Bar = 5.8 cm.

Transcriptome Assembly, Analysis of Differentially Expressed Genes (DEGs), and qRT-PCR Validation
To investigate the molecular mechanisms underlying the response of R. delavayi to waterlogging stress at the transcriptome level, RNA sequencing was performed on leaves from CK, WS30d, and WS30d-R10d plants, based on the stress symptoms of R. delavayi leaves caused by waterlogging stress. Over 20,163,917 clean reads were obtained from each sample, with the lowest mapped ratio being 76.95% in WS30d-R10d-2. The guanine and cytosine/total base (GC content) and the quality of base calling accuracy at 99.9% (Q30) were more than 47.42% and 92.33%, respectively (Table S1). After assembly of the high-quality reads, a total of 97,152 unigenes and 402,244 trancripts with mean lengths of 870 bp and 1705 bp, respectively, were obtained (Table S2).
To determine the variation and similarity of the gene expression profiles among the samples, we conducted principal component analysis (PCA) of all detected genes using normalized Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values. The difference in gene expression among the three samples was statistically significant (p = 0.032), with PC1 and PC2 accounting for 58.5% and 26.4% of the total variation, respectively ( Figure 2A). The PCA plot showed that the data from the three biological replicates were closely clustered together and separated between CK and WS30d and Control (CK), waterlogging stress for 10 days (WS10d), 20 days (WS20d), and 30 days (WS30d) were indicated. Post-waterlogging recovery for 10 days was indicated by WS10d-R10d, WS20d-R10d, WS30d-R10d. Wilting of leaves after WS30d and WS30d-R10d was represented by white and red arrow, respectively. Bar = 5.8 cm.

Transcriptome Assembly, Analysis of Differentially Expressed Genes (DEGs), and qRT-PCR Validation
To investigate the molecular mechanisms underlying the response of R. delavayi to waterlogging stress at the transcriptome level, RNA sequencing was performed on leaves from CK, WS30d, and WS30d-R10d plants, based on the stress symptoms of R. delavayi leaves caused by waterlogging stress. Over 20,163,917 clean reads were obtained from each sample, with the lowest mapped ratio being 76.95% in WS30d-R10d-2. The guanine and cytosine/total base (GC content) and the quality of base calling accuracy at 99.9% (Q30) were more than 47.42% and 92.33%, respectively (Table S1). After assembly of the high-quality reads, a total of 97,152 unigenes and 402,244 trancripts with mean lengths of 870 bp and 1705 bp, respectively, were obtained (Table S2).
To determine the variation and similarity of the gene expression profiles among the samples, we conducted principal component analysis (PCA) of all detected genes using normalized Fragments Per Kilobase of transcript per Million mapped reads (FPKM) values. The difference in gene expression among the three samples was statistically significant (p = 0.032), with PC1 and PC2 accounting for 58.5% and 26.4% of the total variation, respectively ( Figure 2A). The PCA plot showed that the data from the three biological replicates were closely clustered together and separated between CK and WS30d and WS30d-R10d (Figure 2A), demonstrating that the sequencing data from stressed leaves differed significantly from those of the CK leaves.
WS30d-R10d (Figure 2A), demonstrating that the sequencing data from stressed leaves differed significantly from those of the CK leaves. We subsequently identified differentially expressed genes (DEGs) among the WS30d, WS30d-R10d, and CK samples. Comparion of WS30d leaves to the CK leaves revealed 2521 up-regulated and 2774 down-regulated DEGs, while comparison of WS30d-R10d to CK leaves identified 1589 up-regulated and 1840 down-regulated DEGs ( Figure 2B, Tables S3 and S4). DEGs between WS30d and WS30d-R10d were also obtained (Table S5). These results suggested that waterlogging stress and post-waterlogging recovery significantly induced gene transcription in R. delavayi leaves.
To verify the reliability of DEGs identified by RNA sequencing, we randomly selected six genes and quantified their expression levels using quantitative Real-Time PCR (qRT-PCR) analysis. The results demonstrated the consistency between the qRT-PCR data and the RNA sequencing data, thereby confirming the reliability of the DEG screened by RNA sequencing ( Figure S1).

KEGG Pathways by DEGs of Waterlogging Stress Versus Control
To further comprehend the biological functions of the 5295 DEGs induced by 30 days of waterlogging stress, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The up-regulated DEGs were significantly enriched in the "purine metabolism" and "glutathione metabolism" pathways (p < 0.05) ( Figure 3A), while the downregulated DEGs were significantly enriched in "photosynthesis", "photosynthesis-antenna proteins", "fat acid biosynthesis", and "glycosaminoglycan degradation" pathways (p < 0.05) ( Figure 3B). We subsequently identified differentially expressed genes (DEGs) among the WS30d, WS30d-R10d, and CK samples. Comparion of WS30d leaves to the CK leaves revealed 2521 up-regulated and 2774 down-regulated DEGs, while comparison of WS30d-R10d to CK leaves identified 1589 up-regulated and 1840 down-regulated DEGs ( Figure 2B, Tables S3 and S4). DEGs between WS30d and WS30d-R10d were also obtained (Table S5). These results suggested that waterlogging stress and post-waterlogging recovery significantly induced gene transcription in R. delavayi leaves.
To verify the reliability of DEGs identified by RNA sequencing, we randomly selected six genes and quantified their expression levels using quantitative Real-Time PCR (qRT-PCR) analysis. The results demonstrated the consistency between the qRT-PCR data and the RNA sequencing data, thereby confirming the reliability of the DEG screened by RNA sequencing ( Figure S1).

KEGG Pathways by DEGs of Waterlogging Stress Versus Control
To further comprehend the biological functions of the 5295 DEGs induced by 30 days of waterlogging stress, we performed Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. The up-regulated DEGs were significantly enriched in the "purine metabolism" and "glutathione metabolism" pathways (p < 0.05) ( Figure 3A), while the down-regulated DEGs were significantly enriched in "photosynthesis", "photosynthesisantenna proteins", "fat acid biosynthesis", and "glycosaminoglycan degradation" pathways (p < 0.05) ( Figure 3B).
In the "purine metabolism" pathway, the enzymes encoded by the up-regulated DEGs directly or indirectly catalyzed ATP production (Table 1), indicating their crucial role in providing energy. In the "Glutathione metabolism" pathway, 21 up-regulated DEGs were annotated as glutamate-cysteine ligase, glutathione peroxidase, glutathione reductase, and isocitrate dehydrogenase, which participate in the glutathione cycle (Table 1) and facilitate the accumulation of glutathione, as evidenced by metabolome detection (Table 2). Furthermore, four DEGs were up-regulated and annotated as L-ascorbate peroxidase and ribonucleoside-diphosphate reductase (Table 1), suggesting their involvement in the ascorbate cycle.  In the "purine metabolism" pathway, the enzymes encoded by the up-regulated DEGs directly or indirectly catalyzed ATP production (Table 1), indicating their crucia role in providing energy. In the "Glutathione metabolism" pathway, 21 up-regulated DEGs were annotated as glutamate-cysteine ligase, glutathione peroxidase, glutathione reductase, and isocitrate dehydrogenase, which participate in the glutathione cycle (Table  1) and facilitate the accumulation of glutathione, as evidenced by metabolome detection (Table 2). Furthermore, four DEGs were up-regulated and annotated as L-ascorbate pe roxidase and ribonucleoside-diphosphate reductase ( Table 1), suggesting their involve ment in the ascorbate cycle.    In the "photosynthesis" pathway, the down-regulated DEGs were annotated to encode proteins such as cytochrome b6/f complex, electron transport chain components, ATPase, and subunits of photosystem I (PSI) or II (PSII) (as shown in Table 1). Further analysis revealed a significant reduction in CO 2 assimilation rate, stomatal conductance, transpiration rate, and maximum photochemical efficiency by 91.2%, 95.3%, 93.3%, and 8.4%, respectively, in the WS30d leave in comparison with CK ( Figure 4). Additionally, the chlorophyll a and total chlorophyll content in the WS30d leaves decreased by 13.5% and 16.6%, respectively (Figure 4). These observations suggested that excessive light energy might lead to the accumulation of reactive oxygen species (ROS) in the chloroplast. This was further supported by the significant increase in hydrogen peroxide content in WS30d leaves, which was found to be 3.39 times higher than the control (Figure 4). 9 of 21 transpiration rate, and maximum photochemical efficiency by 91.2%, 95.3%, 93.3%, and 8.4%, respectively, in the WS30d leave in comparison with CK ( Figure 4). Additionally, the chlorophyll a and total chlorophyll content in the WS30d leaves decreased by 13.5% and 16.6%, respectively (Figure 4). These observations suggested that excessive light energy might lead to the accumulation of reactive oxygen species (ROS) in the chloroplast. This was further supported by the significant increase in hydrogen peroxide content in WS30d leaves, which was found to be 3.39 times higher than the control (Figure 4). In the "fatty acid biosynthesis" pathway, down-regulated DEGs were annotated as FabF, FabH, FabI, fatty acyl-ACP thioesterase, and long-chain acyl-CoA synthetase (Table  1). These genes are involved in synthesizing of long-chain acyl-CoA. Similarly, in the "glycosaminoglycan degradation" pathway, down-regulated DEGs were annotated as heparanase and hexosaminidase (Table 1). In the "fatty acid biosynthesis" pathway, down-regulated DEGs were annotated as FabF, FabH, FabI, fatty acyl-ACP thioesterase, and long-chain acyl-CoA synthetase ( Table 1). These genes are involved in synthesizing of long-chain acyl-CoA. Similarly, in the "glycosaminoglycan degradation" pathway, down-regulated DEGs were annotated as heparanase and hexosaminidase (Table 1).

Major Transcription Factors Families Active during Waterlogging Waterlogging Recovery
Transcription factors (TFs) play a crucial role in the transcri that occurs in response to abiotic stress. To identify the enriched to waterlogging stress, we screened the DEGs in WS30d leaves. Ou In the "fatty acid metabolism" and "fatty acid biosynthesis" pathways, the proteins encoded by DEGs were consistent with those in the "fatty acid biosynthesis" pathway mentioned above (Table S6). The down-regulated DEGs may impede the synthesis of longchain acyl-CaA. In the "fatty acid elongation" pathway, DEGs were annotated as 3-ketoacyl-CoA synthase, very-long-chain (3R)-3-hydroxyacyl-CoA dehydratase, and very-long-chain 3-oxoacyl-CoA reductase (Table S6). The decreased expression of these DEGs suggested that the synthesis of long-chain fatty acids with long-chain acyl-CaA as a precursor may be inhibited. Long-chain fatty acids are utilized for the biosynthesis of cutin, suberin, and wax. The down-regulated DEGs involved in these pathways may impede the biosynthesis of long-chain esters (Table S6).

Major Transcription Factors Families Active during Waterlogging Stress and Post-Waterlogging Recovery
Transcription factors (TFs) play a crucial role in the transcriptional reprogramming that occurs in response to abiotic stress. To identify the enriched TF families in response to waterlogging stress, we screened the DEGs in WS30d leaves. Our analysis revealed that 20 TF families were encoded by 118 up-regulated DEGs and 117 down-regulated DEGs (Figure 7). Among the up-regulated DEGs, the WRKY family was the most dominant, while the bHLH family was the most dominant among the down-regulated DEGs. Interestingly, the DEGs encoding most TF families, such as MYB and AP2/ERF-ERF, were both up-regulated and down-regulated under waterlogging stress (Figure 7). In the common DEGs between WS30d and WS30d-R10d, the MYB-related and bHLH families were the most dominant TF families encoded by the up-regulated and down-regulated DEGs, respectively ( Figure S2).

Metabolites Accumulations in Response to Waterlogging Stress and Post-Waterlogging Recovery
To investigate the accumulation of metabolites in response to waterlogging stress an post-waterlogging recovery, we utilized the GC-MS platform to determine the metabol levels in CK, WS30d, and WS30d-R10d leaves. To summarize the similarities and diffe ences between these groups, we employed principal component analysis (PCA) and o thogonal projections to latent structures-discriminant analysis (OPLS-DA). Our PCA an ysis revealed that the data from the three duplicate samples were closely clustered, wi PC1 and PC2 accounting for 20.6% and 16.7% of the total variation, respectively (Figu S3A). In the OPLS-DA score plot, component 1 and component 2 accounted for 19.8% an 12.2% of the total variation, respectively ( Figure S3B). We identified a total of 149 putati metabolites and screened for differential metabolites between treatment groups (VIP > and p < 0.05). Compared to CK, WS30d leaves showed significant upregulation in 22 m tabolites, including seven sugars (glucose, sedoheptulose, galactose, sucrose, lyxo galactonic acid, N-Acetyl-beta-D-mannosamine), three organic acids, three lipids, two cohols, two flavonoids, and three others ( Table 2). WS30d-R10d leaves showed upregu tion in six sugars, three lipids, two alcohols and flavonoids, and five others (Table S Interestingly, the common metabolites in WS30d and WS30d-R10d leaves revealed th soluble sugars remained highly accumulated during post-waterlogging recovery (Tab S8).

Waterlogging Stress Inhibited Photosynthesis in R. delavayi Leaves
Gas exchange is usually used to assess the tolerance of plant species to waterloggi [21]. Tolerant species show a slight reduction in their photosynthetic rate [3], where sensitive species experience a significant reduction [22,23]. In the case of R. delavayi, t CO2 assimilation rate in WS30d leaves decreased by 91.2%, suggesting that it is a wate logging-sensitive species (Figure 4A), similar to Quercus petraea and Persea America [23,24]. This decrease in CO2 assimilation rate in WS30d leaves was accompanied by simultaneous reduction in stomatal conductance and transpiration rate ( Figure 4B,C), i dicating that photosynthesis may be affected by stomatal limitation [25]. Although t precise mechanism behind the decrease in stomatal conductance caused by waterloggi

Metabolites Accumulations in Response to Waterlogging Stress and Post-Waterlogging Recovery
To investigate the accumulation of metabolites in response to waterlogging stress and post-waterlogging recovery, we utilized the GC-MS platform to determine the metabolite levels in CK, WS30d, and WS30d-R10d leaves. To summarize the similarities and differences between these groups, we employed principal component analysis (PCA) and orthogonal projections to latent structures-discriminant analysis (OPLS-DA). Our PCA analysis revealed that the data from the three duplicate samples were closely clustered, with PC1 and PC2 accounting for 20.6% and 16.7% of the total variation, respectively ( Figure S3A). In the OPLS-DA score plot, component 1 and component 2 accounted for 19.8% and 12.2% of the total variation, respectively ( Figure S3B). We identified a total of 149 putative metabolites and screened for differential metabolites between treatment groups (VIP > 1, and p < 0.05). Compared to CK, WS30d leaves showed significant upregulation in 22 metabolites, including seven sugars (glucose, sedoheptulose, galactose, sucrose, lyxose, galactonic acid, N-Acetyl-beta-D-mannosamine), three organic acids, three lipids, two alcohols, two flavonoids, and three others ( Table 2). WS30d-R10d leaves showed upregulation in six sugars, three lipids, two alcohols and flavonoids, and five others (Table S7). Interestingly, the common metabolites in WS30d and WS30d-R10d leaves revealed that soluble sugars remained highly accumulated during post-waterlogging recovery (Table S8).

Waterlogging Stress Inhibited Photosynthesis in R. delavayi Leaves
Gas exchange is usually used to assess the tolerance of plant species to waterlogging [21]. Tolerant species show a slight reduction in their photosynthetic rate [3], whereas sensitive species experience a significant reduction [22,23]. In the case of R. delavayi, the CO 2 assimilation rate in WS30d leaves decreased by 91.2%, suggesting that it is a waterloggingsensitive species (Figure 4A), similar to Quercus petraea and Persea Americana [23,24]. This decrease in CO 2 assimilation rate in WS30d leaves was accompanied by a simultaneous reduction in stomatal conductance and transpiration rate ( Figure 4B,C), indicating that photosynthesis may be affected by stomatal limitation [25]. Although the precise mechanism behind the decrease in stomatal conductance caused by waterlogging stress remains unclear, it has been widely reported. In addition, the decrease in maximum photochemical efficiency of PSII (Fv/Fm) in WS30d leaves ( Figure 4D) suggested that photoinhibition had occurred in the stressed leaves [26,27].
The suppression of photosynthesis could also be affected by both chlorophyll content and photosynthetic enzyme activity. In the case of WS30d leaves, there was a decrease of 13.5% and 16.6% in chlorophyll a and total chlorophyll content, respectively (Figure 4), which is consistent with previous studies [3,27]. Additionally, the decrease in chlorophyll a and chlorophyll (a + b) content in stressed leaves indicates that the inhibition of photosynthesis, which results in more damage to PSII, is also more limiting [3]. This, in turn, may lead to the production of excessive ROS. RNA-sequencing analysis of stressed leaves revealed that down-regulated DEGs were enriched in "photosynthesis" and "photosynthesis-antenna proteins" pathways ( Figure 3B). The DEGs annotated in these pathways were found to be involved in PSI, PSII, cytochrome b6-f complex, and ATP synthase ( Table 1), suggesting that the activity of proteins in the photosynthesis pathway was inhibited by waterlogging stress. Although the chlorophyll a and chlorophyll (a + b) content in WS30d-R10d leaves could be restored to the control level, the CO 2 assimilation rate was still reduced (Figure 4). These findings are similar to previous studies on cotton, a woody plant that is sensitive to waterlogging stress [28,29].
Transcription factors play a crucial role in regulating gene expression, particularly in the synthesis of plant secondary metabolites such as anthocyanins, flavonoids, and lignin. MYB and bHLH transcription factors, in particular, have been identified as key regulators of this process [30,31]. In waterlogged leaves, it is posited that the downregulation of MYB and bHLH genes may exert negative regulation on the expression of genes involved in flavonoid biosynthesis, leading to flavonoid accumulation, or alternatively, positively regulate genes involved in lignin biosynthesis, thereby inhibiting lignin synthesis ( Figure 6). Further, WRKY and MYB-related transcription factors are instrumental in regulating ABA biosynthesis [30,32], a process that has been shown to play a significant role in the waterlogged plants [33]. It is hypothesized that the up-regulation of WRKY-and MYB-related transcription factors in waterlogged leaves may regulate ABA biosynthesis, ultimately causing leaf stomatal closure and reduction in stomatal conductance (Figure 4). Additionally, plant hormones induce the expression of ethylene response factors, such as RAP2.3, improving expression levels of downstream genes and enhancing waterlogging or post-waterlogging resistance [32,34]. Based on our findings, it is postulated that the up-regulation of ERF transcription factors may stimulate H 2 O 2 accumulation in the stressed leaves by ethylene biosynthesis (Figure 4). Notably, the observed alterations in the expression levels of both up-regulated and down-regulated ERF family members in conjunction with the varying expression patterns observed in other transcription factors and plant hormone signaling pathways indicates a highly intricate regulatory network governing the plant response to waterlogging stress [17].

Waterlogging Stress Induced Oxidative Stress in R. delavayi Seedlings
The existing body of evidence indicates that short-term waterlogging stress has the potential to generate ROS in plant roots. In leaf tissues, several abiotic stressors can result in a reduction in CO 2 assimilation rate, resulting in the accumulation of excess H 2 O 2 and consequent oxidative stress [35]. In addition, waterlogging stress appears to alter the accumulation of ethylene and provoke the production of H 2 O 2 , as evidenced by the increased H 2 O 2 accumulation in the WS30d leaves ( Figure 4H). Such H 2 O 2 accumulation may disrupt normal metabolic processes by promoting lipid peroxidation, compromising membrane integrity, as well as protein and DNA oxidation [36][37][38]. Analysis of gene expression patterns revealed that the up-regulated DEGs in the stressed leaves were significantly enriched in the "purine metabolism" pathways ( Figure 3A), suggesting that the plants may initiate gene expression mechanisms to repair DNA damage caused by H 2 O 2 accumulation (Table 1).
In the context of ROS accumulation under waterlogging stress, the ascorbate-glutathione cycle plays a critical role in enabling plants to cope with oxidative stress to a certain extent [39]. Notably, the observed increases in ascorbate or glutathione concentrations in the WS30d leaves or WS30d-R10d (Tabels 2 and S9) were found to be consistent with the up-regulated DEGs in the glutathione metabolism pathway ( Figure 3A). This strongly suggests that the accumulation of glutathione or ascorbate could serve as non-enzymatic antioxidants to preserve the balance of H 2 O 2 . However, despite the fact that waterlogged seedlings were able to recover, their H 2 O 2 content remained at a higher level than the post-waterlogging recovery ( Figure 4H). As such, oxidative stress was sustained and even accelerated the damage inflicted on the plants during the recovery, as shown in Figure 1.

Lignin and Cuticle Biosynthesis was Continuously Inhibited during Post-Waterlogging Recovery
Lignin, a phenolic biopolymer, is primarily synthesized through the phenylpropanoid biosynthesis pathway and is typically concentrated in the secondary cell wall of vascular plants. Lignin plays a vital role in providing mechanical support to plant tissues, as well as facilitating the transportation of nutrients or carbohydrates [40]. Notably, the accumulation of lignin is known to enhance cell wall reinforcement, thereby improving resistance to waterlogging stress [41]. Previous studies have shown that waterlogging stress can reduce lignin content by inhibiting the expression of genes involved in lignin biosynthesis [42]. Our transcriptome data analysis revealed that the DEGs in WS30d-R10d leaves were significantly enriched in the phenylpropanoid biosynthesis pathway ( Figure 5D), with their expression being noticeably inhibited (Table S6). These results were strongly aligned with the observed decrease in lignin content in WS30d-R10d leaves ( Figure 6). As such, we hypothesize that waterlogging may reduce the biosynthesis of lignin in the stressed leaves, ultimately leading to leaf wilting ( Figure 1) and decreased nutrient or carbohydrate transportation (Figure 8). A similar result was reported in poplar trees [17].
waterlogging recovery ( Figure 4H). As such, oxidative stress was sustained and even ac-celerated the damage inflicted on the plants during the recovery, as shown in Figure 1.

Lignin and Cuticle Biosynthesis was Continuously Inhibited during Post-Waterlogging Recovery
Lignin, a phenolic biopolymer, is primarily synthesized through the phenylpropanoid biosynthesis pathway and is typically concentrated in the secondary cell wall of vascular plants. Lignin plays a vital role in providing mechanical support to plant tissues, as well as facilitating the transportation of nutrients or carbohydrates [40]. Notably, the accumulation of lignin is known to enhance cell wall reinforcement, thereby improving resistance to waterlogging stress [41]. Previous studies have shown that waterlogging stress can reduce lignin content by inhibiting the expression of genes involved in lignin biosynthesis [42]. Our transcriptome data analysis revealed that the DEGs in WS30d-R10d leaves were significantly enriched in the phenylpropanoid biosynthesis pathway ( Figure  5D), with their expression being noticeably inhibited (Table S6). These results were strongly aligned with the observed decrease in lignin content in WS30d-R10d leaves (Figure 6). As such, we hypothesize that waterlogging may reduce the biosynthesis of lignin in the stressed leaves, ultimately leading to leaf wilting ( Figure 1) and decreased nutrient or carbohydrate transportation (Figure 8). A similar result was reported in poplar trees [17].
The cuticle is a layer of fatty substances that envelops above-ground plant organs, composed mainly of cutin and waxes [43]. As reported by Tellechea-Robles et al. (2019), wetland plants have adapted to waterlogging stress by augmenting the thickness of their cuticle on leaves [44]. In this study, transcriptome analysis revealed that the DEGs were significantly enriched in several pathways during post-waterlogging recovery, namely "fatty acid metabolism", "fatty acid biosynthesis", "fatty acid elongation", and "cutin, suberin, and wax biosynthesis" ( Figure 5D). These pathways may reduce the production of cutin and wax, leading to thinner leaves and a diminished support function.  The cuticle is a layer of fatty substances that envelops above-ground plant organs, composed mainly of cutin and waxes [43]. As reported by Tellechea-Robles et al. (2019), wetland plants have adapted to waterlogging stress by augmenting the thickness of their cuticle on leaves [44]. In this study, transcriptome analysis revealed that the DEGs were significantly enriched in several pathways during post-waterlogging recovery, namely "fatty acid metabolism", "fatty acid biosynthesis", "fatty acid elongation", and "cutin, suberin, and wax biosynthesis" ( Figure 5D). These pathways may reduce the production of cutin and wax, leading to thinner leaves and a diminished support function.

Waterlogging Stress Prevented Transportation of Soluble Sugar from Leaves
Maintaining a stable supply of glycolysis, particularly in roots, is a crucial factor in waterlogging tolerance. In contrast, waterlogging-sensitive plants often fail to uphold sufficient carbohydrate levels, leading to tissue death under stress conditions [21]. Previous studies have established that waterlogging-sensitive species frequently experience substantial disruptions in sugar transportation from phloem to root cells during periods of waterlogging, resulting in increased sugar accumulation in stressed leaves [22,28,45,46]. Our current study confirmed that soluble sugars, such as glucose, sedoheptulose, galactose, and sucrose, were significantly accumulated in waterlogged leaves as compared to CK leaves ( Table 2). In addition, we noted excessive accumulation of soluble sugar in the stressed leaves during post-waterlogging recovery (Tables S9 and S7). These findings suggest that sugar transport in R. delavayi seedlings may have been severely impaired by waterlogging stress, which ultimately resulted in a deficiency of available carbohydrates in the roots [28].

Plant Material and Waterlogging Treatment
Three-year-old Rhododendron delavayi seedlings were cultivated in plastic flowerpots (8 cm in diameter and 8 cm in height), filled with a mixture of nutrient soil and humus soil (1:1 in volume). To induce waterlogging stress, the seedlings were placed within a greenhouse (with a photoperiod 16 h/8 h, temperature 22 • C, light intensity 400 µmol m −2 s −1 , relative humidity 60-70%) at Guizhou Normal University.
The flowerpots with the seedlings were placed in trays with length, width, and height of 55 cm, 30 cm, and 4.5 cm, respectively. The trays were filled with water up to a height of 4 cm, submerging the root area of the R. delavayi seedlings. The water level was sustained at 4 cm by daily additions, creating a waterlogging stress that was imposed for 10 days (WS10d), 20 days (WS20d), and 30 days (WS30d). After the waterlogging stress period, the water within the trays was removed, and the seedlings were allowed to recover for 10 days, marked as WS10d-R10d, WS20d-R10d, and WS30d-R10d, respectively. R. delavayi seedlings before the waterlogging stress were used as the control (CK). Three flowerpots were placed within each tray as a repeat, with three biological repeats set for each treatment.

Measurement of Photosynthesis Parameters and Chlorophyll Fluorescence
To measure photosynthesis parameters, we utilized the portable photosynthesis system (LI-6400, LI-COR Corporate, Lincoln, NE, USA), which was coupled with the 6400-02B chamber (6400-02B, LI-COR Corporate, Lincoln, NE, USA) providing stable and activated intensities. The third young leaf from three individuals was selected and measured between 8.00 and 12.00 a.m. Prior to measurement, the leaves were activated using a light intensity of 1000 µmol quanta m −2 s −1 for 20 s. Subsequently, the light intensity (400 µmol quanta m −2 s −1 ) and CO 2 concentration (ambient CO 2 concentration) were set, and the net CO 2 assimilation rate, stomatal conductance, and transpiration rate were recorded. Additionally, we measured chlorophyll fluorescence using the portable photosynthesis system (LI-6400, LI-COR Corporate, Lincoln, NE, USA) with the 6400-40 chamber. After dark adaptation at night, we conducted measurements before dawn, recording the maximum quantum efficiency of photosystem II (Fv/Fm).

Measurement of Hydrogen Peroxide and Lignin
In this study, we extracted 0.5 g of leaves and subsequently homogenized them with 3 mL of 50 mmol/L phosphate buffer (pH = 6.5). Next, we centrifuged the mixture at 1000 rpm for 5 min at 4 • C, following which we added 1 mL of 0.1% titanium sulfate solution to 3 mL of supernatant and mixed it thoroughly. The solution was then centrifuged at 1000 rpm for 10 min at 4 • C. The resulting precipitate was washed with cooled acetone and subsequently dissolved in 2 mmol/L H 2 SO 4 solution. Subsequently, we measured the absorbance at 410 nm and the hydrogen peroxide content was calculated using the standard curve. Meanwhile, lignin was extracted using an assay kit (BC4200, Solarbio Technology Co., Ltd., Beijing, China) and determined according to the instructions.

RNA Extraction, Library Preparation, and Sequencing
After a 30-day period of exposure to waterlogging stress, the aged leaves of the plant under investigation exhibited signs of wilting ( Figure 1). To gain insight into the effects of this stress, we selected CK, WS30d, and WS30d-R10d leaves and roots for RNA sequencing (RNA-Seq). Total RNA was extracted from the samples using Trizol reagent (Thermo Fisher Scientific, Pleasanton, CA, USA) and evaluated for RNA degradation and contamination using agarose gelelectrophoresis (1%). RNA purity was checked using the Nano Photometer spectrophotometer (Thermo Fisher Scientific, Wilmington, NC, USA), and RNA concentration was measured using the Qubit RNA Assay Kit in Qubit 2.0Fluorometer (Thermo Fisher Scientific, Wilmington, NC, USA). RNA integrity was assessed using the RNA Nano6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, USA). The NEBNext UltraRNA Library Prep Kit for Illumina (New England Biolabs, Inc., Beijing, China) was used to generate sequencing libraries according to the manufacturer's recommendations. However, due to the short roots and lignification of R. delavayi seedlings, RNA extraction from these roots deemed unsuitable for database construction (Table S9), and therefore, RNA-Seq from waterlogged root was abandoned. The index-coded samples were clustered using TruSeq PE Cluster Kit v3-cBot-HS (Illumia) on a Cluster Generation System, following the manufacturer's instructions. Once clusters had been generated, library preparations were sequenced on an Illumina Hiseq 2000 platform, with paired-end reads produced. Each sample yielded more than 6.03 Gb of clean data. Sequencing was completed by the Beijing Biomarker Biotechnology Company (Biomarker biotech, Beijing, China). The RNA-Seq data were obtained in three biological replicates.

Transcriptome Assembly and Differentially Expressed Genes (DEGs) Analysis
Transcriptome assembly was accomplished using the Trinity method [47] with min_kmer_cov set to 2 by default, with all other parameters set to default values. The expression abundance of gene was calculated using the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) method [48]. To analyze DEGs, this was performed using the DESeq R package (1.10.1) between CK, WS30d, and WS30d-R10d. The resulting p values were adjusted using Benjamini and Hochberg's approach to control the false discovery rate (FDR). Genes were considered differentially expressed when FDR was less than 0.01 and fold change (FC) was more than |1.5| between the treatments.

Annotation and Classification of DEGs
To annotate the functions of DEGs, we employed BLAST with an E-value 1 × 10 −10as the cutoff to align these sequences against the NCBI non-redundant database, COG, GO, and KEGG databases. This approach allowed us to assign putative functional annotations to each DEG, based on their sequence similarity to sequences in known databases.

Quantitative Real-Time PCR (qRT-PCR) Validation
To validate the RNA-seq results, we selected six genes from differentially expression genes for quantitative real-time PCR (qRT-PCR) analysis. Total RNA was extracted from the samples using the OmniPlant RNA Kit (CW2598S, Cwbiotech, Beijing, China) and reversetranscribed using the TransScript All-in-One First-Strand cDNA Synthesis SuperMix for qPCR Kit (AT341-01, Transgen Biotech, Beijing, China), following the manufacturer's protocol. The primer sequences used for qRT-PCR, including β-actin as an internal control, were listed in Table S10. The qRT-PCR was performed using a Rotor-Gene Q real-time PCR system (Rotor-Gene, Qiagen, Germany). A total of 10 µL 2 X TransStart Top/TipGreen qPCR SuperMix (AQ141-02, Transgen Biotech, Beijing, China) was added to the reaction mixture, following the manufacturer's instructions. The relative expression levels of the selected genes were normalized to the expression level of the β-actin gene. Cycle threshold values were then used to calculate expression levels using the 2 −∆∆Ct method [49].

Gas Chromatograph Coupled with a Time-of-Flight Mass Spectrometer (GC-TOF-MS) Analysis
The leaves from three individual seedlings were sampled and mixed as a duplicate, and the three duplicate samples were used for GC-MS detection. Metabolites were extracted from leave samples following the method described in [50], with some modifications. Approximately 20 mg of the sample was transferred into a 2 mL EP tubes and extracted with 500 µL of pre-cold extraction liquid (V Methanol :V dH 2 O = 3:1). An aliquot of 10 µL of internal standard (adonitol) was added and mixed by vortexing for 30 s. The sample was then homogenized in a ball mill for 4 min at 35 Hz and subsequently treated with ultrasound for 5 min (incubated in ice water bath). The mixture was then centrifuged at 12,000 rpm at 4 • C for 15 min. The supernatant (100 µL) was transferred to a fresh 1.5 mL EP tube before being dried completely in a vacuum concentrator without heating. Next, the extracts were dissolved in 40 µL of methoxyamine hydrochloride (20 mg mL −1 in pyridine) and then incubated at 80 • C for 30 min. Finally, the samples were derivatized with 60 µL of N,O-bis(trimethylsilyl) trifluoroacetamide (BSTFA) reagent (with 1% TMCS, v/v) at 70 • C for 1.5 h. Quality control (QC) samples (a mixture of all samples to be analyzed) were also processed for detection.
Metabolite detection was conducted using an Agilent 7890 gas chromatography system coupled with a Pegasus HT time-of-flight mass spectrometer (Agilent Technologies, Santa Clara, CA, USA). The analysis was conducted in a randomized order after the addition of fatty acid methylesters (FAMEs). ADB-5MS capillary column coated with 5% diphenyl cross-linked with 95% dimethylpolysiloxane (30 m × 250 µm inner diameter, 0.25 µm film thickness; J&W Scientific, Folsom, CA, USA) was employed. An aliquot of 1 µL of the analyte was injected in splitless mode, while helium was used as a carrier gas. The front inlet purge flow was 3 mL min −1 , and the gas flow rate through the column was 1 mL min −1 . The initial temperature was held at 50 • C for 1 min, then raised to 310 • C at a rate of 10 • C min −1 and held at 310 • C for 8 min. The injection, transfer line, and ion source temperatures were maintained at 280 • C, 280 • C, and 250 • C, respectively. The energy was −70 eV in electron impact mode, and the mass spectrometry data were acquired in full-scan mode with an m/z range of 50-500 at a rate of 12.5 spectra per second after a solvent delay of 6.27 min.

Data Preprocessing and Compound Identification
Chroma TOF 4.3X software and the LECO-Fiehn Rtx5 database (LECO Corporation, Benton Harbor, MI, USA) were used for raw peak extraction, database line filtering and calibration, peak alignment, deconvolution analysis, peak identification, and integration of peak area. Both mass spectrum match and retention index match were considered during metabolite identification. Additionally, peaks detected in less than 50% of the QC samples or for which relative standard deviation (RSD) greater than 30% in the QC samples were removed [50,51]. Metabolites separated by GC-TOF-MS were identified using LECO ChromaTOF 4.3X software and the LECO/Fiehn Rtx5 metabolite mass spectral library by matching the mass spectrum and retention index.

Data Analysis
A three-dimensional data matrix was constructed, consisting of metabolite names (tentatively identified by GC-TOF-MS), sample information (three biological replicates of each treatment), and raw abundance (peak area for each tentatively identified metabolite) This matrix was uploaded to MetaboAnalyst 5.0 and analyzed in accordance with the provided instructions. Raw data were normalized using three categories of normalization: none, log transformation, and auto scaling. Multivariate analyses between CK, WS30d, and WS30d-R10d groups were performed using MetaboAnalyst 5.0, including orthogonal projections to latent structures-discriminant analysis (OPLS-DA) and principal component analysis (PCA). Moreover, the value of the variable importance in the projection (VIP) of the first principal component in OPLS-DA analysis was determined. Metabolites with VIP > 1 and p < 0.05 (student t-test) were considered significantly changed between groups.

Statistical Analysis
The data were presented as mean ± standard deviation (SD) and analyzed using either a one-way analysis of variance (ANOVA) or SNK t-test. Analysis was carried out in at least three replicates, with a p value < 0.05 considered statistically significant. The statistical analysis was performed using SPSS 18.0 for Windows. The difference test of gene expression among the three samples was performed using the R version 4.2.1 vegan package.

Conclusions
We utilized R. delavayi seedlings as a system to elucidate the response mechanism of alpine Rhododendron to waterlogging stress. Waterlogging stress was observed to reduce CO 2 assimilation rate, stomatal conductance, transpiration rate, and maximum photochemical efficiency of PSII in the leave. The stress generated excessive H 2 O 2 , leading to oxidative stress. Additionally, waterlogging stress negatively affected lignin and cuticle biosynthesis, leading to leaf wilting, decreased carbohydrate transportation, and impeded sugar transport from leaves to roots through phloem, which ultimately results in sugar accumulation in the stressed leaves ( Figure 8). Taken together, waterlogging stress has an irreversible effect on the R. delavayi seedlings in both physiological and molecular aspects. Regrettably, our study only shed light on the response mechanism of leaves, prompting further studies to develop into the stress response mechanism of young R. delavayi roots.