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

ABSTRACT

Growth-regulating factors (GRFs) are plant-specific 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 (GA₃) 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 GA₃ signals. The evolutionary patterns and expression profiles of PpGRFs detected in this study increase understanding of the important roles in peach.

KEYWORDS:

• Peach
• growth-regulating factor
• gene family
• UVB
• GA₃
• new shoot

1. 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 (GRF-interacting 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₃ 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₃ 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 narrow-leaf 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 (Liu et al. 2009; 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 (Li et al. 2017). 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 (Chen et al. 2016; Sun et al. 2016; Zhang et al. 2018, 2012; Li et al. 2019a, 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.

2. Methods

2.1. 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 (Li et al. 2017), the optimal UVB treatments were carried out with 1.44 Kj m⁻²d⁻¹ 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 (Li et al. 2017), hormone treatments were conducted by spraying young leaves with 100 µM GA₃ . 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.

2.2. Sequence collection and identification

The four Rosaceae genome sequences were downloaded from the Genome Database for Rosaceae (GDR) (http://www.rosaceae.org/). To identify putative GRF genes from peach, European pear (Pyrus communis), strawberry (Fragaria vesca), and apple (Malus d omestica), the HMM profiles of the GRF domains (PF08879 and PF08880) were downloaded from Pfam (http://pfam.xfam.org/) were used to search the database using HMMER3.0 (http://hmmer.janelia.org/) with the default E-value. SMART (http://smart.embl-heidelberg.de/) and the NCBI CDD (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi) were also used to determine the PF08879 and PF08880 domains of GRF, and sequences lacking similar integrity were removed. Finally, their molecular weights and isoelectric points were predicted using the ExPASy proteomics server database (http://expasy.org/).

2.3. 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 (Chen et al. 2018).

2.4. 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 (Chen et al. 2018; Xie et al. 2018).

2.5. 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).

2.6. 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) (Chen et al. 2016). 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).

2.7. 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.

2.8. 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₃), 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₂Cl₂) 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.

The standard fluids of diiferent hormones (IAA, ABA, GA₃, JA, TZR) were prepared for different concentrations (0.2, 0.5, 1, 2, 5, 10, 20, 50, 200 ng·mL⁻¹) using respective standard dissloved in methanol (0.1% formic acid).

The liquid phase conditions are set as follows: chromatographic column:Agilent ZORBAX SB-C18, reversed-phase chromatographic column, 2.1×150, 3.5 μm; column temperature: 30°C. injection volume:10 μL. Mobile phase A: B = water /0.1% formic acid: methanol. The mass spectrum conditions are set as follows: curtain gas: 15 psi; spray voltage: +4000 V, −3500 V; atomizing gas pressure: 65 psi; auxiliary pressure: 70 psi; atomization temperature: 350°C.

2.9. Target prediction of miR396 and expression analysis following UVB exposure

Target prediction of miR396 followed the rules referring to Li et al. (2017). 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.

3. Results

3.1. Identification of GRFs in peach

