Genome-wide analysis of the GRF gene family and their expression profiling in peach (Prunus persica)

ABSTRACT Growth-regulating factors (GRFs) are plant-speciﬁc transcription factors with vital roles in multiple biological processes. Although GRFs have been identified in various plant species, a comprehensive analysis of GRF genes in peach (Prunus persica) has not yet been reported. In this study, 10 PpGRF genes distributed on 6 chromosomes were identified in peach genome and their properties were analyzed systematically. Expression pattern analysis suggested that most PpGRFs were preferentially expressed in young tissues. Multiple types of cis-elements were observed in PpGRF promoters, and PpGRFs positively respond to ultraviolet B-rays (UVB) and gibberellin (GA) treatments at the transcriptional level. Also, the content of gibberellic acid 3 (GA3) and indole-3-acetic acid (IAA) changed significantly after UVB irradiation. PpGRF 3, 4, 5, 6, 7, 9 and 10 positive responses to UVB and GA3 signals. The evolutionary patterns and expression profiles of PpGRFs detected in this study increase understanding of the important roles in peach.


Introduction
Growth-regulating factors (GRFs) are plant-specific transcription factors (TFs) that play pivotal roles in plant growth and development. The first GRF gene was identified in rice and was shown to mediate gibberellin (GA)-induced stem growth regulation in deepwater rice (Van et al. 2000). The GRF family proteins contain evolutionally conserved QLQ (glutamine, leucine, glutamine; Gln, Leu, Gln) and WRC (tryptophan, arginine, cysteine; Try, Arg, Cys) domains in their N-terminal region ( Van et al. 2000;Kim et al. 2003;Choi et al. 2004). The amino acid sequence of the QLQ domain is similar to the protein-protein interaction domain of SWITCH2/SUCROSE NONFERMENTING2 (SWI/SNF) SWI2/SNF2 (Treich et al. 1995). The QLQ domain is thought to help GRF proteins interact with growth regulators (GRFinteracting factors, GIFs) to regulate plant growth and development (Kim et al. 2003;Chen et al. 2019). The WRC domain is a plant-specific domain capable of binding to DNA and has a nuclear localization signal and a zinc finger motif ( Van et al. 2000;Kim et al. 2003). Also, most GRF proteins have transactivation activity in the C-terminal region and share some less conserved motifs, such as the TQL (threonine, glutamine, leucine; Thr, Gln, Leu), GGPL (glycine, glycine, proline, leucine; Gly, Gly, Pro, Leu), and FFD structure domains (Zhang et al. 2008;Chen et al. 2019).
As development-related TFs, GRFs are strongly expressed in actively growing and developing tissues, involved in root and leaf development, stem elongation, meristem maintenance, cell expansion, floral organogenesis, and seed formation (Kim et al. 2003;Horiguchi et al. 2005;Debernardi et al. 2014;Lee et al. 2015). Leaf and cotyledon sizes are increased in Arabidopsis lines overexpressing AtGRF1, 2. This finding suggests that the AtGRF1, 2 genes have roles in leaf and cotyledon development (Kim et al. 2003). At first, only the Atgrf1/2/3 triple mutant had smaller, narrower leaves and shorter petioles than the wild type, whereas single mutants displayed no changes in phenotype and double mutants displayed only minor phenotypic changes (Kim et al. 2003). Later reports demonstrated that AtGRF4 is required to develop leaves, cotyledons, and shoot apical meristem by constructing a quadruple mutant by crossing it to the grf1/2/3 triple mutant. The quadruple mutant exhibited much smaller leaves and severe fusion of cotyledons (Kim and Lee 2006). Interestingly, in contrast to AtGRF1 and AtGRF2, the overexpression of AtGRF5 showed a larger leaf area due to an increased number of cells (Horiguchi et al. 2005). Also, the biochemical analysis showed that GRF and GIF combine to form a transcription complex in vivo. The GRF-GIF complex might activate its own transcription through a positive feedback regulatory loop, thereby regulating cell proliferation and ultimately controlling leaf size (Horiguchi et al. 2005;Kim and Lee 2006;Lee et al. 2009;Debernardi et al. 2014;Omidbakhshfard et al. 2015). Overexpression of TCPs (AtTCP4) leads to reduced levels of GIF transcripts, suggesting that they act as mediators in the GIF-GRF-miR396 regulatory network (Schommer et al. 2014;Kim and Tsukaya 2015).
GAs are a major class of plant hormones that control many developmental processes. Previous studies demonstrated that GA regulated the transcriptional expression of the GRF genes (Choi et al. 2004;Ma et al. 2017). After GA 3 treatment of celery cabbage leaves, the transcript levels of several Chinese cabbage GRFs (BrGRF2,4,5,7,8,9,11,12,15,and 16) were enhanced (Ma et al. 2017). Similarly, GA 3 treatment enhanced the transcript levels of several rice GRFs (OsGRF 1,2,3,7,8,10,and 12), while the transcript levels of OsGRF9 decreased (Choi et al. 2004). However, few studies have focused on effect of GA on transcript levels of GRFs in peach, making it crucial to explore the relationship between GRFs and GA.
MicroRNAs (miRNAs) are comprised of approximately 21 nucleotide (nt) RNAs that play important regulatory roles in growth and development by targeting mRNAs for cleavage or translational repression. MiR396 is a miRNA conserved among the dicot and monocot plants. It has been shown that miRNA396 directly inhibits the expression of GRFs in various plant species through post-transcriptional regulation (Liang et al. 2014;Chen et al. 2019). The heterologous expression of Ath-miR396a (from Arabidopsis thaliana) in tobacco causes a reduction in the transcription level of the GRF gene in tobacco, accompanied by similar narrowleaf phenotypes such as reduced leaf size. Also, Ath-miR396a overexpression affects flower development relative to wild type, for example, increasing the number of petals, anthers, and carpels and producing shortened (curved) stamens (Yang et al. 2009). Similarly, ectopic expression of Ath-miR396 resulted in changes in tobacco growth and flower development, accompanied by a decrease in the abundance of NtGRF-like (FG137771, FG165999, FG167390, FG194569) transcripts (Yang et al. 2009). Not all GRFs are miR396 target genes. In Arabidopsis, miR396 targets seven GRF genes (AtGRF 1, 2, 3, 4 and AtGRF 7, 8, 9) encoding putative transcription factors with roles in plant leaf growth (Hoe and Tsukaya 2015). Interestingly, overexpression of miR396 also inhibits AtGRF 6 expression, which is not a target of miR396 Rodriguez et al. 2010). MiR396 expression is caused by various abiotic stresses, such as high salinity, low temperature, drought stress, and ultraviolet B (UVB).
UVB is ultraviolet radiation, accounting for about 1.5% of the total solar radiation. Lowdose UVB is an environmental regulator affecting gene expression, cellular and metabolic activities, and growth and development (Brosche and Strid 2003;Jenkins 2009). Several reports have focused on the correlation between UVB and GRFs (Casadevall et al. 2013;Fina et al. 2017). Casadevall et al. (2013) revealed that UVB exposure affects leaf growth by inhibiting cell division in proliferating tissues, a process mediated by miR396 and GRFs in Arabidopsis. Fina et al. (2017) showed that UVB radiation limits leaf growth through mediating GRF activity and GA levels in maize. Soybean plant morphology is substantially altered in response to shade stress (with reduced UVB radiation). At the same time, under shade stress conditions, the GmGRF9 and GmGRF17 expression levels were strongly suppressed, and GmGRF5 was up-regulated (Chen et al. 2019). As a new agricultural form, protected cultivation has been rapidly developed. The protected cultivation of fruit is more economical than traditional cultivation, but it limits the growth of fruit trees, yield, and fruit flavor . Commonly, UVB radiation is insufficient in protected cultivation fruit facilities. A better understanding of UVB radiation and GRFs with other regulatory factors (such as miR396, GA) is crucial for fruit development.
Peach (Prunus persica) is native to China, where they have a cultivation history of more than 3000 years (Zheng et al. 2014). Peach has a relatively small genome and high economic and nutritional value. In 2010, the Genome Sequencing Project of the peach double haploid cultivar'Lovell'was completed, generating 230Mb genome sequence and 202 assembly scaffolds (Verde et al. 2013). Therefore, the peach is considered a useful forest model species for genetic and ecological research (Wang et al. 2018). Many studies have focused on the regulation of Prunus persica genes by TFs, including WRKY, NF-Y, MYB, AP2, TCP, NAC, and SPL TFs through evolutionary analyses Sun et al. 2016;Zhang et al. 2018Zhang et al. , 2012Li et al. 2019aLi et al. , 2019b. However, few studies have focused on analysis of PpGRFs, making it interesting to explore the functions of these genes in peach.
In this study, 10 members of the PpGRF family in peach (Prunus persica) were identified and annotated. We analyzed the phylogenetic relationships among GRF genes, chromosomal location, their protein motifs, and gene structures. Furthermore, we verified gene expression patterns by transcriptome analysis and quantitative real-time PCR (qRT-PCR) in peach and identified several candidate GRF genes closely associated with peach. PpGRFs positively respond to UVB and GA treatments at the transcriptional level. These results in this study increase understanding of the important roles played by these genes in peach.

