Transcriptome profiling reveals molecular mechanisms underlying waterlogging tolerance in daylily

From a historical perspective, daylily (Hemerocallis fulva) has been cultivated in China for over 2,000 years. It not only has high ornamental and practical value, but also carries deep cultural significance, being honored as the Chinese 'Mother Flower'^[1]. As a perennial monocot, daylily belongs to the subfamily Hemerocallidoideae of Asphodelaceae, with 14 species being distributed in the worldwide, 11 of which are native to China. A large quantity of daylily wild species and cultivars have been found in China, which is thus well-known as the distribution center of daylily germplasm resources^[2].

Daylily plants normally possess strong adaptability to external environmental conditions among the most widely used perennial flowers. In recent years, however, the ecological environment has suffered substantial damage. The frequent summer flooding severely threatens the growth and development of daylily plants, posing a serious challenge to urban landscaping. Due to the genetic constraints of plants themselves, and the challenges brought about by environmental changes to their growth adaptability, the quality of introduced daylily cultivars has gradually decreased, with submergence stress emerging as one of the main limiting factors for the conservation and landscaping applications of daylily^[3]. Therefore, the establishment of an efficient evaluation system for daylily submergence tolerance is crucial for the sustainable production of daylily. On the other hand, a high-throughput RNA sequencing (RNA-seq) can be used to elucidate the molecular mechanisms of daylily root responses to waterlogging, thereby laying a solid theoretical foundation for the genetic breeding of daylily.

Water, serving as the primary basis of biological system, is an important ecological factor influencing plant growth and development. Water itself is not harmful to plants, but the hypoxia caused by excess water and accumulation of carbon dioxide (CO₂), and other toxic substances are the major constraining factors for plant life^[4]. Oxygen uptake plays a decisive role in plant survival under hypoxic stress. Plants are partially or completely submerged upon exposure to waterlogging, reducing light intensity, soil nutrient penetration, and diffusion of internal gases. This process may be accompanied by the proliferation of numerous microorganisms and toxic substances, resulting in significant damage to the internal structure and external morphology of plants. The susceptibility of the host to various biotic and abiotic stresses can therefore be increased^[5,6]. Previous reports have indicated that the damage caused by flooding stress on plants can be categorized into five major aspects: reactive oxygen species (ROS) damage, energy depletion hypothesis, toxicity of anaerobic respiration byproducts, cytoplasm acidification, and disruption of endogenous hormone metabolism^[7−10].

Plants have evolved complex physiological and molecular mechanisms to cope with flooding stress during the long-term process of natural environment adaptation^[11−14]. For instance, ethylene (ETH) is a pivotal endogenous hormone in regulating plant responses to waterlogging. It is synthesized by 1-aminocyclopropane-1-carboxylic acid synthase (ACS) and oxidase (ACO). ETH and ROS signaling lead to programmed cell death (PCD) in root tissues, thereby facilitating the formation of aerenchyma tissues involved in gas exchange during the flooding process^[15]. Abscisic acid (ABA) and gibberellins (GAs) are another two important hormones in plant adaptation to waterlogging through the regulation of shoot elongation^[16]. GAs promote internode elongation by reducing growth-inhibiting proteins and starch and loosening cell walls^[17]. ABA modulates the transpiration, germination, dormancy, and adaptation to stress, and especially imposes an inhibiting impact on root elongation^[18,19]. Application of ETH decreases ABA accumulation but increases GA production in the waterlogged plants^[20]. However, the detailed crosstalk between these hormones under submergence still remains elusive.

Additionally, it has been reported that flooding stress results in relatively high expression levels of pyruvate decarboxylases (PDCs), and alcohol dehydrogenases (ADHs), which can promote temporary adenosine triphosphate (ATP) demand through alcohol fermentation^[21−23]. Overexpression of these genes in transgenic plants has improved tolerance to flooding stress compared to wild-type plants^[24,25]. Transcription factors (TFs) play a critical role in specifically regulating the downstream genes in response to stress. It has been reported that the regulation of genes under waterlogging or hypoxic conditions may depend on ethylene response factors (ERFs) of group VII, which serve as the primary inducers of several hypoxia-responsive genes^[26,27]. The submergence 1A (Sub1A) gene in rice confers tolerance to flooding stress by suppressing ETH synthesis, GA response, and carbohydrate consumption^[28]. The barley HvERF2.11 gene regulates the expression of superoxide dismutase 1 (SOD1), peroxidase 1 (POD1), and ACO1 in transgenic Arabidopsis plants by increasing the activity of antioxidant enzymes and ethanol dehydrogenase, thereby enhancing plant tolerance to submergence^[29]. Overexpression of ZmEREB180, a maize ERFVII, promotes adventitious root formation, dynamic ROS balance, and PCD in epidermal cells, thus improving plant survival under flooding stress^[30]. However, the transcriptional regulatory mechanisms of waterlogging tolerance are not fully understood.

