Skip to main content
Download PDF

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 Biologyvolume13, Article number:119(2013)Cite this article

Abstract

Background

As one of the fastest-growing lignocellulose-abundant plants on Earth, bamboos can reach their final height quickly due to the expansion of individual internodes already present in the buds; however, the molecular processes underlying this phenomenon remain unclear. Moso bamboo (Phyllostachys heterocyclacv. Pubescens) internodes from four different developmental stages and three different internodes within the same stage were used in our study to investigate the molecular processes at the transcriptome and post-transcriptome level.

Results

Our anatomical observations indicated the development of culms was dominated by cell division in the initial stages and by cell elongation in the middle and late stages. The four major endogenous hormones appeared to actively promote culm development. Using next-generation sequencing-based RNA-Seq, mRNA and microRNA expression profiling technology, we produced a transcriptome and post-transcriptome in possession of a large fraction of annotated Moso bamboo genes, and provided a molecular basis underlying the phenomenon of sequentially elongated internodes from the base to the top. Several key pathways such as environmental adaptation, signal transduction, translation, transport and many metabolisms were identified as involved in the rapid elongation of bamboo culms.

Conclusions

This is the first report on the temporal and spatial transcriptome and gene expression and microRNA profiling in a developing bamboo culms. In addition to gaining more insight into the unique growth characteristics of bamboo, we provide a good case study to analyze gene, microRNA expression and profiling of non-model plant species using high-throughput short-read sequencing. Also, we demonstrate that the integrated analysis of our multi-omics data, including transcriptome, post-transcriptome, proteome, yield more complete representations and additional biological insights, especially the complex dynamic processes occurring in Moso bamboo culms.

Background

Biofuels have been proposed as alternatives to relieve the problem of greenhouse gas emissions and energy shortages [1,2]。The increasing demand for lignocellulosic biomass for the production of biofuels make it necessary to exploit fast-growing and high-yield wood resources. Thus, it is of great importance to understand the underlying mechanism of growth especially in height, which is highly correlated with biomass yield, for lignocellulose-abundant plants. As one of the fastest growing lignocellulose-abundant plants on Earth, bamboos can reach their final height of 5–20 m in a single growing season of 2–4 months due to the expansion of individual internodes already present in the buds [3]。To better understand the growth characteristics and physical properties of these bamboos, numerous studies have focused on the general mode of growth, anatomical structure, hormone distributions and chemical and physical characteristics of the culms [4–13], and have observed sequentially elongated internodes from the base to the top [14]。Some internode elongation-associated genes such asACO1,CENL1,EUI1,OsGLU1,SNORKELandSSD1have been identified in other plants of Gramineae [15–20], whereas limited molecular information has been provided in the subfamily of bamboos. Despite the sequencing of a set of cDNA [21–26], Expressed Sequence Tags (EST) [27,28], protein expression profiles [29], draft genome [30] , RNA-seq [31] and monoclonal antibody banks [32] in bamboo, the molecular processes underlying rapid internode elongation remain unclear.

Measurements of mRNA and miRNA expression levels, clarity of the regulatory relationships between miRNAs and their corresponding mRNA targets are critical to understanding many pathways and biological systems. With the advent of second-generation sequencing-based technologies such as RNA-Seq and Digital Gene Expression (DGE), it is possible to measure a genome-wide dynamic range of expression in an unbiased manner. These technologies have a high sensitivity and reproducibility compared with existing technologies (e.g. DNA microarrays, cDNA or EST sequencing) [33,34], and will undoubtedly lead to novel insights into plant development and biotic and abiotic stress responses.

To comprehensively understand the molecular processes underlying the rapidly internode elongation ofPhyllostachys heterocyclacv. Pubescens at the whole genome level, we have studied the differently expressed proteins in different development stages using two-dimensional gel electrophoresis (2D-PAGE) technology [29], here, we performed RNA-Seq, global mRNA and miRNA expression profiling on four different developmental stages and three different internodes. In addition to gaining more insight into the unique growth characteristics of bamboo, we provide a good case study to analyze gene expression of non-model plant species using high-throughput short-read sequencing. And also, the correlation and concordance of transcript and post-transcriptional levels, protein and transcript levels will be studied to explain complex dynamic development processes in the rapid growth of developing culms in Moso bamboo.

Results

Anatomical and endogenous hormone variation at different developmental stages of internodes

To better understand the anatomical structure of the culms, transverse and longitudinal sections of culms from the defined development stages were observed under a microscope. Four developmental stages or three developmental internodes were defined according to the different unearthed heights of individual plants (0.05, 1.00, 6.00 or 12.00 m; nominated as G1M–G4M, respectively) or different portions (basal internode, G3B; middle internode, G3M; and top internode, G3T) of the same culm (at the G3M stage), respectively. The number of nuclei declined with culm development (from G1M to G4M stages), and almost disappeared in the later stages (G4M stage). Meanwhile, the length of parenchyma and fiber cells increased continuously with advancing development. Comparison of transverse and longitudinal sections of three developing internodes (G3T, G3M and G3B internodes) in the same developmental stage (at the G3M stage) indicated that the number of nuclei in the tissue increased continuously from basal to apical internodes, while cell length showed the opposite trend (Figure1). Internode elongation is usually correlated with cell division and elongation. The gradually elongated cell length and gradually decreasing numbers of cell nuclei showed that cell division and elongation occurred simultaneously and thus affecting internode elongation: the former was predominant at initial stages, with the latter predominant at late stages. Thus, the culm showed sequential development from basal to apical internodes throughout the growth period.

Figure 1
figure 1

Variations in cell length (μm), occurrence rate of cell nuclei (%), four endogenous hormones concentrations forP. heterocyclaculm at different developmental stages and internodes.Abscisic acid (ABA: ng/g), indole acetic acid (IAA: μg/g), gibberellic acid (GA3: ng/g) and zeatin riboside (ZR: ng/g). Four developmental stages or three developmental internodes were defined according to the different unearthed heights of individual plants (0.05, 1.00, 6.00 or 12.00 m; nominated as G1M–G4M, respectively) or different portions (basal internode, G3B; middle internode, G3M; and top internode, G3T) of the same culm (at the G3M stage), respectively. Each data point represents the mean of three biological replicates and three experimental replicates. Different letters on a column with the same pattern indicate significant differences atP ≤ 0.05 according to the LSD test. Bars represent SE (n > 150 in cell length and occurrence rate of cell nuclei; n = 9 in four endogenous hormones concentrations).

Full size image

The dynamic changes in indole-3-acetic acid (IAA), gibberellic acid (GA3), abscisic acid (ABA) and cytokinin zeatin riboside (ZR) concentrations in four different developmental stages all displayed a unimodal type. The concentration of IAA, GA3and ZR initially gradually increased and then gradually decreased after the G3M stage. In contrast, the ABA concentration decreased initially and then gradually increased after the G2M stage. However, the changes of GA3, ZR and ABA concentrations in three different developmental internodes were basically inconsistent with the above mentioned results, unlike IAA (Figure1). Substantial changes in cellular activities and hormone levels involved in these growth stages prompted us to investigate the molecular processes at the transcriptome level.

Analysis of transcriptome sequencing

Trinity assembly of paired-end sequencing reads

RNA of different types of culms (including four developmental stages and three portions) was pooled to provide a broad gene library associated with internode elongation. A total of 45.9 million reads with length of 180 bp were obtained in the transcriptome sequencing and 154,903 contigs with an average length of 259 bp were obtained. Throughtrinityassembly, we obtained 60,393 unigenes with an average length of 612 bp and a total length of 36.96 Mbp. Of these, we gained 48,888 unigenes (80.95%) with length 200–900 bp, 9248 unigenes (15.31%) with length 1000–2000 bp, and 2257 unigenes (3.74%) with length > 2000 bp (Table1; Additional file1: Figure S1). Sequencing reads were realigned to the unigenes using Trinity [35], allowing up to 2 bp mismatches: about 58.82, 18.16 or 2.01% of the unigenes were realigned by more than 10, 100 or 1000 reads, respectively (Additional file1: Figure S1).

Table 1Length distribution of contigs, scaffolds and unigenes
Full size table

Functional classification by gene ontology (GO) and cluster of orthologous groups (COG)

