Genome-Wide Identi ﬁ cation and Expression Analysis Elucidates the Potential Role of PFK Gene Family in Drought Stress Tolerance and Sugar Metabolism in Cotton

Drought has been identi ﬁ ed as a major threat for global crop production worldwide. Phosphofructokinase (PFK) is vital for sugar metabolism. During phosphorylation, plants have two enzymes: ATP-dependent phosphofructokinase (PFK) and pyrophosphate-dependent fructose-6-phosphate phosphotransferase (PFP). Genome-wide identi ﬁ cation led to the identi ﬁ cation of 80 PFK genes, 26 genes in G. hirsutum and G. barbadense , and 14 genes in G. arboreum and G. raimondii . Phylogenetic, gene structure, and motif analyses showed that PFK genes were grouped into two main categories, namely, PFK and PFP, with 18 and 8 genes in the allotetraploid species and 10 PFK and 4 PFP genes in the diploid species, respectively. Using the RNA-seq expressions of 26 genes from GhPFK , a co-expression network analysis was performed to identify the hub genes. GhPFK04 , GhPFK05, GhPFK09 , GhPFK11 , GhPFK13 , GhPFK14 , and GhPFK17 in leaves and GhPFK02, GhPFK09 , GhPFK11 , GhPFK15 , GhPFK16, and GhPFK17 in root tissues were found as hub genes. RT-qPCR analysis validated the expressions of identi ﬁ ed hub genes. Interestingly, GhPFK11 and GhPFK17 were identi ﬁ ed as common hub genes, and these might be the true candidate genes involved in the drought stress tolerance. In the KEGG enrichment analysis, amino acids such as L-valine, L-histidine, L-glutamine, L-serine, L-homoserine, L-methionine, L-cysteine, and gluconic acid were signi ﬁ cantly upregulated, whereas sugars, mainly fructose-1-phosphate, D-mannitol, D-sorbitol, dulcitol, and lactose, were signi ﬁ cantly downregulated during drought stress. Genome-wide analysis paves the way for a deeper understanding of the PFK genes and establishes the groundwork for future research into PFK ’ s role in enhancing drought stress tolerance and sugar metabolism in cotton. of PFK genes under drought stress. Drought stress induces biosynthesis and upregulation of amino acids production and activates the plant osmo regulation activity. Additionally, PFK involves in sugar metabolism and acts as the energy source for root development during plant drought stress tolerance.


