Two cotton (Gossypium hirsutum L.) lines “CJ1831056 and CJ1831072” were chosen for this study. Seed materials were sourced from the Institute of Cotton Research, Chinese Academy of Agricultural Sciences (ICRCAAS). Cotton seeds were planted in a controlled environment (greenhouse) under natural lighting and temperature (22.2°C/day and 18.2°C/night) and relative humidity ranging from 70% to 85%.

After seed emergence, seedlings were treated with a modified Hoagland nutrient solution containing 5 mM nitrate ( NO 3 − ), as described previously (Sharma et al., 2019). After 30 days of culturing seedlings, plants were screened at a 0-day treatment as a control and at 5-, 10, 15-, and 20-day treatments to determine the ideal time points for waterlogging for subsequent experiments, where pots were kept in 5-cm-deep tap water (waterlogged, WL). After waterlogging stress time-point screening, vigorous physio-morphological changes in cotton leaves, stems, and roots were observed at both 10 and 20 days of stress. Hence, control (CK) experiments at 0 days and treated groups at 10 and 20 days of waterlogging time points were used. The primary roots, lateral roots, and adventitious roots were collected from each plant, frozen separately in liquid nitrogen, and stored at −80°C.

A total of l µg of RNA per sample was used as the input material for the RNA sample preparations. Sequencing libraries were generated using the NEBNext® Ultra™ RNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations. Index codes were added to attribute sequences to each sample. mRNA was purified from total RNA using poly-T oligo-attached magnetic beads. Fragmentation was carried out using divalent cations at elevated temperatures in NEBNext First Strand Synthesis Reaction Buffer (5×). First-strand cDNA was synthesized using random hexamer primers and M-MuLV Reverse Transcriptase. To synthesize the second cDNA strand, the first cDNA strand was combined with DNA Polymerase I and RNaseH. The remaining overhangs were converted into blunt ends via exonuclease and polymerase activities. After adenylation of the 3′ ends of the DNA fragments, the NEBNext Adaptor with a hairpin loop structure was ligated to prepare for hybridization. To select cDNA fragments that were preferentially 250–300 bp in length, the library fragments were purified with the AMPure XP system (Beckman Coulter, Beverly, USA). Then, 3 μl of USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37°C for 15 min, followed by 5 min at 95°C before PCR. PCR was then performed with Phusion High-Fidelity DNA polymerase, universal PCR primers, and Index(X) Primer. Finally, PCR products were purified (AMPure XP system), and library quality was assessed using the Agilent Bioanalyzer 2100 system. The clustering of the index-coded samples was performed on a cBot cluster generation system using the TruSeq PE Cluster Kit v3-cBot-HS (Illumina), according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina NovaSeq platform, and 150-bp paired-end reads were generated.

Reference genome and gene model annotation files were downloaded directly from the Cotton Genome website. The index of the reference genome was built using Hisat2 v2.0.5, and paired-end clean reads were aligned to the reference genome using Hisat2 v2.0.5. Feature Counts v1.5.0-p3 was used to count the read numbers mapped to each gene. The fragments per kilobase of exon per million mapped (FPKM) of each gene was calculated based on the length of the gene and the read count mapped to this gene. Based on the raw count data, differential expression analysis between samples was performed using the DESeq2 R package (1.16.1). The resulting P-values were adjusted using Benjamini and Hochberg’s approach to controlling the false discovery rate. Genes with an adjusted P-value<0.05 found by DESeq2 were assigned as differentially expressed.

Gene Ontology (GO) enrichment analysis of differentially expressed genes was implemented using the cluster Profiler R package, in which gene length bias was corrected. GO terms with a corrected P value of less than 0.05 were considered significantly enriched by differentially expressed genes. KEGG is a database resource for understanding high-level functions and utilities of the biological system, such as the cell, organism, and ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used the cluster Profiler R package to test the statistical enrichment of differentially expressed genes (DEGs) in the KEGG pathways.