Out of 60,393 unigenes, 42,127 (69.75%) showed significant similarity to known proteins in the Nr database and matched 27,408 unique protein accessions. However, 24,647 (40.81%) had BLAST hits in the Swiss-Prot database and matched 7,099 unique protein accessions. Altogether, a total of 27,408 unique protein accessions were identified in BLAST searches. This indicated that our study produced a large fraction of Moso bamboo genes. Analysis of 60,393 non-redundant unigenes using BLASTX and ESTscan software revealed that about 42,483 had reliable coding sequences, and the length of most was > 200 aa when translated into proteins. Comparison with the Nr and Swiss-Prot databases revealed that 42,280 unigenes had good comparability with known gene sequences in existing species.

转让条款32064 unigenes透露209,737 functional terms. Of them, assignments to the cellular component were the majority (91,142, 43.46%), followed by biological processes (79,250, 37.79%) and molecular function (39,345, 18.76%). Under the cellular component category, cell (23,125 unigenes, 25.37%) and organelles (19,280 unigenes, 21.15%) were notably represented. Under the category of biological processes, cellular processes (17,170 unigenes, 21.67%) and metabolic processes (16,813, 21.22%) were the majority. In the molecular function category, 17,710 and 15,853 unigenes were linked to binding and catalytic activity, respectively; it is worthy of mention that a few genes were assigned to nutrient reservoir activity, metallochaperone activity, channel regulator activity and protein tag (Additional file2: Figure S2). These indicated that some important cell activities, biological, cellular and metabolic processes occurred in culms of developing Moso bamboo shoot.

Of 42,127 Nr hits, 13,957 unigenes were assigned to the COG classifications. Of all 25 COG categories, the largest group was the cluster for general function prediction only (4720, 33.82%); followed by transcription (3300, 23.64%); replication, recombination and repair (2694, 19.30%); posttranslational modification, protein turnover and chaperones (2341, 16.77%); signal transduction mechanisms (2340, 16.77%); translation, ribosomal structure and biogenesis (2284, 16.36%); and carbohydrate transport and metabolism (1975, 14.15%). Only a few unigenes were assigned to nuclear structure and extracellular structure. Furthermore, 614 sequences were assigned to secondary metabolite biosynthesis, transport and catabolism (Additional file2: Figure S2).

Functional classification by Kyoto encyclopedia of genes and genomes (KEGG)

The KEGG Pathway database can help to further determine biological functions and interactions of genes. Based on a comparison against the KEGG database using BLASTx withE < 10–5, of the 60,393 unigenes, 23,897 (39.57%) had significant matches and were assigned to 128 KEGG pathways (Additional file3: Table S1). Among them, 5,826 (24.38%) unigenes with enzyme commission (EC) numbers were assigned to the metabolic pathways group. Of the unigenes, 1,982 (8.29%), 1,684 (7.05%), 1,652 (6.91%) and 1,627 (6.81%) were linked to the biosynthesis of secondary metabolites, endocytosis, RNA transport and glycerophospholipid metabolism, respectively.

Digital gene expression (DGE) analysis in rapidly growing culms of bamboo

Tag identification and quantification

A total of 4.2 million raw tags of the mRNAs extracted from Moso bamboo culms were identified by base calling [36,37]。Approximately 0.12, 0.11, 0.14, 0.12, 0.14 and 0.12 million high quality non-redundant clean tags were obtained in G1M–G4M, G3T and G3B groups (Additional file4: Figure S3), respectively. Gene annotation was performed by tag mapping analysis using the 60,393 non-redundant unigenes from RNA-seq-based transcriptome analysis as the reference transcript database. Results showed that 45.34, 50.76, 45.17, 46.12, 45.04 and 52.42% of all distinct tags within different group could be mapped to the entire reference database (sense or anti-sense), respectively (Additional file5: Table S2). The detected sense and anti-sense strands, sense or anti-sense strands mapped by the tags among groups were different (Additional file4: Figure S3). Among the detectable expressed unigenes, 20,849 had successful annotations.

Depth of sampling

Saturation analysis can be used to check whether the number of detected genes continues to increase when sequencing amount (total tag number) increases. When the sequencing amount reached 4 M, the number of detected genes almost ceased increasing (Additional file6: Figure S4) in G1M–G4M, G3T and G3B libraries.

Comparison of gene expression level between libraries

The differences in tag frequencies in the G1M–G4M, G3T and G3B libraries were used to estimate gene expression levels at different development stages or different internodes. The black and white histograms represent transcripts up-regulated or down-regulated in number by more than twofold in every library, respectively. For developmental stages, the highest amount of differential genes (2,410) was detected between G1M and G3M, and the least (1,262) between G1M and G2M. For developmental positions, the highest were observed between G3T and G3B (2,174), followed by G3M vs G3B (1,692), and G3T vs G3M (1,622). The number of down-regulated genes was more than that of up-regulated genes with the exception of G4M vs G3M (Figure2A). Moreover, 5837 genes that met our threshold criteria for at least one of the 15 DGE comparisons were clustered using principal coordinate analysis (PCA) and hierarchical clustering. Six samples were separated into two groups, one was G1M, G2M and G3T group, and another was G3M, G4M and G3B group (Figure2B). The result was consistent with that of hierarchical clustering (Additional file7:图S5)。

Figure 2
figure 2

The histogram showing the number of differentially expressed genes (A) and miRNAs (B) between libraries. Crepresents the length distribution of the small RNA genes in different libraries.Drepresents Principal Coordinates Analysis (PCA) of differentially expressed genes, miRNAs and proteins in different libraries. PCA indicated that six samples were divided into two groups, one was G1M, G2M and G3T, the other was G3M, G4M and G3B. G1M, G2M, G3M and G4M represent four developmental stages in turn. G3T, G3M and G3B represent top, middle and basal internode of the third developmental stage, respectively. The red circle represent mRNAs, green triangle represent miRNAs, and blue square represent proteins.

Full size image

Expression profiles of mRNAs during culms development

The 5,837 genes were further clustered using Short Time-series Expression Miner (STEM) [38]。It identified 11 (1,024 genes) or 6 (923 genes) significant gene profiles (P-Value ≤ 0.01) in development or internode, respectively (Figure3A and3C; Table2; Additional file8: Table S3).

Figure 3
figure 3

Short Time-series Expression Miner clusters of expression profiles.The profiles ordered based on the p-value significance of number (at bottom-left corner) of genes or miRNAs assigned versus expected. Colour square frame denote significant profiles (P-value ≤ 0.01). Each graph displays the mean pattern of expression (black lines) of the profile genes. The number of profiles in each cluster is at the top left corner of each STEM.AandBrepresent the development-specific mRNA and miRNA profiles, respectively.CandDrepresent the internode-specific mRNA and miRNA profiles, respectively. The x-axis represents time (development-specific or internode-specific stages) and they-axis represents log2fold change of gene expression.

Full size image
Table 2Genes and miRNAs grouped into significant cluster profiles by STEM analysis
Full size table

The development-specific gene clusters, STEM profile 9 (177 genes) and profile 42 (177 genes) showed respectively down and up regulated along with the culm development and consistent with the decrease of cell nuclei number and prolong of cell length (cell division) along with clum development and aging. STEM profile 18 (110 genes) and 34 (164 genes) have a reversal in the trend of expression change between development and aging, and profile 18 was consistent with the change of ABA content. Profile 21 (157 genes) and Profile 31 (179 genes) gene clusters have a reversal in G2M stage and no change at later stage, while profile 22 (251 genes) and Profile 30 (148 genes) gene clusters have a reversal in G3M stage and no change at earlier stage. Profile 23 (371 genes) and Profile 29 (298 genes) show up- or down-regulated during the middle stage, and no change at the earlier and later stages. Contrary to profile 23 and 29, Profile 40 contained 95 genes (Figure3A; Additional file8: Table S3).

The internode-specific gene clusters, Profile 2 (171 genes) and 7 (304 genes) represented down-regulated genes, were consistent with the decrease of cell nuclei (cell division) along with internode development. Profile 11 (284 genes), 12 (226 genes), 13 (237 genes) and 15 (239 genes), up-regulation internode-specific genes, were consistent with the prolong of cell length (cell division) along with internode development (Figure3C; Additional file8: Table S3).

miRNAs regulate expression changes in rapidly growing culms of bamboo

测序和数据分析

