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

Drought has been identified 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 identification led to the identification 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 identified hub genes. Interestingly, GhPFK11 and GhPFK17 were identified 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 significantly upregulated, whereas sugars, mainly fructose-1-phosphate, D-mannitol, D-sorbitol, dulcitol, and lactose, were significantly 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.

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 (Zhang et al., 2015). 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: ATP-dependent 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 (Wang et al., 2021), 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.

Materials and Methods

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 (Hu et al., 2015). 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 (Chen et al., 2020). 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 cis-regulatory 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 (Chen et al., 2020). 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.

Results

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.

TABLE 1

TABLE 1. Transcript and physiochemical features of PFK genes in G. hirsutum.

TABLE 2

TABLE 2. Transcript and physiochemical features of PFK genes in G. arboreum.

TABLE 3

TABLE 3. Transcript and physiochemical features of PFK genes in G. raimondii.

TABLE 4

TABLE 4. Transcript and physiochemical features of PFK genes in G. barbadense.

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.

FIGURE 1

FIGURE 1. Phylogenetic tree analysis of PFK genes in G. hirsutum, G. barbadense, G. arboreum, G. raimondii, A. thaliana, and T. cacao. The tree analysis categorized the PFK gene family into two big groups of PFK and PFP which are highlighted in blue and violet colors.

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.

FIGURE 2

FIGURE 2. Gene architecture and motif enrichment analysis of PFK proteins in Gossypium species. I. G. hirsutum; II. G. barbadense; III. G. arboreum; IV. G. raimondii. (A–C) stands for phylogenetic relationship, gene structure, and motif identification in all the Gossypium species, respectively.

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).

FIGURE 3

FIGURE 3. Chromosomal location of the PFK genes. (A) G. hirsutum At subgenome; (B) G. hirsutum Dt subgenome; (C) G. barbadense At subgenome; (D) G. barbadense Dt subgenome; (E) G. arboreum; (F) G. raimondii.

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).

FIGURE 4

FIGURE 4. Gene ontology enrichment analysis using CottonFGD for the PFK gene family. (A) G. hirsutum; (B) G. barbadense; (C) G. arboreum; (D) G. raimondii. Bar graph was done using PrismPad; the scale indicates the number of genes enriched in each function.

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).

FIGURE 5

FIGURE 5. Subcellular localization prediction of the PFK genes of Gossypium species (A) G. hirsutum; (B) G. arboreum; (C) G. raimondii; (D) G. barbadense.

Cis-Regulatory Element Analysis

CottonFGD (www.cottonfgd.org) was used to identify cis-regulatory 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, MeJA-responsiveness, 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.

FIGURE 6

FIGURE 6. Cis-regulatory elements analysis of the PFK gene family. (A) G. hirsutum; (B) G. barbadense; (C) G. arboreum; (D) G. raimondii.

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.

FIGURE 7

FIGURE 7. Protein–protein interaction analysis of PFK hub proteins. (A). Interaction network of GhPFK11 gene; (B). Interaction network of GhPFK17 gene; Colored nodes: query proteins and first shell of interactors; Red lines: gene fusions; Green lines: gene neighborhood; Black lines: co-expression; Blue lines: gene co-occurrence.

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² = 77.5% in the leaf and R² = 64.9% in the root) of both transcriptomes and RT-qPCR analysis results.

FIGURE 8

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.

For gene co-expression network analysis, a total of 26 genes were used from the PFK gene family. We used the RNA-sequencing 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.

FIGURE 9

FIGURE 9. Metabolite enrichment analyses of PFK genes during drought stress. (A). Biosynthesis of amino acids; (B). Carbon metabolism; (C). Galactose metabolism; (D). Fructose and mannose metabolism. Heatmap was visualized by TBtools using the log2 data normalization of the metabolite expression.

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) (Wang et al., 2021).

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 (Wang et al., 2021).

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, light, and anaerobic) elements such as ARE, chs-CMA1a, LAMP-element, 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 (Wang et al., 2021).

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₂ 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 co-expression 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.

FIGURE 10

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.

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

