Search tips
Search criteria 


Logo of scirepAboutEditorial BoardFor AuthorsScientific Reports
Sci Rep. 2016; 6: 27640.
Published online 2016 June 21. doi:  10.1038/srep27640
PMCID: PMC4914925

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.

Bamboo, as one of the large woody grasses on earth, belongs to Bambusoideae within Poaceae that includes rice, maize, wheat and other cereals1. As most important non-timber forest resources and fast-growth plants, a total of 75 genera bamboos including 1,250 species were identified in tropical and subtropical areas, playing a crucial role in forest ecosystems2. As one of the most important non-timber forest resources and forest products, moso bamboo (Phyllostachys edulis) is a large woody bamboo that has ecological, economic and cultural value in Asia, which accounts for ~70% of the total bamboo growth area and 5 billion US dollars of annual forest production in China3. Moso bamboo has broadly usages, such as timber, paper, art ware, and the shoots as delicious food. As an significant application, the timber of moso bamboo has been broadly utilized in bamboo textile industry for its physical properties4. It is an important for increasing farmer income and promoting local economic development in bamboo growth and production area. Moreover, with the publication of the draft genome and genomic database of moso bamboo3,5, moso bamboo has been received increasing attention in recent years.

The high-throughput RNA sequencing (RNA-Seq) represents one of essential next-generation sequencing technologies, which reveal a snapshot of RNA presence and quantity from a genome at a given moment in time6,7. As a cost-effective and a high-throughput tool, deep RNA-Seq obtains a comprehensive set of transcribed regions in the genome and yields an amount of information, such as gene expression, novel splice junctions, novel transcripts, and alternative splicing8. Especially, based on the achievement of the draft genome sequence of moso bamboo3, these sequences will facilitate to interpret functional elements and reveal molecular composition of different bamboo tissues. However, the previous studies of gene expression profiles in bamboo focused on single tissue9,10 or comparative analyses between two tissues11,12. To date, the analysis of multi-tissue transcriptome still remains elusive in bamboo.

Among the transcript factors, the GRAS gene family is important one of plant-specific. GRAS transcription factors have participated in plant growth and development processes, such as in gibberellin acid (GA) signal transduction, axillary shoot meristem formation, root radial pattering, male gametogenesis, and phytochrome A signal transduction13,14,15, as well as in plant disease resistance and abiotic stress responses16. GRAS was defined based on nuclear localization, DNA binding, and transcriptional activation features17,18. Based on considerable protein sequences, GRAS is classified into ten subfamilies and its family members have a variable N-terminal domain (N-domain) and a highly conserved C-terminal domain (GRAS-domain), which are divided into five motifs derived from the significant amino acids: leucine-rich region I (LRI), VHIID, leucine-rich region II (LRII), PFYRE and SAW19,20. Moreover, as the increasing species have their complete genome sequenced, the genome-wide analyses for GRAS family had performed in nearly 30 plant species21, such as Arabidopsis, rice, Chinese cabbage, pine, tomato, poplar, Prunus mume and so on22,23,24,25,26,27,28. These provide the possibility of the genome-wide comparative analysis of GRAS family members for evolutionary analysis. Despite the important roles of these genes in plant growth regulation, no findings on the GRAS genes of bamboo have been reported.