探讨小rna的microrna的组成部分nd the dynamic changes of the miRNAs during the rapid growth of developing culms in Moso bamboo, we purified the culms of G1M-G4M, G3T and G3B from bamboo and sequenced their small RNAs using Solexa high-throughput technology. After removing sequences of low quality, 5′ adaptor contaminants, 3′ adaptor null, insert null, poly A and RNAs smaller than 18nt, we obtained 12,724,584, 14,638,070, 13,279,680, 12,065,802, 12,841,021 and 12,360,293 clean reads, respectively (Additional file9: Table S4). These high-quality small RNAs were used for further analysis.

数以百万计的高质量的小分子rna从the individual libraries, 72.21–95.81% was 20 to 24 nucleotides in length, which is the typical size range for Dicer-derived products. The major component of small RNAs in all libraries was 24 nucleotides long. Throughout clum development, the proportion of 24-nucleotide small RNAs decreased and the 21-nucleotide population increased in G1M-G4M, while G3M contained mostly 21-nucleotide small RNAs (Figure2D).

Prediction of novel miRNA and comparison of miRNA expression level between libraries

In total, we obtained 732 miRNAs and 453 predicted novel miRNAs. Of these novel miRNAs, 26, 61, 101, 57, were expressed in G1M-G4M, and 92 in G3T, 116 in G3B, respectively (Additional file10: Table S5). Differ from the result of differentially expressed genes; the number of differential miRNAs in the initial stages was more than that in the descendent stages, the amount of up-regulated miRNAs was more than that of down-regulated miRNAs with the exception of G1M vs G2M (Figure2B). However, PCA and hierarchical clustering showed the consistency between miRNAs (165 miRNAs) and mRNAs (5837 genes), that is, six samples were separated into two groups, one was G1M, G2M and G3T group, another was G3M, G4M and G3B group, just as protein (213 proteins) clustering result (Figure2C; Additional file7:图S5)。此外,相比之下,microrna的经验ression, higher consistency between mRNAs and proteins was observed.

Expression profiles of miRNAs during culms development

To directly compare the expression patterns of these miRNAs in the different development stages or internodes, four development-specific profiles (54 miRNAs) or one internode-specific profile (101 miRNAs) were identified significantly by STEM cluster, respectively (Figure3B and3D and 3D; Table2; Additional file8: Table S3). Profile 4 (down-regulation, 27 miRNAs) and 25 (up-regulation, 24 miRNAs) continue the development trend. Profile 17 (28 miRNAs) and 18 (25 miRNAs) have a reversal in the trend of expression change between development (Figure3B). The internode-specific significant miRNAs cluster, profile 4 (down-regulation) contained 130 miRNAs (Figure3D).

The integrated analysis between STEM and positive/negative correlation

In general, the mature miRNA negatively regulates their target mRNAs by either translational repression or mRNA degradation[39]。Through target prediction using the Miranda web service, we found 6598 targets for 671 known-miRNAs and 4616 targets for 298 novel-miRNAs (Additional file11: Table S6). Analysis to identify negative correlations between miRNA and mRNA expression was done using an in house R script, and pearson correlation coefficient (cor) was computed for each miRNA and their predicted target mRNAs. Forty-five of 121 miRNA-mRNA pairs combined by 37 unique miRNAs and their 44 target mRNAs presented an anti-correlationship between the expression level of miRNAs and mRNAs (cor < −0.81 and P-value of correlation ≤0.05) (Additional file12: Table S7). Positive correlation between 2-DE staining intensity and gene expression levels determined by RNA-Seq for fast-growing clum in Moso bamboo are summarized in Additional file12: Table S7. Statistical pearson correlation ranged from 0.53–0.93, and eight of 25 mRNA-protein pairs combined by 24 unique corresponding proteins presented a plus correlationship (cor > 0.81 and P-value of correlation ≤0.05).

Based on the analysis of intersection/integration between significant STEM and positive/negative correlation, seventy three pairs were obtained, including 64 genes, 55 miRNAs (34 known- and 21 novel-miRNAs) and 15 proteins. Of the 73 pairs, 34 pairs (27 genes) were significant correlation (P-value of correlation ≤0.05) (Additional file12: Table S7). According to interaction of 64 genes and correlation between miRNA-mRNA and mRNA-protein, a complex network was built in Figure4, 64个基因在核心位置(大喊ow circle). The GO annotation of 64 genes involved in grow and reproduction (especially xylem or phloem histogenesis and formation, leaf morphogenesis), auxin and ABA mediated signaling, cellular process (ribosome biogenesis and cell division or death), metabolic process (protein ubiquitination, photorespiration, glycolysis, proteolysis, glycan and energy metabolism, especially hemicellulose metabolic process), biological regulation (such as regulation of ARF GTPase activity, regulation of anthocyanin biosynthetic process), and response to stimulus (response to cold or cadmium ion). The pathways included metabolism (carbohydrate metabolism, lipid metabolism, nucleotide metabolism, amino acid metabolism, energy metabolism, glycan biosynthesis and metabolism, metabolism of terpenoids and polyketides, biosynthesis of other secondary metabolites), genetic information processing (translation), cellular processes (transport and catabolism), environmental information processing (signal transduction, membrane transport) and environmental adaptation (plant circadian rhythm) (Additional file12: Table S7).

Figure 4
figure 4

A network diagram based on gene interaction, STEM cluster and correlation between mRNA-protein and miRNA-mRNA (left).The negative correlation analysis depend on miRNAs and microRNA-regulated target genes, genes and their encode proteins carry out positive analysis. Gene-gene interaction was compared with rice gene interactions. Red circle represent genes, yellow circle represent 64 target genes, 27 statistically significant genes were labeled in bigger yellow circles, blue square represent proteins, green triangle represent miRNAs. Gene Ontology classification of assembled unigenes (right). Fifty of 64 unigenes with BLAST matches to known proteins were assigned to three main categories: biological process, cellular component and molecular function.

Full size image

Discussion

Here, we report on the first next-generation sequencing study of the temporal and spatial mRNA and miRNA expression profiling of bamboo internodes. A DGE approach enabled identifying those mechanisms related to rapid growth of developing culms in Moso bamboo. STEM, positive and negative correlation analysis were perfectly applied in integrated analysis of multi-omics data. These provided a molecular basis underlying the phenomenon of sequentially elongated internodes from the base to the top.

Temporal and spatial growth patterns of internodes in moso bamboo

The initial, ascending, boosting and terminal stages are involved in growth phases in bamboo culms. Sequentially elongated internodes from the base to the top in bamboo have been reported [14]。In our previous study, the growth stages of different internodes within the same culm were established using morphological observations, endogenous hormones and anatomical structure analysis [29]。However, the molecular basis underlying this phenomenon remains unclear. To determine transcriptional and post-transcriptional regulation changes that occurred during different growth stages of bamboo culms, we analyzed six internodes from different stages (G1M–G4M) or different portions (G3B, G3M and G3T) of the same culm (at the G3M stage). It is generally recognized that different numbers of differentially expressed genes could represent or describe different growth degrees, i.e. the greater the number between stages, the larger the difference between them. Comparisons of gene expression level between the six DGE libraries indicated that the relationship between G1M and G2M was the closest and the farthest between G2M and G3M for four developmental stages. However, for the three spatial positions the relationship between G3T and G3B was the farthest, followed by G3M and G3B, and then G3T and G3M. In addition, PCA and hierarchical clustering analysis based on differentially expressed genes showed that G1M, G2M and G3T were classified in the one group; whereas G3M, G4M and G3B were classified in another group, and cells of the latter were longer than the former. This was consistent with the clustering result from differentially expressed proteins and similar to differentially expressed miRNAs, indicating the consistency in the differentially expressed levels of transcripts, miRNAs and proteins among stages. All our data suggested that G1M, G2M and G3T were involved in the initial or ascending stages; and G3M, G4M and G3B were involved in the boosting or descendent stages. This result was not surprising, because although the whole plants were defined as a specific growth stage (usually defined according to the middle internode), different portions of the same plant may be at different stages resulting from the sequential basal–apical elongation of the culm internodes[14]。When the growth of the whole plant was in boosting stage (G3M), the growth of the base internode (G3B) might be in descendent stage, while the top internode (G3T) might be in the initial or ascending stage. Thus, our study provided a molecular basis underlying the phenomenon of sequentially elongated internodes from the base to the top.

GO and pathways involved in rapidly elongating culms ofP. Heterocycla

