Comprehensive analysis of multi-tissue transcriptome data and the genome-wide investigation of GRAS family in Phyllostachys edulis

GRAS family is one of plant specific transcription factors and plays diverse roles in the regulation of plant growth and development as well as in the plant disease resistance and abiotic stress responses. However, the investigation of GRAS family and multi-tissue gene expression profiles still remains unavailable in bamboo (Phyllostachys edulis). Here, we applied RNA-Seq analysis to monitor global transcriptional changes and investigate expression patterns in the five tissues of Ph. edulis, and analyzed a large-scale transcriptional events and patterns. Moreover, the tissue-specific genes and DEGs in different tissues were detected. For example, DEGs in panicle and leaf tissues were abundant in photosynthesis, glutathione, porphyrin and chlorophyll metabolism, whereas those in shoot and rhizome were majority in glycerophospholipid metabolism. In the portion of Ph. edulis GRAS (PeGRAS) analyses, we performed the analysis of phylogenetic, gene structure, conserved motifs, and analyzed the expression profiles of PeGRASs in response to high light and made a co-expression analysis. Additionally, the expression profiles of PeGRASs were validated using quantitative real-time PCR. Thus, PeGRASs based on dynamics profiles of gene expression is helpful in uncovering the specific biological functions which might be of critical values for bioengineering to improve bamboo breeding in future.


