High-throughput sequencing and characterization of the small RNA transcriptome reveal features of novel and conserved microRNAs in Panax ginseng.

microRNAs (miRNAs) play vital regulatory roles in many organisms through direct cleavage of transcripts, translational repression, or chromatin modification. Identification of miRNAs has been carried out in various plant species. However, no information is available for miRNAs from Panax ginseng, an economically significant medicinal plant species. Using the next generation high-throughput sequencing technology, we obtained 13,326,328 small RNA reads from the roots, stems, leaves and flowers of P. ginseng. Analysis of these small RNAs revealed the existence of a large, diverse and highly complicated small RNA population in P. ginseng. We identified 73 conserved miRNAs, which could be grouped into 33 families, and 28 non-conserved ones belonging to 9 families. Characterization of P. ginseng miRNA precursors revealed many features, such as production of two miRNAs from distinct regions of a precursor, clusters of two precursors in a transcript, and generation of miRNAs from both sense and antisense transcripts. It suggests the complexity of miRNA production in P. ginseng. Using a computational approach, we predicted for the conserved and non-conserved miRNA families 99 and 31 target genes, respectively, of which eight were experimentally validated. Among all predicted targets, only about 20% are conserved among various plant species, whereas the others appear to be non-conserved, indicating the diversity of miRNA functions. Consistently, many miRNAs exhibited tissue-specific expression patterns. Moreover, we identified five dehydration- and ten heat-responsive miRNAs and found the existence of a crosstalk among some of the stress-responsive miRNAs. Our results provide the first clue to the elucidation of miRNA functions in P. ginseng.

Introduction

Panax ginseng C.A. Meyer, belonging to the genus Panax of the family Araliaceae, is an economically significant medicinal plant species and has attracted a lot of attention for the phytomedicinal use since ancient times [1]. Evidence from laboratory and clinical trials showed that the roots of P. ginseng could play significant roles in immune system modulation, anti-stress, anti-cancer and anti-diabetes in human beings [2]–[6]. The major bioactive ingredients of ginseng roots are triterpene saponins, which are known as ginsenosides. To date, more than 30 ginsenosides have been identified [7]. Among them, Rb1, Rb2, Rc, Rd, Re, Rf, Rg1 and Rg2 are most relevant to the pharmacological effects of ginseng roots [7]–[9]. In addition to ginsenosides, ginseng roots also contain some other bioactive compounds, such as polysaccharides, phenolics, and flavonoids [10]–[12].

