Abstract

Background

   Submergence stress is a prevalent abiotic stress affecting plant growth
   and development and can restrict plant cultivation in areas prone to
   flooding. Research on plant submergence stress tolerance has been
   essential in managing plant production under excessive rainfall. Red
   clover (Trifolium pratense L.), a high-quality legume forage, exhibits
   low tolerance to submergence, and long-term submergence can lead to
   root rot and death.

Results

   This study assessed the microstructure, physiological indicators, and
   the key genes and metabolic pathways under submergence stress in the
   root system of red clover HL(Hong Long) and ZY(Zi You) varieties under
   submergence stress at 0 h, 8 h, 24 h, 3 d, and 5 d. Based on 7740
   transcripts identified in the leaves at 0 h, 8 h, and 24 h submergence
   stress, Weighted Gene Co-expression Network Analysis (WGCNA) was
   performed on the differentially expressed genes (DEGs) at 8 h and 24 h.
   Functional annotation of the DEGs in the four key modules was obtained.
   Based on the results, the red clover root system exhibited epidermal
   cell rupture, enlargement and rupture of cortical thin-walled cells,
   thickening of the mid-column, and a significant increase in the number
   of air cavities and air cavity area of aeration tissue with the
   prolongation of submergence stress. The malondialdehyde content,
   relative conductivity, peroxidase, and superoxide dismutase initially
   increased and decreased as submergence stress duration increased. Four
   specific modules (cyan, purple, light cyan, and ivory) closely
   correlated with each stress were identified by WGCNA. The 14 obtained
   Hub genes were functionally annotated, among which six genes, including
   gene51878, gene11315, and gene11848, were involved in glyoxylate and
   dicarboxylic acid metabolism, carbon fixation in photosynthetic
   organisms, carbon metabolism, biosynthesis of pantothenic acid and CoA,
   flavonoid biosynthesis.

Conclusion

   In this study, using WGCNA, the molecular response mechanisms of red
   clover to submergence stress was proposed, and the core genes and
   metabolic pathways in response to submergence stress were obtained,
   providing a valuable data resource at the physiological and molecular
   levels for subsequent studies of submergence stress tolerance in
   plants.

Supplementary Information

   The online version contains supplementary material available at
   10.1186/s12870-024-05804-z.

   Keywords: Trifolium pratense, Submergence stress, Morphological
   structure, Microstructure, Physiological index, WGCNA

Introduction

   Water availability is an essential environmental variable for plant
   growth and development and is an important component of plant
   photosynthesis and organic matter production. According to the
   Intergovernmental Panel on Climate Change (IPCC), annual rainfall will
   significantly increase at high altitudes and in the equatorial Pacific
   at the end of the twenty-first century, and extreme rainfall events
   will be more frequent in low-altitude continental regions and the wet
   tropics [[36]1]. Such regular extreme rainfall events lead to water
   stagnation on the soil surface, resulting in a very low soil oxygen
   content [[37]2]. Rainfall and water table depth affect soil hydrology
   and can lead to plant physiological adaptations significantly affecting
   anatomical features such as root depth and root morphology, as well as
   physiological and metabolic processes [[38]3], and even cause plant
   yield reduction and mortality, thus seriously affecting regional
   ecological and economic development.

   Long-term submergence stress in plants can lead to decreased plant root
   vigor and decreased growth indicators such as root length, root
   biomass, and plant height [[39]4]. Submergence stress in plants can
   also lead to an imbalance in the production and scavenging system of
   reactive oxygen species, resulting in cell membrane lipid peroxidation
   damage [[40]5]. Moreover, PS II reaction centers are damaged when
   plants are subjected to submergence stress, resulting in reduced
   photosynthetic activity and a decrease in chlorophyll content [[41]6,
   [42]7]. Although water is essential for plant growth, submergence
   stress, by disrupting tissue oxygenation, can severely impact normal
   plant growth and development [[43]8] and, thus, is a major constraint
   on the quality and yield of plants.

   The plant root system is the main organ that absorbs water and mineral
   nutrients to support plant growth, and it is the tissue that initially
   senses soil-related stress. Submergence stress can quickly lead to
   hypoxia in the plant root system and inhibit its growth. If the plant
   is completely submerged, the hypoxic stress is further intensified,
   which can cause serious damage and even lead to plant death. Plant
   growth and development are determined by a combination of genetic
   factors and external factors, such as biotic and abiotic stresses
   [[44]9]. During their long-term evolutionary adaptation, plants have
   gradually evolved various morphological structures and molecular
   mechanisms to adapt to adverse environmental conditions [[45]10,
   [46]11]. Flood-tolerant plant varieties play an important role in
   managing extreme precipitation phenomena, reducing the cost of
   large-scale production and operation, expanding planting area, reducing
   yield losses, and improving the utilization of flooded land. Therefore,
   research on plant submergence stress tolerance mechanisms is important
   in breeding cultivars that are less affected by soil flooding.

   Transcriptome sequencing of multiple samples has been widely
   implemented in life science research, and mining biologically
   significant information from the extensive transcriptome data generated
   has become the focus of current transcriptome analyses. Weighted gene
   co-expression network analysis (WGCNA) [[47]12, [48]13] is among the
   most effective analytical approaches for the identification of
   co-expression networks. It is commonly implemented for the
   identification of gene expression correlation modules and to predict
   gene regulatory relationships among genes, and thus identify
   biologically significant target modules and target genes
   [[49]14–[50]16]. More recently, WGCNA approaches have also been applied
   in a wide range of agricultural-related research fields, such as
   genetic marker screening for important traits in plants and animals,
   plant developmental biology, reproductive biology, and pathology.

   Red clover (Trifolium pratense L.) is a perennial cool-season
   herbaceous legume. It is one of the most important forage legume
   cultivated worldwide, including the United Kingdom, the United States,
   New Zealand and China. Red clover has high nutritional and medicinal
   value and is rich in flavonoids, proteins, amino acids, and other
   components [[51]17, [52]18]. It is characterized by good palatability,
   easy digestibility, fast regeneration, and a high biomass yield. It is
   rich in flavonoids and [[53]19] and is an important raw material for
   extracting flavonoids used in women's health products and medicinal
   products [[54]20]. It can be used as a green manure plant and companion
   plant [[55]21] and also as a forage plant with significantly higher
   utilization of protein in silage than alfalfa (Medicago sativa)
   [[56]22]. It prefers warm and humid climates, has poor drought
   tolerance but a significant tolerance to high soil moisture conditions,
   and is suitable for growing in regions with abundant annual
   precipitation and neutral alkaline soils. However, red clover does not
   tolerate submergence, which causes rotting and eventually plant death.
   Nevertheless, large differences in submergence resistance have been
   identified in clover cultivars.

   Extensive studies have been undertaken in recent years on plant (mainly
   Arabidopsis and rice) submergence-associated stresses, such as
   submergence, waterlogging, hypoxia, and anoxia, to identify the
   underlying molecular response mechanisms of plant tolerance to
   submergence, including oxygen perception and signaling. In the present
   study, we investigated the responses of red clover seedlings under
   total submersion and analyzed the root morphology, microstructure, and
   physiological response pathways. We then conducted a WGCNA analysis
   based on the obtained transcriptome data, aiming to elucidate the red
   clover response mechanisms to submergence stress in a multidimensional
   manner and to uncover the core genes and the key metabolic pathways
   involved in submergence stress adaptation and tolerance.

