Two semi-wild accessions of G. hirsutum, namely, Marie-galantie85 and Lattifolium40 that were originated from Mexico and distributed in Guadeloupe and Guatemala and a released variety CRI12 for drought tolerant, sensitive, and regular checks respectively were used for transcriptome and metabolome analysis. Marie-galantie85 was specifically used for functional validation of the genes associated with the key metabolites identified in this study through gene knockout, also referred to as Virus-Induced gene silencing. A variety of CRI12 was also used for drought tolerant, sensitive, and regular checks (Xu et al., 2020).

The seeds of these lines were sterilized and aseptically kept immersed in water for one day at 30°C before being placed for germination in an absorbent paper under 28°C of day and night temperature. They were placed in 16/8 h of alternating light and dark periods (Zhang et al., 2017). The experiment was set up in a greenhouse with three biological replications in a completely random design. The properly germinated seeds on paper rafts were transferred to a hydroponic setup of Hoagland solution (Hoagland and Arnon, 1950).

The VIGS and the wild-type seedlings were subjected to 17% PEG-6000 for PEG-induced drought at the three-leaf stage (Kirungu et al., 2020). PEG-6000 was used to simulate drought stress treatment in various crops. In Triticum aestivum L. (Faisal et al., 2019), in soybean (Basal et al., 2020), in tobacco (Yang et al., 2017), in Stevia rebaudiana (Ahmad et al., 2020), in Egyptian barley cultivars (Hellal et al., 2018), in peanut (Meher et al., 2018) and Medicago sativa L. (Zhang and Shi, 2018) used PEG-6000 induced treatment for drought stress.

The sampling parts of the materials were the roots and leaves. Sampling was done at 0, 24, and 48 h for transcriptome profiling sequence and 0, 3, 6, 9, 12, 24, and 48 for VIGS trials in three replications. A total of 54 samples were done for RNA-sequencing. The samples were then rapidly submerged in liquid nitrogen and kept at −80°C until RNA extraction was completed (Xu et al., 2020).

Phenotypic and physiological data including excised leaf water loss (ELWL), relative leaf water content (RLWC), chlorophyll content, and cell membrane stability (CMS) were recorded from three randomly sampled plants. We took leaf tissue samples for the determination of biochemical analysis. In the biochemical analysis, the antioxidant catalase (CAT), superoxide dismutase (SOD), hydrogen peroxide (H₂O₂), and malondialdehyde (MDA) enzyme activities were measured. All the antioxidants and oxidant enzymes CAT, SOD, H₂O₂, and MDA were evaluated using the Solarbio life sciences kit (www.solarbio.com) using the instructions of the manufacturer. For physiological traits, the following formulas were used to calculate the values.

Note, FW = Fresh weight, SW = Saturated weight, WW = Wilted weight, DW = Dry weight, T1 = Reading 1, T2 = Reading 2, C1 = Control 1, C2 = Control 2.