This research was funded by the National Key R&D Program of China [(2021YFE0101200), PSF/CRP/18thProtocol (07)], the National Natural Science Foundation of China (32171994, 32072023).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2022.922024/full#supplementary-material

References

Aslam, S., Khan, S. H., Ahmed, A., and Dandekar, A. M. (2020). The Tale of Cotton Plant: From Wild Type to Domestication, Leading to its Improvement by Genetic Transformation. Amer. J. Mol. Biol. 10 (2), 91–127. doi:10.4236/ajmb.2020.102008

CrossRef Full Text | Google Scholar

Azhar, M. T., and Rehman, A. (2018). “Overview on Effects of Water Stress on Cotton Plants and Productivity,” in Biochemical, Physiological and Molecular Avenues for Combating Abiotic Stress Tolerance in Plants (Academic Press), 297–316. doi:10.1016/B978-0-12-813066-7.00016-4

CrossRef Full Text | Google Scholar

Bailey, T. L., Boden, M., Buske, F. A., Frith, M., Grant, C. E., Clementi, L., et al. (2009). MEME SUITE: Tools for Motif Discovery and Searching. Nucleic Acids Res. 37 (Suppl. l_2), W202–W208. doi:10.1093/nar/gkp335

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, C., Chen, H., Zhang, Y., Thomas, H. R., Frank, M. H., He, Y., et al. (2020). TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data. Mol. Plant 13 (8), 1194–1202. doi:10.1016/j.molp.2020.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y., Zhang, Q., Hu, W., Zhang, X., Wang, L., Hua, X., et al. (2017). Evolution and Expression of the Fructokinase Gene Family in Saccharum. BMC Genomics 18 (1), 1–15. doi:10.1186/s12864-017-3535-7

PubMed Abstract | CrossRef Full Text | Google Scholar

da Silva, M. D., de Oliveira Silva, R. L., Ferreira Neto, J. R. C., Benko-Iseppon, A. M., and Kido, E. A. (2017). Genotype-Dependent Regulation of Drought-Responsive Genes in Tolerant and Sensitive Sugarcane Cultivars. Gene 633, 17–27. doi:10.1016/j.gene.2017.08.022

PubMed Abstract | CrossRef Full Text | Google Scholar

dos Santos, T. B., Ribas, A. F., de Souza, S. G. H., Budzinski, I. G. F., and Domingues, D. S. (2022). Physiological Responses to Drought, Salinity, and Heat Stress in Plants: A Review. Stresses 2 (1), 113–135. doi:10.3390/stresses2010009

CrossRef Full Text | Google Scholar

Du, Y., Zhao, Q., Chen, L., Yao, X., Zhang, W., Zhang, B., et al. (2020). Effect of Drought Stress on Sugar Metabolism in Leaves and Roots of Soybean Seedlings. Plant Physiol. Biochem. 146, 1–12. doi:10.1016/j.plaphy.2019.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Fang, D. D., Jenkins, J. N., Deng, D. D., McCarty, J. C., Li, P., and Wu, J. (2014). Quantitative Trait Loci Analysis of Fiber Quality Traits Using a Random-Mated Recombinant Inbred Population in Upland Cotton (Gossypium Hirsutum L.) BMC Genomics 15 (397), 1–15. doi:10.1186/1471-2164-15-397

PubMed Abstract | CrossRef Full Text | Google Scholar

Givan, C. V. (1999). Evolving Concepts in Plant Glycolysis: Two Centuries of Progress. Biol. Rev. 74 (3), 277–309. doi:10.1111/j.1469-185X.1999.tb00188.x

CrossRef Full Text | Google Scholar

Granot, D., Kelly, G., Stein, O., and David-Schwartz, R. (2014). Substantial Roles of Hexokinase and Fructokinase in the Effects of Sugars on Plant Physiology and Development. J. Exp. Bot. 65 (3), 809–819. doi:10.1093/jxb/ert400

PubMed Abstract | CrossRef Full Text | Google Scholar