For metabolite extraction, 20 mg of the sample was weighted onto an EP tube after grinding with freeze-drying, and 1000 μl of extract solution (methanol: water = 3: 1, with isotopically labeled internal standard mixture) was added. Then, the samples were homogenized at 35 Hz for 4 min and sonicated for 5 min in an ice water bath. The homogenization and sonication cycles were repeated three times. The samples were then incubated at −40°C for 1 h before being centrifuged at 12,000 rpm (RCF = 13,800g, R = 8.6 cm) for 15 min at 4°C. The resulting supernatant was transferred to a fresh glass vial for analysis. The quality control (QC) sample was prepared by mixing an equal aliquot of supernatants from all samples. LC-MS/MS analyses were performed using an UHPLC system (Vanquish, Thermo Fisher Scientific) with a UPLC HSS T3 column (2.1 mm × 100 mm, 1.8 μm) coupled to a Q Exactive HF-X mass spectrometer (Orbitrap MS, Thermo). To visualize group separation and find significantly changed metabolites, supervised orthogonal projections to latent structures-discriminate analysis (OPLS-DA) were applied. Then, a sevenfold cross validation was performed to calculate the values of R2 and Q2. R2 indicates how well the variation of a variable is explained, and Q2 indicates how well a variable can be predicted. To check the robustness and predictive ability of the OPLS-DA model, 200 permutations were conducted. Afterward, the R2 and Q2 intercept values were obtained. Here, the intercept value of Q2 represents the robustness of the model, the risk of overfitting, and the reliability of the model, which will be better the smaller the intercept value. Furthermore, in the OPLS-DA, the value of variable importance in projection (VIP) of the first principal component was obtained. It summarizes the contribution of each variable to the model. The metabolites with VIP >1 and P< 0.05 (Student t test) were considered significantly changed metabolites. In addition, commercial databases, namely, KEGG (http://www.genome.jp/kegg/) and MetaboAnalyst (http://www.metaboanalyst.ca/), were used for pathway enrichment analysis.

Quantitative real-time PCR (qRT-PCR) was conducted using six genes, namely, XTH9, PER1, RAP2-3, metallothionein (MT1), RbohD, and ADH. Respective qRT-PCR primers are listed in Supplementary Table S6. Total RNA was extracted from upland cotton fiber using TRIzol Reagent. cDNA synthesis was performed using an RT reagent kit (Tiangen, China). qRT-PCR was analyzed in a 20-μl reaction system (including 10 μl SYBR Premix Ex Taq™ II (2×), 2 μl cDNA, and 0.8 μl upstream and downstream primers) and a simple procedure (50°C for 2 min; 40 cycles at 95°C for 30 s, 95°C for 5 s, and 60°C for 20 s; and a final extension at 72°C for 10 min). Three biological repeats were included for each condition, and the GhUBQ gene was used as a control. The relative expression levels were calculated using the 2^-ΔΔCt method.

Morphological changes between CJ1831056 and CJ1831072 genotypes following 0 days of no stress and then 10 and 20 days of waterlogging stress were assessed. At 10 days of waterlogging compared with the control, the two cotton genotypes displayed green leaves and no signs of damage (Figures 1A, B). Formation of hypertrophic lenticel and few adventitious roots on cotton stems and roots were observed at 10 days of stress in both genotypes (Figures 1A, B). Although cotton leaves of both genotypes turned yellow and chlorotic at 20 days of stress, formation of hypertrophic lenticel and adventitious roots was vigorous in CJ1831056 and CJ1831072 (Figures 1A, C).

Here, a total of 18 qualified libraries were established. Clean reads were obtained by removing low-quality reads. Approximately 800 million clean reads were obtained with a total of 120.07 Gb of sequence data. Clean readings at Q20 and Q30 were over 98% and 93%, respectively (Supplementary Table 1). Principal component analysis (PCA) revealed that replicated biological samples significantly clustered together. This shows that biological replicates, consistent samples, and sequencing results were reliable (Figure 1B). Pearson’s correlation coefficient (PC) analysis was performed to check the correlation between all samples (Figure 1C). Results show a higher correlation at 10 and 20 days of replicated samples compared with 0 days. Furthermore, transcriptomics differentially expressed genes revealed 101,599 genes. To effectively analyze and interpret the DEGs, a P-value< 0.05 and |log2 fold change| ≥2 were used to effectively analyze and interpret DEGs as visualized among comparisons (9) (Figure 1D). In general, CJ56T10 vs. CJ56T0 and CJ56T20 vs. CJ56T0 had 41,832 DEGs higher than 8,306 DEGs in CJ72T10 vs. CJ72T0, and CJ72T20 vs. CJ72T0 during waterlogging stress (Figure 1D). This shows that the CJ1831056 genotype has a stronger response to waterlogging than the CJ1831072 genotype at the transcription level. Moreover, 33,547 DEGs were assigned for gene regulation analysis using a Venn diagram. A total of 1,862 upregulated and 843 downregulated DEGs were shared among all four comparisons (Figure 1E). This revealed that upregulated genes were more responsive to waterlogging than downregulated genes. Hence, CJ56T10 vs. CJ56T0, CJ56T20 vs. CJ56T0, CJ72T10 vs. CJ72T0, and CJ72T20 vs. CJ72T0 comparisons were selected for further analysis. For GO analysis, the total annotated number of genes was 43,894, into three ontologies (Supplementary Table 2; Additional File 1). In the biological process category, metabolic (17,784), cellular (15,675), and single-organism processes (10,731) were recorded to have a significantly expressed number of genes. The cellular component category was mainly cell (4,972), cell parts (4,972), and organelles (3,445). Molecular function categories included binding (25,575) and catalytic activity (17,034), as highly represented (Supplementary Table 2). For subsequent analysis, the assembled gene functions were evaluated through a search against the KOG database for functional prediction and classification (Supplementary Table 2).

Pathway analysis of DEGs was conducted by using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway database. Various pathways were induced in cotton roots by waterlogging stress (Supplementary Table 3). The top enriched KEGG pathways among upregulated and downregulated DEGs in CJ1831056 and CJ1831072 are shown in Supplementary Table 3. A higher differentially expressed gene count was recorded in phenylpropanoid biosynthesis, glycolysis, gluconeogenesis, and starch and sucrose metabolism pathways (Figures 2A–D and Supplementary Table 3). The phenylpropanoid biosynthesis pathway, popularly known to regulate under hypoxia and waterlogging conditions (Zhang et al., 2017; Yu et al., 2020), revealed a higher expression of peroxidases (EC: 1.11.1.7) and beta glucosidases (EC: 3.2.1.21) (Figure 2E). This is a clear indication of ROS scavenging and activation of the plant defense system, hormone production, and cell wall lignification (Kongdin et al., 2021). In the glycolysis pathway, phosphofructokinase (EC: 2.7.1.11), alcohol dehydrogenase (EC: 1.1.1.1), and pyruvate decarboxylase-2 (EC: 4.1.1.1) genes were significantly expressed in CJ1831072 than in CJ1831056 (Figure 2F). Starch and sucrose metabolism pathways are known to be induced by waterlogging stress (Loreti et al., 2018). Similarly, starch synthase (EC: 2.4.1.21) SS4, sucrose synthase (EC: 2.4.1.13), chloroplast beta-amylase (EC: 3.2.1.2), and hexokinase (EC: 2.7.1.1) genes significantly responded to waterlogging in both cotton genotypes (Figure 2G).

Identification of ROS and antioxidant and redox-sensing mechanisms plays a key role in waterlogging tolerance mechanisms. The ROS-generating RBOH-D (107905423) gene showed higher upregulation at 20 days of waterlogging in CJ1831056 than in CJ1831072 (Figure 3). Moreover, two CAT genes (107905765 and 107905725) were more upregulated in CJ1831056 than in CJ1831072 at 20 days of waterlogging, compared to the control (Figure 3). Peroxidase DEGs played a role in antioxidant signaling under waterlogging stress in cotton. For instance, 1-Cys peroxiredoxin (107916248) and 1-Cys peroxiredoxin-like (107896588) were upregulated more after 10 and 20 days of waterlogging in both genotypes compared with the control. Four SOD genes were also induced in ROS scavenging in defense against waterlogging stress. Among DEGs involved in SOD, superoxide dismutase (FE) % chloroplastic 2C was more upregulated in CJ1831072 than in CJ1831056 at 10 days of waterlogging, whereas superoxide dismutase (MN) % mitochondrial was more upregulated in CJ1831056 than in CJ1831072. GR genes (107924091 and 107923983), DHAR genes (107925966 and 107959887), L-ascorbate oxidase (AAO) genes (107943178, 107951678, and 121209663), and L-ascorbate peroxidase (APX) genes were also induced by waterlogging in cotton roots (Figure 3).

To further understand the potential tolerant mechanism of cotton toward genes waterlogging stress, highly responsive TFs were identified and analyzed. Interestingly, a considerable number of differentially expressed TFs were recorded. For instance, TFs such as MYB DNA-binding, AP2, NB-ARC, NAM, bHLH, WRKY, LRRNT2, GRAS, and bZIP were induced at 0, 10, and 20 days (Figure 4). However, among the two cotton lines, CJ1831056 genotypes showed higher recorded TFs than CJ1831072 genotypes. Parallel tendencies were also found, as DEGs downregulated were significantly higher than upregulated DEGs (Figure 4). Furthermore, among the TF families, differentially expressed AP2 (2), bZIP (3), bHLH (1), MYB (1), NAM (1), and WRKY (2) were significantly upregulated in both genotypes. On the other hand, two AP2s, one bZIP, two MYBs, one NAM, and one WRKY gene were significantly downregulated (Table 1). However, similar genes of the same TF family showed a differential expression between the two genotypes.

For the UHPLC system, OPLS-DA and PCA were used to reduce dimensionality from metabolomics data groups generated by LC-MS/MS analyses. The ionization source of the QE platform is electrospray ionization, which comprises both positive ion modes (POS) and negative ion modes (NEG). The explanation and predictability values were measured for the first two principal components (PCs) for POS (23.3% and 20%) and NEG (30.6% and 19.7%). In general, S56T10 and S72T10 have strong relationships with quality control (QC). However, the treated groups for 20 days were more significant in POS and NEG (Figure 5A). In addition, the OPLS-DA score chart differences between all groups of samples were significant (within the 95% confidence interval), confirming that the OPLS-DA results were valid (Figure 5B). The separation of control groups from treated groups and variation within treated groups indicate clear differences in metabolite accumulation between samples. Metabolite studies also show more differentially expressed metabolites (DEM) of 171 DEM at 20 days and 52 DEM at 10 days of stress compared to 51 DEM at 0 days of no stress in both POS and NEG metabolite detection coverage (Figure 5C). As shown in Table 2, top 20 significantly accumulated metabolites were selected by their respective fold change values and OPLS-DA methods (Table 2). Zero days of waterlogging stress revealed two significantly upregulated metabolites (dibutyl phthalate 1.6 FC and equol 1.1 FC) and one downregulated metabolite (15-KETE −1.5 FC) (Table 2). At 10 days of waterlogging, phenylpropanoids and polyketides (nobiletin −3.1 FC), organoheterocyclic (urocanic acid −2.3 FC), and benzenoids (sinapyl alcohol −1.8 FC) were significantly downregulated metabolites responsive to waterlogging. In contrast, an organic acid compound, galactaric acid 3.0 FC, was significantly upregulated at 10 days (Table 2). Moreover, 20 days of waterlogging revealed highly significantly upregulated metabolites including lipids and lipid-like molecules (gamma-linolenic acid 1.1 FC and testosterone 2.0 FC), organic acids, and derivatives (L-asparagine 1.5 FC, pyrrolidonecarboxylic acid 1.3 FC, and L-valine 1.1 FC) (Table 2). Conversely, significantly repressed metabolites responsive to 20 days of waterlogging were organoheterocyclic (E,E-trichostachine −3,6FC), nucleosides, nucleotides, analogue compound (adenosine monophosphate), organic acids and derivatives (N-acetyl-L-phenylalanine −1.9 FC and citric acid −1.7 FC), organic oxygen compounds (melibiose −1.2 FC), benzenoid (carmine red −2.0 FC), and phenylpropanoids and polyketides (procyanidin B2 −1.6 FC) (Table 2).

To investigate the relationship of DEGs and differentially accumulated metabolites in cotton under long-term waterlogging, a correlation analysis was performed between DEGs in CJ56T20 vs. CJ56T0 and CJ72T10 vs. CJ72T0 and between POS and NEG highly expressed metabolites with their respective KEEG ID. The correlation of each DEG–metabolite pair was assessed at a correlation coefficient of 0.89. The details of the representative phenylpropanoid DEG–metabolite correlation analysis are listed in Additional file 2. In addition, the correlations of representative DEG–metabolites are visualized in Figure 6. In CJ56T20 vs. CJ56T0, POS metabolites correlated with DEGs; adenosine, L-asparagine, L-valine, and testosterone metabolites had a considerably negative correlation with DEGs, such as PA2, PAL1, PAL2, BGLU42, C4H, CYP98A3, and FAH1. In contrast, melibiose, positively correlated with PRX52 whereas sinapyl alcohol correlated positively with PER64 (Figure 6A). In the NEG CJ56T20 vs. CJ56T0 comparison with metabolites, galactaric acid forms a significant negative correlation with DEGs, such as CYP98A3, BGLU42, BGLU46, PAL1, PA2, PRX52, and RCI3 (Figure 6A). Conversely, hypoxanthine metabolites formed significant positive correlations with CYP98A3, BGLU42, PAL1, PRX52, and RCI3 (Figure 6A). Furthermore, in CJ72T10 vs. CJ72T0, DEGs BGLU40, BGLU46, 107937102, 107936009, 107926676, CAD9, and CAD5 were significantly positively correlated with procyanidin B2 metabolite. However, these DEGs were negatively correlated with adenosine, L-valine, L-asparagine, and testosterone metabolites (Figure 6B). Conversely, RHS19, RCI3, 121215291, 121214674, 121209489, and 121202811 DEGs significantly positively correlated with metabolites, carmine red (Figure 6B).

To verify the accuracy of the DEGs, qRT-PCR analysis was conducted to compare the consistency in the expression of some genes (Wang et al., 2022). Six differentially expressed genes, namely, ADH-class-P-like (107910431), PER1 (107916248), RbohD (107957571), RAP2-3 (107896651), MT1 (107936947), and XTH9 ((107959252), were randomly chosen (Figure 7). As shown in Figure 7, ADH-class-P-like and XTH9 genes had significant expression levels during 10 days of waterlogging stress in CJ1831056 than in CJ1831072, whereas the relative expression levels in RbohD, RAP 2-3, and MT1 were significantly higher in CJ1831072 than in CJ1831056 (Figure 7). In addition, the relative expression levels of PER1 and ADH-class-P-like were significantly higher at 20 days of waterlogging in CJ1831056 than in CJ1831072, whereas in RAP 2-3 and MT1 genes at 20 days of waterlogging, significant expressions were observed in CJ1831072 than CJ1831056 (Supplementary Table 6 and Figure 7). The expression patterns of five out six genes were consistent with RNA-seq data (Figure 7). These results highly confirm the credibility of the RNA-seq results.

Plants’ tolerance against flooding has always been extensively studied. However, the ability of plants to form adventitious roots and hypertrophic lenticel are major keys to tolerance against flooding (Else et al., 2009; Shimamura et al., 2010). Formation of adventitious roots at the hypocotyl internodes or at the stem base during waterlogging stress facilitates gas exchange and nutrition absorption (Qi et al., 2019). To a greater extent, these root tissues usually replace primary roots that die as a result of hypoxic stress, allowing normal growth and development to progress (Eysholdt-Derzsó and Sauter, 2019). In the present investigation, vigorous formation adventitious roots were found on waterlogged cotton root and stems. However, formation was numerous in the CJ1831056 genotype at 20 days of stress (Figure 1A). This illustrates that waterlogged cotton may use root and stem adventitious root formation as an escape tolerance strategy for subsequent respiration. Moreover, hormones play a vital role in adventitious root growth in plants’ tolerance against waterlogging. For example, in cucumber, both auxins and ethylene were found to induce adventitious root growth under waterlogged conditions (Qi et al., 2019). Similarly, our present study showed a higher expression of differentially expressed auxins and ethylene genes in CJ1831056 genotypes in comparison with CJ1831072 genotypes (Figure S1). We hypothesized that both hormones may have a significant role in cotton tolerance against waterlogging stress. A recent report showed that a putative cell-wall-loosening enzyme xyloglucan endotransglucosylase/hydrolase (XTH) may be involved in cell wall metabolism during flooding-induced aerenchyma and adventitious root growth. Although the XTH gene family is thought to play a major role in cell wall remodeling, it is also proven to play a significant role in plants’ tolerance against abiotic stresses such as flooding (Song et al., 2018). For instance in barley, upregulation of xyloglucan endotransglucosylase/hydrolase in roots improves waterlogging resistance (Luan et al., 2018). In cotton, GhXTH1 and GhXTH3 were upregulated in the roots of tolerant waterlogged genotypes (Zhang et al., 2020). Similarly in transgenic soybean, overexpression of the AtXTH31 gene induced more adventitious roots in tolerance against flooding stress at the early seedling growth stage (Song et al., 2018). Clearly, our results showed the expression of four XTH DEGs (XTH9, XTH22, XTH30, and XTH21) (Supplementary Table 4). We hypothesized that XTH-mediated cell wall adjustment may have a role in cotton adaptation to waterlogging stress, along with a useful gene family to develop flooding tolerance lines through molecular transgenic breeding study. Hypertrophic lenticel formation is a sign of tolerance in flooding conditions. Previous reports show that it regulates oxygen transport to the root by acting as a conduit through the cambium layer for gaseous uptake (Le Provost et al., 2016). As shown in Figure 1, 20 days of stress showed vigorous lenticel organs in the CJ1831056 cotton genotype than CJ1831072. This shows that induction of hypertrophic lenticel organs may play a vital role in cotton tolerance against waterlogging stress.

The incidence of abiotic stress activates ROS overproduction. ROS accumulation is usually regulated by the respiratory burst homolog D (RbohD) which enhances ROS production and then acts as a sensory signal for other cells to increase their ROS production (Mansoor et al., 2022). In high quantities, ROS negatively affects plant growth but, when minimized, is capable of inducing defense and hypersensitive responses and protein signaling (Qi et al., 2019; Fardus et al., 2021). In the present study, RhohD (107957571) was induced by waterlogging in both genotypes. However, it was more repressed in CJ1831056 than in CJ1821072 (Figure 3), indicating that ROS production may have played a significant role in waterlogged cotton tolerance, and moreover, CJ1831056 may have a tendency to withstand waterlogging stress by repressing ROS accumulation. A clear signal of high production of harmful ROS substance is proportional to the high production of antioxidant scavenging enzymes. In the present study, antioxidant scavenging enzymes such as PODs, DHARs, SODs, MDHAR, AAOs, CATs, MTs and APXs significantly expressed to reduce the harmful effects of ROS (Figure 3). According to existing reports, the CAT enzyme catalyzes the conversion of the ROS compound (H₂O₂) into water and oxygen (Kaushal et al., 2018). In addition, SOD is popularly known to remove ROS in anaerobic organisms based on the type of metal cofactor, and POD also catalyzes the reduction of hydrogen peroxide and alkyl hydroperoxides against oxidative damage (Nevalainen, 2010; Cheng et al., 2015). Clearly, this present study showed significant expressions of CAT and POD DEGs at 10 and 20 days in CJ1831056 than in CJ1831072, respectively (Figure 3). In addition, MnSOD and FeSOD gene expressions were significant in CJ1831056 than in CJ1831072 at 10 days of waterlogging, but the FeSOD gene was poorly induced at 20 days in both genotypes (Figure 3). We hypothesized that this could contribute to waterlogging resistance in the CJ1831072 genotype at 10 days of waterlogging (Figure 1A). Metallothionein (MT) is a small cysteine-rich protein that plays a role in ROS scavenging in biotic and abiotic stresses (Patankar et al., 2019). Although the role of MTs in waterlogging resistance is minimal, in maize, downregulation of the MT1 gene is thought to play a role in the regulation of aerenchyma formation in the root cortex (Rajhi et al., 2011). Similarly, in the present study, MT1 was highly downregulated in both genotypes, with a higher expression in the CJ1831056 at the designated waterlogging time points (Figure 3). Our findings clearly suggest that MTs may not only play a role in ROS scavenging but also determines the fate of inducible aerenchyma formation in cotton genotypes under waterlogging stress. Sucrose synthase, catalyzing sucrose to UDP-glucose and fructose, plays a fundamental role in providing an adequate sugar supply during anoxic stress, Zeng and colleagues recently demonstrated that waterlogging for 5–10 days reduced sucrose synthase activities, most likely due to changes in gene expression in starch and sucrose metabolism (Zeng et al., 2021). Contrarily, the presence of sucrose in anoxic conditions significantly alleviated meristem death and was capable of restoring root tip viability (Springer et al., 1986). Here, the starch and sucrose pathways were significantly expressed, and related genes were significantly downregulated in both genotypes (Figure 2G and Figure S2B). We speculated that the downregulation of numerous sucrose synthase genes in response to waterlogging may have played a role in root chlorosis and leaf dropping in cotton genotypes (Figure 1A).

Stress-responsive TFs play critical roles in abiotic stress responses and waterlogging stress tolerance in plants including members of the AP2/ERF, MYB, bHLH, NAC, WRKY, and bZIP families (Yoon et al., 2020; Zhang et al., 2022). Similar expressions of these TFs were discovered in this study (Figure 4). N-end pathway TFs are popularly known as oxygen-sensing mechanisms thought to trigger responses that destabilize proteins in oxygen-deprived environments (Perata and Dongen, 2011). Key N-end rule pathway DEGs, such as RAP2.12, RAP2.3, HRE2, PCO1, and PCO2, were upregulated by waterlogging and downregulated by reoxygenation in both tolerant and resistant chrysanthemum cultivars (Zhao et al., 2018). However, in the present study, all RAP2-3 and RAP2-12 DEGs were significantly upregulated at 0 days of no stress but repressed at waterlogging for 10 and 20 days (Figure S2). Based on the present results, we speculated that VII ERFs played a significant role in tolerance, but expression could be relative to the plant species. PCO activity is known as a signaling response to changes in oxygen availability. According to recent reports, during submergence, PCO activity drops, resulting in increased stabilization of ERF-VIIs (Taylor-kearney and Flashman, 2021). However, in this study, PCO DEGs were highly expressed during waterlogging stress but suppressed at no stress in CJ1831056, compared with CJ1831072. We speculated that the higher expression of PCO genes could be a result of ERF-VII instability. Among other TFs, APETALA2 (AP2) domain ABR1-like genes and ethylene-responsive transcription factor ERF071 were significantly induced. ABR1 acts as a repressor of ABA (Chen et al., 2020). However, ERF71 is believed to contribute to tolerance in anoxia stress condition by increasing anaerobic gene expression and adventitious root growth under hypoxia (Licausi et al., 2010; Eysholdt-Derzsó and Sauter, 2017). In this study, both ABR1-like and ERF71 genes were significantly upregulated (Table 1). We hypothesized that both genes may play a significant role in cotton tolerance to waterlogging in ABA repression and adventitious root growth. WRKY TFs play important roles, mainly in the innate immune system of plants. For instance, recent reports have shown that the expression of GhWRKY27 was induced by leaf senescence in early-aging cotton (Gu et al., 2019). A similar result in this study showed that the WRKY27 gene was more expressed in the CJ1831056 than CJ1831072 genotype at 10 days of waterlogging, but at 20 days of waterlogging, expression in the CJ1831056 genotype decreased as expression in the CJ1831072 genotype increased (Table 1). This shows that the WRKY27 senescence gene was significantly enhanced in CJ1831072 compared with CJ1831056, as waterlogging time increased.

Metabolomics studies in cotton genotypes, CJ1831056 and CJ1831072, significantly induced phenylpropanoid biosynthesis, purine metabolism, and galactose metabolism pathways (Supplementary Table 5). Zhang et al. (2021) explained that phenylpropanoid biosynthesis and purine metabolism are the major pathways that respond to abiotic conditions (Zhang et al., 2021). Similar to the current investigation, phenylpropanoid biosynthesis and purine metabolism pathways were highly influenced at the metabolite level of study in responses to waterlogging stress. Moreover, adenosine and sinapyl alcohol metabolites were highly accumulated in these pathways at 10 days of waterlogging (Supplementary Table 3). Sinapyl alcohol metabolite acts as an antioxidant against oxidative injury (Takahama et al., 1996). In the present study, sinapyl alcohol was upregulated in CJ1831072 but downregulated in CJ1831056 (Figure S3B). We speculated that the early accumulation of sinapyl alcohol may play a role in the phenotypic resistance observed at 10 days of stress (Figure 1B). Adenosine, produced from hypoxic and damaged tissues, reduces proinflammatory activities by facilitating anti-inflammatory activities (Ohta, 2016). Adenosine metabolite accumulation observed in this study was upregulated in CJ1831056 but downregulated in CJ1831072 (Figure S3B). First, we hypothesized that accumulation of adenosine was triggered by ROS in both genotypes at 10 days of waterlogging; however, it had high-tolerance responses in CJ1831056 than CJ1831072.

Furthermore, two metabolites, galactaric acid and nobiletin, were highly accumulated metabolites under 10 days of waterlogging, with fold change values of 3.0 and −3.1, respectively (Figures S3B, S4B). Galactaric acid has been shown to be involved in osmoprotection and stress signaling under drought (Marczak and Augustyniak, 2016). Therefore, its upregulation in CJ1831056 and downregulation in CJ1831072 may support the notion that the CJ1831056 genotype was protected from osmotic stress more than CJ1831072 (Figure S4B). Wang et al. (2018b) further reported that nobiletin phytochemicals protect Saccharomyces cerevisiae cells from the oxidative damage imposed by H₂O₂ (Wang et al., 2018b). We speculated that nobiletin metabolites also contributed to antioxidant protection at 10 days of waterlogging in cotton. However, its expression was downregulated in CJ1831056 and upregulated in CJ1831072 (Figure S4B). In response to several kinds of abiotic stressors, plants greatly increase their accumulation of free amino acids and sugar compounds. In this study, valine, leucine, and isoleucine biosynthesis and galactose metabolism pathways were the two main pathways that responded to 20 days of waterlogging in cotton (Supplementary Table 5). Amino acid accumulation is evident in enhancing stress elasticity in plants (Fardus et al., 2021). In this study, L-glutamic acid, L-asparagine, and L-valine metabolites were induced at 20 days of waterlogging (Figures S3C and S4C). Reports have confirmed that exogenous applications of glutamate and glutamic acid enhance drought and salt tolerance in radish and Lens culinaris Medik., respectively (Fardus et al., 2021). These findings suggest that L-glutamic acid may play a role by increasing stress tolerance in plants, indicating that L-glutamic acid metabolite accumulation might regulate “arginine and proline metabolism” and “glutathione” pathways in this study (Supplementary Table 5). Additionally, at 20 days of waterlogging, amino acid metabolites, such as malic acid and citric acid, were both upregulated in CJ1831072 and downregulated in CJ1831056. Other metabolites, such as vanillic acid, gamma-linolenic acid, procyanidin C1, 9-OxoODE, and procyanidin B2, also accumulated significantly at 20 days of waterlogging (Figure S4C).