Genome-wide identification of growth-regulating factors in moso bamboo (Phyllostachys edulis): in silico and experimental analyses

Growth-regulating factor (GRF), a small plant-specific transcription factor (TF) family, is extensively involved in the regulation of growth and developmental processes. However, the GRF family has not been comprehensively studied in moso bamboo (Phyllostachys edulis), a typical non-timber forest member. Here, 18 GRF genes were identified and characterized from the moso bamboo genome, and they clustered into three subfamilies (A, B and C). PeGRF genes were analyzed to determine their gene structures, conserved motifs and promoter. The non-synonymous/synonymous substitution ratios of paralogous and orthologous were less than 1, indicating that the GRF family mainly experienced purifying selection during evolution. According to the analysis of tissue-specific expression patterns, the participation of moso bamboo GRFs might be required during the formation and development of these five tissues. Moreover, PeGRF proteins might be involved in the regulation of plant development in biological processes. The qRT-PCR analysis demonstrated that PeGRF genes played essential roles in combating hormonal stresses and they might be involved in hormone regulation. PeGRF11, a nuclear localized protein as assessed by a subcellular localization assay, could interact with PeGIF3 in yeast and in planta according to yeast two-hybridization and bimolecular fluorescence complementation assays (BiFC) assays. But PeGRF11, as a TF, had no transcriptional activity in yeast. These results provide useful information for future functional research on the GRF genes in moso bamboo.


