Deep sequencing identifies novel and conserved microRNAs in peanuts (Arachis hypogaeaL.)

Peanut has a complex small RNA population

To date, 92,988 peanut ESTs, including 86,724 ESTs from cultivated peanuts and 6,264 ESTs from wild-type peanuts, have been deposited in the NCBI EST database. These sequences are minor compared with the 2,800-Mb genome of the allotetraploid cultivated peanut or even the genome of the diploid wild-type peanut. Previous studies have demonstrated, using computational approaches and EST analysis, that only three conserved miRNAs exist in peanut [34, 38, 41]. With the limited amount of peanut ESTs in the EST database, it is not possible to perform a comprehensive study of peanut miRNAs using only a computational analysis. Experimental cloning and subsequent functional analysis, combined with computational prediction, appears to be the most effective method to identify peanut miRNAs.

Next generation high through-put sequencing, including 454 and Solexa technologies, provides a powerful tool for miRNA cloning. By using the high through-put Solexa sequencing technology, a total of 6,009,541 sequences were obtained from a small RNA library, which was constructed from the cultivated peanut variety Fenghua-1. After removing the low quality sequences and adapter sequences, 4,994,631 sequences were obtained with 3-30 nt in length, among which 4,598,005 sequences ranged from 18-30 nt in length. After further removing tRNAs (437,117), rRNAs (133,410), snRNAs (1,282), and snoRNAs (240), a total of 4,025,956 small RNA sequences were obtained. Although some small RNAs were very high in abundance and present thousands of times in our dataset, the majority of small RNAs were sequenced only a few times. For example, 2,232,910 out of 4,598,005 small RNAs were sequenced only one time in our experiment. This result suggests that (1) the expression of different small RNAs in peanut varies drastically and (2) the small RNA survey in peanut is far from saturated. This also suggests that peanut contains a large and diverse small RNA population.

In peanut, the size of the small RNAs was not evenly distributed (Figure 1). Among these sequences, the number of 24-nt sequences was significantly greater than shorter or longer sequences (Figure 1) and accounted for 45% of the total sequence number. This result was consistent with that of Medigcago [42] and rice [43], as well as Arabidopsis 454 sequencing results [44]. In Arabidopsis, the 24-nt small RNAs accounted for about 60% of its small RNA transcriptome [45]. However, the size distribution differs from wheat and conifer sequences obtained through 454 high through-put sequencing [43, 46] and Chinese yew sequences obtained through Solexa sequencing [47]. In conifer, the fraction of 24-nt RNAs was very small (2.6%) due to the lack of DCL3, the enzyme that matures 24-nt RNAs in angiosperms [43, 48]. In total, 620,060 sequences (13.5%) with 21-nt, which is the typical length of plant mature miRNAs, represented the second highest abundance of sequences in the peanut library.

Identification of conserved peanut miRNAs

To identify conserved miRNAs in peanuts, all small RNA sequences were Blastn searched against the currently known miRNAs in the miRNA database miRBase (March 9, 2009). In total, 1,763 known miRNAs from diverse plant species were utilized in order to identify conserved peanut miRNAs from the small RNA dataset.

Next generation high through-put sequencing provides an alternative way to estimate expression profiles of protein coding genes and/or miRNA genes [44, 46]. Millions of peanut small RNA sequences, generated from Solexa sequencing, allowed us to determine the abundance of various miRNA families and even distinguish between different members of a given family. Interestingly, peanut miRNA families displayed significantly varied abundance from each other. For example, ahy-miR157a, ahy-miR168a, and ahy-miR156a were detected 95,381, 19,898, and 17,058 times respectively (Table 1). In comparison to other plant species, tae-miR169b in wheat and osa-miR169 in rice were the most frequently sequenced miRNAs while miR156 in rice and wheat exhibited low abundance [46]. This suggests a species-specific expression profile for miRNAs. miR156a was also found to be highly expressed in another legume species, Medicago [49]. In Arabidopsis, miR156a, located on chromosome 2 [49], targets 10 mRNAs that code for the squamosa promoter-binding protein (SBP) box, which is involved in leaf morphogenesis [50]. Similar to miR156a, miR157a, which is located on chromosome 1 in Arabidopsis thaliana, was thought to target mRNAs coding for proteins comprising the SBP box [49]. However, the mechanisms, causing the differential expression profile of a same miRNA in different plant species, are unknown. A majority of peanut miRNAs were only sequenced less than 1,000 times, and some rare miRNAs were detected less than 10 times. Compare with the most abundant miRNA ahy-miR157a, their expression level is about 9,500 times lower (Table 1). miRNAs of moderate abundance included ahy-miR157d, ahy-miR164a, ahy-miR166a, ahy-miR166 g, ahy-miR166a, ahy-miR167f, and ahy-miR172a were detected 2,000-10,000 times in the library. The relative abundance of the 22 conserved peanut miRNA families is represented in Figure 2.

