3.1 Identification of GL3 genes and chromosomal distribution
All the retrieved GL3 gene sequences of G. hirsutum, G. barbadense, G. raimondii and G. arboreum from Cotton FGD. were verified by Pfam and InterPro. The number of verified GL3 genes in Gossypium raimondii (GrGL3s), Gossypium arboreum (GaGL3s), Gossypium hirsutum (GhGL3s), and Gossypium barbadense (GbGL3s) were 3, 3, 6, and 6, respectively. The number of GL3 genes in tetraploid cotton (G. hirsutum and G. barbadense) was double than that in diploid cotton (G. raimondii and G. arboreum). Utilizing publicly available information from four cotton genomes, including gene locations and chromosome lengths, distribution of genes across chromosomes was analyzed (Fig. 1a). The chromosomal distribution pattern of GL3 genes exhibited a significant level of similarity among the four cotton varieties. The GL3s were unevenly distributed on the chromosomes. There were six GbGL3s distributed on four chromosomes, three on the At subgenome and three on the Dt subgenome. One GbGL3 was located on chromosome A09 and two on chromosome A11. There was 1 GbGL3 on chromosome D08 and 2 on chromosome D11. A similar trend was found for GhGL3sGb, the only difference was the presence of GhGL3 on chromosome A08 instead of chromosome A09. There was one GaGL3 located on chromosome A08 and two on chromosome A11. There was 1 GrGL3 on chromosome D08 and 2 on chromosome D11 (Fig. 1b). There was no positive correlation between chromosome length and the number of GL3. This again demonstrated that allotetraploid cotton was formed by intergenomic hybridization of A-genome diploids and D-genome diploids, followed by chromosome doubling (Paterson et al., 2012).
3.2 Phylogenetic analysis of the cotton GL3 gene family
As the sequences of Gb_A09G2187 and Gh_A08G2088 were particularly short compared to the other GL3s, no further analysis was performed for these two genes. To further explore the possible evolutionary relationship of GL3 proteins among four cotton species, an unrooted Maximum Likelihood (ML) phylogenetic tree was constructed using the full-length amino acid sequences of the identified genes. The phylogenetic analysis presented in Figure 2 revealed the division of GL3 genes into three primary clades. Among the cotton species mentioned, Clades I, II, and III exhibited 6, 4, and 6 GL3 genes, respectively. Notably, each clade encompassed members from all species under study, with tetraploid cotton displaying nearly twice as many GL3 genes compared to diploid cotton.Gh_D11G0216 (GhDEL61) in classes Ⅰ and Gh_D08G2104 (GhDEL65) in class Ⅱ were mainly localized in the nucleus.
3.3 GL3s protein sequence analysis
The GL3s have one type of structural domains, bHLH_domain (HLH, named IPR011598 in the InterPro and PF00010 in the Pfam). bHLH domain consists of a basic N-terminal region with at least five amino acids in the basic region and a highly conserved HER motif (His5–Glu9–Arg13). The helix-loop-helix (HLH) region at the C-terminal consists of two amphiphilic alpha-helices, mainly composed of hydrophobic amino acids, which are linked by a loop of variable length. Specific recognition of proteins containing bHLH motifs relies on the formation of homo- or heterodimers with other bHLH proteins (Feller et al., 2011; Shangguan et al., 2016). All GL3s possess bHLH-MYC and R2R3-MYB transcription factors N-terminal (bHLH-MYC_N, named PF14215 in the Pfam), which is the N-terminal region of MYB family and MYC transcription factors (Fig. 3). The physicochemical properties are listed in Table 2. The GbGL3 protein was typically composed of a range of 590 to 631 amino acids. Similarly, the GhGL3 protein exhibited a length between 624 and 631 amino acids. On average, GhGL3 was found slightly longer than GbGL3.GbGb. The average MW of GbGL3s was 62339 Da, and those of GhGL3s, GaGL3s, and GrGL3s were 66402.85 Da, 70027.28 Da and 70027.28 Da, respectively. The PI value for the GL3s of the four cotton species ranged from 5.35 to 5.95. The predicted subcellular location of the GL3 proteins suggested that the GL3 proteins were located almost in the nuclear region. Ga11G3131, Gorai.007G099300 and Gb_D11G0963 were located in Cytoplasm. Gh_D11G0951 was located in Chloroplast.
3.4 Gene structure and motif analysis
To analyze intron-exon structure a comparative analysis of intron-exon structures was performed by aligning genomic DNA with full-length cDNA as depicted in (Fig. 4b). All the four cotton species shared similar pattern of gene structure. GhGL3s, GbGL3s, GrGL3s, and GaGL3s possessed 7 ~ 9 exons and could be divided into three categories based on exon number, relative position, and size. GL3s in group Ⅰ had 7 exons, group Ⅱ and group Ⅲ each had 8 exons, except Ga11G3131 with 9 exons. Notably, the intron-exon structures were consistent with the alignment clusters of GL3s, indicating that the exon-intron structure and the phylogenetic relationship among GL3s were highly correlated. MEME software was used to identify conserved motifs proteins. Twenty conserved motifs (designated motif 1 to motif 20) were identified in GhGL3s, GbGL3s, GrGL3s, and GaGL3s, as shown in s (Fig. 4c). Sequence information for each motif from G. hirsutum, G. barbadense, G. raimondii, and G. arboreum is provided in Table 3. It is clear that members of GL3s in the same group have similar composition pattern. The function of some of the conserved motifs remains unclear. Motif 18 and motif 19 belong to group Ⅰ and group Ⅲ; motif 12 and motif 15 belong to group Ⅰ and group Ⅱ; and motif 13 and motif 16 belong to group Ⅲ.
3.5 Cis-regulatory elements analysis of GL3 genes in allotetraploid cotton
To elucidate the underlying function of GL3s in abiotic stress responses, the promoter region within 2 kb upstream of GL3s was analyzed to identify cis-regulatory elements using data from the PlantCARE database. A total of 27 cis-acting elements wre identified in the GL3 promoter regions, primarily divided into two categories, i.e. environmental stress-related elements and hormone-responsive elements (Fig. 5). Within the environmental stress-related elements, it was observed that a minimum of 8 types of cis-acting elements were identified as light-responsive elements.. Moreover, a number ofy anaerobic induction (ARE) associated elements; low-temperature response (LTR), defence and stress response (TC-rich repeats), and drought induction (MBS) elements were also identified. Some cis-elements predicted to be involved in phytohormone responses were found in GaGL3, GrGL3, GhGL3, and GbGL3s. There were 2 GaGL3, 3 GrGL3, 4 GhGL3, 4 GbGL3 genes with at least 1 salicylic acid-responsive element (TCA-element), 2 GaGL3, 3 GrGL3, 3 GhGL3, 2 GbGL3 genes with at least 1 ABA-responsive element (ABRE), 1 GaGL3, 1 GrGL3, 2 GhGL3, 2 GbGL3 genes with 2 MeJA responses (CGTCA-motif & TGACG-motif), and 1 GaGL3, 2 GrGL3, 2 GhGL3, 3 GbGL3 genes with at least 1 gibberellin response element (P box) (Table 4). In this study, these four hormone response elements were identified in G. barbadense. These cis-regulatory elements were found to be associated with abscisic acid (ABA), salicylic acid, gibberellin, and MeJA, suggesting that GL3 genes potential role in hrormone signal transduction. Phytohormones play a crucial role in plant growth and development. Trichome differentiation in plants is also controlled by phytohormones. InA. thaliana, trichome initiation requires multiple plant hormone signaling, including GAs, cytokinins (CKs), jasmonates (JAs), and salicylic acid (SA) (Traw and Bergelson, 2003; An et al., 2011; Qi et al., 2014). These results suggest that GL3 genes may play an important role in phytohormone signaling in both G. hirsutum and G. barbadense.
3.6 Expression pattern Analysis of GL3s
To gain further insights into the potential role of GL3 genes in allotetraploid cotton, the transcriptomic data from Cotton Omics Database was examined to analyze the expression of cotton GL3 genes in 12 tissues, including root, stem, leaf, ovule_0 DPA, ovule_1 DPA, ovule_3 DPA, ovule_5 DPA, ovule_10 DPA, ovule_20 DPA, fibre_10 DPA, fibre_20 DPA, and fibre_25 DPA. According to the heat map (Fig. 6), GhGL3s and GbGL3s were differentially expressed across different tissues. In G. barbadense, Gb_A11G0219 and Gb_D11G0214 showed high transcript levels in roots and stems, Gb_D08G2182 showed high transcript levels in ovule_3 DPA and 5 DPA, and Gb_A11G0932 and Gb_D11G0963 showed low transcript levels in most tissues. In G.hirsutum, Gh_D08G2104 showed high transcript levels in ovule_1 DPA. The expression of other GhGL3 was relatively low in different tissues.
The Promoter cis-acting elements revealed that GbGbGL3s may respond to MeJA, ABA, SA, and GA regulation. To verify these predictions, the expression changes of Gb_D11G0219, Gb_D11G0214, and Gb_D08G2182 were monitored by RT-qPCR at different time points (0, 3, 6, and 24 h) under MeJA, ABA, SA, and GA treatments. Due to the striking similarity in sequences between Gb_A11G0219 and Gb_D11G0214 alleles, designing distinct RT primers to distinguish two genes became a challenging task. Furthermore, the tissue expression patterns of both genes exhibited significant similarities, leading to the speculation of potential functional redundancy. Consequently, universal primers were employed to identify the expression of these two genes. Gb_A110219 and Gb_D08G2182 showed up regulation under ABA and JA treatments, as depicted in (Fig. 7a, b). Conversely, under GA treatment, the two genes expression levels were down-regulated (Fig. 7c); Moreover SA treatment, expression levels of both genes were down-regulated at 3 h, up-regulated at 6 h, and then down-regulated at 24 h (Fig. 7d). According to the statistics, there were no regulatory elements of SA, ABA, and MeJA in the 2000 bp promoter region upstream of Gb_D11G0219, and no regulatory elements of ABA and MeJA in the 2000 bp promoter region upstream of Gb_D11G0214, but the hormone response results depicted that Gb_D11G0219 and Gb_D11G0214 at least respond to GA, SA, and ABA, and it is possible that the ABA and MeJA response elements are located at the distal upstream regions of the gene. On the other hand, Gb_D08G2182 (Fig. 7e-h), did not show any discernible trend in response to hormone treatments, Therefore, our focus can be directed towards analyzing the expression patterns of Gb_D11G0219 and Gb_D11G0214 in the subsequent studies.