P. ginseng is a part shade, low-growing perennial plant native to Asian [7]. It is a tetraploid plant species (2n = 4x = 48) with a large genome. Although the whole genome of P. ginseng has not been decoded, the whole chloroplast genome and a large number of expressed sequence tags (ESTs) have been sequenced [13]–[17]. The cDNA and gene sequences currently available includes 31,741 unique sequences assembled from 217,519 high quality 454 sequencing reads [13], 11,412 ESTs obtained from normal cDNA library, and about 563 other cDNA and gene sequences (http://www.ncbi.nlm.nih.gov/). The genetic transformation system of P. ginseng has been developed and several ginsenoside biosynthesis-related genes have been characterized through this system. For example, overexpression of the squalene synthase gene PgSS1 in adventitious roots of transgenic P. ginseng resulted in a considerable increase of phytosterol and triterpene saponin (ginsenoside) content [4]. RNA interference (RNAi) of the squalene epoxidase gene PgSQE1 and the dammarenediol synthase gene DDS in transgenic P. ginseng led to the reduction of ginsenoside production [18], [19]. It provides great information and tools for studying gene functions and regulatory networks in P. ginseng.

microRNAs (miRNAs) are a class of small endogenous non-coding RNAs. They originate from long primary transcripts (pri-miRNAs), which have stem-loop structures [20]. In plants, miRNAs regulate target genes either through direct cleavage of transcripts [21], [22] or translational repression at post-transcriptional level [23], or, in some case, by methylation at transcriptional level [24]. It has been shown that miRNAs play important regulatory roles in plant development and stress responses [25]–[30]. For example, miR172 regulates the development of flowers by targeting APETALA2 (AP2) and its homologs with two tandem AP2 DNA-binding domains [31]–[33]. miR160 is involved in root development through the regulation of auxin response factors (ARFs) [32]–[34]. An increasing number of stress-responsive miRNAs are being identified [28], [30], [35].

Since the first discovery of plant miRNAs from Arabidopsis, plant miRNAs have been intensely studied using experimental and computational approaches [36], [37]. The number of miRNAs identified has reached to 4011 according to miRBase release 18 [38]. However, there is no information available for miRNAs from P. ginseng. In this study, we constructed a small RNA library of P. ginseng and carried out high-throughput Illumina sequencing. We identified a total of 101 conserved and non-conserved miRNAs from the small RNA dataset. Computational analysis of miRNAs and their precursor sequences revealed some features of P. ginseng miRNAs. Many non-conserved P. ginseng miRNAs showed tissue-specific expression and some of them responded to dehydration and heat treatments. Consistently, many predicted targets of non-conserved P. ginseng miRNAs are associated with metabolism, signal transduction and stress responses.

Results

Features of P. ginseng Small RNA Population

In order to identify small RNAs in P. ginseng, we generated a small RNA library using pooled RNA isolated from roots, stems, leaves and flowers. The library was then sequenced by the high-throughput Illumina sequencing technology [39]. A total of 13,326,328 raw sequence reads were obtained. After removing adaptors, low quality sequences and those smaller than 18 nt, 12,000,591 clean reads with sizes ranged from 18 to 30 nt were collected (Table S1). They were represented by 1,502,664 unique small RNA sequences, of which 114 were sequenced more than 10,000 times, whereas 1,178,394 (78.4%) and 152,436 (10.1%) were sequenced only one and two times, respectively. It indicates that a large, diverse and highly complicated small RNA population exists in P. ginseng. Additionally, 30,240 small RNAs, representing 1,283,459 reads, could be mapped to the chloroplast genome of P. ginseng using the SOAP2 program with no mismatches allowed [16], [40], indicating some small RNAs are probably originated from the chloroplast genome.

Among 12,000,591 clean reads, the majority (71.22%) range from 20 to 24 nt in length, wherein the 21 (20.24%) and 22 (18.03%) nt ones are two major size groups (Figure 1A). The group of 24 nt clean reads also accounts for a relatively high proportion (13.54%), whereas it is significantly less compared with those identified from Arabidopsis [41], Medicago truncatula [42], Brachypodium distachyon [43], Arachis hypogaea [44], Citrus trifoliate [45], Gossypium hirsutum [46], Oryza sativa [47] and Zea may [48], where the 24 nt group is the most abundant. Therefore, we analyzed the pattern of low molecular weight (LMW) RNA from the pooled total RNA using gel electrophoresis as previously described [49]. Unexpectedly, the results of electrophoresis show that the most abundant group of small RNAs is 24 nt (Figure 1B). The inconsistency of results from electrophoresis and high-throughput sequencing was also observed in grape, in which a nearly equal number of 21 nt and 24 nt small RNAs was detected by electrophoresis, whereas the number of 21 nt sequence reads are more than five times of the 24 nt reads [50]. The factors causing the inconsistency remain to be elucidated.

Figure 1

Analysis of P. ginseng small RNA population.

(A) Size distribution of clean reads. (B) Electrophoretic analysis of small RNAs. Low molecular weight (LMW) RNA was resolved on a 15% denaturing polyacrylamide gel and stained with ethidium bromide. The sizes of RNA ladder are indicated. (C) Size distribution of unique small RNA sequences. (D) Percentage of unique sequences starting with an A, C, G, and U in each size group.

Size distribution analysis of the 1,502,664 unique small RNA sequences showed that the group of 24 nt small RNAs was the biggest, which accounted for 38.24% of total unique sequences (Figure 1C). It suggests that 24 nt small RNAs are the most diverse in P. ginseng, although their average of abundance is less compared with those in many other plant species. The groups of 21 and 22 nt have the highest percentage of small RNAs beginning with a U nucleotide, whereas the 24 nt group enriches with sequences having an A at the 5′ end (Figure 1D). It indicates a great proportion of miRNAs and trans-acting small interfering RNAs (ta-siRNAs) in the 21 and 22 nt groups and significant amount of repeat-associated siRNAs (ra-siRNAs) and natural antisense transcript-derived siRNAs (nat-siRNAs) in the 24 nt group, since miRNAs and ta-siRNAs tend to begin with a U and ra-siRNAs and nat-siRNAs prefer an A at the 5′ end [51]–[55]. It is consistent with features of small RNAs from other plant species, such as Arabidopsis and Populus [54], [56]. The diversity of 24 nt ginseng small RNAs may indicate the complexity of P. ginseng genome, although the whole genome has not been decoded.

Identification of Conserved miRNAs

Using the SOAP2 program, we mapped 1727 of 1,502,664 unique small RNA sequences to known miRNA precursors with no more than 2 mismatches [38], [40]. Among them, 661 and 41 unique sequences, representing 305,744 and 406 reads, were mapped to mature miRNA and miRNA* regions, respectively, whereas the other 1025 sequences, representing 126,128 reads, were mapped to other regions of known miRNA precursors. Further alignment of the small RNAs mapped to miRNA regions of known precursors allowed us to group the 661 sequences into 73 miRNA candidates, members of 33 conserved miRNA families (Table S2). Similarly, the 41 unique sequences mapped to miRNA* regions were grouped into 22 miRNAs* putatively produced from 9 gene families, including MIR156, MIR166, MIR167, MIR171, MIR172, MIR390, MIR396, MIR482, and MIR4376 (Table S3).

Among the 33 identified miRNA families, 22 are deeply conserved among many plant species, whereas the other 11 have less evolutionary conservation. For instance, MIR156, MIR166, MIR169, MIR319 and MIR394 have been identified from more than 40 monocots and dicots, whereas MIR403 is found only in several dicots, such as Arabidopsis, P. trichocarpa, Vitis vinifera, and Ricinus communis [30], [38], [57]–[59]. The other 10 miRNA families with less conservation include MIR1439, MIR5072 and MIR5079 that were previously identified from rice only, MIR1509 and MIR1510 found in soybean and Medicago, MIR894 and MIR1024 identified from Physcomitrella patens, MIR1863 found in rice and Norway spruce, MIR4376 found in soybean and tomato, and MIR5139 previously identified from a medicinal plant species, Rehmannia glutinosa [42], [60]–[72]. The identification of less conserved miRNAs indicates the specificity and complexity of P. ginseng miRNAs.

It has been suggested that the expression profiles of miRNA genes can be estimated by high throughput sequencing-based technology [73], [74]. Among the 33 conserved miRNA families, MIR166 was sequenced 228,538 times, which accounted for 74.8% of total conserved miRNA reads, suggesting it is the most abundant miRNA family in the P. ginseng tissues analyzed (Figure 2). The second abundant miRNA family is MIR159, which was sequenced 61,264 times and accounted for 20.0% of total conserved miRNA reads. The other conserved miRNA families are less abundant and each has less than 1.5% of total conserved miRNA reads. Among them, 15 have less than 10 reads in our small RNA data set (Table S2). Compared with other plants, P. ginseng has significantly different expression profiles for some miRNAs. For example, MIR156 is highly expressed in peanut [44], whereas its expression in P. ginseng is low (Figure 2). These results suggest the differential expression of miRNAs in P. ginseng and indicate the specificity of miRNA expression in different plant species.

Figure 2

Number of reads for the conserved miRNA families.

The read numbers of MIR159 and MIR166 are shown.

We mapped the 73 conserved miRNAs to unigenes of P. ginseng using the SOAP2 program and then analyzed secondary structures of the unigenes using the mfold program [40], [75]. As a result, four miRNA precursor sequences were identified. Two of them are precursors of MIR482. The other two generate miR2118 and miR4376 (Figure S1). Further mapping of the small RNAs to miRNA precursors allowed the identification of three additional miRNA* sequences: one for miR482a, one for miR4376, and the other for both miR482b and miR2118 (Table S3). No precursor sequences for other conserved miRNAs were identified. It could be due to the limited number of P. ginseng unigenes.

Identification of Non-conserved miRNAs

In order to identify non-conserved miRNAs as much as possible, we used two computational approaches. First, we mapped small RNAs to 31,371 unigenes of P. ginseng using the SOAP2 program with no mismatches allowed [40] and analyzed the secondary structures of unigenes with small RNA mapped using the mfold program [75]. Next, we used MIREAP software (http://sourceforge.net/projects/mireap/) to predict precursors from P. ginseng unigenes. In total, we identified 25 miRNA precursors. These precursors produced a total of 28 mature miRNAs. We then submitted the precursors and the corresponding mature miRNAs to miRBase for official names [38]]. Based on sequence similarities, the 25 miRNA precursors were grouped into 9 miRNA gene families, although some members of a family generated distinct mature miRNA sequences (Table 1, Figure S1). Except MIR6135 and MIR6140 families, which have 11 and 4 members, respectively, the other 7 non-conserved miRNA gene families have only one or two members. The results showing few members for non-conserved miRNA families have been previously found in other plant species, such as Arabidopsis [76]. In addition to miRNAs, we identified a miRNA* sequence for the non-conserved miRNA, miR6136b, which verified computational prediction of non-conserved miRNAs (Table 1).

Table 1

Non-conserved miRNAs and miRNAs* in Panax ginseng.

miRNA name	miRNA sequence (5′→3′)	Size (nt)	reads	MFEIa		Unigene ID		Length (bp)
miR6135a	TGGTAAGTTGGTCAATTGGC	20	107		1.49		PUT-183a-Panax_ginseng-10400		390
miR6135b	TGGTAAGTTGGTCAATTGGC	20	107		1.48		PUT-183a-Panax_ginseng-1888		3603
miR6135c	GGGTAAGTTGGTCAATTGAC	20	17		1.45		PUT-183a-Panax_ginseng-1888		3603
miR6135d	GGGTAAGTTGGTCAATTGAC	20	17		1.66		PUT-183a-Panax_ginseng-5770		779
miR6135e.1	GGGTAAGTTGGTCAATTGAC	20	17		1.47		PUT-183a-Panax_ginseng-7458		799
miR6135e.2	AATTGACTAATAGAATACTGACAC	24	21		1.47		PUT-183a-Panax_ginseng-7458		799
miR6135f	GAGTAAGTTGGTCAATTGGC	20	4		1.32		PUT-183a-Panax_ginseng-1153		514
miR6135g	GAGTAAGTTGGTCAATTGGC	20	4		1.63		PUT-183a-Panax_ginseng-11661		155
miR6135h	GTAAGTTGGTCAATTGGC	18	5		1.02		PUT-183a-Panax_ginseng-17044		767
miR6135i	AATTGGCCAATAGAATACTGACAC	24	30		1.43		PUT-183a-Panax_ginseng-5770		779
miR6135j	AATTGACTAATAGAATACTGACAC	24	21		1.32		FW1NBNE01BCPRY		474
miR6135k	CGTGTCGATACTGTATTGGT	20	3		1.54		PUT-183a-Panax_ginseng-1587		475
miR6136a.1	TAGACGACGGTTGTATGACCG	21	32		1.40		FW1NBNE01BEGNB		465
miR6136a.2	ACGGGTGAGTAAGATAAGGGGTAT	24	88		1.40		FW1NBNE01BEGNB		465
miR6136b	TCATACAACCGTCGTCTATAC	21	70		1.45		FW1NBNE01BEGNB		465
miR6136b*	TATAAATGATGGTTGTATGAC	21	1		1.45		FW1NBNE01BEGNB		465
miR6137a	ATGAAAATTGTCGCTATAGATC	22	3		2.13		FW1NBNE01A9PXQ		425
miR6137b	ATGAAAATTGTCGCTATAGATC	22	3		1.56		PUT-183a-Panax_ginseng-4023		606
miR6138	TACGTTTGGATTGAAGGAATGAAA	24	13		0.98		PUT-183a-Panax_ginseng-10763		517
miR6139	AAGAATCATTGGGAAGGGAAGAAA	24	9		1.15		FW1NBNE01AHR0J		608
miR6140a	AATGTTTGTAGAATAGTTTGTGTC	24	3		1.38		FW1NBNE01A9GRG		391
miR6140b	AATGTTTGTAGAATAATTTGTGTA	24	2		1.14		PUT-183a-Panax_ginseng-19298		1417
miR6140c	GCTGAGGTGGAGTATGCCACATC	23	5		1.32		FW1NBNE01ATMKX		436
miR6140d	CGTTGATGTGGCATACTTCACC	22	3		1.12		PUT-183a-Panax_ginseng-8684		972
miR6141	TAACTAAATCTGGCCTGTAGCGGA	24	7		1.56		FW1NBNE01B00A0		391
miR6142	GACGATTTTTTGGGCTATGACGAC	24	5		0.93		FW1NBNE01CDN8L		235
miR6143a	AGTACTGTATTGGGCATGAAG	21	13		1.20		FW1NBNE01CELC5		535
miR6143b-5p	ACAATGTCGACACGCAGGCGGAGA	24	2		1.24		PUT-183a-Panax_ginseng-16660		399
miR6143b-3p	CAGCACTGTATTGAACATGAA	21	4		1.24		PUT-183a-Panax_ginseng-16660		399

MFEI, minimal folding free energy index of the hairpin structures.

Special Features of miRNA Genes in P. ginseng

The precursors of miR482a and miR2118 cluster in a unigene, PUT-183a-Panax_ginseng-17125 (Figure 3A). Existence of two precursors is verified by the finding of miR482a* and miR2118* in our small RNA data set. Clustering of miRNAs is a universal phenomenon in animal [77]–[79]. It has also been observed for some plant miRNAs, such as some members of the MIR156 family in wheat and tobacco [80], [81], some members of the MIR166 family in Arabidopsis, rice, wheat and soybean [21], [81], [82], moss miR1219a and miR1219b [83], some member of the MIR395 family in Arabidopsis and wheat [41], [81], switchgrass miR2118a and miR2118b [84], loblolly pine miR950a and miR950b [85], and Pinus pinaster miR1314a and miR1314b [73]. Both miR482a and miR2118 are located in the 3′ arm of hairpin structures and are unidirectional in the unigene (Figure 3A). The presence of two miRNA precursors in a unigene indicates that the generated miRNAs are co-regulated and are probably involved in the network of a cellular process [77]–[79]. Consistently, miR482a and miR2118 show a high sequence similarity (Figure 3B), which enable them to regulate same targets and/or different members of a gene family in P. ginseng.

Figure 3

The MIR482a and MIR2118 precursors cluster in a unigene.

(A) The secondary structure of PUT-183a-Panax_Ginseng-17125 was predicted by the mfold program using the default parameters [75]. Arrows indicate the position and orientation of miRNAs and miRNA*s. (B) Alignment of the mature miR482a and miR2118 sequences. Vertical dashes indicate the Watson-Crick pairing between miR482a and miR2118 sequences.

Several miRNA precursors generate small RNAs from distinct regions (Figure S1). For instance, both miR6136a.1 and miR6136a.2 are derived from FW1NBNE01BEGNB. PUT-183a-Panax_ginseng-7458 produces two putative miRNAs, including miR6135e.1 and miR6135e.2. Similarly, PUT-183a-Panax_ginseng-16660 generates two mature miRNAs, miR6143b-5p and miR6143b-3p. miR6135e.1 and miR6135e.2 are located in the 5′ arm of a hairpin structure, miR6136a.1 and miR6136a.2 are located in the 3′ arm, whereas miR6143b-5p and miR6143b-3p are produced from different arms of a hairpin structure. Generation of multiple miRNAs has previously been observed for several conserved miRNA genes, such as Arabidopsis MIR159, MIR319, MIR163, MIR822, MIR829 and MIR839, rice MIR159 and MIR319 [50], [52], [86], [87]. Both MIR159 and MIR319 precursors generate multiple miRNAs in a phased loop-to-base manner [88], [89]. However, mechanisms for the production of multiple miRNAs from a precursor remain to be elucidated in P. ginseng.

Three unigenes generating non-conserved miRNAs could be bidirectionally transcribed (Table 1, Figure S1). It includes PUT-183a-Panax_ginseng-1888, PUT-183a-Panax_ginseng-5770, and FW1NBNE01BEGNB. PUT-183a-Panax_ginseng-1888 produced miR6135b, whereas its antisense transcript was predicted to generate miR6135c. Similarly, PUT-183a-Panax_ginseng-5770 and FW1NBNE01BEGNB are precursors of miR6135d and two members of MIR6136 family (miR6136a.1 and miR6136a.2), whereas their antisense transcripts produce miR6135i and miR6136b, respectively. Bidirectionally transcribed miRNA precursors were previously reported in various animals and plants, such as Drosophila [90]–[92], soybean [82], and switchgrass [84]. Because the transcription direction of these ginseng unigenes is unknown, it is possible that one or both of the predicted precursors could be authentic.

Conserved and Non-conserved Targets of Conserved miRNAs

Plant miRNAs show perfect or near-perfect complementarities to their targets [93], [94]. It allows an effective prediction of miRNA targets through computation. To predict the targets of miRNAs in P. ginseng, we used psRNATarget, a useful web server, which searches potential targets using an improved iterative parallel Smith-Waterman algorithm and a weighted scoring schema and calculates penalty scores for mismatched patterns in the miRNA:mRNA duplexes within a 20-base sequence window [95], [96]. With the application of penalty score cutoff threshold 0–3, we identified a total of 99 unigenes from 31,371 assembled P. ginseng unigenes to be targets of 28 conserved miRNA families (Table S4). The number of targets for each miRNA family ranges from one to eight. Each unigene is predicted to be the target of a conserved miRNA except the TIR-NBS disease resistance protein gene, FW1NBNE01A3PPX, which is targeted by pgi-miR482 and pgi-miR1510. No targets were found for MIR160, MIR168, MIR319, MIR894, and MIR1024 in the unigene set. It could be due to the limited number of P. ginseng unigenes or the low expression levels of their targets in the P. ginseng materials used for EST library construction.

Twenty six of the 99 predicted targets are conserved among various plant species, suggesting the conserved regulatory roles of some miRNAs in plants (Table S4). Most of the conserved targets are transcription factor genes, including six squamosa promoter-binding-like protein genes (SPLs) regulated by miR156, four miR167-targeted auxin response factor genes (ARFs), three GRAS family transcription factor genes regulated by miR171, two homeobox-leucine zipper protein genes cleaved by miR166, and each of the NAC-domain protein gene, CCAAT-box binding factor gene, AP2 domain-containing protein gene and growth-regulating factor gene targeted by miR164, miR169, miR172, and miR396, respectively. In addition to transcription factor genes, the other seven conserved targets include four F-box family protein genes regulated by miR393 or miR394, an ATP sulfurylase gene cleaved by miR395, a laccase gene targeted by miR397, and a miR2118-regulated TIR-NBS-LRR resistance protein gene. Using the modified 5′ RLM-RACE method, we experimentally validated eight conserved targets of six miRNA gene families in P. ginseng, confirming the functional conservation of some miRNAs (Figure 4, Table S4).

Figure 4

Validation of the predicted mRNA targets.

The cleavage sites were determined by 5′ RLM-RACE. The nucleotide positions of the miRNA complementary sites of unigenes are indicated. The unigene sequence of each complementary site from 5′ to 3′ and the cloned miRNA sequence from 3′ to 5′ are shown. Watson-Crick pairing (vertical dashes) and G:U wobble pairing (circles) are indicated. Vertical arrows indicate the 5′ termini of miRNA-guided cleavage products, as identified by 5′ RLM-RACE, with the frequency of clones shown.

Consistent with previous reports showing species-specific regulatory roles of many conserved miRNAs in various plant species, such as Populus trichocarpa, Pinus taeda and Digitalis purpurea [28], [35], [85], [97], seventy three of the 99 predicted miRNA targets are not conserved in other plants (Table S4). Among them, sixty could be important in plant development and defense responses, whereas the other thirteen are function-unknown. Using the modified 5′ RLM-RACE method, we experimentally validated contig2586, which encodes the Ca2+-transporting ATPase, to be a target of MIR4376. It suggests that conserved miRNAs may play species-specific functions by the cleavage of non-conserved targets.

Targets of Non-conserved miRNAs

Using the computational approach, we predicted 31 targets for seven of nine non-conserved miRNA gene families (Table 2). Many of the predicted targets are associated with metabolism, signal transduction and stress responses, and about one-third of the targets encode proteins with unknown functions. These results are consistent with our previous findings from P. trichocarpa and confirm that a significant fraction of non-conserved miRNAs are part of the regulatory networks associated with specific growth conditions or developmental processes [35]. Interestingly, non-conserved miR6139 and miR1439 conserved between rice and P. ginseng were predicted to target the ZIP metal ion transporter gene (PUT-183a-Panax_ginseng-20083) for cleavage (Table 2 and Table S4), suggesting both miR6139 and miR1439 are probably involved in the homeostasis of metal ion and play a role in plant response to salt stress [98], [99]. Genes to be targeted by more than one miRNAs have been previously observed in other plant species [28]. For instance, five disease resistance protein genes were found to be targets of miR472, miR482 and/or miR1448, and fourteen pentatricopeptide repeat-containing protein (PPR) genes were regulated by both miR475 and miR476 in P. trichocarpa [28]. Because of the capability in cleaving the same target gene, miR6139 and miR1439 are probably members of a regulatory network.

Table 2

Targets of non-conserved miRNAs in Panax ginseng.

miRNA family	Target function	Unigene ID (penalty score)
MIR6135	Zinc finger CCCH-type antiviral protein	FW1NBNE01BOG42(1.5)
	β-Cyanoalanine synthase	PUT-183a-Panax_ginseng-3189(2.5)
	Mapkkk3	PUT-183a-Panax_ginseng-5059(3)
	Unknown	FW1NBNE01CDUDE(1.5), Contig2012(2.5), PUT-183a-Panax_ginseng-2440(3)
MIR6136	Unknown	PUT-183a-Panax_ginseng-1090(3), PUT-183a-Panax_ginseng-11673(3)
MIR6137	Armadillo/beta-catenin-like	PUT-183a-Panax_ginseng-3440(2.5)
	Protein phosphatase	PUT-183a-Panax_ginseng-4655(2.5)
	Repeat-containing protein	PUT-183a-Panax_ginseng-8005(3)
	Unknown	Contig872(3)
MIR6138	Asparagine synthetase	FW1NBNE01BSZXK(2)
MIR6139	Serine/threonine-protein phosphatase	FW1NBNE01BWTI7(2.5)
	ZIP transporter	PUT-183a-Panax_ginseng-20083(2.5)
	3-Hydroxybutyryl-CoA dehydratase	Contig1795(3)
	UDP-sugar transporter	PUT-183a-Panax_ginseng-9068(3)
	Unknown	FW1NBNE01B23TH(2.5), FW1NBNE01AQA9Q(2.5)
MIR6140	Phytanoyl-CoA dioxygenase-like protein	PUT-183a-Panax_ginseng-6406(2.5)
	KH domain-containing protein	PUT-183a-Panax_ginseng-11439(2.5)
	Casein kinase 2 subunit beta	PUT-183a-Panax_ginseng-10659(2.5)
	Villin	PUT-183a-Panax_ginseng-2502(3)
	Mitochondrial chaperonin-60	PUT-183a-Panax_ginseng-18641(3)
	DnaJ-like protein	FW1NBNE01CCFHV(3)
	Phytochrome	FW1NBNE01A44YJ(3)
	Unknown	PUT-183a-Panax_ginseng-19392(3), PUT-183a-Panax_ginseng-13359(3)
MIR6143	Homeobox-leucine zipper protein	Contig703(2.5)
	Acyl-CoA dehydrogenase	FW1NBNE01BR56P(3)
	Unknown	Contig1700(3)

Tissue-specific Expression of miRNAs in P. ginseng

In order to further predict the functions of miRNAs in P. ginseng, we analyzed the expression of mature miRNAs with known precursor sequences, including conserved miR482a/b, miR2118 and miR4376 and all of non-conserved miRNAs, in leaves, stems, flowers, roots, and embryogenic calli using the poly(A) miRNA-based qRT-PCR method as described [28], [35], [100], [101]. This method verifies the specificity of PCR products by single peak of dissociation curves and can readily discriminate the expression of mature miRNAs having as few as only one nucleotide difference [101]. As a result, the majority of analyzed mature miRNAs were specifically PCR-amplified and many of them showed tissue-specific expression patterns (Figure 5). For example, miR6135i, miR6143a and miR6143b-5p exhibited the highest expression in embryogenic calli, indicating their potential roles in embryogenic callus development and/or maintenance. On the contrary, miR4376, miR6135e,2/j, miR6136a.1, miR6136a.2, miR6136b, miR6140a, miR6140c, miR6140d and miR6142 showed low expression in embryogenic calli, whereas their expression was relatively high in stems and/or flowers. Differential expression was also found for many other miRNAs analyzed (Figure 5). It indicates tissue-specific regulation of miRNA expression in P. ginseng.

Figure 5

Tissue-specific expression of miRNAs in P. ginseng.

Fold changes of miRNAs in leaves (L), stems (S), flowers (F), roots (R) of five-year-old P. ginseng plants and embryogenic calli (C) induced from the cotyledons of P. ginseng seeds are shown. miRNAs were analyzed using the poly(T) adaptor RT-PCR method. miRNA levels in leaves were arbitrarily set to 1.

miR482a/b, miR2118 and miR6139 showed very similar expression patterns with the highest in stems, followed by flowers and roots (Figure 5). Among them, miR482a and miR2118 are generated from two precursors clustering in the same unigene PUT-183a-Panax_ginseng-17125 (Figure 3A). Similar expression patterns were previously observed for other clustered miRNAs that regulated same genes or different transcripts encoding functionally related proteins [77], [102], [103]. Consistently, both miR482 and miR2118 were predicted to target genes encoding disease resistance proteins (Table S4). It indicates that miR482a and miR2118 are involved in relevant cellular processes in P. ginseng.

miR6135i and miR6135e.2/j are members of a miRNA family. However, their expression patterns are distinct. miR6135i showed the highest expression in embryogenic calli, the lowest expression in roots, and moderate expression in stems and flowers, whereas miR6135e.2/j exhibited the highest expression in flowers, followed by stems, leaves and roots, and showed the lowest expression in embryogenic calli (Figure 5). It suggests that miR6135i and miR6135e.2/j could be functionally different in P. ginseng, although they belong to the same gene family. Interestingly, the expression patterns were also different between miR6136a.1 and miR6136a.2, two miRNAs generated from distinct regions of a precursor, FW1NBNE01BEGNB (Figure 5). Similarly, distinct expression patterns were observed for miR6143b-5p and miR6143b-3p produced from the unigene PUT-183a-Panax_ginseng-16660 (Figure 5). These results indicate that the level of miRNAs is not only associated with the transcription of primary miRNAs but also related to the maturation process and/or degradation rate of mature miRNAs. The expression pattern of miR6136a.2 is similar to miR6136b, a miRNA generated from the antisense transcript of FW1NBNE01BEGNB (Figure 5). It implies the existence of functional relationship between miR6136a.2 and miR6136b.

Responses of Non-conserved miRNAs to Abiotic Stresses

In vitro cultured embryogenic callus has been proved to be a great model for studying the role of miRNAs in meristem development, embryogenesis and abiotic stress responses in various plant species, such as rice [104], [105], Japanese Larch [106], and sweet orange [107]. To test whether non-conserved miRNAs are involved in plant response to environmental stresses, we analyzed the expression patterns of non-conserved miRNAs in embryogenic calli treated with dehydration and heat, two significant stressors experienced during the growth and development of P. ginseng, using the miRNA-specific ploy(T) adaptor RT-PCR method [101]. As a result, we identified five dehydration- and ten heat-responsive miRNAs, which showed more than two-fold changes between at least two time-points treated with dehydration and heat, respectively (Figures 6 and 7A). Among them, miR6135e.2/j, miR6135i, miR6138, miR6140a and miR6143b-3p responded to heat treatment only, whereas the other five, including miR6136b, miR6135k, miR6139, miR6140d, and miR6141, responded to dehydration and heat. It suggests the existence of a crosstalk among some non-conserved miRNAs in response to dehydration and heat stresses.

Figure 6

Responses of non-conserved miRNAs to dehydration stress.

Fold changes of non-conserved miRNAs in ginseng embryogenic calli treated with dehydration for 0, 1, 3, 6, 12 and 24 h are shown. miRNAs were analyzed using the poly(T) adaptor RT-PCR method. miRNA levels in calli treated for 0 h were arbitrarily set to 1.

Figure 7

Responses of non-conserved miRNAs to heat stress.

Fold changes of non-conserved miRNAs in ginseng embryogenic calli treated with heat (37°C) for 0, 1, 3, 6, 12 and 24 h are shown. miRNAs were analyzed using the poly(T) adaptor RT-PCR method (A) and the stem-loop RT-PCR method (B). miRNA levels in calli treated for 0 h were arbitrarily set to 1.

Due to low expression levels of non-conserved miRNAs, we used the other precise and sensitive miRNA-specific stem-loop RT-PCR method to further verify the expression patterns of ten heat-responsive miRNAs [108], [109]. The results suggest that the expression patterns of miRNAs from stem-loop RT-PCR and poly(T) adaptor RT-PCR are largely consistent (Figure 7B), validating the ploy(T) adaptor RT-PCR method in analysis of miRNA expression. Moreover, we cloned and sequenced some PCR amplicons. All of the miRNAs analyzed were found to be specifically amplified and had integral 3′ ends (Table S5). These results verify the specificity of PCR-amplification and suggest the integrity of miRNAs in the samples analyzed.

Among the five dehydration-responsive miRNAs, miR6315k was down-regulated after stress. The expression of miR6136b, miR6139, miR6140d and miR6141 was fluctuated (Figure 6). Similarly, different expression patterns were observed for ten heat-responsive miRNAs. miR6135i and miR6136b were induced after heat treatment for 6 h. miR6135e.2/j was highly up-regulated at 12 h after stress. However, the expression of miR6135k, miR6138, miR6139 and miR6141 was down-regulated after treatment. The expression of miR6140a, miR6140d and miR6143-3p was fluctuated in response to heat stress (Figure 7). These results imply that the underlying biological functions vary among stress-responsive miRNAs.

Analysis of Other Small RNAs

In order to characterize other small RNAs, we mapped the remaining 1,502,563 unique small RNA sequences to 31,371 P. ginseng EST sequences using the SOAP2 program with no mismatches allowed and then searched for EST sequences with small RNA mapped [40]. Among the 31,371 ESTs, 14,743 have corresponding small RNAs. 77.65% of 14,743 generated less than five small RNAs. About 2% (294 ESTs) produced more than 50 unique small RNAs. It includes 18 with more than 400 unique small RNAs mapped. Among them, eleven produced from rRNAs, three from chloroplast genome, and one from mitochondrion DNA. The other three were predicted to encode ginsenoside biosynthesis-related proteins [110]. It suggests the complexity of small RNAs in P. ginseng. Full-length cDNA cloning and transgenic analysis may help to further characterize the small RNAs.

The availability of whole genome sequence of P. ginseng chloroplast enabled us to further characterize the small RNAs derived from chloroplast [16]. Blast analysis of 1,502,563 unique small RNA sequences against the chloroplast genome showed that 45,971 small RNAs, representing 2,520,635 clean reads, generated hits over the chloroplast genome. Size distribution of clean reads revealed the 22 nt group, which represents 49.69% of the 2,520,635 reads, to be the biggest size group of chloroplast-derived small RNAs (Figure 8A). The reads of the 24 nt group account for 2.11% of the 2,520,635 reads, which is significantly less compared with the percentage of the 24 nt group in the total 12,000,591 clean reads (Figure 1A). Size distribution analysis of the 45,971 unique small RNA sequences showed the 21 nt group and the 22 nt group, which account for 20.33% and 17.32% of the total unique small RNA sequences, respectively, to be two major groups (Figure 8B). Further examination of the chloroplast-derived small RNAs revealed that they were mapped to both sense and antisense orientation in the P. ginseng chloroplast genome (Figure 8C). There are two clusters of small RNAs abundant in the loci annotated as chloroplast rRNAs, which are located in the inverted repeats (IRa and IRb) of chloroplast genome [16]. These results suggest the diversity of chloroplast-derived small RNAs in P. ginseng. Chloroplast-derived small RNAs have been previously identified in various species, such as tobacco [111] and Cucumis [112]. In consistence with our results, two clusters of small RNAs derived from chloroplast rRNAs were also observed in melon, suggesting chloroplast small RNAs may be highly conserved.

Figure 8

Analysis of chloroplast-derived small RNAs.

(A) Size distribution of clean reads. (B) Size distribution of unique sequences. (C) The number of small RNA sequences mapped onto P. ginseng chloroplast genome. X-axis shows nucleotide positions of P. ginseng chloroplast genome. Y-axis shows the number of sRNAs mapped at each position. Black represents sense small RNAs, Red represents antisense small RNAs. The number of antisense sRNAs is shown as negative numbers.

Discussion

The estimated genome of P. ginseng is about 3.12×103 million letters of genetic code. It is 21.5 and 7.43 times of the genome of Arabidopsis and rice, respectively, indicating that the P. ginseng genome could be much more complicated than many other plant genomes [113], [114]. Consistently, we found in this study that a large, diverse and complex small RNA population existed in P. ginseng. Moreover, the population of P. ginseng small RNAs identified is far from exhausted because 88.5% of the total 1,502,664 unique small RNA sequences were sequenced only once or twice. Analysis of the size distribution of P. ginseng small RNAs revealed that the 24 nt group, which mainly contains repeat-associated siRNAs (ra-siRNAs) and transposable element-derived siRNAs (TE-siRNAs), is the most diverse (Figure 1), indicating the existence of a large number of repetitive DNA and transposons in the P. ginseng genome [115]–[117]. Many repetitive elements are located in the heterochromatic regions including centromeres and telomeres. ra-siRNAs generated from heterochromatic regions play significant roles in developing or/and maintaining the repressed chromatin state through DNA methylation and heterochromatic histone modification [55], [118]–[120]. It suggests that the existence of a complex and diverse small RNA population is particularly important to maintain the growth and development of plants, such as P. ginseng, which have a large genome.

From the P. ginseng small RNA data set, we identified a total of 73 conserved and 28 non-conserved miRNAs, which belong to 33 and 9 miRNA families, respectively. However, only four precursors of conserved miRNAs were identified from total 31,371 unigenes, suggesting the EST sequences available are far from sufficient for the identification of miRNA precursors. On the other hand, we were able to identify 25 non-conserved miRNA precursors from the same unigene set. It indicates that the number of non-conserved miRNAs is probably much larger than that of conserved miRNAs in P. ginseng. Among four precursors for conserved miRNAs, two (MIR482a and MIR2118) cluster in a transcript (Figure 3). Both miR482 and miR2118 were predicted to regulate disease resistance proteins (Table S4). Consistently, the expression patterns of miR482 and miR2118 are similar in the tissues analyzed. It suggests that miR482 and miR2118 are involved in defense responses in P. ginseng through relevant regulatory networks.

Among the non-conserved miRNAs, miR6136a.1 and miR6136a.2 are generated from the same precursor FW1NBNE01BEGNB, whereas its antisense transcript produces miR6136b. miR6136a.1 showed the highest expression in flowers, followed by leaves, whereas miR6136a.2 and miR6136b were mainly expressed in stems (Figure 5). The expression of miR6136b was altered in embryogenic calli treated with dehydration and heat, whereas both miR6136a.1 and miR6136a.2 appeared to be not responsive to the stresses (Figures 6 and 7). These results suggest the complexity of miRNA regulatory networks in P. ginseng.

Environmental factors are critical to plant growth and development. Increasing evidence indicates that miRNAs are involved in plant response to various environmental stresses, such as drought, salinity, cold, nutrient starvation, oxidative stress, mechanical stress and UV-B radiation [28], [30], [35], [100], [121]–[132]. In this work, we examined the expression profiles of non-conserved miRNAs under dehydration and heat treatments and identified 5 dehydration- and 10 heat-responsive P. ginseng miRNAs, which appeared to be an incomplete set of dehydration and heat-responsive miRNAs in P. ginseng, considering that the small RNAs analyzed in this study were isolated from plants growing under normal conditions. More stress-responsive miRNAs could be identified if plant tissues treated with dehydration and heat were used. Additionally, we found that a crosstalk could exist among some miRNAs in ginseng response to dehydration and heat stresses. However, the underlying mechanism remains to be elucidated.

Materials and Methods

Plant Materials

P. ginseng C.A. Meyer cv. Shizhu was grown in a field nursery at Kuandian, Liaoning province, China. Leaves, stems, flowers and roots were collected from five-year-old plants and frozen in liquid nitrogen for total RNA extraction and small RNA library construction. No specific permits were required for the described field studies. The location of the field nursery is not privately-owned or protected in any way and the field studies did not involve endangered or protected species. Embryogenic callus was induced as described [133]. Briefly, ginseng seeds were harvested and then stratified for about 6 months. Pericarp was removed before sterilization in 70% ethanol for 1 min and 1% sodium hypochlorite for 20 min. After washing with sterile water, cotyledon explants were dissected out from seeds and cultured on MS media containing 1.0 mg.L-1 2,4-D and 3% sucrose in the dark for one month. Explants with embryogenic calli were then cultured at 24±2°C under a 16/8-h (light/dark) photoperiod with light supplied by white fluorescent tubes at an intensity of 24 µmol.m-2.s-1. Embryogenic calli were maintained by regular 3-week subcultures.

Stress Treatments

For dehydration treatment, embryogenic calli were transferred to a 9 cm Petri dish with a dry sterile filter paper. The Petri dish was sealed with Parafilm and kept in a clean bench for 0, 1, 3, 6, 12 and 24 h at 24°C. For heat treatment, the calli on a callus induction medium were incubated at 37°C in the dark for 0, 1, 3, 6, 12 and 24 h. After each treatment, the calli were collected and blotted briefly with a piece of dry filter paper. Tissues were immediately stored in liquid nitrogen until use.

Small RNA Library Construction and High-throughput Sequencing

Total RNA was extracted from leaves, stems, flowers and roots of P. ginseng using the TRIzol Reagent (Invitrogen, USA). Equal amounts (5 µg) of total RNA from each tissue were pooled. The quality and quantity of total RNA were analyzed using Agilent 2100 (Figure S2). Small RNA library was prepared as described [134]. Briefly, 10–30 nt small RNAs were purified from a 15% denaturing polyacrylamide gel and then ligated with the 5′ and 3′ adapters. After being reverse-transcribed by Superscript II reverse transcriptase (Invitrogen, USA), small RNAs were amplified by PCR. High-throughput sequencing of small RNAs was performed using an Illumina G1 Genome Analyzer at BGI-Shenzhen, China.

Electrophoretic Analysis of Small RNAs

Low molecular weight (LMW) RNA from 200 µg pooled total RNA was resolved on a 15% denaturing polyacrylamide gel and stained with ethidium bromide as described by Lu et al. [49]. The pictures were taken with the Bio-Rad Universal Hood II Gel Imager.

Data Sets

The 20621 PlantGDB-assembled unique transcripts (PUTs) of P. ginseng were downloaded from PlantGDB (http://www.plantgdb.org, April 25th, 2011) and further assembled with singletons identified previously by EGassembler using the default parameters (http://egassembler.hgc.jp/) [13], [135]. The resulting 31371 unigenes were then used for small RNA mapping. Known mature miRNAs and miRNA precursors were downloaded from miRBase for the analysis of conserved miRNAs in P. ginseng [38].

Computational Identification of Small RNAs

Adapters were removed and low-quantity sequences were filtered from raw small RNA sequence data using the PHRED and CROSS MATCH programs [136]. The 18–30 nt clean small RNAs were mapped to known miRNA precursors using the SOAP2 program with no more than 2 mismatches allowed [40]. Small RNAs with matched precursors were then checked for mature miRNAs and miRNA* manually. For identification of non-conserved miRNAs, two approaches were used. First, small RNAs were mapped to 31371 unigenes of P. ginseng using the SOAP2 program with no mismatches allowed [40]. Unigenes with small RNA mapped were analyzed for secondary structures using the mfold program [75]. Second, small RNAs were mapped to P. ginseng unigenes and miRNAs were predicted using the MIREAP software (http://sourceforge.net/projects/mireap/) with the following parameters applied. It includes minimal miRNA length (18), maximal miRNA length (26), minimal miRNA (reference) length (20), maximal miRNA reference length (24), uniquencess of miRNA (280), maximal energy (-30), minimal space (5), maximal space (450), minimal mature pair (16), maximal mature bulge (3), maximal duplex asymmetry (4), and flank sequence length (20). The resulting hairpin structures from both approaches were then manually checked. The criteria applied to manual miRNA annotation includes (1) at least 2 reads for miRNAs, (2) a single-stranded, hairpin precursor with the free energy (dG) value less than -30, (3) four or fewer mismatched miRNA bases suggested by Meyers et al (2008) [137], (4) no more than two asymmetric bulges within the miRNA/miRNA* duplex, and (5) at least 0.85 of the minimal folding free energy index (MFEI) calculated as described [138].

Computational Prediction, Annotation and Experimental Validation of miRNA Targets

miRNA targets were predicted from the assembled P. ginseng unigene set by the psRNATarget server using default parameters [95]. Unigenes were annotated by blastx analysis against the Nr protein database using default parameters (http://www.ncbi.nlm.nih.gov/). The transcriptional direction of unigenes was determined based on the direction of predicted proteins. Unigenes without a hit in the Nr protein database were not considered as miRNA targets because the transcriptional direction was unknown without further experiments. Experimental validation of miRNA targets were carried out on mRNAs isolated from leaves or pooled tissues of leaves, stems, flowers and roots of P. ginseng using the modified 5′ RNA ligase-mediated (RLM)-RACE method as described previously [35]. The nesting and nested primers used for PCR were listed in Table S6.

Quantitative Real-time Reverse Transcription-PCR (qRT-PCR) Analysis of miRNAs

Total RNA was isolated from leaves, stems, flowers, roots and calli using the TRIzol Reagent (Invitrogen, USA) and digested with RNase-free DNase I (Promega, USA) to remove genomic DNA. Two qRT-PCR methods, including the poly(T) adaptor RT-PCR method [101] and the stem-loop RT-PCR method [108], [109], were used for miRNA quantification. The poly(T) adaptor RT-PCR method used about 1 µg DNase-treated total RNA as the starting material, which was polyadenylated using the Poly(A) Tailing kit (Ambion, USA) as described previously [101] and then reverse-transcribed into single-strand cDNA in a 20 µl reaction mix with SuperScript III Reverse Transcriptase (Invitrogen, USA) and 0.5 µg Poly(T) adapter (5′-GCGAGCACAGAATTAATACGACTCACTATAGG(T)12VN-3′, Ambion, USA). qRT-PCRs were performed in triplicates using the Bio-Rad CFX96 Real-Time PCR System C1000 Thermal Cycler (Bio-Rad, USA). Each PCR reaction was carried out in a final volume of 20 µl containing 10 µl 2×SYBR Premix Ex Taq II (TaKaRa, Japan), 0.25 µM each of the miRNA specific forward primer and the universal reverse primer (Table S7), and the cDNA reverse-transcribed from 100 pg total RNA using the following conditions: 95°C for 20 sec, 40 cycles of 95°C for 15 sec, 60°C for 15 sec and 72°C for 15 sec. 5.8 S rRNA was used as an endogenous reference as described previously [28], [35], [100], [101]. Standard deviations were calculated from three PCR replicates. The specificity of PCR products was verified based on single peak of dissociation curves. The relative abundance of miRNAs was determined using the 2-ΔΔCq method, where Cq represents the threshold cycle [139].

For the stem-loop RT-PCR method, about 100 ng DNA-free total RNA was hybridized with a miRNA-specific stem-loop RT primer. The hybridized miRNA molecules were then reverse-transcribed into cDNA as described [108], [109]. cDNA was diluted and then PCR-amplified in triplicates using the Bio-Rad CFX96 Real-Time PCR System C1000 Thermal Cycler (Bio-Rad, USA). Each PCR reaction was performed in a final volume of 20 µl containing 10 µl 2×SYBR Premix Ex Taq II (TaKaRa, Japan), 0.25 µM each of the miRNA specific forward primer and the universal reverse primer (Table S8), and the cDNA reverse-transcribed from 10 ng total RNA using the following conditions: 95°C for 20 sec, 40 cycles of 95°C for 10 sec, 60°C for 10 sec and 72°C for 6 sec. The endogenous reference, 5.8 S rRNA, was analyzed on the cDNA template converted from total RNA using random primers (Invitrogen, USA). In order to verify the specificity of PCR-amplification, some amplicons were cloned and sequenced.

Predicted hairpin structures of miRNA precursors. Mature miRNA sequences are indicated in red. miRNA* sequences are indicated in blue.

(DOC)

Click here for additional data file.

The quality and quantity of total RNA used for small RNA library construction. One µl of total RNA were diluted with 1 µl of RNase-free water and then denatured at 70°C for two min. Total RNA was analyzed using Agilent 2100. The results showed that the sample had a RNA integrity number (RIN) of 8.1 and a 28s/18s rRNA ratio of 1.8, suggesting the RNA samples were not degraded. The concentration of original total RNA is 302 ng/µl, which is good for small RNA library construction.

(DOC)

Click here for additional data file.

Statistics of high-throughput sequencing results of small RNAs.

(DOC)

Click here for additional data file.

Conserved miRNAs in .

(DOC)

Click here for additional data file.

miRNAs* putatively produced from conserved miRNA gene families in .

(DOC)

Click here for additional data file.

Targets of conserved miRNAs in .

(DOC)

Click here for additional data file.

Sequencing analysis of miRNA amplicons from stem-loop RT-PCR.

(DOC)

Click here for additional data file.

Primers used for 5′ RLM-RACE.

(DOC)

Click here for additional data file.

Primers used for poly(T) adaptor RT-PCR.

(DOC)

Click here for additional data file.

Primers used for stem-loop RT-PCR.

(DOC)

Click here for additional data file.