Genome-wide analysis of growth-regulating factors (GRFs) in Triticum aestivum

Identification of GRF genes in T. aestivum, T. urartu, T. dicoccoides and Ae. tauschii

Genome-wide data for Triticum aestivum (IWGSC v1.1), Triticum urartu (v1.43), Triticum dicoccum (v1.0.43), and Aegilops tauschii (v4.0.43) were downloaded from Ensembl Plants database (http://plants.ensembl.org/index.html). First, the Hidden Markov Model (HMM) of WRC (PF08879) and QLQ (PF08880) domains were obtained from PFAM (http://pfam.xfam.org/) and used as query sequences for HMMER3.0 (http://hmmer.org/download.html) searching (e-value ≤1e⁻¹⁰). Second, download of known GRF protein sequence were used as query sequences, including 9 GRFs from Arabidopsis (Berardini et al., 2015), 14 GRFs from Zea mays (Andorf et al., 2015), and 12 GRFs from Oryza sativa (Ouyang et al., 2006). They were then used as query sequences for BLASTp searching the wheat database. The first uncurated protein sequence list was genereted by e-values lower than 1 × 10⁻¹⁰. Next, we combined the results and deleted the redundant sequences. Finally, predicted proteins were considered as GRFs only if they contained QLQ and WRC conserved domains verified by NCBI CDD (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) and SMART (http://smart.embl-heidelberg.de/) (Letunic & Bork, 2017).

Characterization of TaGRFs proteins

ExPASy server10 (https://web.expasy.org/compute_pi/) was used to predict the amino acid length, molecular weight (MW), isoelectric point (pI), stability, and grand average of hydropathicity (GRAVY) for TaGRFs proteins (Li et al., 2018); and subcellular localization prediction was carried out by Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/) (Chou & Shen, 2010).

Chromosomal location and gene duplication of TaGRFs

The wheat genome GFF3 gene annotation file was from the wheat database IWGSC v1.1 (https://wheat-urgi.versailles.inra.fr/Seq-Repository/Assemblies). Gene structure annotations of TaGRFs were extracted from the GFF3 file. The start and end location information of the TaGRFs in the corresponding chromosomes was used to draw the physical map by the software MapInspect (Fang et al., 2019). The orthologous genes from wheat and its subgenome donor were identified by the common tool “all against all BLAST search”. The cutoff values (e-value <10⁻¹⁰, identity >80%) were used to assure the reliability of the orthologues. Then we used Multiple Collinearity Scan toolkit (MCScanX) to depict their homology relationships (Wang et al., 2012). Synteny diagrams were generated using the R package “circlize”. Gene duplication events were divided into tandem duplication events and segmental duplication events. Tandem duplication events were determined by the following evaluation criteria: (1) length of the aligned sequence >80% but of each sequence, (2) identity >80%, (3) threshold ≤10⁻¹⁰, (4) only one duplication can be recognised when genes are tightly linked; and (5) intergenic distance is less than 25 kb. When genes passed the criteria for (1), (2), and (3), but were on a different chromosome, they were deemed to be segmental duplications (Fang et al., 2020; Jiang et al., 2020).

Analysis of TaGRFs gene structures and motifs

According to the TaGRFs annotation information, GSDS2.0 (http://gsds.gao-lab.org/Gsds_about.php) was used to produce TaGRFs genetic structure (Hu et al., 2017). The MEME v4.9.1 (http://meme-suite.org/) was used to identify conserved TaGRF protein motifs (Zheng et al., 2017). The trained parameters were applied as follows: each sequence may contain any number of nonoverlapping occurrences of each motif, up to 20 different motifs, and a motif width range of 6 to 50 amino acids (aa). These motif patterns were drawn using TBtools software (Chen et al., 2020). The annotations of those predicted motifs were analyzed by SMART (http://coot.embl-heidelberg.de/SMART/) (Letunic & Bork, 2017). Multiple amino acid sequences were aligned using DNAMAN6.0 (Lynnon Biosoft).

Phylogenetic analyses of TaGRFs

The 99 protein sequences (9 AtGRFs, 12 OsGRFs, 14 ZmGRFs, 30 TaGRFs, 6 TuGRFs, 10 AeGRFs, and 18 TdGRFs) were conducted multiple comparisons by using ClustalW2 software (Thompson, Higgins & Gibson, 1994). Then, the phylogenetic relationships were inferred using the Neighbor-Joining (NJ) method with bootstrap analysis for 1,000 repetitions by MEGA7.0 (Kumar, Stecher & Tamura, 2016). Finally, the midpoint rooted base tree was drawn using Interactive Tree of Life (IToL, v4, http://itol.embl.de) (Letunic & Bork, 2019).

Cis-acting elements analysis of TaGRFs

The PlantCARE (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used to predict cis-acting elements in the regions 1,500 bp upstream of 30 TaGRFs start codons (Lescot et al., 2002). The predicted results were organized and displayed by the R package “pheatmap” (Jiang et al., 2019; Zhu et al., 2019b).

Gene Ontology annotation in TaGRF family genes

The functional annotation of GRF sequences and the analysis of annotation data were performed using Blast2GO (http://www.blast2go.com) (Conesa et al., 2005). The full-length amino acid sequence of the TaGRF proteins were uploaded to the original program, drawn and annotated. The program provides the output defining three categories of GO classification namely biological processes, cellular components, and molecular functions.

Multiple conditional transcriptome analysis of TaGRFs

The multiple transcriptome data were downloaded from the Wheat Expression Browser (http://www.wheat-expression.com/) (Ramírez-González et al., 2018); and the heat maps of TaGRFs were generated using the R package “pheatmap”.

Growth and stress treatment of wheat seedlings

Seeds of Emai 170 (a hexaploid common wheat cultivars) were sterilized on the surface with 1% hydrogen peroxide, rinsed thoroughly with distilled water, and germinated in an incubator at 25 °C for 2 days (He et al., 2020; Ma et al., 2016). According to the reported method, the seedlings were transferred to 1/2 strength Hoagland nutrient solution and cultured in continuous ventilation (Yin et al., 2019; Zhu, Gong & Yin, 2019). After five days (when wheat seeding reached the stage of one heart and one leaf), 85.5 mM NaCl and 82.5 mM mannitol were applied to seedings. Every two days, 2 M KOH or 0.4 M H₂SO₄ was used to adjust the pH of culture solution to 6.0. During the application, the plants were grown at 16 h/8 h (day/night) and 25 °C. Leaves and roots were collected at 2 h, 4 h, 8 h, 12 h, 24 h, 96 h and 144 h after treatments. Three biological repeats are included for each treatment. Finally, the samples were immediately frozen with liquid nitrogen and stored at −80 °C.

RNA isolation and qRT-PCR analysis

According to the manufacturer’s instructions, total RNA of samples were extracted by TRizol reagent (Invitrogen, USA) and cleansing DNA with DNaseI (TaKaRa, USA). The first cDNA was reverse-transcribed from RNA by RevertAid Reverse Transcriptase (Vazyme, China). Gene-specific primers were designed using Primer 5.0; and the ADP-ribosylation factor Ta2291 (F: GCTCTCCAACAACATTGCCAAC, R: GCTTCTGCCTGTCACATACGC) was used as an internal reference gene for qRT-PCR analysis (Paolacci et al., 2009). The qRT-PCR reaction system and protocol were carried out as manufacturer’s instructions for SYBR® (Vazyme, China). For each sample, settings included three technical replicates. Relative gene expression level was calculated using the 2^−ΔΔCt method (Yin et al., 2018b).

Identification and analysis of wheat GRF transcription factor gene family members

For identification of GRF TF genes in wheat, both BLAST and Hidden Markov Model (HMM) searches were performed. The 35 known GRF proteins (Table S1), including Arabidopsis (9), maize (14), rice (12) as the query sequences to conduct BLASTp against the wheat reference genome IWGSCv1.1. Using the HMM of WRC (PF08879) and QLQ (PF08880) domains were used as the query sequences for HMMER3.0 searching. The candidate proteins were verified by NCBI CDD and SMART Online Tools to determine that the TaGRF contained both WRC and QLQ domains. Finally, a total of 30 TaGRFs were identified from the wheat genome. We named wheat GRF genes (TaGRFs) according to the naming rule of Schilling et al. (2020); the corresponding gene IDs are shown in Table 1. Using the same method, we identified 6, 10, and 18 GRFs from T. urartu, Ae. tauschii and T. dicoccoides, respectively (Table S2). The deduced polypeptides ranged in length from 206 (TaGRF2-2A) to 611 (TaGRF4-4B) amino acids, with the predicted molecular weights ranging between 21.6 to 64.2 kDa. Their isoelectric points ranged from 4.72 (TaGRF1-2B) to 10.23 (TaGRF2-2A). Their instability parameters were between 41.46 (TaGRF7-6A) to 66.25 (TaGRF12-7B). Their average hydrophilicity coefficient ranged from 0.269 (TaGRF2-2A) to 0.882 (TaGRF12-7B) (Table 1). Subcellular localization predictions showed that all GRF proteins except TaGRF1-2A and TaGRF2-2A were localized only in the nucleus, while TaGRF1-2A was located in the chloroplast, cytoplasm and nucleus, and TaGRF2-2A was located in the chloroplast and nucleus.

According to the phylogenetic relationships (Fig. 1), the 30 TaGRFs could be divided into four sub-categories (Group I to IV). Group I consisted of a single member, TaGRF7-6A. Group II included TaGRF1-2A, TaGRF2-2A, TaGRF1-2B, TaGRF2-2B, TaGRF1-2D, TaGRF2-2D, TaGRF8-6A, TaGRF8-6B, and TaGRF8-6D. Group III consisted of TaGRF4-4A, TaGRF5-4A, TaGRF4-4B, TaGRF5-4D, and TaGRF4-4D. The remaining TaGRF genes were classified in Group IV.

The TaGRF gene structure map showed that all wheat GRF gene members contain 1 to 4 introns, with the majority having 2 to 3 introns (Fig. 1, and Table S3). There were 2 to 5 exons, with most TaGRF having 2 to 4 exons. The exon number of TaGRF genes within same group were relatively consistent.

Conservative motif analysis indicated that TaGRF protein domains are highly conserved among the 30 members. Each member contains only two structural domains: WRC (Motif 1) and QLQ (Motif 2) (Fig. 1). Lengths and the most matching sequences of 20 motifs were shown in Table S4.

In order to further analyze the conservation degree of QLQ and WRC domains in TaGRFs, we performed multiple sequence alignment of these two domains. The results indicate that, as highlighted in Fig. 1, the QLQ and WRC motifs are highly conserved. The N-terminal QLQ motif was conserved with one Leu and two Gln in all the TaGRF proteins. The WRC motif was also highly conserved with one Trp, Arg, and Cys in each of the TaGRF proteins. A zinc finger motif (CCCH) was also found within the WRC domain in all TaGRF proteins (Fig. 2).

Chromosome localization of wheat TaGRF genes

Based on the GFF3 genome reference files, the chromosome map of TaGRF genes was generated using MapInspect software (Fig. 3, and Table S3). The three sub-genomes A, B, and D contained 15, 11, and 10 TaGRFs, respectively. But the TaGRFs are not uniformly distributed among chromosomes (chromosome 2, 9; chromosome 4, 6; chromosome 6, 10; and chromosome 7, 5). However, distribution range of genes in different group was diverse. Members of TAGRF group II are only distributed on chromosome 4.

Phylogenetic analysis of GRF transcription factor family members in wheat, rice, maize, Arabidopsis, T. urartu, Ae. tauschii and T. dicoccoides

The phylogenetic analysis of wheat (30), rice (12), maize (14), Arabidopsis (9), T. urartu (6), Ae. tauschii (10) and T. dicoccoides (18) GRFs showed that 99 GRFs could be divided into 4 sub-categories (Group I to IV) (Fig. 4). Group I contained only AtGRF7 and AtGRF8 of Arabidopsis, TaGRF7-6A of wheat. The Group II included 9 TaGRFs, 3 AtGRFs, 2 ZmGRFs, 3 OsGRFs, 3 AeGRFS, 6 TdGRFs, and 3 TuGRFs. The group III consisted of 5 TaGRFs, 2 AtGRFs, 3 ZmGRFs, 4 OsGRFs, 2 AeGRFs, and 4 TdGRFs. The Group IV included, 15 TaGRFs, 2 AtGRFs, 5 OsGRFs, 9 ZmGRFs, 5 AeGRFs, 8 TdGRFs, and 3 TuGRFs. In addition, we found that in Group II - IV, TaGRF has a closer phylogenetic relationship with TuGRF, AeGRF and TdGRF, followed by OsGRF and ZmGRF. and relatively distant from the AtGRF.

TaGRF gene promoter Cis- element analysis

Analysis of the cis-elements in the promoter sequence was important for understanding the regulatory functions of genes. The cis-acting element analysis was performed in Plant-CARE by using upstream sequences (1.5 kb) of TaGRF genes extracted from the wheat genome. The detailed information including function and location were displayed in Table S5. The results showed that all 30 TaGRF genes contained several TATA boxes and CAAT boxes, indicating that TaGRF genes can be normally transcribed. When focusing on the cis-acting elements associated with wheat growth and development, hormonal and stress responses, it can be seen in Fig. 5 the TaGRF gene promoter contains a large number of cis-elements, with the largest number found in TaGRF9-6D having 19 cis-elements, and TaGRF3-2A with the least, containing only 8 cis-elements. There were several different light-related elements in these cis-elements, such as AE-box, Box 4, I-box, C-box, Sp1, circadian, CAG-motif, 3-AF1 binding site, LAMP-element, TCT-motif, GATT-motif, ATCT-motif, and Gap-box. This suggests that the GRF gene family may play a role in light response.

In addition, a large number of responsive hormones and stress-related cis-elements were found in the promoter region of the TaGRFs, including auxin (11 TGA-elements), gibberellin (8 GARE-motifs and 6 P-boxs), jasmonic acid methyl ester (65 TGACG-motifs), abscisic acid (72 ABREs) and other hormone response components, as well as anaerobic induction (25 AREs), drought (19 MBSs) and low temperature (3 LTRs) and other stress response cis-elements. This suggests a potential role for the wheat GRF family in wheat growth and development and in a variety of hormones and stress.

Gene Ontology annotation in TaGRF family genes

The GO item analysis was performed using Blast2Go and the results indicated the putative participation of 30 TaGRF proteins in diverse biological processes (Fig. 6, and Table S6). Total ten different GO items of biological processes were defined. Majority of the TaGRFs were predicted to function in ‘regulation of transcription, DNA-templated (GO: 0006355)’ (76.67%), followed by ‘response to deep water (GO: 0030912)’ (20%) and ‘response to gibberellin (GO: 0009739)’ (20%). Molecular function prediction showed that about 76.67% of the TaGRFs were evidenced to participation of ‘ATP binding (GO: 0005524)’ and ‘hydrolase activity, acting on acid anhydrides, in phosphorus-containing anhydrides (GO: 0016818)’. Cellular localization prediction indicated that the majority of TaGRF proteins (80%) were localized in the nucleus (Fig. 6).

Homologous gene pairs and synteny analysis

Gramineae evolved 50–70 million years ago, and the Pooideae subfamily, which includes barley and wheat, evolved about 20 million years ago (Inda et al., 2008; Peng, Sun & Nevo, 2011). Obviously, common wheat has a complicated evolutionary history, and its ancestors’ origin is affected by many factors, but research shows that wheat has two major polyploid evolutionary events (Ling et al., 2013). Homology reflects the phylogeny of a species, so it can be used to transfer annotations for one known gene to another newly sequenced genome. In order to further infer the evolutionary origin and homology of the wheat GRF family, Sixty-four GRFs were identified from T. aestivum (30 TaGRFs), T. urartu (6 TuGRFs), T. dicoccoides (18 TdGRFs) and Ae. tauschii (10 AeGRFs) using a computer-based method (Fig. 7A, and Table S7). There were no paralogous gene pairs in Ae. tauschii and T. urartu, and 21 and 5 paralogous gene pairs in T. aestivum and T. dicoccoides, respectively. Among them, 18 orthologous gene pairs were identified from T. aestivum and Ae. tauschii, 5 orthologous gene pairs were identified from T. dicoccoides and Ae. tauschii, 5 orthologous gene pairs were identified from T. urartu and Ae. tauschii, 20 orthologous gene pairs were identified from T. aestivum and T. dicoccoides, 12 orthologous gene pairs were identified from T. aestivum and T. urartu, and 3 orthologous gene pairs were identified from T. urartu and T. dicoccoides. Given phylogenetic analyses and homology results of four wheat species, it was speculated that 8 TaGRFs (TaGRF2-2B, TaGRF3-2B, TaGRF4-4A, TaGRF4-4B, TaGRF6-4A, TaGRF8-6B, TaGRF9-6A and TaGRF9-6B) originated from T. dicoccoides, 9 TaGRFs (TaGRF2-2D, TaGRF3-2D, TaGRF4-4D, TaGRF5-4D, TaGRF8-6D, TaGRF9-6D, TaGRF10-6D, TaGRF11-7D and TaGRF12-7D) from Ae. tauschii, and 5 TaGRFs (TaGRF3-2A, TaGRF8-6A, TaGRF8-6B, TaGRF10-6A and TaGRF12-7A) from T. urartu.

Gene expression pattern analyses of TaGRFs

For multigene families, analysis of gene expression patterns often provides useful clues for determining gene function. Transcriptome data from growth and abiotic stresses were downloaded from Wheat Expression Browser to examine their expression patterns. The results showed that 28 TaGRF genes (except TaGRF2-2A and TaGRF7-6A) were expressed in different tissues or under different stress treatments (Fig. 8, and Table S8). TaGRF1-2A and TaGRF1-2D were highly expressed in various tissues. TaGRF5-4A and TaGRF5-4D had the highest expression in shoot tip meristem. Most of the genes were expressed in shoot tip meristems more significantly than other tissues. About half of the TaGRF genes were expressed under NaCl treatment, and TaGRF4-4A, TaGRF4-4B and TaGRF4-4D were significantly expressed under NaCl treatment. It was speculated that TaGRF family members may play important roles in the development of wheat shoot tip meristems, and TaGRF4 may play an important role in wheat salt tolerance.

Quantitative-real time PCR analysis

To further understand the potential role of the TaGRF genes in abiotic stresses (NaCl and mannitol), qRT-PCR was used to analyze the expression pattern of TaGRFs. Based on transcriptome analysis, we selected 14 TaGRFs for qRT-PCR. In the two treatments of this study, the expression of all 14 TaGRFs differed from the control, although their degree of difference was often substantial (Fig. 9).

After treatment with NaCl, the expression of TaGRF4-4A, TaGRF4-4D, TaGRF5-4A, TaGRF10-6A, TaGRF10-6B and TaGRF10-6D were higher in treated leaves than in the control 96 h after treatment. Expression of the genes TaGRF1-2B, TaGRF1-2D, TaGRF3-2D, and TaGRF9-6D were lower in the leaves 96 h after treatment than in the control group. However, in the roots, the expression of TaGRF1-2B and TaGRF3-2D were much higher than in the control 96 h after treatment, while the expression of TaGRF4-4A, TaGRF4-4D, TaGRF10-6B, and TaGRF10-6D were much lower than in the control 96 h after treatment. Among these, the expression trends of TaGRF1-2B, TaGRF3-2D, TaGRF4-4A, TaGRF4-4D, TaGRF10-6B, and TaGRF10-6D were completely reversed in the roots compared to the leaves 96 h after treatment.

After treatment with mannitol, TaGRF1-2B, TaGRF6-4A, TaGRF10-6A and TaGRF9-6D were expressed higher in the roots than in the control group 24 h and 96 h after treatment. However, in the leaves, only TaGRF6-4A, TaGRF10-6A and TaGRF9-6D was expressed higher than the control group at 24 h and 96 h after treatment. The expression level of TaGRF1-2B did not change compared with the control group. TaGRF1-2B and TaGRF3-2B were significantly down-regulated in leaves at 2 and 4 h after treatment.