Grover, C. E., Pan, M., Yuan, D., Arick, M. A., Hu, G., Brase, L., et al. (2020). The Gossypium Longicalyx Genome as a Resource for Cotton Breeding and Evolution. G3 Genes Genom. Genet. 10 (5), 1457–1467. doi:10.1534/g3.120.401050

PubMed Abstract | CrossRef Full Text | Google Scholar

Hafeez, A., Gě, Q., Zhāng, Q., Lǐ, J., Gōng, J., Liú, R., et al. (2021). Multi-Responses of O-Methyltransferase Genes to Salt Stress and Fiber Development of Gossypium Species. BMC Plant Biol. 21 (1), 1–17. doi:10.1186/s12870-020-02786-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Haghighi, M., Saadat, S., and Abbey, L. (2020). Effect of Exogenous Amino Acids Application on Growth and Nutritional Value of Cabbage under Drought Stress. Sci. Hortic. 272, 109561. doi:10.1016/j.scienta.2020.109561

CrossRef Full Text | Google Scholar

Hu, B., Jin, J., Guo, A.-Y., Zhang, H., Luo, J., and Gao, G. (2015). GSDS 2.0: An Upgraded Gene Feature Visualization Server. Bioinformatics 31 (8), 1296–1297. doi:10.1093/bioinformatics/btu817

PubMed Abstract | CrossRef Full Text | Google Scholar

Kato-Noguchi, H. (2002). The Catalytic Direction of Pyrophosphate: Fructose 6‐Phosphate 1‐Phosphotransferase in Rice Coleoptiles in Anoxia. Physiol. Plant. 116 (3), 345–350. doi:10.1046/J.0031-9317.2002.00002.X

CrossRef Full Text | Google Scholar

Kaur, H., Manna, M., Thakur, T., Gautam, V., and Salvi, P. (2021). Imperative Role of Sugar Signaling and Transport during Drought Stress Responses in Plants. Physiol. Plant. 171 (4), 833–848. doi:10.1111/ppl.13364

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamaoui, M., Jemo, M., Datla, R., and Bekkaoui, F. (2018). Heat and Drought Stresses in Crops and Approaches for Their Mitigation. Front. Chem. 6, 26. doi:10.3389/fchem.2018.00026

PubMed Abstract | CrossRef Full Text | Google Scholar

Lescot, M., Déhais, P., Thijs, G., Marchal, K., Moreau, Y., Van de Peer, Y., et al. (2002). PlantCARE, a Database of Plant Cis-Acting Regulatory Elements and a Portal to Tools for In Silico Analysis of Promoter Sequences. Nucleic Acids Res. 30 (1), 325–327. doi:10.1093/nar/30.1.325

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Fan, G., Lu, C., Xiao, G., Zou, C., Kohel, R. J., et al. (2015). Genome Sequence of Cultivated Upland Cotton (Gossypium Hirsutum TM-1) Provides Insights into Genome Evolution. Nat. Biotechnol. 33 (5), 524–530. doi:10.1038/nbt.3208

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, Z., Liu, Z., Wei, Y., Liu, Y., Xing, L., Liu, M., et al. (2021). Genome-Wide Identification of the MIOX Gene Family and Their Expression Profile in Cotton Development and Response to Abiotic Stress. PLoS One 16 (7), e0254111. doi:10.1371/journal.pone.0254111

PubMed Abstract | CrossRef Full Text | Google Scholar

Lü, H., Li, J., Huang, Y., Zhang, M., Zhang, S., and Wu, J. (2019). Genome-Wide Identification, Expression and Functional Analysis of the Phosphofructokinase Gene Family in Chinese White Pear (Pyrus Bretschneideri). Gene 702, 133–142. doi:10.1016/j.gene.2019.03.005

PubMed Abstract | CrossRef Full Text | Google Scholar

Mayes, P. A., and Bender, D. A. (2003). “Glycolysis and the Oxidation of Pyruvate,” in A LANGE Medical Book, 136. McGraw-Hill, Medical Publishing Division. doi:10.13140/RG.2.2.10309.78567

CrossRef Full Text | Google Scholar