Materials and methods

Plant growth, submergence stress application, and sample collection

   In this experiment, the flood-tolerant variety HL (Hong Long) and
   sensitive variety ZY (Zi You) were obtained from Beijing Best Grass
   Industry. The experiment was conducted at Southwest University in
   2021–2022. Seeds from both varieties were sterilized with 10% H[2]O[2]
   for 10 min during the pre-emergence period. Afterward, they were rinsed
   five times with sterile water to disinfect them. Subsequently, the
   seeds were placed in Petri dishes on a moist double layer of filter
   paper and incubated in a germination chamber at a constant temperature
   of 24°C for 4–7 days. When the cotyledons were fully expanded,
   seedlings with vigorous and uniform growth were selected and
   transplanted into pots (15.0 cm in diameter and 13.5 cm in height)
   containing a volume ratio of nutrient soil: vermiculite: perlite of
   3:1:1. Fifteen seedlings were transplanted in each pot. The plants were
   placed in an illuminated plant growth chamber at 22/15°C (day/night),
   80% relative humidity, a 10,000 lx light setting, and a 16/8 h
   photoperiod (day/night) and watered every two days. Submergence stress
   was applied when the leaves were about 2 cm long and V-shaped white
   spots appeared on the leaf surface. Plants with the same overall growth
   were selected and placed in a water tank (length 80 cm × width 57
   cm × height 50 cm), and water(sterilized water) was added to completely
   submerge the plants in water and the water level was kept 5 cm above
   the plants.

   After the submergence stress treatment initiation, the leaves of HL and
   ZY clover seedlings were cut at 0 (CK, No submerged stress treatment),
   8, and 24 h of stress treatment (three biological replicates at each
   time point). They were then immediately snap-frozen in liquid nitrogen
   and stored in a freezer at -80℃ for subsequent transcriptome sequencing
   analysis. The roots were also sampled under the same submergence stress
   conditions at 0 h, 8 h, 24 h, 3 d, and 5 d and stored in a refrigerator
   at 4℃ for microstructure and physiological index analysis.

Histological analyses

   Roots sampled at each time point of submergence stress were scanned
   using a Microtek ScanMaker i800 plus root scanner. The obtained images
   were processed with the WinRHIZO software to obtain data on phenotypic
   traits such as root length, diameter, surface area, and number of root
   tips. Paraffin sections of red clover roots under submergence stress
   treatment were prepared with the following steps: Extraction and
   fixation; dehydration and clearing of the tissue; wax immersion and
   embedding of the tissue; sectioning, spreading; dewaxing rehydration;
   and staining. ions, the tissue sections were microscopically scanned
   and imaged using a PANNORAMIC panoramic section scanner (3DHISTECH,
   Hungary). Subsequently, observations were made at a magnification range
   from 1–400 × using CaseViewer 2.4 software. Target areas of selected
   tissues were imaged at 100 × to observe the structural changes in the
   root system under different submergence stresses of each variety.

Measurement of physiological indices

   Physiological indices of the HL and ZY clover roots were determined at
   each time point of submergence stress. The electrolyte leakage method
   was used to measure the relative conductivity of the root system. The
   root malondialdehyde (MDA) content was determined using the
   thiobarbituric acid reactive substances (TBARS) assay [[57]23].
   Antioxidant enzyme activities (peroxidase, POD; superoxide dismutase,
   SOD) were quantitatively determined using an ELISA-based protocol
   (TECAN, Tech Spark 10M, Zurich, Switzerland) [[58]24]. Three biological
   replicates were assessed per treatment.