Short Time-series expression data/experiments provide a wealth of information regarding the global view of the dynamic networks that are activated in dynamic biological processes [38,40]。Many studies show that a comprehensive understanding of the control of gene expression will require quantitative information at all levels, from DNA variants through to differential stability of the products or their regulative factors [41–43]。Our temporal or static omics data, including transcriptome, post-transcriptome and translatome/proteome, play an even bigger part in molecular mechanisms study of dynamic bamboo’ culm development. By integrating such data with their positive/negative correlation analysis data sets we obtain the complete set of gene/protein or miRNA that are activated, their dynamics and interactions, the role that different gene/protein or miRNA have in the process and the differences in processes within and between different growth stages. At last the core genes focused on 64 genes and their correlative 55 miRNAs and 15 proteins were obtained, the GO and pathway analysis revealed that genes involved in response to stimulus, cell cycle, cellular component, regulation of plant hormone levels, many metabolism and signal transduction played major roles in persistently elongating culms ofP. heterocycla. Meanwhile, with the development of bamboo culm, gradually declined cell nuclei, gradually elongated parenchyma and fiber cells, dynamic changes of a unimodal-type endogenous hormones concentration were observed in our study. According to these, we proposed some scenarios involved in culm elongation (Figure5). During the growth season, some environmental cues such as sufficient moisture, appropriate temperature, fertile soil and low light drive plant circadian rhythm, stimulate cell differentiation, division and growth, and biosynthesis of plant hormones, especially auxin and ABA, through some signaling sensor and transduction pathways. Furthermore, the cells take up energy or nutrients, possibly supplied by the mother bamboo rhizome and root system, in the cellular path through glycolysis, photorespiration, glycan and energy metabolism to maintain normal metabolic activity, which promotes cell division, elongation or death, lignin and cellulose deposition, eventually leading to growth of the culm.

Figure 5
figure 5

A scenario of bamboo’ culm elongation and development at molecular, physiological, cellular and tissue level.During the growth season, some environmental cues such as sufficient moisture, appropriate temperature, fertile soil and low light drive plant circadian rhythm, stimulate cell differentiation, division and growth, and biosynthesis of plant hormones, especially auxin and ABA, through some signaling sensor and transduction pathways. And then the cells take up energy or nutrients through glycolysis, photorespiration, glycan and energy metabolism to maintain normal metabolic activity, which promotes cell division, elongation or death, lignin and cellulose deposition, eventually leading to elongation and development of the Moso bamboo’ culm.

Full size image

Stimulus response of an organism is usually due to some external condition or event. Many genes/proteins involved in response to stimulus have been reported in the developing maize rachis [44] and fruit development of sweet orange [45]。Seventeen unigenes (Additional file12: Table S7) involved in response to stimulus, such as carbohydrate, cytokinin, auxin or ABA stimulus, cold, blue light, cadmium ion, gravitropism were detected during the fast-growing stage of bamboo, indicating the important role of environmental cues in bamboos’ rapid growth. Signal transduction occurs when an extracellular signaling molecule activates a cell surface receptor. In turn, this receptor will alter intracellular molecules with a response. The hormone-mediated signaling pathway was one signal transduction found to be regulated in fast-growing culms of bamboo. It comprised auxin mediated signaling pathway (one unigene) and abscisic acid mediated signaling pathway (one unigene). The data indicated some genes in diverse hormone (auxin, ABA, ethylene and GAs) signaling pathways are regulated in different trends or stages in rice seed development [46]。The phosphatidylinositol signaling system was also found to be regulated in fast-growing culms of bamboo. It comprised inositol or phosphatidylinositol kinase activity (one unigenes) and nucleoside-triphosphatase activity (five unigenes). Some studies also illustrated that the phosphatidylinositol signaling system played an important role in plant development, such as storage root formation [47], maize growth [48] and rice seedling germination [49]。Our study will help to expand knowledge of the function of the phosphatidylinositol signaling system in plant growth.

Plant hormones are major regulators of plant growth and development. Endogenous hormones as well as environmental factors can regulate progression through the cell cycle: shaping the plant and affecting seed growth, time of flowering, sex of flowers, and senescence of leaves and fruit [50,51]。Auxins, the first identified plant hormones (and commonly IAA), are compounds that positively influence cell enlargement, bud formation and root initiation. They affect cell elongation by altering cell wall plasticity. They stimulate the cambium (a subtype of meristem cells) to divide, and in stems cause secondary xylem to differentiate. Some studies have shown that auxin promotes proliferation activities in the cambium during secondary growth. For example, it was indicated that auxin maxima are established and genes associated with auxin signaling are up-regulated in the actively proliferating cambial cells in hybrid aspen [52]。By contrast, general reduction in polar auxin transport and down-regulation of auxin signaling genes were observed during the phase of cambium dormancy [52,53]。Consistently, disruption in auxin signaling/responsiveness led to detrimental effects on cambial activities, reducing secondary growth in hybrid aspen plants [54,55]。Similarly, several other plant hormones, in particular cytokinin, ethylene and gibberellins are known to regulate cambial development [56,57]。In our studies, four major endogenous hormones (IAA, GA3, ZR and ABA) appeared to strongly influence the cell elongation phase. In the present study, four unigenes involved in auxin efflux and abscisic acid mediated signaling pathway, ethylene biosynthetic process, response to cytokinin stimulus, were identified (Additional file12: Table S7). This indicates that in the bamboo growth season, environmental cues might stimulate cell division, growth or death through regulating hormone transport of influencing the expression of hormone response factors.

A cellular component may be an anatomical structure (e.g. rough endoplasmic reticulum or nucleus) or a gene product group (e.g. ribosome, proteasome or protein dimer). During secondary vascular development in poplar, the majority of the integral proteins identified were plasma membrane proteins [58]。Some studies have suggested that up-regulation of plasma membrane aquaporins improves the photosynthetic activity and growth of Eucalyptus trees [59]。In our study, numerous cell nuclei were found in parenchyma and fiber cells – an obvious sign of cell division – indicating the presence of meristematic tissue. As the culm developed, the number of nuclei declined, until there were almost no detectable nuclei during the late stages of development. Meanwhile, the length of parenchyma and fiber cells increased continuously during development [29]。Many genes involved with the plasma membrane (12 unigenes), organelles (38 unigenes) such as the nucleolus (two unigenes), organelle envelope (six unigenes), chloroplast thylakoid membrane (two unigenes), cell wall (three unigenes), and Golgi membrane (one unigenes), cytosol (10 genes) and apoplast (five genes) influenced the fast growth of the culm (Additional file12: Table S7). Most of them peaked at the ascending or boosting stages. This suggests that cell division is essential to rapid elongation of the culms during their early developmental stages. Remarkably, six genes involving in biosynthesis of lignin and cellulose/hemicellulose, the major cell wall components, were detected, it suggests that lignocellulose is accumulated along with the development of bamboo’ culm.

代谢过程是必要的。一个麻木er of excellent studies have indicated that metabolic processes were among the most regulated during developing stages of plants, and these genes were involved in many biological processes: e.g. glycan biosynthesis and metabolism, and the metabolisms of carbohydrates, lipids, amino acids and energy [26,44,45,60,61]。In developing bamboo culm in the present study, the dominant metabolic processes (27 unigenes) included glycan biosynthesis and metabolism, the metabolisms of carbohydrates, lipids, amino acids, nucleotides, energy and biosynthesis of other secondary metabolites (Additional file12: Table S7). These metabolic processes can provide the energy and components for DNA replication, signal transduction, cellular growth, hormone biosynthesis, cellulose and lignin biosynthesis, which are essential for the rapid culm elongation under conditions of light deficiency.

Conclusions

To the best of our knowledge, this study is the first exploration to characterize the temporal and spatial transcriptome, miRNA and mRNA expression profiling in a developing bamboo culm. We provided a molecular basis underlying the phenomenon of sequentially elongated internodes from the base to the top. Several key pathways such as environmental adaptation, signal transduction, translation, transport and many metabolism, and GO terms such as hormone-mediated signaling, cell growth, division and differentiation, primary shoot apical meristem specification, xylem and phloem pattern formation, protein modification, response to stimuli, metabolic process and biological regulation were identified as involved in rapid growth of bamboo culms. In addition to gaining more insight into the unique growth characteristics of bamboo, we provide a good case study to analyze mRNA and miRNA expression of a non-model plant species using high-throughput short-read sequencing. Cluster concordance of protein, transcript and post-transcriptional levels are obvious. Above all, we demonstrate that the integrated analysis of our multi-omics data, including transcriptome, post-transcriptome/miRNAs and translatome/proteome, yield more complete representations and additional biological insights, especially the complex dynamic processes occurring in the rapid growth of developing culms in Moso bamboo.