Next generation high through-put sequencing technology also provides a method for distinguishing and measuring miRNA sequences with only a few nucleotide changes. Based on the results from the Solexa sequencing, different family members displayed drastically different expression levels. For example, the abundance of miR156 family varied from 4 read (ahy-miR156f) to 17,058 reads (ahy-miR156a) in the deep sequencing. This was also the case for some other miRNA families, such as ahy-miR164 (from 1 read to 4,116 reads) and ahy-miR166 (from 1 read to 9577 reads). The existence of a dominant member in a miRNA family may suggest that the regulatory role of this family was performed by the dominant member at the developmental time when the samples were collected for RNA extraction. Abundance comparisons of different members in one miRNA family, during various growth conditions or specific developmental stages, may provide valuable information on the role that miRNAs play in plant growth. Expression levels of two members of the ahy-miR159 family (ahy-miR159a and ahy-miR159b) were similar and were detected 66 and 41 times, respectively (Table 1).

Identification of novel peanut miRNAs

Besides these 14 identified novel candidate miRNAs, we also discovered two small RNAs, with 701 and 159 reads in our small RNA dataset, which correspond to Phaseolus vugaris legume-specific miRS1 and miR2118. These two miRNAs were able to detected in peanut by northern blot analysis [51]. Interestingly, the expression of miR2118 has previously been shown to be induced in Phaseolus vugaris by abiotic stress, especially drought and ABA treatment [51]. We did not include these two sequences in the list of novel peanut miRNAs because we could not find their precursor sequences in the current databases. In addition to miRS1 and miR2118, we also found the third small RNA with 137 reads in our dataset that had only one mismatch with Phaseolus vugaris miR159.2. A fourth 21-nt small RNA with 729 reads was also identified in our dataset, which had 4 mismatches and one nucleotide missing to compare with Phaseolus vugaris miR482*.

Based on the number of detection times and sequences in the small RNA library, novel peanut miRNAs displayed lower expression levels compared to the majority of conserved families. The low abundance of novel miRNAs might suggest a specific role for these miRNAs under various growth conditions, in specific tissues, or during developmental stages. The library enriched only small RNAs that play a role during early seedling stages under normal growth conditions. Whether these low-abundant miRNAs are expressed at higher levels in other tissues and organs, such as flowers, gynophores, pods, or seeds, or whether they are regulated by biotic or abiotic stress, remains to be investigated. Future experiments would provide more insight into the function of these miRNAs.

Validation of peanut miRNAs

Stem-loop qRT-PCR is a reliable method for detecting and measuring the expression levels of miRNAs. The stem-loop primers increase the sensitivity of the reactions such that this method can significantly distinguish two miRNAs with only one single nucleotide change [52]. In this study, we adopted this technique to validate and measure the expression of 4 novel miRNAs (miRn1, miRn2 and miRn2*, miRn3, and miRn4) as well as 5 conserved miRNAs (miR156, miR157, miR162, miR172, and miR396). All of these miRNAs were identified in peanut by Solexa sequencing. The qRT-PCR results demonstrate that all tested miRNAs, and one miRNA*, are expressed in peanut leaves (Figure 3). However, the expression levels of the different miRNAs varied.

The results of the qRT-PCR reaction show that conserved miRNAs are expressed in peanut. Based on the threshold cycle (C_T), miR172 and miR156 were highly expressed with C_T values of 19.6 ± 3.5 and 20.5 ± 5.3, respectively. In one of our previous studies, we also found that miR172 is highly expressed in cotton leaves [53]. Other studies have shown that conserved miR172 and miR156 play very important roles in plant growth and development [41]. miR156 is involved in Arabidopsis leaf development by negatively regulating the Squamosa-promoter binding protein (SBP) [38, 42]. miR172 controls flower development by regulating the expression of apetala2 (ap2) in Arabidopsis [4, 43] and glossy 15 in maize [44]. Aberrant expression of miR156 and miR172 in plants disrupts normal leaf and flower development. Compared with miR156 and miR172, the expression levels of miR157 and miR162 are moderate while the expression of miR396 is low. The expression patterns of these miRNAs appear to be related to their function.