Differentially Expressed Gene (DEG) Screening and Weighted Gene Co-expression
Network Analysis (WGCNA)

   In order to obtain transcript information under submergence stress for
   the red clover HL and ZY varieties, total RNA was extracted from
   seedling leaves of both varieties at 0, 8, and 24 h under submergence
   stress using the Trizol kit (Invitrogen, USA). The extracted RNA
   concentration and purity were assessed using Nanodrop 2000 (Thermo
   Scientific, USA). The RNA Nano 6000 assay kit for the Agilent
   Bioanalyzer 2100 system (Agilent Technologies, USA) was used to assess
   the RNA integrity to ensure that high RNA-quality samples were used for
   transcriptome sequencing. Subsequently, cDNA sequencing libraries were
   constructed from the samples that passed the quality control. The
   library concentration was accurately quantified using q-PCR. The cDNA
   libraries were pooled based on the target sequencing depth and coverage
   and then sequenced using the Illumina high-throughput sequencing
   platform based on the sequencing by synthesis (SBS) technology. Before
   data analysis, the raw data were filtered, and high-quality clean data
   were obtained after strict quality control, such as removing the
   adapter sequences from the reads and removing low-quality reads. The
   obtained high-quality clean data were aligned and mapped to the
   reference clover genome based on the HISAT2 [[59]25] alignment
   software. The reference genome is available for download at:
   [60]ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/583/005/GCA_00058300
   5.2_Tp1.0/GCA_000583005.2_Tp1.0_genomic.fna.gz. After mapping and
   assembly, the reads were quantitatively analyzed using String Tie.
   Differential expression analysis between treatment groups was performed
   using DESeq2. Fold Change ≥ 1.5 and P-value < 0.01 were used as the
   screening criteria to compare the gene expression levels in the leaves
   of seedlings of the two varieties at different periods of submergence
   stress and to identify DEGs between the comparison groups.

   Transcripts with an absolute value of FPKM > 1, obtained from the HL
   and ZY clover varieties under submergence stress treatments at 0 h, 8,
   and 24 h, were screened for differential expression. Co-expression
   networks of differentially expressed genes were constructed using the R
   software and the WGCNA package for gene expression profile matrices
   [[61]12]. Firstly, the gene network was assumed to obey a scale-free
   distribution, and the two-by-two correlation coefficients between genes
   were calculated. The Pearson correlation coefficient for each pair of
   gene i and gene j was Sij, and the gene co-expression correlation
   matrix S = [Sij] was constructed using Sij =|cor(i,j)|, where S denotes
   the similarity matrix. Then, the threshold value of the correlation
   coefficient between the specified genes was determined using the power
   exponential adjacency function aij = power (Sij, β) =| Sij |β, where
   aij is used as a measure of the relationship between gene i and gene j.
   According to the principle of scale-free networks, the parameter
   weighting factor β of the adjacency function was defined. Subsequently,
   a topological overlap dissimilarity measure (dissTOM) was performed.
   The matrix S was first converted into an adjacency matrix A = [aij],
   and subsequently, the adjacency matrix was converted into a topological
   overlap matrix (TOM), TOM = direct correlation + indirect correlation.
   In order to eliminate the error caused by pseudo-association, diss
   TOM = 1- TOM was used to calculate the dissimilarity between genes and
   obtain the dissimilarity matrix. Finally, the hierarchical clustering
   tree between each two genes was constructed based on the dissimilarity
   coefficients, and the resulting clustering tree was cut using Dynamic
   tree cut (DTC) to identify gene co-expression modules. To merge modules
   with similar expression and identify a greater number of significant
   modules, the module similarity threshold was set to 0.25, and the
   minimum number of genes per gene network module was set to 30. This
   process merged genes with similar expression patterns on the same
   branch, with each branch representing a co-expression module and each
   color corresponding to genes on the clustering tree belonging to the
   same module. Each module eigengene (ME) and the correlation of the
   eigenvectors of the resulting modules with the treatments were
   calculated. Modules with a correlation greater than 0.8 were selected
   as the modules that were significantly correlated with each treatment.
   After identifying the significantly correlated modules, Gene Ontology
   (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) annotations and
   enrichment analyses were performed for the genes in the module,
   finally, the co-expression network was visualized using Cytoscape
   (version 3.9.0). Nodes with weight > 0.5 in each module were selected
   for gene interaction network construction, and the degree was
   calculated using the Centiscape 2.2 plug-ins, which displays node sizes
   according to degree size, and the three or four genes with the largest
   degree were selected as Hub genes.

qRT-PCR validation

   Eight co-expressed DEGs were randomly selected, as well as Hub genes
   identified from WGCNA. To validate the transcriptome sequencing
   results, quantitative real-time PCR (qRT-PCR) analysis was performed on
   the selected genes using the GAPDH gene as the internal reference gene
   [[62]26]. The primers used are listed in (Table [63]S1). The
   PrimeScript ™ RT reagent Kit with the gDNA Eraser kit (Takara, Dalian,
   China) was used to reverse transcribe RNA into cDNA. qRT-PCR was
   performed with the TB Green Premix Ex TaqTM II kit (Takara, Dalian,
   China), and the reaction system is shown in (Table S2). Three technical
   replicates were performed for each target gene. PCR amplification was
   performed using a three-step method: pre-denaturation at 95°C for 30 s;
   40 cycles of 95°C for 15 s, 60°C for 1 min; melt curve analysis at 95°C
   for 15 s, 60°C for 1 min, 95°C for 15 s. Finally, the relative
   expression of genes was calculated using the 2^−ΔΔCt method [[64]27].

Statistical analysis

   Microsoft Excel 2016 and Prism 10.0 Software (GraphPad, USA) were used
   to organize and summarize the data. Statistical significance between
   treatment and control groups was analyzed by one-way analysis of
   variance (ANOVA) via SPSS 26.0 (SPSS Institute Inc, USA). Duncan's
   method was used to determine significant differences between treatments
   (P < 0.05), different letters indicate significant differences.

Results

Clover morphological features and changes under submergence stress

   Red clover plants became elongated, and both the aboveground parts and
   root system grew after 8 and 24 h of submergence stress (Fig. [65]1a).
   At 3 d of submergence stress, leaves of both HL and ZY clover varieties
   showed signs of decay, the tips of their root system gradually turned
   black, and the length of the root system was gradually reduced. At 5 d
   of submergence, the aboveground tissues of ZY were severely decayed,
   and the number of leaves was reduced due to leaf abscission. At 24 h of
   submergence, the root length of the HL and ZY varieties reached their
   maximum values of 63.10 mm and 53.55 mm, respectively. HL had
   significantly longer root than ZY and showed a highly significant
   difference at 3 and 5 d (Fig. [66]1b). Both HL and ZY exhibited a trend
   of an initially increasing and then decreasing root surface area, with
   maximum values of 7.12 mm^2 and 7.23 mm^2 at 24 h of submergence
   stress. However, no significant differences were observed between the
   varieties (Fig. [67]1c). The HL and ZY root diameters reached a maximum
   value of 0.36 mm and 0.44 mm at 24 h of submergence stress. ZY showed
   the most significant differences at 24 h and 3 d, having a 23.06% and
   13.43% greater root diameter compared to HL (Fig. [68]1d). Moreover, HL
   had a higher number of root tips compared to ZY at all time points of
   submergence stress, with significant differences between varieties at
   24 h and 3 d of submergence stress (Fig. [69]1e). Taken together, these
   results suggest that the different submergence symptoms expressed in HL
   and ZY root morphology at different time points of submergence stress
   may be due to the differences in submergence tolerance between the two
   varieties due to genotypic differences.