Methods

Plant materials

Culm samples from a natural population ofP. heterocyclawere harvested in spring 2009 at the Dagangshan Forest Ecosystem Research Station in Jiangxi Province of China (27°30′–27°50′N, 114°30′–114°45′E), whereP. heterocyclais typically distributed in the subtropical zone. To ensure that there were absolute minimal differences in the growth conditions, all culms were sampled within the distribution of 50 × 50 m2. To ensure that equivalent staged samples and developmental stages were taken as replicates, individual plants were sampled with the same height at each stage. Based on previous observations [14,29] and our hormonal and anatomical results, four developmental stages or three developmental internodes were defined according to the different unearthed heights of individual plants (0.05, 1.00, 6.00 or 12.00 m; nominated as G1M–G4M, respectively) or different portions (basal internode, G3B; middle internode, G3M; and top internode, G3T) of the same culm (at the G3M stage), respectively. The samples, each with three biological replicates (three trees), were cut from tissue inside the epidermis located in the basal part of each middle internode for different stages or for different internodes at the same stage as mentioned above. The samples were immediately frozen in liquid nitrogen and kept at −80°C until analysis.

Light microscopy and ELISA

Approximately 0.5-cm3samples of culm tissue were fixed in FAA (Formalin-Acetic-Alcohol) buffer and exhausted with an aspirator pump. Subsequently, serial transverse and longitudinal sections (8 μm thick) from paraffin-embedded tissue were sequentially stained with safranin and fast green. Finally, these sections were observed with a Zeiss Axiophot light microscope (Zeiss, Jena, Germany). Digital images were analyzed using Zeiss Axio Vision 4.6.

The method for extraction, purification and quantification of endogenous plant hormones IAA, GA3, ABA and ZR were modified from the description of Wang [62]。激素的ELISA试剂盒用于估计levels came from China Agricultural University (Beijing, China). One-way ANOVA and least significant difference (LSD) tests were used to check for differences in cell morphology and hormone content.P ≤ 0.05 was considered statistically significant.

Sample preparation for RNA-seq and transcriptome analysis

