Identification, characterization, and expression profiles of the GASA genes in cotton

Abstract

Background

GASA (Giberellic Acid Stimulated in Arabidopsis) gene family plays a crucial role in the phytohormone signaling pathway, growth and development, and stress responses in plants. Many GASA homologs have been identified in various plants. Nevertheless, little is known about these proteins in cotton.

Results

In the current study, we identified 19, 17, 25, 33, and 38 GASA genes via genome-wide analyses of Gossypium herbaceum, G. arboreum, G. raimondii, G. barbadense, and G. hirsutum, respectively, and performed comprehensive bioinformatics and expression analyses. According to our results, 132 GASA proteins shared similar protein structures and were classified into four groups based on the phylogenetic tree. A synteny analysis suggested that segmental duplication was a key driver in the expansion of the GASA gene family. Meanwhile, the cis-element and protein interaction analyses indicated that GhGASA proteins play a significant role in the hormone responses. Transcriptomic and qRT-PCR (Quantitative real time-polymerase chain reaction) analyses revealed diverse expression profiles of the GhGASA genes in different organs under abiotic stresses, indicating that some GhGASA genes possibly participate in fiber development and abiotic-stress responses.

Conclusions

The GASA genes in cotton were systematically identified and analyzed for the first time in this paper, and it suggested that the GASA genes are important to the development and growth of cotton. These results will support future exploration of the functions of GASA genes in cotton.

Introduction

GASA (Giberellic acid stimulated in Arabidopsis) is a CRP (cysteine-rich peptide) protein (Silverstein et al. 2010) and a low-molecular-weight peptide that is widespread in the plant world. Numerous studies have illustrated that GASA proteins are important for plant growth and development, but much remains unknown. The GAST1 gene was the first identified member of the GASA gene family, found in the gib1 tomato mutant (Shi et al. 1992). Previous studies on the GASA gene family have been conducted in A. thaliana (Aubert et al. 1998; Roxrud et al. 2007), Moso Bamboo (Hou et al. 2018), apple (Fan et al. 2017), rice (Muhammad et al. 2019), wheat (Lv et al. 2018), Glycine max (Ahmad et al. 2019), Solanum tuberosum (Nahirñak and Rivarola 2016), grapevine (Ahmad et al. 2020), and Sorghum bicolor (Filiz and Kurt 2020). Almost all GASA protein domains have a C-terminal region of approximately 60 amino acids consisting of 12 cysteine residues (Aubert et al. 1998). Studies have shown that missing or alterating crucial peptide residues may lead to a loss of biological function of GASA domain (Sun et al. 2013). This indicates that the evolution of cysteine-rich regions is relatively conservative, and that these regions possesse important biological activities (Sun et al. 2013).

Several studies have shown that GASA proteins influence processes of phytohormone response and plant growth and development, including signal transduction, seed development, floral and root development, and flowering time. For instance, in rice, OsGSR1 can influence brassinosteroid signaling by interacting with DIM/DWF1 (DIMINUTO/DWARF1), a brassinosteroid synthetase (Wang et al. 2010), and OsGASA1/7 genes are important in the brassinolide (BR) and gibberellin (GA) signaling networks (Lee et al. 2015; Boonpa et al. 2018). The expression levels of several MdGASA genes were up- or downregulated under GA and 6-benzylaminopurine (6-BA) treatments in apple (Fan et al. 2017). In Arabidopsis, AtGASA2/3/5 and AtGASA14 are involved in flower induction and participate in the abscisic acid (ABA) signaling response (Roxrud et al. 2007). Moreover, maize ZmGSL genes exert different functions in GA₃-regulated root growth (Zimmermann et al. 2010). Furthermore, the expression of potato GASA genes is affected by GA and ABA hormones, and is highly induced during the seedling stage and floral formation (Nahirñak and Rivarola 2016). Notably, TaGASR7 is associated with grain length in wheat (Dong et al. 2014), and OsGASA genes increase grain size and length in rice (Muhammad et al. 2019). VvGASA2 and VvGASA7 are involved in seed development in grapevines (Vitis vinifera L.) (Ahmad et al. 2020). These results show that GASA genes have a great potential for increasing the crop yield and fruit breeding. However, many reports have indicated that various GASA gene family members serve antagnistic functions. For example, AtGASA4 facilitates flowering while AtGASA5 delays or inhibits flowering in Arabidopsis (Aubert et al. 1998; Roxrud et al. 2007; Zhang et al. 2009), and GEG (Gerbera homolog of GAST1 gene) delays petal development while PRGL (proline-rich and GASA-like) promotes petal development in Gerbera hybrida (Kotilainen et al. 1999; Peng et al. 2010).