To isolate RNA samples from both leaf and root tissues for RNA-Seq analysis, Trizol R reagent (Invitrogen, Waltham, MA, USA) was used. Using an agarose gel electrophoresis and a NanoDrop 2000 spectrophotometer, the quality and concentration of RNA were assessed (Thermo Fisher Scientific, USA) (Yu et al., 2017). The Agilent 2100 Bioanalyzer RNA Nanochip was used to complete a more precise RNA quantification (Agilent Technologies, Waldbronn, Germany) (Shi et al., 2018). The cDNA fragments were then cleaned using a Qia Quick PCR extraction kit that was corrected at the ends with added poly(A) and ligated to Illumina sequencing adapters. Gene Denovo Biotechnology Co. used the Illumina HiSeqTM 2500 to amplify and sequence PCR products. The size of the ligation products was determined using agarose gel electrophoresis from Gene Denovo Biotechnology Co (Guangzhou, China). PCR amplified and sequenced the PCR using the Illumina HiSeqTM 2500. Tophat2 (V.2.0.13) was used to extract clean reads from raw reads, which were then associated to the reference genome of G. hirsutum (http://cottonfgd.org) to generate mapped reads (Kim et al., 2013). Cufflinks v2.2.1 was used to analyze gene expression and to calculate the differences between the treatment and control fragments per kilobase of transcript per million fragments (FPKM) values. Readings from the transcriptome were measured and adjusted to fragments per kilobase of transcript per million fragments linked to the reference genome (Trapnell et al., 2014).

R software (http://www.r-project.org/) was used to find differentially expressed genes between samples. The false discovery rate (FDR) and log2 Fold Change thresholds for differentially expressed genes (DEGs) were set at < .05 and > 1, respectively. The DEGs were then submitted to a GO function and KEGG pathway enrichment analysis. First, through the AgriGO web tool (http://systemsbiology.cau.edu.cn/agriGOv2/), all DEGs were mapped to GO (Gene and Consortium, 2000). Pathway enrichment analysis revealed considerably enhanced metabolic transduction pathways when DEGs were compared to the entire genomic background.

To transfer genes into co-expressed modules, weighted gene co-expression network analysis (WGCNA) was done in the R package (Zhang and Horvath, 2005; Langfelder and Horvath, 2008). Before generating an adjacency matrix, the FPKM values were normalized. WGCNA software was used to import the phenotypic data, and a correlation-based relationship between accessions, time points, and gene modules was done using the default settings. The WCGNA software was used to transform the adjacency matrix into TOM (topological overlap matrix). Transcripts with comparable expression patterns were clustered into one module once the network was created, and eigengenes were determined for these modules. Using Cytoscape's default parameters, each module's genes were exported.

Leaf and root tissues of the three semi-wild cotton species Marie-galantie85, Lattifolium40, and Upland cotton at 0-, 24-, and 48-h time points were used as phenotype data in the WGCNA analysis. Thirteen unique gene modules were recognized from RNA Seq data by entering WGCNA using the FPKM values. Coexpression network analysis of the three highly and positively correlated to drought modules namely brown, yellow, and pink was established by Cytoscape v.3.7.2 (Shannon et al., 2003). Cotton functional genomics database website (www.cottonfgd.org/) was used to download the protein sequence of the candidate gene for G. hirsutum and G. arboreum, while G.raimondii, A. thaliana, Theobroma cacao, Brassica rapa, Glycine max, Medicago truncatula, Oryza sativa, Triticum aestivum, Zea mays, and Sorghum bicolor were downloaded from phytozome (https://phytozome.jgi.doe.gov). The ClustalX tool was used to align the full-length ALDH protein sequences (Larkin et al., 2007) and MEGA 7 was used to create the phylogenetic tree using the neighbor-joining method with 1000 bootstrap replications (Tamura et al., 2011).

The tender leaves and roots of cotton at the three-leaf seedling stage were used for metabolite profiling by the liquid chromatography-mass spectrometry (LC-MS) approach. To detect reproducibility under the same treatment, the sample extracts were mixed to prepare a quality control sample (QC). During the LC-MS analysis, QC samples are arranged by mixing sample extracts to analyze the recurrence ability of samples under the same processing method. In the process of instrument analysis, a QC sample is usually inserted into every 10 test analysis samples to investigate the repeatability of the analysis process. Three biological replications were maintained in sampling and metabolome analysis (Zhou et al., 2012).

RNA sequencing of G. hirsutum yielded a substantially upregulated gene, which was then transformed into Arabidopsis thaliana (Colombia-0). The PCR study with forward (CGGCCATTTAAATAGTGGATTCGG) and reverse (GCCACCATCTATCCTCAACGA) pair of primer sequences of Gh_A06G1257 synthesized from Invitrogen, Beijing, China, verified the pWM101-35S: Gh_A06G1257 build in Agrobacterium tumefaciens GV3101. The floral dip approach was used to transform wild-type A. thaliana plants.

The infiltration media was prepared as prescribed previously (Lu et al., 2019). The seedlings were moved to a growing environment with a temperature of 25°C and a 16 light/8 dark hours cycle from the selection medium at three-leaf stages. The little plastic pots were filled with a 1:1 mix of vermiculite and humus (Sadau et al., 2021).

The first-generation seeds were collected after the seedlings from generation T0 had been grown to set seeds (T1). True lines were found by estimating the antibiotic-selectable marker segregation ratio of 3:1 after T1 seeds were sown in antibiotic media. Only the lines with a 100% success rate were chosen for T3 generation growth. After RT-qPCR, T3 homozygous lines were chosen from a T2 generation. Three of the six ALDH transgenic lines that were successfully transformed (OE-1, OE-9, and OE-10) were chosen from a T2 generation. The RT-qPCR prototype was Gh_A06G1257 forward primer sequence and Gh_A06G1257 reverse primer sequence with complete complementary DNA (cDNA). T3 homozygous generation was used to conduct the phenotypic studies.

The plants were subjected to drought stress treatments after 21 days of growth. Samples from leaves were collected from transgenic lines and control. The pots were then watered with water containing 15% PEG-6000 treatment was used to alleviate drought tension. After 8 days, physiological and phenotypic traits were observed. All measurements were replicated three times biologically and three times technically (Li et al., 2019).

Drought simulated stress conditions were used to assess the germination percentage and root length of transgenic lines and wild form. OE-1, OE-9, OE-10, and wild-type seeds were sown in.5 MS plates complemented with 0, 100, 200, and 300 mM mannitol concentrations to simulate drought. After 10 days, the germination rate was assessed. Transgenic and wild-type seeds were seeded for 6 days in 0.5 MS media before being moved to 0.5 MS supplemented with different amounts of 0, 100, 200, and 300 mM mannitol for the root length assay. All measurements were replicated three times biologically and three times technically (Li et al., 2019).

The functional characterization of substantially upregulated gene expression in G. hirsutum was investigated using this approach. A gene Gh_A06G1257 with 296 bp fragments was transformed by Forward (CTGTGAGTAAGGTTACCGAATTCTCTAGAGAGATGTGGAATCCCCTTGGAA) and Reverse (TCGAGACGCGTGAGCTCGGTACCGGATCCACCTTCGAACTCCCCGTGA) primer sequences into pTRV vector using enzymes XbaI and BamHI to develop a 35S promoter-driven pTRV2: ALDH7B4. The promoter cells of A. tumefaciens LBA4404 were deformed with the recombinant vector by freezing and thawing (Velásquez et al., 2009). Wild, pTRV: 00 were used positively with Phytoene desaturase (PDS) as a negative control.

To validate RNA-Seq data, DEGs were verified by RT-qPCR using (Schmittgen and Livak, 2008) procedures. Primers were created using NCBI (https://www.ncbi.nlm.nih.gov) and gene sequences from G. hirsutum. Quantitative PCR was performed using an SYBR® Green PCR Master Mix Kit (Applied Biosystems, Foster City, CA, USA) and an ABI-7900 system. The 2-^ΔΔCT was used to determine the relative gene expression (Schmittgen and Livak, 2008). Internally, GhActin was employed as a reference. The RT-qPCR experiment was carried out in three biological and technical replications (Chen et al., 2015).

Chenomx NMR Suite 7.7 was used to do the multivariate statistical analysis, principal component analysis (PCA), partial least squares determinant analysis (PLS-DA), and orthogonal projections to latent structures (OPLS-DA). The variance in the data matrix was presented using PCA. PLS-DA can maximize the differentiation between groups, which helps find differential metabolites. OPLS-DA data was used to examine subsequent model tests and differential metabolite screening. Fisher's exact test in R was used to calculate statistical significance. The Benjamini-Hochberg correction was used to adjust the false recovery rate (FDR) in the transcriptome analysis (Noble, 2009). Coexpression network analysis was established by Cytoscape v.3.7.2 (Shannon et al., 2003). The SPSS program was used to perform both t-tests and ANOVA. Significant differences were declared at the p < 0.05 probability level (Worley and Powers, 2015).

Wilting was observed in the plants under drought stress with reduced evaporation. There were no obvious differences in the phenotype of three semi-wild cotton lines in response to drought stress (Figure 1I). Clear and significant differences (p < 0.05) were recorded in the case of physiological and biochemical parameters (Figure 1II). Chlorophyll and relative water contents in MG85 and CRI12 were higher as compared to LT40. A lower ion leakage with a higher excised leaf water loss was observed in LT40 as compared to CRI12 and MG85 (Figure 1III).

Higher activities of CAT and SOD were observed in MG85 and CRI12 under drought stress conditions, however lower CAT and SOD were observed in the case of LT40. Moreover, we also measured the H₂O₂ and MDA activities. We found out that they were higher in LT40 as compared to MG85 and CRI12. Thus, suggesting that MG85 and CRI12 were more resistant to drought as compared to LT40 (Figure 1IV).

We selected semi-wild cotton lines, namely, Marie-galantie85, Lattifolium40, and Upland cotton for the transcriptome sequencing in response to drought treatment. After sequencing, 54 transcriptome libraries were obtained. For the accuracy of the subsequent analysis of the new data, reads with lower, linker sequences, sequences having >10% unknown N were deleted. After filtering the new data, we finally acquired high-quality clean reads of 2.66 Gigabytes in all the 54 libraries (Table 1). The average length of each clean-read was 150 bp. We, therefore, acquired over 95% high-quality clean reads. The Q20 and Q30 base percentages were all above 96 and 90%, respectively, and the GC content was higher than 43% (Table 1). Bowtie2 was used for the screening of high-quality reads and then comparing them with G. hirsutum reference genome by TopHat2 software. The samples were tested for Pairwise repeatability and the correlation coefficient was all above 84%, indicating that the samples were reproducible.

We found that the upregulated genes were higher in number as compared to the downregulated genes in both the leaves and roots. Specifically, TRL0 vs. TRL24 (23701 vs. 12610), TGL0 vs. TGL24 (21903 vs. 1118), and TSL0 vs. TSL48 (11618 vs. 9632) in leaves. Whereas TRR0 vs. TRR48 (19935 vs. 16030), TGR0 vs. TGR48 (19598 vs. 14874), and TSR0 vs. TSR48 (19031 vs. 13490) were highly upregulated and downregulated genes in roots (Figure 2II).

The reason for higher DGEs number in roots at different time points of drought application could be that they are the first to be in contact with PEG-6000, thus triggering the drought response. The amount of differently expressed genes in root tissues grew over time, owing to the roots' direct interaction with PEG-6000, which resulted in more sustained drought response in the roots. Intersection and union analysis of differentially expressed genes at different time points in the same material was performed by applying Venny online software (Figure 2I).

The expression patterns of genes change with conditions such as specific environments and time. The trend analysis of genes that were differentially expressed was performed by STEM software for the expression of leaf and root regulatory genes at different time points. The differentially expressed genes of the experimental materials were divided into 8 gene expression trends. According to the standard of (p < 0.05), the gene expression trend showed the same significant trend (trend 1, trend 6, and trend 7) with time in the three materials in both leaves and roots (Table 2).

According to the principal component analysis of the six quality control samples gathered together, the dispersion is very small, indicating that the stability of the instrument is better in the analysis of metabolites, and the results of metabolite analysis are reliable. Thirty-six samples were divided into four components by the main components PCA1 (39.9%) and PCA2 (12.8%). The main component PCA1 (39.9%) separated the leaves and root samples, indicating that under drought stress, leaves and roots respond differently (Figure 2III). The main component PCA2 (12.8%) separates the samples at different treatment time points of the material, and the difference in the roots is small. It is indicated that the changes in leaves are more obvious at different treatment time points of drought stress. This may be due to the sensitivity of leaf tissues to stress as compared to root tissues. LC-MS analysis was helpful to identify 445 metabolites totally (Supplementary Table 4). According to the analysis of the differentially expressed metabolites, the highest upregulated metabolites were recorded in the leaves at MGL0 vs. MGL48 (67 metabolites), MRL0 vs. MRL48 (61 metabolites), and MSL0 vs. MSL48 (56 metabolites), respectively, whereas the highest downregulation of metabolites was found in roots from MRR0 vs. MRR48 (54 metabolites), MSR0 vs. MRR48 (47 metabolites) and MSR0 vs. MSR48 (36 metabolites) successively (Figure 2IV).

The identified differentially expressed genes in both the roots and leaves tissues were mainly enriched in the molecular functions, cell components, and the biological processes of functional categories (Supplementary Figure 1). In MG85, the leaf tissues were enriched in biological and molecular functions, whereas the root tissues were enriched only in biological processes. In CRI12, both tissues are enriched in two functions, namely, biological and molecular, and biological and cellular, respectively. Differently from the above lines LT40 enriched in all the three GO functions biological, cellular, and molecular in leaf tissues only, there is no significant enrichment in the root tissues (Supplementary Figure 1). In the biological process, the most important DEGs were involved in metabolic processes (GO:0008152), cellular processes (GO:0009987), biogenesis (GO:0071840), and single organism processes (GO:0044699), etc. On the other hand, in cellular components, DEGs were mainly involved in cells (GO:0005623), cell parts (GO:0044464), organelles (GO:0043226), macromolecular complex (GO:0032991), and extracellular region (GO:0005576). The molecular functions mainly included binding (GO:0005488) and catalytic activity (GO:0003824). Among them, we found that most of the DEGs were enriched in the metabolic processes in the tissues of different materials. Thus revealing that primary metabolites and secondary metabolites play an essential part in response to drought stress.

Here we used the FPKM values of commonly expressed DEGs to perform a weighted gene coexpression network analysis for the identification of hub genes associated with drought stress tolerance. A cluster dendograms was generated to see the number of modules along with DEGs, a network heatmap was also generated from the genes located in each of the identified module (Figures 3A,B). Using WGCNA, thirteen different modules were detected. Out of these thirteen modules, three modules (Yellow, brown, and pink) were found to be highly and positively associated with our phenotype. The yellow module contains 1,568 genes, the brown module contains 1,234 genes whereas 765 genes were present in the pink module. The yellow module has significant associations with MG85 0 h, CRI12 0 h, and LT40 0 h in root with r = 0.02 and p = 0.54, respectively. Brown modules have significant correlations with MG85 0 h, CRI12 0 h, and LT40 0 h in leaf with p = 0.02, p = 0.02 and p = 0.03 with r² values of 0.55, 0.55, and 0.52, successively. Pink modules have significant associations with MG85 0 h and CRI12 0 h with r² = 0.42, 0.44 and p = 0.49 and 0.42, respectively (Figure 3C). Cytoscape v.3.7.2 was used for the network visualization.

WGCNA helped us to import the top 30 genes for each of the significant modules to build the coexpression networks. The gene lists were further used for the network visualizations. We used Cytoscape software for this purpose. To identify the hub genes, a built-in Cytoscape extension by the name of “cytohubba” was used (Shannon et al., 2003). Five hub genes from the yellow module include Gh_A07G0563, Gh_D05G0221, Gh_A05G3716, Gh_D12G1438, and Gh_D05G0697, the brown module contains three hub genes (Gh_A06G1257, Gh_A06G1256, and Gh_D06G1578), and the pink module has five hub genes, namely, Gh_A02G1616, Gh_D12G2599, Gh_D07G2232, Gh_A02G0527, and Gh_D07G0629. The most interesting thing is that all the thirteen hub genes belong to the same ALDH family. Furthermore, the expression profiles of all hub genes in different plant parts under drought stress were validated via RT-qPCR. Results from RT-qPCR suggest that the Gh_A06G1257 gene, having the highest expression under drought stress, might be the true candidate responsible for drought stress tolerance. We further validate this gene on a functional basis (Figure 4).

Two distinct clusters were formed as a result of phylogenetic analysis. The candidate gene Gh_A06G1257 together with the two Gossypium species together with A. thaliana, T. cacao, and B. rapa were grouped in one cluster. G. max, M. truncatula, O. sativa, T. aestivum, Z. mays, and S. bicolor were also grouped in the second cluster.

We performed PCR and RT-qPCR analysis to see the expression level of six positive lines. OE-1, OE-9, and OE-10 overexpressed lines conformed with the highest expression level (Figure 8I). Seedlings of overexpressed lines were exposed to drought for studying the functions of Gh_A06G1257 during drought stress. The wild-type and transgenic seedlings grew and germinated uniformly in the MS medium. On the other hand, the MS medium containing 300 mM mannitol seedlings of OE lines (OE-1, OE-9, and OE-10) showed better germination, grew faster than wild-type seedlings, and demonstrated a drought-resistant phenotype. Increased mannitol concentration resulted in a steady decrease in the rate of germination. In normal conditions, germination was calculated to be >90% in both wild and transgenic lines Figures 5A,B. However, under drought stress treatments, wild-type seedlings showed a decrease in germination rate to <25%, while a significantly higher germination rate of 60% was observed in transgenic lines in 300 mM mannitol. Under mannitol treatment, longer roots (13 mm) were observed in the lines with over-expressed ALDH gene while in wildtype a shorter root length of 6 mm, suggesting overexpressed lines are more resistant to drought stress (Figures 5C,D).

In response to drought stress, we compared the physiological responses of the three transgenic lines with the wildtype. Under normal conditions, there were no significant differences between the transgenic lines and the wild type in any of the measurements, but when the plants were stressed, the transgenic lines showed greater stress tolerance than the wild type. RLWC was much higher in transgenic lines than in control lines (Figure 6I). Wild-type plants displayed more stress-induced ion leakage when they were stressed by drought. Under regulated conditions, the ELWL did not differ statistically between wild-type and transgenic plants. However, under stress conditions, the wild-type leaves lost more water than the transgenic plant leaves Figure 6II. The fact that the transformed gene improved drought stress tolerance in the transgenic Arabidopsis lines was demonstrated by the increased degree of tolerance between the transgenic lines.

Activities of oxidants and antioxidants were measured in the transgenic and wild-type plants after applying drought treatment. Under control circumstances, there were no significant variations in SOD and CAT activity between transgenic lines and the wild type. However, when seedlings were exposed to drought stress after 8 days, the transgenic lines OE-1, OE-9, and OE-11 showed significant changes from the wild type, with increased SOD and CAT activity, lower MDA, and lower H₂O₂ levels (Figures 7A–D). Our results suggest that Gh_A06G1257 is playing a critical role in drought stress tolerance.

We choose four abiotic responsive genes ABF4, SOS1, RAB18, and RD22 to measure the expressions by qRT-PCR in over-expressed and wild-type seedlings for a better understanding of the role Gh_A06G1257 is playing when plants are exposed to abiotic stress. Results from the RT-qPCR results revealed that upon exposure of transgenic lines to drought, the expressions of ABF4, SOS1, RAB18, and RD22 were recorded higher in the transgenic lines (Figure 8II). Thus, suggesting that Gh_A06G1257 has a key role in coping with drought stress.

Significant variations in the physio-morphological traits of cotton were observed under drought conditions. The VIGS plants, positive controlled plants, and the wild types significantly (p < 0.05) differs in their response to drought stress in both morphological and physiological traits (Figure 9). Plant height and root length showed no significant differences, while VIGS plants had a considerably lower shoot and root fresh weights than wild-type pes and positive controls. Moreover, VIGS plants have a higher excised leaf water loss, and ion leakage but significantly lower relative leaf water contents than positive control and wild-type plants under drought stress.

We observed significant differences (p < 0.001) in the oxidants and antioxidants activities among VIGS and wild-type plants. Overall, lower antioxidant activities and higher oxidant activities were measured in the case of VIGS plants (Figure 10II). CAT and SOD contents before treatment were almost the same, while after being treated with PEG-6000, the wild plants showed an increase over the VIGS plants. An increase in the contents of the oxidant enzymes activities showed that the VIGS plants experience higher oxidative stress than wild-type and positive controls. Activities of MDA and H₂O₂ were similar in normal conditions but recorded higher under drought stress mainly in the VIGS plants. Thus, our results suggest that the knockdown of the Gh_A06G1257 gene has a significant influence on coping with drought stress.

Samples from the leaves and roots of wild-type, positive control, and VIGS plants were collected. To clarify the role of the ALDH gene (Gh_A06G1257) under drought, we collected samples from different tissues to check the expression. The expression of Gh_A06G1257 in VIGS plants was lower as compared to the wild-type and the positive control plants. Thus, indicating that Gh_A06G1257 is playing a critical role in coping with drought (Figure 10IIE,F).

In the KEGG pathway annotation, the metabolites were enriched in metabolism, genetic information processing, and environmental processing information groups. A total of 25 genes playing key roles in different pathways have been identified from metabolome KEGG analysis (Supplementary Table 1B). Within this, the ALDH genes were enriched in amino acids metabolism, carbohydrates and lipids metabolism, and the replication and repair activities. The metabolism enrichment analysis found that the major enriched pathways during drought stress were arginine and proline metabolisms, tryptophan metabolism, valine and leucine degradation, and lysine degradation (Figure 11). In general, higher numbers of metabolites were upregulated in the leaves than roots under drought. In arginine and proline metabolism pathway out of nine metabolites agmatine and spermidine, in tryptophan metabolism out of ten metabolites, tryptamine, N-Acetyl-5-hydroxytryptamine, and kynurenic acid, in valine and leucine degradation mainly L-valine showed upregulation in leaf tissue. On the contrary, L-ascorbic acid in glycine, serine, and threonine metabolism, glutaric acid and N6-Acetyl-L-Lysine in lysine degradation, L-valine in valine and leucine degradation, and L-proline in arginine and proline metabolism showed upregulation in root tissues during drought stress. The rest of the metabolites showed a similar trend in both tissues.