To identify the potential submergence-tolerant daylily cultivars and elucidate their molecular mechanisms in response to waterlogging, 23 daylily cultivars and species were selected to observe their phenotypic changes under submergence and an effective submergence tolerance evaluation system established. Through RNA-seq, the molecular mechanisms of daylily root response to waterlogging were preliminarily elucidated, and the variation in gene transcription between tolerant and sensitive daylily cultivars was thus investigated.

All the experimental materials were introduced from American floral company to Dezhou city, Shandong province, China. All of them are internationally recognized as excellent popular cultivars with different plant types, petal shapes, flower colors, fragrances, and ploidy levels, thereby possessing great ornamental and economic significance. The evaluation of their sensitivities to waterlogging stress will be helpful for growers to understand the suitable planting areas in China. They were previously introduced into the botanical garden of Shanghai Institute of Technology (Shanghai, China). A total of 23 daylily cultivars and species were planted in the field, with a consistent growth cycle and a growth period of three years. In the autumn of November, the plants were uniformly pruned (leaving 5–8 cm of the above-ground portion and root systems), and then planted in 2-gallon flower pots with an equal amount of soil mix (perlite : vermiculite : turf = 3:3:4). After winter natural low temperature-induced dormancy release, they were transferred to a growth chamber at 25/20 °C day/night temperature with a 16/8 h light/dark photoperiod and relative humidity of 75% to ensure the identical growth environment. In the late spring of May, when the vigorous growth season of daylilies occurs, the seedlings with similar growth status were selected for the treatment with submergence stress.

Submergence stress treatment

Flower pots with daylily plants were placed in plastic boxes, then submerged at a water level 5–6 cm above the soil surface. The water was replenished at regular intervals to maintain a stable water level, with the submergence stress lasting for 10 d. After waterlogging for 10 d, each daylily cultivar was subjected to recovery for 40 d, with normal maintenance. Ten independently biological replicates were used for each treatment. In the process of waterlogging, the morphological characteristics of the above-ground parts of the daylily plants were observed and recorded, including leaf shape, leaf color, stem color, and morphological changes (the height of newly grown seedlings after 40 d of recovery). The damage index upon exposure to submergence was categorized into five levels (Table 1). After determining the grading range of each indicator, the total score was calculated according to a grading standard. The submergence tolerance of various daylily cultivars or species was comprehensively evaluated based on the established scheme.

RNA extraction and transcriptome sequencing

Based on the submergence tolerance evaluation of daylily plants, the highly submergence-tolerant 'Here Lies Butch', and highly submergence-sensitive 'Bonibrae Sharky', both of which are tetraploid, were selected as experimental materials for RNA-seq. Their root samples were harvested at 0 and 24 h of waterlogging (the root samples of 'Here Lies Butch' at 0 and 24 h were referred to as B1 and B5; the root samples of 'Bonibrae Sharky' at 0 and 24 h were referred to as S1 and S5). Total RNA from the roots was extracted using the Trizol method (Vazyme Biotech, Nanjing, China). The RNA concentration and quality, based on OD260/OD280 and OD260/OD230 values, were measured using a microplate reader. The RNA integrity was assessed using 1.5% agarose gel electrophoresis with GelRed stained (Biotium, Hayward, CA, USA). The cDNA libraries were then constructed using the VAHTS Universal V6 RNA-seq Library Prep Kit for Illumina according to the manufacturer's instructions. After quality control of the constructed libraries through the Agilent 2100 Bioanalyzer, the sequencing was performed using the Illumina HiSeqX Ten sequencer (OE Biotech, Shanghai, China).

De novo assembly and differentially expressed gene analysis

The adaptor, low-quality, and unknown nucleotides of raw reads were removed to generate clean reads, which were de novo assembled through Trinity^[27]. The assembled gene fragments were further enriched through a threshold of length ≥ 300 bp by TIGR Gene Indices clustering tools^[31]. The resulting products were considered as unigenes in the transcriptome database. They were annotated against the public protein databases, including NR, Swiss-Prot, COG, and Kyoto Encyclopedia of Genes and Genomes (KEGG). The fragments per kilobase of transcript per million mapped reads (FPKM) were determined based on the read counts of unigenes using the bowtie2 and eXpress^[32]. FPKM values represent the transcript abundances of unigenes through Cuffinks. The differentially expressed genes (DEGs) were identified using DESeq2 according to a cutoff of fold change ≥ 2.0 and false discovery rate (FDR) ≤ 0.05. Gene Ontology (GO) classification of DEGs was performed using Blast2GO and the WEGO program^[33,34]. To clarify the metabolic pathways associated with DEGs, the unigenes were mapped to the KEGG database^[35]. Statistical significance of enriched pathways was examined using a cutoff of FDR ≤ 0.05.

Reverse transcription-quantitative PCR assay

Under the treatment with submergence, the daylily roots at 0 and 24 h were taken as experimental materials. RNA extraction was conducted as mentioned above, and the cDNA samples were synthesized using a HiScript II Q RT SuperMix for qPCR Kit (Novogene Biotech, Jiangsu, China). Reverse transcription-quantitative PCR (RT-qPCR) was performed to validate the accuracy of the transcriptome data, using the CFX Connect system (BIO-RAD, Hercules, CA, USA) in combination with the SYBR Premix Ex Taq II (Takara, Otsu, Shiga, Japan). The daylily gene TAP42 interacting protein of 41 kDa (TIP41), serves as an internal control. Specific primers were designed using Primer 5.0 software (Supplementary Table S1). The formula based on 2^−ΔΔCᴛ was used to calculate relative expression levels. A minimum of three biological replicates were used for each analysis.

