Glyceraldehyde-3-phosphate dehydrogenase Gh_GAPDH9 is associated with drought resistance in Gossypium hirsutum

Plant material

Based on previous research, drought-resistant (Acala1517-08 (which originated as a single-plant selection derived from a cross between B7636, an unreleased breeding line, and LA 887, which is moderately resistant to root-knot nematodes) and Tashigan7 hao (which was bred in the former Soviet Union and preserved by the Xinjiang Bazhou Institute of Agricultural Sciences, numbered ZM-06734)) and drought-sensitive (Xinluzhong36 (which was selected and bred by the Bayingolin Autonomous Prefecture Agricultural Science Research and Xinjiang registration number is 20100108) and Xinluzao26 (which is a new variety cultivated by Xinjiang Tianhe Seed Industry Co., Ltd. in 2006, derived from the variant of Xinlu Zao No. 8)) upland cotton materials were selected, all of which were conserved and provided by the Germplasm Innovation Laboratory of Xinjiang Agricultural University (Fan et al., 2021). The cotton seeds of the four materials were sterilized by soaking in hydrogen peroxide for 4 h, washed in sterile water and soaked for 24 h to germinate the seeds. Then, germination boxes that were 18 cm long, 12 cm wide and 6 cm high were selected, and 50 seeds were placed in each box. The boxes were then placed in a constant-temperature (25 °C) incubator with a light intensity of 1200 lx or more and a relative humidity of 60–80%, set to 8 h of darkness and 16 h of light, and incubation was continued until the cotyledons expanded. The cotton plants were then transferred to white plastic pots that were 60 cm long, 45 cm wide and 20 cm high and cultivated in 50% Hoagland’s nutrient solution. The nutrient solution was changed once every 2 d (to ensure that the nutrient solution did not pass by the roots), three true leaves were allowed to grow, and processing was started. Before treatment, the water in each hydroponic pot was replaced with clean, fresh water, in which incubation was continued for 2 d. Thereafter, the water in each hydroponic pot was replaced with a 15% solution of PEG6000 configured to simulate drought treatment, and root, stem and leaf tissues were collected after 0, 0.5, 1, 3, 6, 12 and 24 h of treatment and stored frozen in a −80 °C refrigerator. Three biological replicates of each treatment were performed. Columbia wild-type Arabidopsis thaliana (WT) was used as receptor material for transfecting the GAPDH9 gene.

Identification and bioinformatics analysis of the GAPDH gene family in cotton