Plant material and treatments
The peach cultivars (Prunus persica var. nectarine Zhongyou No.5) used in this study were collected from Tai'an, China (36°18'N, 117°13'E). All trees were cultivated for 7 years under standard horticultural practices and were fully productive. The following two treatments were set up in this study: (1) UVB treament, (2) GA treatment. As our previous research described , the optimal UVB treatments were carried out with 1.44 Kj m −2 d −1 UVB radiation. UVB was provided by the dedicated UV lamp of 40W, 297 nm (NanjingKazhi), hanging at 1.5 m above the plants. UV light was kept on from 9:30 am to 10:30 am in this treatment. The leaves within the range of 80-120 cm below the lamp were selected. And as our previous research described , hormone treatments were conducted by spraying young leaves with 100 µM GA 3 . Leaf samples were collected at 3, 6, 9, and 24 h post-treatment. At each time point, we randomly collected 10 young leaves, which were immediately frozen in liquid nitrogen and stored at −80°C until use. In the tissue expression pattern experiment, we collected samples of five tissues, including root, stem, leaves, shoot tip, and hypocotyl (the hypocotyl of the cotyledons of nectarine seeds after germination), which were immediately frozen in liquid nitrogen and stored at −80°C until use.

Phylogenetic analysis and sequence alignment
To investigate the phylogenetic relationship of the GRF proteins between model plants and Rosaceae species, full-length GRF protein sequences of Arabidopsis, soybean, rice, poplar and four Rosaceae species (Prunus persica, Pyrus communis, Fragaria vesca and Malus domestica) were downloaded to construct a phylogenetic tree. Multiple sequence alignments of GRFs were analyzed by Clustal W. Then MEGA 7.0 software was employed to construct phylogenetic trees using the neighbor-joining (NJ) method and the bootstrap analysis (1000 times) (Tamura et al. 2013). The conservative motifs were predicted by the MEME online program (http:// meme-suite.org/tools/meme), in which the maximum value of the motif was set to 10. The intron-exon arrangements were determined by aligning its cDNA sequence with the genomic sequence by the TBtools software package. Then, the TBtools software package was used to integrate these analyses ).

Chromosomal location and collinearity analysis
MCScanX (http://chibba.pgml.uga.edu/mcscan2/) was used to analyze the WGD/segmental, tandem, proximal, and dispersed duplications to investigate the mechanisms mediating GRF gene family evolution. According to the results of MCScanX, the Ka, Ks and Ka/Ks of duplicated genes were calculated using the KaKs_Calculator 2.0 software (Zhang et al. 2006). Synteny analysis among the four Rosaceae genomes was conducted locally using the TBtools software package to map the genes to chromosomal locations Xie et al. 2018).

Cis-elements in the promoters of the PpGRFs analysis
The promoter sequences of the PpGRFs were identified by searching the P. persica database, and the cis-elements in the promoters were predicted using PlantCARE (Lescot et al. 2020).

RNA extraction and quantitative RT-PCR analysis
Total RNA was isolated from various tissues using the RNeasy Plus Mini Kit (Qiagen, Valencia, CA, USA) according to the manufacturer's instructions. RNA quality and quantity were determined by agarose gel electrophoresis and a NanoPhotometer P-360 (Implen, Munich, Germany) . To ensure that the RNA was free of DNA, three duplicate samples of total RNA (approximately 3 mg each) were reverse-transcribed with the PrimeScript RT reagent kit and gDNA Eraser (Perfect Real Time) (Takara Biotechnology, Dalian, China). The primers are listed in Additional File 1 (Supplementary Table S1) using Ppactin gene as the reference and designed using Primer Premier 6.22. RT-qPCR was performed using a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, California, USA) with the TB Green® Premix Ex Taq™ (TaKaRa Biotechnology). The protocol consisted of predenaturation at 95°C for 30 s, followed by 40 cycles of denaturation at 95°C for 5 s, primer annealing/extension at 60-63°C for 30 s, and detection of the fluorescent signal. Primer specificity was verified by melting curve analysis, subjecting the samples to 95°C for 10 s, 60°C for 5 s, and 95°C for 0.5 s. Each sample was analyzed in triplicate. Then, the 2 −△△ct method was used to determine expression levels relative to the reference genes (Livak and Schmittgen 2001).

Determination of growth index
The measurement of growth index started on June 7 and was measured once per month until the peach tree stopped growing and enters the dormant induction period on October 7. During the measurement, five peach trees with uniform growth in each experimental group were selected. The length and thickness of the new shoots and the plant height were repeatedly measured, and the average was calculated.

Determination of hormone contents
Using liquid chromatography-mass spectrometry (LC-MS) to determine the different hormones [indoleacetic acid (IAA), abscisic acid (ABA), gibberellic acid 3(GA 3 ), jasmonic acid (JA), trans-zeatin riboside (TZR)] contents of the stem tip tissues under UVB treatment and CK were taken on August 7 (vigorous shoot growth) and stored at −80°C until use (Mikiko et al. 2009). The ground vigorous shoot growth samples (0.8 g/replication) were transferred to a centrifuge tube containing 8 mL isopropyl alcohol/hydrochloric acid (IPA/HCl) buffer shaken for 30 min at 4°C. 15 ml dichloromethane (CH 2 Cl 2 ) was added shaken for 30 min at 4°C. The mixture was centrifuged at 13000 rpm for 5 min at 4°C to get organic phase. Then the organic phase was dried with a nitrogen blower in dark. 400 μl methanol (0.1% formic acid) was added to dissolve the above organic phase. Finally the mixture was filtered with 0.22 μm organic mesofiltration membrane for tests of hormone contents using LC-MS.
2.9. Target prediction of miR396 and expression analysis following UVB exposure Target prediction of miR396 followed the rules referring to . Target sites were searched using peach genome information. To better understand the roles of miRNAs in peach under low-dose UVB treatment, the potential target functions were annotated using the Gene Ontology and KEGG pathway database.

Phylogenetic and structural analyses of PpGRF genes
To gain an insight into the evolutionary relationships, a neighbor-joining phylogenetic tree was constructed among GRFs from Arabidopsis (9 members), soybean (22 members), rice (12 members), poplar (19 members), and four Rosaceae species (46 members, as mentioned above in Table 1). As shown in Figure 1, 108 GRFs from the different species were divided into six subgroups (I-VI). The majority [∼23% (25/108)] of GRF proteins belonged to clade I, 21 to clade II (without AtGRF), 24 to clade III, 15 to clade IV (without OsGRF), 13 to Clade V (without OsGRF), and 10 to clade VI (without GmGRF) (Figure 1). No species-specific clades were discovered, and every subgroup contained PpGRFs, same as other three Rosaceae species (Pyrus communis, Fragaria vesca, Malus domestica). The phylogenetic tree showed that the PpGRFs showed a closer relationship with PcGRFs, FvGRFs, and MdGRFs, which suggests that GRF proteins have been conserved during evolution in the four Rosaceae species.
To better explore the evolutionary relationships in peach and other three Rosaceae species and predict the function of GRF proteins, we investigated the exon-intron patterns and motif characteristics of GRFs. We found some structural features among clades or subclades according to the alignment and motif results. As shown in Figure 2, a total of 10 conserved motifs of GRFs were found using the MEME online software. Each subgroup has two to nine conserved motifs, and the motif composition is similar within subgroups. Some motifs appeared in only certain specific subgroups, which possibly contribute to functional diversity. For example, motif 6 is unique to subgroup I, and motif 5 is unique to subgroup II, while motifs 3, 4, 7, 9, and 10 are concentrated in group III. Motifs 3, 4, and 9 are specific to group IV. Motifs 3 and 4 are specific to group VI. We also investigated GRF gene structures to further characterize their evolutionary trajectory. PpGRFs contained conserved QLQ and WRC domains in their Nterminal regions, which was also found in the other species. However, group V has two WRC domains. These results indicate that they contain different numbers of exons, varying from 3 to 12. The gene structures were similar or the same in each subgroup. Overall, the gene structures and motif characteristics strongly supported the phylogenetic relationships of GRFs between peach and other three Rosaceae species.

Chromosomal localization and collinearity analyses of PpGRF genes
The distribution of PpGRF genes is random as seen in Figure  3. PpGRF2, 3 and 4 were on distributed on chromosomes 2, PpGRF6, 7 on chromosomes 5, PpGRF9,10 on chromosomes 7, PpGRF1 on chromosomes 1, PpGRF8 on chromosomes 6 ( Figure 3). Chromosome 4, 8 do not have a PpGRF gene. In peach, some chromosomes do not have a GRF gene, such as chromosome 4 and chromosome 8. This phenomenon was also found in Pyrus communis, Fragaria vesca, Malus domestica, which suggested that the GRF genes were distributed unevenly in chromosomes.
Subsequently, to further infer the phylogenetic relationship between peach and other Rosaceae species, we analyzed the collinear relationships between peach and three Rosaceae species (Figure 3). We found that 9 FvGRF genes were collinear with PpGRF genes, followed by MdGRFs (14) and PcGRFs (10). Moreover, the collinear relationship in the peach genome was also examined to elucidate the evolution and origin dynamics of PpGRF genes. Three pairs of PpGRF genes (PpGRF01/08, PpGRF03/07, and PpbZIP05/10) have a collinear relationship and were generated by whole-genome duplication (WGD) or segmental duplication (Figure 3 and Additional File 2: Supplementary Table S2). Furthermore, to determine the selection constraints on the duplicated PpGRF genes, we calculated the non-synonymous/synonymous substitution ratio (Ka/Ks) of each pair of duplicated genes. The Ka/Ks ratios of most PpGRF pairs were less than 1, indicating that these PpGRFs had undergone purifying selection processes (Additional File 2: Supplementary Table S2).

Analysis of cis-elements in the promoter sequences of PpGRF genes
To explore the involvement of the PpGRF genes in light responsiveness, hormone signaling pathways, plant growth and development, their promoter sequences were analyzed using PlantCARE software. All cis-elements identified in the promoter regions of the PpGRF gene family members are shown in Figure 4. There are at least six commonly occurring light-responsive elements (LREs): AE-box, ATCT-motif, G-box, GT1-motif, GATA-motif and I-box, which have been demonstrated to be essential for the regulation of light-mediated transcriptional activity. The results indicated that all 10 PpGRF promoter regions contained two or more LREs among which G-box is the most abundant element (Additional File 3: Supplementary Table S3). In this region, we identified many other important potentials-elements, such as GARE motif and P-box (GA responsive element), a CAT-box (cis-acting regulatory element related to meristem expression), AuxRR-core, and TGA-element (auxin-responsive element), an RY element (seed-specific regulation) and a TC-rich repeats (defense and stress responsiveness). GARE, P-box, CATbox, AuxRR-core, TGA-element, RY element and TCrich repeats were distributed within 3, 3, 7, 2, 3, 1 and 5 PpGRF promoter regions, respectively (Additional File 4: Supplementary Fig S1).

Expression pattern of GRF genes in peach
Although we identified GRF genes in the peach genome, the functions of these genes remain largely unknown. In general, GRFs were highly expressed in growing tissues in Arabidopsis. Next, we investigated the expression patterns of PpGRFs in root, stem, leaves, the shoot tip and hypocotyl. Except for PpGRF8 (data not shown), these genes were differentially expressed in different tissues ( Figure 5). Most of these genes were more expressed, especially in the shoot tip and root ( Figure 5). For example, the expression levels of PpGRF 2, 3, 4, 5, 6, 7, 9, and 10 were the highest in the shoot tip. However, PpGRF1 showed relatively strong expression in the root. Of the analyzed genes, PpGRF3, 9and 10 were more preferentially expressed in the shoot tip, more than 20 times greater than that in other tissues. To further investigate the potential functions of PpGRFs, we also examined the effects of GA 3 and UVB on the expression of the PpGRF gene ( Figure 6). More than half of the total genes were up-regulated in response to abiotic stresses, with PpGRF1, 4, 5, 6, 7, and 10 being up-regulated after  GA 3 treatment; PpGRF2, 3, 4, 5, 6, and 7 were up-regulated after UVB treatment. On the other hand, several PpGRF genes (including PpGRF2, 6 and 7) were suddenly downregulated at 9 h after GA 3 treatment ( Figure 6). Also, PpGRF5, 9, and 10 were suddenly down-regulated at 6 h after UVB treatment. Based on the expression patterns of PpGRFs, we propose that they likely have multiple functions in regulating growth and responsiveness to abiotic stresses.

UVB effect on peach leaf growth
The peach trees were treated with 1.44 Kj m −2 d −1 UVB radiation during the whole growing period. It can be seen from Figure 7

Hormone profiles after UVB exposure in peach
Previous studies have demonstrated that GRFs play important roles in regulating leaf growth (Kim et al. 2003;Horiguchi et al. 2005;Gonzalez et al. 2012;Debernardi et al. 2014;Lee et al. 2015). The shoot tip is the growth point of the new shoot, and its hormone content determines the growth of the new shoot. At this facility, the new shoots of peach grow vigorously from June to September, so in this experiment, the shoot tip tissue was taken on August 7 and frozen for hormone content analyses. The hormone content analyses are shown in Figure 8: Under UVB treatment, the content of IAA varies greatly, and the IAA content of the shoot tip tissue of the CK group (1.978917 pmol/g) is 4.9 times that of the UVB treatment group (0.40339 pmol/g). At the same time, the GA 3 content was increased in the UVB treatment group (0.59663pmol/g) relative to the CK group (0.22767 pmol/g). In our hormone levels analysis, IAA content showed a significant difference (P < 0.01). The GA 3 content also showed a significant difference (P < 0.05). Correlation analysis showed that there was significant correlation among the IAA and GA 3 (CK treatment: r = 0.9979; UVB treatment: r = −0.9979). The content of abscisic acid (ABA), jasmonic acid (JA) and trans-Zeatin-riboside (TZR) were not significantly different when comparing the CK group. Therefore, we speculate that the UVB-mediated growth inhibition of the new shoot is closely related to the IAA and GA 3 contents in the shoot tip.

Prediction of the miR396 target site and how its expression changes following UVB exposure
Differently to Arabidopsis, all GRFs in peach have a sequence that is partially complementary to miR396, which is located at the WRC domain. To further explore the evolutionary relationship among GRFs, we analyzed the WRC structures of all GRFs in P. persica. According to Figure 9, we found that all PpGRFs contain this part, which is complementary to the miR396, share a bulge at position 7 (counting from the 5 end of the miRNA). From this point of view, it is to be expected that GRFs in P. persica are highly conserved. We have identified three differentially expressed miRNAs predicted to be responsive to low-dose UVB in P. persica.
To identify the expression profiles of miRNA396, we utilized our transcriptome data derived from Illumina RNA-Seq reads generated and analyzed by . We found that, in peach, miR396 was down-regulated after UVB irradiation, but this did not reach statistical significance. Also, we isolated another mirR319 from P. persica, which regulates TCP4, induces miR396, and represses GRF activity, which was also non-significantly down-regulated following UVB irradiation (Additional File 5: Supplementary  Table S4).

Discussion
GRFs are plant-specific TFs that play important roles in plant growth and development, especially in regulating organ size Omidbakhshfard et al. 2015;Li et al. 2016). The GRF gene family in model plants (Arabidopsis, rice, soybean, and poplar) has been identified and described in detail, but little is known about the GRFs in P. persica (Kim et al. 2003;Choi et al. 2004;Cao et al. 2016;Chen et al. 2019). To better understand the characteristics and functions of GRFs in P. persica, we studied PpGRF family members by bioinformatics assay with the other three Rosaceae plants. Meanwhile, by performing RT-qPCR analysis, plant morphology analysis and hormone content determination to analyze the function of PpGRF. In the study, we identified 46 GRF members in four representative Rosaceae plants, 10 members of the GRF family in Prunus persica, Pyrus communis, Fragaria vesca, respectively and 16 members of the GRF family in Malus domestica (Table 1) and found that the number of GRFs in apples is relatively large, perhaps due to a genome-wide replication event in the apple genome (Velasco et al. 2010;Cao et al. 2013). Phylogenetic analysis showed that the four species of Rosaceae GRF could be divided into six groups (I-VI) (Figure 1), which is consistent with the classification of Arabidopsis, soybean, and rice GRFs (Kim et al. 2003;Choi et al. 2004;Chen et al. 2019). The homologous pairs of GRF proteins in the four Rosaceae plants are more common. According to the topology of the phylogenetic tree, we found that some ancestral GRFs existed before the evolutionary divergence of the Rosaceae. By analyzing the GRF domains, we found that all of the speculated GRFs contain QLQ and WRC domains, and the VI subfamily contains two WRC domains. Similar results also exist in the Arabidopsis genome but not in the rice genome, suggesting that a major expansion occurred after the eudicots diverged from the monocots and before the separation of Rosaceae (Additional File 6, Supplementary Table S5). The structures of GRFs belonging to different subgroups in intron and exon structures showed low homology. GRFs of the same subfamily have the same or similar genetic structure, especially the number and length of exons. Also, based on our MEME analysis, each subfamily contains distinct motif arrangements. The differences in these characteristics of subfamilies suggest that GRF members are functionally diverse, which also supports the classifications used in this study (Kim et al. 2003;Choi et al. 2004;Cao et al. 2016;Chen et al. 2019). The expansion of the GRF family occurs mainly through gene duplication, especially large-scale duplication (WGD or segment duplication) Chen et al. 2019). In this study, three pairs of PpGRF genes were distributed in duplication blocks, which means large-scale, repeats that promote the amplification of the GRF gene family in peach. We also found that the Ka/Ks ratios of all the PpGRFs gene pairs were less than 1, indicating that they have experienced strong purifying selection.
We further explored the expression patterns of PpGRFs. First, RT-qPCR was used to detect the expression of PpGRFs in various tissues. There are some important parallels in the spatial expression patterns of PpGRF genes and those from other plant species. Several studies have reported that, with the aging of organs, the transcription levels of GRFs decrease. The expression levels of PpGRF2,3,4,5,6,7,9,and 10 in the shoot tip meristem were higher than those in mature tissues, suggesting that these genes play important roles in elongation of the new shoot. PpGRF1 was more highly expressed in the root than in the other tissues, indicating that PpGRF1 is also important in root growth. These results showed that most of PpGRF genes are expressed in actively growing and developing tissues, in agreement with previous studies (Kim et al. 2003;Lee et al. 2015;Chen et al. 2019), indicating that PpGRFs play critical roles in plant growth and development. A well-documented effect of UVB exposure on plants is the reduction of biomass, and miR396/GRFs exert their effect on organ growth at least in part by regulating the expression of genes involved in GA metabolism Hewezi et al. 2012;Vercruyssen et al. 2015). Our study demonstrates that PpGRFs can widely respond to UVB and GA 3 treatments at transcriptional levels, with almost half of the members show induction and a few show repression.
The excessive vigor of fruit trees will lead to new shoots, which is not conducive to the storage of assimilates, flower bud differentiation, or the balance between vegetative and reproductive growth (Chen 2009). This is one of the critical  prevention and control goals in fruit tree cultivation management technology. Here we found that during the growth period of the new shoot, UVB supplementary light treatment can significantly inhibit the growth of the new shoot but has little effect on the changes in the thickness of the new shoot. Perhaps UVB inhibits the proliferation of cells in developing leaves. The shoot tip is the growth point of the extension of the new shoot, and its hormone content determines the growth of the new shoot. GA stimulates cellular expansion and induces the expression of the PpGRFs, these genes appear to function in cellular expansion, leading to an increase in leaf size and so on. IAA promotes cell elongation and required for primary elongation in peach. UVB treatment can significantly reduce the concentration of IAA and increase the concentration of GA 3 in the shoot tips. The concentration of IAA had significant negative correlation with the concentration of GA 3 under UVB radiation. Therefore, we infer that the levels of IAA and GA 3 in the shoot tips inhibit the growth of the new shoot to reach equilibrium and thus mediate tree vigor. Additional studies are needed to explore this phenomen in more detail. Studies in maize indicate that the inhibitory effect of UVBs on growth is mediated by GRF, some of which might be regulated by miR396, thereby modifying GA levels in the leaf growth zone (Casati and Walbot 2003;Casadevall et al. 2013;Fina et al. 2017). To further verify the above hypothesis, we tested the expression of miR396. Under UVB treatment, miR396 maintained a low expression level in the function leaves and was slightly down-regulated. Unlike in Arabidopsis, miR396 was expressed at low levels in the leaf primordia and then steadily accumulated with the development of the leaf in Arabidopsis (Rodriguez et al. 2010). At the same time, mirR319 expression was also less expressed (although not reaching statistical significance), which maybe result from detecting the miRNA in the mature leaves. As previous reports showed that most GRF genes are strongly expressed in actively growing and developing tissues, such as shoot tip, root and flower bud, but weakly in mature stem and leaf tissues (Kim et al. 2003;Omidbakhshfard et al. 2015).
Notably, PpGRF1, 2, 3, 5, and 9 showed a response peak after 9 h of GA 3 treatment, inconsisitent with the response peak after UVB treatment. The analysis of cis-elements in the PpGRF promoter sequences found that some members contained multiple types of light-responsive elements and GA-responsive elements, implying the involvement of PpGRFs in different light responses and hormone signaling pathways. This shows that there are differences in how GRFs participate in these two signaling pathways. There is much evidence indicating the role of GRFs in various biological processes through interacting with different clients, such as hormone pathways, stomatal behavior, reactive oxygen species balance, and ion transport He et al. 2017;Kaundal et al. 2017;Liu et al. 2017;Yang et al. 2017). In Arabidopsis and maize, UVB radiation inhibits leaf growth by decreasing the expression of GRFs (Casati and Walbot 2003;Casadevall et al. 2013;Fina et al. 2017). This ultimately leads to a reduced expression of GA biosynthetic genes and increased levels of catabolic transcripts (Casati and Walbot 2003;Rodriguez et al. 2010;Casadevall et al. 2013;Fina et al. 2017). In addition, the PpGRF promoter regions containing multiple types of cis-elements, which is useful for exploring the upstream genes of the PpGRFs involved in regulating growth and responsiveness to abiotic stresses. In analysis of cis-elements in the PpGRF promoter sequences, we found that all PpGRFs contained light-responsive elements, at the same time, UVB treatment led to changes of all these genes. In Arabidopsis, all AtGRFs contained light-responsive elements according to the analysis of cis-elements (Additional file 7 Supplementary Fig  S2). Casadevall et al. (2013) reported that UVB radiation modulates the expression of GRF1, 2 and 3. Future research should explore the veiled functions of GRFs with regard to UVB radiation in plants. Interestingly, we found that PpGRF2, 3, 4, and 10 have GA-responsive elements, and the expression analysis of these genes showed response to GA treatment. GA-responsive elements have not been detected in PpGRF1, 5, 6, 7, 8, and 9, but these genes still responded to GA treatment. We propose that this result maybe due to some new GA-responsive elements has not been discovered in these genes or other upstream TFs respond to GA treatment regulate the expression of these genes (PpGRF1,5,6,7,8,and 9). Overall, further study of the function of PpGRFs can reveal the behaviors of PpGRFs during regulating tree growth and development, providing useful information for further identification of the biological functions of PpGRFs.