Principle component and interaction network analysis

The correlations between genes were calculated according to Pearson's correlation coefficient (r) through Excel 2016. The gene expression heatmaps based on Log₂-transformation FPKM values were generated using the TBtools^[36]. Principal component analysis (PCA) was conducted using SIMCA-P software (Umetrics, Umea, Sweden). The interaction network was built in the R environment (www.r-project.org) with a coefficient of R ≥ 0.5 or R ≤ −0.5. The interaction results were visualized using the Cytoscape tool^[37]. Within the network, relatively larger graphs (circles or triangles) with more nodes indicate more significant correlations, and conversely the smaller ones with less nodes suggest less significant correlations. The hub genes refer to the interacting genes with close relationships.

Statistical analysis

At least three independently biological replicates were used for each experiment in this study. Statistical significance of the data among samples was evaluated through Duncan's multiple range test at p < 0.05 and denoted by various letters. The JMP software (SAS Institute Inc., Cary, NC, USA) was employed for determining the significance of difference.

To investigate the morphological responses of daylily plants to waterlogging, five daylily cultivars were chosen for submergence treatment. The phenotypic changes in daylily leaves became evident after 4 d of submergence. From 0 to 2 d of submergence, all the plants tested remained healthy with the green leaves extending outward and submerged stems maintaining the natural color. Subsequently, the symptoms of daylily leaves became evident after 4 d of submergence. The middle and lower leaves of plants began to curl and droop. Under submergence for 8 d, almost half or over half of the lower leaves turned yellow. Finally, most of the plant leaves showed chlorosis and withering at 10 d of submergence, and the submerged stems turned white and released unpleasant odors. Different daylily cultivars exhibited variable symptoms under waterlogging. As the ability of the plant to resume growth post-submergence stress is a crucial indicator of tolerance, we further treated the waterlogged daylily plants with 40 d of normal maintenance. Consistent with the waterlogging symptoms, the relatively tolerant daylily cultivars, such as 'Here Lies Butch', 'Barbara Barnes', 'Alexander Hay', and 'Face of the Stars', displayed new and vigorous shoot growth after recovery (Fig. 1). These results indicate that the submergence stress can remarkably change the phenotypes of daylily plants, with the tolerance of different daylily cultivars varied.

Figure 1. 

Representative symptoms of different daylily cultivars at different days after waterlogging treatment. The symptoms were recorded every two d post treatment. After waterlogging for 10 d, each daylily cultivar was subjected to recovery (Rec.) for 40 d of normal maintenance.

Establishment of a submergence tolerance evaluation system for daylily plants

Based on the observed phenotypic changes of daylily plants under submergence, a number of indicators, including leaf shape, leaf color, stem color, and morphological changes after 40 d of recovery, were selected to assess submergence tolerance. The scoring range and criteria for each indicator were determined following the previously described standards^[38]. The morphological traits of above-ground parts of daylily were observed and recorded during 10 d of continuous submergence and 40 d of recovery post-submergence stress. Total scores of four indicators were calculated without weighting to determine the five-level submergence tolerance.

Next, the evaluation system was used to examine the waterlogging tolerance of one wild daylily species, H. fulva, and 22 cultivars. The tolerance evaluation revealed five highly submergence-tolerant (1) cultivars ('Darla Anita', 'Chesapeake Crablegs', 'Here Lies Butch', 'Mary Reed', and 'Stella de Oro'), two submergence-tolerant (2) cultivars ('Happy Hooligan' and 'Barbara Barnes'), and nine medium submergence-tolerant (3) cultivars and species ('Rosy Rhino', 'Alexander Hay', 'Big City Eye', H. fulva, 'Autumn Red', 'Bela Lugosi', 'Chicago Apache', 'Strawberry Candy', and 'Chicago Fire'), three medium submergence-sensitive (4) cultivars ('Face of Stars', 'Water Dragon', and 'Pardon Me'), and four highly submergence-sensitive (5) cultivars ('Bonibrae Sharky', 'Get Jiggy', 'Chicago Sunrise', and 'Orange Crush'). Different daylily cultivars exhibited significant variations in response to submergence stress in terms of appearance and morphology (Table 2).

Transcriptome sequencing analysis of daylily roots under waterlogging