组织样本来自18个茎(6开发鹿es or internodes, each with three biological replicates) were collected for RNA preparation. Total RNA was isolated using a TRIzol reagent (Gibco BRL) following the manufacturer’s instructions, and then beads with Oligo (dT) were used to isolate poly (A) mRNA. After fragmentation, random hexamer-primer was used to synthesize the cDNA. Short cDNA fragments were purified and followed by end reparation, adenylation and sequencing adapter ligation. After PCR amplification of the selected fragment, one pooling library was sequenced using Illumina HiSeq™ 2000. The raw sequencing datasets of transcriptome were deposited in the Sequence Read Archive (SRA) of NCBI (http://www.ncbi.nlm.nih.gov/sra) under accession number SRX329521.

Transcriptomede novoassembly was carried out with the short reads assembling program – Trinity [63] : it first combined reads with certain length of overlap to form contigs, and then reads are mapped back to contigs; with paired-end reads, it is able to detect contigs from the same transcript as well as the distances between these contigs. Finally, Trinity connects the contigs, and gets sequences that cannot be extended on either end. Such sequences are defined as Unigenes. After clustering, the unigenes will be divided to two classes. One is clusters with the prefix CL, and another is singletons with the prefix Unigene. Finally, BLASTx alignment (E < 10–5) between unigenes and protein databases such as non-redundant protein (Nr), Nucleotide (Nt), Swiss-Prot and KEGG, then COG was performed, and the best aligning results used to determine sequence direction of unigenes. When a unigene was not aligned in any of the above databases, ESTScan [64] was used to predict its coding regions and sequence direction. On the basis of Nr annotation, the Blast2GO and the WEGO program were used to perform GO functional classification. COG and the KEGG Pathway database were also used to further determine the biological functions and interactions of unigenes.

DGE-tag profiling

DGE sample preparation was performed using Illumina Gene Expression Sample Prep Kit. Of total RNA, 6 μg was purified using Oligo (dT) magnetic beads, and then Oligo (dT) was used as a primer to synthesize the first and second-strand cDNA. The 5′-ends of tags were generated byNlaIII. The fragments apart from the 3′-cDNA fragments connected to Oligo (dT) beads were washed away and the Illumina adaptor 1 was ligated to the sticky 5′-end of the digested bead-bound cDNA fragments. The junction of Illumina adaptor 1 and CATG site was cut byMmeI. After removing 3′-fragments by magnetic bead precipitation, Illumina adaptor 2 was ligated to the 3′-ends of tags, acquiring tags with different adaptors of both ends to form a tag library. After 15 cycles of linear PCR amplification, purification by 6% TBE PAGE gel electrophoresis, and denaturation, the single-chain molecules were fixed onto the Illumina Sequencing Chip (flowcell). After cluster generation, sequencing with the method of sequencing by synthesis was performed using Illumina HiSeq™ 2000. The raw reads were submitted to the SRA of NCBI under the accession number SRX329613.

Aligning DGE tags to the reference transcriptome data set and identification of differentially expressed genes

Clean tags and count number of six DGE libraries from G1M–G4M, G3B and G3T groups were collected and summarized using custom Bioperl scripts. All tags were mapped to the reference transcriptome generated by RNA-seq. To monitor mapping events on both strands, both sense and complementary anti-sense sequences were included in the mapping process. Only perfect matches over the entire 21-bp length of the 17-bp tag plus the 4-bpNlaIII recognition site were allowed. This study was limited to tags that mapped to ORFs (open reading frames) only and cannot show tags that mapped to mRNA with long 3′-UTRs (untranslated regions). Rigorous algorithms were developed to identify differentially expressed genes between two groups [65]。The correlation of the detected count numbers between parallel libraries was statistically assessed using Pearson’s correlation. In addition to theP-value, false discovery rate (FDR) was used in DGE to determine the threshold ofP-value in multiple test and analysis [66]。In this research,P-value ≤ 0.01, FDR ≤ 0.001, and the absolute value of log2Fold-change ≥1 or ≤ −1 were used as threshold to judge the significance of gene expression differences.

Small RNAs preparation and analysis

The experiment process of miRNA was described in Hafner et al. [67]。Briefly, low molecular-weight RNA was isolated and sequenced following the Small RNA Sample Preparation Protocol (Illumina). The 50nt sequence tags from HiSeq sequencing will go through the data cleaning first, which includes getting rid of the low quality tags and several kinds of contaminants. Length distribution of clean tags is then summarized. The raw reads are available in the NCBI SRA database under the accession number SRX330460. Afterwards, the standard bioinformatics analysis will annotate the clean tags into different categories and take those which cannot be annotated to any category to predict the novel miRNA. The characteristic hairpin structure of miRNA precursor was used to predict novel miRNA and software Mireap (http://sourceforge.net/projects/mireap/files/mireap/) was developed to predict novel miRNA by exploring the secondary structure, the Dicer cleavage site and the minimum free energy of the unannotated small RNA tags which could be mapped to genome.

Correlation analysis

Analysis to identify negative correlations between miRNA and mRNA expression was done using an inhouse R script. Briefly, normalized miRNA and mRNA data were sample-matched for all samples with both miRNA and mRNA sequencing data. Then for each miRNA, Pearson correlation coefficients were computed for its predicted target mRNAs, and a contingency table was created for all mRNAs, which was used to assess the enrichment level of the negative correlated mRNAs (correlation < 0 and P-value of correlation ≤ 0.05) within predicted targets of the intended miRNA using Fisher’s exact test. The method of positive correlation analysis between mRNAs and corresponding proteins was the same as that of negative correlations.

References

  1. 法雷尔AE、Plevin RJ特纳BT,琼斯广告,奥黑尔M, Kammen DM: Ethanol can contribute to energy and environmental goals. Science. 2006, 311: 506-10.1126/science.1121416.

    PubMedCASArticleGoogle Scholar

  2. Salas Fernandez MG, Becraft PW, Yin Y, Lübberstedt T: From dwarves to giants? plant height manipulation for biomass yield. Trends plant sci. 2009, 14: 454-461. 10.1016/j.tplants.2009.06.005.

    PubMedCASArticleGoogle Scholar

  3. Magel E, Kruse S, Lütje G, Liese W: Soluble carbohydrates and acid invertases involved in the rapid growth of developing culms in Sasa palmata (Bean) Camus. Bamboo Science and Culture. 2005, 19: 23-29.

    Google Scholar

  4. Zheng YS, Hong W, Qiu EF, Chen LG: Hormone content and distribution inPhyllostachys heterocyclacv.Pubescensduring period of shoot emergence. Scientia Silvae Sinicae. 1998, 34: 100-104.

    Google Scholar

  5. He XQ, Suzuki K, Kitamura S, Lin JX, Cui KM, Itoh T: Toward understanding the different function of two types of parenchyma cells in bamboo culms. Plant cell physiol. 2002, 43: 186-195. 10.1093/pcp/pcf027.

    PubMedCASArticleGoogle Scholar

  6. Li XB, Shupe TF, Peter GF, Hse CY, Eberhardt TL: Chemical changes with maturation of the bamboo species phyllostachys pubescens. J Trop For Sci. 2007, 19: 6-

    Google Scholar

  7. Yu Y, Tian G, Wang H, Fei B, Wang G: Mechanical characterization of single bamboo fibers with nanoindentation and microtensile technique. Holzforschung. 2011, 65: 113-119.

    CASGoogle Scholar

  8. 日元TM,李JS:比较碳seque地上stration between moso bamboo (Phyllostachys heterocycla) and china fir (Cunninghamia lanceolata) forests based on the allometric model. For Ecol Manage. 2011, 261: 995-1002. 10.1016/j.foreco.2010.12.015.

    ArticleGoogle Scholar

  9. Wang XQ, Ren HQ, Zhang B, Fei BH, Burgert I: Cell wall structure and formation of maturing fibres of moso bamboo (Phyllostachys pubescens) increase buckling resistance. J R Soc Interface. 2011, 9: 988-996.

    PubMedPubMed CentralArticleGoogle Scholar

  10. 温家宝GS,张LY,张RM,曹,周通用、黄H, Wong MH: Temporal and spatial dynamics of carbon fixation by moso bamboo (Phyllostachys pubescens) in Subtropical China. Bot Rev. 2011, 77: 271-277. 10.1007/s12229-011-9068-x.

    PubMedPubMed CentralArticleGoogle Scholar

  11. Lee C, Chin T: Comparative anatomical studies of some Chinese bamboos. Acta Bot Sin. 1960, 9: 76-95.

    Google Scholar

  12. Murphy R, Alvin K: Variation in fibre wall structure in bamboo. IAWA Bull(NS). 1992, 13: 403-410.

    ArticleGoogle Scholar

  13. Lin J, He X, Hu Y, Kuang T, Ceulemans R: Lignification and lignin heterogeneity for various age classes of bamboo (Phyllostachys pubescens) stems. Physiol plantarum. 2002, 114: 296-302. 10.1034/j.1399-3054.2002.1140216.x.

    CASArticleGoogle Scholar

  14. Jiang ZH: Bamboo and rattan in the world. Shen Yang, China: LiaoNing Science and Technology Bublishing House: 2002.

    Google Scholar

  15. Luo A, Qian Q, Yin H, Liu X, Yin C, Lan Y, Tang J, Tang Z, Cao S, Wang X: EUI1, Encoding a putative cytochrome P450 monooxygenase, regulates internode elongation by modulating gibberellin responses in rice. Plant cell physiol. 2006, 47: 181-

    PubMedCASArticleGoogle Scholar

  16. Zhou HL, He SJ, Cao YR, Chen T, Du BX, Chu CC, Zhang JS, Chen SY: OsGLU1, A putative membrane-bound endo-1, 4-β-D-glucanase from rice, affects plant internode elongation. Plant mol biol. 2006, 60: 137-151. 10.1007/s11103-005-2972-x.

    PubMedCASArticleGoogle Scholar

  17. Ruonala R, Rinne PLH, Kangasjärvi J, Van Der Schoot C: CENL1 Expression in the rib meristem affects stem elongation and the transition to dormancy in populus. Plant Cell. 2008, 20: 59-74. 10.1105/tpc.107.056721.

    PubMedCASPubMed CentralArticleGoogle Scholar

  18. Hattori Y, Nagai K, Furukawa S, Song XJ, Kawano R, Sakakibara H, Wu J, Matsumoto T, Yoshimura A, Kitano H: The ethylene response factors SNORKEL1 and SNORKEL2 allow rice to adapt to deep water. Nature. 2009, 460: 1026-1030. 10.1038/nature08258.

    PubMedCASArticleGoogle Scholar

  19. Asano K, Miyao A, Hirochika H, Kitano H, Matsuoka M, Ashikari M: SSD1, Which encodes a plant-specific novel protein, controls plant elongation by regulating cell division in rice. Proc Jpn Acad Ser B Phys Biol Sci. 2010, 86: 265-273. 10.2183/pjab.86.265.

    PubMedCASPubMed CentralArticleGoogle Scholar

  20. Iwamoto M, BABA-ASAI A, Kiyota S, Hara N, Takano M: ACO1, a gene for aminocyclopropane - arboxylate oxidase: effects on internode elongation at the heading stage in rice. Plant Cell Environ. 2010, 33: 805-815.

    PubMedCASGoogle Scholar

  21. Chiu WB, Lin CH, Chang CJ, Hsieh MH, Wang AY: Molecular characterization and expression of four cDNAs encoding sucrose synthase from green bamboo bambusa oldhamii. New Phytol. 2006, 170: 53-63. 10.1111/j.1469-8137.2005.01638.x.

    PubMedCASArticleGoogle Scholar

  22. Hsieh CW, Liu LK, Yeh SH, Chen CF, Lin HI, Sung HY, Wang AY: Molecular cloning and functional identification of invertase isozymes from green bamboo bambusa oldhamii. J Agric Food chem. 2006, 54: 3101-3107. 10.1021/jf052711s.

    PubMedCASArticleGoogle Scholar

  23. Chen TH, Huang YC, Yang CS, Yang CC, Wang AY, Sung HY: Insights into the catalytic properties of bamboo vacuolar invertase through mutational analysis of active site residues. Phytochemistry. 2009, 70: 25-31. 10.1016/j.phytochem.2008.10.004.

    PubMedCASArticleGoogle Scholar

  24. Chen CY, Hsieh MH, Yang CC, Lin CS, Wang AY: Analysis of the cellulose synthase genes associated with primary cell wall synthesis in bambusa oldhamii. Phytochemistry. 2010, 71: 1270-1279. 10.1016/j.phytochem.2010.05.011.

    PubMedCASArticleGoogle Scholar

  25. Peng ZH, Lu TT, Li LB, Liu XH, Gao ZM, Hu T, Yang XW, Feng Q, Guan JP, Weng QJ: Genome-wide characterization of the biggest grass, bamboo, based on 10,608 putative full-length cDNA sequences. BMC Plant Biol. 2010, 10: 116-10.1186/1471-2229-10-116.

    PubMedPubMed CentralArticleGoogle Scholar

  26. Gao ZM, Li CL, Peng ZH: Generation and analysis of expressed sequence tags from a normalized cDNA library of young leaf from Ma bamboo (Dendrocalamus latiflorus Munro). Plant cell rep. 2011, 30: 2045-2057. 10.1007/s00299-011-1112-0.

    PubMedCASArticleGoogle Scholar

  27. 林XC Chow泰,陈HH,刘CC,周SJ,黄BL, Kuo CI, Wen CK, Huang LC, Fang W: Understanding bamboo flowering based on large-scale analysis of expressed sequence tags. Genet Mol Res. 2010, 9: 1085-1093. 10.4238/vol9-2gmr804.

    PubMedCASArticleGoogle Scholar

  28. 高周M,杨P, P,唐DQ:识别differentially expressed sequence tags in rapidly elongating Phyllostachys pubescens internodes by suppressive subtractive hybridization. Plant Mol Biol Rep. 2011, 29: 224-231. 10.1007/s11105-010-0222-0.

    ArticleGoogle Scholar

  29. Cui K, He CY, Zhang JG, Duan AG, Zeng YF: Temporal and spatial profiling of internode elongation-associated protein expression in rapidly growing culms of bamboo. J Proteome Res. 2012, 11: 2492-2507. 10.1021/pr2011878.

    PubMedCASArticleGoogle Scholar

  30. Peng Z, Lu Y, Li L, Zhao Q, Feng Q, Gao Z, Lu H, Hu T, Yao N, Liu K: The draft genome of the fast-growing non-timber forest species moso bamboo (Phyllostachys heterocycla). Nat Genet. 2013, 45: 456-461. 10.1038/ng.2569.

    PubMedCASArticleGoogle Scholar

  31. R: Transcriptome sequencing and de novo analysis for ma bamboo (Dendrocalamus latiflorus Munro) using the illumina platform. PloS one. 2012, 7: e46766-10.1371/journal.pone.0046766.

    PubMedCASPubMed CentralArticleGoogle Scholar

  32. Wu YJ, Chen HM, Wu TT, Wu JS, Chu RM, Juang RH: Preparation of monoclonal antibody bank against whole water-soluble proteins from rapid-rowing bamboo shoots. Proteomics. 2006, 6: 5898-5902. 10.1002/pmic.200600278.

    PubMedCASArticleGoogle Scholar

  33. Wang Z, Gerstein M, Snyder M: RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009, 10: 57-63. 10.1038/nrg2484.

    PubMedCASPubMed CentralArticleGoogle Scholar

  34. Haas BJ, Zody MC: Advancing RNA-seq analysis. Nat Biotechnol. 2010, 28: 421-423. 10.1038/nbt0510-421.

    PubMedCASArticleGoogle Scholar

  35. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q: Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature biotechnol. 2011, 29: 644-652. 10.1038/nbt.1883.

    CASArticleGoogle Scholar

  36. Ewing B, Green P: Base-calling of automated sequencer traces using Phred II. error probabilities. Genome res. 1998, 8: 186-194.

    PubMedCASArticleGoogle Scholar

  37. Ewing B, Hillier LD, Wendl MC, Green P: Base-calling of automated sequencer traces using Phred. I. Accuracy assessment. Genome res. 1998, 8: 175-185. 10.1101/gr.8.3.175.

    PubMedCASArticleGoogle Scholar

  38. Ernst J, Nau GJ, Bar-Joseph Z: Clustering short time series gene expression data. Bioinformatics. 2005, 21: i159-i168. 10.1093/bioinformatics/bti1022.

    PubMedCASArticleGoogle Scholar

  39. van Rooij E, Olson EN: MicroRNAs: powerful new regulators of heart disease and provocative therapeutic targets. J Clin Investig. 2007, 117: 2369-2376. 10.1172/JCI33099.

    PubMedCASPubMed CentralArticleGoogle Scholar

  40. Bar-Joseph Z, Gitter A, Simon I: Studying and modelling dynamic biological processes using time-series gene expression data. Nat Rev Genet. 2012, 13: 552-564. 10.1038/nrg3244.

    PubMedCASArticleGoogle Scholar

  41. Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, Chen W, Selbach M: Global quantification of mammalian gene expression control. Nature. 2011, 473: 337-342. 10.1038/nature10098.

    PubMedArticleGoogle Scholar

  42. Somel M, Guo S, Fu N, Yan Z, Hu HY, Xu Y, Yuan Y, Ning Z, Hu Y, Menzel C: MicroRNA, mRNA, and protein expression link development and aging in human and macaque brain. Genome res. 2010, 20: 1207-1218. 10.1101/gr.106849.110.

    PubMedCASPubMed CentralArticleGoogle Scholar

  43. Zhu M, Yi M, Kim CH, Deng C, Li Y, Medina D, Stephens R, Green J: Integrated miRNA and mRNA expression profiling of mouse mammary tumor models identifies miRNA signatures associated with mammary tumor lineage. Genome biol. 2011, 12: R77-10.1186/gb-2011-12-8-r77.

    PubMedCASPubMed CentralArticleGoogle Scholar

  44. Pechanova O, Pechan T, Ozkan S, McCarthy FM, Williams WP, Luthe DS: Proteome profile of the developing maize (Zea mays L.) rachis. Proteomics. 2010, 10: 3051-3055. 10.1002/pmic.200900833.

    PubMedCASArticleGoogle Scholar

  45. Liu Q, Zhu A, Chai L, Zhou W, Yu K, Ding J, Xu J, Deng X: Transcriptome analysis of a spontaneous mutant in sweet orange [Citrus sinensis (L.) Osbeck] during fruit development. J Exp Bot. 2009, 60: 801-10.1093/jxb/ern329.

    PubMedCASPubMed CentralArticleGoogle Scholar

  46. Xue L-J, Zhang J-J, Xue H-W: Genome-wide analysis of the complex transcriptional networks of rice developing seeds. PLoS One. 2012, 7: e31081-10.1371/journal.pone.0031081.

    PubMedCASPubMed CentralArticleGoogle Scholar

  47. Yang J, An D, Zhang P: Expression profiling of cassava storage roots reveals an active process of glycolysis/GluconeogenesisFree access. J Integr Plant Biol. 2011, 53: 193-211. 10.1111/j.1744-7909.2010.01018.x.

    PubMedCASArticleGoogle Scholar

  48. Zheng J, Fu J, Gou M, Huai J, Liu Y, Jian M, Huang Q, Guo X, Dong Z, Wang H: Genome-wide transcriptome analysis of two maize inbred lines under drought stress. Plant mol biol. 2010, 72: 407-421. 10.1007/s11103-009-9579-6.

    PubMedCASArticleGoogle Scholar

  49. Ya HY, Chen QF, Wang WD, Qin GY, Jiao Z: Pathway clustering of the Up-regulated genes in the stimulated rice seedlings that germinating from the seeds implanted by the Low-energy nitrogen Ion. Procedia Environmental Sciences. 2011, 8: 679-686.

    CASArticleGoogle Scholar

  50. Teale WD, Paponov IA, Palme K: Auxin in action: signalling, transport and the control of plant growth and development. Nat Rev Mol Cell Biol. 2006, 7: 847-859. 10.1038/nrm2020.

    PubMedCASArticleGoogle Scholar

  51. Sorefan K, Girin T, Liljegren SJ, Ljung K, Robles P, Galván-Ampudia CS, Offringa R, Friml J, Yanofsky MF, Østergaard L: A regulated auxin minimum is required for seed dispersal in Arabidopsis. Nature. 2009, 459: 583-586. 10.1038/nature07875.

    PubMedCASArticleGoogle Scholar

  52. Moyle R, Schrader J, Stenberg A, Olsson O, Saxena S, Sandberg G, Bhalerao RP: Environmental and auxin regulation of wood formation involves members of the Aux/IAA gene family in hybrid aspen. Plant J. 2002, 31: 675-685. 10.1046/j.1365-313X.2002.01386.x.

    PubMedCASArticleGoogle Scholar

  53. Schrader J, Baba K, May S, Palme K, Bennett M, Bhalerao R, Sandberg G: Polar auxin transport in the wood-forming tissues of hybrid aspen is under simultaneous control of developmental and environmental signals. Proc Natl Acad Sci USA. 2003, 100: 10096-10.1073/pnas.1633693100.

    PubMedCASPubMed CentralArticleGoogle Scholar

  54. Nilsson J, Karlberg A, Antti H, Lopez-Vernaza M, Mellerowicz E, Perrot-Rechenmann C, Sandberg G, Bhalerao RP: Dissecting the molecular basis of the regulation of wood formation by auxin in hybrid aspen. Plant Cell. 2008, 20: 843-855. 10.1105/tpc.107.055798.

    PubMedCASPubMed CentralArticleGoogle Scholar

  55. Zhou J, Sebastian J, Lee JY: Signaling and gene regulatory programs in plant vascular stem cells.Genesis. 2011, 49: 885-904.

    Google Scholar

  56. Love J, Björklund S, Vahala J, Hertzberg M, Kangasjärvi J, Sundberg B: Ethylene is an endogenous stimulator of cell division in the cambial meristem of populus. Proc Natl Acad Sci USA. 2009, 106: 5984-5989. 10.1073/pnas.0811660106.

    PubMedCASPubMed CentralArticleGoogle Scholar

  57. Nieminen K, Immanen J, Laxell M, Kauppinen L, Tarkowski P, Dolezal K, Tähtiharju S, Elo A, Decourteix M, Ljung K: Cytokinin signaling regulates cambial development in poplar. Proc Natl Acad Sci USA. 2008, 105: 20032-20037. 10.1073/pnas.0805617106.

    PubMedCASPubMed CentralArticleGoogle Scholar

  58. Bi T, Li L: Characterization of the plasma membrane proteins and receptor-like kinases associated with secondary vascular differentiation in poplar. Plant mol biol. 2011, 76: 97-115. 10.1007/s11103-011-9771-3.

    PubMedCASPubMed CentralArticleGoogle Scholar

  59. Tsuchihira A, Hanba YT, Kato N, Doi T, Kawazu T, Maeshima M: Effect of overexpression of radish plasma membrane aquaporins on water-use efficiency, photosynthesis and growth of eucalyptus trees. Tree physiol. 2010, 30: 417-430. 10.1093/treephys/tpp127.

    PubMedCASArticleGoogle Scholar

  60. Dharmawardhana P, Brunner A, Strauss S: Genome-wide transcriptome analysis of the transition from primary to secondary stem development in populus trichocarpa. BMC genomics. 2010, 11: 150-10.1186/1471-2164-11-150.

    PubMedPubMed CentralArticleGoogle Scholar

  61. Li YZ, Pan YH, Sun CB, Dong HT, Luo XL, Wang ZQ, Tang JL, Chen B: An ordered EST catalogue and gene expression profiles of cassava (Manihot esculenta) at key growth stages. Plant mol biol. 2010, 74: 573-590. 10.1007/s11103-010-9698-0.

    PubMedCASArticleGoogle Scholar

  62. Wang Y, Du ST, Li LL, Huang LD, Fang P, Lin XY, Zhang YS, Wang HL: Effect of CO2 elevation on root growth and its relationship with indole acetic acid and ethylene in tomato seedlings. Pedosphere. 2009, 19: 570-576. 10.1016/S1002-0160(09)60151-X.

    CASArticleGoogle Scholar

  63. Li R, Zhu H, Ruan J, Qian W, Fang X, Shi Z, Li Y, Li S, Shan G, Kristiansen K: De novo assembly of human genomes with massively parallel short read sequencing. Genome res. 2009, 20: 265-272.

    PubMedArticleGoogle Scholar

  64. Iseli C, Jongeneel CV, Bucher P: ESTScan: a program for detecting, evaluating, and reconstructing potential coding regions in EST sequences. Proc Int Conf Intell Syst Mol Biol. 1999, 138-148.

    Google Scholar

  65. Audic S, Claverie JM: The significance of digital gene expression profiles. Genome Res. 1997, 7: 986-995.

    PubMedCASGoogle Scholar

  66. Benjamini Y, Yekutieli D: The control of the false discovery rate in multiple testing under dependency. The annals of statistics. 2001, 29: 1165-1188. 10.1214/aos/1013699998.

    ArticleGoogle Scholar

  67. Hafner M, Landgraf P, Ludwig J, Rice A, Ojo T, Lin C, Holoch D, Lim C, Tuschl T: Identification of microRNAs and other small regulatory RNAs using cDNA library sequencing. Methods. 2008, 44: 3-12. 10.1016/j.ymeth.2007.09.009.

    PubMedCASPubMed CentralArticleGoogle Scholar

Download references

Acknowledgements

We acknowledge the Beijing Genomics Institute at Shenzhen for its assistance in original data processing and related bioinformatics analysis. We would like to thank Zhong-fu Deng and Jian-jiu Liu for their help in sample collection. This work was supported by grants from the National Natural Science Foundation of China (30972332 and 31100454).

Author information

Affiliations

  1. State key laboratory of tree genetics and breeding, Research Institute of Forestry, Chinese Academy of Forestry, Beijing, 100091, China

    Cai-yun He, Kai Cui, Jian-guo Zhang, Ai-guo Duan & Yan-fei Zeng

  2. Key Laboratory of Tree Breeding and Cultivation, State Forestry Administration, Research Institute of Forestry, Chinese Academy of Forestry, Beijing, 100091, China

    Cai-yun He, Kai Cui, Jian-guo Zhang, Ai-guo Duan & Yan-fei Zeng

  3. Research Institute of Resources Insects, Chinese Academy of Forestry, Kunming, 650224, China

    Kai Cui

Authors
  1. Cai-yun He
    View author publications

    You can also search for this author inPubMedGoogle Scholar

  2. Kai Cui
    View author publications

    You can also search for this author inPubMedGoogle Scholar

  3. Jian-guo Zhang
    View author publications

    You can also search for this author inPubMedGoogle Scholar

  4. Ai-guo Duan
    View author publications

    You can also search for this author inPubMedGoogle Scholar

  5. Yan-fei Zeng
    View author publications

    You can also search for this author inPubMedGoogle Scholar

Corresponding author

Correspondence toJian-guo Zhang.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

JGZ designed and supervised the study. CYH designed experiment and analyzed the data as a whole and wrote the manuscript. KC performed proteome analysis. AGD and YFZ coordinated character investigation and collected bamboo samples. All authors read and approved the final manuscript

Electronic supplementary material

12870_2013_1333_MOESM1_ESM.jpeg

Additional file 1: Figure S1: Quality and coverage evaluation of assembled unigenes. Distribution of unique-mapped reads of the assembled unigenes. A: Length distribution of unigene. B: Read distribution of unigene. (JPEG 684 KB)

12870_2013_1333_MOESM2_ESM.jpeg

Additional file 2: Figure S2: Gene Ontology (right) and clusters of orthologous groups (COG) (left) classification of assembled unigenes. A total of 32,064 unigenes with BLAST matches to known proteins were assigned to three main categories: biological process, cellular component and molecular function. Out of 42,127 Nr hits, 13,957 sequences were assigned to 25 COG classifications. (JPEG 1 MB)

Additional file 3: Table S1: List of 125 KEGG pathways for unigenes. (XLS 34 KB)

12870_2013_1333_MOESM4_ESM.jpeg

Additional file 4: Figure S3: Distribution of tags and gene expression among G1M-G4M, G3T and G3B groups. The distribution of tags matches that of gene expression among groups. Furthermore, an increase in tags or gene expression is accompanied by a decrease in the frequencies of tags or genes expression. (JPEG 1 MB)

Additional file 5: Table S2: Summary of tag mapping in DGE analysis for G1M–G4M, G3T and G3B groups. (XLS 102 KB)

Additional file 6: Figure S4: Sequencing saturation evaluation of G1M-G4M, G3T and G3B libraries. (JPEG 660 KB)

12870_2013_1333_MOESM7_ESM.jpeg

Additional file 7: Figure S5: Hierarchical clustering of differentially expressed genes, miRNAs and proteins. A total of 213 differentially expressed proteins (DEPs) were detected by two-dimensional gel electrophoresis (2-DE) and identified by matrix-assisted laser desorption/ionization time-of-flight/time-of-flight mass spectrometry (MALDI-TOF/TOF MS). G1M, G2M, G3M and G4M represent four developmental stages in turn. G3T, G3M and G3B represent top, middle and basal internode of the third developmental stage, respectively. The color ranges from green to red for the down-regulated and up-regulated genes, respectively. (JPEG 3 MB)

12870_2013_1333_MOESM8_ESM.xls

Additional file 8: Table S3: A. The significant STEM clusters of mRNAs. Eleven (1,024 genes) development-specific and six (923 genes) internode-specific gene profiles (P-Value ≤ 0.01) were obtained separately. B: The significant STEM profiles of miRNAs. Four (54 miRNAs) development-specific and one (101 miRNAs) internode-specific miRNA profiles were identified separately. (XLS 408 KB)

Additional file 9: Table S4: Statistics of small RNA sequences from the individual libraries. (XLSX 13 KB)

Additional file 10: Table S5: The list of know-miRNA and novel-miRNA in each libraries. (XLSX 114 KB)

Additional file 11: Table S6: The list of miRNAs and their predicted target genes. (XLS 2 MB)

12870_2013_1333_MOESM12_ESM.xls

Additional file 12: Table S7: A: The negative correlation pairs between miRNAs and mRNAs, the positive correlation pairs between mRNAs and proteins. B: The statistical data pairs from the integrated analysis of STEM, positive/negative correlation. C: The GO annotation and enrichment of 64 core genes. (XLS 228 KB)

Authors’ original submitted files for images

Rights and permissions

Open AccessThis article is published under license to BioMed Central Ltd. This is an Open Access article is distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Reprints and Permissions

About this article

Cite this article

He, Cy., Cui, K., Zhang, Jg.et al.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 Biol13,119 (2013). https://doi.org/10.1186/1471-2229-13-119

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:https://doi.org/10.1186/1471-2229-13-119

Keywords

  • Digital Gene Expression
  • Zeatin Riboside
  • Bamboo Culm
  • Moso Bamboo
  • Digital Gene Expression Library

BMC Plant Biology

ISSN: 1471-2229

Contact us