Fig. 1.

   [70]Fig. 1
   [71]Open in a new tab

   Effects of submergence stress on morphological indicators of red
   clover. a Morphological changes. b root length. c root surface area. d
   root diameter.(e) number of root tips

   Different lowercase letters indicate significant differences (P < 0.05)
   between different submergence stress treatment time points of the same
   variety; * indicates significant differences (P < 0.05) between
   different varieties of the same submergence stress treatment time
   points, ** indicates highly significant differences (P < 0.01) between
   different varieties of the same submergence stress treatment time
   points, the same below.

Microstructure of the root system under submergence stress

   Root anatomical properties represent the most direct manifestation of
   the level of root development and are intricately linked to
   physiological functions. In this study, it was found that submergence
   stress significantly impacted the microstructure of the red clover root
   system. Based on root transverse section observations (Fig. [72]2a),
   the increase in submergence duration resulted in an increased
   thickening of the mid-columns of both HL and ZY varieties. The number
   of conduits exhibited an increasing trend in both HL and ZY at all time
   points of submergence stress, showing highly significant differences at
   8 h and 24 h, with 50% and 15.56% more conduits in HL than ZY
   (Fig. [73]2b). However, the conduit diameter did not significantly
   differ (Fig. [74]2c). HL and ZY were significantly thickened at both 8
   and 24 h of submergence stress and reached a maximum value of 0.32 mm
   and 0.26 mm at 24 h. At 3 and 5d of submergence stress, a significant
   reduction of the cortex thickness of ZY was observed, with the most
   significant change at 5 d, with the cortex thickness measured at 0.16
   mm (Fig. [75]2d). The vascular column diameters of both HL and ZY
   showed a gradual increase with the increase in the duration of
   submergence stress, with significant differences observed at each time
   point. The vascular column diameter of HL was significantly longer
   compared to ZY at 24 h and 3 d (Fig. [76]2e). In conclusion, red clover
   may adapt to short-term submergence mainly by increasing the number of
   root conduits, cortex thickness, and vascular column diameter at the
   microstructural level.

Fig. 2.

   [77]Fig. 2
   [78]Open in a new tab

   Effects of submergence stress on the microstructure of red clover root
   systems. a Microstructure of the root system. b number of conduits. C
   conduit diameter. d cortical thickness. E diameter of vascular string

Physiological indices of the root system under submergence stress

   The relative conductivity, malondialdehyde content, peroxidase, and
   superoxide dismutase activities of HL and ZY roots exhibited a trend of
   initially increasing and then decreasing with the increase in
   submergence stress duration. Among them, the relative conductivity
   reached its maximum values at 24 h, with ZY consistently exhibiting a
   higher relative conductivity compared to HL (Fig. [79]3a), which was
   significantly higher at 5 d of submergence (P < 0.01). This indicated
   that submergence stress had a more significant effect on membrane
   permeability in ZY. Malondialdehyde content reached its maximum values
   at 24 h of submergence stress in HL. At the later time points, it was
   lower compared to ZY, indicating that the degree of membrane lipid
   peroxidation was differentially affected in the two varieties. It also
   suggested that HL was more sensitive to submergence stress in the
   initial submergence stage (Fig. [80]3b). The SOD activity in the root
   systems of HL and ZY did not exhibit a significant trend at all time
   points of submergence stress (Fig. [81]3c). On the other hand, the POD
   activity in the root systems of HL and ZY tended to increase with the
   increase of submergence stress duration and eventually decreased at
   later points in time. Specifically, at 8 h of submergence stress, the
   POD activity in the root systems of both HL and ZY was significantly
   increased, reaching its maximum value, and then showed a decreasing
   trend (Fig. [82]3d). Correlation analyses of SOD, POD, relative
   conductivity and MDA with phenotypic metrics revealed more significant
   correlations between the metrics in ZY ((Fig. [83]S1b)) compared to HL
   (Fig. [84]S1a).

Fig. 3.

   [85]Fig. 3
   [86]Open in a new tab

   Effect of submergence stress on (a) relative conductivity, (b)
   malondialdehyde, (c) superoxide dismutase(There was no significant
   correlation among treatments.) and (d) peroxidase of red clover roots

Identification of differentially expressed genes

   The raw sequencing data from 18 samples generated by Illumina's
   high-throughput sequencing platform were subjected to strict quality
   control. They yielded a total of 124.53Gb of high-quality clean data,
   with Q30 percentages of 93.58% and above. The clean reads of each
   sample were compared with the designated reference genome for sequence
   alignment, with the alignment coverage ranging from 77.54% to 83.50%
   (Table S3). Gene expression analysis was then performed based on the
   processed and aligned sequencing reads. Differentially expressed genes
   (DEGs) were identified based on their expression in different samples
   and were used for subsequent analyses.

   To obtain the DEGs between the HL and ZY varieties under submergence
   stress at 8 and 24 h, different comparison groups were designated
   compared to the control: H8h vs. 0h, H24h vs. 0h, Z8h vs. 0h, Z24h vs.
   0h. Differential expression analysis was performed using DESeq2, and
   Fold Change ≥ 1.5 and P-value < 0.01 to identify the DEGs. The
   screening criteria were used to determine the number of DEGs within
   each comparison group. The results are shown in (Fig. S2, Fig. [87]4a)
   (Table [88]1). The number of down-regulated genes was greater than the
   number of up-regulated genes in both varieties after submergence
   stress, and the number of differentially expressed genes increased
   after 24 h of submergence stress.