The GRF Hidden Markov Model (HMM) configuration file (PF08879 and PF08880) was used to identify GRF members. Furthermore, the SMART tool (http://smart.embl-heidelberg.de/) was used to verify the candidate genes’ presence of two characteristic conserved domains (QLQ and WRC). After manually removing the incomplete sequences from Rosaceae database, we identified 10 PpGRF, 10PcGRF, 10 FvGRF and 16MdGRF protein sequences (Table 1). The PpGRF were unevenly distributed across the chromosomes, similar to other three Rosaceae species (Pyrus communis, Fragaria vesca, Malus domestica). The 10 predicted GRF proteins ranged from 331(PpGRF4) to 612 (PpGRF8) amino acid residues in length, and their relative molecular mass varied from 36.23 (PpGRF4) to 66.01(PpGRF8) kDa. The isoelectric point (pI) values range from 6.88 (PpGRF8) to 9.47 (PpGRF2), with 1 member having pI value <7, and 9 members having pI values >7 (Table 1).

Table 1. Annotation of GRF genes in Prunus persica and three Rosaceae species (Pyrus communis, Fragaria vesca and Malus domestica).

Annotation of GRF genes in Prunus persica,Pyrus communis, Fragaria vesca and Malus domestica
Species name	Gene name	Gene ID	Chr locous	5′ End	3′ End	Sequence Length(aa)	Chain	Molecular weight	Isoelectric point
Prunus persica	PpGRF1	Prupe.1G419000	Pp01	36302980	36306676	591	+	63720.55	8.8
	PpGRF2	Prupe.2G148100	Pp02	20403542	20406569	444	+	47925.72	9.47
	PpGRF3	Prupe.2G317500	Pp02	29721181	29724355	546	+	59301.99	7.16
	PpGRF4	Prupe.2G329200	Pp02	30235010	30239972	331	–	36227.06	8.23
	PpGRF5	Prupe.3G041800	Pp03	2977871	2980492	526	+	59104.78	9.03
	PpGRF6	Prupe.5G021400	Pp05	2337057	2342732	336	+	36948.34	8.92
	PpGRF7	Prupe.5G131600	Pp05	12672215	12674510	557	+	61570.07	8.09
	PpGRF8	Prupe.6G217000	Pp06	22407811	22411275	612	+	66013.13	6.88
	PpGRF9	Prupe.7G115500	Pp07	14076370	14079828	390	–	42539.21	8.7
	PpGRF10	Prupe.7G201000	Pp07	18750960	18753565	379	+	42814.15	8.58
Malus domestica	MdGRF1	MD02G1088800	Chr02	6989543	6991135	229	–	25790.98	9.65
	MdGRF2	MD02G1189500	Chr02	17428479	17432083	593	–	64281.31	7.65
	MdGRF3	MD02G1201000	Chr02	19847667	19850955	467	–	49847.81	9.53
	MdGRF4	MD04G1215900	Chr04	29871031	29873763	485	–	52187.74	8.97
	MdGRF5	MD06G1014700	Chr06	1876356	1881403	306	–	33821.7	8.57
	MdGRF6	MD06G1121000	Chr06	26143211	26145619	481	+	52738.11	8.28
	MdGRF7	MD07G1297800	Chr07	35757203	35760130	519	+	55964.31	6.87
	MdGRF8	MD09G1183300	Chr09	15790752	15792881	325	+	36261.78	8.78
	MdGRF9	MD12G1020500	Chr12	2035411	2038156	394	+	42772.49	7.31
	MdGRF10	MD14G1017600	Chr14	1652199	1655992	397	+	43120.8	8.83
	MdGRF11	MD14G1137600	Chr14	22492159	22495911	477	+	52294.42	8.58
	MdGRF12	MD15G1065000	Chr15	4531603	4535131	579	+	62492.89	9.14
	MdGRF13	MD15G1216000	Chr15	17431154	17433559	369	–	41693.06	8.96
	MdGRF14	MD15G1300900	Chr15	28847105	28850233	600	–	64887.1	7.66
	MdGRF15	MD00G1001000	Chr00	238417	244764	340	–	37988.66	8.63
	MdGRF16	MD00G1142400	Chr00	30993049	30995600	365	–	41389.59	8.83
Pyrus communis	PcGRF1	pycom02g07120	Chr2	4973472	4975168	365	–	41196.57	8.69
	PcGRF2	pycom02g16430	Chr2	13577598	13580385	465	–	49430.38	9.61
	PcGRF3	pycom06g01200	Chr6	1236492	1240784	329	–	36010.15	8.72
	PcGRF4	pycom07g27020	Chr7	26868100	26870390	528	+	56812.23	7.17
	PcGRF5	pycom08g06450	Chr8	5037777	5040611	570	+	61254.46	7.75
	PcGRF6	pycom09g10100	Chr9	8213622	8215108	329	+	36481.02	8.43
	PcGRF7	pycom14g01600	Chr14	1191231	1193776	395	+	43060.74	9.03
	PcGRF8	pycom15g19160	Chr15	13652357	13654338	378	–	42200.61	9.68
	PcGRF9	pycom15g26320	Chr15	21175651	21178440	361	–	64675.67	6.24
	PcGRF10	pycom16g25430	Chr16	27332996	27338421	187	–	20698.61	9.9
Fragaria vesca	FvGRF1	FvH4_1g08440	Fvb1	4473841	4475847	369	–	41542.65	8.4
	FvGRF2	FvH4_1g18220	Fvb1	10597643	10600896	592	–	64067.9	6.09
	FvGRF3	FvH4_2g32670	Fvb2	24704625	24708301	580	+	62689.86	9.09
	FvGRF4	FvH4_5g10750	Fvb5	6093251	6095578	474	–	52627.9	8.99
	FvGRF5	FvH4_5g23900	Fvb5	15201347	15204913	330	–	36795.22	9.12
	FvGRF6	FvH4_6g18130	Fvb6	11941323	11943517	367	+	40258.77	8.75
	FvGRF7	FvH4_6g34730	Fvb6	27401421	27403313	328	–	37536.34	8.83
	FvGRF8	FvH4_7g12130	Fvb7	11046464	11049139	445	+	48302.07	9.25
	FvGRF9	FvH4_7g32780	Fvb7	23498439	23500857	543	+	58521.69	8.47
	FvGRF10	FvH4_7g34250	Fvb7	24166270	24173059	843	–	94268.61	7.58

3.2. 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.

Figure 1. Phylogenetic analysis of GRF from Arabidopsis, oryza, populus, glycine max, peach, european pear, strawberry and apple. Using the Neighbor-Joining method to construct the phylogenetic tree. Six subfamilies are shown in different colors.

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 N-terminal 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.

Figure 2. Phylogenetic relationships, gene structures and motif structure in GRF genes from Prunus persica and three Rosaceae plants (Pyrus communis, Fragaria vesca and Malus domestica). (a) The phylogenetic tree was constructed based on the full-length sequences of GRF proteins using MEGA7.0. (b) Schematic representation of the 10 conserved motifs in GRF proteins. Members of same group are arranged sequentially according to phylogenetic classification. The conserved motifs labeled with different colors boxes. The length of the protein can be estimated using the scale at the bottom. (c) Exon-intron structures of the GRF genes. Green boxes indicate untranslated 5′- and 3′-regions; yellow boxes indicate exons; and black lines indicate introns. The QLQ and WRC domains are highlighted by blue and red boxes. The number indicates the phases of the corresponding introns. The length of the nucleic acid can be estimated using the scale at the bottom.

3.3. 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.

Figure 3. Gene location and collinearity analysis of the GRF genes between Prunus persica and three Rosaceae plants (Pyrus communis, Fragaria vesca and Malus domestica). Gray bars in the background indicate the collinear blocks within the peach and other plant genomes. Different genomes are indicated by group-specific color. Syntenic gene pairs are highlighted by lines in different colors.

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).