To elucidate the molecular mechanisms of waterlogging tolerance in daylily, the highly tolerant cultivar 'Here Lies Butch' (B), and highly sensitive cultivar 'Bonibrae Sharky' (S) were used as experimental materials for RNA-seq, and 12 cDNA libraries were then constructed from root samples at 0 h (B1 and S1) and 24 h (B5 and S5) after submergence, with three libraries for each sample. The original data of each sample showed a Q30 distribution ranging from 92.13% to 92.82%, with an average GC content of 46.88%. A total of 87,573 unigenes were assembled (Supplementary Table S2), with a total length of 69,826,487 bp and an average length of 797 bp. Alignment of clean reads with the unigenes resulted in a total mapping rate of 86.37% to 87.99% and a unique mapping rate of 65.31% to 66.60% (Supplementary Table S3). These results indicate that the sequencing data is specific to daylily and suitable for subsequent data analysis. PCA showed that PC1 and PC2 explained 66% and 21% of the variance, respectively, suggesting highly significant variations among four samples (Fig. 2a). The correlation test showed that the correlation coefficient is approximately 1.0 for three replicates of each sample, revealing a close genetic distance among them (Fig. 2b). The cluster analysis based on the sample-to-sample distance showed that three replicates of each sample clustered together (Fig. 2c), which is in accordance with the above results.

Figure 2. 

Transcriptome analysis of daylily roots under waterlogging treatment. (a) Principal component analysis (PCA) of daylily root transcriptome data for each sample. B1 and B5, 'Here Lies Butch' root samples at 0 and 24 h after waterlogging; S1 and S5, 'Bonibrae Sharky' root samples at 0 and 24 h after waterlogging. (b) The correlation test among samples in daylily root transcriptome data. (c) Cluster analysis of all samples based on sample-to-sample distances. (d) Number of differentially expressed genes (DEGs) for different pairwise comparisons. (e), (f) Venn diagrams of DEGs for different pairwise comparisons.

Compared to S1, a total of 9,297 DEGs were observed in S5, including 4,702 up-regulated genes and 4,595 down-regulated genes, and 8,769 DEGs were found in B5 compared to B1 (3,922 up-regulated genes and 4,847 down-regulated genes). There were 15,426 DEGs (7,856 up-regulated genes and 7,570 down-regulated genes) between the B1 vs the S1 group, and 15,936 DEGs (6,768 up-regulated genes and 9,168 down-regulated genes) between B5 vs S5 group (Fig. 2d; Supplementary Tables S4 and S5). A total of 31,527 DEGs were present in the four comparison groups, and 319 DEGs were shared by the four pairwise comparison groups (Fig. 2e). Specifically, 1,370 DEGs were shared by the pairwise comparisons of B5 vs B1, and S5 vs S1 in terms of down-regulation. For the up-regulation, the number of shared DEGs between B5 vs B1, and S5 vs S1 was 1,334 (Fig. 2f).

To verify the reliability of transcriptome data, nine DEGs were selected for further examination using the RT-qPCR method. These DEGs included glutathione S-transferase (GST), gretchen hagen 3 (GH3), aldolase (ALDO), β-glucosidase B (BGIB), hexokinase (HK), phosphoglycerate mutase (PGM), ADH1, ADH5, and E3.2.1.21. The final quantitative results of DEGs were consistent with the transcriptome data, confirming the reliability of the RNA-seq (Supplementary Fig. S1).

Functional classification and pathway analysis of DEGs

A total of 28,615 DEGs were annotated to the GO database, with all DEGs classified into three functional categories: biological process (BP), cellular component (CC), and molecular function (MF). Most of the DEGs in the BP category were assigned to cellular processes, metabolic processes, and regulation of biological processes. The three most abundant subcategories in the CC classification were cell, cell part, and organelle. For the MF category, the three most abundant structural components were binding activity, catalytic activity, and transporter activity (Fig. 3a).

Figure 3. 

GO and KEGG analysis of daylily root transcriptome data in response to waterlogging stress. (a) The GO functional classification of differentially expressed genes (DEGs) in all pairwise comparisons. (b) KEGG pathway enrichment of DEGs between different groups.

Additionally, 10,392 unigenes were annotated to 199 pathways based on KEGG, involving metabolism, genetic information processing, environmental information processing, and cellular processes. During submergence from 0 to 24 h, a total of 53 significantly enriched metabolic pathways were identified in the S5 vs S1 group, while 27 pathways were identified in the B5 vs B1 group (Fig. 3b). The enrichment of more metabolic pathways in the S5 vs S1 group indicates that the highly submergence-sensitive 'Bonibrae Sharky' mobilized a greater number of metabolic pathways to resist submergence. Additionally, 17 metabolic pathways showed significant differences in both the B5 vs B1, and S5 vs S1 groups, suggesting their important roles in submergence tolerance of daylily. Among them, the pathways such as ascorbate and aldarate metabolism (ko00053), starch and sucrose metabolism (ko00500), plant hormone signal transduction (ko04075), glutathione metabolism (ko00480), arginine and proline metabolism (ko00330), and glycolysis/gluconeogenesis (ko00010) were highlighted.

Totally 21 and 30 metabolic pathways were significantly enriched in the B1 vs S1, and B5 vs S5 groups, respectively (Fig. 3b). Eight metabolic pathways consistently showed significant differences in both the B1 vs S1, and B5 vs S5 groups. Among them, nitrogen metabolism (ko00910), phenylalanine metabolism (ko00360), phenylpropanoid biosynthesis (ko00940), ribosome (ko03010), cutin, suberine, and wax biosynthesis (ko00073), and tropane, piperidine, and pyridine alkaloid biosynthesis (ko00960) are related to submergence stress. The results suggest that these metabolic pathways exhibit strong activities in 'Here Lies Butch' responding to waterlogging.