Fig. 4.

   [89]Fig. 4
   [90]Open in a new tab

   a Venn diagram of differentially expressed genes in each treatment
   group under submergence stress. b Soft threshold (power) selection
   diagram. c Module hierarchical clustering tree diagram. d Module gene
   clustering heat map. e Heat map of correlation between samples and
   modules

Table 1.

   Statistics of differentially expressed genes in different treatment
   comparison groups
   Variety Comparison groups Total Up    Down  Unique Common
   HL      H8h vs 0 h        5,060 2,440 2,620 1036   4024
           H24h vs 0 h       9,016 4,277 4,739 4992
   ZY      Z8h vs 0 h        5,040 2,244 2,796 1491   3549
           Z24h vs 0 h       6,570 2,922 3,648 3021
   [91]Open in a new tab

Construction and analysis of weighted gene co-expression networks

   WGCNA implements a clustering algorithm for constructing gene
   co-expression networks based on differential expression data. In this
   study, a total of 7740 transcripts with FPKM absolute values > 1 were
   selected as differentially expressed between the two varieties of red
   clover ZY and HL under submergence stress at 0, 8, and 24 h to
   construct gene co-expression matrices. Co-expression networks were
   constructed utilizing the R software and the WGCNA package. The degree
   of gene co-expression in different modules was determined, and the
   association between gene expression modules and each treatment was
   explored to identify target genes and their respective networks.

   Initially, a soft threshold was determined. To measure whether genes
   exhibited similar expression patterns and to ensure that the
   connections between genes obey the scale-free network distribution, the
   square of the correlation coefficient (R^2) was set to 0.8, above which
   the expression patterns are considered similar. Higher R^2 values
   indicate that the network approximates a scale-free distribution. In
   this study, β = 14 (R^2 > 0.8) was chosen to construct a weighted gene
   co-expression network (Fig. [92]4b). Next, module hierarchical
   clustering was performed. The hierarchical clustering tree among genes
   was constructed based on the TOM value among each gene node, and the
   gene co-expression modules were divided, as shown in (Fig. [93]4c).
   Each color in the figure represents genes on the corresponding
   clustering tree that belong to the same module, and the number of
   different colors represents the number of modules. 13 gene
   co-expression modules were eventually obtained. The number of genes in
   each module is shown in (Table S4). Finally, module gene-clustering was
   performed. The H-clust function of the R package was used to construct
   a module-gene clustering heat map for the 13 modules obtained above
   (Fig. [94]4d).

   In general, if the correlation between a module and a sample is
   significantly higher than other modules, it indicates that this
   specific module may have the strongest association with that treatment.
   Therefore, modules with module-sample correlation coefficient R^2 > 0.8
   and significance P < 0.05 were selected for analysis in this study.
   Based on this information, the module-sample correlation heat map was
   constructed, as shown in (Fig. [95]4e) below. Based on the results, the
   cyan and purple modules had very high correlations with samples H8 and
   H24, respectively (correlation coefficients were 0.96 and 0.99 with
   P-values of 0.003 and 6e-05). This suggests that genes in the cyan and
   purple modules may be involved in the response of the HL variety to
   submergence stress at 8 and 24 h. It is worth noting that the light
   cyan module was strongly correlated with sample H0 (correlation
   coefficient 0.94, P-value 0.006) and exhibited a negative correlation
   with other samples of HL as the submergence stress duration increased,
   suggesting that the genes in the light cyan module may have a negative
   regulatory effect on the response of HL to submergence stress. However,
   none of the above modules showed a significant correlation among the ZY
   samples. The ivory and light green modules in ZY exhibited significant
   correlations with Z8 and Z24, respectively (correlation coefficients of
   0.92 and 0.83, P-values of 0.009 and 0.04). Similarly, these two
   modules did not exhibit any correlation with the HL samples. This
   suggests that there may be significant differences in gene expression
   regulation between the submergence-tolerant and sensitive varieties HL
   and ZY in response to submergence stress.

GO enrichment analysis in key gene modules

   The four modules, cyan, purple, light cyan, and ivory, with the highest
   correlation coefficients with the samples and the strongest
   significance, were selected for GO enrichment analysis. The 211 DEGs in
   the cyan module were annotated into 40 GO terms, containing 15
   biological process GO terms, 14 cellular component GO terms, and 11
   molecular function GO terms (Fig. S3a). The 215 DEGs in the purple
   module were annotated into 37 GO terms, containing 15 biological
   process GO terms, 15 cellular component GO terms, and 7 molecular
   function GO terms (Fig. S3b). The 258 DEGs in the light cyan module
   were annotated into 33 GO terms, containing 12 biological process GO
   terms, 11 cellular component GO terms, and 10 molecular function GO
   terms (Fig. S3c). The 249 DEGs in the ivory module were annotated into
   30 GO terms, containing 13 biological process GO terms, 12 cellular
   component GO terms and 5 molecular function GO terms (Fig. S3d). The
   genes in the four modules with the strongest correlations with the
   samples were enriched to varying degrees in the GO categories redox
   process, enzyme activities such as oxidoreductases and hydrolases,
   endomembrane system, cell cortex, cell wall, and terpene metabolism,
   suggesting the involvement of these pathways in the responses of red
   clover to the submergence stress. Meanwhile, the light cyan module
   genes, which were significantly correlated with the H0 samples that
   corresponded to the control, were significantly enriched in
   carbohydrate metabolism, fatty acid metabolism, and lipid transport.
   This module exhibited a negative correlation with H8 and H24 samples of
   submergence stress. Therefore, we hypothesized that submergence stress
   might affect carbohydrate metabolism, fatty acid metabolism, and lipid
   transport processes in red clover.