Conclusions
In this study, we identified 10 PpGRFs from the peach genome and investigated their phylogenetic classification, protein motifs, and gene structures. By further analysis, PpGRFs belonging to the same subgroup had similar gene structures and motif compositions. Chromosomal localization and collinearity analyses suggested that large-scale duplication contributed to the expansion of the PpGRF family. Physiological experiments and expression analyses have revealed the involvement of PpGRFs under UVB and GA 3 treatment and identified some important candidates for regulating the elongation of the new shoot. Additional studies are needed to explore this hypothesis in greater detail. This systematic study gave a basic understanding of GRFmediated signal cascades regulating peach tree development and abiotic stress response. These results will provide a foundation for future research on the GRFs in peach. Figure 9. Prediction of the miR396 target site. Scheme representing the GRF genes. The interaction of the ten GRFs from peach with miR396 is shown. QLQ and WRC indicate the conserved domains that define the GRF family.

Acknowledgement
We thank all the colleagues in our laboratory for providing useful discussions and technical assistance. We are very grateful to TopEdit (www.topeditsci.com) for English language editing of this manuscript. We sincerely thank Dr. Chengjie Chen share of TBtool software package.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Notes on contributor
Li Liu, research assistant, research interest: cultivation physiology and molecular biology of fruit tree, formation and regulation of fruit quality, mineral nutrition of fruit tree, nitrogen nutrition, calcium signaling.
Xiu-jie Li, associate professor, research interest:cultivation physiology of fruit tree, mineral nutrition of fruit tree, formation and regulation of fruit quality.
Bo Li, professor, research interest: cultivation physiology and molecular biology of fruit tree, formation and regulation of fruit quality, mineral nutrition of fruit tree, nitrogen nutrition, calcium signaling.
Ming-yue Sun, teaching assistant, research interest: cultivation physiology and molecular biology of fruit tree, bioinformatics of fruit tree, photomorphogenesis.
Shao-xuan Li, associate professor, research interest: cultivation physiology and molecular biology of fruit tree, bioinformatics of fruit tree, photomorphogenesis.

Availability of data and materials
The four Rosaceae species genome datasets (Prunus persica, Pyrus communis, Fragaria vesca and Malus domestica) used for identifying the GRF genes in this study were downloaded from the GDR (http: //www. rosaceae.org). The other species sequences (Arabidopsis thaliana, Glycine max, Oryza sativa and Populus trichocarpa) were downloaded from the plantTFDB (http: // planttfdb. cbi. pku. edu. cn). The RNA-seq data were obtained from our laboratory. The datasets supporting the conclusions of this article were included with in the article and its Supplementary files.

Ethics approval and consent to participate
All the peach materials used in the experiment were collected from Facilities Fruit Tree Laboratory, Shandong Agricultural University. These plant materials are widely used all over the world. and no permits are required for the collection of plant samples. This article did not contain any studies with human participants or animals performed by any of the authors, and did not involve any endangered or protected species.