Transcriptomic and metabolomic analyses reveal the potential mechanism of waterlogging resistance in cotton (Gossypium hirsutum L.)

Introduction Cotton (Gossypium hirsutum L.) is susceptible to long-term waterlogging stress; however, genomic information of cotton response mechanisms toward long days of waterlogging is quite elusive. Methods Here, we combined the transcriptome and metabolome expression level changes in cotton roots after 10 and 20 days of waterlogging stress treatment pertaining to potential resistance mechanisms in two cotton genotypes. Results and discussion Numerous adventitious roots and hypertrophic lenticels were induced in CJ1831056 and CJ1831072. Transcriptome analysis revealed 101,599 differentially expressed genes in cotton roots with higher gene expression after 20 days of stress. Reactive oxygen species (ROS) generating genes, antioxidant enzyme genes, and transcription factor genes (AP2, MYB, WRKY, and bZIP) were highly responsive to waterlogging stress among the two genotypes. Metabolomics results showed higher expressions of stress-resistant metabolites sinapyl alcohol, L-glutamic acid, galactaric acid, glucose 1-phosphate, L-valine, L-asparagine, and melibiose in CJ1831056 than CJ1831072. Differentially expressed metabolites (adenosine, galactaric acid, sinapyl alcohol, L-valine, L-asparagine, and melibiose) significantly correlated with the differentially expressed PRX52, PER1, PER64, and BGLU11 transcripts. This investigation reveals genes for targeted genetic engineering to improve waterlogging stress resistance to enhance abiotic stress regulatory mechanisms in cotton at the transcript and metabolic levels of study.


Introduction
Waterlogging is a major type of flooding problem, depending on the water level encountered by plant species (Jia et al., 2021). It primarily causes oxygen deficiency in soil in plant roots, resulting in morphological, physiological, molecular, and anatomical changes in plant tissues that affect plant growth and development. Typically, major damage from waterlogging stress is similar but speciesspecific. However, very common intrinsic disruptive substances and reactive oxygen species, such as hydrogen peroxides (H 2 O 2 ) (Xu et al., 2013) and superoxide anion radical (O ·− 2 ) (Wang et al., 2018a) and malondialdehyde (MDA) (Tian et al., 2019), usually accumulate under oxidative conditions posed by waterlogging. These substances are detrimental to plant growth and development because they reduce plant photosynthesis by causing leaf chlorosis, root death, and a reduction in chlorophyll content and chloroplast damage (Parent et al., 2008;Broughton et al., 2015;Zhang et al., 2022). In contrast, when these substances are activated during waterlogging, plants also use defensive mechanisms and counteract them by adapting to anaerobic conditions. Adaptive mechanisms are normally characterized morphologically, biochemical, anatomically, and molecularly, with the aim of discovering resistance mechanisms for plant survival. Popularly reported morphological defensive mechanisms mainly include adventitious root formation (Broughton et al., 2015), hypertrophic lenticel formation (Le Provost et al., 2016), aerenchyma air space formation (Chavez-Arias et al., 2019), radial oxygen loss spaces (Dodd et al., 2013), and formation of apical meristem (Blokhina et al., 2001). Biochemical defensive mechanisms consist of the production of enzymatic scavengers such as peroxidase (POD), superoxide dismutase (SOD), dehydroascorbate reductase (DHAR), catalase (CAT), monodehydroascorbate reductase, (MDHAR), glutathione reductase (GR), and non-enzymatic scavengers, such as proline (Zhang et al., 2016;Chavez-Arias et al., 2019;Hasanuzzaman et al., 2020;Park et al., 2020). In addition, hormone biosynthesis, such as ethylene abscisic acid, gibberellins, and auxin biosynthesis, has been reported to help in plant signaling and resistance to waterlogging stress (Coutinho et al., 2018;Zhao et al., 2018).
Cotton is a prevalent crop due to its chief source of high-quality fiber, making it an important asset to mankind. It is largely cultivated in more than 80 countries worldwide, including China, India, the USA, and Pakistan (Zahra et al., 2019). However, it is commonly known to be susceptible to waterlogging stress (Bange et al., 2004). Waterlogging stress damages cotton plant growth and development, as well as nutrient uptake (Dodd et al., 2013), resulting in yield reduction (Cao et al., 2012;Najeeb et al., 2015). Existing research on cotton tolerance mechanisms has proven that cotton adapts to waterlogging using three primary mechanisms: escape, quiescence adaptation, and self-regulation and compensation . The plant uses the escape strategy under short-term waterlogging stress by increasing adventitious root growth, producing aerenchyma, and accelerating stem elongation (Hussain et al., 2014). Self-regulation and compensation strategies include indeterminate growth habits, compensatory abilities, and the acceleration of plant growth and development after stress recovery (Sadras, 1995;Zhang et al., 2017). Finally, the quiescent adaptation strategy involves changes in the activity of protective enzyme systems, changes in hormone concentration and distribution, differential expression of genes , and hydrogen peroxide H 2 O 2 signaling (Zhao et al., 2018). In China, areas closer to the Yangtze River and the Yellow River along the middle and lower plains have high potential for supporting plantations of fast-growing plants in the future. However, these places are highly flooded (Gong et al., 2007). Moreover, in the context of China's global climate change records, flooding is a frequently occurring abiotic stress, in addition to drought, which poses yearly threats to agricultural productivity (Qian et al., 2020). Furthermore, an extreme number of reports on cotton tolerance mechanisms have centered on escape strategies, self-regulatory strategies, and compensational waterlogging stress adaptability strategies (Yin et al., 2010). However, research on quiescence adaptability, which involves more molecular studies by coupling transcriptomics and metabolomics resistance mechanisms, is minimal. In addition, there is lesser information of cotton-tolerant responses to long-term waterlogging. Therefore, the objectives of this study were to observe the morphological changes in "CJ1831056" and "CJ1831072" cotton genotypes during longterm waterlogging stress and primarily the potential tolerance mechanisms of cotton at the transcriptome and metabolome profiles in response to waterlogging stress. To attain this goal, we first revealed induced waterlogging-responsive genes through the analysis of differentially expressed genes. Subsequently, we classified differentially expressed metabolites during long-term waterlogging using metabolomics analysis. Our results can provide a system-level context for advanced studies on the potential mechanism of waterlogging resistance in cotton and further contribute to the breeding of waterlogging-tolerant lines.