3.4. 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, CAT-box, AuxRR-core, TGA-element, RY element and TC-rich repeats were distributed within 3, 3, 7, 2, 3, 1 and 5 PpGRF promoter regions, respectively (Additional File 4: Supplementary Fig S1).

Figure 4. Schematic model of six cis-elements in the promoter sequences of PpGRF genes. Members of GRF are arranged sequentially according to phylogenetic classification. The length of the nucleic acid can be estimated using the scale at the bottom.

3.5. 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₃ 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₃ 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 down-regulated at 9 h after GA₃ 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.

Figure 5. Expression profiles of PpGRF genes in sixs peach tissues. Root (R), stem (S), leaves (L), shoot tip (ST) and hypocotyl (H). The x-axis represents different peach organs. The y-axis represents the relative expression levels of PpGRF genes compared with actin. The expressions of PpGRF genes in root were normalized to 1. Error bars represent standard deviations of the means of three independent replicates. Different letter denote statistical significance according to one-way ANOVA and Duncan’s test (P < 0.05).

Figure 6. Relative expression analysis of the PpGRFs under UVB and GA₃ treatment. The x-axis represents different time. The y-axis represents the relative expression levels of PpGRF genes compared with actin. Data presented are the means ± SD of three replicates. Different letters denote statistical significance according to one-way ANOVA and Duncan’s test (P < 0.05).

3.6. UVB effect on peach leaf growth

The peach trees were treated with 1.44 Kj m⁻² d⁻¹ UVB radiation during the whole growing period. It can be seen from Figure 7 that during the growth period of the new shoot, UVB supplementary light treatment can inhibit the elongation of the new shoot. After harvesting for 120 days (September 7), the new shoot of the UVB group stopped lengthening, while the CK group continued to increase slowly. The average new shoot length (114.08 cm) was about 1.5 times that of the UVB group (69.45 cm). The thickness of the new shoot in the CK group increased by 0.3 cm, while the thickness of the new shoot in the treatment group increased by 0.29 cm; the difference was not statistically significant. Based on this finding, this dose of UVB treatment had little effect on the thickness of the new shoot.

Figure 7. Effects of UVB treatment on the growth of peach new shoots. a. new shoot length, b. new shoot thickness. The data indicate variation after harvest from 7 June 2019–7 October 2019. Values are means (n = 3) and error bars represent the standard deviation.

3.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₃ 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₃ content also showed a significant difference (P < 0.05). Correlation analysis showed that there was significant correlation among the IAA and GA₃ (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₃ contents in the shoot tip.

Figure 8. Contents of different hormones in stem apex. a: IAA; b: ABA; c: GA₃; d: JA; e: TZR. Data is presented as the means ± SD, n = 3. Different asterisks above bars denote statistical significance according to one-way ANOVA and Duncan’s test.

3.8. 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 Li et al. (2017). 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).

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.

4. Discussion

GRFs are plant-specific TFs that play important roles in plant growth and development, especially in regulating organ size (Debernardi et al. 2014; 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) (Cao et al. 2016; 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 (Liu et al. 2009; Hewezi et al. 2012; Vercruyssen et al. 2015). Our study demonstrates that PpGRFs can widely respond to UVB and GA₃ 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₃ in the shoot tips. The concentration of IAA had significant negative correlation with the concentration of GA₃ under UVB radiation. Therefore, we infer that the levels of IAA and GA₃ 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₃ 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 (Chen et al. 2017; 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.

5. 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₃ 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 GRF-mediated signal cascades regulating peach tree development and abiotic stress response. These results will provide a foundation for future research on the GRFs in peach.

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.

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).

Correction Statement

This article has been corrected with minor changes. These changes do not impact the academic content of the article.

Additional information

Funding

This work was supported by Key Laboratory of Biology and Genetic Improvement of Horticultural Crops (Germplasm Resources Utilization), Ministry of Agriculture, P.R.China: [grant number NYZS202101]; Shandong fruit industry technology system: [grant number SDAIT-06-09].

Notes on contributors

Li Liu

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

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

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

Ming-yue Sun, teaching assistant, research interest: cultivation physiology and molecular biology of fruit tree, bioinformatics of fruit tree, photomorphogenesis.

Shao-xuan Li

Shao-xuan Li, associate professor, research interest: cultivation physiology and molecular biology of fruit tree, bioinformatics of fruit tree, photomorphogenesis.