The genomic and proteomic data for G. arboreum (CRI, version 1.0), G. raimondii (JGI, version 2.0), G. hirsutum (HUA, version 1.0) and G. barbadense (ZJU, version 1.1) were obtained from the COTTONGEN (http://www.cottongen.org/) database (Yu et al., 2014). The Cryptomarkov model of the structural domain of the GAPDH genes (PF00044 and PF02800) and HMMER software (http://hmmer.org/) were used for analysis (Potter et al., 2018). Using known GAPDH homologous sequences, machine learning in the algorithm generates a hidden Markov model, which we then use to predict homologous gene sequences in cotton. Such a hidden Markov model is actually a file on the Pfam database with the suffix “hmm” (PF00044 and PF02800) (Potter et al., 2018). After the removal of redundancy, all candidate genes were subjected to structural domain validation in the NCBI-CDD (https://www.ncbi.nlm.nih.gov/cdd/) database (Lu et al., 2020). The number of amino acid residues, relative molecular masses, and theoretical isoelectric points of the upland cotton GAPDH proteins were calculated using ExPASy (http://cn.ExPASy.org/) (Mariethoz et al., 2018) and Plant-mPLoc online software (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) to predict the subcellular localization of upland cotton GAPDH proteins (Chou & Shen, 2010).

Phylogenetic and collinear analysis of the GAPDH gene family in cotton

To further understand the evolution and differentiation of the GAPDH genes in terrestrial cotton, the GAPDH protein sequences of different cotton species were compared using Clustal W in MEGA 7 software for multiple sequence alignment, and then a phylogenetic tree was constructed using the proximity method based on the results of the alignment, where the bootstrap value was set to 1,000 (Kumar, Stecher & Tamura, 2016). The resulting phylogenetic trees were embellished using the online tool Evolview (https://evolgenius.info/). Local BLAST searches were performed using GAPDH protein sequences from the four cotton species with the e value set to 1e-10. The results of the comparisons were filtered and analyzed using MCScanX software (Wang et al., 2012b), and the results were plotted using CIRCOS software (Zhang, Meltzer & Davis, 2013).

Chromosomal localization, gene structure and motif analysis of the GAPDH gene family in upland cotton

To further understand the composition of the terrestrial cotton GAPDH genes and their locations on the chromosomes, information on the chromosomal locations of upland cotton GAPDH gene family members was extracted using TBtools software version 1.120 (Chen et al., 2020), and the chromosomal locations of upland cotton GAPDH genes were mapped using MapChart software version 2.32 (Voorrips, 2002). An evolutionary tree of upland cotton GAPDH was constructed by MEGA 7 software to obtain NWK files (Kumar, Stecher & Tamura, 2016). Motif analysis was performed with the MEME program (setting the number of functional domains to 10, http://meme-suite.org/) (Bailey et al., 2009). All of the above results were visualized using TBtools software version 1.120 (Chen et al., 2020).

Analysis of the cis-acting elements in upland cotton GAPDH gene promoters

Transcription factors control transcription by identifying specific DNA sequences in gene promoter regions to guide genome expression and thus regulate plant function. It is generally believed that the 2,000 bp region upstream of a gene is the promoter region of that gene. The 2,000 bp DNA sequences upstream of the GAPDH genes were extracted as promoter sequences based on the positional information of the upland cotton GAPDH genes using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) to predict possible cis-acting elements, which were visualized using TBtools software (Chen et al., 2020).

RNA-seq analysis

Transcriptomic data from different tissues (calycle, leaf, petal, pistil, root, stamen, stem and torus) of upland cotton and data obtained after drought, cold, heat and salt stress treatments were downloaded from the NCBI SAR (Sequence Read Archive) database (genome sequencing project accession: PRJNA532694). The data were filtered and subjected to quality control using fastp software, and the resulting clean data were used for subsequent analysis. The upland cotton (TM-1) genome (https://www.cottongen.org/species/Gossypium_hirsutum/HAU-AD1_genome_v1.0_v1.1) was used as the reference genome for alignment using HISAT2 software, and String Tie was used to compare the postalignment sequences for expression quantification (Hu et al., 2019). Different RNAs may have different lengths, and the longer the length, the more reads there are. When each RNA is divided by its length, the relative expression of different genes in the same sample can be compared. It’s important to note that various samples may have been sequenced at different sequencing depths, and the deeper the sequencing depth, the greater the number of corresponding reads. Normalizing the results by the number of reads per million in their respective libraries provides a robust measure of the relative expression of the same gene in two distinct samples. The FPKM (fragments per kilobase exon per million fragments mapped) value is the number of reads per kilobase mapped to one exon per million reads in the map. The FPKM method was used to assess gene expression. Heatmaps were created using TBtools software (Chen et al., 2020).

qRT–PCR

The qRT‒PCR procedure was performed according to Li et al. (2022). Total RNA was isolated by referring to the Polysaccharide Polyphenol Total RNA Extraction Kit instructions (DP411; Tiangen, Beijing, China). cDNA was synthesized using a FastKing RT Kit (KR116; Tiangen, Beijing, China). Specific primers were designed using the NCBI Primer Tool, and the primer sequences are shown in Table S1. qRT‒PCR was performed using the ABI 7500 platform and PerfectStart Green qPCR SuperMix (AQ601; TransGen Biotech, Beijing, China). The reaction program was a thermal cycling program at 94 °C for 30 s followed by 40 cycles of 95 °C for 5 s, 57 °C for 5 s, and 72 °C for 34 s (three biological replicates). The internal reference gene was GhUBQ7, and the results were calculated using the 2^−ΔΔCt method (Zhao et al., 2021).

VIGS and overexpression

VIGS technology produces RNAi directly in contemporary plants, resulting in significant time and cost savings. It can avoid the drawbacks of traditional transgenic systems and can downregulate the expression of deletion-lethal genes. A 409 bp fragment of GhGAPDH9 was amplified from the cDNA of Acala1517-08 at the EcoR and XbaI sites. pTRV2-GhGAPDH9 was successfully constructed and transformed into Agrobacterium tumefaciens GV3101. The TRV vector-carrying Agrobacterium cells were mixed and injected into the cotyledons of 7-d-old seedlings of the drought-tolerant material Acala1517-08, which was then incubated at 22–25 °C under a 16 h/8 h light/dark cycle after 24 h of protection from light. The silencing efficiency was verified by qRT‒PCR in pTRV2-GhGAPDH9 and control pTRV2-injected plants after the appearance of bleaching in leaves injected with pTRV2-CLA. After the confirmation of silencing, watering was stopped until phenotypic differences appeared between silenced and control cotton. Under drought conditions, the changes in malondialdehyde (MDA), proline (Pro) and chlorophyll (Chl) contents were important indicators to judge the drought resistance of plants. The physiological indicators of pTRV2-GhGAPDH9 and control pTRV2, including MDA, Pro and Chl, were subsequently determined according to Fan et al. (2021). Overexpression can cause the gene expression product to exceed normal levels, enabling observation of the biological behavior of the gene and thus elucidating the function of the gene. A 1,224 bp fragment of GhGAPDH9 was amplified from the cDNA of Acala1517-08 at the SacI sites, and pCAMBIA3301-GhGAPDH9 was successfully constructed and transformed into Agrobacterium tumefaciens GV3101 (Yu et al., 2020). Wild-type Arabidopsis and GhGAPDH9-overexpressing Arabidopsis plants were grown and transplanted, and overexpression of the GhGAPDH9 gene in pure strains was determined using qRT‒PCR. The plants were subjected to natural drought treatment (from the time of treatment, watering was stopped, simulating natural drought conditions), and the phenotypes were observed periodically after the drought treatment and photographed and recorded according to the degree of leaf wilt in Arabidopsis. All primers used in this experiment are listed in Table S1.

Statistical analysis

Excel 2019 was used to collate physiological indicators and qRT‒PCR data, SPSS 26.0 software (SPSS Inc., Chicago, IL, USA) was used to analyze the data by ANOVA, and R language was used for graphing (R Core Team, 2020).

Identification of the GAPDH gene family in cotton

Based on the protein functional domains PF00044 and PF02800, a total of 84 GAPDH proteins encoded by four cotton genomes (G. arboreum, G. raimondii, G. hirsutum and G. barbadense) were identified, and the results of the search were validated in the NCBI-CDD database. The copy number changes in the GAPDH gene family during the evolution of cotton were systematically investigated. The G. arboreum, G. raimondii, G. hirsutum and G. barbadense genomes encoded 14, 14, 27 and 29 GAPDH proteins, respectively, suggesting that there has been no gene loss or chromosomal rearrangement in the GAPDH gene family as a result of chromosome doubling or the evolution of cotton. We named the 27 upland cotton GAPDH1 to Gh_GAPDH27 proteins according to the order of their chromosomal positions. The open reading frames (ORFs) of the upland cotton GAPDH family genes are 858–1,368 bp in length and encode proteins containing 285–455 amino acid residues, which is less variable than the ORF sequence. The relative molecular masses ranged from 31.13 to 48.65 kDa, and the theoretical isoelectric points ranged from 6.32 to 10.18, indicating that the physicochemical properties of the GAPDH family genes did not differ significantly. Subcellular localization analysis of the proteins showed that eight localized to the chloroplasts, 14 to the cytoplasm and five to the cytoplasm/mitochondria (Table 1).

To investigate the genomic distribution of the upland cotton GAPDH genes on chromosomes, we studied the chromosomal locations of the Gh_GAPDH family members. The results showed that the 27 Gh_GAPDH genes were distributed on 17 upland cotton chromosomes (Fig. 1), with subgroup A containing 15 GAPDH genes and subgroup D containing 12 GAPDH genes. Previous studies have suggested that G. arboreum and G. raimondii are donor species of the G. hirsutum A and D subgenomes, respectively, although the number of GAPDH genes in the A subgenome of upland cotton is one greater than the number of GAPDH genes in G. arboreum, and the number of GAPDH genes in the D subgenome is two less than the number of GAPDH genes in G. raimondii (Hu et al., 2019). This difference suggests that the GAPDH genes may have been relatively conserved during the evolution of upland cotton and that no tandem duplication or segmental duplication events have occurred in this species, which further suggests that the loss of a particular GAPDH gene may have occurred during evolution of the upland cotton GAPDH genes.

Evolutionary analysis of the GAPDH genes in cotton

To further understand the evolutionary relationships of the upland cotton GAPDH gene family, we constructed an evolutionary tree using the protein sequences of the GAPDH genes from G. arboreum, G. raimondii, G. hirsutum and G. barbadense. In addition, 84 cotton GAPDH protein sequences were classified into three groups, GAPA/B, GAPC and GAPCp, according to isoform classification (Marri et al., 2005; Anoman et al., 2015; Mirko et al., 2013), which were located in Group 1, Group 3 and Group 2, respectively. Upland cotton contained eight GAPA/B genes, 14 GAPC genes and five GAPCp genes (Fig. 2). The number of GAPDH genes in each subgroup was essentially twice as high in G. hirsutum and G. barbadense as in G. arboreum and G. raimondii. This result was consistent with those of the previous analysis and in line with the evolutionary relationships of cotton, which suggests that GAPDH family genes have remained relatively conserved during the evolution of cotton. Although there were relatively few members in Group 2, they have existed throughout the evolution of cotton, suggesting that they may have played an important role in the development of cotton and the response to adverse stress conditions.

Based on gene numbers, chromosome positions and phylogenetic tree analysis, GAPDH was found to be conserved during the evolution of cotton. To investigate the evolutionary relationships of the cotton GAPDH family in depth, we selected upland cotton as the core species and constructed a covariance matrix of G. hirsutum with G. arboreum, G. raimondii and G. barbadense using MCScanX software (Fig. 3). G. hirsutum shared 135 collinear gene pairs with G. arboreum, 123 collinear gene pairs with G. raimondii, and 271 collinear gene pairs with G. barbadense, and G. hirsutum itself exhibited 230 collinear gene pairs. These results suggest that upland cotton GAPDH has remained relatively conserved during the evolution of cotton.

Evolutionary tree, gene structure and motif analysis of the GAPDH genes in upland cotton

Evolutionary tree, gene structure and motif analyses based on the full-length sequences, CDSs and protein sequences of the upland cotton GAPDH genes were carried out (Fig. 4). The upland cotton GAPDH members were divided into three subgroups based on the results of evolutionary tree analysis. This finding suggests that the main difference between Groups 1 and 2 may be due to the difference in intron length between them, while the difference in Group 1 relative to Groups 2 and 3 is mainly due to the different positions of the motif. Most members of the same subgroup showed similar motifs, lengths and structures, indicating functional similarity within each group. In Group 1, the amino acid sequences of the four upland cotton GAPB members (Gh_GAPDH2, Gh_GAPDH6, Gh_GAPDH16 and Gh_GAPDH19) were longer than those of the four GAPA members (Gh_GAPDH9, Gh_GAPDH11, Gh_GAPDH21 and Gh_GAPDH23). In Group 2, Gh_GAPDH1 contained 2 motifs 7, and Gh_GAPDH26 lacked motif 4 and had no UTR. Each subgroup of protein sequences was highly conserved, but there were significant differences between the groups, especially between the sequences of Group 1 and those of Groups 2 and 3.

Analysis of cis-acting elements in the promoters of upland cotton GAPDH genes

Transcription factors (TFs) can regulate plant growth and development and responses to adverse stresses by regulating gene expression (Inoue & Horimoto, 2017). To further investigate the possible transcriptional regulatory mechanisms of the upland cotton GAPDH genes, the cis-acting elements in the 2,000 bp promoter sequences upstream of the 27 Gh_GAPDH genes were analyzed. The results showed that these cis-acting elements were mainly divided into hormone and stress response elements (Fig. 5). Most of the Gh_GAPDH gene promoters contained abscisic acid, methyl jasmonate (MeJA), gibberellin, low temperature, growth hormone and v-myb avian myeloblastosis viral oncogene homolog (MYB) response elements. Each upland cotton GAPDH gene promoter contained different numbers and types of cis-acting elements. The largest number of cis-acting elements in the Gh_GAPDH gene were related to ABA and MeJA, which suggests that the upland cotton GAPDH family may be mainly dependent on ABA and MeJA.

Different tissue-specific expression patterns of GAPDH genes in upland cotton

The expression levels of different GAPDH genes in different tissues of upland cotton (calycle, leaf, petal, pistil, root, stamen, stem and torus) were analyzed (Fig. 6). The results showed that Gh_GAPDH2, 6, 9, 11, 19, 21, 23 and 24 were highly expressed in the calycle; Gh_GAPDH7, 8 and 18 were highly expressed in the petal; Gh_GAPDH13, 20 and 25 were highly expressed in the root; Gh_GAPDH10, 15, 22 and 27 were highly expressed in the stamen; Gh_GAPDH4 and 17 were highly expressed in the stem; Gh_GAPDH1, 4 and 17 were highly expressed in the leaf; and Gh_GAPDH2, 9, 16 and 21 were highly expressed in the torus. The above results suggest that the gene expression patterns of members of the Gh_GAPDH gene family are not only specific but are also related to complex functions.

Roles of the upland cotton GAPDH genes in gene expression under different abiotic stresses

Previous studies have shown that GAPDH can play a role in the response to abiotic stresses through nonglycolytic functions (Munoz et al., 2009, 2010; Guo et al., 2012; Zhao et al., 2023). Thus, it is important to investigate the expression pattern of the GAPDH genes under different stresses. We therefore analyzed the expression patterns of the GAPDH family genes after exposure to cold, heat, salt and drought stresses (Fig. 7). The expression levels of 11 genes (Gh_GAPDH2, 6, 7, 9, 11, 12, 16, 18, 19, 21 and 23) were higher than those of the other genes after cold treatment, so these 11 GhGAPDH genes may be more responsive to cold stress (Fig. 7A). After heat treatment, 10 genes (Gh_GAPDH6, 7, 9, 11, 12, 16, 18, 19, 21 and 23) tended to be expressed at high levels, but the expression levels of these 10 Gh_GAPDH genes gradually decreased as the heat treatment time increased (Fig. 7B). After salt stress treatment, nine genes (Gh_GAPDH1, 3, 4, 5, 14, 15, 17, 26 and 27) were expressed at low levels on the first branch (Fig. 7C-1), and 11 genes (Gh_GAPDH2, 6, 7, 9, 11, 12, 16, 18, 19, 21 and 23) were expressed at high levels on the third branch (Fig. 7C-3). These 11 highly expressed Gh_GAPDH genes showed a response to salt invasion at 1 h, when their expression levels reached a maximum, but their expression levels began to decline within 6 h as the duration of salt treatment increased (Fig. 7C). However, after drought stress (Fig. 7D), changes in the expression of all but six genes (Gh_GAPDH8, 15, 16, 17, 18 and 26) were observed. The expression of 11 genes (Gh_GAPDH1, 2, 4, 9, 10, 11, 19, 21, 22, 25 and 27) was upregulated and reached a maximum in the late stage of stress, indicating that these 11 upland cotton GAPDH genes are induced by drought stress and may accordingly play a role in the response of upland cotton to drought stress.

In brief, four genes (Gh_GAPDH9, Gh_GAPDH11, Gh_GAPDH19 and Gh_GAPDH21) were highly expressed simultaneously under cold, heat, salt and drought-induced abiotic stresses, indicating significant variability in gene expression following different stress treatments. Among these genes, the expression levels of the Gh_GAPDH11, Gh_GAPDH19 and Gh_GAPDH21 genes reached their maxima in the late stage of drought stress, while the Gh_GAPDH9 gene was upregulated from 6 h of drought stress onward, and its maximum expression occurred at 12 h. This result suggests that the Gh_GAPDH9 gene may play an important role in the early defense response when upland cotton is subjected to drought stress. Based on the above analysis, we predicted that the Gh_GAPDH9 gene might be involved in drought resistance in cotton.

Expression analysis of the Gh_GAPDH9 gene under drought stress

To further analyze the expression pattern of the Gh_GAPDH9 gene, two drought-resistant materials (Acala1517-08 and Tashigan7 hao) and two drought-sensitive materials (Xinluzhong36 and Xinluzao26) were selected and subjected to qRT‒PCR analysis. We examined the relative expression levels of the Gh_GAPDH9 gene at seven time intervals after 15% PEG6000 treatment (Fig. 8A). The results showed that the expression levels of drought-resistant and drought-sensitive materials were significantly different at 0.5, 1 and 3 h. The expression levels in the drought-resistant materials Acala1517-08 and Tashigan7 hao showed a continuous increase from 0 to 3 h and decreased after 6 h of stress. The expression of this gene was downregulated in the drought-sensitive materials Xinluzhong36 and Xinluzao26 from 0.5 to 24 h, and it was barely expressed at the beginning of the drought treatment.

The Gh_GAPDH9 gene was expressed at different levels in different tissues, and understanding its expression in different tissues provides a basis for better elucidating the molecular mechanism of the gene. The expression of the Gh_GAPDH9 gene was greatest in the drought-tolerant material at 3 h of drought stress, and it is hypothesized that cotton may respond to 15% PEG6000 stress after 3 h of exposure. Therefore, we also examined the involvement of the Gh_GAPDH9 gene in the response to drought in other tissues of cotton (Fig. 8B). The expression levels of Gh_GAPDH9 were highly significantly different only in leaves (p < 0.01), showing highly significant increases in the leaves of both drought-resistant materials (Acala1517-08 and Tashigan7 hao) and drought-sensitive materials (Xinluzhong36 and Xinluzao26) after PEG drought stress.

Functional analysis of the Gh_GAPDH9 gene in response to drought stress

In this experiment, we constructed a pTRV2-GhGAPDH9 VIGS vector. At 10 d after cotton cotyledon injection, pTRV2-CLA plants showed bleaching, confirming the success of the silencing system in this assay (Fig. 9A). The expression of Gh_GAPDH9 was significantly lower in the silenced plants than in the control plants, indicating that the Gh_GAPDH9 gene was effectively silenced (Fig. 9B). The control and silenced plants were then divided into two groups, one subjected to normal irrigation treatment (Control) and one subjected to natural drought treatment (Drought). The plants were photographed and their characteristics were recorded before drought (0 d) and during drought treatment (15 d). The results showed that there was no significant difference between the silenced plants and the control plants before stress treatment, and after stress treatment, the GhGAPDH9-silenced plants showed more severe leaf wilting than the control plants (Fig. 9C).

To analyze the physiological mechanism of Gh_GAPDH9 gene silencing in response to drought stress, we examined the MDA, Pro, chlorophyll, Chla and Chlb contents of silenced GhGAPDH9 and control plants before and after drought treatment. The results showed that the contents of these five indicators did not differ significantly in silenced Gh_GAPDH9 and control plants before stress treatment, whereas after stress treatment, silenced Gh_GAPDH9 plants showed significantly higher MDA content and significantly lower Pro, chlorophyll and Chlb contents compared with the control plants (Fig. 9D). To further determine the function of Gh_GAPDH9, we transferred it to A. thaliana. The expression of Gh_GAPDH9 in transgenic strains was also confirmed by qRT‒PCR (Fig. 9E). Compared with wild-type (WT), the leaf wilting and dryness of transgenic strains under drought stress were lower, indicating that Gh_GAPDH9 is a positive regulator of drought resistance (Fig. 9F).