Plant material
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%.

Waterlogging stress treatment and sampling
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 physiomorphological 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.

RNA-seq library preparation
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 ml of USER Enzyme (NEB, USA) was used with size-selected, adaptorligated 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.

RNA-seq analysis
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 pairedend 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.

GO and KEGG enrichment analyses
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.

Metabolite extraction, LC-MS/MS analysis, and metabolite quantification
For metabolite extraction, 20 mg of the sample was weighted onto an EP tube after grinding with freeze-drying, and 1000 ml 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 mm) 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.

Validation of RNA-seq by qRT-PCR
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-ml reaction system (including 10 ml SYBR Premix Ex Taq ™ II (2×), 2 ml cDNA, and 0.8 ml 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 -DDCt method.

Morphological variations of the two cotton genotypes to waterlogging stress
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).

Transcriptome analysis
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 singleorganism 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).

Analysis of DEGs for TFs induced by waterlogging
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.

Metabolome analysis of cotton root under waterlogging stress
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 (   Heatmap analysis of DEGs involved in ROS and antioxidant signaling. The bar represents the scale of the expression levels of each gene (log2 FPKM) in each sample, as indicated by red/blue rectangles. Red rectangles represent the high expression of genes, and blue rectangles represent the low expression.

DEGs validation by qRT-PCR
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 RNAseq results.

Morphological tolerant mechanisms in waterlogged cotton
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 . 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 . Similarly in transgenic soybean, overexpression of the AtXTH31 gene induced more adventitious roots in tolerance against flooding stress at the early seedling growth stage . Clearly, our results showed the expression of four XTH DEGs (XTH9, XTH22, XTH30, and XTH21)  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. qRT-PCR verification of some key genes with RNA-seq expression. Bars represent qRT-PCR expression, and zigzags represent FPKM expression. Each bar represents the mean ± SD. Asterisk *, ***, and **** represent the significant differences (P< 0.1 to<0.05 two-way ANOVA with Tukey's HSD post-hoc test), whereas ns indicates no significance, respectively. 56T0/72T0, 56T10/721T0, and 56T20/72T20 denote 0 days of no waterlogging, 10 days of waterlogging, and 20 days of waterlogging, respectively.

FIGURE 6
Spearman correlation analysis of DEGs in the phenylpropanoid biosynthesis pathway-accumulated POS and NEG metabolite modes. (A) Correlation between DEGs in CJ56T20 vs. CJ56T0 and metabolites in POS and NEG modes. (B) Correlation between DEGs in CJ72T10 vs. CJ72T0 and metabolites in POS and NEG modes. Each square of the heat map indicates a correlation coefficient score resulting from a Spearman correlation coefficient of 0.89. Red indicates a positive correlation, whereas blue indicates a negative correlation. An asterisk * indicates a significant correlation.

Mechanism of ROS and antioxidant tolerance in waterlogged cotton
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 2 O 2 ) 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 cysteinerich 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).

