Identification, Classification, and Expression Analysis of GRAS Gene Family in Malus domestica

GRAS genes encode plant-specific transcription factors that play important roles in plant growth and development. However, little is known about the GRAS gene family in apple. In this study, 127 GRAS genes were identified in the apple (Malus domestica Borkh.) genome and named MdGRAS1 to MdGRAS127 according to their chromosomal locations. The chemical characteristics, gene structures and evolutionary relationships of the MdGRAS genes were investigated. The 127 MdGRAS genes could be grouped into eight subfamilies based on their structural features and phylogenetic relationships. Further analysis of gene structures, segmental and tandem duplication, gene phylogeny and tissue-specific expression with ArrayExpress database indicated their diversification in quantity, structure and function. We further examined the expression pattern of MdGRAS genes during apple flower induction with transcriptome sequencing. Eight higher MdGRAS (MdGRAS6, 26, 28, 44, 53, 64, 107, and 122) genes were surfaced. Further quantitative reverse transcription PCR indicated that the candidate eight genes showed distinct expression patterns among different tissues (leaves, stems, flowers, buds, and fruits). The transcription levels of eight genes were also investigated with various flowering related treatments (GA3, 6-BA, and sucrose) and different flowering varieties (Yanfu No. 6 and Nagafu No. 2). They all were affected by flowering-related circumstance and showed different expression level. Changes in response to these hormone or sugar related treatments indicated their potential involvement during apple flower induction. Taken together, our results provide rich resources for studying GRAS genes and their potential clues in genetic improvement of apple flowering, which enriches biological theories of GRAS genes in apple and their involvement in flower induction of fruit trees.