Four novel miRNAs and one miRNA*, all identified by Solexa sequencing, were validated by qRT-PCR. The expression levels of the miRNAs differed from one another in peanut leaves. miRn2 and miRn1 were expressed much higher than other tested peanut-specific miRNAs with a C_T value of 21.2 ± 1.0 and 24.6 ± 3.2, respectively. The expression levels are much lower for miRn3 and miRn2* with C_T values of 37.9 ± 1.8 and 33.1 ± 4.2. However, more studies need to be performed to elucidate the function that these miRNAS have on the growth and development of peanut.

Target prediction of peanut miRNAs

To better understand the functions of the newly identified species-specific as well as conserved peanut miRNAs, putative targets of these miRNAs were predicted using the described criteria and methods. The target genes of thirteen conserved and seven novel peanut miRNA families were predicted. Transcription factors, including GRAS family transcription factor, nuclear transcription factor Y subunit and NAC1 were predicted to be potential targets of peanut miRNAs. Furthermore, genes directly involved in protein synthesis, e.g., ribosomal protein S12, were targets of peanut miRNAs. A previous study indicates that auxin signaling is regulated by miRNAs [18]; our current result is consistent with this study and the auxin signaling F-box 3 is a potential target of peanut miR393. Resveratrol synthase, NAM (no apical meristem)-like protein, growth regulator factor 5, basic blue copper protein, endonuclease, a protein kinase, transport inhibitor response 1 and a disease resistance response protein were also predicted to be potential targets of identified peanut miRNAs (Additional file 2).

Plant materials

Peanuts (Arachis hypogaea L. cv. Fenghua-1) were grown in a growth chamber, with a light intensity of 3000 lx, at a relative humidity of 75%, and 26/20°C day/night temperatures. Leaves, stems, and roots from 14-day-old seedlings were collected and immediately stored in liquid nitrogen until total RNA extraction.

RNA extraction and miRNA cloning

Total RNA was isolated from leaves and roots using Trizol agent (TaKaRa, Dalian, China), according to the manufacturer's instructions. Total RNA was isolated from stems using a modified CTAB method with isopropanol instead of lithium chloride for RNA precipitation [54]. Briefly, one gram of stem tissue was ground to a fine powder using liquid nitrogen and mixed thoroughly with 5 ml of pre-warmed (65°C) extraction buffer (2% CTAB, 2% PVP, 0.1 M Tris-HCl, 2.0 M NaCl, 25 mM EDTA, 2% beta-mercaptoethanol, pH 8.0). The mixture was incubated at 65°C for 5 min and shaken three individual times during the incubation period. After a brief cooling of the mixture, 2.5 ml of chloroform and 2.5 ml of isopropanol were added. The mixture was vortexed for 1 min and then centrifuged at 12000 rpm for 15 min at 4°C. After DNase treatment of the extract, RNA was precipitated at room temperature (25°C) for 10 min using an equal volume of isopropanol. The RNA was resuspended in an equal volume of phenol:chloroform:isopropanol (25:24:1), and then resuspended again with an equal volume of chloroform:isopropanol (24:1). A total of 1/10 volume of 3 M NaOAC (pH 5.2) and 2.5 volumes of cold ethanol were added to precipitate the RNA overnight at -20°C.

To identify as many tissue- or developmental-specific miRNAs as possible, we pooled the total RNAs from leaf, stem, and root samples in an equal fraction ratio. miRNA cloning was performed as described previously by Sunkar and Zhu [21]. Briefly, 0.5 M NaCl and 10% PEG8000 were used to precipitate and enrich RNAs with low molecular weight. Next, a 15% polyacrylamide denaturing gel was employed to separate the low-molecular weight RNA. During gel electrophoresis, about 100 μg of total RNA was applied to the gel and two labeled RNA oligonucleotides, approximately 18 and 26 nt in length, were used as size standards. After gel electrophoresis, small RNAs with 18-26 nt were excised from the gel and eluted with 0.4 M NaCl overnight at 4°C. The RNA was dephosphorylated using alkaline phosphatase (New England Biolabs, Beijing China) and recovered by ethanol precipitation. The isolated small RNAs were then sequentially ligated to RNA/DNA chimeric oligonucleotide adapters, reversely transcribed, and amplified by PCR. Finally, Solexa sequencing technology was employed to sequence the small RNAs from pooled peanut samples (BGI, Beijing China).

Identification of conserved and peanut-specific miRNAs