Differentially expressed TF responses in waterlogged cotton roots
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 . 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.

Responsive metabolites in waterlogged cotton root
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 antiinflammatory 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 2 O 2 (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- Molecular and morphological adaptation strategies of cotton in response to waterlogging stress. Molecular adaption denoted transcriptomic expression of genes (ADH: alcohol dehydrogenase, PDC: pyruvate decarboxylase, XTH: xyloglucanendo-transglucosylase/hydrolase, RBOHD: respiratory burst oxidase homolog protein D, CAT: catalase, SOD: superoxide dismutase, POD: peroxidase, MT1: metallothionein 1, TFs: transcription factors, IAA: indole-3-acetic acid, Eth: ethylene and ABA: abscisic acid) and metabolic expression of metabolites (organic acid, phenylpropanoid, organic oxygen, and nucleotide analogue). Morphological adaptation represents adventitious roots and hypertrophic lenticels formation. 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 Lglutamic 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).

Conclusion
The combination of metabolomic and transcriptomic analyses revealed the potential mechanism of long-term waterlogging resistance between the cotton genotypes CJ1831056 and CJ1831072. Transcriptomic analysis with RNA-seq revealed significant genes such as PER1, PRX52, PER64, ADH, PDC, MT1, XTH, and SUS to be the key genes highly contributing to waterlogging resistance and highly expressed in CJ1831056 than in CJ1831072. Moreover, innately resistant TFs, including WRKY, AP2/ERF, and MYB, were significantly induced in waterlogged CJ8131056 genotype roots, which hypothetically show signs of tolerance in response to waterlogging. Contrarily, metabolomics studies with metabolites revealed significant variations between the two genotypes after 10 to 20 days of waterlogging duration. Stresstolerant metabolites, including sinapyl alcohol, adenosine, Lglutamic acid, galactaric acid, nobiletin, and L-asparagine, were significantly expressed in CJ1831056 than in CJ1831072 in response to long-term waterlogging. Based on these findings, a schematic model of cotton adaptable mechanisms to long-term waterlogging durations (Figure 8) is presented. These analyses provide comprehensive knowledge on the molecular understanding of cotton's response to long-term waterlogging, providing effective guidance for future breeding techniques in cotton production.

Data availability statement
The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA911339.

Funding
This research was supported by the Key Research and Development Plan of Anhui province (No. 202004e11020002), and this work was partially supported by the Anhui Provincial College Students Innovation and Entrepreneurship Training Program (S202010364227). The funding bodies were not involved in the design of the study. at different positions represent the relative expression amount of the metabolites at the corresponding positions. Red indicates high expression of the substance, and blue indicates low expression.

SUPPLEMENTARY FIGURE 4
Analysis of accumulated metabolites at 20 days of waterlogging in cotton waterlogged roots. (A) Metabolites' absolute value of the fold change (waterlogging/Control) at 20 days in cotton roots. (B) The bubble plot represents a metabolic pathway. (C) Hierarchical cluster analysis. The abscissa represents the different experimental groups, the ordinate represents the comparative metabolites of the group, and the color blocks at different positions represent the relative expression amount of the metabolites at the corresponding positions. Red indicates high expression of the substance, and blue indicates low expression.