Here, we identified a large number of expressed genes in deeply sequencing pool based on RNA-Seq data from the five tissues (leaf, rhizome, root, shoot and panicle) of moso bamboo and further analyzed the bamboo GRAS family with three bamboo leaf samples under high light, both using the Illumina HiSeq 2000 sequencing platform. The tissue-specific genes and DEGs in different tissues were detected, the expression patterns of genes involved in mainly biological processes were further analyzed. Moreover, combined with RNA-Seq data, we applied a considerable analysis of bamboo GRAS members, including sequence phylogeny, gene structure, conserved motifs as well as performed the analysis of PeGRAS expression profiles in response to high light and a co-expression analysis. These results maybe build a foundation for dissecting the gene expression profiles in different tissues of bamboo, provide the evidence of PeGRASs in response to high light, and offer a genome-wide insight into GRAS gene function, which would be helpful in selecting genes for future bamboo breeding via genetic engineering.


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 BambooGDB5 (Bamboo Genome Database, to estimate the profile of expressed genes in each library. The widely used software of TopHat29 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.

Table 1
Summary of RNA-Seq data based on the different tissues collected from moso bamboo.

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

Table 2
Expressed genes detected only in the single tissue.

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

Figure 1
Cluster analysis of gene expression in moso bamboo.

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.

Figure 2
Expression patterns of 28 different families/subfamilies or groups in central metabolic pathways, energy metabolism, transcription factors, and GRAS subfamilies and the statistics of differential expression genes.

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 C2H2, 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 log2 (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.

Figure 3
The log2 (fold change) distribution analysis of genes compared with root and other tissues and the average expression profiles of transcription factor families.

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 PlantTFDB31 (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 cabbage20,22,23. A total of 59 genes encoding GRAS proteins were identified based on BambooGDB (Table 3). In order to compare with the number of GRAS in other species, we selected 75 plants including monocots and dicots (see Supplementary Table S8), totally 3,961 GRAS proteins were analyzed, in which the largest number was in Populus trichocarpa (151). The number of GRAS members in moso bamboo was more than that of Brachypodium distachyon (48), Triticum aestivum (56), A. thaliana (37), Cucumis sativus (43), Brassica rapa (48), Solanum lycopersicum (54), Prunus mume (46) and Vitis vinifera (43), whereas it was less than other species, such as Zea mays (104), Sorghum bicolor (86), and Malus domestica (127).

Table 3
The information of GRAS transcription factors.

Phylogenetic analysis of GRAS

To genome-widely explore the phylogenetic relationship for GRAS TFs among the four plants, including Ph. edulis, O. sativa subsp. Indica, A. thaliana, and B. distachyon (Table 3), we applied Neighbor-Joining method to construct the phylogenetic tree. According to the previous report32, 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.

Figure 4
Phylogenetic analysis of GRAS and conserved motifs in GRAS proteins.

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 GRAS subfamilies with more than 13 members were HAM (Ph. edulis), AtPAT1 (Ph. edulis), and LISCL (B. distachyon).

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

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-W20. 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 treatments16. 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.225,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.


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 genes35. 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 bamboo11,12,36. Moreover, we also defined 270 genes of Arabidopsis thaliana37 and 415 genes of Oryza sativa38 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 lancea39, Oryza longistaminata40 and Zingiber officinale41. 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 data42. 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).

Figure 5
Validation of RNA-Seq data by quantitative real time RT-PCR and expression analysis of GRAS genes in different tissues of moso bamboo.

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 pepper43,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 longiflorum45. 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.


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 provided46. 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 described48. Then, the five pooled libraries from five tissues were 2 × 100 bp by Illumina HiSeqTM 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 Trimmomatic48 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 software49. Then, as the reference, genome sequences and annotation of moso bamboo was downloaded from BambooGDB ( The filtered sequences were mapped to the reference of bamboo genome using TopHat29. 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 Cufflinks50. 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 reads51. 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 studies52,53. Moreover, in the analysis of functional and structural annotation, GO enrichment was carried out using Ontologizer54.

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 ( In addition, the exon-intron structures of genes were performed with GSDS ( The information of GRAS proteins of Arabidopsis and rice were referred to the previous study32. The GRAS protein sequences in Arabidopsis, rice and B. distachyon were obtained from The Arabidopsis Information Resource database (, Rice Genome Annotation Protect ( and PlantTFDB (, respectively. Moreover, all GRAS protein sequences were aligned with ClustalW. The phylogenetic tree was constructed by Neighbor-Joining method of MEGA 6 ( The MEME program (version 4.9.1) was used to predict the motifs in 59 GRAS proteins in moso bamboo (

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 ddH2O. 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[big up triangle, open][big up triangle, open]Ct method58 using NTB as reference gene59.

Additional Information

How to cite this article: Zhao, H. et al. Comprehensive analysis of multi-tissue transcriptome data and the genome-wide investigation of GRAS family in Phyllostachys edulis. Sci. Rep. 6, 27640; doi: 10.1038/srep27640 (2016).

Supplementary Material

Supplementary Information:
Supplementary Dataset 1:


The work was supported by the Special Fund for Forest Scientific Research in the Public Welfare from State Forestry Administration of China [No. 201504106], the National Science Foundation of China [No. 31400557], the National Sci-Tech Support Plan [No. 2015BAD04B01 and No. 2015BAD04B03] and Fundamental Research Funds for International Center for Bamboo and Rattan [No. 1632015008].


Author Contributions H.Z. and Z.G. designed the experiments; H.Z. and L.D. performed bioinformatics analysis; Y.L., L.L. and L.W. made the experimental validation; H.Z., L.D. and Z.G. analyzed experimental data; L.D., H.S. and Z.L. contributed reagents and materials; H.Z., L.D. and Z.G. wrote the paper. All authors read and approved the final manuscript.


  • Han J. et al. Diversity of culturable bacteria isolated from root domains of moso bamboo (Phyllostachys edulis). Microb Ecol 58, 363–373, doi: (2009).10.1007/s00248-009-9491-2 [PubMed] [Cross Ref]
  • Yuming Y., Kanglin W., Shengji P. & Jiming H. Bamboo diversity and traditional uses in Yunnan, China. Mt Res Dev 24, 157–165 (2004).
  • Peng Z. et al. The draft genome of the fast-growing non-timber forest species moso bamboo (Phyllostachys heterocycla). Nat Genet 45, 456–461, 461e1–2, doi: (2013).10.1038/ng.2569 [PubMed] [Cross Ref]
  • Li X. Physical, chemical, and mechanical properties of bamboo and its utilization potential for fiberboard manufacturing. Beijing Forestry Univ (2004).
  • Zhao H. et al. BambooGDB: a bamboo genome database with functional annotation and an analysis platform. Database 2014, bau006, doi: (2014).10.1093/database/bau006 [PMC free article] [PubMed] [Cross Ref]
  • Metzker M. L. Sequencing technologies - the next generation. Nat Rev Genet 11, 31–46, doi: (2010).10.1038/nrg2626 [PubMed] [Cross Ref]
  • Schuster S. C. Next-generation sequencing transforms today’s biology. Nat Methods 5, 16–18, doi: (2008).10.1038/nmeth1156 [PubMed] [Cross Ref]
  • Wang Z., Gerstein M. & Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet 10, 57–63, doi: (2009).10.1038/nrg2484 [PMC free article] [PubMed] [Cross Ref]
  • Peng Z. et al. Genome-wide characterization of the biggest grass, bamboo, based on 10,608 putative full-length cDNA sequences. BMC Plant Biol 10, 116, doi: (2010).10.1186/1471-2229-10-116 [PMC free article] [PubMed] [Cross Ref]
  • Gao Z. M., Li C. L. & Peng Z. H. Generation and analysis of expressed sequence tags from a normalized cDNA library of young leaf from Ma bamboo (Dendrocalamus latiflorus Munro). Plant Cell Rep 30, 2045–2057, doi: (2011).10.1007/s00299-011-1112-0 [PubMed] [Cross Ref]
  • Gao J. et al. Characterization of the floral transcriptome of Moso bamboo (Phyllostachys edulis) at different flowering developmental stages by transcriptome sequencing and RNA-seq analysis. PloS one 9, e98910, doi: (2014).10.1371/journal.pone.0098910 [PMC free article] [PubMed] [Cross Ref]
  • Peng Z. et al. Transcriptome sequencing and analysis of the fast growing shoots of moso bamboo (Phyllostachys edulis). PloS one 8, e78944, doi: (2013).10.1371/journal.pone.0078944 [PMC free article] [PubMed] [Cross Ref]
  • Dill A., Jung H.-S. & Sun T.-p. The DELLA motif is essential for gibberellin-induced degradation of RGA. Proc Natl Acad Sci USA 98, 14162–14167 (2001). [PubMed]
  • Greb T. et al. Molecular analysis of the LATERAL SUPPRESSOR gene in Arabidopsis reveals a conserved control mechanism for axillary meristem formation. Genes Dev 17, 1175–1187, doi: (2003).10.1101/gad.260703 [PubMed] [Cross Ref]
  • Lee M. H. et al. Large-scale analysis of the GRAS gene family in Arabidopsis thaliana. Plant Mol Biol 67, 659–670, doi: (2008).10.1007/s11103-008-9345-1 [PubMed] [Cross Ref]
  • Bolle C. “Functional aspects of GRAS family proteins”, in Plant Transcription Factors, Evolutionary, Structural, and Functional Aspects. ed. Gonzalez D. H. (Cambridge: Elsevier), 295–311 (2016).
  • Itoh H., Ueguchi-Tanaka M., Sato Y., Ashikari M. & Matsuoka M. The gibberellin signaling pathway is regulated by the appearance and disappearance of SLENDER RICE1 in nuclei. Plant Cell 14, 57–70 (2002). [PubMed]
  • Morohashi K., Minami M., Takase H., Hotta Y. & Hiratsuka K. Isolation and characterization of a novel GRAS gene that regulates meiosis-associated gene expression. J Biol Chem 278, 20865–20873, doi: (2003).10.1074/jbc.M301712200 [PubMed] [Cross Ref]
  • Pysh L. D., Wysocka-Diller J. W., Camilleri C., Bouchez D. & Benfey P. N. The GRAS gene family in Arabidopsis: sequence characterization and basic expression analysis of the SCARECROW-LIKE genes. Plant J 18, 111–119 (1999). [PubMed]
  • Bolle C. The role of GRAS proteins in plant signal transduction and development. Planta 218, 683–692, doi: (2004).10.1007/s00425-004-1203-z [PubMed] [Cross Ref]
  • Hirsch S. & Oldroyd G. E. D. GRAS-domain transcription factors that regulate plant development. Plant Signal Behav 4, 698–700 (2009). [PMC free article] [PubMed]
  • Song X. M. et al. Genome-wide analysis of the GRAS gene family in Chinese cabbage (Brassica rapa ssp. pekinensis). Genomics 103, 135–146, doi: (2014).10.1016/j.ygeno.2013.12.004 [PubMed] [Cross Ref]
  • Tian C., Wan P., Sun S., Li J. & Chen M. Genome-wide analysis of the GRAS gene family in rice and Arabidopsis. Plant Mol Biol 54, 519–532, doi: (2004).10.1023/B:PLAN.0000038256.89809.57 [PubMed] [Cross Ref]
  • Abarca D. et al. The GRAS gene family in pine: transcript expression patterns associated with the maturation-related decline of competence to form adventitious roots. BMC Plant Biol 14, 354, doi: (2014).10.1186/s12870-014-0354-8 [PMC free article] [PubMed] [Cross Ref]
  • Grimplet J., Agudelo-Romero P., Teixeira R. T., Martinez-Zapater J. M. & Fortes A. M. Structural and functional analysis of the GRAS gene family in grapevine indicates a role of GRAS proteins in the control of development and stress responses. Front Plant Sci 7, 353, doi: (2016).10.3389/fpls.2016.00353 [PMC free article] [PubMed] [Cross Ref]
  • Huang W., Xian Z., Kang X., Tang N. & Li Z. Genome-wide identification, phylogeny and expression analysis of GRAS gene family in tomato. BMC Plant Biol 15, 209, doi: (2015).10.1186/s12870-015-0590-6 [PMC free article] [PubMed] [Cross Ref]
  • Lu J., Wang T., Xu Z., Sun L. & Zhang Q. Genome-wide analysis of the GRAS gene family in Prunus mume. Mol Genet Genomics 290, 303–317, doi: (2015).10.1007/s00438-014-0918-1 [PubMed] [Cross Ref]
  • Ma H. S., Liang D., Shuai P., Xia X. L. & Yin W. L. The salt- and drought-inducible poplar GRAS protein SCL7 confers salt and drought tolerance in Arabidopsis thaliana. J Exp Bot 61, 4011–4019, doi: (2010).10.1093/jxb/erq217 [PMC free article] [PubMed] [Cross Ref]
  • Trapnell C., Pachter L. & Salzberg S. L. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 25, 1105–1111, doi: (2009).10.1093/bioinformatics/btp120 [PMC free article] [PubMed] [Cross Ref]
  • Saeed A. I. et al. TM4: a free, open-source system for microarray data management and analysis. Biotechniques 34, 374–378 (2003). [PubMed]
  • Jin J., Zhang H., Kong L., Gao G. & Luo J. PlantTFDB 3.0: a portal for the functional and evolutionary study of plant transcription factors. Nucleic Acids Res 42, D1182–1187, doi: (2014).10.1093/nar/gkt1016 [PMC free article] [PubMed] [Cross Ref]
  • Sun X. et al. A functionally required unfoldome from the plant kingdom: intrinsically disordered N-terminal domains of GRAS proteins are involved in molecular recognition during plant development. Plant Mol Biol 77, 205–223, doi: (2011).10.1007/s11103-011-9803-z [PubMed] [Cross Ref]
  • Wu N. et al. Unusual tandem expansion and positive selection in subgroups of the plant GRAS transcription factor superfamily. BMC Plant Biol 14, 373, doi: (2014).10.1186/s12870-014-0373-5 [PMC free article] [PubMed] [Cross Ref]
  • Obayashi T. & Kinoshita K. Rank of correlation coefficient as a comparable measure for biological significance of gene coexpression. DNA research: an international journal for rapid publication of reports on genes and genomes 16, 249–260, doi: (2009).10.1093/dnares/dsp016 [PMC free article] [PubMed] [Cross Ref]
  • O’Rourke J. A. et al. An RNA-Seq based gene expression atlas of the common bean. BMC Genomics 15, 866, doi: (2014).10.1186/1471-2164-15-866 [PMC free article] [PubMed] [Cross Ref]
  • He C. Y., Cui K., Zhang J. G., Duan A. G. & Zeng Y. F. Next-generation sequencing-based mRNA and microRNA expression profiling analysis revealed pathways involved in the rapid growth of developing culms in Moso bamboo. BMC Plant Biol 13, 119, doi: (2013).10.1186/1471-2229-13-119 [PMC free article] [PubMed] [Cross Ref]
  • Arabidopsis Genome, I. Analysis of the genome sequence of the flowering plant Arabidopsis thaliana. Nature 408, 796–815, doi: 10.1038/35048692 (2000). [PubMed]
  • International Rice Genome Sequencing, P. The map-based sequence of the rice genome. Nature 436, 793–800, doi: 10.1038/nature03895 (2005). [PubMed]
  • Huang Q. et al. Differential gene expression between leaf and rhizome in Atractylodes lancea: A comparative transcriptome analysis. Front Plant Sci 7, 348, doi: (2016).10.3389/fpls.2016.00348 [PMC free article] [PubMed] [Cross Ref]
  • He R. et al. A systems-wide comparison of red rice (Oryza longistaminata) tissues identifies rhizome specific genes and proteins that are targets for cultivated rice improvement. BMC Plant Biol 14, 46–46, doi: (2014).10.1186/1471-2229-14-46 [PMC free article] [PubMed] [Cross Ref]
  • Koo H. J. et al. Ginger and turmeric expressed sequence tags identify signature genes for rhizome identity and development and the biosynthesis of curcuminoids, gingerols and terpenoids. BMC Plant Biol 13, 27–27, doi: (2013).10.1186/1471-2229-13-27 [PMC free article] [PubMed] [Cross Ref]
  • Jin J. et al. An Arabidopsis transcriptional regulatory map reveals distinct functional and evolutionary features of novel transcriptionfactors. Mol Biol Evol 32, 1767–1773, doi: (2015).10.1093/molbev/msv058 [PMC free article] [PubMed] [Cross Ref]
  • Engstrom E. M. et al. Arabidopsis homologs of the petunia hairy meristem gene are required for maintenance of shoot and root indeterminacy. Plant Physiol 155, 735–750 (2011). [PubMed]
  • David-Schwartz R., Borovsky Y., Zemach H. & Paran I. CaHAM is autoregulated and regulates CaSTM expression and is required for shoot apical meristem organization in pepper. Plant Sci 203–204, 8–16, doi: (2013).10.1016/j.plantsci.2012.12.011 [PubMed] [Cross Ref]
  • Morohashi K., Minami M., Takase H., Hotta Y. & Hiratsuka K. Isolation and characterization of a novel GRAS gene that regulates meiosis-associated gene expression. J Biol Chem 278, 20865–20873 (2003). [PubMed]
  • Zhao H. et al. Transcriptome and comparative gene expression analysis of Phyllostachys edulis in response to high light. BMC Plant Biol 16, 34, doi: (2016).10.1186/s12870-016-0720-9 [PMC free article] [PubMed] [Cross Ref]
  • Dasgupta P. et al. Nicotine inhibits apoptosis induced by chemotherapeutic drugs by up-regulating XIAP and survivin. Proc Natl Acad Sci USA 103, 6332–6337, doi: (2006).10.1073/pnas.0509313103 [PubMed] [Cross Ref]
  • Severin A. J. et al. RNA-Seq Atlas of Glycine max: a guide to the soybean transcriptome. BMC Plant Biol 10, 160, doi: (2010).10.1186/1471-2229-10-160 [PMC free article] [PubMed] [Cross Ref]
  • Haas B. J. et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis. Nat Protoc 8, 1494–1512, doi: (2013).10.1038/nprot.2013.084 [PMC free article] [PubMed] [Cross Ref]
  • Trapnell C. et al. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc 7, 562–578, doi: (2012).10.1038/nprot.2012.016 [PMC free article] [PubMed] [Cross Ref]
  • Mortazavi A., Williams B. A., McCue K., Schaeffer L. & Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods 5, 621–628, doi: (2008).10.1038/nmeth.1226 [PubMed] [Cross Ref]
  • Chen J., Chen B. & Zhang D. Transcript profiling of Populus tomentosa genes in normal, tension, and opposite wood by RNA-seq. BMC Genomics 16, 164, doi: (2015).10.1186/s12864-015-1390-y [PMC free article] [PubMed] [Cross Ref]
  • Xiong D. et al. Deep mRNA sequencing reveals stage-specific transcriptome alterations during microsclerotia development in the smoke tree vascular wilt pathogen, Verticillium dahliae. BMC Genomics 15, 324, doi: (2014).10.1186/1471-2164-15-324 [PMC free article] [PubMed] [Cross Ref]
  • Bauer S., Grossmann S., Vingron M. & Robinson P. N. Ontologizer 2.0–a multifunctional tool for GO term enrichment analysis and data exploration. Bioinformatics 24, 1650–1651, doi: (2008).10.1093/bioinformatics/btn250 [PubMed] [Cross Ref]
  • Hu B. et al. GSDS 2.0: an upgraded gene feature visualization server. Bioinformatics, doi: (2014).10.1093/bioinformatics/btu817 [PMC free article] [PubMed] [Cross Ref]
  • Tamura K., Stecher G., Peterson D., Filipski A. & Kumar S. MEGA6: Molecular evolutionary genetics analysis version 6.0. Mol Biol Evol 30, 2725–2729, doi: (2013).10.1093/molbev/mst197 [PMC free article] [PubMed] [Cross Ref]
  • Bailey T. L. et al. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res 37, W202–208, doi: (2009).10.1093/nar/gkp335 [PMC free article] [PubMed] [Cross Ref]
  • Livak K. J. & Schmittgen T. D. Analysis of relative gene expression data using real-time quantitative PCR and the 2−ΔΔCT method. Methods 25, 402–408 (2001). [PubMed]
  • Fan C. et al. Selection of reference genes for quantitative real-time PCR in bamboo (Phyllostachys edulis). PloS one 8, e56573, doi: (2013).10.1371/journal.pone.0056573 [PMC free article] [PubMed] [Cross Ref]

Articles from Scientific Reports are provided here courtesy of Nature Publishing Group