Mehari, T. G., Xu, Y., Magwanga, R. O., Umer, M. J., Shiraku, M. L., Hou, Y., et al. (2021). Identification and Functional Characterization of Gh_D01G0514 (GhNAC072) Transcription Factor in Response to Drought Stress Tolerance in Cotton. Plant Physiol. Biochem. 166, 361–375. doi:10.1016/j.plaphy.2021.05.050

PubMed Abstract | CrossRef Full Text | Google Scholar

Mehari, T. G., Xu, Y., Umer, M. J., Shiraku, M. L., Hou, Y., Wang, Y., et al. (2021). Multi-omics Based Identification and Functional Characterization of Gh_A06G1257 Proves its Potential Role in Drought Stress Tolerance in Gossypium Hirsutum. Front. Plant Sci. 12, 746771. doi:10.3389/fpls.2021.746771

PubMed Abstract | CrossRef Full Text | Google Scholar

Mustroph, A., Sonnewald, U., and Biemelt, S. (2007). Characterisation of the ATP-Dependent Phosphofructokinase Gene Family from Arabidopsis Thaliana. FEBS Lett. 581 (13), 2401–2410. doi:10.1016/j.febslet.2007.04.060

PubMed Abstract | CrossRef Full Text | Google Scholar

Nazir, M. F., Jia, Y., Ahmed, H., He, S., Iqbal, M. S., Sarfraz, Z., et al. (2020). Genomic Insight into Differentiation and Selection Sweeps in the Improvement of Upland Cotton. Plants 9 (6), 711. doi:10.3390/plants9060711

PubMed Abstract | CrossRef Full Text | Google Scholar

Ohlendieck, K. (2010). Proteomics of Skeletal Muscle Glycolysis. Biochim. Biophys. Acta 1804 (11), 2089–2101. doi:10.1016/j.bbapap.2010.08.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Plaxton, W. C. (1996). The Organization and Regulation of Plant Glycolysis. Annu. Rev. Plant. Physiol. Plant. Mol. Biol. 47 (1), 185–214. doi:10.1146/annurev.arplant.47.1.185

PubMed Abstract | CrossRef Full Text | Google Scholar

Rehman, A., Atif, R. M., Azhar, M. T., Peng, Z., Li, H., Qin, G., et al. (2021). Genome Wide Identification, Classification and Functional Characterization of Heat Shock Transcription Factors in Cultivated and Ancestral Cottons (Gossypium spp.). Int. J. Biol. Macromol. 182, 1507–1527. doi:10.1016/j.ijbiomac.2021.05.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Riggs, J. W., Cavales, P. C., Chapiro, S. M., and Callis, J. (2017). Identification and Biochemical Characterization of the Fructokinase Gene Family in Arabidopsis Thaliana. BMC plant Biol. 17 (1), 1–18. doi:10.1186/s12870-017-1031-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Ros, S., and Schulze, A. (2013). Balancing Glycolytic Flux: The Role of 6-Phosphofructo-2-Kinase/Fructose 2, 6-Bisphosphatases in Cancer Metabolism. Cancer Metab. 1 (1), 1–10. doi:10.1186/2049-3002-1-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Schmittgen, T. D., and Livak, K. J. (2008). Analyzing Real-Time PCR Data by the Comparative CT Method. Nat. Protoc. 3 (6), 1101–1108. doi:10.1038/nprot.2008.73

PubMed Abstract | CrossRef Full Text | Google Scholar

Stein, O., Damari‐Weissler, H., Secchi, F., Rachamilevitch, S., German, M. A., Yeselson, Y., et al. (2016). The Tomato Plastidic Fructokinase Sl FRK 3 Plays a Role in Xylem Development. New Phytol. 209 (4), 1484–1495. doi:10.1111/nph.13705

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, J., Li, L., Pang, C., Wei, H., Wang, C., Song, M., et al. (2016). Two Genomic Regions Associated with Fiber Quality Traits in Chinese Upland Cotton under Apparent Breeding Selection. Sci. Rep. 6 (1), 38496. doi:10.1038/srep38496

PubMed Abstract | CrossRef Full Text | Google Scholar