The raw sequences were processed using PHRED and CROSS MATCH programs as previously reported [21, 55]. After removing the vector sequences, trimmed sequences longer than 17 nt were used for further analyses. First, rRNA, tRNA, snRNA, and snoRNA, as well as those containing the polyA tail, were removed from the small RNA sequences and the remaining sequences were compared against rice and Arabidopsis ncRNAs deposited in the NCBI Genbank database and Rfam8.0 database. Then, the unique small RNA sequences were used to do a Blastn search against the miRNA database, miRBase 13.0, in order to identify conserved miRNAs in peanuts. Only perfectly matched sequences were considered to be conserved miRNAs. To study potential miRNA precursor sequences, we used the identified peanut mature miRNA sequences to do Blastn searches against peanut ESTs in NCBI. Non-coding ESTs, which met previously described criteria [56], were then considered to be miRNA precursors. Specifically, dominant, mature sequences residing in the stem region of the stem-loop structure and ranging between 20-22 nt with a maximum free-folding energy of -25 kcal mol^-1 were considered. A maximum of six unpaired nucleotides between the miRNA and miRNA* was allowed. The distance between the miRNA and miRNA* ranged between 5 and 240-nt. After removing the conserved miRNA sequences, the rest of the small RNA sequences were used to perform Blastn searches against peanut ESTs in order to obtain precursor sequences for novel potential miRNAs. The selected EST sequences were then folded into a secondary structure using an RNA-folding program mFold. If a perfect stem-loop structure was formed, the small RNA sequence was sit at one arm of the stem as well as other criteria were followed, this small RNA was consisted as one novel peanut miRNA.

miRNA validation

The RT-PCR and qRT-PCR reactions, for validating and detecting peanut miRNAs, were followed using the same protocols as our previous report [37, 53]. Briefly, miRNA reverse transcription reactions were performed in 200 μL PCR tubes, each containing a total of 20 μL of reaction solution. Each reaction solution contained 1000 ng of total leaf RNAs, 3.33 U/μL MultiScribe reverse transcriptase, 1× reverse transcription buffer, 0.25 mM each of dNTPs, and 0.25 U/μL RNase inhibitor; sterilized RNase-free water was used to adjust the total volume of the reverse transcription reaction to 20 μL. The miRNA reverse transcription reactions were incubated in an Eppendorf Mastercycler (Eppendorf North America, Westbury, NY). The RT-PCR temperature program was adjusted to run for 30 min at 16°C, 30 min at 42°C, 5 min at 85°C, and then 4°C until future use. For each miRNA, three biological replicates were performed. After reverse transcription, the products of each reaction were diluted 10 times to avoid potential primer interference in the following qRT-PCR reaction.

Quantitative real time PCR was performed using the TaqMan^® microRNA Assay kit (Foster City, CA) on an Applied Biosystems 7300 Sequence Detection System (Foster City, CA). Each reaction consisted of 3 μL of product from the diluted reverse transcription reaction, 2 μL of 20× TaqMan MicroRNA Assay primers (forward and reverse), 12.5 μL of 2× TaqMan Universal PCR Master Mix, and 7.5 μL of nuclease-free water. The reactions were incubated in a 96-well plate at 95°C for 10 min, followed by 45 cycles of 95°C for 15s and 60°C for 60s. After the reactions were completed, the threshold was manually set and the threshold cycle (C_T) was automatically recorded. The C_T is defined as the fractional cycle number at which the fluorescence signal passes the fixed threshold [52]. All reactions were run in two replicates for each sample.

Target gene prediction

The potential targets of peanut miRNAs were predicted using the psRNATarget program http://bioinfo3.noble.org/psRNATarget/ with default parameters. Newly identified peanut miRNA sequences were used as custom miRNA sequences; Arachis transcript/genomic library (EST, GSS, and nucleotide databases) were used as custom plant databases.

All predicted target genes were evaluated by scoring, and the criteria used were as follows: each G:U wobble pairing was assigned 0.5 points, each indel was assigned 2.0 points, and all other non-canonical Watson-Crick pairings were assigned 1.0 points each. The total score for an alignment was calculated based on 20 nt. When the query was longer than 20 nt, scores for all possible consecutive 20 nt subsequences were computed, and the minimum score was considered the total score for the query-subject alignment. Because targets complementary to the miRNA 5' end appear to be critical, mismatches other than G:U wobbles at positions 2-7 at the 5' end were further penalized by 0.5 points in the final score [57]. Sequences were considered to be miRNA targets if the total score was less than 3.0 points.

Once potential target mRNA sequences were obtained, redundant sequences were removed using the 'contig express' feature of the Vector NTI program. Blastx was performed using the target sequence and the NCBI database to predict functions of potential targets.