Identification and Characterization of MicroRNAs in the Leaf of Ma Bamboo (Dendrocalamus latiflorus) by Deep Sequencing

MicroRNAs (miRNAs), a class of non-coding small endogenous RNAs of approximately 22 nucleotides, regulate gene expression at the post-transcriptional levels by targeting mRNAs for degradation or by inhibiting protein translation. Thousands of miRNAs have been identified in many species. However, there is no information available concerning miRNAs in ma bamboo (Dendrocalamus latiflorus), one of the most important non-timber forest products, which has essential ecological roles in forests. To identify miRNAs in D. latiflorus, a small RNA library was constructed from leaf tissues. Using next generation high-throughput sequencing technology and bioinformatics analysis, we obtained 11,513,607 raw sequence reads and identified 84 conserved miRNAs (54 mature miRNAs and 30 star miRNAs) belonging to 17 families, and 81 novel miRNAs (76 mature miRNAs and five star miRNAs) in D. latiflorus. One hundred and sixty-two potential targets were identified for the 81 novel bamboo miRNAs. Several targets for the novel miRNAs are transcription factors that play important roles in plant development. Among the novel miRNAs, 30 were selected and their expression profiles in response to different light conditions were validated by qRT-PCR. This study provides the first large-scale cloning and characterization of miRNAs in D. latiflorus. Eighty-four conserved and 81 novel miRNAs were identified in D. latiflorus. Our results present a broad survey of bamboo miRNAs based on experimental and bioinformatics analysis. Although it will be necessary to validate the functions of miRNAs by further experimental research, these results represent a starting point for future research on D. latiflorus and related species.

Bamboo, one of the most important non-timber forest products among the world's plant and forest resources, is widely distributed in the tropical and subtropical areas under fluctuating light conditions. Most bamboos are fast growing, reaching their full height and diameter within a single growth season, which indicates that bamboos may possess unique carbon assimilation mechanisms in their leaves. D. latiflorus is an evergreen species locally known as 'tropical giant bamboo', which forms abundant forests in southern China and southeast Asia, and is a valuable natural resource used as food, building material and other human consumption [27]. Consequently, D. latiflorus is an obvious choice for an initial study of miRNAs in bamboo.
Using high-throughput sequencing and bioinformatics analysis, we identified 84 conserved miRNAs belonging to 17 families, and 81 novel miRNAs from more than 11 million raw sequence reads generated from a small RNA library of ma bamboo leaf. One hundred and sixty-two potential targets were identified for the 81 novel bamboo miRNAs. In addition, to confirm the novel predicted miRNAs in bamboo leaf tissues, the expression profiles of 30 novel miRNAs under different light conditions were validated by qRT-PCR. These results will lay the foundation for understanding miRNA-based regulation during D. latiflorus development.