Su, W., Raza, A., Gao, A., Jia, Z., Zhang, Y., Hussain, M. A., et al. (2021). Genome-Wide Analysis and Expression Profile of Superoxide Dismutase (SOD) Gene Family in Rapeseed (Brassica Napus L.) under Different Hormones and Abiotic Stress Conditions. Antioxidants 10 (8), 1182. doi:10.3390/antiox10081182

PubMed Abstract | CrossRef Full Text | Google Scholar

Tamura, K., Peterson, D., Peterson, N., Stecher, G., Nei, M., and Kumar, S. (2011). MEGA5: Molecular Evolutionary Genetics Analysis Using Maximum Likelihood, Evolutionary Distance, and Maximum Parsimony Methods. Mol. Biol. Evol. 28 (10), 2731–2739. doi:10.1093/molbev/msr121

PubMed Abstract | CrossRef Full Text | Google Scholar

Trovato, M., Funck, D., Forlani, G., Okumoto, S., and Amir, R. (2021). Editorial: Amino Acids in Plants: Regulation and Functions in Development and Stress Defense. Front. Plant Sci. 12, 772810. doi:10.3389/fpls.2021.772810

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, H., Zhao, P., Shen, X., Xia, Z., Zhou, X., Chen, X., et al. (2021). Genome-Wide Survey of the Phosphofructokinase Family in Cassava and Functional Characterization in Response to Oxygen-Deficient Stress. BMC Plant Biol. 21 (1), 1–15. doi:10.1186/s12870-021-03139-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, M., Tu, L., Yuan, D., Zhu, D., Shen, C., Li, J., et al. (2019). Reference Genome Sequences of Two Cultivated Allotetraploid Cottons, Gossypium Hirsutum and Gossypium Barbadense. Nat. Genet. 51 (2), 224–229. doi:10.1038/s41588-018-0282-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Wendel, J. F., and Grover, C. E. (2015). Taxonomy and Evolution of the Cotton Genus, Gossypium. Cotton 57, 25–44. doi:10.2134/agronmonogr57.2013.0020

CrossRef Full Text | Google Scholar

Xu, W., Cui, K., Xu, A., Nie, L., Huang, J., and Peng, S. (2015). Drought Stress Condition Increases Root to Shoot Ratio via Alteration of Carbohydrate Partitioning and Enzymatic Activity in Rice Seedlings. Acta Physiol. Plant 37 (9), 1–11. doi:10.1007/s11738-014-1760-0

CrossRef Full Text | Google Scholar

Yao, K., and Wu, Y. Y. (2016). Phosphofructokinase and Glucose-6-Phosphate Dehydrogenase in Response to Drought and Bicarbonate Stress at Transcriptional and Functional Levels in Mulberry. Russ. J. Plant Physiol. 63 (2), 235–242. doi:10.1134/S102144371602014X

CrossRef Full Text | Google Scholar

Zahid, Z., Khan, M. K. R., Hameed, A., Akhtar, M., Ditta, A., Hassan, H. M., et al. (2021). Dissection of Drought Tolerance in Upland Cotton through Morpho-Physiological and Biochemical Traits at Seedling Stage. Front. Plant Sci. 12, 260. doi:10.3389/fpls.2021.627107

CrossRef Full Text | Google Scholar

Zhang, H.-B., Li, Y., Wang, B., and Chee, P. W. (2008). Recent Advances in Cotton Genomics. Int. J. Plant Genomics 2008, 742304. doi:10.1155/2008/742304

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, T., Hu, Y., Jiang, W., Fang, L., Guan, X., Chen, J., et al. (2015). Sequencing of Allotetraploid Cotton (Gossypium Hirsutum L. Acc. TM-1) Provides a Resource for Fiber Improvement. Nat. Biotechnol. 33 (5), 531–537. doi:10.1038/nbt.3207

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhu, T., Liang, C., Meng, Z., Sun, G., Meng, Z., Guo, S., et al. (2017). CottonFGD: An Integrated Functional Genomics Database for Cotton. BMC Plant Biol. 17 (1), 1–9. doi:10.1186/s12870-017-1039-x

PubMed Abstract | CrossRef Full Text | Google Scholar