In the B5 and S5 groups, 30 metabolic pathways were significantly enriched. These pathways involve a number of amino acid metabolism, energy metabolism, and biosynthesis pathways related to submergence stress. Among these, the pathways such as carbon fixation in photosynthetic organisms (ko00710), glycolysis/gluconeogenesis (ko00010), peroxisome (ko04146), biosynthesis of amino acids (ko01230), tryptophan metabolism (ko00380), limonene and pinene degradation (ko00903), pyruvate metabolism (ko00620), and glutathione metabolism (ko00480) responded to submergence stress (Fig. 3b). This may explain why 'Here Lies Butch' exhibited stronger tolerance to waterlogging in comparison to 'Bonibrae Sharky'.

Identification of DEGs from various pathways associated with waterlogging

To understand the responses of daylily roots to waterlogging, the DEGs involved in different metabolic pathways, as previously described^[16], related to waterlogging were identified. A large number of DEGs were found in the pathways of ETH biosynthesis, ROS production and scavenging, alcoholic fermentation, ABA biosynthesis, GA biosynthesis, and PCD (Fig. 4). In detail, compared to S1 and S5, most of DEGs involved in ETH biosynthesis pathway and ROS production and scavenging pathway exhibited decreased transcription in B1 and B5, especially ACS1, respiratory burst oxidase homolog D (RbohD), RbohFb, SOD1, SOD2, SOD4, catalase 1a (CAT1a), and CAT2. This implies that the highly tolerant cultivar 'Here Lies Butch' may have less ETH and ROS generation under waterlogging than the highly sensitive cultivar 'Bonibrae Sharky'. Conversely, the B1 and B5 samples showed remarkably increased transcription of more DEGs in the alcoholic fermentation pathway, including lactate dehydrogenase 1 (LDH1), LDH2, ADH1a, and ADH7, the ABA biosynthesis pathway, including zeaxanthin epoxidase 2 (ZEP2), short-chain dehydrogenases/reductase 1 (SDR1), SDR2, and SDR3b, and the PCD pathway, including ataxia telangiectasia and Rad3‐related (ATR), ataxia telangiectasia mutated 2 (ATM2), and caspase (CASP), in comparison to S1 and S5. Among them, the alcoholic fermentation and PCD have been revealed to participate in energy production, adventitious root growth, and aerenchyma formation, and therefore they are probably the essential factors for the high tolerance of 'Here Lies Butch'. Besides, some critical GA biosynthesis-associated DEGs, including ent-copalyl diphosphate synthase (CPS), ent-kaurenoic acid oxidase 1a (KAO1a), and KAO2a, in B1 and B5 showed higher transcript abundances than those in S1 and S5, whereas the other three GA pathway DEGs KAO1b, KAO2b, and GA20-oxidase 2 (GA20ox2) displayed lower expression levels (Fig. 4). This indicates that GA may play a complicated role in response to submergence stress between the tolerant and sensitive daylily cultivars.

Figure 4. 

Expression profiles of differentially expressed genes in various pathways associated with waterlogging stress responses through transcriptome analysis. These differentially expressed genes (DEGs) are indicated by italics. Different colors in the heatmap represent variable transcription ranged from low to high levels for each transcript. Transcript levels of DEGs were Log₂-transformed. Solid lines ending with arrows indicate positive regulation, whereas the ones ending with short perpendicular lines represent negative regulation.

Identification of differentially expressed TFs related to waterlogging responses

To dissect the regulatory mechanisms of waterlogging tolerance in daylily roots, the TFs with differential expression in daylily transcriptome data was further analyzed. For the pairwise comparison of B1/S1 and B5/S5, 189 and 204 differentially expressed TFs, belonging to 15 distinct TF families, were identified, respectively (Supplementary Tables S6, S7). The highest number of TFs were found to be the members of the ZFP family, followed by WRKY, bHLH, MYB, NAC, and AP2/ERF families (Fig. 5a). A set of candidate TFs with highly significant differential transcription were selected for further investigation. According to RNA-seq data, transcript abundances of ERF71, aintegumenta-like 5 (AIL5), bHLH113, bHLH144, MYB1, MYB12, NAC76, WRKY65, WRKY75, and ZFP6 in B1 and B5 were markedly higher than those in S1 and S5. The data demonstrate that these TFs possibly function as positive regulators of daylily tolerance to waterlogging. The remaining TFs, except for MYB39, showed relatively higher expression levels in S1 and S5 (Fig. 5b), implying their possible negative roles in regulating waterlogging tolerance.

Figure 5. 

Analysis of differentially expressed transcription factors in waterlogging-treated daylily roots through transcriptome analysis. (a) Number of differentially expressed transcription factors in response to waterlogging stress. (b) Transcript levels of a number of transcription factors from ERF, ARF, bHLH, MADS, MYB, NAC, WRKY, and ZFP families, in various samples under waterlogging treatment. Error bars indicate standard deviation (SD) of the means from three biological replicates. Different letters above the bars suggest statistical significance, which is determined by Duncan's multiple range test at p < 0.05.