INTRODUCTION
Bamboo is an important non-timber forest product worldwide, that financially supports nearly 250 million people (Peng et al., 2013). Moso bamboo, a widely cultivated bamboo species in China, possesses the greatest ecological, economic and cultural values Peng et al., 2013). In China, moso bamboo is widely used for timber, paper, art wares and food (Wu et al., 2016), and its annual forest production has been estimated at approximately five billion US dollars .
Transcription factors (TFs) play essential roles in regulating plants growth and development, including metabolism, differentiation and reproduction. As DNA-binding proteins, TFs regulate gene expression by activating or inhibiting the mRNA transcription of multiple target genes (Filiz, Koç & Tombuloğlu, 2014). To date, many TFs have been identified in moso bamboo and other plants, including the TCP (Chai et al., 2017;Liu et al., 2018), PHD Wang et al., 2015), MYB (Yang, Dai & Zhang, 2012;Cao et al., 2016b;Yang et al., 2019), NAC (Peng et al., 2015;Wang et al., 2016), WRKY (Cai et al., 2017;Li et al., 2017), HD-ZIP Zhao et al., 2014) families, in which their members respond to hormones and participate in plant growth and development. Growth-regulating factor (GRF), a small plant-specific TF family, has been identified to contain two typical domains in the N-terminal region, QLQ (Gln, Leu, Gln; IPR014978) and WRC (Trp, Arg, Cys; IPR014977) (Cao et al., 2016a;Choi, Kim & Kende, 2004;Omidbakhshfard et al., 2015). The QLQ conserved region is also found in the yeast SWI2/SNF2 protein, which interacts with SNF11 to form a complex that participates in chromatin remodeling Zhang et al., 2018). This protein-protein interaction domain, combines with the SNH domain of the GRF-interacting factor (GIF) family, which leads to the formation of the GRF-GIF complex involved in transcriptional activation (Kim, Choi & Kende, 2003). The WRC domain, a plant-specific domain, contains a putative zinc finger structural motif, which is thought to be involved in nuclear localization and DNA binding (Baloglu, 2014;Zhang et al., 2008). Furthermore, less-conserved motifs including FFD (Phe, Phe, Asp) and TQL (Thr, Gln, Leu) are remained in the C-termini of several GRF proteins (Kim, Choi & Kende, 2003;Ma et al., 2017).
In plants, functional analyses of GRFs have been reported, including in root, stem, and leaf development, flower and seed formation and environmental stress tolerance (Omidbakhshfard et al., 2015). OsGRF1, the first reported GRF member, may regulate gibberellic acid (GA)-induced stem elongation and transcriptional activity (Van der Knaap, Kim & Kende, 2000). The GRF proteins also participate in plant morphogenetic periods. In Arabidopsis, AtGRF1 and AtGRF2 overexpression lines show larger leaves and cotyledons, whereas the triple AtGRF1/2/3 mutants produce smaller leaves and cotyledons based on the regulation of cell proliferation and cell expansion (Kim, Choi & Kende, 2003;Kim & Lee, 2006). AtGRF4 has an important effect on the embryonic development of cotyledons and shoot apical meristem, and it is also involved in leaf development, owing to changes in cell numbers (Kim & Lee, 2006). Furthermore, OsGRF4 may regulate panicle traits through its influence on the expression of the CKX5 and CKX1 genes, which are both cytokinin dehydrogenase precursor genes . The overexpression of Brassica rapa L. GRF8 in Arabidopsis causes enlarged leaf organs owing to increased cell numbers . GIFs form functional complexes with GRFs by interacting with QLQ domains, which contribute to the regulation of plant growth and development (He et al., 2017;Kim & Tsukaya, 2015). For example, AtGRF5 cooperate with AtGIF1 to positively regulate the development of leaf primordia (Horiguchi, Kim & Tsukaya, 2005). Additionally OsmiR396d-OsGRF-OsGIF1 plays an important role in the development of inflorescences, which is beneficial to increase the rice yield . Similarly, Li et al. (2016) found that the miR396c-OsGRF4-OsGIF1 module effects grain size determination, which may increase the rice yield. Interestingly, ZmGRF10 plays a negative regulatory role in leaf size and plant height, but it may act as a fine-tuner in maintaining the regulatory functions of GRF-GIF complexes during cell proliferation (Wu et al., 2014). Moreover, GRF genes can also respond to abiotic stresses in plants.
AtGRF7 represses the expression of osmotic stress-responsive genes (including DREB2A) under normal conditions. However, the atgrf7 mutant line enhances tolerances to salt and drought stresses compared with the wild-type and AtGRF7-overexpressing lines .
Here, we identified 18 moso bamboo GRF genes and investigated their phylogenetic relationships, motif organizations, gene structures, evolutionary divergence, promoter elements and expression profiles. Furthermore, we identified a possible functional gene, PeGRF11, which was performed subcellular localization and transcriptional activation assessments. It was then characterized using yeast two-hybridization and bimolecular fluorescence complementation assays (BiFC) assays.
Phylogenetic tree, motif distribution and gene structure A multiple sequence alignment was conducted using ClustalX 2.11 software with the default settings (Tamura, Nei & Kumar, 2004) (File S1). Based on the alignment, a phylogenetic tree was constructed using MEGA 6.0 software with the Neighbor-joining method and a bootstrap analysis of 1,000 replicates (Tamura et al., 2013). At the same time, we used RAxML software to construct a phylogenetic tree, and used the maximum likelihood method (ML) to analyze the data (Stamatakis, 2006). The PeGRF gene structures were analyzed using the online gene structure display server (GSDS 2.0, http://gsds.cbi. pku.edu.cn/) by comparing the genomic DNA sequences and corresponding coding sequences (CDSs) of individual genes (Guo et al., 2007). In addition, using the multiple em for motif elicitation (MEME) program (http://meme-suite.org/tools/meme), with the default parameters of a maximum number of 20 motifs and an optimum width of 6-200 residues, the motifs of PeGRF proteins were predicted.

Cis-acting element analysis and gene ontology annotation
The promoter sequences, 2,000 bp regions upstream of the PeGRF translational start sites were downloaded from the bamboo genome database. The cis-regulatory elements from promoter sequences were analyzed using PlantCARE software (http://bioinformatics.psb. ugent.be/webtools/plantcare/html/) (Higo et al., 1999).

Non-synonymous (Ka)/synonymous (Ks) substitution ratio analysis of homologous pairs
A BLASTN algorithm-based method was used in all-against-all nucleotide sequence similarity searches of three species (Altschul et al., 1997). To identify putative orthologs between two different species (A and B), each sequence from species A was searched against all sequences from species B, and then each sequence from species B was searched against all sequences from species A using BLASTN. If each of two sequences was the best hit to the other and they aligned over more than 300 bp, then they were defined as orthologs . When two PeGRF sequences were aligned over 300 bp and shared at least a 40% identity, the sequences were defined as paralogs. The results were used to calculate the Ka and Ks substitution rates by DnaSP version 5.10.1 (Librado & Rozas, 2009;Rozas, 2009). In addition, a sliding window analysis was performed with the following parameters: window length 150 bp and step size nine bp. The divergence time of the gene pairs was calculated using the formula T = Ks/2λ, where λ = 6.5 × 10 −9 (Peng et al., 2013;Wang et al., 2017).

Plant materials, growth conditions and hormone treatments
A total of 8-week-old moso bamboo seedlings, which grown in the greenhouse under long-day conditions of 16 h light and a continuous temperature of 25 ± 2 C, were used for RNA isolation. The seedlings leaves were sprayed independently with 100 µM gibberellic acid (GA3), 100 µM methyl jasmonate (MeJA), and 100 µM abscisic acid (ABA) solutions and were collected at five time points (1, 3, 6, 12 and 24 h) after the hormone treatment (Khatun et al., 2017;Wu et al., 2016;Wang et al., 2017). The control groups (0 h) were unprocessed leaves before hormone treatments. The treatments were repeated three times and all samples were immediately frozen in liquid nitrogen and stored at −80 C for further use.

Expression profiling of moso bamboo GRF genes
The expression profile data were collected from the NCBI Short Read Archive database (https://trace.ncbi.nlm.nih.gov/Traces/sra/?study=ERP001341). The unprocessed BioProject ERP001341 RNA sequencing reads were trimmed to remove low quality base-calls (Q < 20) and clean adaptor sequences with the pipeline Fastq (Zhang et al., 2015). The paired clean reads were mapped to the Phyllostachys edulis reference genome using the pipeline tophat2, and cufflink, with defaults parameters, was used to detect differentially expressed genes (Trapnell et al., 2012). A heatmap of moso bamboo GRF genes was constructed using the Heatmapper Plus tool (http://bar.utoronto.ca/ntools/ cgi-bin/ntools_heatmapper_plus.cgi) (Toufighi et al., 2005), which included five different tissues or developmental stages.

RNA extraction and qRT-PCR analysis
Total RNA was extracted using an RNAprep Pure Plant Kit (Tiangen, Beijing, China), and the quality was checked by 1.5% agarose gel electrophoresis. Then, the extracted RNAs were reverse transcribed to cDNAs using a Prime-ScriptRT Reagent Kit (TaKaRa, Kusatsu, Japan) in accordance with the manufacturer's instructions. Gene-specific primers were designed using Primer Premier 5.0, and used for qRT-PCR, which was performed as follows: 95 C for 30 s; followed by 40 cycles of 95 C for 10 s, 55 C for 15 s and 72 C for 10 s. The tonoplast intrinsic protein 41 gene acted as an internal control .

Subcellular localization and transcriptional activation
The full-length CDS of PeGRF11 was cloned using the designed primers (Table S1) and was identified by sequencing. Subsequently, the amplified PCR products were digested with Xba I and BamH I restriction enzymes and independently inserted into the pCAMBIA1305-GFP vector containing CaMV35S promoter. These constructs were transformed into Agrobacterium tumefaciens EHA105. Then, Agrobacterium independently harboring the 35S::PeGRF11::GFP and 35S::GFP control vectors were infiltrated into Nicotiana benthamiana leaves using the injection method. The transient expression was examined using a fluorescence microscope (CarlZeiss LSM710; Jena, Germany) after transfected plants had been maintained in the dark for 36-48h.
The cDNA sequences of PeGRF11 were cloned into the pGBKT7 vector using EcoR I and BamH I restriction sites. The co-transform vector (pGBKT7-53 and pGADT7-T) and the empty vector (pGBKT7) were used as positive and negative controls, respectively. The plasmids pGBKT7::PeGRF11, pGBKT7-53+pGADT7-T and the empty vector (pGBKT7) carrying the GAL4 DNA-binding domain and TRP1 nutritional marker were independently transformed into yeast strain Y2HGold. The specific primers were listed in Table S1.

Yeast two-hybridization and BiFC assays
PeGRF11 was cloned into the pGBKT7 bait vector, and the CDS of PeGIF3 was cloned into the pGADT7 vector. Primer sequences were listed in Table S1. Recombinant plasmids were cotransformed into yeast strain AH109 and then plated on selective media (SD/−Leu/−Trp and SD/−Leu/−Trp/−His/−Ade/X-a-GAL) to screen for positive clones (Zeng et al., 2018).
To obtained expression vectors (pFGC-N-YFP and pFGC-C-YFP), eYFP fragments coding for the N-terminal 173 aa and C-terminal 155 aa were cloned into pFGC5941 . The CDSs of PeGRF11 and PeGIF3 were introduced into pFGC-N-YFP and pFGC-C-YFP, respectively, resulting in N-terminal in-frame fusions with N-YFP and C-terminal in-frame fusions with C-YFP, respectively (Table S1). The vector constructs were transformed into Agrobacterium tumefaciens and injected into the leaves of tobacco (N. benthamiana) (Waadt & Kudla, 2008). Images of fluorescence and 4,6-diamidino-2-phenylindole staining of transfected plants were taken using a confocal laser scanning microscope after a 36-48 h dark treatment.

Statistical analysis
Statistical Product and Service Solutions (SPSS) software was used for analyzing the statistically significant differences. All values were presented as means ± standard deviations of three replicates, and significance levels were determined as ÃÃ P < 0.01 and Ã P < 0.05.

Identification and sequence analyses of moso bamboo GRF genes
A BLASTP analysis-based search was performed against the moso bamboo genome database using the consensus protein sequences from the Hidden Markov Model profile as a query. Gene identifier of the GRFs in the NCGR database was provided (Table 1). For easily research, we named the genes PeGRF1-18 in accordance with their scaffold positions ( Table 1). The CDS lengths of the PeGRFs varied slightly: the longest was PeGRF12 (1,656 bp) and the shortest was PeGRF11 (513 bp) (Table 1; Files S2 and S3). The encoded PeGRF proteins ranged from 170 (PeGRF11) to 551 (PeGRF12) amino acids in length, with an average molecular weight of 39,250.5 Da, and the isoelectric points of which varied from 4.92 to 9.57 (Table 1). The details of the GRF genes in rice, maize and Brachypodium distachyon were listed in Table S2.
To further investigate the characteristics of PeGRFs, 18 putative proteins were aligned using DNAMAN software. All PeGRF proteins contained the highly conserved QLQ and WRC domains in their N-terminal regions (Figs. 1A and 1B). Moreover, zinc finger motifs were found within the WRC amino acid stretch (Fig. 1A). Additionally, the FFD and TQL motifs, two smaller stretches of amino acid residues, were found in the C-terminal regions of PeGRF2, -3, -4, -5, -10 and -14 (Fig. 1C).

Phylogenetic analysis of GRF proteins
To evaluate the evolutionary relationships of GRFs between moso bamboo and other species including rice, maize and Brachypodium distachyon, the phylogenetic tree of 54 GRF motif-containing proteins was constructed based on the multiple sequence alignment. As shown in Fig. 2, GRF proteins were divided into three subfamilies on the basis of previous studies in rice and maize (Zhang et al., 2008). Subsequently, we performed a ML phylogenetic tree and found that the GRF genes were divided into the same three subfamilies consistent with the neighbor-joining tree (Fig. S1). Moso bamboo GRFs were relatively evenly distributed among the three subfamilies (A, B and C) (Fig. S2).

Conserved motifs and gene structure analyses of moso bamboo GRF genes
Conserved motifs and intron/exon arrangements were analyzed to further explore the structural diversity of moso bamboo GRFs. In total, 20 distinct motifs in PeGRFs were predicted by the MEME program and were annotated using Pfam online tools. The details of 20 conserved domains, best possible matches and motif lengths were listed in Table S3. Motif 1 and 2 were respectively encoded in WRC and QLQ domains, which were two specific conserved regions of GRF proteins, but the remaining 18 motifs were not. Motif 2 was present in all PeGRFs, except PeGRF18, which was similar to BrGRF . Interestingly, PeGRF genes of the same subgroup possessed similar motifs. For example, most members from subfamily C contained motifs 1, 2, 3, 6, 8 and 9, and they maintained relatively stable structure. The subfamily A also contained similar motif types (Fig. 3).
Subsequently, the intron/exon organizations of PeGRFs were analyzed using the GSDS website to deduce the possible structural evolution. The 18 PeGRFs had different numbers of exons, ranging from two to eight (Table 1). Most of PeGRF (PeGRF2, -3, -6, -8,  -9, -13 and -15) genes had three exons, followed by four (four out of 18) and three (three out of 18) exons (Fig. 4). The intron/exon structures of the 10 paralogous GRF gene pairs were further investigated, and six pairs had been altered. Among them, PeGRF5, -11 and -14 had one less intron compared with PeGRF14, -15 and -4, respectively. However, other homologous pairs had greater differences in the numbers of introns. Thus, the gain and loss of introns during the evolutionary period might result in the structural diversity found among PeGRF genes (Cao et al., 2016a).

Evolutionary patterns of GRF genes
Selective pressure on CDSs were assessed using Ka/Ks ratios, in which Ka/Ks ratios greater than 1, equal to 1 and less than 1 indicated positive, neutral and negative or purifying selection, respectively (Juretic et al., 2005). To further understand the evolutionary selection pressure on the GRF genes in three plant species, 10 pairs of paralogs in moso bamboo (Pe-Pe), 22 orthologous pairs between moso bamboo and rice (Pe-Os), and 20 orthologous pairs between moso bamboo and Brachypodium distachyon (Pe-Bd) were screened by bidirectional best-hit analyses ( Table 2; Table S4). The Ka/Ks ratios of all homologous gene pairs were less than 1, except PeGRF10/BdGRF7 (Fig. 5), indicating that the GRF genes of these species had mainly undergone strong purifying selection.
Remarkably, six paralogous gene pairs of moso bamboo had Ka/Ks ratios of less than 0.5, and the maximum ratio was only 0.716. The Ka/Ks ratios of 31 orthologous pairs  Consequently, we hypothesized that the moso bamboo GRF family underwent a duplication event after its divergence from rice and Brachypodium distachyon (Peng et al., 2013). Sliding-window analyses were also performed to more clearly understand the Ka/Ks ratios of homologous pairs (Figs. S3, S4 and S5).

Gene ontology annotation
Gene ontology functional annotation analysis indicated that ATP binding (GO:0005524) and protein binding (GO:0005515) were enriched in the molecular function category, and the PeGRFs showed protein binding components, which was related to the GRF family protein binding properties. The classification of cellular component showed that all PeGRF proteins were classified as nuclear-related components (GO:0005634), indicating that the moso bamboo GRFs could act as nuclear localization TFs. Additionally, the PeGRFs were involved in 10 different biological processes. In the biological process category, developmental process (GO:0032502) and regulation of transcription, DNA-templated (GO:0006355) were the most abundant terms, both greater than 60% ( Fig. 6; Table S5).

Promoter analysis
The promoter regions 2,000 bp upstream of 17 PeGRF genes were surveyed and this helped us better understand the potential functions and regulatory mechanisms of GRF genes in moso bamboo (Table S6). PeGRF18 was excluded because it did not contain a promoter sequence. The cis-elements of PeGRFs were divided into three major functional categories, plant growth and development, phytohormone responsive, and abiotic or biotic stress (Fig. 7). The Skn-1 and GCN4 motifs, involved in endosperm expression, were found 49 and 14 times, respectively, and accounted for 58% of the plant growth and development-related motifs. The CAT-box and CCGTCC-box cis-elements, which were involved in meristem expression, were detected in the promoters of the PeGRFs, and circadian control elements were also found. Most of PeGRF genes contained MeJA-responsiveness elements, CGTCA and TGACG motifs, which appeared 82 times, representing 56% of the hormone-related cis-acting elements. The GARE and P-box motifs, GA-responsive elements, were observed in five and seven PeGRF genes, respectively, which represented 14% of the hormone-responsive elements of PeGRFs. Moreover, ABA-responsive element (ABRE, motif IIb and CE1), SA-responsive element (TCA-element), and uxin-responsive element (TGA-element, AuxRR-core, and TGA-box) were presented in 10, nine and eight PeGRF genes, respectively. In addition, some abiotic and biotic stress-related cis-acting elements were also identified in the putative promoter regions of the PeGRF genes. Heat-stress-responsiveness elements and low-temperature-responsiveness elements (ARE and GC motifs) appeared 20 and 39 times, accounting for 21% and 40% of the stress-related cis-acting elements, respectively. However, cis-acting elements related to fungal elicitor-responsive (Box-W1) and wound-responsive (WUN motif) elements were only identified in four and one genes, respectively.
Expression analysis of PeGRF genes exposed to hormone treatments Plant hormones play vital roles in plant growth and development. Furthermore, OsGRF1, discovered 18 years ago, was predicted to be a Gibberellin 3 (GA3)-induced gene (Van der Knaap, Kim & Kende, 2000). To investigate whether moso bamboo GRF genes had roles in response to phytohormones, qRT-PCR was performed on the 18 PeGRF genes in moso bamboo leaves 0, 1, 3, 6, 12 and 24 h after hormone treatments. The specific primers for moso bamboo GRF genes were listed in Table S8. The transcriptional levels of the 18 PeGRF genes exhibited different degrees of change after hormone treatments (Figs. 9, 10 and 11). A total of 13 PeGRF genes showed up-regulated expression levels in response to the GA3 treatment. In particular, the expression levels of 12 genes peaked at 6 h, with the exception of PeGRF18 (Fig. 9). PeGRF1, -9 and -12 exhibited strongly up-regulated levels at five time points compared with at 0 h (Figs. 9A, 9I and 9L). Additionally, PeGRF13 and PeGRF14 were respectively up-regulated by 10-fold and fivefold at 6 h compared with at 0 h and were then down-regulated (Figs. 9M and 9N). PeGRF6 was only down-regulated at 12 h by the GA3 treatment and was up-regulated by ninefold at 6 h ( Fig. 9F). In contrast, the expression levels of three members (PeGRF2, -11 and -16) were prominently down-regulated after the GA3 treatment (Figs. 9B, 9K and 9P).
Moso bamboo GRF genes showed different responses following the ABA treatment (Fig. 11), with 12 genes being significant up-regulated more than fourfold compared with

Subcellular localization of PeGRF11
Most TFs localize to the nucleus and play important roles in regulating target genes. PeGRF11 exhibited high expression levels in moso bamboo leaves and was significantly down-regulated after hormone treatment. The full-length CDS of PeGRF11 was amplified by specific primers (Table S1) and cloned into the pCAMBIA1305-GFP vector. Then, Agrobacterium tumefaciens independently harboring 35S::PeGRF11::GFP and the 35S::GFP control vector were infiltrated into N. benthamiana leaves. The 35S::PeGRF11:: GFP fusion proteins were localized in the nucleus by GFP signals, while the expression of 35S::GFP alone was detected throughout the whole cell (Fig. 12). Transcriptional activation of PeGRF11 To research the transcriptional activity of PeGRF11, the constructed pGBKT7::PeGRF11, positive control pGBKT7-53+pGADT7-T and empty pGBKT7 vectors were independently transformed into yeast strain Y2HGold, and these strains were inoculated onto SD/−Trp and SD/−Trp/−His/−Ade/X-a-GAL media, respectively. The three transformants containing pGBKT7::PeGRF11, pGBKT7-53+pGADT7-T (positive control) and pGBKT7 (negative control) vectors grew well on the SD/−Trp medium. On the SD/−Trp/ −His/−Ade/X-a-GAL medium, the positive control group (pGBKT7-53+pGADT7-T) exhibited visible blue colonies with satisfactory growth states, while the pGBKT7:: PeGRF11 and negative control transformants were unable to survive (Fig. 13). Thus, PeGRF11 did not possess transcriptional activity in yeast strains.

PeGRF11 interactions with PeGIF3 as cotranscription factors
Growth-regulating factor proteins interact with GIFs through the QLQ domain, and the GRF-GIF complex regulates leaf growth and morphology (Wu, Wang & Zhuang, 2017).
To gain insights into the regulatory mechanisms of the GRF genes in moso bamboo, we performed yeast two-hybrid and BiFC assays. As depicted in Fig. 14A, BD-PeGRF11 and AD-PeGIF3 co-transformed yeast (AH109) cells grew well and turned blue on SD/−Leu/−Trp/−His/−Ade/X-a-GAL selective medium, similar to the positive control (pGBKT7-53+pGADT7-T), suggesting that PeGRF11 interacted with PeGIF3 in yeast. Furthermore, the constructed PeGRF11-N-YFP and PeGIF3-C-YFP vectors were co-transformed into tobacco (N. benthamiana) leaves for the BiFC assays. YFP signals were clearly distributed in the nuclear compartment based on staining with 4,6-diamidino-2-phenylindole (Fig. 14B). In contrast, YFP signals were not detected in the two negative controls (PeGRF11-N-YFP co-expressed with unfused C-YFP and PeGIF3-C-YFP co-expressed with unfused N-YFP). Thus, PeGRF11 physically interacted with PeGIF3 in yeast and in planta.
The genome-wide identification of GRF gene family has been reported in various plant species, including rice (Choi, Kim & Kende, 2004), maize (Zhang et al., 2008), Arabidopsis (Kim, Choi & Kende, 2003), Chinese cabbage , Chinese pear (Cao et al., 2016a) and tomato (Khatun et al., 2017). Nevertheless, little was known of the GRF gene family in moso bamboo. Here, 18 GRF genes containing the highly conserved QLQ and WRC domains in their N-terminal regions were identified from moso bamboo genome (Table 1; Fig. 1). The 54 GRF genes from rice, maize, Brachypodium distachyon and moso bamboo clustered into three major subfamilies based on the phylogenetic analysis, with moso bamboo GRF genes distributed across all three subfamilies (Fig. 2). Moso bamboo GRF proteins clustered onto relatively short branches with rice or Brachypodium distachyon, comparing with maize, which indicated that the evolutionary relationships between the moso bamboo GRF family in and those of rice and Brachypodium were closer than with the GRF family of maize (Peng et al., 2013). Thus, the GRF genes, which might share a common ancestor, had gone in distinct evolutionary directions among the different lineages. Furthermore, motif compositions and exon/intron structures among members of the same subfamily possessed certain similar features. In contrast, GRF genes from different subfamilies did not share consistent features. These results were in accordance with previous studies in tobacco (Zhang et al., 2018) and Chinese pear (Cao et al., 2016a), which further supported the close evolutionary relationships within GRF subfamilies and the dependability of our phylogenetic analysis. Additionally, the similar characteristics within the subfamilies might indicate that the moso bamboo GRF members had similar functions in plant growth and development.

Molecular characterization of GRFs
The QLQ domain, a protein-binding region, could combine with GIFs to form a complex involving in the regulation of plant growth and development. AtGRF1 interacted with AtGIF1 to regulate the growth and shape of leaves and petals . The interplay between maize GRF10 and GIF2 formed a complex that might participate in the plant morphology formation process (Wu et al., 2014). PeGRF11, a homolog of ZmGRF10, was cloned, and we isolated PeGIF3, which was high consistency with ZmGIF2, to investigate this hypothesis in moso bamboo (Fig. S6). The yeast two-hybridization and BiFC assays revealed that the PeGRF11 and PeGIF3 proteins interacted with each other to form a complex that might be involved in the regulatory network of moso bamboo (Fig. 14). The WRC domain served as a nuclear localization signal (Omidbakhshfard et al., 2015). Nuclear localization signals had been confirmed for several GRF members, such as OsGRF4 and OsGRF10 (Kuijt et al., 2014;Sun et al., 2016). In GO annotation analysis, the nuclear-related components were abundant in the cellular component category, which was related to the functional domain of the GRF transcription factor. PeGRF11 contained a WRC domain in its N-terminal region (Fig. 1), and its fusion protein was visibly expressed in the nucleus, which further verified this feature (Fig. 12).
Additionally, the C-terminal regions acted as transactivation domains, as shown by the transactivation activity assay of AtGRF1 and AtGRF2 fragments . ZmGRF10, defined as endogenous truncated gene owing to the loss of the entire C-terminal region, showed no transcriptional activity in yeast (Wu et al., 2014). Similarly, PeGRF11 was not transcriptionally activated in yeast probably owing to the lack of a C-terminal domain ( Fig. 13; Fig. S6). Thus, we hypothesized that the loss of the C-terminal region of PeGRF11 might occur during a gene duplication event.

Potential functions of GRFs in moso bamboo
The expression profiles of PeGRF gene in different tissues provided valuable clues for predicting gene functions during plant growth and development. GO annotation analysis indicated that the moso bamboo GRF genes involved in different biological processes, including leaf development (GO:0048366), seed development (GO:0048316), developmental process (GO:0032502), etc. These suggested that the moso bamboo GRF genes might have functional differentiation. GRF genes showed high expression levels in relatively actively growing tissues, such as germinating seeds and buds (Horiguchi, Kim & Tsukaya, 2005;Kim, Choi & Kende, 2003;Kim et al., 2012;Zhang et al., 2018). Based on a microarray data analysis, PeGRF genes were expressed more highly in moso bamboo leaf and early panicle where cell proliferation vigorously occurred in. Similar phenomenon occurred in tobacco (Zhang et al., 2018) and Brassica napus . Arabidopsis overexpressing AtGRF5 exhibited larger leaves than normal and could enhance cell proliferation in leaf primordia (Horiguchi, Kim & Tsukaya, 2005). However, the overexpression of ZmGRF10 resulted in smaller leaf size and shorter plant heights because of the decrease in cell proliferation (Wu et al., 2014). Therefore, GRFs might regulate leaf size by enhancing or decreasing cell proliferation. Furthermore, PeGRF11, a homolog of ZmGRF10 in moso bamboo, was highly expressed in leaves, which might affect biological pathways by involving in plant growth and development.
Plant hormones, which are metabolic products in plants, are involved in plant growth and development. GRFs played roles in responding to plant hormones (Khatun et al., 2017; Van der Knaap, Kim & Kende, 2000). Seven OsGRFs were strongly up-regulated owing to the application of GA3 (Choi, Kim & Kende, 2004). However, AtGRF gene expression levels did not remarkably change after GA3 treatments (Kim, Choi & Kende, 2003). NtabGRF genes significantly respond to various hormones, including GA3, ABA, IAA, BR and 6-BA (Zhang et al., 2018). Most Camellia sinensis GRF genes were up-regulated or maintained when treated with salicylic acid or IAA, whereas only one gene was up-regulated after a GA3 treatment (Wu, Wang & Zhuang, 2017). Here, the expression levels of nine PeGRF genes showed remarkable changes after GA3 treatments; six members were significantly up-regulated to fourfold or more compared with at 0 h, and three members were significantly down-regulated to less than 0.5-fold (Fig. 9). Similarly, PeGRF genes had significant responses to MeJA and ABA treatments (Figs. 10 and 11). Thus, moso bamboo GRF members might be necessary for regulatory hormone feedback mechanisms. Promoter analysis also provided important information on potential functions and gene expression changes in response to hormone treatments. GARE and P-box motifs, which were GA-responsive elements, were identified in the promoter regions of 10 PeGRFs. MeJA-responsiveness elements (including CGTCA and TGACG motifs) were positioned in the promoter regions of moso bamboo GRF genes. In addition, several PeGRFs contained ABRE, IIb and CE1 motifs, which were involved in the ABA responses during the growth process. Therefore, we hypothesized that GRFs might have different roles in the regulation hormonal feedback mechanisms and physiological processes.

CONCLUSIONS
In this study, we identified 18 PeGRF genes from moso bamboo genome, which were further classified into three subfamilies based on phylogenetic relationship analysis. Analyses of motif compositions and exon/intron structures showed that the most members had conserved structural features within the same subfamily. The Ka/Ks ratios of the paralogous (Pe-Pe) and orthologous (Pe-Os, Pe-Bd) gene pairs were less than 1, indicating that the GRF families of these species had undergone strong purifying selection. We also performed the qRT-PCR analysis of moso bamboo GRF genes under hormone (GA3, MeJA and ABA) treatments, which displayed that the PeGRF gene had significant responses to hormone treatment. Besides, promoter analysis of PeGRFs also suggested that they might be necessary for regulatory hormone feedback mechanisms. The GO annotation analysis showed biological process for the PeGRF regulation in plant development. The expression profiles of moso bamboo GRFs were explored by transcriptome data, from which we screened PeGRF11 which had high expression levels in various tissues, and performed subcellular localization and transactivation activity, yeast two-hybridization and BiFC assays. The results revealed valuable clues to the functions of PeGRFs in moso bamboo growth and development.

GRF
Growth-regulating factor GIF GRF-interacting factor BiFC bimolecular fluorescence complementation assays NCBI National Center of Biotechnology Information GA3 Gibberellin 3 MeJA Methyl Jasmonate ABA abscisic acid Ks number of synonymous substitutions per synonymous site Ka number of non-synonymous substitutions per non-synonymous site CDS coding sequence qRT-PCR quantitative real-time PCR.