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

Identification of the GRF gene family in moso bamboo

The GRF protein sequences of rice were collected from the Rice Genome Annotation Project (http://rice.plantbiology.msu.edu/) (Choi, Kim & Kende, 2004), and the genomic data of both maize and Brachypodium distachyon were downloaded from the Phytozome database (http://www.phytozome.net/) (Filiz, Koç & Tombuloğlu, 2014; Zhang et al., 2008). The moso bamboo genomic data were obtained through the National Center For Gene Research (NCGR) (http://server.ncgr.ac.cn/bamboo/index.php) to identify PeGRF genes. Based on Hidden Markov Model profiles, the QLQ (PF08880) and WRC (PF08879) domains from the Pfam database (http://pfam.janelia.org/) were employed to search against the moso bamboo genome database (Choi, Kim & Kende, 2004; Wheeler & Eddy, 2013). The SMART software was used to screen for candidate protein sequences containing the known conserved QLQ and WRC domains (Finn et al., 2016). Subsequently, the physicochemical parameters of the PeGRF proteins, such as the predicted number of amino acids, molecular weight and isoelectric point, were calculated using ExPASy (http://www.expasy.ch/tools/pi_tool.html). Additionally, the sequences of the PeGRF proteins were calculated using the DNAMAN (http://dnaman.software.informer.com/) algorithm.

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

The protein sequences of PeGRFs were annotated using the Blast2GO program (https://www.blast2go.com/) (Conesa & Götz, 2008), and gene ontology (GO) terms were divided into three categories: molecular function, cellular component and biological process.

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 (Guo et al., 2014). 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⁻⁹ (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 (Fan et al., 2013).

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-α-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 (Hu et al., 2013). 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 (Wang et al., 2014). 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 (Pe-Os and Pe-Bd) were less than 0.5. In addition, the relative Ks values varied between 0.086 and 0.716 for paralogous pairs, deducing that the divergence times of moso bamboo GRF genes were concentrated between 7.18 and 36.65 million years ago. For orthologous genes, Pe-Os and Pe-Bd, the Ks values were centrally distributed between 0.2–0.45 and 0.2–0.4, respectively, which indicated that moso bamboo diverged from rice and Brachypodium distachyon 15.38–34.62 and 15.38–30.77 million years ago, respectively. 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 profiles of PeGRF genes

The GRF family was involved in the development of plant organs (Khatun et al., 2017; Wang et al., 2014; Wu, Wang & Zhuang, 2017). Different PeGRF genes were detected in different tissues (Fig. 8; Table S7). Ten moso bamboo genes (PeGRF1, -2, -3, -6, -7, -8, -9, -11, -15 and -18) exhibited high expression levels in five tissues or developmental stages, suggesting that these genes were most likely involved in the formation of these tissues. Subsequently, four paralogous pairs (PeGRF4/PeGRF5, PeGRF4/PeGRF14, PeGRF8/PeGRF13 and PeGRF5/PeGRF14) showed distinct expression patterns, implying that duplicated genes’ functions might have differentiated. Notably, all moso bamboo GRF genes were highly expressed in leaf, with the exception of PeGRF13, potentially indicating that the GRF gene family played an important role in periods of leaf growth and development.

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

In total, nine of 18 PeGRF genes (PeGRF1, -5, -6, -9, -10, -12, -13, -14 and -17) were strongly up-regulated over fourfold by the MeJA treatment, while three genes (PeGRF2, -11 and -16) were down-regulated to less than 0.5-fold compared with the control group (0 h) (Fig. 10).

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 the control (0 h). Notably, PeGRF11 and PeGRF15 had comparatively lower transcriptional levels after the ABA treatment.

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

GRFs in moso bamboo

As TFs regulating plant growth and development, GRFs have been studied since the discovery of OsGRF1 (Omidbakhshfard et al., 2015; Van der Knaap, Kim & Kende, 2000). 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 (Wang et al., 2014), 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 (Kim & Kende, 2004). 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 (Kim & Kende, 2004). 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 (Ma et al., 2017). 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.