Key gene module pathway enrichment analysis

   To further explore the functions of genes in the identified modules and
   to gain a deeper understanding of the red clover response mechanisms to
   submergence stress, we used the Cluster Profiler software to identify
   the enriched KEGG pathways among genes in cyan, purple, light cyan, and
   ivory modules at P < 0.05 (Table S5).

   There were 139 genes annotated to 52 metabolic pathways in the cyan
   module, among which the significantly enriched pathways were flavonoid
   biosynthesis, glyoxylate and dicarboxylate metabolism, fatty acid
   elongation, and carbon metabolism (Fig. [96]5a). 140 genes were
   annotated to 69 metabolic pathways in the purple module, with
   significant enrichment of pathways such as homologous recombination,
   alanine, aspartate and glutamate metabolism, and phenylpropanoid
   biosynthesis, nucleotide excision repair (Fig. [97]5b); 187 genes were
   annotated to 59 metabolic pathways in the lightcyan module, among which
   the significantly enriched pathways were glycosphingolipid
   biosynthesis, other glycan degradation, sphingolipid metabolism,
   glycosaminoglycan degradation, fatty acid degradation, fatty acid
   biosynthesis, fatty acid metabolism, starch and sucrose metabolism
   (Fig. [98]5c). Finally,106 genes were annotated to 54 metabolic
   pathways in the ivory module, which were significantly enriched in the
   carotenoid biosynthesis, sesquiterpenoid and triterpenoid biosynthesis,
   pantothenate and CoA biosynthesis pathways (Fig. [99]5d).

Fig. 5.

   [100]Fig. 5
   [101]Open in a new tab

   Bubble diagram of KEGG enrichment analysis under (a) cyan, (b) purple,
   (c) lightcyan and (d) ivory module

   The results suggest that the red clover responses to submergence stress
   might involve the pathways of flavonoid biosynthesis, glyoxylate and
   dicarboxylic acid metabolism, carbon metabolism, amino acid, and
   terpenoid biosynthesis. Notably, genes in the light cyan module,
   significantly associated with the unflooded stress H0 samples, were
   significantly enriched in fatty acid degradation, fatty acid
   metabolism, and sugar metabolism pathways. Moreover, the genes in this
   module exhibited a negative correlation with the samples corresponding
   to submergence stress. Thus, consistent with our suggestion above,
   submergence stress may have inhibited fatty acid metabolism and starch
   and sucrose biosynthesis.

Construction and visualization of key gene interaction networks

   The Cytoscape software was used to map the gene interactions in the
   target modules, and the nodes with weight > 0.5 in each module were
   selected for subsequent analysis. The degree was calculated using the
   Centiscape 2.2 plug-ins, the size of the nodes was displayed according
   to the degree size. And the first three or four genes with the largest
   degree were selected as Hub genes. The Hub genes obtained from the gene
   network interaction analysis in the cyan module were gene53873,
   gene51878, gene51609, and gene64012 (Fig. [102]6a), in the purple
   module were gene48428, gene45795, and gene39564 (Fig. [103]6b), in the
   lightcyan module were gene4283, gene38670, and gene37375
   (Fig. [104]6c), and in the ivory module were gene29857, gene11848,
   gene11315, gene59154 (Fig. [105]6d).

Fig. 6.

   [106]Fig. 6
   [107]Open in a new tab

   Map of gene interaction network in the target module. a Cyan module
   network interworking diagram. b Purple module network interworking
   diagram. c Lightcyan module network interworking diagram. d Ivory
   module network interworking diagram

   The functional annotation of the above 14 Hub genes is shown in
   (Table [108]2). The genes gene11315, gene11848, gene29857, gene51878,
   gene45795, gene4283, based on their annotation, were assigned to
   functions such as RNA degradation, glyoxylate and dicarboxylic acid
   metabolism, carbon fixation in photosynthetic organisms, carbon
   metabolism, pantothenic acid and CoA biosynthesis,
   circadian-vegetative, flavonoid biosynthesis, keratin folinic acid and
   wax biosynthesis, and degradation of other sugars across nine metabolic
   pathways. These six genes may have key functions in the response of red
   clover to submergence stress. We suggest that red clover may activate
   carbon metabolism and material translocation in the early stage of
   submergence stress. It may also respond to the hypoxic environment by
   promoting the synthesis of corresponding enzymes and flavonoids, etc.
   The plants tend to senesce with the prolongation of stress supporting
   the enrichment of RNA degradation, which is one of the key markers of
   plant senescence.

Table 2.

   Functional annotation of Hub genes in the key modules
   Module Gene ID KEGG pathway Nr annotation
   ivory gene11315 RNA degradation hypothetical protein L195_g011316
   gene11848 Glyoxylate and dicarboxylate metabolism ribulose bisphosphate
   carboxylase small chain 3A
   Carbon fixation in photosynthetic organisms
   cyan

   gene29857

   gene51878

   Carbon metabolism

   Pantothenate and CoA biosynthesis

   Circadian rhythm – plant

   Flavonoid biosynthesis

   nodulin 6

   hypothetical protein TSUD_324920

   purple

   light cyan

   gene45795

   gene4283

   Cutin, suberin, and wax biosynthesis

   Other glycan degradation

   protease inhibitor, partial

   cytosolic endo-beta-N-acetylglucosaminidase-like protein
   [109]Open in a new tab

qRT-PCR validation

   In order to verify the validity of the transcriptome sequencing data,
   we performed qRT-PCR expression analysis on the eight co-expressed DEGs
   identified. As determined by the correlation analysis for these eight
   genes, the R^2 between qRT-PCR and RNA-Seq in HL was 0.8719
   (Fig. [110]7a) and 0.7911 in ZY (Fig. [111]7b). This indicated that the
   qRT-PCR results were highly consistent with the expression levels from
   the RNA-seq sequencing results, indicating that the data obtained by
   transcriptome sequencing were highly reliable.