In other researches, some GASA proteins have been shown to have important impacts on biotic- and abiotic-stress resistance. GsGASA1 isolated from Glycine soja was shown to be involved in inhibiting root growth via the accumulation of DELLA proteins under chronic cold stress (Li et al. 2011). Arabidopsis GASA4 (Ko et al. 2007) and GASA5 play a role in response to heat stress (Zhang and Wang 2011), and the overexpression of Beechnut FsGASA4 in Arabidopsis improves tolerance to abiotic stress (Alonso-Ramírez et al. 2009). Another report noted that GASA14 can modulate the accumulation of reactive oxygen species (ROS) to regulate abiotic-stress resistance (Sun et al. 2013). In disease resistance, HbGASA genes in rubber plants show high expression levels when faced with fungi (An et al. 2018). Indeed, snakin1, − 2, and − 3 in Solanum tuberosum responded similarly when inoculated with bacterial fungal pathogens (Nahirñak and Rivarola 2016). Meanwhile, Snakin-1 and Snakin-2 proteins exert antibacterial activity in potato (Segura et al. 1999; Berrocal-Lobo et al. 2002; Almasia et al. 2008; Kovalskaya and Hammond 2009; Balaji and Smart 2012). Moreover, CaSn, which is a GASA/snakin protein, is involved in halting the intrusion of root-knot nematodes (Meloidogyne spp.) in pepper (Mao et al. 2011).

Cotton is an important fiber crop and is the main raw material used to produce lint, along with cooking oil and other products (Zhang et al. 2008). In recent years, many cotton gene families have been reported and analyzed, including INDETERMINATE DOMAIN genes (Ali et al. 2019), plant U-box genes (Lu et al. 2020), DNA demethylases (Yang et al. 2019b), OSCA (Hyperosmolality-gated calcium-permeable channels) gene (Yang et al. 2019a), and VEFS-box (Ma et al. 2020). However, there has been little research on GASA genes. An earlier study identified GhGASL1–GhGASL7 from the cDNA libraries of upland cotton and suggested that GhGASL genes are involved in fiber development and GA signaling responses (Liu et al. 2013). Overall, most GASA genes in cotton are yet to be identified and researched. Here, we identified the GASA genes in five cotton species using the latest genomic database and bioinformatics methods, and illustrated their chromosomal locations, chemical characterization, phylogenetic relationships, gene structures, synteny, cis-elements, and expression patterns, predicted their potential functions. Our findings will further expand the available information on GASA gene family members and provide a theoretical basis for further studies on the functions of GASA genes in cotton.

Materials and methods

Identification of the GASA proteins in Gossypium spp.

