To obtain accurate information on the GAD family, several datasets and multiple steps were used to search for the sequences. The genome files and protein sequences of Gossypium hirsutum (G. hirsutum) (ZJU), Gossypium barbadense (G. barbadense) (ZJU), Gossypium arboreum (G. arboreum) (CRI), and Gossypium raimondii (G. raimondii) (JGI) were downloaded from the Cotton Functional Genomics Database CottonFGD (https://cottonfgd.org/) (Zhu et al., 2017). Genome data of other seven species Arabidopsis thaliana (A. thaliana), Vitis vinifera (V. vinifera), Populus trichocarpa (P. trichocarpa), Theobroma cacao (T. cacao), Glycine max (G. max), Oryza sativa (O. sativa) and Zea mays (Z. mays) were obtained from the Ensembl Plants database (http://plants.ensembl.org/i ndex.html) (Kumar et al., 2016; Wang et al., 2020). BLAST (Basic Local Alignment Search Tool) was downloaded from NCBI (https://www.ncbi.nlm.nih.gov/). The online website Softberry (http://www.softberry.com/) was used to predict subcellular localization.

A total of five published Arabidopsis GAD sequences as the queries. The BLAST program was used to identify all candidate cotton GADs (E-value < e⁻⁶). The Hidden Markov Model (HMM) profile of the Pyridoxal_deC (PF01694 in Pfam) was downloaded and used in local searches of the datasets, all the possible members of the GAD gene family were retrieved using hmmer (version 3.3.1) (http://www.hmmer.org/). Then, the common id of genes obtained by the two methods was selected as the candidate genes. To further confirm these genes, these sequences were further verified via CD-Search Tool (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi) and Simple Modular Architecture Research Tool (SMART) (http://smart.embl-heidelberg.de/). Manually deleting sequences that do not belong to the conserved binding domain and contain incomplete C and N terminals.

The identified 10 GAD family genes sequences of upland cotton were used as probes, BLAST program was used to identify the GAD family genes in the other 7 species.

The full-length amino acid sequences of 11 plant species including G. hirsutum, G. barbadense, G. arboreum, G. raimondii, A. thaliana, V. vinifera, P. trichocarpa, T. cacao, G. max, O. sativa, and Z. mays encoded by GAD genes were aligned with the ClustalW program with the default settings, and then manually adjusted in MEGA7.0. Subsequently, the neighbor-joining (NJ) tree was constructed with 1000 bootstrap replicates using the Poisson substitution (p-distance) model with default parameters in MEGA7.0 (Kumar et al., 2016). The website EvolView (https://www.evolgenius.info/evolview) was used to decorate the obtained phylogenetic tree.

The chromosomal locations of G. hirsutum, G. barbadense, G. arboreum, and G. raimondii were plotted using TBtools software (Chen et al., 2020). The reference genome GFF3 files were downloaded from CottonFGD.

To investigate the collinearity and to analyze the syntenic relationship among the GAD family of four cotton species, the complete genome sequences of these cotton species along with genome annotation files were subjected to the MCScanX tool (Wang et al., 2012). The collinear and homologous chromosomal regions among 4 cotton species were visualized using the advanced Circos package in TB tools. Gene duplication was assessed through MCScanX. To visualize duplicated regions in the 4 species of cotton, lines were drawn between duplicated genes in Circos using TBtools (Chen et al., 2020).

To investigate the selection pressure experienced by GAD duplicated gene pairs from 4 cotton species, the rates of synonymous (Ks) and non-synonymous (Ka) substitutions along with their ratios were calculated by Ka/Ks calculator in TBtools.

We used the website MEME (http://meme-suite.org/tools/meme) to predict gene motifs, the parameters were as follows: the maximum number of motifs was 15, and the rest parameters were set by default (Bailey et al., 2009). The file of the structure domain was obtained from the CD-Search Tool. The software TBtools was used to draw the association analysis diagram of the evolutionary relationship, gene structure, domain, and motifs composition of genes.

The 2000 bp DNA sequence of the upstream region of GhGADs was obtained from the CottonFGD database (http://www.cottonfgd.org/) (Zhu et al., 2017). The predicted cis-acting elements related to abiotic stresses and plant hormones in promotor regions of the GhGADs were obtained from the PlantCARE website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) for further analysis. We used RNA-Seq data (PRJNA490626) from NCBI (National Center for Biotechnology Information) (https://www.ncbi.nlm.nih.gov/) to analyze the expression level (FPKM) of GhGADs under cold (4°C), heat (37°C), salt (0.4 M NaCl), and PEG (200 g/L) stress (Hu et al., 2019). RNA-Seq data (GSE126671) (Han et al., 2019) from NCBI to analyze the FPKM of GhGADs in different tissues under 4 mM Cd²⁺ treatment for 9 h. The heatmaps were drawn based on the FPKM of GhGADs. Finally, TBtools software was used to draw a picture containing an evolutionary tree, cis-acting elements, and a heatmap of expression levels for visual observation.

The GhGAD protein interaction network was analyzed through the STRING database (https://string-db.org/) (Zhang et al., 2022). On the basis of A. thaliana orthologs to predict the interaction of GhGAD family genes with other genes in cotton.

Han 242 was used as a material with better Cd²⁺ tolerance (Han et al., 2019). We planted Han 242 in the sand, The cultural conditions were 28°C/16 h of light and 25°C/8 h of dark cycle culture. Cotton plants were soaked in 4 mM Cd²⁺ solution at the three-leaf-one-heart stage, take the root at various periods (0, 3, 6, 9, 12, 15 h) and place in liquid nitrogen (Wang et al., 2021).

We extracted RNA and reverse transcription into cDNA as a template for qRT-PCR. The primers for qRT-PCR of GhGADs are designed on NCBI (Supplementary Table S1). According to the manufacturers protocol using TransStart Top Green qPCR Supermix (TransGene Biotech Co., LTD, Beijing, China), the qRT-PCR experiment was performed on the Bio-Rad 7500 fast fluorescence quantitative PCR platform and the experiment was carried out in three independent replicates. Actin (AY305733) was used as an internal reference gene, and the relative expression level of GhGADs was calculated using 2^−ΔΔCt, then the significance analysis was carried out in SPSS software.

To verify the function of the GAD genes, we selected a highly expressed gene GhGAD6 (GH_D01G1621). The fragments of 300 bp were designed by SGN-VIGS (https://vigs.solgenomics.net/). The fragment was ligated into the pYL156 vector. The recombinant vector was transformed into Agrobacterium tumefaciens GV3101. We injected GV3101 bacterial solution carrying control pYL156 (empty vector), pYL156:GhGAD6, pYL156:PDS (positive control), and pYL192 (helper vector) into the cotyledons of Han 242. After 24 h of dark treatment, cotton was grown in an incubator with 25°C/16 h of light and 23°C/8 h of dark cycle culture (Fan et al., 2022). Subsequent treatment with Cd²⁺ is described in 2.9.

The SPAD-502 PLUS measuring instrument was used to detect the chlorophylⅡ content in the leaves (Konica Minolta (China) Investment Ltd). After Cd²⁺ treatment, 0.1 g of sample powder mixed with at least 20 cotton plants were taken to determine the superoxide dismutase (SOD) activity by the SOD activity detection kit (Nanjing Jiancheng Bioengineering Institute, A001-3-1).

H₂O₂ was detected by diaminobenzidine (DAB) staining as described previously (Ji et al., 2020). Two leaves per plant were taken from three randomly selected plants of pYL156 and pYL156:GhGAD6 lines under Cd²⁺ stress. The leaves were placed in 1 g/L DAB staining solution and treated in the dark at 28°C for 12 h, then add 95% ethanol to decolor the leaves. A deep brown polymerization product represented the reaction between DAB and H₂O₂.

29 sequences were retrieved from 4 Gosspium, 10, 9, 5, and 5 putative GAD proteins and were detected by genome-wide identification analysis in G. hirsutum, G. barbadense, G. arboreum, and G. raimondii, respectively. The open reading frame (ORF) of all 10 GhGADs ranges from 1437 (GhGAD3) to 1512 (GhGAD8) bp. The encoded proteins range from 478 (GhGAD3) to 503 (GhGAD8) amino acids, with pI varying from 5.747 (GhGAD4) to 7.14 (GhGAD5) and MWs varying from 54.229 (GhGAD3) to 57.323 (GhGAD8) kDa. Their exons are 5 or 6 (Supplementary Table S2).

Select T. cacao that is closely related to Gossypium, and A. thaliana, O. sativa, Z. mays, V. vinifera, G. max, and P. trichocarpa that are more studied in plants. GAD family genes were identified in 7 other species, 5 in A. thaliana, 5 in O. sativa, 5 in Z. mays, 4 in V. vinifera, 8 in G. max, 2 in T. cacao, and 9 in P. trichocarpa. Then we renamed the GADs based on their location on their chromosome (Supplementary Table S3). The two tetraploid cotton species G. hirsutum and G. barbadense had twice the number of GAD family genes as the two diploid G. arboreum and G. raimondii. The GAD family genes of tetraploid cotton are significantly more than other species, indicating that cotton had undergone a large-scale expansion during its evolution.

To understand the evolutionary relationship of the GAD family, we utilized 67 protein sequences to build a phylogenetic tree from G. hirsutum, G. barbadense, G. arboreum, G. raimondii, A. thaliana, V. vinifera, P. trichocarpa, T. cacao, G. max, O. sativa, and Z. mays (Figure 1). GADs can be divided into four clades based on sequence similarity, tree topology, gene structural characteristics, and motifs in each sequence (Figure 3). The results showed that the GADs clade Ⅲ had the largest number (36), of which 6 were GhGAD, clade Ⅰ and clade Ⅱ had 10 and 14 genes, respectively, each containing 2 GhGAD family genes. Clade Ⅳ contains 7 genes, namely GmGAD5, GmGAD6, GmGAD7, VitGAD2, ZmGAD1, OsGAD2, and OsGAD4, which were less close to the genes in clade Ⅰ, Ⅱ, Ⅲ. AtGAD only exists in clade Ⅰ, Ⅲ, indicating that GADs in clade Ⅱ have different functions. CcGAD is close to the cotton GAD branch, which showed cocoa and cotton are closely related and originated from the same ancestor, consistent with previous studies (Li et al., 2014). T. cacao contains only two GADs, revealing that the important role of GADs in evolution has been amplified.

To gain a more intuitive understanding of the distribution of genes on chromosomes, we constructed physical maps of the chromosome distributions of GAD gene family members in four cotton species (Figure 2). Chromosome location analysis showed that the chromosomal locations were unevenly distributed, GADs of G. hirsutum and G. barbadense were distributed on chromosomes 1, 3, 9, and 12 of At and chromosomes 1, 2, 9, and 12 of Dt, respectively (Table 1). The chromosome distribution of GADs in G. arboreum was consistent with the At, but the locations were different. The distribution on the chromosome of G. raimondii was different from the Dt, which implied the GADs rearrangement occurred in the process of tetraploid. A GAD gene is missing at the end of chromosome 12 in Dt of G. barbadense, compared with that of G. hirsutum. We supposed that may be due to the loss of the G. barbadense gene during evolution or incomplete genome assembly.

We analyzed evolutionary relationships, motifs, domains, exons, and introns to study the conserved structure of GAD family genes (Figure 3). All GAD family genes had Pyridoxal_deC (Type II pyridoxal phosphate-dependent decarboxylase) domain, which can bind to PLP to achieve its catalytic function (Figure 3C).

They were classified according to the tree topology of the evolutionary tree (Figure 3A). The distribution patterns of exons and introns correlate with their biological functions, and their arrangement can be used to analyze evolutionary associations between members of different gene families. Interestingly, the first introns of GADs in clade Ⅱ (GhGAD5 and GhGAD10) were significantly longer than the others (Figure 3D). GAD family gene motifs were relatively consistent. GADs in clade Ⅰ (GhGAD3 and GhGAD8) lacked motif 13 at the N-terminus and motifs 12 and 14 at the C-terminus. GADs in clade Ⅱ lacked motif 12 at the C-terminus (Figure 3B). It was speculated that functional changes had occurred during evolution.

Gene duplication events are considered to play an important role in the amplification of gene families. To explore the amplification mechanism of the GAD gene family, by comparing the genomes of Ga-Ga, Ga-Gb, Ga-Gh, Gb-Gb, Gb-Gr, Gb-Gh, Gr-Gr, Gr-Ga, and Gh-Gh, a total of 114 homologous gene pairs were identified. There were 21, 21, 2, and 4 duplication GAD gene pairs identified in Gh-Gh, Gb-Gb, Ga-Ga, and Gr-Gr, these 48 paralogous gene pairs were predicted as segmental duplications according to the chromosomal location (Figure 4). There were 66 GAD gene pairs that underwent whole genome duplication (WGD), the numbers were 21, 22, and 23 in Ga-Gb, Ga-Gh, and Gb-Gr, respectively. From these results, we presumed that segmental duplication and WGD were the main reason for the evolution of the GAD gene from diploid to tetraploid.

During evolution, duplicated gene pairs may also deviate from their original functions, eventually leading to neofunctionalization (loss of original function), subfunctionalization (a division of original function), and neofunctionalization (gain of new function). To investigate the driving forces of the GAD family gene during evolution, we calculated non-synonymous substitution (Ka) and synonymous substitution (Ks) values for 164 repeated gene pairs from four Gossypium (Figure 5). The selection pressure of duplicate gene pairs can be inferred according to the ratio of Ka/Ks. It is generally believed that Ka/Ks = 1 indicates neutral selection (pseudogene), Ka/Ks < 1 indicates purification or negative selection (purification selection), and Ka/Ks > 1 indicates positive selection. There are 164 duplicate gene pairs in the GAD family genes in the four Gossypium, including Ga-Ga, Ga-Gb, Ga-Gr, Gb-Gb, Gb-Gr, Gh-Ga, Gh-Gb, Gh-Gh, Gh-Gr, and Gr-Gr. 153 (97%) duplicate gene pairs Ka/Ks ratio are between < 0.5, 5 (3%) duplicate gene pairs Ka/Ks ratio are between 0.5 and 0.99, they are GaGAD2-GbGAD2, GaGAD2-GhGAD2, GaGAD5-GbGAD4, GaGAD5-GhGAD4, GrGAD4-GbGAD9, indicating that GAD is evolving slowly and had strong purifying selection pressure with the limited functional divergence that occurred after segmental duplications and WGD (Table 2).

Through promoter analysis, we understand the response of GhGAD to hormonal and abiotic stress, which was conducive to further analysis of the regulatory network. GhGAD is related to plant hormones (ABA, MeJA, GA, IAA, SA) and various stresses (low temperature, drought, hypoxia, defense, and stress responsiveness) (Figure 6B, Supplementary Table S4). GhGAD10 contains an MYB binding site that regulates flavonoid synthesis. GAD contains many GA and MeJA regulatory elements. GhGAD5 and GhGAD10 increased at 1–6 h and decreased at 12 h after salt, heat, and drought treatments, which may be involved in the regulation of salt, heat, and drought (Figure 6C). GhGAD2 and GhGAD7 were significantly up-regulated after cold stress, and GhGAD4 and GhGAD9 were up-regulated by PEG induction.

To analyze the expression patterns of the GAD family in different tissues of cotton, we used transcriptome data to analyze the FPKM values of cotton 8 tissues (roots, stems, leaves, torus, petal, stamen, pistil, and calycle). The result showed that GADs were expressed in various tissues, and the overall expression is the highest in stamens (Figure 7). GhGAD5 and GhGAD10 were highly expressed in all tissues, and their expression was significantly higher than that of other GAD family genes, which may be necessary to maintain the normal life activities of cotton. GhGAD3 and GhGAD8 were specifically expressed in stamens, we speculate that they may play an important role in the development of stamens. GhGAD4 was highly expressed in stems, GhGAD9 was highly expressed in petals, and GAD orthologous gene pairs showed different expression patterns. GAD gene family is highly expressed in stamens (Figure 7B). These results suggest that GAD had tissue-specific expression under normal growth conditions.

To investigate GhGADs responses to abiotic stress, especially Cd²⁺, we used qRT-PCR to study the expression changes of the GAD gene in upland cotton roots under Cd²⁺ (Figure 8). The expression levels of GhGADs in clade Ⅲ (GhGAD1, GhGAD2, GhGAD4, GhGAD6, GhGAD7, and GhGAD9) increased significantly after 3 h under Cd²⁺ stress and continued high expression thereafter. The clade Ⅱ GhGAD genes (GhGAD5 and GhGAD10) decreased significantly after stress. The clade Ⅰ GhGAD genes (GhGAD3 and GhGAD8) did not change significantly. The results revealed that different clades of GAD family genes had different expression patterns in response to Cd²⁺ stress.

The expression patterns of the GhGADs in clade Ⅲ were all up-regulated in roots under Cd²⁺ treatment (Figure 9). GhGAD1, GhGAD2, GhGAD6, and GhGAD7 were significantly increased in the stem. Only GhGAD6 responded to Cd²⁺ stress in roots, stems, and leaves. Therefore, we selected GhGAD6 for further analysis.

Based on the homologous gene ATGAD1 in Arabidopsis with the highest homology to GhGAD6, an interaction network was constructed using the STRING database to analyze the function of the GAD protein. (Figure 10). ATGAD1 interacts with 4-aminobutyrate pyruvate transaminase (POP2), succinate-semialdehyde dehydrogenase (ALDH5F1), delta1-pyrroline-5-carboxylate dehydrogenase (ALDH12A1), glutamate dehydrogenase (GDH), and glutamate synthase (GLT) proteins. By analyzing the KEGG pathway of GhGAD6 from the transcriptome data, we found that GhGAD6 was mainly involved in Alanine, aspartate, and glutamate metabolism (ko00250), GAD synthesized GABA through POP2, ALDH5F1 converted GABA to succinate, and then entered the TCA cycle. ALDH12A1, GLT, and GDH were involved in the synthesis of glutamate. We speculated that GAD interacted with these proteins to respond to Cd²⁺ stress by regulating the content of glutamate and GABA.

We used VIGS experiments to verify the role of GhGAD6 under Cd²⁺ stress. The pYL156:PDS exhibited obvious chlorosis, and the relative expression level of GhGAD6 was determined by qRT-PCR, which showed a 70% decrease in pYL156:GhGAD6 than pYL156 (Figure 11B), indicating that it had a good silencing effect. After Cd²⁺ stress, cotton showed blackening of stems and veins, wilting of leaves, and pYL156:GhGAD6 was more severe than pYL156 (Figure 11A). The DAB staining showed the same result, more brown was produced at the veins of pYL156:GhGAD6 than pYL156, indicating that H₂O₂ was accumulated in pYL156:GhGAD6 after Cd²⁺ stress (Figure 11C). It can be seen that both pYL156 and pYL156:GhGAD6 are up-regulated of GhGAD6 after Cd²⁺ stress, and the up-regulation range is 3–4 times, but pYL156:GhGAD6 is still significantly lower than pYL156. Furthermore, the chlorophylⅡ and SOD content of pYL156:GhGAD6 also increased significantly after Cd²⁺ stress (Figures 11D,E).

We detected the expression of GABA transaminase (GABA-TP), GABA transporter (GAT) in the GABA shunt, and GABA receptor protein aluminum-activated malate transporter (ALMT) in plants (Figure 11F). It was found that the expression of GABA-TP, GAT, and ALMT9 decreased significantly after silencing GhGAD6, which indicated that the GABA content was reduced. After Cd²⁺ stress, the expressions of GABA-TP and GAT were up-regulated, and the up-regulated expression of GABA shunt genes coped with Cd²⁺ stress.