Results
Summary of the bamboo transcriptome. To effectively expand our understanding of a global gene expression profile in bamboo, high-throughput RNA-Seq using Illumina HiSeq 2000 was performed on the five tissues of moso bamboo, including leaf (LF), rhizome (RH), root (RT), shoot (SH), and panicle (PN). The statistics shown approximately 668 million reads (~65 Gb) raw data were produced, with an average of 134 million reads (~13 Gb) per tissue. Before short reads alignment, the adaptor sequences and low quality reads were trimmed during sequence preprocessing. Lastly, about 602 million reads (~59 Gb) were acquired as high quality reads. As a significant step for RNA-Seq analysis, clean reads were aligned to the reference genome from BambooGDB 5 (Bamboo Genome Database, www.bamboogdb.org) to estimate the profile of expressed genes in each library. The widely used software of TopHat 29 was employed for sequence alignment. In total, the output demonstrated that ~560.4 million reads, accounted for up to 93.09% of total clean reads, were mapped to the bamboo genome (Table 1). Approximately 6.91% of total reads did not align to the reference genome, likely indicating either setting strictly parameters of alignment, sequencing/assembling errors, gaps in the current genome, or alternative splicing that exists in the reference genome. Moreover, 70.20%, 38.70% and 54.38% of total reads were mapped as perfect match, unique match and multiple position matches, respectively. Analysis of expressed genes. RNA-Seq data could be used in the quantitative analysis of gene expression levels, determined by Fragments Per Kilobase of gene per Million mapped fragments (FPKM). To detect expressed genes in each tissue, we investigated the distribution of gene expression values among the five tissues (see Supplementary Fig. S1). The statistics of expressed genes revealed that the genes of FPKM > 0 accounted for ~91% of the total annotated genes and ~43% of expressed genes were considered as low expression values (0 < FPKM < 1). However, the number of genes with the moderate expression values (1 < FPKM ≤ 100) and high expression values (FPKM > 100) accounted for ~48% of the total annotated genes. Ultimately, the genes with FPKM ≥ 1, ranging from 14,596 (45.63%) to 16,023 (50.09%) in the five tissues, were considered as expressed genes in this study.
Relying on the comparative analysis of expression profile in five tissues, 26,328 (82.31%) and 15,488 (48.42%) genes of the total 31,987 annotated genes in bamboo genome were expressed at a statistically significant value (FPKM ≥ 1) among the different tissues (marked as among-tissue) and within an individual (marked as within-tissue), respectively. Moreover, 209 genes were not expressed (FPKM = 0) in all tissues and 27 expressed genes were detected only in one tissue ( Table 2). The majority of expressed genes in the five tissues provided a significant resource for the dynamic profiles of gene expression in bamboo. Furthermore, an enrichment analysis of Gene Ontology terms was performed using all bamboo genes as the background to explore conservatively biological functions for expression genes in within-tissue (see Supplementary Table S1). Then, the results illustrated that some expressed genes were highly enriched in the processes related to "translation (GO: 0006412)", "catabolic component (GO: 0005575)" and "structural molecule activity (GO: 0005198)".
Change of expressed genes in the different tissues. We performed the clustering affinity search technique (CAST) 30 to elucidate dynamic transcriptome changes in the different tissues of moso bamboo. The results showed 15,488 expressed genes in within-tissue were clustered into 298 groups, with gene numbers within the clusters ranging from 1 to 3835 (Fig. 1a). As shown in Fig. 1b, the cluster less than 50 genes was dominated and accounted for 77.93% of all clusters, while only 3 clusters had more than 1000 genes, in which the number of genes accounted for more than 33% of total expressed genes. On the other hand, 14% clusters (42 clusters) had more than 100 genes. Of these, the top number of genes appeared in the cluster of 100-500 genes (Fig. 1c).
According to the cluster analysis results, some groups of expressed genes shared differentially expressed patterns. Thus, GO enrichment was analyzed in the top 3 clusters to increase our understanding of gene functions. The result indicated that some key biological processes were concentrated in the top 3 clusters, respectively ( Fig. 1d and see Supplementary Table S2-4). For instance, "translation (GO: 0006412)" and "biosynthetic process (GO: 0009058)" were enriched in the Cluster 1. The enrichment of the Cluster 2 contained "regulation of nucleobase-containing compound metabolic process (GO: 0019219)" and "nitrogen compound metabolic process (GO: 0051171)". "Organic substance catabolic process (GO: 1901575)" and "catabolic process (GO: 0009056)" were relatively abundant in the Cluster 3.
Expression patterns in mainly biological processes. To comprehensively unveil expression patterns of gene families or groups, we had selected and analyzed 28 different families or groups based on the pathway predicted by BambooGDB in moso bamboo. The results indicated different families or groups possessing diverse expression patterns. In Fig. 2a, the five groups could be approximately classified into two types based on similar expression patterns. One type represented TCA cycle and glycolysis cycle. Another type composed of pentose phosphate pathway (PPP), glycolysis/glucogeogensis and fatty pathway, in which the maximum and minimum expression emerged in PN and SH. Another feature was that the expression values had a rapid rise and sharply drop in LF to PN and PN to SH, respectively. In the energy metabolism ( Fig. 2b), some significant changes of expression value were undetected in carbon fixation pathway, oxidative phosphorylation, methane metabolism, nitrogen metabolism and sulfur metabolism. Nevertheless, the maximum value of photosynthesis pathway appeared in PN and all the genes in other tissues were obviously down-regulated relative to the panicle sample, indicating the photosynthesis pathway played a vital role in panicles during energy metabolism.
Based on the classification of TF family, eight TF families were selected to analyze the expression patterns (Fig. 2c). The results illustrated that the four TF families, including MYB, NAC, M-type and C 2 H 2 , shared similar expression patterns, in which maximum expression values were found in PN. Whereas the expression values of other TF families in PN, containing GATA, MIKC, GRAS and AP2, were relatively lower than those of other tissues. The expression value of MIKC exhibited an obviously change in PN, indicating that MIKC may play specific roles in panicles. As the subfamilies of GRAS family (Fig. 2d), the expression value of DLT, DELLA, AtSCR and AtSHR had great variations in different tissues and shared similar expression patterns, in which maximum and minimum expression value appeared in LF and PN, respectively. However, the expression value of AtPAT1 and LISCL were higher in PN than in other tissues.
Analyzing differentially expressed genes. Differentially expressed genes (DEGs) were identified based on pair-wise comparison between different tissues, utilizing the following cutoff: fold change ≥ 2 and q-value < 0.05. Thus, the result showed that 3,038 genes were defined as DEGs, accounted for ~0.46% (3,038/665,297) of all pair-wise comparisons ( Fig. 2e and see Supplementary Table S5). As shown in Fig. 2f, the distribution of DEGs depicted 1,607 up-regulated genes and 1,431 down-regulated genes by the paired comparison. The up-regulated genes accounted for 52.90% of all DEGs, in which the set of DEGs was majority in PN and SH. However, the down-regulated DEGs were mainly distributed in RT and PN. Moreover, a total of 173 genes (85 up-regulated genes and 88 down-regulated genes) were identified as unique DGEs in all sets of DGEs (see Supplementary  Table S6). To deep investigate significant expression genes, we applied the log 2 (fold change) distribution analysis by compared with RT vs. other tissues (Fig. 3a). The result, thus, indicated the number of up-and down-regulated DEGs was majorly distributed in 5-10 fold-changes, showing the large part of DEGs within RT vs. other tissues had obvious differences. After filtering duplicated genes, 942 and 517 genes were identified as the up-and down-regulated genes in RT vs. other tissues, respectively. Of all 78 genes involved in photosynthesis metabolism, 55 genes were detected as the up-regulated genes. Except the photosynthesis metabolism, some genes in responses to stresses were also found, such as PH01000000G3800 (heat stress transcription factor B-1), PH01000081G0140 (heat stress transcription factor), PH01000174G0590 (heat stress transcription factor), PH01000239G0570 (abscisic stress-ripening), PH01000944G0560 (stress responsive A/B Barrel domain containing protein), and PH01004112G0160 (stress responsive protein). In the down-regulated genes, 14 genes involved in phenylalanine metabolism were mainly enriched, followed by plant hormone signal transduction and flavonoid biosynthesis. For instance, the genes of PH01002904G0190 and PH01002231G0280 in phenylalanine metabolism were found in each pair-wise comparison. In total, these DGEs maybe perform crucial functions in bamboo growth and development.