Plant material, RNA isolation and small RNA highthroughput sequencing
Cutting seedlings of ma bamboo (D. latiflorus) were potted in our laboratory under a regime of 16 h light and 8 h darkness at 25°C, with a light intensity of 200 μmol·m -2 ·s -1 and a relative humidity of 75%. As experimental materials, we chose the third piece of new functional leaf (blade tissue only) from the top of the branch, which could be considered a juvenile leaf. The leaves were collected from 2-year-old cuttings and quickly frozen in liquid nitrogen. Total RNA was isolated from leaf tissues using the Trizol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer's instructions.
The small RNA library construction for ma bamboo and Solexa sequencing were carried out at BGI-Shenzhen (Shenzhen, China) using the standard Solexa protocol [28]. Briefly, small RNAs of 15-30 nt in length were first isolated from the total RNA through 15% TBE urea denaturing polyacrylamide gels. Subsequently, 5′ and 3′ RNA adaptors were ligated to these small RNAs, followed by reverse transcription into cDNAs. These cDNAs were amplified by PCR and subjected to Solexa sequencing. After removing low quality reads and trimming adapter sequences, small RNAs ranging from 18-30 nt were collected and used for further analyses. Finally, the selected clean reads were analyzed by BLAST against the Rfam database (http://rfam.sanger.ac.uk/) [29] and the GenBank non-coding RNA database (http:// www.ncbi.nlm.nih.gov/) to discard rRNA, tRNA, snRNA and other ncRNA sequences. In addition, sequencing data within this study were firstly uploaded to NIH Short Read Archive (accession number: SRX347876).

Data content
To obtain and analyze miRNAs from the leaf of ma bamboo, three types of data were used in this study (1). MiRNA data. The miRNA data from the leaf of ma bamboo was obtained by next generation sequencing, and all the reference sequences of mature miRNAs, star miRNAs and their precursors (pre-miRNAs) were downloaded from miRBase, Release 19.0 (http://www.mirbase.org/index.shtml) [13]. (2) Genome data of moso bamboo (Phyllostachys heterocycla var. pubescens) [30]. Although the whole genome sequence of ma bamboo is not available, phylogenomic analyses suggested that ma bamboo and moso bamboo had the closest relationship, with high sequence similarity [31]. For this reason, the sequenced genome of moso bamboo from our previous study was used as the genome to analyze ma bamboo miRNA data (3). Expressed sequence tag (EST) and mRNA data from ma bamboo. The approach outlined in (2) would inevitably produce some biases and some miRNAs would not be identified. Therefore, to fill the gap, ESTs and other mRNA data from ma bamboo were obtained from the NCBI database site and used to identify additional miRNAs.

Prediction of conserved miRNAs, novel miRNAs and potential miRNA targets in ma bamboo
As miRNA precursors have a characteristic fold-back structure, 150 nt of the sequence flanking the genomic sequences of small RNAs (sRNAs) was extracted and used to predict miRNAs. For analysis of conserved miRNAs in ma bamboo, unique sRNAs were aligned with plant mature miRNAs in miRBase Release 19.0. First, after rigorous screening, all retained sequences with three or more copies were considered as potential miRNAs. Second, sRNA sequences with no more than four mismatched bases were selected by BLAST searching against miRBase. Third, the remaining 15-26 nt reads were used to map the genome of moso bamboo and ESTs of ma bamboo using the BLASTN program. Sequences with a tolerance of two mismatches were retained for miRNA prediction. RNAfold (http:// www.tbi.univie.ac.at/RNA/) [32] was used for secondary structure prediction (hairpin prediction) of individual mapped miRNAs, using the default folding conditions to identify the known conserved miRNAs in ma bamboo. Finally, sequences that were not identical to the conserved miRNAs were termed novel miRNAs.
Potential target sequences for the newly identified miRNAs were predicted using the psRNATarget program (http:// plantgrn.noble.org/psRNATarget/) with default parameters. Newly identified miRNA sequences for ma bamboo were used as custom miRNA sequences, while coding sequences for moso bamboo, as well as EST and mRNA databases for ma bamboo, were used as custom plant databases, respectively. All predicted target genes were evaluated by the scoring system and criteria defined in a previous report [33]. Sequences with a total score less than 3.0 were identified as miRNA targets.

Expression analysis of novel miRNAs by qRT-PCR
The stem-loop qRT-PCR method [3] was used to detect the expression levels of novel miRNAs in D. latiflorus. Forward primers were specifically designed for each individual miRNA, as detailed in a previous method [3]: six nucleotides of the 3' end of the stem-loop RT primer were complementary to the 3' end of the mature miRNA, and the sequence 5' -CTCAACTGGTGTCGTGGAGTC -3' was used as the universal reverse primer [34]. U6 snRNA was used as an internal control [35]. More detailed information is supplied in File S1.
The qRT-PCR was conducted using a SYBR Green I Master Kit (Roche, Germany) on a LightCycler ® 480 Real-Time PCR System (Roche). The final volume was 20 μl, containing 7.5 μl 2×SYBR Premix Ex Taq, 0.3 μl of each primer (10 μM), 2 μl of cDNA and 7.2 μl of nuclease-free water. The amplification was carried out as follows: initial denaturation at 95°C for 10 min, followed by 41 cycles at 95°C for 10 s, 55°C for 20 s, and 72°C for 10 s. The melting curves were adjusted as 95°C for 5 s and 55°C for 1 min and then cooled to 40°C for 30 s [36]. All reactions were repeated three times.
For each condition, the qRT-PCR experiments were performed as biological triplicates and expression levels were normalized according to that of the internal control. The relative value of the gene expression was calculated using the 2 -ΔΔCt method [37]. Statistical tests were performed on the qRT-PCR data using SPSS (Statistical Product and Service Solutions) 18.0 software. Error bars representing the standard deviation were derived from the three experiments in triplicate.

Results and Discussion
Overview of small RNA library sequencing Deep sequencing of the small RNA library from D. latiflorus leaf tissues produced 11,513,607 raw sequence reads. Low quality sequences, adapters and small sequences shorter than 16 nucleotides (nt) were removed, leaving 10,593,305 clean reads and 6,320,379 unique sequences. After further removal of unannotated small RNAs and non-coding RNAs, such as tRNAs, rRNAs, siRNAs, snRNAs, snoRNAs and other noncoding RNAs, 910,151 miRNA sequences, accounting for 8.59% of the total sRNA, were identified ( Table 1).
The significant feature of the size profile permitted miRNAs to be distinguished from other small RNAs. The sRNAs length distribution (10-30 nt) of the original reads demonstrated that the most abundant reads were those of 20-24 nt in length (Figure 1), which were consistent with sRNAs with known function [19]. The most abundant class was those of 24 nt, which was also highly consistent with those small RNA of other Poaceae plants based on Solexa sequencing technology, such as rice [38], barley [39] and wheat [40].

Identification of conserved miRNAs in ma bamboo
Conserved miRNAs have been found in many plant species and play significant roles in plant development and stress responses [4]. To identify the conserved miRNAs in D. latiflorus, the sRNA library was searched using BLASTN for unique mature plant miRNA sequences in miRBase Version 19.0, in which miRNAs of bamboo had not been identified. Following a set of strict filtering criteria and further sequence analysis, 54 conserved miRNAs and 30 conserved star miRNAs were identified from D. latiflorus, which belonged to 17 families (Table 2).
Compared with the conserved miRNA families in plants reported by Cuperus et al [41], 15 conserved miRNA families, accounting for 88% of the conserved miRNA families in ma bamboo, were represented in ma bamboo. However, of the 29 conserved miRNA families in ehrhartoideae [41], 14 conserved miRNA families were not represented in ma bamboo. This may have resulted from the data in miRBase database being a different version, may represent the particular evolution and function of miRNAs in ma bamboo, or may have resulted from the inevitable biases in processes such as the preparation of samples, sequencing and analysis. Further experiments using increased amounts of data are required to verify the conserved miRNA families in ma bamboo. In addition, based on the moso bamboo genome, 13 conserved miRNA families were predicted and compared with those experimentally derived from ma bamboo. This analysis indicated that eight conserved miRNA families, accounting for 62% of all families predicted in moso bamboo, were represented in ma bamboo. This may be explained by the fact that they are two different bamboo species with distinct biological and evolutionary features. For example, the number of chromosomes is 68 and 48 in ma bamboo and moso bamboo, respectively.
MiRNAs with high sequencing frequencies have been shown to play fundamental and essential regulatory functions in maintaining biological processes. Therefore, the read counts for known miRNA families were analyzed ( Figure 2). The most abundant miRNA family (594,744 reads) was miR168, which was represented five to 41 times more frequently than other relatively high abundance miRNAs, including miR156/157, miR535, miR165/166 and miR167, whose total abundances ranged from 14,448 to 109,638 reads. Moreover, the top three The number of members in the differently conserved miRNA families was also analyzed. Four families (miR156/157, miR165, miR167 and miR535) contained multiple members, with nine, six, five and six members respectively. Five families, miR390, miR396, miR399, miR827 and miR1878, only had one member.
The bamboo flowering cycle can take up to 120 years and also involves infrequent and unpredictable flowering events, as well as peculiar monocarpic behavior, e.g., flowering once before culm death [26]; therefore, the mechanism of flowering, as a unique characteristic of bamboo, has received much research interest. Studies in other plants indicated that miRNAs play important roles in flower development [42,43]. The transition from juvenile to adult phase is controlled by miR156 and miR172. MiR172 also influences floral organ identity, as evidenced by failure of carpel abortion from the male inflorescence [42]. In addition, extensive analysis of certain monocotyledonous species, such as Brachypodium distachyon, Oryza sativa, Sorghum bicolor, and Zea mays, revealed more than three members in each miR172 family. Therefore, absence of miR172 may contribute to bamboo's special regulatory mechanism of flower development. Other miRNAs (miR164, miR167, miR169 and miR319) with functions during different stages of flower development were detected in D. latiflorus.
As seen in Figure 3, to explore the evolutionary roles of these conserved miRNAs, deep analyses focused on extensive comparisons against known conserved miRNAs in other plant species, including Picea abies, Pinus taeda, Physcomitrella patens, Selaginella moellendorffii, Arabidopsis thaliana, Brassica napus, Ricinus communis, Medicago truncatula, Citrus sinensis, Vitis vinifera, O. sativa, S. bicolor and Z. mays. Based on BLAST searches and sequences analysis, some miRNA families of D. latiflorus (miR156/157, miR160, miR165/166, miR171, miR319, miR390 and miR396) are highly conserved and ancient. For example, miR156/157 is a family of endogenous miRNAs with a relatively high expression level in the juvenile phase of many plants. The level of miR156/157 gradually decreases with plant age [44,45].
However, some miRNA families (miR528, miR535, miR827, miR1318 and miR1878) are less evolutionarily conserved and are involved in regulation of diverse physiological processes. For example, the miR528-target recognition site is only present within monocot genes, while all eudicot genes orthologous to SsCBP1 from Arabidopsis, poplar, grape, and soybean genomes completely lack the miR528-target recognition site, which supports the view that miR528 is a monocot-specific miRNA [46]. Another example is miR535, which was predicted to target a gene encoding a brassinosteroid signaling positive regulator protein in Physocmitrella patens [47]. It was reported that mature miR535 may be considered as a divergent member of the miR156/157 family, and targeting of SPL genes may be a vestigial function of miR535, which may have been performed more efficiently in the course of evolution by other members of the family [48].
The remaining five miRNA families, including miR164, miR167, miR168, miR169 and miR399, are homologous in eudicotyledons and monocotyledons, indicating these miRNA  families may be recent and are more sensitive to certain stresses. For example, miR168 regulates loop ARGONAUTE1 (AGO1) homeostasis, which is crucial for the regulation of gene expression and plant development [49]. Another report demonstrated that NFYA5 contained a target site of miR169, which was downregulated by drought stress through an ABAdependent pathway [50,51]. MiR399 is upregulated in phosphate starvation and its target gene, encoding a ubiquitinconjugating E2 enzyme, is downregulated in A. thaliana [50,51].

Identification of novel miRNAs in ma bamboo
Previous studies have shown that each species has speciesspecific miRNAs [17,18,24,25,33,39,[52][53][54]; therefore, ma bamboo is also likely to have unique miRNAs. In addition, a distinct feature of miRNAs is the ability of their pre-miRNA sequences to adopt the canonical stem-loop hairpin structure. After excluding sRNA reads homologous to known miRNAs and other non-coding sRNAs, the remaining 18-24 nt sRNAs were subjected to secondary structure analysis of their precursors using RNA fold software. Seventy-six novel miRNAs and five novel star miRNAs were identified as possible miRNA candidates in ma bamboo ( Table 3). ESTs of ma bamboo were used to predict miRNAs; however, no novel miRNAs were  Moreover, the length distributions of plant pre-miRNAs are more heterogeneous than animal pre-miRNAs. The shortest pre-miRNA is 53 nt in length (miRBase ID: ath-MIR5645b) and   the longest in 938 nt in length (miRBase ID: aly-MIR858) in miRBASE Version 19. In ma bamboo, the distribution of pre-miRNAs is 54-96 nt, which is consistent with the pre-miRNAs distributions in other plants [17,25,39]. In addition, the length distribution of novel miRNA ranged from 18-24 nt in length. Forty-eight (45.3%) of the novel miRNAs belong to the 24nucleotide class, representing the most abundant novel miRNAs. Among the novel miRNAs, dla-miRC1 had the highest expression in our data, with 4,801 reads. Based on their frequencies and sequences in the small RNA library, although the expression levels of these candidates ranged from thousands of reads to single reads, in general, novel miRNA candidates of D. latiflorus showed lower expression compared with most of the conserved families. The low abundance of novel miRNAs might indicate that these miRNAs play a specific role in certain tissues or developmental stages, and may be considered young miRNAs in terms of evolution [41]. This sRNA library was only generated from leaf tissues, and future experiments will be carried out to determine whether these lowabundant miRNAs are expressed at higher levels in other organs, such as flowers and seeds, or whether they are regulated by certain stresses.

Prediction of miRNA targets in ma bamboo
The identification of miRNA targets using bioinformatics approaches is an essential step to further understand the regulatory function of miRNAs [25,33,55]. Given the lack of genomic data for ma bamboo, the CDS from the genome of the closely related species moso bamboo and ESTs of ma bamboo were used to predict new miRNA target genes, using the criteria described in the material and methods. We identified 176 potential targets, including 139 targets from moso bamboo in File S2 and 37 targets from ma bamboo (Table 4).
Among the novel miRNAs, dla-miRC5 has 26 putative target genes with different functions, which indicates that dla-miRC5 might be involved in regulating the expression of multiple genes in ma bamboo. Some of the novel miRNAs target transcription factor (TF) genes that have been confirmed to play key roles in plant development. The targets of dla-miRC1 are auxin response factor (ARF) TF genes. ARF can regulate auxin response genes by binding specifically to cis-elements of their promoters to affect the developmental process of plants [55]. To date, miR160 and miR167 have been demonstrated to inhibit ARF10, ARF16, ARF17 [8,56] and ARF6, ARF8 [57], respectively; however, dla-miRC1 is barely homologous to  miR160 and miR167, indicating that it might be another specific suppressor of ARF genes in ma bamboo. In addition, an ERF gene (belonging to AP2/EREBP family) is predicted to be a target of dla-miRC37, a MYB TF gene is a putative target of dla-miRC56, and both dla-miRC9 and dla-miRC16 have one target gene belonging to the BHLH family of TFs. Scarecrow, of the GRAS TF gene family, is targeted by dla-miRC33. All the targets of these novel miRNAs were TF genes that function in regulating the development of plants [58][59][60][61], which indicated these miRNAs might play important roles in ma bamboo development and stress responses.
Some novel miRNAs, including dla-miRC2, dla-miRC5, dla-miRC13, dla-miRC35 and dla-miC45, are predicted to target genes related to chloroplast synthesis and photosynthesis in the leaf, which is the most important photosynthetic organ of plants. This result suggested that the miRNAs might be leafspecific and be involved in the process of regulating chloroplast synthesis and photosynthesis.

Expression profiles of novel miRNAs in response to light
The expression patterns of miRNAs could provide clues to their functions [22]. As an efficient and sensitive method for detecting gene expression [62], stem-loop qRT-PCR was developed based on the common qRT-PCR method, with advantages such as increased sensitivity and high accuracy. Therefore, stem-loop qRT-PCR has been widely employed to distinguish two miRNAs with small differences [3].
Among the novel predicted miRNAs, the top 30 novel miRNAs (according to the number of reads) were selected and primers were designed based on their highly specific stem-loop reverse sequences to experimentally identify their expression. As shown in File S3, the result of stem-loop qRT-PCR showed that each of the selected miRNAs had a highly specific dissolution curve, which indicated that the primers could be used for further analysis. Moreover, as shown in Figure 4, the in-depth analysis demonstrated the expression of the miRNAs under high light stress conditions; the number of miRNAs that were upregulated and downregulated was 10 and 16, respectively. The expressions of dla-miRC18, dla-miRC27-5p and dla-miRC27-3p were significantly increased under high light (P<0.01), among which dla-miRC18 was upregulated to 15 times the level of the control. However, the expression of dla-miRC1, dla-miRC19 and dla-miRC28 was significantly downregulated (P<0.01). Under dark condition, the expressions of dla-miRC1, dla-miRC22, dla-miRC25, dla-miRC27-5p and dla-miR27-3p were upregulated significantly (P<0.01), while the expressions of dla-miRC5, dla-miRC17 and dla-miRC29 were downregulated significantly (P<0.01).
This analysis also indicated that the expressions of dla-miRC1 and dla-miRC19 were downregulated significantly under high light and upregulated markedly under dark conditions (P<0.01), indicating that they were inhibited by light. However, dla-miRC29 was upregulated under high light and downregulated under dark, indicating it was light-induced expression. MiR27-5p and miR27-3p were cleaved from the same precursor and had similar expression profiles, being upregulated under both high light and dark, which indicated they might have the same promoter and cis-elements, and synergistic expression characteristics. The expression results indicated that these miRNAs were affected by light, and might play vital roles in the regulation of genes involved in light signal transduction and light stress.
Bamboo converts light energy into chemical energy through photosynthesis, which is one of the necessary processes supplying carbohydrates for the rapid expansion of cells. As the first comparative identification of miRNAs among leaves treated with distinct light conditions, this analysis showed that bamboo may possess a unique light regulation mechanism involving miRNAs, although its function and mechanism are currently unknown.

Conclusions
We identified 81 novel miRNAs and 84 conserved miRNAs belonging to 17 families in the leaf of ma bamboo using highthroughput sequencing and bioinformatics analysis. The results of qRT-PCR indicated the miRNAs might regulate the expression of genes involved in photosynthesis, which acts as a key metabolic pathway in the fast growth of bamboo. These miRNAs will add to the growing database of new miRNAs and lay the foundation for further understanding of miRNA function in the regulation of ma bamboo development and other biological characteristics. These miRNAs identified in the leaf of ma bamboo provide new opportunities for future functional genome research in bamboo and other related species.

Supporting Information
File S1. Primers used in novel miRNAs qRT-PCR.