Interaction of differentially expressed TFs with genes in various pathways under waterlogging

To investigate the possible interaction between the above TFs and genes in various pathways, PCA was carried out based on their expression profiles. The results showed that some genes in various waterlogging-associated pathways were correlated with different TFs. In the ETH biosynthesis pathway, specifically, ACS1 was pooled together with NAC100, ZFP62, ZFP63, ARF16, MYB4, MYB39, and ZFP53, while ACO1b was pooled together with bHLH49, ZFP19, ERF1, WRKY17, and WRKY24 (Fig. 6a). This suggests a closer transcriptional relationship between them. For the genes involved in ROS production and scavenging, RbohD, CAT2, and SOD2/SOD4 were grouped with MYB4/ZFP63, bHLH149/MADS3/bHLH49, and MADS2/MADS23, respectively (Fig. 6b). ADH1b, a key gene in alcoholic fermentation pathway, was next to MYB39, ZFP53, ZFP63, ARF16, NAC100, and MYB4 (Fig. 6c). The ABA biosynthetic gene SDR2 had a close distance with ZFP62/WRKY65/bHLH113/AIL5, and SDR3a/ZEP1 with NAC100/ARF16/ZFP63/MYB4/ZFP53/ MYB39 (Fig. 6d), which were also pooled together with the GA biosynthetic genes KAO1b, KAO2b, and GA20ox2 (Fig. 6e). Moreover, the PCD-associated gene ATM2 displayed a high correlation with bHLH113, NAC76, WRKY65, and AIL5 (Fig. 6f).

Figure 6. 

Principle component analysis of the selected differentially expressed transcription factors with the genes in various pathways associated with waterlogging tolerance. The functional genes in various pathways of (a) ethylene biosynthesis, (b) reactive oxygen species (ROS) production and scavenging, (c) alcoholic fermentation, (d) abscisic acid (ABA) biosynthesis, (e) gibberellin (GA) biosynthesis, and (f) programmed cell death, are denoted by red solid circles, while the selected transcription factors from various families are marked by green solid circles. The genes with relatively close correlation distance are indicated by dashed circles.

To further study the correlation of candidate TFs with various waterlogging-related genes, an expression data-based interaction network was then built. The network revealed the high correlation between ACS1 and ZFP62/NAC100, ACO1b and ERF1/ZFP19/bHLH49, RbohD and MYB4/ZFP63, SOD2/SOD4 and MADS23, and ADH1b and MYB39 (Fig. 7). In contrast, the other TFs and pathway genes displayed relatively large correlation distance, representing no remarkable transcriptional link between them. The results indicate that these genes in response to waterlogging may be specifically modulated by the corresponding strongly correlated TFs.

Figure 7. 

An interaction network between the selected differentially expressed transcription factors with the genes in various pathways associated with waterlogging tolerance. Grey solid or dashed lines indicate positive or negative regulation, respectively. Larger circles or triangles with more nodes indicate more significant correlations, whereas smaller ones with less nodes represent the opposite correlations. The genes with relatively close correlation distance are indicated by dashed circles.

The primary influence of submergence stress on plants is hypoxia, resulting from insufficient oxygen in water to meet the normal growth and development of plants^[39]. As plants are inherently immobile, their resilience to submergence-induced hypoxia is normally enhanced through physiological and biochemical responses such as the activation of antioxidant enzyme systems and modification of external morphology^[40,41]. The waterlogging tolerance evaluation of daylily cultivars based on the external morphological changes is a prerequisite for correlating alterations in antioxidant enzymes and osmolytes with waterlogging tolerance. Considering the morphological characteristics and growth states of daylily, this study established an evaluation system for submergence tolerance by exposing plants to artificial submergence stress. Significant phenotypic changes, including leaf shape, leaf color, stem color, and morphological changes after 40 d of recovery were observed in 23 tested daylily cultivars and species, such as the highly tolerant cultivar 'Here Lies Butch' and highly sensitive cultivar 'Bonibrae Sharky', during submergence treatment. It suggests that phenotypic variation can serve as a crucial indicator for assessing plant resistance traits. The evaluation outcomes showed variable tolerance grades of daylily plants to waterlogging, which would be helpful to guide the selection and cultivation of daylily cultivars in different rainfall areas. Nonetheless, as the qualitative descriptions, the phenotypic traits seem not to be sufficient for the establishment of a perfect evaluation system. The incorporation of quantitative indicators, such as chlorophyll content and biomass measurements, is quite essential to improve the objectivity and reliability of the evaluation.

In this study, the post-waterlogging recovery duration (40 d) of daylily plants was determined based on a previous report in coffee plants^[42]. It was also found that the recovery duration was less than 40 d in many plant species post waterlogging, such as 10 d of recovery in the Rhododendron tree^[43], 14 d of recovery in wheat^[44], and 15 d of recovery in forage grasses^[45]. These reports indicate that 40 d of recovery after waterlogging in daylily plants is sufficient to confirm the restoration of normal growth.