Analysis of transcription factors (TFs).
Relying on the annotation of protein function and protein structure in BambooGDB, we identified 1,768 TFs classified into 54 families in moso bamboo. The result was consistent with the prediction by PlantTFDB 31 (see Supplementary Table S7). As shown in Fig. 3b, the average expression value of TF families in all tissues displayed various expression profiles. For example, the average expression of WRKY in shoot was much higher than that of other tissues, while one of WRKY was relatively lower expression in panicle. Moreover, the expression value of HRT-like was undetected in the five tissues and YABBY was not found in LF.
As a plant-specific transcription factor family, GRAS played critical roles in multifarious processes and was regarded as a large family in Arabidopsis, rice, and Chinese cabbage 20 (Table 3), we applied Neighbor-Joining method to construct the phylogenetic tree. According to the previous report 32 , GRAS family was divided into ten subfamilies, containing AtLAS, AtSCL4/7, HAM, AtSCR, DLT, AtSCL3, DELLA, AtPAT1, AtSHR, and LISCL (Fig. 4a). However, the number of GRAS proteins in Ph. edulis (6), O. sativa (8) and B. distachyon (7) was unavailable in the ten subfamilies, respectively, indicating these TFs may belong to a new subfamily. Based on the phylogenetic tree, the member of GRAS in different subfamilies has various characteristics. Among the ten subfamilies, numerous members were concentered in HAM and AtPAT1, while fewer one was AtSCL4/7 in moso bamboo. The comparative result showed that the abundant members were in HAM, AtPAT1, and LISCL subfamilies, whereas the few members in AtLAS, AtSCL4/7, and DLT subfamilies. Moreover, the Analysis of gene structure and conserved motifs in bamboo GRAS. Gene structural diversity is a mechanism for evolution of multi-gene families. For the further studies on the structural diversity of GRAS genes, we analyzed the coding sequences with relevant genome sequences of each GRAS gene in moso bamboo. A detailed diagram of the exon-intron structures was shown in Fig. 4b except for PeGRAS-58. The result demonstrated that the members in different subfamilies had different numbers of intron. For example, HAM subfamily contained one or no introns with exception of PeGRAS-11 including two introns. Whereas, it was diverse in AtPAT1 subfamily that the largest number of intron was four. Most GRAS members had a similar exon-intron structure, indicating that these conserved intron structures may be necessary for the regulation of GRAS gene expression 33 .
GRAS proteins have the five distinct sequence motifs in the C-terminal domain, including LRI, VHIID, LRII, PFYRE and SAW. Among these, the LRI and LRII in the GRAS proteins are two leucine-rich zones. In the most   cases, LRI contains a putative nuclear localization signal. A VHIID motif flanked by LRI and LRII presented in all members of GRAS proteins. The previous studies suggested that LRI, VHIID, and LRII motifs played a critical role in binding between GRAS proteins and their DNA or protein partners. A PFYRE motif consists of three units: P, FY, and RE, and the SAW motif is characterized by three pairs of conserved residues: R-E, W-G, and W-W 20 . Since amino-acid sequence patterns having various biological significances, ten GRAS protein sequences of moso bamboo were chosen to verify the five conserved motifs (see Supplementary Fig. S2).
To further analyze motif compositions of GRAS in moso bamboo, ten conserved motifs were predicted by MEME program (see Supplementary Fig. S3). Most of the GRAS proteins in the same subfamily had conservative motifs, suggesting that they own similar functions. Most of the GRAS proteins in moso bamboo had the ten conserved motifs including VHIID, PFYRE, SAW and so on. However, several members had incomplete motifs. For instance, PeGRAS-18 had only one motif. Thus, the similarity of gene structures and motifs of GRAS proteins provided evidences for the analysis of phylogenetic tree. Moreover, the diversity of motifs compositions between the subfamilies suggest that the functions of GRAS in moso bamboo are diverse.
Expression profiles of PeGRASs in response to high light. The studies shown many GRAS proteins have been participated in biotic/abiotic stresses, showing various response under different environmental conditions and treatments 16 . To deeply investigate PeGRASs in response to high light and provide the valuable evidences for GRAS genes in response to abiotic stresses, we have performed high-throughput RNA-Seq for three samples of moso bamboo leaves treated with high light (1200 μ mol·m −2 ·s −1 ) for 0 h (CK), 0.5 h (0.5H), and 8 h (8H). As shown in Supplementary Fig. S4(a), the cluster results of expression value shown 12 PeGRASs were detected high expression (FPKM > 10 in each sample) and 21 PeGRASs were in lower expression value (FPKM <1 in a sample), suggesting that theses low expressed gene may also perform the essential functions. In addition, with the increased time of high light, PeGRASs continued to be expressed, the down-regulated PeGRASs were abundant in 0.5H and the up-regulated ones were concentrated in CK and 8H, indicating that the expression of PeGRAS may be first repressed in short-term high light, then the expression reach the normal value. Moreover, the results of CAST clustering in Supplementary Fig. S4(b) indicated all PeGRASs were grouped into 6 groups based on expression pattern. Though identifying different expressed patterns for PeGRASs in response to high light, many PeGRASs (~63%) shared one pattern (Cluster 5). These evidences may indicate that PeGRASs exhibit various modules in response to high light, but mainly response pattern was one module, i.e. expression value first declined in initial stage of high light, then ascended during the increasing time of high light. The result of CAST clustering was consistent with that of expression value. Therefore, these findings will help us expand the knowledge for GRAS genes in response to high light stress.
We had performed a co-expression analysis using the relative values of expression in above RNA-Seq data, centered on the average expression, based on a Pearson Correlation Coefficient (PPC) with the cutoff of 0.2 25,34 . The clustered genes were considered under same pattern possessing similar regulatory elements. Thus, the results shown that 56 genes, including 3 PeGRASs, were in the co-expressed gene set, which involved in some TFs, plant hormone signal transduction, ribosome metabolism, and transporters. The annotations of the co-expressed gene set were provided in Supplementary Table S9.