INTRODUCTION
Transcription factors function as trans-acting factors, combining with specific cis-elements in eukaryotic gene promoter regions to regulate plant growth and development.
Recent studies have revealed that the N-termini of GRAS proteins play important roles in their specific functions. The N-termini contain intrinsically disordered regions (IDRs); thus, GRAS proteins are also known as intrinsically disordered proteins (IDPs) (Sun et al., 2011). Several studies have focused on the structures and functions of the IDRs in the N-termini of GRAS proteins (Sun et al., 2010(Sun et al., , 2011(Sun et al., , 2012. Different subfamilies have different types and numbers of molecular recognition features in their N-termini (Sun et al., 2010(Sun et al., , 2011(Sun et al., , 2012. For instance, the DELLA subfamily has three molecular recognition features in its IDR at the N terminus: DELLA, VHYNP, and L(R/K) XI. Protein motifs at the C-terminus include the LHR I, HIID, LHRII, PFYRE, and SAW motifs (Pysh et al., 1999).
The functions of GRAS family genes have been studied in recent years; there is a large body of evidence that the GRAS genes family play important roles in various biological processes in many different plant species. They are involved in diverse processes related to root meristem development, axillary bud outgrowth, phytochrome signaling pathways, nodular signal transduction, and stress adaptation, especially in the gibberellins (GA)-signaling pathway (Ait-ali et al., 2003;Greb et al., 2003;Heckmann et al., 2006;Torres-Galea et al., 2006;Koizumi et al., 2012;Park et al., 2013). DELLA proteins are negative regulators of the GA signal, and have been shown to repress growth in many different plant species including Arabidopsis (AtGAI, AtRGA, AtRGL1, AtRGL2, and AtRGL3), O. sativa (OsSLR1), Zea mays (ZmD8), Malus domestica (MdRGL2a), and Vitis vinifera (VvGAI) (Foster et al., 2007;Vandenbussche et al., 2007;Aleman et al., 2008;Dai and Xue, 2010). Various studies have shown that the GA-GID1-DELLA complex and ubiquitin ligase E3 play important roles in GA signaling, and that SCR functions in the development of roots and aboveground organs in Arabidopsis (Zhang et al., 2011;Koizumi et al., 2012). PhHAM (Petunia hybrid), a petunia hairy meristem gene, was shown to be necessary to maintain the apical meristem in petunia, and the petunia recessive mutant phham (ham-B4281) produced fewer flowers than wild type (Stuurman et al., 2002). SCL13 was shown to positively regulate phyB, and the Abbreviations: DAFB, days after full blossom; ES, early stage; MS, middle stage; LS, later stage; GA 3 , gibberellic acid; 6-BA, 6-benzylaminopurine. scl13 mutant showed impaired sensitivity in Arabidopsis (Torres-Galea et al., 2006). Several studies have shown that some GRAS genes are up-regulated under stress conditions in Arabidopsis, rice and tobacco (Day et al., 2003;Torres-Galea et al., 2006;Czikkel and Maxwell, 2007). Among all the biological processes, their flowering related characterizations and functions have been researched in some model plants. On the one hand, several GRAS genes are associated with the flowering phenotype. For example, in tomato, Ls (Later Suppressor), defined as another GRAS gene, its mutant ls showed lower numbers of inflorescence (Schumacher et al., 1999;Bolle, 2004). In Arabidopsis, the LS subfamily gene AtLAS and the SHR subfamily gene AtSHR were involved in shoot apical meristem. The Ls subfamily gene AtSCL26 and the HAM subfamily genes AtSCL6 and AtSCL27 were highly expressed in flowers (Bolle, 2004); the rga and gai mutant also showed earlier flowering in Arabidopsis (Galvao et al., 2012). On the other hand, GRAS genes can integrate or regulate flowering related genes such as SPL (SQUAMOSA PROMOTER BINDING PROTEIN-LIKE), LFY (LEAFY), and AP1 (APETALA1) as well as FT (FLOWERING LOCUS T). They all acted positive roles and intergrated various signals in flower induction, of which FT and SPL affected flowering time; while LFY and AP1 belonged to flower meristem identify genes. Several GRAS genes can be recruited and integrated with SPL proteins to regulate LFY and AP1 expression in Arabidopsis (Yamaguchi et al., 2014). Additionally, the PAT and SCL subfamilies of GRAS genes were both reported involved in photosignal in many plants, it influenced the expression of phyA (phytochrome A), a phytochrome A signal transduction gene, which caused to the change of the circadian clock and expression of FT and other flowering genes (Bolle et al., 2000;Torres-Galea et al., 2006;Li et al., 2014). However, little is known about the GRAS family of transcription factors in woody fruit trees, such as its expression in response to flower induction in apple.
Apple is one of the most widely cultivated and economically important fruit trees in the temperate regions of the world. In apple production, flower induction and flower buds formation sustain two growing seasons. Flower induction is a key stage, and always limits fruit yield (Guitton et al., 2012). Apple (Malus domestica Borkh.) cv. "Fuji" is one of the most popular cultivars, and it accounts for 65% planting areas in China. "Fuji" has trouble in flowering induction. Thus, characteristic the molecular mechanism of apple flower induction is necessary. Flower induction was affected by environmental and internal factors (Guitton et al., 2012;Xing et al., 2014Xing et al., , 2015Fan et al., 2016). Sugar and hormones were involved in flower induction (Guitton et al., 2012;Xing et al., 2014Xing et al., , 2015Fan et al., 2016;Li et al., 2016;Zhang et al., 2016). Among all kinds of hormones, including auxin, cytokinins, abscisic acid, GA and other hormones. GA was defined the most important hormones in regulating flower induction and GA pathway was recognized as one of the flowering pathways (Guitton et al., 2012;Porri et al., 2012;Li et al., 2016;Zhang et al., 2016). 6-BA and sucrose, which can promote flower bud differentiation, were also important and widely used (Xing et al., 2015;Li et al., 2016). Several floweringrelated genes families have been identified and characterized in the apple genome, including SPL, IDD, and MADX-box Kumar et al., 2016;Fan et al., 2017). GRAS genes were reported to partially or entirely involve in flower induction in various plants (Hynes et al., 2003;Foster et al., 2007;Vandenbussche et al., 2007;Galvao et al., 2012;Yamaguchi et al., 2014). Little is known about MdGRAS genes and their potential involvement in apple flower induction. The recently published genome sequence of apple (Velasco et al., 2010) has made it possible to characterize the structures, functions, evolution of apple GRAS genes and their response to flower induction. Moreover, it is also available to utilize candidate MdGRAS genes for further genetic improvement aiming at flowering induction.
In this study, we performed a genome-wide survey to identify GRAS genes in the apple genome. A systematic analysis including gene classification, gene characterization, gene structures and gene phylogenies were determined. We further investigated their expression levels of all the 127 putative MdGRAS genes among different tissues and varieties during apple growth and development with ArrayExpress database. Moreover, eight MdGRAS genes were further analyzed.
We investigated their expression levels of the candidate eight genes among different tissues (stems, leaves, flowers, fruits, and buds), different flowering-related treatments (GA 3 , 6-BA and sugar ) and different flowering varieties (Yanfu No. 6 and Nagafu No. 2). Changes in response to these flowering-related treatments indicated their potential involvement in apple flower induction. Our results provided basic clues for further analyses of MdGRAS in apple growth and development.

Identification of GRAS Genes in Apple
Hidden Markov Model (HMM) searches (Finn et al., 2011) were performed to identify GRAS genes in the apple genome. The HMM profiles (PF03514) were downloaded from the Pfam database (http://pfam.sanger.ac.uk). Additionally, AtGRAS protein sequences were used as queries for BLASTP search against the Apple Genome Database (https://www.rosaceae.org/) with an e-value < 1e-4 to identify further possible members of GRAS gene as previous investigation (Li et al., 2011;Jia et al., 2013;Cui et al., 2015;Xu J. N. et al., 2016). Arabidopsis GRAS genes were downloaded from TAIR (http://www.arabidopsis.org/). Finally, genes with the same accession number were removed. All the identified GRAS genes were further manually checked for the GRAS domain.

Gene Structures and Locations Determination
Annotations for candidate GRAS genes were obtained from the Genome Database for Rosaceae (GDR) (http://www.rosaceae.org), and chromosome locations were retrieved from these annotations. MapDraw was used to show the accurate locations of the genes on each chromosome. The online software Gene Structure Display Server (GSDS: http://gsds.cbi.pku.edu.cn) was used to generate the exon-intron structures of all the GRAS genes.

Chemical Characteristics, Elements in the Promoters and Phylogenetic Analysis
Prosite ExPaSy (http://web.expasy.org/protparam/) was used to predict the potential chemical characteristics of the MdGRAS genes. A phylogenetic tree was generated using MEGA 6.06 with the Maximum Likelihood (ML) method and 1000 bootstrap replications. The upstream of 1.5 Kb were used for cis-elements in the promoters of the candidate MdGRAS genes. PlantCARE software (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) was used for searching regulatory elements.

Tandem Duplication and Synteny Analysis
Tandem duplication and synteny analysis were conducted using Circos v. 0.63 (http://circos.ca/). Tandem duplication of the MdGRAS genes was identified according to their physical locations on individual chromosomes in the apple genome. Adjacent homologous GRAS genes on the same apple chromosome with no more than one intervening gene were considered tandem duplicates. The Plant Genome Duplication Database (http://chibba.agtec.uga.edu/duplication/) was used to identify syntenic blocks.

Plant Material
In this study, 72 uniform 6-year-old "Fuji"/T337/Malus robusta Rehd. apple trees were randomly divided into three groups: 18 were treated with GA 3 , 18 were treated with sucrose, 18 were treated with 6-BA and 18 were sprayed with water as a control. These tress were grown at the experimental orchard of Northwest A& F University in Yangling (108 • 04 ′ E, 34 • 16 ′ N), China. And the orchard was well managed. Each group consisted of three blocks, with three replicates. The experiment was conducted from 30 to 70 days after full bloom (DAFB) in 2015. (1) GA 3 treatment was performed as described by Zhang et al. (2016) with a slight modification: 700 mg L −1 GA 3 (Sigma, Deisenhofen, Germany) was sprayed once on a clear morning at 30 DAFB (May 9). (2) 300 mg L −1 6-BA (Sigma, Deisenhofen, Germany) was also sprayed on a clear morning at 30 DAFB (May 9). (3)For sugar treatment, trees were sprayed twice on clear mornings at 30 and 37 DAFB (May 9 and May 16) with 15,000 mg L −1 and 20,000 mg L −1 sucrose. All treatments were performed on the whole plant and were applied with a low-pressure hand-wand sprayer. Terminal buds on current-year spurs (<5 cm), which were chosen according to our previous studies (Xing et al., 2014(Xing et al., , 2015Li et al., 2016;Zhang et al., 2016), were collected into liquid nitrogen at 30, 40, 50, 60, and 70 DAFB, and then stored at −80 • C for use in gene expression analyses.
Different organs were also collected for tissue-specific expression analysis. Flowers were collected at full blossom on April 9 in 2015. Stems were collected from new shoots with a diameter of 2-3 mm. Mature leaves were collected from the adjacent buds. Fruits were collected with a diameter of 2-3 cm. All tissues were immediately frozen in liquid nitrogen and stored at −80 • C for gene expression analysis. Additionally, the online ArrayExpress database (https://www.ebi.ac.uk/arrayexpress/, E-GEOD-42873) was also employed to investigate their expression patterns among different apple varieties and tissues.

Flowering Rate Analysis
Flowering rate was also investigated among different exogenous treatment and two different flowering varieties. One hundred twenty terminal buds on short shoots which were less than 5 cm and no fruits on them were tagged randomly at 15 DAFB. The number of flowers was counted of the tagged shorts as finally flowering rate. And it was investigated at full blossom on April 10 in 2016.

RNA Extraction and cDNA Synthesis
Total RNA was extracted from plant tissue samples using the cetyltrimethyl ammonium bromide (CTAB) method with slight modifications (Gambino et al., 2008). Briefly, 900 µL extraction buffer (2% CTAB, 2.5% PVP-40, 2 M NaCl, 100 mM Tris-HCl [pH 8.0], 25 mM EDTA [pH 8.0], and 2% β-mercaptoethanol) was preheated at 65 • C and added to 2-mL microcentrifuge tubes just before use. Samples containing 200 mg of bud tissue stored at −80 • C were ground to a powder, added to the tubes, and mixed with extraction buffer. After shaking and inverting each tube vigorously for 5 min and incubating at 65 • C for 30 min, an equal volume of chloroform:isoamyl alcohol (24:1, v/v) was added. Each tube was shaken and inverted vigorously and then centrifuged at 12,000 × g for 10 min at 4 • C. The supernatant was collected into a new tube and re-extracted with an equal volume of chloroform:isoamyl alcohol (24:1, v/v). The resulting supernatant was then transferred into a new 2-mL tube and LiCl (3 M final concentration) was added. The mixture was incubated at −20 • C for 4 h and the RNA was selectively pelleted by LiCl after centrifugation at 18,000 × g for 20 min at 4 • C. The pellet was resuspended in 500 µL of SSTE buffer (10 mM Tris-HCl [pH 8.0], 1 mM EDTA [pH 8.0], 1% SDS, and 1 M NaCl) preheated to 65 • C and an equal volume of chloroform:isoamyl alcohol. The mixture was then centrifuged at 12,000 × g for 10 min at 4 • C. The supernatant was transferred to a new microcentrifuge tube, and the RNA was precipitated with 2.5 volumes of cold ethanol at −80 • C for at least 30 min and centrifuged at 1,200 × g for 20 min at 4 • C. Finally, the pellets were washed with 70% ethanol and resuspended in diethylpyrocarbonate-treated water.
Total RNA integrity was verified by running the samples on 2% agarose gels. First-strand cDNA was synthesized from 1 µg of total RNA using a PrimeScript RT Reagent kit with gDNA Eraser (Takara Bio, Shiga, Japan) following the manufacturer's instructions.

Expression Analysis by qRT-PCR
The expression levels of all identified MdGRAS genes were analyzed by qRT-PCR using primer pairs designed with Primer 5.0 ( Table S1). The real-time PCR assay mix (20 µL) consisted of 2-µL cDNA samples (diluted 1:8), 10 µL 2× SYBR Premix ExTaq II (Takara Bio), 0.8 µL of each primer (10 µM) ( Table 1), and 6.4 µL distilled deionized H 2 O. Each PCR assay was run on an iCycler iQ Real Time PCR Detection System (Bio-Rad) with an initial denaturation at 95 • C for 3 min, followed by 40 cycles at 94 • C for 15 s, 62 • C for 20 s, and 72 • C for 20 s. The resulting fragments were subjected to melting-curve analysis to verify the presence of gene-specific PCR products. The meltingcurve analysis was performed directly after real-time PCR and included an initial step of 94 • C for 15 s, followed by a constant increase from 60 • to 95 • C at a 2% ramp rate. The apple EF-1α gene (GenBank accession no. DQ341381) was used as an internal control to normalize all mRNA expression levels. The 2 − Ct method was used to calculate the relative amount of template present in each PCR amplification mixture (Livak and Schmittgen, 2001).

Statistical Analysis
Data were analyzed using variance (ANOVA) at the 5% level with the SPSS 11.5 software package (SPSS, Chicago, IL, USA).

Identification and Characterization Analysis of MdGRAS Genes
To identify members of the MdGRAS gene family, both HMM profiles and BLASTP were used with the default parameters. After manual checking, a total of 127 candidate MdGRAS genes were identified ( Table 1). The genes were named according to their chromosomal locations. The 127 candidate MdGRAS genes were distributed on 16 chromosomes of the apple genome (Figure 1). Chromosome 11 contained the most MdGRAS genes (15.75%), followed by chromosome 3 (11.81%), while chromosome 16 contained no MdGRAS genes (0%) ( Figure S1). Additionally, four genes (MDP0000950387, MDP0000319342, MDP0000640034, and MDP0000640034) could not be located on any of the chromosomes and were named MdGRAS124 to MdGRAS127. Only two genes were distributed on each of chromosomes 7 and 8. The open reading frame (ORF) lengths of the MdGRAS genes ranged from 333 to 5,787 bp, encoding peptides ranging from 110 to 1,928 amino acids. We used ExPaSy to predict the characterization of MdGRAS proteins. The predicted molecular weights of the GRAS proteins ranged from 12.50 to 218.53 kDa. The grand average of hydropathicity was less than 0 for all GRAS proteins. The instability index of most of the 127 MdGRAS proteins was greater than 40, except for MdGRAS22, MdGRAS33, MdGRAS41, MdGRAS42, MdGRAS43, MdGRAS53, and MdGRAS65, indicating that these proteins were unstable.

Phylogenetic and Gene Structure of the GRAS Gene Family
Since Arabidopsis is the most popular model species and the functions of some GRAS genes have been well characterized. Phylogenetic analysis of GRAS genes was investigated in apple Frontiers in Physiology | www.frontiersin.org   and Arabidopsis. A total of 160 GRAS sequences were collected and analyzed. They were classified into eight distinct groups according to the phylogenetic tree: DELLA, SCL3, SCR, Ls, LISCL, PAT1, SHR, and HAM (Figure 2). 14 GRAS proteins were located in the DELLA subfamily, including five AtGRAS (AtGAI, AtRGA, AtRGL1, AtRGL2, and AtRGL3) and nine MdGRAS proteins. However, MdGRAS37 did not contain the integral DELLA and TVHYNP domains. MdGRAS37 was classified into the DELLA subfamily because of its high homology with the other proteins. To examine the gene structures, we analyzed the exon-intron structures of the 127 MdGRAS genes according to their genome sequences ( Figure S2). Thirteen genes (MdGRAS100, MdGRAS101, MdGRAS25, MdGRAS7, MdGRAS102, MdGRAS2, MdGRAS33, MdGRAS51, MdGRAS80, MdGRAS45, MdGRAS29, MdGRAS1, and MdGRAS3) had two or more introns. Other genes had few or no introns.

Expansion Patterns of the MdGRAS Gene Family
Segmental and tandem duplications provide information on the expansion of the gene family. Here, we used the Circos program to evaluate the segmental and tandem duplications of the MdGRAS genes. Tandem duplications were identified from adjacent homologs on a single chromosome, while segmental duplications were identified from homologs on different chromosomes. More than 15 pairs of MdGRAS genes, such as MdGRAS122/53, MdGRAS116/50, and MdGRAS115/49, were located in duplicated genomic regions (Figure 3). All were located in different chromosomes. Chromosomes 11 and 3 had the most duplication, which could partly explain the larger numbers of MdGRAS genes on chromosomes 11 and 3 ( Figure S1). Chromosomes 1, 6, 7, 8, 13, and 14 did not contain any duplicated genes. In summary, segmental and tandem

Comparative Analysis of the GRAS Gene Families in Different Species
To further analyze the evolutionary relationships and expansion patterns of plant GRAS genes, we compared the published GRAS genes in different species (Table S2). The number of GRAS genes was the largest in apple and was close to that in Populus. However, the apple genome is twice the size of the Populus genome. The number of GRAS genes was also similar between Chinese cabbage (48) and P. mume (46), and between grape (52) and tomato (53). However, their genome sizes are also very different (Chinese cabbage: 485.0 Mb vs. P. mume: 280.0 Mb; grape: 487.0 Mb vs. tomato: 760.0 Mb). These results suggested that the number of plant GRAS genes was partly related to genome size.
As Arabidopsis is a model plant species, the functions of the GRAS genes in Arabidopsis have been well researched. We generated a comparative GRAS synteny map between apple and Arabidopsis. Many pairs of syntenic orthologous genes were matched, such as AtSCL32-MdGRAS108, AtRGA-MdGRAS114, and AtPAT1-MdGRAS63 (Figure 4). These genes were detected on 14 chromosomes with duplicated regions, including 4 chromosomes in Arabidopsis and 10 chromosomes in apple.

Tissue-Specific Analysis of MdGRAS Expression Using ArrayExpress Data
To elucidate the potential roles and functions of the MdGRAS genes in apple, we downloaded expression profile data for different tissues from the ArrayExpress database (E-GEOD-42873). Eight organs and samples with two biological replicates were represented in this expression array, including seedlings, roots, stems, flowers, fruit, and leaves ( Figure 5). The GRAS genes showed different expression patterns among different apple varieties and tissues. In the "M74" apple genotype, most of the MdGRAS genes were expressed at higher levels in the flower than in the fruit. Moreover, the expression levels of the MdGRAS genes were higher in ripe fruit (at harvest) than in unripe fruit, indicating their important roles in fruit ripening. Most MdGRAS genes were expressed at lower levels in seedlings, except for MdGRAS100, MdGRAS101, MdGRAS126, MdGRAS18, MdGRAS79, and MdGRAS27. In "M14" leaves, MdGRAS85 showed the highest expression. However, other GRAS genes were hardly detected in leaves. Moreover, most GRAS showed lower levels in stems except MdGRAS26, MdGRAS50, and MdGRAS64.

QRT-PCR Analysis of the Candidate Flowering Related MdGRAS Genes in Different Tissues in "Nagafu No. 2"
We selected eight higher expression genes (MdGRAS6, MdGRAS26, MdGRAS28, MdGRAS44, MdGRAS53, MdGRAS64, MdGRAS107, and MdGRAS122) from our transcriptome data for further analysis (Figure S4, Table S3). These eight genes were from different GRAS subfamilies, including the LS (MdGRAS6), SHR (MdGRAS26), SCL (MdGRAS28), PAT1 (MdGRAS44), DELLA (MdGRAS53, MdGRAS107, and MdGRAS122) and LISCL (MdGRAS64) subfamily. We employed qRT-PCR to investigate their expression patterns in different tissues (stems, leaves, flowers, fruits and buds). All candidate MdGRAS genes showed constitutive expression patterns in the five tissues, except that the levels of MdGRAS26 and MdGRAS28 were low in flowers (Figure 6). The transcription levels of MdGRAS6, MdGRAS26, MdGRAS44, MdGRAS53, MdGRAS64, MdGRAS107, and MdGRAS122 were high in buds and leaves, which indicated their important roles in mediating flower induction. Additionally, MdGRAS107 and MdGRAS28 were all showed lower levels than other candidate genes among the five tissues.

QRT-PCR Analysis of the Candidate MdGRAS Genes in Different Apple Varieties
To examine the expression of the MdGRAS genes in two varieties with contrasting flowering rates, we selected the easyflowering variety "Yanfu No. 6" and the difficult-flowering "Nagafu No. 2." As shown in Table S4, "Yanfu No. 6" has a 1.5fold higher flowering rate than "Nagafu No. 2." We examined the expression patterns in these varieties using qRT-PCR (Figure 7). All candidate genes showed different expression patterns during flower induction. MdGRAS6 expression was higher in the buds of "Yanfu No. 6" during all stages. MdGRAS44 expression was higher in "Yanfu No. 6" except at 80 DAFB, while MdGRAS122 expression was higher in "Nagafu No. 2" except at 70 DAFB. MdGRAS26 expression was lower at the first stage in "Nagafu No. 2, " but was higher in all subsequent stages. The transcription level of MdGRAS107 was higher in "Yanfu No. 6" except at 30 DAFB.

QRT-PCR Analysis of the Candidate MdGRAS Genes in Response to Treatments Affecting Flowering
Previous studies have shown the exogenous treatments (GA, 6-BA, and sucrose) can alter flowering rates, among which GA always decreases the flowering rate, while 6-BA and sugar increase flowering rates (Xing et al., 2015;Li et al., 2016;Zhang et al., 2016). We used these kinds of treatments to further investigate the expression of the candidate MdGRAS genes during flowering. As shown in Table S4, exogenous GA 3 treatment inhibited flowering rate while sugar and 6-BA promoted flowering. The genes displayed various expression patterns at different time points (Figure 8). Some of them were down-regulated, while others were up-regulated. MdGRAS6 was up-regulated after GA 3 treatment except at 40 DAFB and after 6-BA treatment except at 70 DAFB, but was down-regulated after sugar treatment except at 30 DAFB. MdGRAS26 showed similar changes after GA 3 treatment, and was up-regulated except at 40 DAFB. However, the level of MdGRAS26 was downregulated after 6-BA and sugar treatment except at 30 DAFB. MdGRAS28 and MdGRAS44 were down-regulated after GA 3 treatment except that MdGRAS28 was up-regulated at 30 and 50 DAFB. MdGRAS28 was down-regulated after sugar treatment at all time points. The expression of MdGRAS53 was promoted by exogenous GA 3 except at 40 and 80 DAFB, but was inhibited by exogenous 6-BA and sugar at most time points. MdGRAS64 expression was up-regulated after sugar treatment, but was not consistent in response to GA and 6-BA. MdGRAS107 was downregulated after GA 3 treatment except at 30 and 50 DAFB, but was up-regulated by 6-BA except at 30 and 70 DAFB. MdGRAS107 was down-regulated after sugar treatment except at 30 DAFB. The transcription level of MdGRAS122 was lower in 6-BA treated buds than control except at 40 DAFB, and it was down-regulated after sugar treatment at all time points. In order to further confirm the results, we also investigated their expression patterns in GA-treated buds in 2014 and 6-BA treated buds in 2013 based on our previous research Zhang et al., 2016). Samples were collected in the early stage (ES), middle stage (MS), and later stage (LS) of flower induction. As shown Figures S5, S6, their expression patterns in 2014 and 2015 were similar with our results.

Analysis Related cis-Elements in the Candidate MdGRAS Genes
To further analyze the potential roles of MdGRAS genes in response to various circumstances, a 1.5 kb promoter FIGURE 6 | Expression of candidate MdGRAS genes in different tissues. Each value represents the mean ± standard error of three replicates. Means followed by small letters (a, b, and c) are significantly different at a level of 0.05. of candidate MdGRAS genes were used. The related ciselements were identified as previous study (Fan et al., 2017).The candidate MdGRAS genes all shared large light-responsive boxes, followed by stress ( Figure S7). Additionally, hormones-related cis-elements were also found in all the candidate MdGRAS genes, including MeJA, salicylic acid, gibberellin, auxin and ethylene. These identified motifs in the promoter of each gene indicated that MdGRAS may be regulated by various cis-elements within the promoter during growth.

DISCUSSION
Genome-wide identification of GRAS genes has been reported in Arabidopsis, P. mume, Populus, rice, grape, castor bean, and Chinese cabbage (Tian et al., 2004;Song et al., 2014;Huang et al., 2015;Lu et al., 2015;Grimplet et al., 2016;Xu W. et al., 2016), but little in known in apple. In this study, we identified MdGRAS genes in the apple genome (Velasco et al., 2010), and analyzed their phylogenetic relationships, gene characteristics, gene structures, expansion history, and syntenic relationships. We also analyzed their tissue-specific expression patterns and their expression profiles in response to flower induction. Our comprehensive study provides a basic for further investigation of GRAS genes in apple, which can also have potential values in genetic improvement of flower induction in apple as well as some other relative species.

Phylogenesis of Apple GRAS Genes
We identified 127 putative MdGRAS genes in apple, a larger number than that reported for other plant species such as P. mume, Populus, rice, grape, and Chinese cabbage (Tian et al., 2004;Song et al., 2014;Huang et al., 2015;Lu et al., 2015;Grimplet et al., 2016;Xu W. et al., 2016). The larger number of GRAS genes in apple might be because its genome (881.3 Mb) is larger than those of P. mume (280.0 Mb), rice (372.0 Mb), grape (487.0 Mb), and Brassica rapa (283.3 Mb). The MdGRAS genes were located on almost all of the chromosomes except chromosome 16. The structures of most MdGRAS genes were similar, and most had few introns, similar to those in P. mume and Populus (Tian et al., 2004;Song et al., 2014). These findings indicated that the structures of GRAS genes are highly conserved among different plant species. A previous study identified six genes in the DELLA subfamily in apple, based on data in the EST database (Foster et al., 2007). In our study, nine putative genes belonging to the DELLA subfamily were identified from the apple genome  according to the phylogenetic tree (Figure 2). However, MdGRAS37, which was classified into the DELLA subfamily, lacked the conversed DELLA domain and other features of the DELLA subfamily ( Figure S3), including gene length and structural domains. MdGRAS37 was probably placed in the DELLA clade because of its high sequence similarity with Arabidopsis sequences of the DELLA subfamily (Lu et al., 2015). The DELLA domain of MdGRAS37 may have gradually degenerated during the evolutionary process or its sequence may be partly similar with DELLA proteins. However, this needs to be further researched.

Expansion and Synteny of GRAS Genes
Gene duplications can be tandem, segmental, or arise from whole genome duplication. Duplicated genes are known to play important roles in plant growth and development.
To date, various kinds of genes encoding transcription factors have been identified, and most of them have shown gene duplications, including AP2, MADS, SPL, and DOF (Zahn et al., 2005;Moreno-Risueno et al., 2007;Shigyo et al., 2007;Kumar et al., 2016). Gene duplication is important for gene family expansion (Day et al., 2003;Cannon et al., 2004;Leister, 2004). The expansion and duplication of GRAS genes have also FIGURE 7 | Expression of candidate MdGRAS genes in apple buds in the easy-flowering variety "Yanfu No. 6" and different-flowering "Nagafu No. 2." Samples were collected at 30, 40, 50, 60, 70, and 80 days after full bloom (DAFB). Each value represents the mean ± standard error of three replicates. Means followed by small letters (a, b) are significantly different at a level of 0.05. been reported in Arabidopsis, P. mume, and Chinese cabbage (Liu and Widmer, 2014;Song et al., 2014;Lu et al., 2015). It was reported that 85% of AtGRAS genes have undergone segmental duplications, including two tandemly duplicated genes (AtSCL33 and AtSCL34). Here, we used the Circos program to identify tandemly duplicated genes in the apple genome. Synteny analysis of the duplicated blocks in the apple genome showed that some MdGRAS genes were putative segmental duplicates (Figure 2). The present result was consistent with previous study in Arabidopsis, which GRAS also occurred gene duplication in apple genome. Plant GRAS genes were involved in various biological processes, and the expansion of the MdGRAS might contribute to their diversity functions in regulation plant growth and developemt. Additionally, it was reported that a genomewide duplication (more than 50 million years ago) contributed to the expansion of the apple chromosome number, which led to a change to 17 chromosomes from nine chromosomes (Velasco et al., 2010). So the whole genome duplications promoted the expansion of GRAS genes in apple. This expanded by segmental duplication or whole genome duplications of GRAS genes in apple genome was also proved in grape and other species (Grimplet et al., 2016). Totally, the duplication and expansion of GRAS genes have played important roles in gene evolution, and contributed to their varied structures and functions.
To date, little functional research has been performed on rosaceous GRAS genes. Comparison with homologous GRAS genes of the model plant Arabidopsis can help us further analyze the MdGRAS genes. Comparisons of genome sequences among different species can be valuable for the reconstruction of individual gene families (Koonin, 2005). Such genome comparisons can be used to transfer annotations from a betterstudied taxon (whose genome structure and functions have been elucidated) to a less-studied taxon (Lyons et al., 2008). Based on this, we inferred the functions of the less-studied MdGRAS genes from the better-understood AtGRAS genes. We used the Circos program to compare apple and Arabidopsis. Most GRAS genes were detected in syntenic genomic regions (Figure 4), allowing us to predict their functions. Currently, most AtGRAS genes are well understood, including AtSCL3 (Zhang et al., 2011), AtSHR (Levesque et al., 2006), AtSCL13 (Torres-Galea et al., 2006, and AtSLR (Aleman et al., 2008). Thus, we can analyze the potential functions of MdGRAS genes based on their homologs in Arabidopsis. However, functions predicted in this way needs to be experimentally verified.

Expression Profiles of MdGRAS Genes in Different Tissues
To reveal their potential functions in apple growth and development, we employed ArrayExpress data to analyze the expression of the MdGRAS genes in different tissues and organs. These potential investigations of each GRAS can be useful to access its possible function. Their various expression patterns were consistent with their diverse gene locations, lengths and sequences. As shown in Figure 5, the expression patterns of the GRAS genes differed among different tissues and varieties, which is consistent with other species (Tian et al., 2004;Song et al., 2014;Huang et al., 2015;Lu et al., 2015;Grimplet et al., 2016;Xu W. et al., 2016). For example, tomato SlGRAS18 and SlGRAS38 transcript levels were high at the breaker and redripening stages, and their expression was regulated by RIN (Fujisawa et al., 2012(Fujisawa et al., , 2013. In another study, transcripts of 14 SlGRAS genes were detected during fruit ripening (Huang et al., 2015). Our analyses of GRAS expression profiles confirmed that most of the MdGRAS genes had higher transcript levels in fruit than in other apple tissues, indicating their important roles in fruit ripening. MdGRAS9 and MdGRAS85 expression levels FIGURE 8 | Expression of candidate MdGRAS genes in apple buds treated with GA 3 , 6-BA, and sugar. Samples were collected at 30, 40, 50, 60, 70, and 80 days after full bloom (DAFB). Each value represents the mean ± standard error of three replicates. Means followed by small letters (a, b, and c) are significantly different at a level of 0.05.
were extremely high in "M20" during the harvest stage. These results suggested that MdGRAS9 and MdGRAS85 play important roles during fruit development; this was consistent with previous study. In Arabidopsis, three AtGRAS genes (SCL6, SCL22, and SCL27) were shown to be involved in leaf development. In the apple genotype "M14, " MdGRAS85 transcripts were detected at high level in leaves, suggesting that this gene may share similar functions with AtSCL6, AtSCL22, and AtSCL27 in leaf development, which was consistent with their close evolutionary relationships to MdGRAS85 (Figure 2). However, MdGRAS85 did not show high level in leaves of "M49"; the different expressions of MdGRAS85 may be explained by their different genotypes of "M14" and "M49." The high expression patterns of MdGRAS126, MdGRAS18, and MdGRAS79 in seeds were also consistent with PmGRAS20 and PmGRAS36 expression, and suggested that these genes may have important roles in seed development (Lu et al., 2015). These different expression patterns of GRAS indicated their various functions in apple development; further functional analysis needs to be confirmed.
We also analyzed the expression of the candidate MdGRAS genes in different tissues (stem, leaf, flower, fruit, and bud) in "Nagafu No. 2" with qRT-PCR. The MdGRAS genes all showed different expression patterns among the five tissues. Leaves and buds have been used to research flower induction in Arabidopsis and other plants (Galvao et al., 2012;Porri et al., 2012;Xing et al., 2015;Fan et al., 2016). The candidate genes MdGRAS6, MdGRAS26, MdGRAS44, MdGRAS53, MdGRAS64, and MdGRAS122 mostly showed high expression in buds and leaves, which proved their important roles in flower induction. A previous study showed that DELLA has important roles in regulating internode elongation (Bolle, 2004). We found that the expression of MdGRAS107, a DELLA subfamily member, was also higher expressed in the apple stem, which indicates it might also have functions in regulating internode elongation.

Characterization of MdGRAS Expression During Flower Induction
To our best knowledge, the relationship between GRAS and flower induction was received less attention, apart from MdGRL2a. Flower induction is associated with several environmental and internal factors, including temperature, hormones and photosynthesis (Galvao et al., 2012;Porri et al., 2012;Xing et al., 2015Xing et al., , 2016. Hormones and sugar are of the most important factors affected apple flowering, the relationship about GRAS involved in hormones or sugar remains unknown. Up to now, GA was considered to be the most associated with flower induction; 6-BA and sugar were also involved in flower induction. And they were all well applied in orchard cultivation. Additionally, our previous studies showed that exogenous treatments can alter flowering rates; GA 3 can inhibit flowering, while 6-BA and sugar can increase flowering rates (Xing et al., 2015;Li et al., 2016;Zhang et al., 2016). These exogenous treatments mainly alter the nutrition or hormone balance and influence the expression of related genes, which inhibits (GA) or promotes (sugar and 6-BA) flowering (Xing et al., 2015;Li et al., 2016;Zhang et al., 2016). Based on this, we investigated the candidate MdGRAS genes in response to these floweringrelated treatments (Figure 8). With their multiple functions, the candidate MdGRAS genes showed different expression patterns after the treatments, and were down-or up-regulated at most time points. "Yanfu No. 6, " a bud mutant of "Nagafu No. 2, " has a higher proportion of spurs, shorter internodes and a higher flowering rate than "Nagafu No. 2" ( Table S4). The differential expression of the candidate MdGRAS genes in these two varieties with contrasting flowering rates indicated their important roles in flower induction (Figure 7).
In recent years, although the roles of GRAS genes in regulating plant growth and development have become better understood, little is known about the roles of GRAS genes in regulating flower induction in fruit trees. Several studies have shown that GRAS genes are important in bud growth and development. Most DELLA subfamily gene regulate flower induction by the GA pathway, and it can influence the expression of some flowering control genes (Bolle, 2004;Galvao et al., 2012;Porri et al., 2012;Xing et al., 2014Xing et al., , 2015Huang et al., 2015). However, genes from other subfamilies were also associated with flower induction. AtGRAS genes promoted the expression of FT and the della mutant always showed earlier flowering (Galvao et al., 2012). In Arabidopsis, DELLA proteins were shown to inhibit flower induction, and interact with SPL to control flower induction (Shigyo et al., 2007). Previous studies showed that the expression of one-fifth of the candidate genes was increased during flower development in P. mume, and DELLA gene expression was also higher in tomato buds, indicating their important roles in flower induction (Huang et al., 2015;Lu et al., 2015). In our study, we selected eight important MdGRAS genes from transcriptome data. They belong to different subfamilies, including LS (MdGRAS6), SHR (MdGRAS26), SCL (MdGRAS28), PAT1 (MdGRAS44), DELLA (MdGRAS53, MdGRAS107, and MdGRAS122) and LISCL (MdGRAS64). Our result suggested that these candidate MdGRAS genes from different subfamilies were important and they all have potential involvement in flower induction. Additionally, functions of their homologous genes have been described and verified in other species, which provided more evidence for their important response to flower induction (Fu et al., 2001;Hynes et al., 2003;Petty et al., 2003;Wang et al., 2010;Galvao et al., 2012). They can regulate flower induction directly and indirectly. In Arabidopsis, AtLAS and AtSHR were involved in shoot apical meristem. AtSCL6, AtSCL26, and AtSCL27 were highly expressed in flowers (Bolle, 2004); the Arabidopsis rga and gai mutant showed earlier flowering (Galvao et al., 2012). Additionally, GRAS genes in other species were all indicated their potential involvement in flowering induction, such as tomato (SlGRAS11 and SlGRAS18), castor beans (27568.m000253, 28650.m000187, 28650.m000190, and 28650.m000191), and grape (VviSHR1 and VviPAT6) (Huang et al., 2015;Grimplet et al., 2016;Xu W. et al., 2016). Our candidate MdGRAS genes were homologous with these wellknown GRAS genes in other plants, which is in agreement with their important functions in flower induction. However, further functional research needs to be confirmed about these MdGRAS genes. Their different expression patterns during flower induction indicated they were involved in GA-mediated, 6-BA mediated or sugar-mediated flowering in apple trees. Previous transgenic Arabidopsis of MdRGL2a (one of the DELLA members in apple) showed delayed flowering phenotype, which proved their functions in regulating flowering. Additionally, the predication of cis-elements within their promoters indicated their potential roles in response to hormone mediated signal during apple growth and development ( Figure S7). However, with the large number of GRAS gene in Arabidopsis and apple genome, we cannot capture more genes aiming at flower induction except the candidate eight MdGRAS. The other left MdGRAS genes need to be further confirmed. And more GRAS genes need to be functionally researched.
All in all, extensive cross-talk about GRAS genes has been identified in herbaceous and woody plants; extensive research effort has also been devoted to characterization GRAS genes in plant growth and development. However, the GRAS genes of woody economically important tree species have received less attention. In this study, we identified 127 putative GRAS genes from the apple genome, which were located on 16 chromosomes. We classified these genes and performed systematic phylogenetic, structural and synteny analyses. We also analyzed their expression in different tissues (stems, leaves, flowers, fruits and buds). The spatiotemporal expression patterns indicated they potentially have multiple functions in regulating plant growth and development. Additionally, some candidate MdGRAS genes were investigated under treatments with several flowering-related compounds, which indicated the candidate MdGRAS genes are involved in the complex flower induction process. These results may provide useful strategies for further improvement of flower induction in apple.

AUTHOR CONTRIBUTIONS
MH, SF, and DZ conceived and designed the experiment. SF, CG, and MZ performed the experiment. SF, HW, YL, and YS analyzed the data. MH and SF wrote the manuscript.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fphys. 2017.00253/full#supplementary-material