Discovery and Comparative Profiling of microRNAs in Representative Monopodial Bamboo (Phyllostachys edulis) and Sympodial Bamboo (Dendrocalamus latiflorus)

Background According to the growth pattern of bamboo, sympodial bamboo and monopodial bamboo are considered as two mainly kinds of bamboo. They have different phenotypes and different characteristics in developmental stage. Much attention had been paid on the study of bamboo cultivation, processing, physiology, biochemistry and molecular biology, which had made great progresses in the last decade, especially for the highlighted achievement of the bamboo genomics. However, there is no information available on concerning comparative profiling of miRNAs between sympodial bamboo and monopodial bamboo, which might play important roles in the regulation of bamboo development. Methodology/Principal Findings We identified the profiles of small RNAs using leaf tissues from one sympodial bamboo i.e. moso bamboo (Phyllostachys edulis) and another monopodial bamboo i.e. ma bamboo (Dendrocalamus latiflorus). The result showed that there were 19,295,759 and 11,513,888 raw sequence reads, in which 92 and 69 conserved miRNAs, as well as 95 and 62 novel miRNAs were identified in moso bamboo and ma bamboo, respectively. The ratio of high conserved miRNA families in ma bamboo is more than that in moso bamboo. In addition, a total of 49 and 106 potential targets were predicted in moso bamboo and ma bamboo, respectively, in which several targets for novel miRNAs are transcription factors that play important roles in plant development. More importantly, annotation of differentially expressed target genes was performed based on the analysis of pathway and gene ontology terms enrichment. Conclusions/Significance This study provides the first large-scale sight of discovery and comparative characterization of miRNAomes between two representative bamboos belonged to sympodial bamboo and monopodial bamboo, respectively. Although it will be necessary to validate the function of miRNAs through more experimental research in further, these results lay a foundation for unraveling the miRNA-mediated molecular processes in different kinds of bamboo.


Introduction/Background
MicroRNAs (miRNAs) are endogenous, non-coding small RNAs of approximately 22 nucleotides that regulate the flow of genetic information by controlling the translation or stability of mRNAs. It has been estimated that miRNAs account for ,1% of predicted genes in higher eukaryotic genomes, and that up to 10-30% of genes may be regulated by miRNAs, which exist widely in animals, plants and some viruses. Since the discovery of the first miRNA in Caenorhabditis elegans, a large number of miRNAs have been identified by experimental cloning or bioinformatics methods with the development of biotechnology. More importantly, with the increasing development of high-throughput technology, the integrated analysis based on miRNA sequences in the evolutionary process provides a significant number of evidences that miRNAs evolved independently in different strains and generally were conserved in the process of evolution [1,2]. However, some data and studies associated function with evolution for bamboo miRNAs have been rare till now.
As a tribe of flowering and evergreen perennial monocots, classified in the subfamily Bambusoideae within the grass family Poaceae that includes rice, maize, wheat and other cereals [3], bamboo represents one of the fastest-growing plants and one of the most important non-timber forest resources in the world [4,5]. More importantly, bamboo is of notable economic and environmental significance, due to being used for versatile raw products. About 2.5 billion people depend on bamboo economically, and the international trade in bamboo amounts to over 2.5 billion US dollars per year [6].
Bamboo has rhizomes and roots underground that lives exclusively in forest and grows into large woody culms. There are two mainly patterns of bamboo for their growth, namely as sympodial bamboo and monopodial bamboo. Monopodial bamboo grows fast, through their roots or rhizome, which can spread widely underground and send up new culms to break through the surface. Therefore, this kind was named as running bamboo vividly. Monopodial bamboo ranges in all shapes and sizes, from a modest 30 cm (i.e. pygmy bamboo) to a mighty 18-20 m in height (i.e. moso bamboo). As one of significant bamboo in monopodial bamboo, moso bamboo (Phyllostachys edulis) is the most important economic one with many advantages such as fast growing rate, high yield, extensive use, short crucial period formation and strong regeneration capacity.
Sympodial bamboo grows like a tussock, and their new clums emerged around the outside edge of clump, which was aptly named as clumping bamboo. Due to unique root and rhizome structures, sympodial bamboo produces a tight cluster of culms and spreads only a couple centimeters out from the base each season. They range in height from less than 20 cm to more than 30 m, depending on different species of bamboo. Distinguished from monopodial bamboo, sympodial bamboo have some special features, such as lower spreading rate and u-shaped rhizome making new culms next to original plant. As one of the most popular and valuable bamboo species in sympodial bamboo, ma bamboo (Dendrocalamus latiflorus) fascinated people because of its evergreen color and delicious shoot, which was widely cultivated in southern China.
To facilitate omics studies in bamboo and help unveil the miRNAs feature of bamboo, small-RNA libraries of leaf tissues from moso bamboo and ma bamboo, respectively, were constructed and sequenced. In an attempt to identify new conserved miRNAs and novel miRNAs in bamboo, a total of 7,450 known miRNAs from 70 plants species from miRBase served as a reference data set. Moreover, the potential target genes of novel miRNAs were predicted based on the high-throughput technology and bioinformatics analysis. The profiling of micro-RNAs from moso bamboo and ma bamboo were compared. Ultimately, these results will play the necessary role in understanding the regulation of miRNAs in the development of monopodial bamboo and sympodial bamboo.