INTRODUCTION
Cotton is an appealing model for investigating polyploid origins, evolution, and domestication. Diversity in Gossypium increased by transoceanic, long-distance dispersal and broad hybridization among lineages that are currently geographically separated. Gossypium hirsutum L. (AD1) and Gossypium barbadense L. (AD2) are two cultivated tetraploids that resulted through transoceanic hybridization . They are developed separately and are domesticated in different parts of the world. G. hirsutum, which has been domesticated, shows greater adaptability and high production, whereas G. barbadense produces uniquely high-quality fibers (Wang et al., 2019). G. hirsutum is the main economical essential fiber crop, the primary source of renewable textile fibers, rich in protein and oilseed production. It was one of the first genetically modified crops to be widely used, and human-mediated breeding has resulted in current upland varieties with increased production and fiber quality (Su et al., 2016;Nazir et al., 2020).
Drought stress is a frequent abiotic stress that restricts crop growth and yield around the globe (Yao and Wu, 2016). In plants, drought triggers a complex set of molecular responses that starts with stress detection, progresses through a signal transduction cascade, and ends with physio-morphological changes at the cellular level, such as stomatal closure, cellular respiration activation, and inhibition of cell growth and photosynthesis (Lamaoui et al., 2018;dos Santos et al., 2022). Plants can also generate and accumulate certain metabolites that are specifically engaged in stress tolerance when they are subjected to drought (da Silva et al., 2017).
Glycolysis is a metabolic activity in which enzymes convert glucose into two molecules of pyruvate in the availability of oxygen or two molecules of lactate without oxygen. Anaerobic glycolysis, the latter pathway, is thought to be the first natural method to create adenosine triphosphate (ATP) (Plaxton, 1996). Because mitochondria are absent in few cells, like mature red blood cells, glycolysis is the sole way to produce ATP. It is a cytoplasmic method for converting glucose into two three-carbon molecules, which generates energy. Phosphorylation of glucose with the enzyme hexokinase traps glucose. All cells in the body use it to generate energy. In aerobic situations, glycolysis creates pyruvate, while in anaerobic ones, it produces lactate. To produce extra energy, pyruvate enters the Krebs cycle (Givan, 1999;Ohlendieck, 2010).
The glycolysis and pentose phosphate pathways have a significant impact on crop tolerance and abiotic stress. Phosphofructokinase (PFK) is a range-limiting enzyme in the carbon-flow controlling pathways of biological activities (Yao and Wu, 2016). The phosphorylation of fructose-6-phosphate to fructose-1,6-bisphosphate, a crucial regulatory step, is catalyzed by PFK, a glycolytic enzyme. The transfer of a phosphoryl group from ATP, which is mediated by enzymes, is a vital step in numerous biological activities (Ros and Schulze, 2013).
In crops, phosphofructokinase is vital for sugar metabolism. During phosphorylation, plants have two enzymes such as: ATPdependent phosphofructokinase (PFK) and pyrophosphate-dependent fructose-6-phosphate phosphotransferase (PFP) (Lü et al., 2019). PFK gene families have been characterized and a number of genes have been identified in some crops, such as eleven genes in Arabidopsis thaliana (Mustroph et al., 2007), fifteen PFK genes in rice (Kato-Noguchi, 2002), fourteen genes in white pear (Lü et al., 2019), thirteen genes in cassava , and studies were also performed in spinach and Saccharum (Chen et al., 2017).
The sugar required for root growth and metabolism is transported from leaves, enhancing the sucrose transport from leaves to roots is conducive to maintaining root growth under drought stress (Xu et al., 2015). In line with this, a higher root/ shoot (R/S) ratio was pronounced under drought stress in soybeans. It increased the contents of soluble sugar and sucrose in the leaves, but decreased starch content; in the roots, all of these parameters were increased. This may be related to the enhanced carbohydrate metabolism activity under drought stress, including notable changes in the activities of sugar metabolism enzymes and the expression levels of genes in soybeans (Du et al., 2020). Sugar metabolism is an important process for root development under drought stress (Xu et al., 2015;Du et al., 2020).
Drought stress has a negative impact on cotton yield and productivity. Like other abiotic stresses, it is one of the most significant environmental elements influencing cotton production and growth, as well as fiber quality, due to inadequate cellulose cell formation in the bolls (Azhar and Rehman, 2018). Drought has put cotton growers' future prospects in jeopardy; thus, finding a solution to this challenge is critical (Zahid et al., 2021).
In this research, using published genome as a reference, we performed genome-wide identification and expression analyses on the four cotton species. Phylogenetic tree, gene structure, motif analysis, chromosomal location, gene ontology, subcellular localization, protein-protein interaction, and promoter analysis along with the RNA-seq validation by RT-qPCR analysis were performed to analyze the evolution and potential of the PFK gene family under drought stress tolerance and sugar metabolism. The identified drought stress response genes, as well as the PFK gene family in general, will be used to better understand genomic architecture and functional structure, as well as to characterize the PFK gene family in cotton drought tolerance and sugar metabolism studies.

Planting Materials and Stress Treatments
The seeds of three cotton species were obtained from the Institute of Cotton Research, Chinese Academy of Agricultural Sciences (Mehari et al., 2021a). The seeds were soaked in water overnight. They were sown for about 7 days in a filter paper until good germination and then transferred to a greenhouse with 16 h of light and 8 h of darkness. When the plants grew to the three true leaves stage, the seedlings were subjected to PEG-6000. The leaf and root tissues at 0 h, 24 h, and 48 h were collected. The samples were promptly frozen in liquid nitrogen and kept at −80°C until the RNA extraction was completed (Mehari et al., 2021b).

