Key wheat GRF genes constraining wheat tillering of mutant dmc

Tillering is a key agronomy trait for wheat (Triticum aestivum L.) production. Previously, we have reported a dwarf-monoculm wheat mutant (dmc) obtained from cultivar Guomai 301 (wild type, WT), and found growth regulating factors (GRFs) playing important roles in regulating wheat tillering. This study is to systematically investigate the roles of all the wheat GRFs (T. aestivum GRFs, TaGRFs) in regulating tillering, and screen out the key regulators. A total of 30 TaGRFs were identified and their physicochemical properties, gene structures, conserved domains, phylogenetic relationships and tissue expression profiles were analyzed. The expression levels of all the TaGRFs were significantly lower in dmc than those in WT at early tillering stage, and the abnormal expressions of TaGRF2-7(A, B, D), TaGRF5-7D, TaGRF10-6(A, B, D) and TaGRF11-2A were major causes constraining the tillering of dmc. The transcriptions of TaGRFs were significantly affected by exogenous indole acetic acid (IAA) and gibberellin acid (GA3) applications, which suggested that TaGRFs as well as IAA, GA signaling were involved in controlling wheat tillering. This study provided valuable clues for functional characterization of GRF genes in wheat.


INTRODUCTION
Transcription factors are the most important regulators in plants. They are involved in various biological processes such as plant growth and development, metabolism, reproduction and differentiation (Chen & Cao, 2015;Kim & Tsukaya, 2015). Up to now, more than 60 transcription factors have been found in plants (Zhao et al., 2019). Growth-regulating factor (GRF) is a kind of plant-specific transcription factor, which regulates plant growth and development, plant cell size, participates in chloroplast proliferation and abiotic stress response (Chen & Rajewsky, 2007;Omidbakhshfard et al., 2015;Ruan et al., 2018). The most prominent feature of the GRF proteins is that there are two conserved domains in the N-terminal region, namely QLQ (Gln, Leu, Gln) and WRC (Trp, Arg, Cys) domains (Kim & Kende, 2004).
In recent years, many studies have demonstrated that GRFs play a variety of important regulatory roles in plant growth and development (Ma et al., 2017). GRFs participate in the growth and development of root (Bao et al., 2014), stem (Kim & Lee, 2006), leaf (Horiguchi, Kim & Tsukaya, 2005;Liu et al., 2012Liu et al., , 2015 and flower (Liang et al., 2014;Lee et al., 2015), and play an important role in signal transduction and stress response of plants (Liu et al., 2008;Hewezi et al., 2012;Casadevall et al., 2013). OsGRF4 controls grain shape, panicle length and seed shattering in rice (Sun et al., 2016). However, there is almost no evidence of the GRFs regulating tiller development in wheat. A recent study indicates that TaGRFs may be involved in regulating wheat (Triticum aestivum L.) tillering (Wang et al., 2019a).
Common wheat is one of the most important food crops worldwide. Previously, we have reported a dwarf-monoculm wheat mutant (dmc) derived from Guomai 301 (He et al., 2018), and found that wheat miR396b (tae-miR396b) and several of its target TaGRFs are involved in regulating wheat tillering in dmc . This provides a good opportunity to design experiments exploring the molecular mechanisms of wheat tiller development. Here, we attempt to thoroughly investigate the expression profiles of all TaGRFs in Guomai 301 and mutant dmc under normal growth and development condition, as well as exogenous indole acetic acid (IAA) and gibberellin acid (GA) applications. The results provided a theoretical base for further functional characterization of GRF genes in wheat.

Plant materials and growth conditions
The mutant dmc was obtained according to the method described in the previous article . All the plant materials were planted in our experimental field (34 25′ N, 115 39′ E, 49 m a.s.l.). The field experiments were carried out in a completely randomized design as described by Li et al. (2019). Fertilizer and weed management were similar to that of wheat breeding .

IAA and GA treatments
For hormone treatments, the wheat seedlings at the early three-leaf stage were sprayed with 10 −5 mol L −1 IAA solution, and the wheat seedlings at the three-leaf stage were sprayed with 2 × 10 −4 mol L −1 GA solution on the leaves until all the leaves were wet.
The tillers of seedlings were sampled at 0 (untreated control), 1 and 2 h after hormone treatments. The RNAs of all the samples were immediately extracted.

Tiller sample preparation
Three bulks of tiller samples were prepared separately at the three-leaf stage (WT1, dmc1; sampling date: 15 November 2018), the over-winter stage (WT2, dmc2; sampling date: 6 January 2019) and the rising to jointing stage (WT3, dmc3; sampling date: 16 February 2019) for RNA extraction. The RNAs were used for qRT-PCR analysis. Wheat tillering had been completed at the rising to jointing stage. All the samples had three biological replicates.

Identification of TaGRFs
The genome sequences and protein sequences of latest wheat genome assembly version IWGSC refseqv1.1 (http://plants.ensembl.org/) was used to identify wheat GRF family members. The longest transcript sequence corresponding to each candidate gene was selected as the final sequence.
The Hidden Markov Model (HMM) of WRC (PF08879) and QLQ (PF08880) domains were obtained from the Pfam website (http://pfam.xfam.org/) and were used to identify the conserved domains of wheat GRF proteins with HMMER software (https://www.ebi.ac. uk/Tools/hmmer/). Besides, the protein sequences of Arabidopsis GRF family members and rice GRF family members were downloaded from PlantTFDB (http://planttfdb.cbi. pku.edu.cn/), and these sequences were used as input sequences to BLAST in the wheat protein database. All output protein sequences with e-value ≤ 1 × 10 −10 were collected, removing the redundant sequences. The domains of the candidate GRF family members were confirmed by searching Pfam website and the Conserved Domain Search at NCBI (https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi). The TaGRF genes were named mainly to refer to the annotation information from the UniProt database (https://www.uniprot.org/) and their chromosomal locations.
The lengths of amino acid sequences, molecular weights, pIs and other characteristics of wheat GRF proteins were analyzed using the online analysis software of EXPASY website (https://www.expasy.org/). All identified TaGRF genes were located on chromosomes, and the gene map was drawn using TBtools software (Chen et al., 2020).

Construction of the phylogenetic tree
Amino acid sequences of the GRF proteins in wheat, Arabidopsis and rice were used to conduct multi-sequence alignment. Based on the results of sequence alignment, the phylogenetic tree was constructed using the Neighbor-Joining method in MEGA software (https://www.megasoftware.net/). The check parameter bootstrap value was 1,000, and the default value of the system was used for other parameters.

Conserved motif and gene structure analysis of TaGRFs
The online website MEME (http://meme-suite.org/) was used to analyze the motifs of TaGRFs protein sequences. Parameters were set as follows: the motif discovery mode was classic mode, the site distribution was Zero or One Occurrence Per Sequence (zoops), the maximum number of motifs to find was 10, and other parameters were default. The DNA and cDNA sequences corresponding to predicted proteins from the wheat genome database were downloaded from the wheat genome database. Then, TBtools software was used to draw motif distribution and the structure distribution map of TaGRFs, and the Logo diagram of amino acid conservation was drawn on the WEBLOGO website.

Analysis of the cis-acting elements in TaGRF promoters
The 2,000 bp upstream sequences before transcription start positions of TaGRFs were extracted from wheat genome sequence, and the cis-acting elements were analyzed using the PlantCARE database (http://bioinformatics.psb.ugent.be/webtools/plantcare/html).

Gene duplication and synteny analysis of TaGRFs
Multiple Collinearity Scan toolkit (MCScanX) program with the default parameters was adopted to analyze the gene duplication events of TaGRFs. To exhibit the synteny relationship of the orthologous GRF genes obtained from wheat and other selected species, the syntenic analysis maps were constructed using the Dual Synteny Plotter software (https://github.com/CJ-Chen/TBtools). According to the results of MCScanX, the nonsynonymous substitution rate (Ka) and synonymous substitution rate (Ks) of duplicated GRF genes were calculated by KaKs Calculator 2.0 software. The Ka/Ks ratios for TaGRF genes were used to assess the selection pressure on duplicated genes and Ka/Ks ratio >1, <1, or =1 indicated positive, negative, or neutral evolution, respectively.

Tissue specific expression analysis of TaGRFs
The raw gene expression data were downloaded from the Wheat Expression Browser (http://www.wheat-expression.com/) (Table S1). A total of eight mRNA-seq data from wheat cultivar Chinese Spring and Azhurnaya were analyzed. The samples were prepared from four tissues, including spike (anthesis stage), grain (milk grain stage), leaves/ shoots (flag leaf stage) and roots (flag leaf stage). TBtools software was used to draw the expression heat map of TaGRFs in wheat. The expression profiles were performed based on the transcripts per million (TPM) values of TaGRF genes. The Gene Ontology database (GO, geneontology.org) was used for functional annotation of the GRF genes. GO annotations were mapped according to molecular functions, biological processes and cellular components (Table S2).

RNA-seq and gene expression analysis
The samples for RNA-seq analysis were prepared at the three-leaf stage to four-leaf stage in 2017. Two super bulk samples of the mutant dmc (T01, T02 and T03) and WT (T04, T05 and T06) with three biological replicates were prepared. The transcript abundance of TaGRF genes was calculated as fragments per kilobase of exon model per million mapped reads (FPKM) (He et al., 2018). The expression profiles of all the TaGRF genes in Guomai 301 and mutant dmc were analyzed based on the RNA-seq data in this study (Table S3). The heatmaps were drawn by HemI1.0 based on the transformed data of log 2 (FPKM + 1) values (Deng et al., 2014).

qRT-PCR
qRT-PCR was performed as described previously by An et al. (2019). The primers of TaGRFs were designed using Primer-Blast of NCBI website (https://www.ncbi.nlm.nih. gov/tools/primer-blast/). All the primer sequences were listed in Table S4. The β-actin gene was used as an internal control and each reaction was performed with triplicates. The relative expression of TaGRFs was calculated by 2 −DDCt methods (Livak & Schmittgen, 2001).

Typical traits of Guomai 301 and mutant dmc
The plant height and tiller number between mutant dmc and Guomai 301 were significantly different at heading stage (Fig. 1A). The average plant height and tiller number of Guomai 301 were 64.9 cm and 20.33, the average plant height of mutant dmc was 48.0 cm and most dmc plants had no tiller. At the three-leaf stage (Fig. 1B), the tiller buds of WT formed two primary tillers (PTs) at the base of the main culm (Fig. 1E). Meanwhile, only one tiny protuberance formed at the main culm base of dmc ( Fig. 1H). At the over-winter stage (Fig. 1C), the tiller number of WT had reached 11-14, which mainly consisted of PTs and secondary tillers (Fig. 1F), while there were only two tiny tiller

Identification of TaGRFs
A total of 30 TaGRF genes (including 13 homoeologous groups) ( Table 1) were finally identified from the wheat genome using nine Arabidopsis GRF and 12 rice GRF genes as input sequences. The lengths of TaGRF proteins ranged from 211 to 611 aa and the mean value was 359 aa; the molecular weight ranged from 22.3 to 64.3 kDa and the average value was 38.8 kDa. The TaGRF6-4B was the heaviest and TaGRF10-6B was the lightest. The isoelectric point (pI) was between 4.72 and 9.90, with an average value of 7.9, and the TaGRF10-6D was the highest and TaGRF11-2B was the lowest. Among them, seven TaGRFs with pIs less than seven were slightly acidic, while the remaining 23 TaGRFs with pIs greater than seven were slightly alkaline. These data indicated that most TaGRFs contained more basic amino acids. The fat solubility of TaGRFs was between 45.83 and 70.6, with an average value of 57.056. Their hydrophobic indices ranged from −0.882 to −0.35, indicating that TaGRFs had good hydrophilicity (Grand average of hydropathicity, GRAVY). The 30 TaGRFs were located on 12 of the 21 wheat chromosomes. Each of the 12 chromosomes had 1-4 TaGRFs. Among them, chromosome 6A had four TaGRFs, 4B and 7B had one TaGRF and 4D, 7A and 7D had two TaGRFs, respectively (Fig. 2).

Phylogenetic tree of wheat GRF proteins
In order to better understand the phylogenetic relationship between TaGRFs and those in other plant species. The phylogenetic tree was constructed (Fig. 3), including 30 TaGRFs, nine AtGRFs and 12 OsGRFs. According to the branches of the phylogenetic tree, the 51 GRFs were clustered into four classes: Class I, II, III and IV. There were nine TaGRFs in Class I, six TaGRFs in Class II, five TaGRFs in Class III and 10 TaGRFs in Class IV. Phylogenetic analysis indicated that wheat GRFs were more closely related to rice GRFs than Arabidopsis GRFs.

Structures and conserved motifs of wheat GRF genes
Based on the gene structural information, the phylogenetic tree of TaGRFs was built (Fig. 4A). A total of 10 conservative motifs were predicted in wheat GRF proteins (Fig. 4B). Motif 1 was the WRC domain and motif 2 was the QLQ domain, both were highly ). In addition, most TaGRFs also contained motif 3. The exon-intron structure diagram of TaGRFs showed that their genomic DNA sequence lengths were significantly different (Fig. 4C). The longest was TaGRF4-6A, its length was about 7,000 bp, and the shortest was about 1,000 bp. The exon number of TaGRFs was 2-5, and the exon numbers of TaGRFs were closely related to their classes in the phylogenetic tree of TaGRFs (Fig. 3). For example, each of the six TaGRFs in Class II had two exons. Among the 30 TaGRFs, 25 TaGRFs had both 5′-UTR (untranslated region) and 3′-UTR, two TaGRFs only contained 3′-UTR, one TaGRF only contained 5′-UTR, and the remaining two TaGRFs had no UTRs.

Cis-acting elements in the promoters of TaGRFs and GO annotation
Among the members of the wheat GRF gene family, the promoter sequences of TaGRF1-6B, TaGRF10-6A, TaGRF11-2A, TaGRF11-2B and TaGRF11-2D contained a large number of 'N', so they hadn't been analyzed (Fig. 5). Except for a large number of CAAT-box and TATA-box elements, there were also a large number of cis-acting regulatory elements related to growth and development, hormones and stress responsiveness in the promoter sequences of TaGRFs, including some cis-acting elements involving in auxin (AuxRR-core, TGA-element), gibberellin (P-box), light response (ACE, AEbox and G-box), methyl jasmonate reaction (CGTCA-motif), salicylic acid reaction (TCA-element), circadian rhythm control (circadian), regulation of gliadin metabolism (O2-site), and regulation of meristematic expression (CAT-box and CCGTCC motif). These cis-acting elements implied TaGRFs played various potential roles in regulating wheat growth and development, and response to hormones and stresses. Gene Ontology annotation revealed that the TaGRF genes are mainly involved in regulation of transcription, DNA-templated (GO:0006355) and response to deep water (GO:0030912). Both TaGRF1 and TaGRF2 are involved in response to gibberellin (GO:0009739) (Table S2).

Synteny analysis of TaGRFs
A total of 39 TaGRF gene pairs were probably segmental duplicated genes, including 21 pairs (24 genes) of homoeologous genes, and they distributed on different chromosomes ( Fig. 6; Table S5). According to the methodology of Holub (2001), there were no tandem duplication genes in wheat GRF gene family, and most TaGRFs were generated by large-scale repeated events or segmental duplications. It implied that the segmental duplication played an important role in the evolution of the TaGRFs. Three comparative syntenic maps of wheat were constructed compared with three representative plant species, including Arabidopsis, rice and maize (Fig. 7). Only one wheat GRF gene, TaGRF6-4B, had a syntenic gene in Arabidopsis (AtGRF1) (Fig. 7A). The numbers of syntenic gene pairs between wheat and rice (Fig. 7B), wheat and maize (Fig. 7C) were 40 and 51 (Table S6). It indicated that the evolutions of GRF genes in wheat, rice and maize were similar.
To better understand the evolutionary constraints acting on TaGRFs, we calculated the Ka/Ks ratios of these syntenic gene pairs. The results showed that the Ka/Ks ratios of most TaGRF syntenic gene pairs were less than 1 (Table S7), suggesting that the TaGRFs might have undergone purifying selection processes.

The expression patterns of TaGRFs in different tissues
The expression profiles of TaGRFs in spike, grain, leaf/shoot and root showed that about half of 30 TaGRFs had similar expression patterns in Chinese Spring and Azhurnaya, and the others expressed differentially in the four tissues of Chinese Spring (Fig. 8A). The expression patterns of most TaGRFs in Azhurnaya were similar to that in Chinese Spring. However, the expression profiles of TaGRF1, TaGRF6 and TaGRF9 were significantly different between Azhurnaya and Chinese Spring (Fig. 8B). The expression levels of all the three homoeologous genes of TaGRF11-2(A, B, D) in the four tissues were high, while the expressions of TaGRF2, TaGRF4, TaGRF10 and TaGRF12 were not detected or their levels were very low, implying they played no important roles in the development of the four tissues. TaGRF1, TaGRF5, TaGRF6 and TaGRF9 expressed highly in spike and root, implying their roles in regulating spike and root development. Tissue-specific expression profiles of some TaGRFs between Azhurnaya and Chinese Spring were slightly different, which was probably caused by genotype, sample, and experiment variations. These results can be used as a reference for functional studies of TaGRFs.

Expressions of TaGRFs in TPs of the mutant dmc
The expression patterns of all 30 TaGRFs were investigated using available transcriptome data (Fig. 8C). Among the 30 TaGRFs, TaGRF1-6D and TaGRF5-7A had not been detected, TaGRF10-6D was almost not expressed, and TaGRF12-2A expressed very lowly in all detected samples, these four genes probably were not necessary for wheat tiller development (Fig. 8C). TaGRF1-6A (FPKM > 20 in dmc), TaGRF9-4A and TaGRF11-2A (FPKM > 24 in dmc) expressed relatively highly indicated their basic functions in wheat tiller development. All TaGRFs expressed at a reduced level in dmc than in WT. Compared to WT, the expression level of TaGRF6-4A decreased by 66.2%, TaGRF12-2B decreased by 57.8%, and other 15 TaGRFs decreased by more than 40% in dmc. This indicated that the expressions of all the TaGRFs played positive roles at early tillering stage, the constrained tillering of the dmc was associated with the lower expression levels of TaGRFs. qRT-PCR was performed to analyze the expression patterns of TaGRFs in the TPs of WT and mutant dmc at three tiller developmental stages (Fig. 9). The 30 TaGRF genes had various expression patterns at three tillering stages. All the 30 TaGRFs expressed at a reduced level in dmc at the three-leaf stage, which was consistent with the transcriptome data. A total of 20 TaGRFs expressed highly in dmc at the over-winter stage and the rising to jointing stage, while TaGRF12-2A and TaGRF12-2B expressed lowly in dmc at the rising to jointing stage. Among these genes, TaGRF1-6D, TaGRF2-7D, TaGRF4-6A, TaGRF4-6D, TaGRF5-7A, TaGRF6-4B, TaGRF6-4D and TaGRF12-2D expressed significantly differentially in dmc and WT at wheat tillering stages, implying their important roles in the abnormal tillering in dmc. The expression levels of TaGRF2-7A, TaGRF2-7B, TaGRF5-4A, TaGRF5-7D, TaGRF9-4A, TaGRF9-4D, TaGRF10-6A and TaGRF11-2A were always lower in dmc at two tillering stages, the expression levels of TaGRF2-7A, TaGRF2-7B and TaGRF11-2A in dmc were significantly different from that in WT. TaGRF9-4D and TaGRF11-2A expressed lowly in dmc compared to that in WT at each stage, implying they might affect the tillering of dmc. There were 15 TaGRFs expressed at an upward trend in dmc. TaGRF2-6A, TaGRF9-4A and TaGRF9-4D expressed at an upward trend in WT, implying their roles in regulating tiller development.
TaGRF10-6D and TaGRF12-2B expressed at a downward trend in dmc. The expression patterns of TaGRF genes were complex. In summary, the abnormal expressions of TaGRF2-7(A, B, D), TaGRF5-7D, TaGRF10-6(A, B, D) and TaGRF11-2A were major causes for constraining the tillering of dmc.

Expression patterns of TaGRFs in response to IAA and GA treatments
A total of 20 representative TaGRFs were selected from the 30 TaGRF genes to investigate whether their expressions were affected by IAA and GA treatments. These TaGRF genes were carefully selected based on the cis-acting elements in their promoters (Fig. 5) and their expression levels in tillers of dmc (Fig. 9). The expressions of 20 TaGRFs were affected by IAA treatment in various degrees (Fig. 10). The expressions of 9 TaGRFs in WT and 2 TaGRFs (TaGRF1-6D and TaGRF10-6B) in dmc were continuously significantly repressed by IAA treatment, their expressions decreased by more than 50%. Among these genes, there was a cis-acting element involving in auxin (AuxRR-core, TGA-element) in the 2,000 bp upstream sequences before transcription start positions of TaGRF1-6D, TaGRF3-2A, TaGRF4-6B, TaGRF6-4A, TaGRF9-4D and TaGRF12-2D. The expression TaGRF9-4A and TaGRF10-6D in dmc were significantly up-regulated by IAA treatment, which suggested these TaGRFs played key roles in regulating wheat tillering. The expressions of the 20 TaGRFs were also affected by GA treatment in various degrees (Fig. 11). It was worth mentioning that most of the homoeologous genes showed similar expression patterns in response to GA treatment, such as TaGRF4, TaGRF5, TaGRF6, TaGRF9 and TaGRF10. The expressions of TaGRF1-6D, TaGRF6-4B and TaGRF9-4D in WT were hardly affected by GA treatment, while they were significantly up-regulated in dmc at 2 h after GA treatment. There were 10 TaGRFs in dmc were significantly repressed at 1 h and significantly up-regulated at 2 h after GA treatment. The expressions of all the detected TaGRFs were significantly up-regulated in dmc after GA treatment. There were 1-3 cis-acting elements (indicated in the parentheses behind the gene symbols) involving in gibberellin (P-box) in the 2,000 bp upstream sequences before transcription start positions of TaGRF5-7A (1), TaGRF6-4A (3), TaGRF6-4B (2), TaGRF6-4D (2), TaGRF9-4D (1) and TaGRF10-6D (1), respectively. The expressions of all TaGRFs were affected by both IAA and GA treatments indicated that phytohormone IAA and GA were involved in regulating wheat tillering.

Characteristics and evolution of wheat GRF gene family members
Plant GRF proteins typically have two conserved domains, QLQ and WRC, in the N-terminal regions. The QLQ domain interacts with GRF-interacting factors (GIFs) to form transcription activating factor that participates in the biological processes of plant growth and development (Kim & Kende, 2004). The WRC domain interacts with the cis-acting regions of its downstream genes and plays the biological function of GRF. The genome-wide analyses of GRF gene families have been widely carried out in various plant species whose genomes have been sequenced (Cao et al., 2016). In this study, a total of 30 wheat GRF proteins with the two conserved domains, QLQ and WRC, were identified. Phylogenetic and gene structure analyses showed that most of the GRF genes in the same subfamily had similar exon/intron structures, which provided clues to the evolutionary relationships of plant GRFs (Hu & Liu, 2011). These data indicated that the GRF genes with similar structures have similar evolution histories and functions (Babenko et al., 2004;Roy & Penny, 2007). Here, we found wheat homoeologous genes of TaGRF3-2(A, B, D) and TaGRF10-6(A, B, D) were highly homologous to OsGRF3 and OsGRF10. The distributions of the introns and exons in TaGRF genes were similar to those of rice GRF genes (Cao et al., 2016), suggesting that TaGRF genes and OsGRF genes might have similar functions. The gene duplication events are usually derived from the polyploidization or tandem and segmental duplication, segmental duplication occurs most frequently in plants because most plant species are diploidized polyploids (Moore & Purugganan, 2005;Adams & Wendel, 2005). No tandem repeat genes were found in TaGRFs, there were 24 homoeologous genes in 27 segmental duplication GRF genes, and not every GRF had three homoeologous genes on the homoeologous chromosomes A, B and D. It had been reported that some homoeologous genes could be lost during the polyploidization of the genome (Lynch & Force, 2000). Hence, it was highly possible that the segmental duplication played an important role in the expansion of the TaGRFs. A large number of cis-acting elements related to growth and development, hormones and stress regulation were found in the promoter regions of TaGRF genes, which implied their various functions.

Various functions of TaGRFs
Plant GRFs regulate growth and development. AtGRF5 stimulates Arabidopsis chloroplast division, photosynthesis and leaf longevity (Vercruyssen et al., 2015). Overexpression of ZmGRF10 can reduce leaf size and plant height of maize . OsGRF1 not only regulates growth at the juvenile stage, but also regulates heading in rice (Luo et al., 2005). AtGRF4 and AtGRF5 are involved in the development of leaf size and leaf senescence in Arabidopsis (Horiguchi, Kim & Tsukaya, 2005;Kim & Lee, 2006). As the target genes of miR396, the expressions of GRF genes are post-transcriptionally down-regulated by miR396 (Hewezi & Baum, 2012;Lee & Kim, 2014). There are six GRF genes (AtGRF1-AtGRF3 and AtGRF7-AtGRF9) are regulated by miR396 (miR396a and miR396b) in Arabidopsis (Liu et al., 2009). Down-regulated expressions of AtGRF genes by miR396 overexpression resulted in Arabidopsis with enhanced leaf adaxial-abaxial defects, narrower leaves, abnormal development of the pistil and compromises the shoot meristem, etc (Jones-Rhoades, Bartel & Bartel, 2006;Rodriguez et al., 2010Rodriguez et al., , 2015Wang et al., 2011). Overexpression of Ptc-miR396 from Populus trichocarpa in tobacco resulted in altered plant growth and flower development (Baucher et al., 2013). miR396 plays an important role in controlling carpel number and pistil development via regulation of the GRF/GIF complex (Liang et al., 2014). OsmiR396d regulates the expressions of OsGRF genes and participates in the regulation of flower organ development in rice (Liu et al., 2014).
The functional diversity of genes can be predicted by their tissue-specific expression patterns. GRF genes not only have different distributions and molecular structures, but also have different tissue-specific expression patterns. Studies have shown that TaGRF genes are highly expressed in active development tissues or organs, and they are relatively lowly expressed in mature tissues or organs (Choi, Kim & Kende, 2004;Wu, Wang & Zhuang, 2017). TaGRFs express highly in young spikes and roots, which is consistent with the putative function based on the meristem-related cis-acting elements in the promoter regions of TaGRFs. AtGRF1 and AtGRF3 are the most highly expressed genes in Arabidopsis roots, which are essential for root development (Szakasits et al., 2009;Hewezi et al., 2012). Overexpression of AtMIR396a decreased the transcript levels of AtGRF genes and resulted in a shorter root phenotype (Bao et al., 2014). TaGRF11 and its homoeologous genes (TaGRF11-2A, 2B, 2D) were highly expressed in various tissues, indicating they played essential roles in wheat growth and development (Fig. 8). TaGRF1-6A, TaGRF6-4A and their homoeologous genes expressed highly in roots suggesting their functions in root development (Fig. 8). The data indicated that TaGRFs functioned differently during wheat development.
The key TaGRFs involved in tiller development Plant miR396/GRF is a conservative plant growth regulation module, but there is no definite research report on its involvement in plant branching and tillering (Liebsch & Palatnik, 2020). Overexpressions of OsGRF3 and OsGRF10 reduce formations of tillers and internodes in rice (Kuijt et al., 2014). Therefore, we speculated that the homologous genes of TaGRF3 and TaGRF10 might be related to the growth and development of tiller in wheat. The tiller ability of the transgenic wheat overexpression miR396 is significantly decreased (Song, 2016). Our miRNome and transcriptome integrative analysis about the mutant dmc and WT found that the highly expressed tae-miR396b (T. aestivum microRNA396b) down-regulated the expressions of many TaGRFs (TaGRF1, TaGRF2, TaGRF5, TaGRF9, TaGRF10 and TaGRF12) in dmc during tillering (He et al., 2018;An et al., 2019). It demonstrated that the miR396/GRF regulatory module played a key role in wheat tiller development. Compared with the WT, the expression levels of all TaGRFs in dmc were significantly decreased at early tillering stage, which is positively related to the phenotype of dmc (Fig. 9). All the 30 TaGRFs have different expression patterns in WT and dmc, but only those significantly differentially expressed TaGRFs in TPs are the key tiller development regulators. In this case, TaGRF1-6D, TaGRF2, TaGRF4, TaGRF5-7A, TaGRF11-2B and TaGRF12-2D were significantly differentially expressed at early tillering stage, indicating their important roles in regulating tiller numbers in wheat. In summary, TaGRF2-7(A, B, D), TaGRF5-7D, TaGRF10-6(A, B, D) and TaGRF11-2A were key regulators in wheat tiller growth and development, especially in dmc.

Exogenous IAA and GA affect the expressions of TaGRFs
There is considerable evidence that GRFs play significant roles in regulating plant growth and development, and play a variety of regulatory roles in signal transduction and stress response of plants. There are a large number of cis-acting elements related to hormones in TaGRF promoters, including those related to IAA (AuxRR-core, TGAelement) and GA (P-box). IAA regulates plant morphological formation, such as plant tropism growth, root development and tiller formation (Li et al., 2000;Ren, Dai & Liu, 2012). OsmiR167a represses its targets OsARF12 (auxin response factor, ARF), OsARF17 and OsARF25, to control rice tiller angle by fine-tuning auxin asymmetric distribution in shoots (Li et al., 2020). GA participates in regulation of many physiological processes in wheat including the spike development, stem elongation, plant height and stress responses (Zhang et al., 2007;Pavlista, Santra & Baltensperger, 2013). KNOX (KNOTTED1-like homeobox) proteins contribute to the regulation of meristem maintenance by negatively regulating the production of GA (Sakamoto et al., 2001;Bolduc & Hake, 2009). Repression the activities of KNOX genes is a conserved function of GRFs (Kuijt et al., 2014). These results suggested that the GRF genes positively regulated the production of GAs, and GAs in turn up-regulated GRF gene expressions . OsGRF1 may play a regulatory role in GA-induced stem elongation (Van der Knaap, Kim & Kende, 2000). GA treatment increases the expressions of some OsGRFs (OsGRF1, 2, 3, 7, 8, 10 and 12) in rice, while the expression of OsGRF9 decreases (Bao et al., 2014). Interactions between exogenous cytokinin and nitrogen treatment can regulate tiller bud growth in winter wheat . There are many evidences of IAA and GA affect tiller developments (Lee, Foster & Morgan, 1998;Frantz et al., 2004;Liu et al., 2011;Assuero et al., 2012;Cai et al., 2013). In previous studies, we found that the contents of IAA were significantly lower, but the contents of GA were significantly higher in the tiller tissues of dmc . The results of qRT-PCR also confirmed that the expressions of TaGRFs were significantly affected by IAA and GA applications (Figs. 10 and 11). According to these data, it was considered that TaGRFs as well as IAA and GA signaling were involved in regulating wheat tiller development.

CONCLUSIONS
A total of 30 TaGRF transcription factors with both typical QLQ and WRC conserved domains were identified in wheat. TaGRF genes distribute on 12 of the 21 wheat chromosomes, and their promoter regions have a large number of cis-acting elements related to plant growth and development, hormone signal pathway and stress response. The expression levels of all the TaGRF genes are significantly lower in dmc at the early tillering stage (three-leaf stage), and the expressions also are significantly affected by exogenous IAA and GA. TaGRF1-6A, TaGRF2-7D, TaGRF9-4A and TaGRF11 Yumei Jiang performed the experiments, authored or reviewed drafts of the paper, contributed to the field experiments and maintenance of the wheat accessions, and approved the final draft. Zhixin Jiao performed the experiments, prepared figures and/or tables, helped with sowing and sample preparation, and approved the final draft. Huijuan Li performed the experiments, prepared figures and/or tables, and approved the final draft. Ting Wang performed the experiments, prepared figures and/or tables, helped with sowing and sample preparation, and approved the final draft. Peipei Zhang performed the experiments, prepared figures and/or tables, helped with sowing and sample preparation, and approved the final draft. Mengyao Han performed the experiments, prepared figures and/or tables, helped with sowing and sample preparation, and approved the final draft. Lei Li conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft. Hongjie Liu performed the experiments, authored or reviewed drafts of the paper, selected the mutant, and approved the final draft. Qiaoyun Li performed the experiments, authored or reviewed drafts of the paper, contributed to the field experiments and maintenance of the wheat accessions, and approved the final draft. Jishan Niu conceived and designed the experiments, authored or reviewed drafts of the paper, and approved the final draft.

DNA Deposition
The following information was supplied regarding the deposition of DNA sequences: Raw sequencing reads are available at NCBI under BioProject ID PRJNA670838.

Data Availability
The following information was supplied regarding data availability: The raw measurements for qRT-PCR are available in the Supplemental File.

Supplemental Information
Supplemental information for this article can be found online at http://dx.doi.org/10.7717/ peerj.11235#supplemental-information.