Discussion
Due to the complexity of reference genome as well as limitations of sequencing and alignment methods, it was unavoidable that some reads aligned multi-position. These were ultimately mapped into one position of reference genome randomly, causing some biases in the evaluation of gene expression levels. To obtain better alignment results, the parameters of TopHat were properly set based on genomic architecture and library construction, such as inner distance and standard deviation. In this study, we calculated the two parameters via Trinity and Picard software in each library (see method) (Table 1). Through a comparative analysis of TopHat with/without the two parameters (see Supplementary Table S10 and Fig. S5), the alignment result showed the total number of unique match was increased by 37,787,739 reads (~6.27% of total reads), whereas the average multi-position match was reduced by 19,240,798 reads (~3.19% of total reads). Multi-position matches in the read mapping maybe connect to ancient whole genome duplication. Moreover, a total of the mapping reads was increased by ~3.08%, the unmapped reads was reduced by ~3.08%. In sum, the alignment results were improved by accurately set the two parameters, i.e. inner distance and standard deviation. It contributed to yield more reliable data for further analyses such as expression profile and alternative splicing.
The previous study considered genes exhibiting stable expression profiles in tissues as potential housekeeping genes 35 . Thus, based on the hypothesis, we identified 422 genes with co-variance values less than 0.3 as potential housekeeping genes in different tissues (see Supplementary Table S11) combined with available RNA-Seq data of moso bamboo 11,12,36 . Moreover, we also defined 270 genes of Arabidopsis thaliana 37 and 415 genes of Oryza sativa 38 as the reciprocal best genes with candidate housekeeping genes in moso bamboo, respectively. The results illustrated the pathways contained the candidate housekeeping genes were mainly concentrated in the central pathways (carbohydrate metabolism, energy metabolism, lipid metabolism), transport, translation, and signal metabolism. Then, the dataset of the potential housekeeping genes may be valuable for experimental studies on gene function and regulation in future.
Tissue specific genes in different tissues indicated they may all be functionally related to the tissue growth, development, and metabolic processes. A total of 27 tissue specific genes were identified, including 4 in shoot, 6 in root, 5 in panicle and 12 in leaf. The leaf specific genes were detected from several gene families, including carbonic anhydrase, potassium transporter, MATE efflux family protein and RING-H2 finger protein. These genes might be involved in carbon assimilation, ion transport and stress response. However, no rhizome specific genes were found in this study, which was different from other plants with rhizome such as Atractylodes lancea 39 , Oryza longistaminata 40 and Zingiber officinale 41 . The rhizome specific genes were expected to play important roles in rhizome growth and development, thus the tissue-specific studies need to further investigate.
In the DEGs, the up-regulated genes of PN and LF were mainly concentrated in photosynthesis, glutathione, porphyrin and chlorophyll metabolism. PH01002424G0210 was identified as a DEG in PN vs. SH. Based on the annotation and experimental evidence from a reciprocal best gene (AT2G38310) in A. thaliana, PH01002424G0210 may be considered as abscisic acid sensors, classified in the pathway of plant hormone signal transduction and mainly involved in the development of panicle. Whereas the DGEs in SH and RH were found in glycerophospholipid metabolism, which maybe relate to fast-growth by accelerating cell division and cell growth in the development of shoot and rhizome. PH01002370G0230, whose annotation was transferring acyl groups other than amino-acyl groups, was defined as a unique DEG appeared in SH vs. RT. As a reciprocal best gene of PH01002370G0230, the gene of AT2G26640 in A. thaliana was annotated as function associated with the responses to cold, light and osmotic stress located in membrane indicating PH01002370G0230 located in membrane maybe play an essential role of nutrients transferring and responding to cold and light stresses in the development stage of shoot. Integrated with the growth and development of different bamboo tissues and the DEGs, we predicted that the DEGs might be involved in the regulation of material transport and metabolism in different tissues of bamboo.
Moreover, we had identified a TF network of moso bamboo based on a high-confidence Arabidopsis transcriptional regulatory data 42 . As shown in Supplementary Table S12 and Fig. S6, the results indicated that a total of 97 TFs from 21 families were identified in moso bamboo TFs network. The 6 interactions were found between MYB and NAC, suggesting that both TF families have close connection in transcriptional regulation. Moreover, more self-regulated interactions were detected in NAC and WRKY, suggesting that the TFs in these families may have diverse functions.
In addition, the result demonstrated that the homologous genes encoding AtLAS were not found in moso bamboo by analyzing phylogenetic tree and aligning genomic sequence. However, GRAS members in other three species were detected and clustered into the AtLAS subfamily, which inferred that moso bamboo also has the member of AtLAS subfamily. To further seek the homologous genes, we applied de novo transcriptomic assembling to construct the transcriptomes using Trinity (see Supplementary Table S13). The result showed the homology sequences of AtLAS were still unavailable in transcriptional level. Therefore, we intended to draw the conclusion that the homologous genes of AtLAS maybe lost in moso bamboo.
To experimentally verify the expressed genes based on RNA-Seq, quantitative real time PCR (qRT-PCR) assay was performed using independently collected rhizome tissue, which was in the same developmental stage as those used for the RNA-Seq analysis. Among the 26 selected genes randomly, the Pearson Correlation Coefficient (PCC) based on expression values between qRT-PCR and RNA-Seq data was 0.9491, showing that the RNA-Seq data were highly reliable (Fig. 5a).
In addition, we used RNA-Seq data to analyze different expression levels of GRAS genes in the five tissues of moso bamboo, indicating that their expression value statistically enhanced or repressed in the different tissues (Fig. 5b). Among the 59 genes, most of them were expressed in the different tissues. Some PeGRASs were higher expression in some tissues. For instance, PeGRAS-13 had a high transcript level in RT while PeGRAS-2 and PeGRAS-14 were mostly expressed in SH and RT. All of them were clustered into the HAM subfamily. It was reported that the subfamily was associated with the regulation of the shoot development in A. thaliana, Petunia hybrid and pepper 43,44 . Several PeGRASs show the feature of tissue-specific expressions. For instance, PeGRAS-54 was only expressed highly in LF and SH, while PeGRAS-13 was detected in PN, RT, and SH. The genes from PeGRAS-44 to PeGRAS-53, belonging to LISCL subfamily, were mostly detected in PN. LISCL regulates meiosis-associated genes in Liliuum longiflorum 45 . These indicated that the GRAS genes maybe involve in different mechanisms during plant growth and development.
Taken together, we applied the high-throughput RNA sequencing to investigate the global transcriptomic dynamics and DEGs in five Ph. edulis tissues. Moreover, a total of 1,768 TFs were identified and classified into 54 families in moso bamboo, among which 59 GRAS members were further analyzed, such as sequence phylogeny, gene structure, conserved motifs and gene expression profiles validated by qRT-PCR. We also analyzed the expression profiles of PeGRASs in response to high light and a co-expression analysis. Therefore, these results may provide a key basis for further experimental research on the function of genes and TFs in regulating the growth and development of bamboo, and may help to select the genetic factors valuable for future bamboo breeding by engineering modification.