For RNA-seq analysis, the highly tolerant cultivar 'Here Lies Butch' (B) and the highly sensitive cultivar 'Bonibrae Sharky' (S) were used as experimental materials. It is worth mentioning that three cultivars 'Bonibrae Sharky', 'Get Jiggy', and 'Chicago Sunrise', showed the same lower score of 6. In comparison to the other two cultivars, 'Bonibrae Sharky' is an excellent tetraploid plant with rich gene dosage and strong physiological stability. It features a high number of flower buds, multiple branches, and especially highly double flowers that are pink with a double edge of dark purple and gold teeth. It has very high application potential in the international flower market and is therefore a relatively ideal material for studying the sensitivity of daylily plants to submergence stress.

After 24 h of submergence stress, 17 metabolic pathways were significantly enriched in both cultivars 'Here Lies Butch' and 'Bonibrae Sharky', indicating their involvement in the molecular responses of daylily roots to waterlogging. Furthermore, 8,769 DEGs were identified in 'Here Lies Butch' roots and 9,297 DEGs in 'Bonibrae Sharky' roots. This suggests that the roots of the two daylily cultivars have distinct responsive mechanisms to short-term submergence threat. In the comparison between the B5 and S5 groups, a total of 30 metabolic pathways were significantly enriched, with these pathways being unique to 'Here Lies Butch' or exhibiting higher metabolic activities compared to 'Bonibrae Sharky'. These pathways involve various functions, including energy metabolism, amino acid metabolism, antioxidative processes, and others. Energy metabolism-related pathways can provide daylily plants with sufficient energy required under submergence stress^[46], while amino acid metabolism-related pathways play a multifaceted role in plant tolerance to waterlogging. They assist plants in adapting to adverse conditions by regulating metabolism, antioxidation, and protein synthesis, thereby facilitating maintain cellular homeostasis and physiological functions^[47]. Antioxidation-related metabolic pathways contribute to the maintenance of a dynamic ROS balance to alleviate oxidative damage^[48,49]. These explanations, as provided above, may clarify why the cultivar 'Here Lies Butch' exhibits stronger submergence tolerance when compared to 'Bonibrae Sharky'.

For the enrichment of DEGs, we mainly focused on a few crucial pathways closely related to waterlogging response, including ETH biosynthesis, ROS production and scavenging, alcoholic fermentation, ABA biosynthesis, GA biosynthesis, and PCD. A general trend at the transcript level showed that the ETH and ROS pathways were relatively attenuated but the pathways of alcoholic fermentation, ABA, and PCD were promoted in the highly tolerant 'Here Lies Butch' (Fig. 4). A controversial finding showed that waterlogging treatment strongly increases ETH production in tolerant genotype (126) of rapeseed compared to sensitive genotype (85)^[50]. Exogenous ETH application reduces the sensitivity of maize seedlings to waterlogging by affecting root architecture and aerenchyma formation^[51]. Conversely, the treatment with anti-ETH agent aminoethoxyvinylglycine induces waterlogging tolerance in cotton, leading to decreased fruit abscission and yield losses^[52]. We thus hypothesize that ETH probably plays multifaceted roles in submergence tolerance relying on its application period, plant species, or organs. Moreover, the positive role of GAs in waterlogging tolerance has been previously reported by inducing internode elongation, photosynthetic rate, and accumulation of total phenolics and flavonoids^[53]. However, no clear transcriptional variation trend of GA biosynthetic genes was found through RNA-seq analysis (Fig. 4). Due to the unique morphological characteristics of daylily plants, like short stems, the function of GAs in the responses of daylily plants to waterlogging may differ from the other plant species. A detailed investigation of GAs' role in submergence stress needs to be performed in subsequent studies.

Under submergence-elicited hypoxia, the ROS pathway can be greatly activated to generate superoxide anion radical (O₂⁻) and hydrogen peroxide (H₂O₂). When the ROS levels exceeds the plant's clearance capability, it leads to ROS overproduction and lipid peroxidation, ultimately resulting in plant damage such as root rot and leaf wilting^[54,55]. Therefore, the key to plant survival under submergence conditions lie in maintaining a balance between ROS production and scavenging^[56]. During submergence stress, antioxidative metabolic pathways such as peroxisome (ko04146), glutathione metabolism (ko00480), and ascorbate and aldarate metabolism (ko00053) were enriched. Particularly, the peroxisome (ko04146) pathway includes a number of DEGs associated with CAT and SOD. However, many genes encoding CAT and SOD are down-regulated after 24 h of submergence stress compared to 0 h. It implies that the ability of daylily roots to cope with oxidative stress decreases at later stages, leading to ineffective clearance of ROS and exacerbation of cellular oxidation. Alternatively, similar to cucumber, daylily roots when exposed to submergence may not rely on CAT and SOD for the removal of ROS species^[57]. Additionally, plants can adjust their internal substance metabolism by accumulating soluble sugars, prolines, and soluble proteins to enhance their resilience to waterlogging conditions, thus mitigating the adverse effects of submergence on plants^[58]. Proline serves as a compatible osmotic regulator, stabilizer of proteins and subcellular structures, scavenger of ROS, and balancer of redox equilibrium, thereby playing a pivotal role in plant defense against waterlogging^[59]. In the metabolic pathway of arginine and proline metabolism (ko00330), a total of 32 DEGs were identified. The up-regulation and down-regulation of these genes facilitate the accumulation of prolines in daylily roots, contributing to the maintenance of cell turgor pressure and plant protection^[60,61].

