Genome-wide identification and expression profile of GhGRF gene family in Gossypium hirsutum L.

Identification of GRF genes

The genomic and protein sequences of G. hirsutum (NAU, v1.1; ZJU, v 1.0; JGI, v1.1) were downloaded from CottonFGD (https://cottonfgd.org/about/download.html) (Zhu et al., 2017), NCBI (https://www.ncbi.nlm.Nih.gov/genome/), and Phytozome 13 (https://phytozome.jgi.doe.gov/pz/portal.html). GRF protein sequences of rice were downloaded from NCBI GenBank, and the protein sequences of AtGRFs were acquired from TAIR 10 (http://www.arabidopsis.org). AtGRF sequences were used as query sequences to blast against the G. hirsutum protein database to search homologous candidate sequences in the GRF gene family with default parameters. The candidate sequences were submitted to the Batch CD-search tool in NCBI and ScanProsite tool in the ExPASy prosite database (https://prosite.expasy.org/ scanprosite/) to detect QLQ and WRC domains, which are the structural features to distinguish GRF family members from other gene families. The GRF family genes were determined by predicting the conserved domain of the candidate sequence.

Gene structure and protein motifs of GhGRFs

We used the online software Gene Structure Display Server (GSDS) for gene structure analysis (http://gsds.cbi.pku.edu.cn/) (Hu et al., 2015), as previously described (Wang et al., 2021). The CDS sequences were aligned with the reverse complementary sequence of miRNA396a, which was acquired from previously published research papers (Hewezi et al., 2012; Yang et al., 2009). Motif analysis of the full-length protein sequence was performed using the MEME online software (Bailey et al., 2015) (https://meme-suite.org/meme/) as previously described (Qanmber et al., 2019b), according to the following parameters: classic mode, zero or one occurrence per sequence, and motif number 10.

Biophysical properties and sequence logos of GhGRF genes

Various physical and chemical parameters, such as the number of amino acids, MW (molecular weight), pI (isoelectric point), and GRAVY (grand average of hydropathicity) of GhGRF gene coding protein sequences were determined using the ExPASy ProtParam tool (http://us.expasy.org/tools/protpara-m.html). The GRF protein sequences were aligned using Muscle in MEGA 5.2 (Tamura et al., 2011), and the results were analyzed with the WEBLOG online program (http://weblogo.plusone.com/create.cgi) (Crooks et al., 2004) as previously described (Ali et al., 2020; Qanmber et al., 2019a) to generate the sequence logos of conserved amino acid residues.

Classification and phylogenetic analyses of GRF genes

Classification and phylogenetic analyses were performed using MEGA 5.2. A neighbor-joining tree of G. hirsutum, O. sativa, and Arabidopsis was constructed using full-length protein sequences according to the following parameters: Poisson correction, 95% partial deletion, and bootstrap analysis with 1000 replicates, as previously described (Xiao et al., 2018; Yu et al., 2018). Another phylogenetic tree was built using maximum likelihood (ML) methods (default parameters and JTT+G as amino acid substitution models) to confirm the authentication of the phylogenetic tree.

Promoter cis-elements and expression pattern analyses of GhGRF genes

For cis-element analysis, 2.5-kb upstream promoter regions were downloaded from CottonFGD and submitted to the plantCARE database (https://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002), as previously described (Qanmber et al., 2019a; Wu et al., 2021). The FPKM value of GhGRF gene expression was downloaded from CottonFGD to draw heatmaps and demonstrate the expression profiles of GhGRF genes when developing ovules and fibers, mature organs (calycle, leaf, petal, pistil, root, stamen, stem, and torus), germinating seeds, cotyledons, roots, and seedlings treated with abiotic stresses and hormonal treatments.

Plant material and Real-time quantitative PCR

For expression pattern analysis, tissue samples were collected under field conditions. Cotton plants growing under normal conditions were placed into a container with a predetermined concentration of BL, IAA, JA, PEG, and NaCl. The samples were collected at different time points (0, 0.5, 1, 3, and 5 h), immediately frozen into liquid nitrogen, and stored at −80 °C for subsequent RNA extraction and qRT-PCR analysis. We used the Trizol (Invitrogen, Carlsbad, CA, USA) method to extract total RNA from each cotton sample, as previously described. cDNA was synthesized using a PrimeScript 1st strand cDNA Synthesis Kit (Takara), according to the manufacturer’s instructions. RT-qPCR was performed using SYBR Green Master Mix (Takara) on an IQ5 Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) (Wang et al., 2014b). RT qPCR experiments were performed using the actin gene as an internal reference gene. The PCR program consisted of an initial polymerase activation at 95 °C for 1 min, followed by 40 cycles of 95 °C for 10 s and 60 °C for 30 s. Upon completion of the PCR, a dissociation curve was generated to confirm the specificity of the product and avoid generating primer dimers.

Identification and naming of GhGRF genes

To identify the GRF gene family, candidate sequences lacking both domains (QLQ and WRC) or with only one of the two domains were discarded. Finally, 35 GRF genes were identified (Fig. 1). To further enhance the credibility of the results, two additional blasts were performed against NAU and ZJU protein sequence data of G. hirsutum. The 42 or 41 candidates were also submitted to the Batch CD-search tool in NCBI, and 35 GhGRF family members were confirmed (Figs. S1 and S2). Multiple sequence alignment was performed using JGI (35), N AU (33), and ZJU (35) GhGRF deduced protein sequences. The gene names and their corresponding IDs are displayed in Table 1, and the nomenclature is described in detail in later paragraphs.

Conserved amino acid residues within QLQ and WRC domains

The GRF gene family members were characterized by the presence of both conserved QLQ and WRC domains in the N-terminal end (Kim, Choi & Kende, 2003; Van der Knaap, Kim & Kende, 2000). To further explore the conservative properties of the QLQ and WRC domains of GhGRFs, multiple sequence alignment was performed to depict sequence logos of QLQ and WRC domains for G. hirsutum, Arabidopsis, and rice. The conserved domain sequence logos showed that most amino acid residues distributed in two domains were highly conserved among three plant species (Fig. 2). Amino acid residues include Q (6), E (9), Q (13), Y (19), V (26), P (27), and L (30) in the QLQ domain, while EP (3-4), RC (6-7), RRTDGKKWRC (9-17), KYC (26-28), H (31) and R (38) in the WRC domain were highly conserved.

Phylogenetic analysis and nomenclature of GRF genes

To better understand the phylogenetic relationships between GRFs in G. hirsutum, rice, and Arabidopsis, a neighbor-joining phylogenetic tree was generated using 56 GRF proteins. GRF genes from the three plant species were classified into five groups (G1 to G5) (Fig. 3). To validate the phylogenetic tree constructed using the NJ method, the maximum likelihood (ML) method was used to build another phylogenetic tree. Results indicated that 56 GRFs were also divided into five clades based on the ML tree. The topology of the ML tree was slightly different from the NJ tree, but the composition in each clade was the same, indicating that the NJ tree could be used for further analysis.

The number of GRF genes differed between the clades of the NJ tree. There were 18, 17, six, five, and 10 members in G1, G2, G3, G4, and G5 clades, respectively. The G1, G2, and G5 groups had the highest coverage of all three species, including the monocot (rice) and dicots (G. hirsutum and Arabidopsis). In contrast, the G3 and G4 groups contained GRF genes only from dicots without no GRF member from the monocot. The proportion of GRF family members from three different plants in each clade differed. In the G1 group, the number of AtGRF, OsGRF, and GhGRF members were two, four, and 12, respectively, while in the G2 group, the number of AtGRFs, OsGRFs, and GhGRFs were two, five, and ten, respectively. Similarly, in the G3 group, there were two AtGRFs and four GhGRFs members; the G4 group has two AtGRFs and three GhGRFs genes. G5 group contained one AtGRF, three OsGRFs, and six GhGRFs. These results indicate that GRF families of different plants experienced different evolutionary processes, and that cotton and Arabidopsis thaliana are dicotyledons and have a relatively close genetic relationship.

In the phylogenetic tree, names were assigned to each GhGRF gene based on their phylogenetic relationship and homologous sequence consistency (Fig. 3). Gh is used as a prefix before the GRF name of upland cotton, At for Arabidopsis GRFs. We proceeded according to the following rules of nomenclature: first, GhGRF genes were named after their orthologs in Arabidopsis, which were part of the same clade of the phylogenetic tree; second, if there were more than one orthologous counterpart of Arabidopsis, the GhGRF genes were named after the Arabidopsis homologous gene possessing the highest sequence consistency with the cotton genes; third, “a”, “b”, and “c” were appended to gene names as suffixes to distinguish the GRF genes sharing the common Arabidopsis orthologs based on the sequence consistency order from high to low. Finally, At or Dt was attached respectively to the end of each gene name to distinguish the sub-genomes in tetraploid cotton in which the gene was located; “t” means tetraploid.

Biophysical properties and chromosomal distribution of GhGRFs

Physical and biochemical properties include the distribution of genes on the chromosome, gene ID, length of coding sequence (CDS), number of amino acids (protein length), molecular weight (MW), isoelectric point (pI), and grand average of hydropathicity (GRAVY) of GhGRF genes and their deduced protein is displayed in Table 1. The CDS length of GhGRF genes ranged from 498 bp (GhGRF9b_At) to 1,833 bp (GhGRF1c_At), and the number of amino acids of their corresponding deduced proteins ranged from 165/17.68 to 610/65.38. Moreover, pI values ranged from 5.89 (GhGRF8_At) to 9.65 (GhGRF9c_At). In addition, the lowest (−0.957) and highest (−0.373) GRAVY values were computed for the GhGRF5a_Dt and GhGRF1c_Dt, respectively. We observed that GRAVY scores for all GhGRF genes were negative, indicating that all GhGRF proteins were hydrophilic; however, the degree of hydrophilicity varied greatly.

Annotation of the GhGRF genome sequences (GFF3 profile) helped us to predict the locations of GhGRF genes on chromosomes (Fig. 4). The 35 GRF genes were unevenly distributed among 20 chromosomes, with 17 genes on nine chromosomes of the At-subgenome and 18 genes on 11 chromosomes of the Dt-subgenome. There were three genes on chromosomes At02, At05, At12, Dt05, and Dt12. For At10, At13, Dt03, Dt10, and Dt13, there were two genes on each chromosome. Ten chromosomes (At06, At08, At09, At11, Dt02, D04, Dt06, Dt108, Dt09, and Dt11) have one GhGRF gene. However, six chromosomes (At01, At03-04, At07, Dt01, and Dt07) have no GhGRF gene. Four pairs of homologous chromosomes, including At05/Dt05, At06/Dt06, At08/Dt08, and At13/Dt13, have a similar distribution of GhGRF genes. However, the gene distribution on three chromosome pairs (At02/Dt02 to At04/Dt04) differed. We inferred that gene loss occurred from the At04 homologous chromosome (GhGRF3b_Dt) to the Dt04 chromosome and unidirectional translocation of a chromosome segment containing the two genes, GhGRF6c-At and GhGRF1b_At, from At03 to At02 occurred (Fig. 4). These results are consistent with those of previous studies (Yang et al., 2017).

Analysis of gene structure and deduced protein motif

To explore the exon-intron structure of GhGRF genes, their coding and genomic sequences were aligned and the number and positions of exons and introns were detected (Fig. 5B). Results indicated that the structural patterns of genes were conserved between Arabidopsis and upland cotton, and the majority of GhGRF and AtGRF genes possessed 2-3/3-4 introns/exons, respectively. The MEME web server was used to analyze the types and distributions of motifs in GhGRF deduced proteins, and ten kinds of motifs were detected (Fig. 5C). Motif 1, characterized by the WRC domain, was predicted in 34 out of 35 GhGRFs. Motif 2, characterized by the QLQ domain, was present in all GhGRFs. In the C-terminus, there were many conserved motifs, including motifs 3-5, 7, and 10. Moreover, motifs 6, 8, and 9 were primarily located in the N-terminal end close to WRC and QLQ of G1 members. This conserved exon-intron and motif structure pattern supports the evolutionary relationship of GRF genes in upland cotton and Arabidopsis.

To investigate whether microRNA396 inhibits the regulation of GhGRF genes on plant growth and development, we tried to find fragments in CDS complementary to microRNA396a through multiple sequence alignment using ClustalW in MEGA5.2. Except for GhGRF9b (which lost the second half of the gene), all other gene mRNA sequences contained highly conserved fragments reverse complementary to microRNA396a (Fig. 6). These results suggested that microRNA396a is essential for controlling the function of GhGRF genes when regulating plant growth and development by decreasing their transcripts.

Analysis of cis-regulatory elements

To analyze cis-regulatory elements, 1,500 bp upstream sequence before the transcriptional start codon was retrieved, and along with the 5′-UTR sequence, was submitted to the PlantCARE database. More than 20 types of cis-elements were identified in GhGRF promoter regions and divided into four categories based on their functions: (1) Phytohormone responsive elements (ABA, MeJA, auxin, GA, and SA); (2) environmental or stress-responsive elements (light, defense, low-temperature, anaerobic conditions, and drought); (3) cis-elements related to growth and development (meristem, palisade mesophyll cells, circadian control, cell cycles, and endosperm); (4) cis-elements related to the regulation of bioactive compound metabolism (zein, phytochrome, flavonoid biosynthesis); (Table 2). The distribution of these elements in the promoters regions of GhGRF was not uniform, indicating that GRF genes became functionally differentiated during germline evolution. The diversity of promoter cis-elements demonstrated that the GRF genes play an essential role in the growth and development of cotton.

Expression patterns of GhGRF genes

To analyze the expression patterns of GhGRF genes, publicly available transcriptomic data was retrieved and used to describe the expression levels of GhGRF genes in mature organs, developing tissues, and abiotic stress conditions (Fig. 7). The expression patterns of GhGRF genes in different tissues indicated that most GhGRF genes transcripts developed high amounts of ovules and fibers, except for two genes (GhGRF1f_Dt and GhGRF6a_Dt) whose expression was barely detected during ovule and fiber development (Fig. 7A). More than half of the GhGRF genes were highly expressed during seed germination and the root and cotyledon development processes (Fig. 7B). As shown in Fig. 7C, 27 genes exhibited high expression in mature organs (pistil). The transcript level of 2–5 genes was high for other organs such as calycle, leaf, petal, root, stamen, stem, and torus. Only one gene (GhGRF9b_Dt) indicated a high transcript level in the torus.

We next observed the expression levels of GhGRF genes against abiotic stresses (cold, heat, salt, and PEG) using transcriptomic data. GhGRF1a_At and GhGRF9b_Dt exhibited downregulated expression under the exposure of cold and PEG, while the expression levels of these two genes were upregulated under heat and salt stress. Moreover, the expression of GhGRF1c _At was high at 6 h and 12 h of cold and heat treatment, while GhGRF3a_At was up-regulated only in response to cold. However, GhGRF8_At was down-regulated in response to cold, heat, and salt.

To validate these results, we performed RT-qPCR analysis of 35 GhGRF genes using cotton tissues (root, stem, leaf, flower, pistil, stamen, ovule (−2, 0, 1, 3, 5, 10, 15, 20, 25), fiber (10, 15, 20, 25), abiotic stress (NaCl and PEG), and hormonal treatments (IAA, BL, and GA). Interestingly, tissue-specific expression patterns were similar to the previously published RNA - seq data. Most GhGRF genes (GhGRF1aAt, GhGRF1aDt GhGRF1bAt, GhGRF1bDt, GhGRF1cAt, GhGRF1cDt, GhGRF1dAt, GhGRF1dDt, GhGRF1eAt, GhGRF1eDt, GhGRF1fAt, GhGRF1fDt, GhGRF3aAt, GhGRF3aDt, GhGRF5aAt, GhGRF5aDt, GhGRF5bAt, GhGRF5bDt, GhGRF6aAt, GhGRF6aDt, GhGRF6cAt, GhGRF6cDt, GhGRF8At, GhGRF8Dt, GhGRF9cAt, and GhGRF9cDt) exhibited preferentially high expression during ovule and fiber development, indicating that they are potentially involved in ovule and fiber development (Fig. 8). The tissue expression trend of most genes is consistent when qRTPCR with transcriptome data, including GhGRF1dDt, GhGRF1eAt, GhGRF5aAt, GhGRF5aDt, GhGRF6aDt, GhGRF6bAt, GhGRF7At, and GhGRF7Dt.

Expression pattern analysis for abiotic and hormonal stresses indicated that all GhGRF experienced both up- and down-regulated expression at different time points. The expression level of GhGRF1aAt, GhGRF1aDt, GhGRF1bAt, GhGRF1bDt, GhGRF1dAt, GhGRF1dDt, GhGRF1eAt, GhGRF1eDt, GhGRF1fAt, GhGRF1fDt, GhGRF3aAt, GhGRF3a Dt, GhGRF3bDt, GhGRF5aAt, GhGRF5aDt, GhGRF5bAt, GhGRF6aAt, GhGRF6aDt, GhGRF6bAt, GhGRF6bDt, GhGRF6cAt, GhGRF6cDt, GhGRF8At, GhGRF9aDt, GhGRF9cAt, GhGRF9cDt was significantly higher for all time points under all abiotic stresses (Fig. 9), suggesting that these GhGRF genes could improve plant resistance to abiotic stresses and are potential candidate genes for further study. The transcript level of GhGRF7At decreased at all time points with NaCl and PEG treatment compared to the control. The transcription level of other genes decreased at 0.5, 1, and 3 h after NaCl and PEG treatment, including GhGRF1cAt and GhGRF1cDt. For hormonal treatments, including IAA, BL, and GA, some genes (GhGRF1cAt, GhGRF1cDt, GhGRF5aAt) had relatively high expression levels at all time points. This highlights the important role they play in hormone signaling pathways (Fig. 10). Other GhGRF genes exhibited moderate to low expression when exposed to IAA, BL, and GA treatments.