First, we used the AtGASA protein sequences as queries to search in the cotton genome database (https://cottonfgd.org/; https://www.cottongen.org/) (Zhu et al. 2017; Yu et al. 2014), using the BLAST Similarity Search tool. Second, proteins containing the GASA domain (PF02704) in the Gossypium spp. database were detected via the hidden Markov model (HMM). The complete GASA domains of the protein sequences were examined using Batch CD-Search tools (https://www.ncbi.nlm.nih.gov/Structure/cdd/cdd.shtml), SMART (http://smart.embl.de/smart/set_mode.cgi?NORMAL=1#), and Pfam (http://pfam.xfam.org/). Finally, conserved domains of all candidate protein sequences were verified via manually removing incomplete domain sequences.

Physicochemical characterization, chromosomal location, and sequence alignments

The ExPASy program was used to predict the characteristics of the GASA protein (Artimo et al. 2012). The subcellular locations of the GASA proteins were predicted using an online tool (http://cello.life.nctu.edu.tw/) (Yu et al. 2010). Information on the chromosomal locations of the GASAs was acquired from the Cottongen database (https://www.cottongen.org/), and the gene locations were plotted on the chromosomes using MapChart software (Voorrips 2002) and further embellished using Adobe Illustrator CC 2018. ClustalW in the SeaView software was used to align the GASA protein sequences, while the sequence logos were further created using the WebLogo platform (Crooks et al. 2004), and GeneDoc software was used to color the conserved amino acids. The 3D structures of the GASA proteins were predicted using PHYRE server v2.0 (http://www.sbg.bio.ic.ac.uk/phyre2/html/page.cgi?id=index, which is an online tool.

Phylogenetic tree, gene structure, and promoter region analyses

Phylogenetic trees were constructed using the neighbor-joining method under 1 000 replicates bootstrap test by MEGA 7.0 software (Kumar et al. 2016). The online GSDS tool (http://gsds.cbi.pku.edu.cn) was used to visualize the exon–intron structures of the GASA genes. Prediction of the cis-acting elements was performed using the PlantCARE online tool on the upstream sequence of each gene (1.5 KB) (Lescot et al. 2002), and the image were ultimately created using TBtools software (Chen et al. 2020).

Synteny analysis and protein interactions of GASA genes in Gossypium hirsutum

The synteny and collinearity of GASA genes were determined using the Multiple Collinearity Scan toolkit (MCScanX) (Wang et al. 2012a). We then used the Circos software to visualize the syntenic relationship of the duplicated genes (Krzywinski et al. 2009). The protein interaction database is only available for G. raimondii, which was obtained from the String website (https://string-db.org/), then the CDS sequence of the GASA genes in upland cotton was blast against CDS sequence of G. raimondii to obtain corresponding protein interactions of upland cotton. Related protein information was based on an available database (https://www.uniprot.org/). The Cytoscape software was used to map the protein–protein interactions (Shannon et al. 2003).

Plant materials and treatments

Gossypium hirsutum cv. TM-1 was planted in field in Anyang, Henan, China. The roots, stems, leaves, anthers, and petals of the plants were harvested for RNA extraction. The ovules of the flowers were collected in liquid nitrogen for RNA extraction 3 days before flowering, as well as at 0 days, 1 day, 3 days, and 5 days after flowering. Seedlings with two cotyledons were moved to a light incubator at 4 °C. The whole plants were collected at 1 h, 3 h, 6 h, 12 h, and 24 h in the 4 °C treatment. The control seedlings were planted under consistent natural conditions. We immediately froze all samples in the liquid nitrogen and stored them at − 80 °C. Three independent biological experiments were conducted for each sample.

Transcriptome data analysis and quantitative RT-PCR (qRT-PCR) for GASA genes in Gossypium hirsutum

To analyze the expression profiles of GhGASA genes in different tissues and at different developmental stages, we retrieved the FPKM (fragments-per-kilobase-per-million) expression data from CottonFGD (https://www.cottonfgd.org) (Zhu et al. 2017). The TBtools software (Chen et al. 2020) was used to assemble transcripts for each RNA-seq analysis. Heatmap charts were generated based on the FPKM values. The gene-specific primers used in qRT-PCR (Quantitative real time-polymerase chain reaction) were listed in Additional file 1: Table S1. The PCR conditions were as follows: preheating for 30 s at 95 °C; 35 cycles for 30 s at 95 °C, 60 s at 60 °C, and 10 s at 72 °C for 10 s, and, finally, 2 min at 72 °C for extension. The relative gene expression levels were calculated according to the 2^−ΔΔC_T method (Schmittgen and Livak 2008).

Results

Identification and protein features of GASA genes in Gossypium spp.

To identify the GASA genes in Gossypium spp., we used the Blast by Similarity Search Tool and the Hidden Markov Model (HMM) alignment search on the protein sequence databases of G. herbaceum, G. arboreum, G. raimondii, G. barbadense, and G. hirsutum, respectively. Then, the Batch CD-Search tools, Pfam and SMART were used to verify the obtained sequences. Finally, the numbers of GASA genes confirmed in G. herbaceum (GheGASAs), G. arboreum (GarGASAs), G. raimondii (GrGASAs), G. barbadense (GbGASAs), and G. hirsutum (GhGASAs) were 19, 17, 25, 33, and 38, respectively. All genes were named based on their locations on the chromosomes. Meanwhile, previously reported GhGASL genes were renamed according to their homologs. The predicted sequences of GASA genes were deposited in GenBank (accession numbers MW429289-MW429326, MW450702-MW450776).

Detailed information about cotton GASA genes and proteins characteristics are listed in Table 1 and Additional file 2: Table S2, including the length of gene, molecular weight (MW) of the protein, isoelectric point (pI) of the protein, subcellular localization of the protein, length of the CDS sequence. The analysis results showed that the GarGASA17 was the shortest protein, in G. arboreum, while GbGASA19 was the longest, in G. barbadense. The average MW of GhGASAs was 11 776.78 Da, and the average MW of GheGASAs, GarGASAs, GrGASAs, and GbGASAs were 15 972.93 Da, 11 349.69 Da, 12 350.80 Da, and 14 071.10 Da, respectively. The average isoelectric point (pI) of GheGASAs, GarGASAs, GrGASAs GhGASAs and GbGASAs were 8.65, 8.65, 8.85, 8.68, and 8.70, respectively, and they showed similar statistical results. The pI values of the five cotton species’ proteins ranged from 6.44 to 10.60. The pI of GarGASA11, GheGASA12, and GhGASA32 were the lowest, 6.44, while that of GrGASA3 was the highest, 10.60. The function of a protein is usually closely linked to its subcellular positions. The predicted subcellular locations suggested that most of the GASA proteins in cotton are located in the extracellular. Only a few proteins are located in nuclear and plasma membrane.

Table 1 Information of the GASA genes in G. hirsutum

Based on the predicted 3D protein structures, α helices, β sheets, extended strands, and random coils were found in GhGASA proteins. According to our analysis, all GhGASA proteins have coils which is a flexible structure, at least three α helices, while β sheets are uncommon (Additional file 3: Fig. S1).

Sequence analysis, chromosomal localization, and phylogenetic tree of cotton GASA genes

To verify the GASA domain, the full GASA protein sequences were aligned using SeaView and GeneDoc software. All the GASA proteins shared similar structures: an N-terminal signal peptide containing 18 to 30 amino acids, a variable peptide region, and a highly conserved C-terminal region consisting of about 60 amino acids with 12 cysteine residues (Fig. 1, Additional file 4: Fig. S2, Additional file 5: Fig. S3). Based on the publicly available information on the cotton genome, chromosomal maps were drawn for the 132 GASA genes from G. herbaceum, G. arboretum, G. barbadense, G. hirsutum, and G. raimondii using MapChart software (Additional file 6: Fig. S4, Additional file 7: Fig. S5). In total, 131 GASAs were dispersed across the chromosomes, except for one GASA gene of G. arboretum that was located in the unmapped scaffold. Notably, longer chromosomes did not necessarily contain more members of the GASA gene family, which indicates that the number of GASA genes on each chromosome is independent of its length. This result showed that the GASA protein-coding genes of G. hirsutum, as in other plants (Fan et al. 2017; Muhammad et al. 2019; Ahmad et al. 2020), are not evenly distributed on 26 chromosomes. To study the evolutionary relationships of the GASA gene family in cotton, a phylogenetic tree was constructed using the neighbor-joining method for 173 GASA proteins from G. arboreum (17 GarGASAs), G. herbaceum (19 GheGASAs), G. raimondii (25 GrGASAs), G. hirsutum (38 GhGASAs), G. barbadense (33 GbGASAs), Malus domestica (26 MdGASAs), and A. thaliana (15 AtGASAs) (Fig. 2). These results showed that all GASA genes could be divided into four groups (G1, G2, G3, and G4) based on topological consistency.

Fig. 1

Features of the GhGASA protein domains. a GhGASA protein structure analysis; the 12 cysteines are shown in the red column. b Multiple sequence alignments of the GhGASA conserved-domain. c Sequence logo of the GhGASA conserved-domain. Alphabets in each sequence represents amino acids

Fig. 2

Phylogenetic relationships among the GASA proteins. The four subgroups of GASAs are shown using four different colors

Gene structure analysis of GASA members in G. hirsutum

To further study the evolutionary relationships of GhGASA proteins, a phylogenetic tree was reconstructed (Additional file 8: Fig. S6A). The coding sequences within the corresponding genomic sequences were aligned for a comparative analysis of the intron–exon structures (Additional file 8: Fig. S6B). The structures of GASA genes were divided into four groups. Each member of G4 contained two exons, while the number of exons in G1, G2, and G3 ranged from 2 to 4, 2 to 3, and 1 to 4, showing that the structure of G4 group members is more conserved than that of other subgroups.

Synteny analysis of GASA genes in G. hirsutum

BLAST and MCScanX were used for genomic synteny analysis to further investigate the expansion mechanism of the GASA gene family in G. hirsutum. As shown in Additional file 9: Fig. S7, two genes were tandemly duplicated (GhGASA5 and GhGASA6), while 20 pairs of segmental duplication were observed in the 25 genes. Thus, compared with tandem duplication, segmental duplication has played a dominant role in the evolution of GASA gene family. To explain the evolutionary relationship of GASA genes, the synteny analysis was conducted using the genomes of different cotton species (Fig. 3). G. herbaceum (A₁) and G. arboreum (A₂), G. arboreum (A₂) and G. hirsutum A-subgenomes (GH_At) had 15, 15 collinear gene pairs, respectively (Fig. 3a). We found that most chromosomes showed evidence of whole-genome duplication during evolution, while a few chromosomes (A01, A02) underwent a segmental swap. This indicates that the A-genomes were similar and stable across different cotton species. As shown in Fig. 3b, there were 12 and 16 collinear gene pairs in the G. hirsutum D-subgenome (GH_Dt) and G. raimondii (D₁), and G. barbadense D-subgenome (Gb_Dt) and G. raimondii (D₁). In addition, the collinearity relationships of D-genome of between tetraploid cotton and G. raimondii were irregular, except for the chromosome D03 and D13, which experienced genome duplication. However, the D-genome of different cotton species remained highly consistent in the collinearity blocks, which is consistent with the fact that the ancestor of the D-genome of tetraploid cotton was originated from G. raimondii (Wang et al. 2012b; Paterson et al. 2012; Li et al. 2015).

Fig. 3

Distribution and gene duplications of GASA genes. a Synteny analysis of GASA genes between the G. herbaceum (A₁), G. arboreum (A₂), and G. hirsutum A subgenome (GH_At). b Synteny analysis of GASA genes between G. barbadense D subgenome (Gb_Dt), G. raimondii (D₁) and G. hirsutum D subgenome (GH_Dt). Purple lines and blue lines highlight the syntenic GASA gene pairs

Promoter region analysis of GhGASA genes

Cis-element contribute to the regulation of gene expression. To study the regulatory mechanisms of GhGASA genes, we selected the 1.5 kb upstream of the start codon according to the G. hirsutum genome sequence and analyzed the potential cis-elements (Additional file 10: Fig. S8, Additional file 11: Fig. S9). Cis-elements related to phytohormones were detected in 38 GhGASA gene promoters, including TGACG-motif, TATC-box. P-box, ABRE, TGA-element, AuxRR-core, ERE and TCA-element. Among these, most cis-elements were related to abscisic acid, gibberellin, and salicylic acid. In most genes, cis-elements associated with different types of stress response (LTR and STRE) and disease resistance (TC-rich repeats, W-box, and WUN-motif) were found. The cis-elements involved in zein metabolism regulation (O2-site) and seed-specific regulation (RY-element) were also found in these parts of the promoters. Moreover, cis-elements associated with anaerobic respiration (ARE) and light response (GT1-motif, TCT-motif, Box 4, GATA, G-BOX, and AE-box) were detected in the promoters of most of the GhGASA genes. The results for these cis-acting regulatory elements indicated that GhGASA genes may be regulated by a variety of plant hormones and respond to various stresses in upland cotton.

Expression profiles of GhGASAs in different tissues of upland cotton

To explore the possible biological functions of GhGASAs, published transcriptomic data were downloaded from CottonFGD (https://www.cottonfgd.org) to analyze the expression of 38 GhGASA genes in different tissues of cotton (Additional file 12: Fig. S10). Most of the GASA genes in G. hirsutum were found to have no expression or extremely low expression in multiple tissues, including anther, filament, pistil, bract, sepal, petal, torus, root, leaf, stem, fiber, and ovule. GhGASA9 and GhGASA31 were similar, with expression observed in almost all tissues and showed higher transcriptional levels than those in other tissues in pistil, petal, and torus. And the transcriptional levels of GhGASA14 in leaf, GhGASA4, GhGASA25, and GhGASA26 in stem; and GhGASA8 and GhGASA23 in root and stem were significantly higher than respective expression levels in other tissues. Moreover, the GhGASA38 gene were only expressed in the root. These results show that the above-mentioned GhGASA genes may have special functions in different tissues. According to Additional file 12: Fig. S10, the expression levels of GhGASA3, GhGASA23, and GhGASA24 continuously increased from − 3 to 5 days post anthesis (DPA) fiber, and the corresponding expression levels of 3 and 5 DPA fiber were far higher than − 3 and 1 DPA fiber, suggesting that GhGASA3, GhGASA23, and GhGASA24 play an important role in the early elongation of cotton fiber cell. Moreover, GhGASA3, GhGASA24 in 10 DPA_fiber, 15 DPA_fiber, and 20 DPA_fiber was highly expressed, GhGASA19 and GhGASA23 in 10 DPA_fiber, 15 DPA_fiber were highly expressed, suggesting that the above-mentioned genes have an important influence on the development of fiber.

In order to verify the transcriptome data, the GhGASA4, GhGASA8, GhGASA9, GhGASA13, GhGASA14, and GhGASA26 genes were selected to analyze the expression patterns in anther, stem, root, leaf, and petal by qRT-PCR (Fig. 4a). The results showed that four genes (GhGASA4, GhGASA8, GhGASA13, and GhGASA26) displayed higher expression levels in stem and lower expression levels in anther. GhGASA9 in petal, GhGASA13 in root, GhGASA14 in leaf showed higher transcription levels. In addition, the GhGASA3, GhGASA8, GhGASA13, GhGASA16, GhGASA23, and GhGASA24 genes were used to analyze the expression patterns in the ovule (− 3, 0, 1, 3, 5 DPA) by qRT-PCR (Fig. 4b). GhGASA8 and GhGASA23 showed continuous increasing expression levels from − 3 to 3 DPA ovule. GhGASA16 showed continuous increasing expression from − 3 to 1 DPA ovule, and then decreased from 1 to 5 DPA ovule. The expression level of GhGASA3 and GhGASA24 in 3 DPA ovule showed higher expression. GhGASA13 displayed differential expression levels from − 3 to 5 DPA ovules, whereas the highest expression level was in 1 DPA ovule. These results were consistent with the RNA-seq data (Zhu et al. 2017), indicating that the data from CottonFGD are reliable.

Fig. 4

The expression profiles of 38 GhGASA genes by RT-qPCR. a Expression levels of GhGASA4, GhGASA8, GhGASA9, GhGASA13, GhGASA14, and GhGASA26 in anther, petal, leaf, stem, and root. b Expression levels of GhGASA3, GhGASA8, GhGASA13, GhGASA16, GhGASA23, and GhGASA24 in − 3, 0, 1, 3, and 5 DPA ovule

Expression profiles of GhGASAs under abiotic stress in upland cotton

Transcriptomic data downloaded from CottonFGD (https://www.cottonfgd.org) were used to analyze the expression profiles of GhGASA genes under abiotic stress to study their potential functions in G. hirsutum. As shown in Fig. 5a, more than half of the genes had relatively no or very low expression levels under different stress conditions in cotton leaves. The GhGASA9 and GhGASA31 genes were similar in response to multiple stress treatments, especially the 6 h cold treatment, 6 h heat treatment, 6 h polyethylene glycol (PEG) treatment, and 1 h and 6 h salt treatments, in which presenting higher expression levels of GhGASA9 and GhGASA31 thanin the control treatment. The expression of GhGASA25 was downregulated under 3 h cold, 3 h PEG, and 3 h salt treatment and upregulated under 1 h PEG and 1 h salt treatment. The expression levels of GhGASA4 and GhGASA26 were similar, both showed high expression under 3 h and 12 h cold treatments, and had significantly increased expression levels from 3 h to 6 h under salt treatments. Moreover, the GhGASA12 under 12 h PEG treatment and GhGASA18 and GhGASA37 under 3 h cold and 6 h salt treatments presented higher expression levels than the control treatment. These results indicate that these genes  may play roles in plant resistance to abiotic stress and are potential candidate genes for further research.

Fig. 5

a Expression profiles of GhGASAs in cotton leaves under abiotic stress. b Expression analysis of GhGASA4, GhGASA9, GhGASA18, and GhGASA25 in cotton leaves under cold treatment based on qRT-PCR

Based on previous researches and transcriptome data, GhGASA4, GhGASA9, GhGASA18, and GhGASA25 were selected for qRT-PCR (Fig. 5b). As shown in Fig. 5b, the expression levels of GhGASA4 and GhGASA18 significantly increased after 3 h of cold treatment and significantly decreased over the next three hours. GhGASA25 was poorly expressed under the cold treatment, showing a significant downregulation level at 3 h. GhGASA9 showed differential expression levels under cold treatment, whereas downregulated at 1 h, 3 h, and 12 h and upregulated at 6 h and 24 h.

Protein interactions of GhGASA genes

To predict the functions of GhGASA genes, we mapped the protein interaction network. Meanwhile, the available database (https://www.uniprot.org/) was utilized. Based on subcellular localizations, GhGASA13 and GhGASA35 were selected to the analysis of protein interaction by Cytoscape software (Fig. 6). Our results suggest that the GhGASA proteins may interact with several other proteins, including the TPR repeat-containing thioredoxin TTL1, the ATPase ASNA1 homolog, the ATPase GET3, the DnaJ homolog subfamily C member 7, the serine/threonine protein kinase STY8, and the serine/threonine protein kinase BRI1-like 2.

Fig. 6

Interaction network of GhGASA proteins using the Cytoscape software

Discussion

GASA genes were shown to be involved in plant growth, development, and physiological processes, including seed development, floral organ development, flowering time, phytohormone signaling transduction, stress response, defense responses against pathogens, and plant reproductive processes. Although some GASA genes have been identified and analyzed, GASA genes have not been systematically studied in cotton. The genomes of two allotetraploid species (Hu et al. 2019) and three diploid species (Du et al. 2018; Huang et al. 2020) were improved upon or published, which provided more precise information for the whole-genome analysis of cotton gene family. In this study, 19, 17, 25, 33, and 38 candidate GASA genes were systematically identified in G. herbaceum, G. arboretum, G. raimondii, G. barbadense, and G. hirsutum, respectively.

Expansion analysis of GASA genes in cotton

According to previous reports, allotetraploid cotton was originated from the intergenomic hybridization of the A-genome diploid species and the D-genome diploid species, followed by chromosomal doubling (Paterson et al. 2012). The number of most gene family members in the tetraploid species is close to the sum of the gene family members in the diploid progenitor. However, in this study, the number of GASA genes in tetraploid and diploid cotton did not present the expected relationship. Recent studies showed that all existing A-genomes may have originated from a common ancestor, A₀ (Huang et al. 2020). The number of GASA genes of G. herbaceum (A₁), G. arboreum (A₂), G. barbadense A-subgenome and G. hirsutum A-subgenome hasis 19, 17, 17 and 18, respectively. Moreover, most chromosomes of the A-genome in the cotton species had undergone whole-genome duplication, according to the synteny analysis, indicating tha GASA genes in A-genomes have been highly conserved during evolution. This demonstrated, again, that an unknown A₀ is the possible common origin of the A-genomes. Past research has shown that the D-genomes of the allotetraploid cotton species are derived from the genome of G. raimondii (Wang et al. 2012b; Paterson et al. 2012; Li et al. 2015). However, the number of GASA genes in the G. raimondii (D₁) and G. hirsutum D-subgenome and the G. barbadense D-subgenome was 25, 20, and 16, respectively. Moreover, the collinear relationships of the D-genome in cotton species were shown to be unordered. Hence, we argue that the loss of some GASA genes in the D-subgenomes during evolution may be the reason why the GASA genes of tetraploid cotton are not the sum of the GASA genes of two diploid cotton.

GhGASA genes are possibly involved in plant hormone responses

GA is an important phytohormone that influences multiple developmental and physiological processes (Sun and Gubler 2004; Li et al. 2010), as well as stress responses (Alonso-Ramírez et al. 2009; Qin et al. 2011; Magome et al. 2008; Zhu et al. 2012). As a class of low-molecular-weight signal peptides, GASA proteins play an important role in the GA signaling pathways (Silverstein et al. 2010), thereby affecting plant growth. The functions of many of GASA proteins have been illustrated in plants. In this study, the cis-acting elements of GhGASA genes were classified and analyzed to explore their potential functions (Additional file 10: Fig. S8, Additional file 11: Fig. S9). The results of the cis-acting element analysis indicated that GhGASA genes might be related to many plant hormone responses (abscisic acid, GA, salicylic acid, MeJA, and auxin), as well as stress responses. Meanwhile, we constructed a protein–protein interaction network to predict the functions of the GhGASA proteins (Fig. 6). The results showed that GhGASA13 and GhGASA35 proteins could interact with the TPR repeat-containing thioredoxin TTL1 that participates in osmotic stress and abscisic acid (ABA) responses and be a positive regulator of ABA signals during germination and seedling development under stresses (Rosado et al. 2006; Lakhssassi et al. 2012). Previous studies have demonstrated that GASA proteins can respond to plant hormone signaling. In Arabidopsis, AtGASA2/3/5 and AtGASA14 respond to ABA signaling (Roxrud et al. 2007). In rice, OsGSR1 can influence brassinosteroid signaling (Wang et al. 2010). Moreover, OsGASA1/7 proteins have been suggested to play a role in BR and GA signaling (Lee et al. 2015; Boonpa et al. 2018). Our analysis suggest that GhGASA proteins may be involved in phytohormone responses.

Expression of GhGASA genes in different tissues provide the reference for further research

The expression profiles of genes in different tissues of plants can provide important clues for the identification of their functional properties. To reveal these expression patterns more intuitively, we drew an summary chart to present the expression locations of GhGASA genes in cotton based on the transcriptome data (Fig. 7). As shown, GhGASA23 and GhGASA24 showed high expression levels during the development of the fiber stages. In practice, GhGASL3, the orthologous gene of GhGASA24, and GhGASL5, the orthologous gene of GhGASA23, may be involved in early stage of fiber elongation, according to previous findings (Liu et al. 2013). Actually, qRT-PCR analysis of GhGASA23 and GhGASA24 again hinted that GhGASA23 and GhGASA24 play a positive role in the early elongation of cotton fiber cell. These results indicate that GhGASA23 and GhGASA24 genes may be involved in cotton fiber development, provide the reference for the in-depth analysis of cotton fiber development related genes. Analysis of transcriptome data and qRT-PCR also found other fiber-related genes, including GhGASA3, GhGASA8, GhGASA14, and GhGASA19. Previous studies confirmed that GASA proteins are important for regulating floral development in other species (Kotilainen et al. 1999; Peng et al. 2010). However, little is known about their potential roles in floral development in cotton. Here, we investigated the GhGASA genes expression patterns in floral organs and found that GhGASA9, GhGASA31, and GhGASA18 are highly expressed in most floral organs, including the anther, filament, pistil, bract, sepal, petal, and torus, with the highest expression levels observed in the petal. These suggests GhGASA genes play an important role in the floral development in cotton. Meanwhile, we found that GhGASA38 was only expressed in the root, showing that it may have special functions in root development. GhGASA4, GhGASA8 and GhGASA26 in stem, GhGASA14 in leaf with high expression level by qRT-PCR analysis, suggesting theroles in stem or leaf development.

Fig. 7

Summary chart for cotton and the expression locations of the GhGASA genes. Different colors indicate different tissues. The degree of GhGASA expression is represented by a graduated color

GhGASA genes are involved in abiotic stress response

Research has confirmed that GASA genes are closely related to abiotic-stress responses. For example, the overexpression of GASA4 in Arabidopsis enhanced resistance to heat stress (Ko et al. 2007), abiotic stress affects the expression of the GASA5 gene (Zhang and Wang 2011), and the different expression levels of OsGASA1/8/10 show a response to salt stress (Muhammad et al. 2019). In our study, GhGASA9 and GhGASA31 presented high expression levels under abiotic stress, especially under heat treatment. GhGASA4 and GhGASA26 expression was also increased under salt stress. The expression of GhGASA25 was significantly downregulated under cold stress but upregulated under salt treatment. These results show that GhGASAs (GhGASA4, GhGASA9, GhGASA25, GhGASA26 and GhGASA31) may be involved in stress response. Meanwhile, qRT-PCR analysis of GhGASA4, GhGASA9, GhGASA18, and GhGASA25 confirmed GhGASA genes involved in cold stress response. Especially, GhGASA4 and GhGASA18 were induced intensely in reply to cold environment. Xinjiang province is the primary cotton-growing areas in China. Cold damage on cotton happens widely due to the late spring coldness of Xinjiang, resulting in yield decline. So, it is important to improve the cold-resistance of cotton. This study provides reference for further research on cotton's defense against low temperature.

Conclusion

Cotton is the most important crop for the production of renewable fiber. In this paper, a systematic whole-genome analysis of the GASA gene family was performed for the first time in five cotton species. As a result, 132 GASA genes were identified and subjected to protein structural, chromosomal localization, and phylogenetic characterizations. Further focusing on the gene structures, synteny analyses, tissue-specific expression patterns, expression profiles under abiotic stress, and protein interactions of GhGASA genes enriched our understanding of the GASA gene family of cotton. Finally, we found several GhGASA genes that may be related to fiber development and stress response. The results also suggest that the GhGASA proteins are closely related to phytohormone signaling. These findings provide valuable information for further study into the functions of these candidate GASA proteins.