PFK Gene Identification in Cotton
To identify the PFK genes in different cotton genomes, the PF00365 number was used to search the cotton sequence database, and G. hirsutum from NAU assembly, G. barbadense from HAU assembly, Gossypium arboreum from CRI assembly, and Gossypium raimondii from JGI assembly were filtered from CottonFGD (http://www. cottonfgd.org/search/) database (Zhu et al., 2017).

Phylogenetic Analysis of PFK Genes
The ClustalX2 software was used to align the protein sequences of all the identified PFK genes from the four cotton species, that is, G. hirsutum, G. arboreum, G. barbadense, and G. raimondii as well as A. thaliana and Theobroma cacao using default parameters. The phylogenetic analysis of PFK from all species was constructed using the numerous sequence alignments imported and displayed in MEGA 7.0 using the neighbor joining approach. For statistical reliability, tree nodes were calculated using the Bootstrap method with 1,000 repeats (Tamura et al., 2011).
Gene Architecture and Chromosomal Location of PFK Genes GSDS 2.0 (https://gsds.cbi.pku.edu.cn/) was used to illustrate the gene structure of PFK genes using Genomic DNA and CDS sequences of each species . The CottonFGD genome annotation files were used to identify the chromosomal locations of PFK genes in G. hirsutum, G. arboreum, G.
barbadense, and G. raimondii. This information was also used to construct chromosomal mapping which was then displayed by TBtools . The PFK protein sequences of the four cotton species were submitted to the online motif and domain identification software MEME (http://meme-suite.org/) to determine the conserved domains found in the proteins. The motif search was carried out using a total count of 10 motifs. The MAST tool was used to show the protein database for the identified motifs (Bailey et al., 2009).

Gene Ontology and Subcellular Localization Prediction
The Gene ontology analysis was performed by using the data downloaded from CottonFGD (www.cottonfgd.org) and the figure was visualized using the graphPad prism program (Zhu et al., 2017). The protein sequences of PFK genes were uploaded to the WoLF PSORT online website (https://wolfpsort.hgc.jp) for predictions of subcellular localization.

Cis-Regulatory and Protein-Protein Interaction Analysis of PFK Genes
The PlantCARE (http://bioinformatics.psb.ugent.be/webtools/ plantcare/html/) online tool was used to search putative cisregulatory elements in the PFK genes' promoter regions up to 2,000 bp upstream, which were retrieved from the CottonFGD database (Lescot et al., 2002). The TBtools software was used to create a circular heatmap of cis elements . The amino acid sequences of the common hub genes GhPFK11 and GhPFK17 were used as query sequences to obtain the protein-protein interaction network by using the STRING website (https://cn.string-db.org/).

Expression Analysis of PFK Genes and RT-qPCR Assay
To study the expression of PFK genes, the gene expression profiles were obtained from previously published RNA-Seq data (NCBI accession: PRJNA663204) in the leaf and root tissues of different G. hirsutum races under drought stress. The FPKM values were used to present the PFK expression levels. RNA was extracted by TianGEN kit following the protocol guidelines and was reverse transcribed into cDNA by cDNA synthesis SuperMix for qPCR (one step for gDNA removal). Finally, the SYBR Green SuperMix kit was used to perform quantitative RT-qPCR as directed in the instruction manual. Three technical and biological replications were used in the experiment. The 2 −ΔΔCt approach was used to calculate the relative expression data (Schmittgen and Livak, 2008). The primers used in the experiment were designed via NCBI (https://www.ncbi.nlm.nih.gov/) and shown in the Supplementary Table S1.

PFK Gene Family Identification in Cotton Species
The CottonFGD database (www.cottonfgd.org) was used to identify the PFK encoding genes in both diploid and tetraploid cotton species. The four cotton species were found to have a total of 80 genes, with 26 genes in G. hirsutum and G. barbadense and 14 genes in the diploid species of G. arboreum and G. raimondii. The CDS length of the cotton species ranges between 372 and 1854 bp, 1302 and 1851 bp, 372 and 1851 bp, and 840 and 1854 bp in G, hirsutum, G. barbadense, G. arboreum, and G, raimondii, respectively. Similarly, 299-617 aa, 299-617 aa, 299-617 aa, and 299-617 aa protein length was observed in the above four cotton species successively. Correspondingly, molecular weight distributes between 14.192-67.555 kDa and 47.844-67.38 kDa were recorded in G. hirsutum and G. barbadense, while 14.208-67.41 kDa and 30.493-67.568 kDa were recorded in both G. arboreum and G. raimondii successively. In case of molecular charge, similar distributions were observed, −0.5-15, −1-19, −1.5-18, and −3.5-19.5, in G. hirsutum, G. arboreum, G. raimondii, and G. barbadense, correspondingly. The GRAVY nature of the four cotton species was low and negative (Tables 1-4). This implies that the PFK proteins are hydrophilic.

Phylogenetic Analysis of PFK Proteins
The evolution analysis and patterns of cotton PFK proteins were exhibited by building a phylogenetic tree in G. hirsutum, G. barbadense, G. arboreum, G. raimondii, A. thaliana, and T. cacao. The phylogenetic tree of PFK proteins was composed of two main categories, namely, PFK and PFP group (Figure 1). To construct the tree, ClustalX was used to align 97 protein sequences from G. hirsutum, G. barbadense, G. arboreum, G. raimondii, A. thaliana, and T. cacao. From this, 18 PFK and 8 PFP proteins each from G. hirsutum and G. barbadense, 10 PFK, and 4 PFP proteins from G. arboreum and G. raimondii, 7 PFK with 2 PFP from A. thaliana and 6 PFK with 2 PFP from T. cacao were used to build the tree via MEGA 7.0.

Genetic Architecture and Motif Analysis
The genetic architecture of PFK proteins was demonstrated by studying the exon-intron structural distribution in the four cotton species. The lowest exon number was recorded in genes GaPFK01 and GaPFK05 in G. arboreum, GrPFK06 in G. raimondii, GhPFK03, GhPFK05, and GhPFK14 in G. hirsutum, and GbPFK12 in G. barbadense with 2 exons while the highest number of exons "18" was present in GaPFPA1 in G. arboreum, GrPFPA1 in G. raimondii, GhPFPA1 and GhPFPA2 in G. hirsutum, GbPFPA1 and GbPFPA2 in G. barbadense. In relative to gene length, 10,005 bp, 8536 bp, 8008 bp, and 8016 bp were observed in GhPFK18, GbPFK18, GaPFK10, and GrPFK10 in the four species, respectively ( Figure 2A). This showed us the lowest exon numbers in PFK gene group and the highest exon numbers which harbors the PFP gene group. By analyzing the protein sequences of the PFK gene family using the MEME online tool, 10 conserved motifs were predicted. The identified motifs varied between 6 and 10 in G. hirsutum, 3 and 10 in G. arboreum, 3 and 10 in G. raimondii, and 5 and 10 in G. barbadense. The lowest number of motifs were observed in GaPFK01 with motif numbers of 4, 3, and 1 and GrPFPB1 with 8, 5, and 9 motifs. The PFP subgroup has lower number of motifs between 3 and 6, with motif number 1 as a common motif ( Figure 2B). Generally, the FPK proteins have highly conserved motifs and as a result, proteins with comparable structures are likely to have similar functional tasks.

Chromosomal Mapping of PFK Genes
The distribution of PFK genes on chromosomes varies between the tetraploid and diploid cotton species. In the allotetraploid species of G. hirsutum and G. barbadense, At01, At03, At04, At08, At09, At12 from At subgenome and Dt01, Dt04, Dt08, Dt09, and Dt12 from Dt subgenome had no genes on any chromosome. The gene distribution was between 1 and 4 genes among the rest 15 chromosomes with two scaffolds on the G. hirsutum At subgenome (Figure 3). Similarly, in the diploid cotton species of G. arboreum and G. raimondii, A04, A08, A09, A12 and D02, D04, D06, D08, D12 chromosomes did not possess any PFK genes. The PFK genes were irregularly distributed among the 9 and 8 chromosomes of G. arboreum and G. raimondii, successively between 1 and 4 genes ( Figures 3A-F).

Gene Ontology Enrichment Analysis
Gene ontology enrichment studies were performed on the genes of G. hirsutum, G. barbadense, G. arboreum, and G. raimondii to better realize the functional annotations of PFK family genes in Gossypium species. Gene ontology analysis proved their glycolytic process (79 genes) and fructose-6-phosphate metabolic process (77 genes) as biological function and 6-phosphofructokinase activity (80 genes), diphosphate-fructose-6-phosphate 1-phosphotransferase activity (23 genes), and ATP binding (73 genes) in the molecular function category in all PFK genes of the Gossypium species (Figure 4).

Subcellular Localization Prediction of PFK Proteins
According to the online subcellular localization prediction in the WoLF PSORT (https://wolfpsort.hgc.jp/) tool, the four Gossypium species were expected to be localized in various cell sections (Supplementary Table S2), primarily in the chloroplast, cytoplasm, cytoskeleton, and nucleus. A small number of genes were also found in the mitochondria, plastids, and perisome cells ( Figure 5).

Cis-Regulatory Element Analysis
CottonFGD (www.cottonfgd.org) was used to identify cisregulatory components in the 2,000 base pair of the 5′ upstream region. The findings revealed that cis-acting components fall into numerous categories such as defense and stress responsiveness, abscisic acid responsiveness, light responsiveness, anaerobic induction elements, MeJAresponsiveness, hormone-responsive elements, low-temperature responsiveness, and in drought inducibility elements were identified in the promoter regions of both the tetraploid and diploid cotton species (Figure 6, Supplementary Table S3). In addition, the promoter regions of G. hirsutum, G. barbadense, G. arboreum, and G. raimondii contained five elements important in drought and stress reactivity, including TC-rich repeats, MBS, ABRE, ARE, and WUN-motif.

Protein-Protein Interaction Analysis of PFK Genes
Protein-protein interaction analysis was used to further identify the potential biological activities of PFK hub genes, and 10 putative interactors were discovered (Figure 7). Notably, several sugar metabolism-related proteins that play a key role in glycolysis and gluconeogenesis (FBA4, FBA7, and FBA8) were found. AT5G42740 and PGI1 encode glucose-6-phosphate and plastid phospho-glucose isomerase which encode for defense response and accumulation of starch in root cap cells. TIM and TPI proteins also encode plastidic triose phosphate and triosephosphate isomerase enzymes.

RT-qPCR and Gene Co-Expression Network Analysis
Transcriptome analysis was performed in G. hirsutum races to observe the expression potential of PFK genes in response to drought stress tolerance. The results of the RNA-seq revealed that in both leaf and root tissues the genes of PFK family showed higher expressions in response to drought stress. GhPFK02, GhPFK04, GhPFK05, GhPFK09, GhPFK11, GhPFK13, GhPFK14, and GhPFK17 in the leaf tissue, GhPFK02, GhPFK06, GhPFK09, GhPFK11, GhPFK15, GhPFK16, and GhPFK17 in the root tissue showed consistent upregulation in all the three cotton races at all time points (Figures 8A,C). To validate the transcriptome results, RT-qPCR analysis was performed in 26 PFK genes in similar tissues and time points as the transcriptome. The RT-qPCR analysis verified a similar trend of expression with high correlation (R 2 = 77.5% in the leaf and R 2 = 64.9% in the root) of both transcriptomes and RT-qPCR analysis results. For gene co-expression network analysis, a total of 26 genes were used from the PFK gene family. We used the RNAsequencing data from both leaf and root tissues of G. hirsutum races. A correlation analysis of these genes was performed via a correlation calculator. After evaluating the correlation of the genes, the co-expression network analysis was established by cytoscape v.3.7.2 (Shannon et al., 2003). In both leaf and root tissues, 26 nodes (genes) and 325 edges (correlation) were observed. In leaf tissue expressions, there were 191 positive and 134 negative correlations and in root tissue expressions 174 positive and 151 negative correlations were identified. Results indicated that GhPFK04, GhPFK05, GhPFK09, GhPFK11, GhPFK13, GhPFK14, and GhPFK17 were found as hub genes in leaves. Similarly, in roots, GhPFK02, GhPFK09, GhPFK11, GhPFK15, GhPFK16, and GhPFK17 were found as hub genes via the cytohubba analysis (Figures 8B,D. Supplementary  Table S4). Surprisingly we found that GhPFK11 and GhPFK17 are common hub genes in both leaves and roots.

Metabolite Enrichment Analysis
KEGG enrichment analysis was performed in the G. hirsutum races to identify the significant metabolites that are crucial for drought stress tolerance. Biosynthesis of amino acids, carbon metabolism, fructose and mannose metabolism, and galactose metabolism were all shown to be significantly enriched for the PFK genes (Supplementary Table S5). In the biosynthesis of amino acids and carbon metabolism, most of the amino acids like L-valine, L-histidine, L-glutamine, L-serine, L-homoserine, L-methionine, L-cysteine, and gluconic acid were significantly upregulated, whereas in fructose and mannose, metabolism and galactose metabolism sugars mainly fructose-1-phosphate, D-mannitol, D-sorbitol, dulcitol, and lactose were significantly downregulated under drought stress (Figure 9). In summary, amino acids were upregulated, whereas sugars were downregulated under drought stress treatment.

DISCUSSION
Cotton (genus Gossypium) is a valuable crop that produces the majority of the world's natural fiber. It is a significant crop that has seen substantial production increases during the last century (Aslam et al., 2020;Grover et al., 2020). Gossypium contains more than 50 wild species that are related to cultivated cottons. These species are divided into groups of nearly related species arranged "A-G″ and "K" for diploids and "AD" for tetraploids, and they could provide biotic and abiotic stress tolerance. (Wendel and Grover, 2015;Li et al., 2021). In this study, we used G. arboreum and G. raimondii from the diploid species category of A and D genomes and G. barbadense and G. hirsutum, belonging to the AD-genome tetraploid category (Zhang et al., 2008). G. barbadense and G. hirsutum now account for >95% of the world's cotton production (Fang et al., 2014;Li et al., 2015).
In the four Gossypium species, we identified a total of 80 PFK genes with a majority of PFK (56 genes) and 24 genes from PFP subfamily. The PFP gene subfamily was also categorized as PFPA (Alpha) with 6 genes and PFPB (Beta) with 18 genes. Similarly, in other researched crops, the gene family was divided between PFK and PFP groups. The PFK gene family has eleven members in A. thaliana: four belongs to AtPFP and seven to AtPFK. Rice has fifteen PFK genes, five of which belong to OsPFP and ten to OsPFK. The white pear (Pyrus bretschneideri) has fourteen members in the PFK family, ten of which are PbPFK and four of which are PbPFP. PFK is a crucial protein in plant growth and development with a wide range of functions, with 13 members of the PFK family in cassava (Manihot esculenta Crantz) . According to a subcellular localization analysis of PFK genes, five MePFKs in chloroplast, two in cytoplasm, and four MePFPs were predicted to be localized at the cytoplasm in cassava (M. esculenta Crantz). MePFK03 and MePFPA1, which are expected to be found in the chloroplast and cytoplasm, successively, were chosen to generate GFP fusion proteins to support the abovementioned conclusion. GFP images revealed that they were present in the chloroplast and cytoplasm, successively .
The discovery of cis-regulatory elements unveiled that the PFK genes contain several essential cis-regulatory elements that protect cells from harmful environmental stimuli. Box 4, CGTCA-motif, DRE1, DRE core, ERE, GARE-motif, G-Box, GATA-motif, GT1-motif, LTR, MBS, MYB-like sequence, MYB recognition site, MYB, MYC, STRE, W-box, WRE3, and WUN-motif were some of the vital regulatory elements identified in the PFK family. OMT genes contain W-box, MYB, MYC, DRE, ABRE, G-Box, and MBS. The W-box regulates gene expression and binds WRKY transcription factors. Plants need WRKY TFs to protect themselves from drought, chilling, wounding, heat, and salinity stress (Hafeez et al., 2021).
Cis-regulatory elements belong to numerous classifications such as stress responsive, hormone responsive, and light responsive. They were encoded in the promoter regions of G. hirsutum, G. arboreum, G. raimondii, and G. barbadense. Furthermore, four abiotic stress-responsive promoter elements were found in two diploid and two tetraploid cotton species, including LTRs, TC-rich repeats, MBSs, and WUN motifs confirming their role in abiotic stress tolerance (Rehman et al., 2021). ABRE, CGTCA-motif, TCA-element, TGACG-motif, GARE-motif, TATC-box, P-box, and other phytohormone linked abscisic acid, gibberellin, auxin, methyl jasmonate (MeJA), and salicylic acid were identified in the SOD family of Brassica napus in a similar study. Stress response (drought, cold, FIGURE 8 | RNA-seq, RT-qPCR, and Co-expression network analysis of PFK genes. (A). RNA-seq and RT-qPCR expression analysis of the PFK genes in the leaf tissue of G. hirsutum races; (B). Co-expression analysis and hub gene identification from leaf expression; (C). RNA-seq and RT-qPCR expression analysis of the PFK genes in the root tissue of G. hirsutum races; (D). Co-expression analysis and hub gene identification from leaf expression. Genes in the center of the network analysis with red color are hub genes.
Frontiers in Genetics | www.frontiersin.org June 2022 | Volume 13 | Article 922024 light, and anaerobic) elements such as ARE, chs-CMA1a, LAMPelement, LTR, MBS, TCT-motif, G-box, GT1-motif, and MBS were also found (Su et al., 2021). Drought is one such negative environmental cue that hinders photosynthetic carbon fixation and impacts sugar transport by reducing the cellular osmotic potential. Sugar transporters are a group of proteins that aid in the transit of sugar within cells. The influx/outflow of different sugars and their metabolite intermediates that enable plant growth and development is determined by these transporter proteins. Sugar distribution across the cellular and subcellular compartments is reprogrammed as a result of abiotic stress, particularly drought stress-induced damage (Kaur et al., 2021).
Glycolysis is a metabolic process that converts glucose into pyruvate in the availability of oxygen or lactate molecules without oxygen via enzyme-catalyzed processes. Anaerobic glycolysis, the latter pathway, is thought to be the first natural method to create adenosine triphosphate (ATP) (Mayes and Bender, 2003). It is a crucial metabolic mechanism in all living things. Glycolysis helps cells adjust to abiotic conditions including lack of oxygen, cold, and drought by providing energy for cell and sugar metabolism .
Sugars are plants' basic products, and the hexose sugars glucose and fructose provide the starting material for almost all metabolic pathways in crops. Hexose sugars are necessary to be phosphorylated before they can be digested. Hexokinases and fructokinases are the only two families of enzymes found in plants that catalyze the necessary permanent phosphorylation of glucose and fructose. These enzymes appear to coordinate sugar synthesis with water, sunlight, nutrients, and CO 2 absorption (Granot et al., 2014). In plants, fructose was formed as a breakdown product of both primary sucrose-cleaving enzymes. A fructokinase phosphorylates fructose before it enters metabolism. The phosphofructokinase B family, which includes known FRKs, is a heterogeneous group of carbohydrate/purine kinases (Riggs et al., 2017). The primary step of fructose phosphorylation in plants is catalyzed by two types of fructokinases: cytosolic and plastidic fructokinases (Stein et al., 2016). Fructokinases have an impact on long-lasting developmental processes such as sugar, nutrients, and water transport (Granot et al., 2014).
Metabolite enrichment analysis was used to directly investigate a group of functionally related metabolites that are significantly enriched. Biosynthesis of amino acids, carbon metabolism, fructose and mannose metabolism, and galactose metabolism were significantly enriched during drought stress. Many amino acids response to abiotic stressors have been discovered to play vital functions in plant growth. Amino acids also serve as precursors for a variety of primary and secondary metabolites, and they play an important role in human nutrition (Trovato et al., 2021). In both normally watered and drought-stressed plants, exogenous addition of amino acids can improve nutritional traits like phenol, total protein, proline, and several important amino acids like glutamic acid, glutamine, and asparagine. By reducing the negative effects of drought stress on cabbage, amino acids were found to be efficient in boosting physiological and nutraceutical parameters. Thus, the use of amino acids proved successful in reducing the impacts of drought stress and boosting nutritional value in drought stress situations (Haghighi et al., 2020). Overall, The above results widen our understanding about the role of Gossypium PFK genes on drought stress tolerance and sugar metabolism.

CONCLUSIONS
The glycolysis and pentose phosphate pathways have a significant impact on crop tolerance under abiotic stress. Genome-wide identification, co-expression, RT-qPCR profiling, and metabolite enrichment analysis indicated that PFK gene family has potential in drought stress tolerance by involving in biosynthesis of amino acids and sugar metabolism activity ( Figure 10). Important amino acids like L-valine, L-glutamine, L-methionine, and gluconic acid were significantly upregulated, while sugars like D-mannitol, D-sorbitol, and dulcitol were significantly downregulated during drought stress. The coexpression and RT-qPCR results showed that GhPFK04, GhPFK09, GhPFK13, GhPFK14, GhPFK16, and GhPFK17 are hub genes and consistently upregulated in both leaf and root tissues during the stress condition. GhPFK11 and GhPFK17 were identified as common hub genes and these might be the true candidate genes involved in the drought stress tolerance in G. hirsutum. The current study lays the groundwork for the importance of PFK gene family in drought stress tolerance and sugar metabolism in cotton and opens doors to functional characterization of hub genes. Therefore, further functional validation of the hub genes is needed in order to better understand their involvement in the drought stress tolerance and sugar metabolism mechanisms at molecular and genetic levels.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
TGM: conceptualization, performing the experiment, and writing the original draft. YX: conceptualization and methodology. MJU: writing review and editing. FH: data analysis. YH, ZZ, and XC: resources, software, and validation. BW, KW, and FL: providing funding for the project and supervision. All authors contributed to the article and approved the submitted version. FUNDING FIGURE 10 | A graphic model for the role of PFK genes under drought stress. Drought stress induces biosynthesis and upregulation of amino acids production and activates the plant osmo regulation activity. Additionally, PFK involves in sugar metabolism and acts as the energy source for root development during plant drought stress tolerance.