Fig. 7.

   [112]Fig. 7
   [113]Open in a new tab

   Correlation between RNA-Seq and qRT-PCR gene expression levels in (a)
   HL and (b) ZY. c Hub gene expression level qRT-PCR validation

   Further validation of the above six genes (gene11315, gene11848,
   gene29857, gene51878, gene45795, and gene4283) revealed that gene29857,
   gene11848, and gene51878 were expressed in all samples, and the
   expression levels of these genes in different treatment groups were
   different (Fig. [114]7c). It further suggests that genes gene29857,
   gene11848, and gene51878 may be candidates that play key roles in the
   regulation of red clover responses to submergence stress by influencing
   pantothenic acid and CoA biosynthesis, glyoxylate and dicarboxylic acid
   metabolism, carbon fixation in photosynthetic organisms, carbon
   metabolism, circadian-plant, and flavonoid biosynthetic pathways.

Discussion

   The plant root system is the main organ that absorbs water and mineral
   nutrients to maintain normal plant growth and is also the first organ
   that senses soil-related stress conditions. Its growth and development
   are determined by genetic factors, external biotic stresses, and
   abiotic stresses and their interactions [[115]19]. When plant roots
   absorb water from the soil, it is first transported through cortical
   tissues and then upward through the xylem. Some studies have shown that
   the structure and properties of cortical tissues, the thickness of the
   inner and outer cortical cell walls, and the number and diameter of
   xylem vessels are altered according to the degree of stress [[116]28].
   Submergence stress also leads to disruption of the dynamic homeostasis
   mechanism of reactive oxygen species in plants, resulting in metabolic
   disorders and electron leakage in the cells. This results in the
   excessive accumulation of superoxide anion radicals (O^2ˉ), hydroxyl
   radicals, hydrogen peroxide (H[2]O[2]), nitric oxide, and other ROS,
   causing damage to the plant membrane system and cellular oxidation
   [[117]29]. Among them, the synthesis of antioxidant compounds not only
   scavenges ROS in plants but, at the same time, synergistically
   increases the activity of related enzymes to improve plants' tolerance
   to submergence stress. Studies have shown that the increase in SOD,
   CAT, and APX activities in rice [[118]30], Nicotiana tabacum L.
   [[119]31], and Triticum aestivum L [[120]32]. under submergence stress
   plays an important role in submergence tolerance. However, as
   submergence stress conditions are prolonged or get more severe, the
   antioxidant protection system is severely damaged, leading to the
   continuous accumulation of ROS. It also leads to a gradual increase in
   the content of malondialdehyde, a cellular membrane lipid peroxidation
   product, which will trigger and further exacerbate membrane lipid
   peroxidation, causing severe damage to cell membrane function and
   eventually leading to plant death.

   Weighted gene co-expression network analysis (WGCNA) is an efficient
   clustering analysis method to identify biologically important gene
   co-expression modules and obtain candidate genes by correlating them
   with target traits [[121]33, [122]34]. In this study, based on the
   transcriptome sequencing results from leaves of the flood-tolerant "HL"
   and sensitive "ZY" red clover varieties under 0, 8, and 24 h of
   submergence stress, 7740 transcripts were selected for gene
   co-expression network construction and co-expression module
   identification. GO functional enrichment and KEGG functional annotation
   of the genes in the identified expression modules (cyan, purple, light
   cyan, and ivory) revealed that red clover mainly engaged the flavonoid
   biosynthesis, glyoxylate and dicarboxylic acid metabolism, carbon
   metabolism, fatty acid metabolism, circadian-vegetative, pantothenic
   acid and CoA biosynthesis and other metabolic pathways in response to
   submergence stress.

Roles of flavonoid biosynthetic pathways in submergence stress

   Flavonoids are an important polyphenolic secondary metabolite involved
   in plant growth and defense responses with key functions and properties
   such as antibacterial, antioxidant, and free radical scavenging
   [[123]35]. It has been demonstrated that plants can respond and adapt
   to biotic and abiotic stresses by the accumulation of compounds such as
   phenols and flavonoids [[124]36, [125]37]. A previous study [[126]38]
   on phenolic and other secondary metabolite accumulation under water
   stress in Chrysanthemum morifolium L. revealed significant differences
   in DPPH antioxidant activity, and the accumulation of phenols,
   flavonoids, and anthocyanins between control and water-deprived plants.
   The content of compounds such as quercetin and apigenin increased with
   increasing water stress as measured by high-performance liquid
   chromatography analysis. It has also been suggested that polyphenolic
   compound accumulation under water stress may be related to stress
   intensity, the plant growth stage, or plant species and that water
   submergence stress may promote the biosynthesis of compounds such as
   flavonoids [[127]39]. Studies in Vitis vinifera L. [[128]40], Triticum
   aestivum L. [[129]41], and Reaumuria soongorica L [[130]42]. revealed a
   higher accumulation of these compounds under submergence stress. The
   results of this study showed that genes in both the cyan and ivory
   modules (significantly correlated with samples H8 and Z8, respectively)
   were enriched in the flavonoid biosynthetic pathway, which was
   significantly enriched in the cyan module (P < 0.05), with the Hub gene
   gene 51878 in this module was also annotated in this pathway. However,
   the ivory module did not exhibit a significant enrichment in the
   flavonoid biosynthetic pathway. In response to submergence stress,
   flavonoid biosynthesis in the flood-tolerant variety HL was altered
   more significantly than in the sensitive variety ZY, suggesting
   genotypic differences in the flavonoid biosynthesis pathway in response
   to submergence stress in red clover [[131]43] also found that flavonoid
   content was affected by water deficit and varietal interaction effects.
   Corresponding to the physiological levels, the relative conductivity of
   HL at all time points of submergence stress was less than that of ZY,
   and the activities of SOD and POD were greater than that of ZY.
   moreover, it is more certain that there is varietal variability in the
   response of red clover to submergence stress. The phenotypic and
   molecular level analyses in this study further confirmed that, HL has
   stronger tolerance to submergence stress.