Methods
The material collection of moso bamboo tissues. The five tissues in moso bamboo were collected, containing leaf (LF), rhizome (RH), root (RT), shoot (SH), and panicle (PN), and 3 individual samples were gathered. Of these, four vegetative tissues (LF, RH, RT, and SH) were picked in spring from the city of Lin'an, Zhejiang Province of China (119°26′ 55.0′ ′ E, 30°19′ 13.4′ ′ N, 480 meter in elevation) and one reproductive tissue (PN) was gathered in autumn from the city of Guilin, Guangxi Province of China (110°31′ 20.2′ ′ E, 25°10′ 42.7′ ′ N, 216 meter in elevation). Moreover, to provide and expand the evidences of the GRAS genes in response to abiotic stress, we have collected three samples of moso bamboo seedlings under same high light intensities (1200 μ mol·m −2 ·s −1 ) for different treated times, including 0.5 hour (0.5H) and 8 hours (8H), and 0 hour as control (CK). The moso bamboo seedlings were potted in our laboratory under long-day conditions (16 h light/8 h dark) at 25 °C, with a light intensity of 200 μ mol·m −2 ·s −1 before high light treatment. The high light was provided by cool white fluorescent tubes. The third leaf on the top of seedlings were selected for the RNA sequencing. The detailed information was provided 46 . All samples were collected and quickly frozen in liquid nitrogen. Furthermore, it took two years to look for the floral tissues in main bamboo growth regions of China due to infrequent and unpredictable flowering events of bamboo. It is hard to collect three individual samples of panicle during the same stage in the given time. Therefore, the above situation brought about a consequence that replicated samples for panicle were unavailable by now.
RNA isolation, cDNA library construction, and sequencing. The total RNA was isolated from five tissues using TRIZOL Reagent Solution (Invitrogen, Carlsbad, CA, USA) on the basis of the manufacturer's instructions. The purity and concentration were detected using a NanoDrop 2000 spectrophotometer. Reverse transcription was conducted with Reverse Transcription System (Promage, USA) 47 . The extracted RNA was treated with RNase-free DNase I for 30 min at 37 °C to remove the residual DNA. cDNA library construction and normalization were performed as previously described 48 . Then, the five pooled libraries from five tissues were 2 × 100 bp by Illumina HiSeq TM 2000 platform (Illumina, San Diego, CA, USA).
De novo assembly of the clean reads, mapping reads to the reference bamboo genome, and estimation of gene expression. Firstly, adaptor sequences and low quality sequences were trimmed using Trimmomatic 48 in the preprocessing. Secondly, to accurately calculate two parameters of TopHat, expected inner distance and standard deviation, we performed de novo assembling of the clean reads by Trinity software 49 . Then, as the reference, genome sequences and annotation of moso bamboo was downloaded from BambooGDB (www. bamboogdb.org) 5 . The filtered sequences were mapped to the reference of bamboo genome using TopHat 29 . Based on insert size calculated by Trinity and the feature of genome, some parameters were optimized. For example, max-intron-length was 2000, min-intron-length was 30, expected inner distance was listed in Table 1. The remaining parameters were utilized as default. Subsequently, the aligned read files were processed by Cufflinks 50 . After reads were assembled into transcripts, their abundance was estimated and normalized using the numbers of reads per kilobase of exon sequence in a gene per million mapped reads 51 . Owing to without replicated samples, we had applied cuffdiff program (one of programs in cufflinks) via setting the parameter of dispersion-method as 'blind' to detect differentially expressed genes. This configuration was suitable for a single biological replicate and was utilized in the previous studies 52,53 . Moreover, in the analysis of functional and structural annotation, GO enrichment was carried out using Ontologizer 54 .
Accession numbers. All sequence data for five tissues (LF, RH, RT, SH, and PN) and three samples under high light from this article have been deposited in the Short Read Archive (SRA) at the NCBI database under the following accession numbers: SRX082501, SRX082502, SRX082503, SRX082504, SRX082507, SRX082508, SRX082509, SRX082510, SRX082511, SRR2035212, SRR2035263, and SRR2035327.
Analysis of GRAS gene structure and construction of phylogenetic tree. The CDS and genome sequences of moso bamboo were downloaded from BambooGDB (http://www.bamboobdb.org). In addition, the exon-intron structures of genes were performed with GSDS (http://gsds.cbi.pku.edu.cn) 55 . The information of GRAS proteins of Arabidopsis and rice were referred to the previous study 32 . The GRAS protein sequences in Arabidopsis, rice and B. distachyon were obtained from The Arabidopsis Information Resource database (http://www.arabidopsis.org), Rice Genome Annotation Protect (http://rice.plantbiology.msu.edu) and PlantTFDB (http://planttfdb.cbi.pku.edu.cn), respectively. Moreover, all GRAS protein sequences were aligned with ClustalW. The phylogenetic tree was constructed by Neighbor-Joining method of MEGA 6 (http://www. megasoftware.net) 56 . The MEME program (version 4.9.1) was used to predict the motifs in 59 GRAS proteins in moso bamboo (http://meme.nbcr.net) 57 .
Expression patterns of GRAS transcription factor in moso bamboo. The primers of GRAS genes from moso bamboo were designed using Primer Premier 5.0 (see Supplementary Table S14). All primers were tested with rTaq (TaKaRa, Japan). The qRT-PCR reactions were performed with a qTOWER2.2 Real-Time PCR System (Analytik Jena, Germany) using Roche LightCycler 480 SYBR Green I Master. The qRT-PCR procedure consisted of 95 °C for 10 min and 50 cycles of 95 °C for 20 s, 60 °C for 10 s. The reaction volume was 10 μ L and contained 5.0 μ L 2 × SYBR Green I Master Mix, 0.8 μ L cDNA, 0.2 μ L forward primer and reverse primer each (5 μ M), and 3.8 μ L ddH 2 O. All reactions were repeated three times. For each condition, the qRT-PCR experiments were performed as biological triplicates. The relative gene expression level was calculated with the 2 −△△Ct method 58 using NTB as reference gene 59 .