Plant material and RNA isolation
Divided seedlings of moso bamboo and ma bamboo were potted in our laboratory under a regime of 16 h light and 8 h darkness at 25uC, with a light intensity of 200 mmol?m 22 ?s 21 and a relative humidity of 75%. The leaf tissues were used as materials for this study. We chose the third piece of new functional leaf (blade tissue only) from the top of the branch, which could be considered as a juvenile leaf. The leaves were collected and quickly frozen in liquid nitrogen. Total RNA was isolated from leaf tissues using the Trizol reagent (Invitrogen, Carlsbad, CA, and USA), according to the manufacturer's instructions strictly.

Small RNA high-throughput sequencing
The small RNA library construction and Solexa sequencing were carried out at BGI-Shenzhen (Shenzhen, China) using the standard small RNA sample preparation protocol (Illumia) as described by Hafner et al. [7]. Briefly, small RNAs of 15-30 nt in length were first isolated from the total RNA through 15% TBE urea denaturing polyacrylamide gels, which were ligated to the small RNA-59 adaptor, and then 39 adaptor was ligated to the small RNA-59 adaptor, followed by reverse transcription into cDNAs. These cDNAs were amplified by PCR and subjected to Solexa sequencing.

Sequence data analysis
After removing low quality reads and trimming adapter sequences, initial reads based on Solexa sequencing were processed by summarizing data production, evaluating sequencing quality, calculating the length distribution of small RNA reads. Moreover, small RNAs ranging from 18-30 nt were collected and used for further analyses. The raw reads are available in the NCBI SRA database under the accession number SRX480448. For analyzing conserved and novel miRNA, firstly, the clean reads, which mapped to exons or other non-coding RNAs such as rRNA, snRNA, snoRNA and tRNA, are filtered by BLAST against the coding sequences of moso bamboo genome, the Rfam database (http://rfam.sanger.ac.uk/) [8] and the GenBank non-coding RNA database (http://www.ncbi.nlm.nih.gov/) to discard coding gene, non-coding RNA sequences. Subsequently, unique sRNAs were aligned with plant mature miRNAs in miRBase Release 20. After rigorous screening, sRNA sequences without more than three mismatched bases were selected by BLAST searching against miRBase. Then, the remaining reads were used to map the genome of moso bamboo. Sequences with a tolerance of two mismatches were retained for miRNA prediction. RNAfold (http://www.tbi.univie.ac.at/RNA/) [9] was used for secondary structure prediction (hairpin prediction) of individual mapped miRNAs, using the default folding conditions to identify conserved miRNAs in bamboo. Lastly, sequences without defined as the conserved miRNAs were termed as novel miRNAs.
RPKM (Reads Per Kilobase per Million) is a method of quantifying gene expression from RNA sequencing data by normalizing for total read length and the number of sequencing reads [10]. Based on the same principle, RPKM is suitable the analysis of miRNA as well. The equation for calculating RPKM is in the following. In addition, to reasonably compare the miRNA abundance in libraries of moso bamboo and ma bamboo, the value of RPKM were calculated to seek the normalization of the difference between two samples. As a unique exact match, the information that a read mapped only once in the genome with no mismatch was also provided.
Potential target sequences for the novel miRNAs were predicted on the basis of the psRNATarget program (http://plantgrn.noble. org/psRNATarget/) with default parameters [11]. As the requirement of program, the newly identified miRNA sequences were used as custom miRNA sequences, while coding sequences in moso bamboo genome were used as custom plant databases. All predicted target genes were evaluated by the scoring system and criteria defined in a previous report [12].
To obtain custom enrichments, Gene Ontology (GO) terms enrichment was calculated by Ontologizer [13], which is a java application that can be used to perform statistical analysis for enrichment of GO terms in set of genes or proteins based on the one-sided Fisher's exact test, the novel parent-child method, and topology-based algorithms. As the inputs for Ontologizer, the four files were needed, including one GO file with the format of open biomedical ontologies, one list of all gene IDs of moso bamboo, one list of gene IDs that will be analyzed and one association file of moso bamboo. As parameter settings, term for term analysis with model-based gene set was used. GO terms analysis was performed separately on the set of target genes in sample of moso bamboo and ma bamboo. Besides, GO categories with p-values less than 0.05 were considered significant. Lastly, the biological interpretation of the target genes of differential miRNAs was completed using KEGG (Kyoto Encyclopedia of Genes and Genomes) pathway analysis.

Expression analysis of miRNAs by qRT-PCR
To detect the expression level of miRNA in bamboo, stem-loop qRT-PCR was employed in this study. Primers were designed specifically based on individual miRNA using previous method [14]. U6 snRNA was selected as internal control [15]. More detailed information of primers was supported in File S1.
On the basis of leaf samples in moso bamboo and ma bamboo, cDNAs were synthesized from total RNA with the miRNA-specific stem-loop RT primer according to reference [16]. Subsequently, qRT-PCR was carried out using a SYBR Green I Master Kit (Roche, Germany) on a QTOWER2.2 Real-Time PCR System (Analytik Jena). The final volume was 10 ml, containing 5.0 ml 26SYBR Premix Ex Taq, 0.2 ml of each primer (10 mM), 0.8 ml of cDNA and 3.8 ml of nuclease-free water. All reactions were repeated three times and the amplification was conducted as follows: initial denaturation at 95uC for 10 min, followed by 50 cycles at 95uC for 10 s, and 62uC for 10 s. There was three biological experiments. Expression levels were normalized to that of the internal control, and the relative value was calculated using previous method [17].

MiRNA ioslation by PCR and sequencing
To prove the existence of miRNAs in bamboo, stem-loop primers were also directly used to amplify the mature miRNA sequences using specific templates from moso bamboo and ma bamboo respectively. Moreover, the primers designed for precursors (File S2) were used to amplify with DNA templates of moso bamboo and ma bamboo respectively. All the amplified products were sequenced and analyzed.

Summary of small RNA library sequencing
Based on Solexa sequencing of small RNA libraries from moso bamboo and ma bamboo (File S3), 19,295,759 and 11,513,888 raw reads were produced, respectively. After discarding low quality sequences and collapsed raw data, there were 8,089,795 (moso bamboo) and 10,593,305 (ma bamboo) clean reads as well as 2,653,717 (moso bamboo) and 6,320,379 (ma bamboo) unique sequences. Subsequently, the sequences, which mapped to coding sequences of moso bamboo genome and non-coding RNAs such as tRNAs, rRNAs, siRNAs, snRNA, and snoRNA, were removed. The statistics of small RNA library sequencing were shown in Table 1. As an important feature of the size profile, the distribution of small RNA length was summarized in Figure 1. The general distribution profile of moso bamboo was similar to that of ma bamboo. The most of small RNAs were ranged from 20 nt to 24 nt in length, which were highly consistent with those of small RNAs for known function [18]. As a maximum peak, small RNA with 24 nt in length was also highly agreed with those small RNAs in other Poaceae plants based on Solexa sequencing technology [19][20][21].

Identification of conserved miRNAs
We identified 92 and 69 conserved miRNAs, belonged to 53 families and 35 families for moso bamboo and ma bamboo, respectively, based on sequence alignment by BLAST against the latest miRBase v20, in which miRNAs of bamboo were no available, and a series of strict filtering criteria (File S4).
Normalized analysis of miRNA sequence between moso bamboo and ma bamboo were employed. The number of mapped reads, the length of mapped reads and the number of total mapped reads were considered to make it possible to compare RNA expression levels from different samples. In addition, the previously reports shown miRNAs with high sequencing frequencies played some key regulatory functions in maintaining biological process [22]. Consequently, the read counts for known miRNA families were evaluated in two samples. The result indicated small RNA families of moso bamboo were mapped to numerous conserved miRNA families, while those of ma bamboo were only concentrated in some conserved miRNA families. The reads with RPKM.5000 in two samples were focused on the abundance of conserved miRNA families, which accounted for a high proportion (,99%) of the total conserved miRNA both in ma bamboo and moso bamboo. Apart from high abundant miRNA families, one must bear in mind that, there are a large number of low abundant miRNA families, indicating that these miRNAs were expressed at a low level. Another feature is that some conserved miRNAs were identified only in one sample. For example, miR396, miR397, miR1432 and miR7748 were identified only in moso bamboo, wherever miR170 was only in ma bamboo. This may be explained by the fact that the two kinds of bamboo had their own diversely biological and developmental features in the long progress of evolution. However, the previously study on the miRNA of ma bamboo demonstrated that miR396 family was identified [12]. It was may be caused by tissue and process specific miRNA in two different samples of ma bamboo [23,24]. Taken together, as specific conserved miRNAs in moso bamboo, miR397, miR1432 and miR7748 were focused on using experimental validations further.
The number of members in different conserved miRNA families was also analyzed. The number of conserved miRNA families with more than 3 members was seven both in moso bamboo and ma bamboo. In addition, those with 3 members or less in both of moso bamboo and ma bamboo were accounted more than half of the total miRNA families.

Evolution of miRNA in bamboo
Due to the important role of miRNA in regulation of plant gene expression, researchers have focused on the prediction, identification and functional analysis of miRNA in bamboo. However, the evolution of miRNA in bamboo has not been reported. To investigate the evolutionary roles of the conserved miRNAs, comparisons against conserved miRNAs in plants described by Cuperus et al [25] were carried out. The outcome showed almost of highly conserved miRNA families of bamboo were found (File S5), which further proved the previous arguments that miRNAs were highly conserved in sequences [26]. In addition, the previous study indicated that miRNA families from miR156 to miR408 were identified as high conserved miRNAs in the plant kingdom [27,28]. These miRNAs play essential and conserved functions in plant development, such as flower and leaf development. MiRNAs families behind miR408 are low conserved or non-conserved, which may play roles in more species-specific characteristics in plant growth and development. These miRNA families were considered as the new generation in evolutional process. As shown in Figure 2, comparing conservation for miRNA families between two samples, a total of RPKM ranged from miR156 to miR408 accounted for 69.1%, 81.0% in moso bamboo and ma bamboo, respectively. Moreover, there is higher ratio of low conserved miRNAs families in moso bamboo (31%) than that of ma bamboo (19%). It was indicated that the implementation in regulatory processes, which played more essential roles, might need more involvement of miRNA in ma bamboo leaf than those in moso bamboo.
The previously study demonstrated that miRNA sequences have to complement with their target gene transcripts to carry out their functions, which suggested that miRNA genes have to coevolve with their target genes [29]. MiRNA genes would amplify through duplication events similar to those that drive the evolution of protein-coding genes: tandem gene duplications, segmental duplications, and chromosomal duplications or polyploidization [30,31]. In this study, the different of proportion in low conserved miRNAs figures out more miRNAs with low conserved will be need to regulate newly emerging genes, which help moso bamboo Table 1. Summary of small RNAs sequences from moso bamboo and ma bamboo. adapt to the changing environments better than ma bamboo does [32,33]. The rhizome of monopodial bamboo can spread laterally grown in soil, and also be differentiated into bamboo away from the mother plant, while sympodial bamboo grows in clusters within a relatively small range. Therefore, monopodial bamboo requires more self-regulation to adapt to the environment than sympodial bamboo, which might be consistent with the generation of lower conserved miRNAs families. Studies have shown that the expression of miRNAs is closely related with the developmental stages of plant. Based on the analysis of the previous study on the miRNAs of moso bamboo [34], the analysis based on percentages of conserved miRNA families in different development of internode demonstrated that highly conserved miRNAs may play more essential roles in earlier development of moso bamboo due to some stress response, while non-conserved miRNAs with low expression occupied a dominant position in maintaining the regulatory function during mature stage of moso bamboo ( Figure 3). Moreover, the results based on percentages of conserved miRNA families in different positions of same development stage illustrated that highly conserved miRNAs from different positions of same tissue were possibly stable in one stage.

Category
To further understand the feature of bamboo, the related studies on other species of bamboo will be properly carried out.

Identification of novel miRNAs
According to the analysis, we predicted 95 and 62 candidates of novel miRNAs with more than 3 reads in moso bamboo and ma bamboo, respectively. The mature sequences of novel miRNA with RPKM.300 in two samples were shown in Table 2. The comprehensive information on novel mature miRNAs was provided in File S6.
The number of novel miRNAs and unique mapped reads of ma bamboo were less than those of moso bamboo, probably because    the genome of moso bamboo was analyzed as required genome. Although both of moso bamboo and ma bamboo belonged to Bambusoideae, phylogenetic analyses between these had few study for the lack of the genome information of ma bamboo. Therefore, using the genome of moso bamboo for analyzing miRNA of ma bamboo would inevitably produce some biases. To address the gap, the whole genome sequencing of ma bamboo need be implemented as soon as possible.
Among the novel miRNAs, phe-miRC1-1 with 106143.04 RPKM, more than 3-folds of the first-highest in ma bamboo, had the highest expression level in moso bamboo. On the basis of their frequencies and sequences in the small RNA libraries, although the expression levels of these candidates ranged from thousands of RPKM to lower, in general, novel miRNA candidates showed lower expression level compared with most of the conserved families. The low abundance of novel miRNAs demonstrated that these miRNAs might play a specific and essential role in certain tissues or developmental stages, and were considered as young miRNAs in terms of evolution [25].

Prediction of miRNA targets
The acquisition of knowledge for miRNA target genes contributed to insight into the range of miRNA regulation and unveiling detailed descriptions of miRNA-target interactions. In addition, as an essential step to identify miRNA targets, the methods based on bioinformatics were widely used [11]. Herein, on the basis of the criteria described in the material and methods, a total of 49 and 106 target genes were identified in moso bamboo and ma bamboo, respectively. The target genes of novel miRNAs with RPKM.300 in two samples were shown in Table 3. The comprehensive information on target genes was provided in File S7.
The target genes of the novel miRNAs from bamboo mainly concentrated on transcription factor, signal regulation and functional protein, among which were a considerable number of SPL genes, indicating that these miRNAs might be involved in regulating the expression of SPL genes. SPL genes family is a group of structurally diverse genes encoding putative transcription factors found in photosynthetic organisms [35,36]. The distinguishing feature of the SPL gene family is the SBP-box encoding a conserved protein domain of 76 amino acids in length. Recent studies indicated that SBP-box genes played central roles in controlling certain homeostatic processes, as well as in overcoming stresses [37,38]. In the previous study of Arabidopsis, several targeted SPL genes act redundantly in controlling the juvenile-toadult phase transition, which renders plants sensitive to photoperiodic induction of flowering [39,40]. Therefore, the miRNAs targeted SPL genes might play important roles in bamboo leaves, which response to environmental factors especially for light induction.

Annotation of differentially expressed target genes
Functional annotations were carried out to investigate which processes are differentially regulated by miRNA in moso bamboo and ma bamboo. To infer significant associations between the sets of targets predicted by individual miRNAs and specific biochemical pathways, a total of 20 different metabolic pathways were found, including 5 metabolism, 8 cellular processes, 7 genetic information processing (Table 4), which indicated that miRNA played essential roles in various biological process.
Besides pathway analysis, another functional annotation was performed via the analysis of GO, which is the de facto standard in gene functionality description and is used widely in functional annotation and enrichment analysis [41]. Genes from a large amount of organisms have been annotated to GO terms. The widespread applications based on GO terms enrichment analysis were the identification of annotation-enriched GO terms in a list of genes with some similar characteristics in biology. These terms were often considered as some representations with the outstanding biological features of the genes in the study of set [13]. We performed GO analysis on the result of prediction of miRNA targets using model-based gene set analysis (MGSA), which analyzes all GO terms at once by embedding them in a Bayesian network in order to provide high-level, summarized views of core biological processes.
In order to identify significant target genes, the GO term for target genes with p-value ,0.05 were selected and shown in Figure 4. GO term enrichment analysis identified a total of 20 and 32 GO terms with high significant in moso bamboo and ma bamboo, respectively. The most highly ranked term of the target genes with higher significant in moso bamboo was DNA binding (GO: 0003677), which included 17 GO terms from the total 100 GO terms in the set of the moso bamboo. On the other hand, the most highly ranked term for the genes with significant in ma bamboo was cellular component (GO: 0005575) including 19 ones of the total 69 GO terms in the set of ma bamboo. Moreover, there were 4 GO terms overlapped (GO: 0003677, GO: 0005575, GO: 0043231 and GO: 0043227) in the total of significantly GO terms, of which 3 GO terms were distributed in cellular component part (GO: 0005575, GO: 0043231 and GO: 0043227). Furthermore, more highly ranked terms for the targets genes with higher expression were both in biological process, based on the significant analysis of GO terms between moso bamboo and ma bamboo. For example, on account of the category of GO terms, the ratio of biological process, cellular component and molecular function in moso bamboo was 40%, 30% and 30% of the total of GO terms, respectively, while the corresponded ratio were change to 65.6%, 12.5% and 21.9% in ma bamboo. Higher proportion of biological process in ma bamboo possibly reflected the significance of target genes, which might play more essential role than those in moso bamboo. Validation of specific conserved miRNAs by qRT-PCR and sequencing Specific primers designed for mature miRNA were used to validate the expression of miR396, miR397, miR1432, and miR7748 in leaves of both moso bamboo and ma bamboo. The result of qRT-PCR demonstrated that miR396, miR397, miR1432, and miR7748 were expressed differently in leaf of moso bamboo, among which miR396 was the highest one, followed by miR1432 and miR397, and miR7748 was the lowest one ( Figure 5).
While that in ma bamboo, miR397 was the highest one, followed by miR1432 and miR396, miR7748 was not detectable (File S8). To prove the existence of miRNAs in ma bamboo, stemloop primers were also directly used to amplify the mature miRNA sequences using specific templates from moso bamboo and ma bamboo, respectively, and the amplified products were sequenced. The result showed the sequences of mature miRNA in moso bamboo were all consistent with the sequences of miR396, miR397, miR1432, and miR7748 correspondingly. However, only three pairs of primers had specific amplification products for ma bamboo, one is the sequence of miR396,consistented with that of moso bamboo, while the other sequences were mismatched with miR397 and miR1432, respectively (File S9), which indicated the qRT-PCR result of miR397 and miR1432 for ma bamboo is false positive.
Moreover, the primers designed for precursors of miR396, miR397, miR1432, and miR7748 were used to amplify with DNA templates of moso bamboo and ma bamboo, respectively. The sequencing result demonstrated the sequences of moso bamboo were matched with the genomic sequences and could form the stem-loop structures (File S10). However, there were no specific PCR products from ma bamboo. In order to identify the targets of these miRNAs, we had screened the published ma bamboo data of transcriptome [42] and ESTs [43] in NCBI and bambooGDB database [44] using the targets of miR396, miR397, miR1432, and miR7748, respectively, and no matched sequences was found for them, which indicated that the target genes in moso bamboo and ma bamboo are not comparable. On the basis of these results, miR397, miR1432, and miR7748 were specific conserved in the leaf sample of moso bamboo.

Conclusions
To the best of our knowledge, this work presents the firstly comparative profiling of microRNAs in representative monopodial bamboo (moso bamboo) and sympodial bamboo (ma bamboo), based on high-throughput sequencing analysis of small RNA. A large number of conserved miRNAs, novel miRNAs and target genes were identified, which are worthy of further investigation to improve our understanding of the regulatory mechanisms of miRNA in biological processes in bamboo. Target prediction of differential miRNAs indicated that SPL genes are mainly regulated by miRNAs. Functional annotation of genes targeted by differential miRNAs revealed that high significant biological processes are transcription regulation and photosynthesis. The expression of novel miRNAs and target genes involved in developmental stage of bamboo were further validated. Experiments showed that miR397, miR1432, and miR7748 were specific conserved in the leaf sample of moso bamboo. Taken together, the comparison between moso bamboo and ma bamboo indicated that monopodial bamboo and sympodial bamboo may be share some different miRNAs and theirs target genes in order to better adapt the development in different stages and stress response in their diverse course of evolution. Therefore, this study should lay an important foundation for future genetic or genomic studies on bamboo and help to fill the gap in comparative profiling of microRNAs.

Supporting Information
File S1 Primers for qRT-PCR and PCR of miRNAs.