Role of carbon metabolic pathways in submergence stress

   Carbon metabolism plays a decisive role in biological growth and amino
   acid biosynthesis [[132]44, [133]45]. Iwmer [[134]46] found that during
   induced and targeted metabolic modifications in plants, cellular
   equilibrium is disrupted, carbon metabolism is inhibited, leading to
   excess carbon sources, and cells synthesize amino acids in large
   quantities or initiate other non-related metabolic pathways, resulting
   in metabolic inefficiencies Yalage et al. [[135]47] investigated the
   combined effect of carbon concentration and incubation temperature on
   the biosynthesis of VOCs was investigated using a Box-Behnken
   experimental design. In this study, the carbon metabolism pathway was
   significantly enriched in the cyan module based on KEGG enrichment
   analysis. On the other hand, the biosynthetic pathways of amino acids
   such as alanine, aspartate, and glutamate were enriched in the purple
   module. Therefore, we speculate that submergence stress may have
   disturbed the metabolic balance in HL, causing carbon accumulation and
   excess carbon sources in the pre-submergence period, thus inducing the
   synthesis of organic compounds such as amino acids to maintain the
   energy supply required for red clover growth under submergence stress.

Role of the glyoxylate pathway under submergence stress

   The glyoxylate pathway, also known as the glyoxylate cycle, is a
   process that is only present in plants and microorganisms functioning
   in lipid-to-sugar conversion. It plays an important role in plant
   growth and development and resistance to stress conditions and has been
   relatively well-studied in rice [[136]48]. Among the three catalase
   isozymes isolated from tobacco (Nicotiana tabacum), CAT1 is mainly
   responsible for scavenging H[2]O[2] produced by photorespiration, CAT2
   scavenges H[2]O[2] produced by oxidative stress, and CAT3 mainly
   scavenges H[2]O[2] produced by the glyoxylate cycle [[137]49]. In the
   present study, the glyoxylate and dicarboxylic acid metabolism pathways
   were enriched among the genes of both cyan and purple modules, and they
   were significantly enriched in the cyan module, suggesting that these
   pathways play an important role in the pre-submergence stress of HL.
   The light cyan module genes were significantly enriched in fatty acid
   degradation, fatty acid biosynthesis, fatty acid metabolism, starch,
   and sucrose metabolism. Notably, the light cyan module was strongly
   correlated with sample H0. We speculate that submergence stress may
   have induced the expression of genes such as gene11848 in the
   glyoxylate and dicarboxylic acid metabolism pathways in red clover,
   thus enhancing the conversion of lipids to sugars to adapt to
   short-term submergence stress.

   Fatty acid biosynthesis is one of the fundamental metabolic processes
   in plants, and regulating their metabolism can enhance the organism's
   resilience [[138]50, [139]51]. Unsaturated fatty acids are components
   essential for plant membrane function. It has been demonstrated that
   specific ratios of saturated and unsaturated fatty acids in membrane
   lipids are required for photosynthetic thermal stability and adaptation
   to high temperatures and that inter-organ transfer of fatty acids is
   involved in mediating temperature-induced membrane changes and
   adaptation to high temperature [[140]52, [141]53]. When the relevant
   enzymes that catalyze the biosynthesis of unsaturated fatty acids are
   mutated, the amount of unsaturated fatty acids in plants is reduced,
   leading to reduced cold resistance [[142]54]. The fatty acid
   biosynthesis and degradation pathways were significantly enriched in
   the light cyan module in this study, suggesting that fatty acid
   metabolism may be involved in red clover responses to submergence
   stress.

Conclusion

   Based on WGCNA, we comprehensively assessed the transcriptome of the
   red clover varieties ZY and HL under submergence stress at 0, 8, and 24
   h, with a total of 7740 DEGs identified between the genotypes and
   submergence stress time points. Four specific modules tightly
   associated with submergence stress adaptation of each variety were
   identified using weighted gene co-expression network analysis (cyan,
   purple, light cyan, and ivory modules). The 14 obtained Hub genes were
   functionally annotated, of which six genes, including gene51878,
   gene11315, and gene11848, were annotated to RNA degradation,
   glyoxylate, and dicarboxylic acid metabolism, carbon fixation in
   photosynthetic organisms, carbon metabolism, biosynthesis of
   pantothenic acid and CoA, circadian rhythm—plants, flavonoid
   biosynthesis, terpenoid biosynthesis, and degradation pathways of other
   sugars, providing significant insights on the core genes involved in
   red clover responses to submergence stress. In this study, the
   molecular response mechanism of red clover adaptation to submergence
   stress was revealed for the first time, and the core genes and
   metabolic pathways were identified, providing a valuable data resource
   for subsequent studies of submergence stress tolerance in red clover
   and other plants.

Supplementary Information

   [143]12870_2024_5804_MOESM1_ESM.docx^ (857.4KB, docx)

   Additional file 1: Supplementary Material 1. Fig. S1 Correlation
   analysis plot for each physiological trait of HL(a) and ZY(b). *
   Indicates significant differences between traits (P < 0.05). ** denotes
   highly significant differences between traits (P < 0.01). Supplementary
   Material 2. Fig. S2 Volcano plot of differentially expressed genes in
   each treatment group under submergence stress. (a) H8h vs 0h (b) H24h
   vs 0h (c) Z8h vs 0h (d) Z24h vs 0h. Supplementary Material 3. Fig. S3
   Histogram of GO enrichment analysis under (a)cyan, (b)purple, (c)
   lightcyan and (d) ivory module. Supplementary Material 4: Table S1
   Primers of differentially expressed genes. Supplementary Material 5.
   Table S2 PCR reaction system. Supplementary Material 6. Table S3 Basic
   statistics for quality control of sequencing reads and comparison with
   reference genomes. Supplementary Material 7. Table S4 Statistics of
   differentially expressed genes in different differential groups.
   Supplementary Material 8. Table S5 Metabolic pathway statistics for
   P < 0.05 in the target module.

Acknowledgements