The energy crisis induced by hypoxia under submergence conditions is a critical factor adversely affecting plant growth and development^[62]. Hypoxic challenge during submergence stress leads to a shift in the energy production of submerged tissues from oxidative phosphorylation (a high-efficiency mode) to glycolysis and fermentation (a low-efficiency mode)^[63]. Additionally, pyruvic acid, the end product of glycolysis, can be converted to acetaldehyde by PDC and further transformed into ethanol by ADH, thus generating nicotinamide adenine dinucleotide (NAD+) and providing substrates for glycolysis^[64,65]. However, the breakdown of glucose by the glycolytic pathway can yield ATP, posing a risk of energy crisis for plants under prolonged submergence stress. Therefore, plants also induce amylase to degrade starch into glucose, facilitating the mobilization of energy reserves. The metabolic pathways of starch and sucrose metabolism (ko00500), pyruvate metabolism (ko00620), and glycolysis/gluconeogenesis (ko00010) serve as the biochemical foundation for energy production in plants, thereby enhancing the submergence tolerance of daylily^[66,67].

Furthermore, plant hormones play a crucial role in regulating plant growth and development, integrating various signaling pathways under abiotic stress, and directly or indirectly participating in hormone homeostasis regulation. In this study, the metabolic pathway of plant hormone signal transduction (ko04075) was enriched with numerous DEGs. These DEGs are predicted to play significant roles in the induction of defensive gene expression and the synthesis of secondary metabolites. They are particularly vital in the defense against environmental stressors. In particular, we found some DEGs in the pathways of ETH, ABA, and GA biosynthesis (Fig. 4).

Previous research has indicated that ETH accumulates significantly under submergence and is relatively unaffected by other environmental factors, making it a crucial element in waterlogging tolerance. Normally, gaseous ETH emitted by plants dissipates quickly into the surroundings. However, when faced with submergence, the ETH synthesis pathway is activated, resulting in the rapid accumulation of endogenous ETH within a short timeframe. This accumulation promotes the development of adventitious roots and stomata in plants^[68−70]. It has been reported that reduced concentrations of indole-3-acetic acid (IAA), a kind of auxin, stimulate the synthesis of ETH in maize roots, consequently prompting the development of aerenchyma tissues^[71]. GAs serve as a key regulator for plants to adapt to submergence stress, influencing carbohydrate metabolism, internode elongation, and the formation of adventitious roots^[72]. In most cases, GAs induce the formation of aerenchyma tissues and adventitious roots by promoting ETH synthesis. Additionally, ABA acts as a negative regulatory factor in the adaptive mechanisms such as internode elongation and adventitious root formation under submergence stress^[73]. It was found that the tricarboxylic acid (TCA) cycle of ABA-deficient tomato plants was inhibited under flooding stress, thereby decreasing the synthesis of ATP. In rice, ABA antagonizes ETH to function as a repressor of GA and ETH signaling^[74,75]. Furthermore, brassinosteroid (BR), another plant hormone, has been found to interact with the Sub1A gene in rice, influencing GA signaling and regulating internode elongation under flooding conditions^[76]. A crosstalk among these hormones and the pathway genes probably exists to participate in the daylily root responses to waterlogging, which requires further investigation in future work.

It is noteworthy that some differentially expressed TFs belonging to various families were identified from daylily root transcriptome data under submergence stress, and their interaction with genes in various waterlogging-related pathways was further analyzed (Figs. 5−7). A cluster of TFs were pooled together with some key structural transcripts in the pathways. These TFs, such as ERF1, ZFP19, bHLH49, ZFP62, NAC100, MYB4, MYB39, and others, are presumed to participate in waterlogging responses by specifically modulating the downstream transcripts with close correlations. It has been suggested that ERF1 contributes to auxin accumulation and ETH-induced lateral root growth inhibition under waterlogging stress in Arabidopsis^[77]. MYB39 positively regulates Arabidopsis root suberization^[78], which is believed to be implicated in the plant tolerance to waterlogging. These findings largely support the possible important roles of the identified TFs above in waterlogging response. The functions of the candidate TFs and their direct binding to the promoters of the highly correlated structural transcripts should be further investigated.

In this study, 23 daylily cultivars and species were collected for the characterization of submergence tolerance. An effective tolerance evaluation system for daylily plants upon exposure to submergence was established. Transcriptome sequencing was employed to analyze the molecular mechanisms of daylily root responses to submergence stress. The data demonstrated that the two daylily cultivars 'Here Lies Butch' and 'Bonibrae Sharky' exhibited variable expression patterns to waterlogging. A number of DEGs associated with hormone signaling, antioxidant systems, energy metabolism, and other metabolic pathways were identified, therefore providing a crucial basis for the functional investigation of candidate genes and future molecular genetic breeding of waterlogging-tolerant daylily.