Abstract

To gain insights into a possible role of microRNAs in myeloproliferative neoplasms, we performed short RNA massive sequencing and extensive bioinformatic analysis in the JAK2V617F-mutated SET2 cell line. Overall, 652 known mature miRNAs were detected, of which 21 were highly expressed, thus being responsible of most of miRNA-mediated gene repression. microRNA putative targets were enriched in specific signaling pathways, providing information about cell activities under massive posttranscriptional regulation. The majority of miRNAs were mixtures of sequence variants, called isomiRs, mainly because of alternative, noncanonical processing of hairpin precursors. We also identified 78 novel miRNAs (miRNA*) derived from known hairpin precursors. Both major and minor (*) forms of miRNAs were expressed concurrently from half of expressed hairpins, highlighting the relevance of miRNA* and the complexity of strand selection bias regulation. Finally, we discovered that SET2 cells express a number of miRNA-offset RNAs (moRNAs), short RNAs derived from genomic regions flanking mature miRNAs. We provide novel data about the possible origin of moRNAs, although their functional role remains to be elucidated. Overall, this study shed light on the complexity of microRNA-mediated gene regulation in SET2 cells and represents the basis for future studies in JAK2V617F-mutated cellular models.

Introduction

The JAK2V617F mutation, which occurs in most patients with polycythemia vera and more than or equal to 60% of essential thrombocythemia and primary myelofibrosis,1  is considered integral to the pathogenesis of myeloproliferative neoplasms (MPNs), although additional, antecedent mutations are required for an MPN to develop.2-4  The disease can be reproduced in mice expressing the JAK2V617F allele.5  Cells harboring the JAK2V617F mutation display autonomous activation of several cell signaling pathways, particularly JAK/STAT, and proliferate and mature in a cytokine-independent manner.6  Recent information highlighted that deranged epigenetic gene regulation7,8  and abnormal expression of microRNAs9-11  also contribute to the pathogenesis of MPNs.

MicroRNAs (miRNAs) constitute a large class of noncoding RNAs that act as posttranscriptional regulators of gene expression. Because dysregulation of miRNA expression plays a critical role in genetic and multifactorial disorders and most, if not all, human cancers,12  miRNAs are considered prominent tumor markers, oncogenes (oncomiRNAs), relevant targets for therapy, or therapeutic agents depending on the specific context.13  We and others have identified abnormally expressed microRNAs in CD34+ and blood cells of MPN patients, some of which have a JAK2-dependent pattern of expression.9-11  One of these, miR-16, is overexpressed in polycythemia vera CD34+ cells and contributes to the expansion of erythropoiesis.14 

There is an intrinsic complexity of the microRNA system that derives both from the large number of coding RNAs that represent possible targets of a single mature microRNA and the variety of mature different miRNA molecules that originate from a single genetic locus. Both for intragenic and intergenic (single or clustered) miRNAs, primary miRNAs are transcribed, edited, and cleaved in the nucleus to generate hairpin precursors that are exported to the cytoplasm where maturation occurs.15  Two different mature miRNA sequences can derive from the same precursor hairpin: a major form (eg, hsa-let-7c), which is the most stable and prevalent, and a minor form (eg, hsa-let-7c*). The pre-miRNA hairpin is cleaved in the cytoplasm by Dicer to produce miRNA duplex (∼ 22 nt),16  which is incorporated in the RNA-induced silencing complex (RISC). The RISC recognizes the duplex, unwinds it, selects the guide miRNA strand (while degrading the passenger strand) and mediates recognition and silencing of target coding RNAs. Both miRNA forms can be expressed and are functionally effective, but sometimes the “minor” form prevails.17-19  Major and minor forms are associated with different sets of target coding mRNAs, contributing uniquely to the regulation of cell activities.

Furthermore, deep sequencing technologies (RNAseq) led to the discovery of a great number of novel miRNAs in the last years,20  including in normal and malignant blood cells.21,22  It has also been shown in animals (from Drosophila melanogaster to humans) and plants that mature miRNA sequences are not invariant; rather, they constitute a mixture of miRNA variants: isoforms heterogeneous in 5′- and 3′ ends or with nontemplate 3′ additions of nucleotides, which are collectively called isomiRs.23,24  These variants, all potentially active, target only partially overlapping sets of coding mRNAs. Finally, a novel class of miRNA-related RNAs, micro-RNA offset RNAs (moRNAs), has been recently identified by RNAseq.25 

moRNAs were first reported in the simple chordate ascidian Ciona intestinalis26  as approximately 20-nt-long RNAs that originate predominantly from the 5′-arm of the pre-miRNA, regardless of the position of the major miRNA, suggesting that moRNA and miRNA biogenesis may be linked but not necessarily interdependent. The expression of moRNAs in C intestinalis is developmental regulated, and their abundance sometimes exceeds the corresponding mature miRNA. Initially, C intestinalis moRNAs were considered byproducts of an miRNA processing machinery with intrinsic properties.26  However, shortly thereafter, Langenberger et al reported moRNAs, derived from 78 miRNA loci, as being weakly expressed in human prefrontal cortex,27  whereas Taft et al showed that moRNAs were enriched in the nucleus in the human leukemia cell line THP-1 and almost exclusively derived from the 5′-arm, regardless of the position of the prevalent miRNA.28  Some moRNAs were even more expressed than the prevalent miRNA of the same locus. Thereafter, moRNAs (eg, hsa-miR-410 5′moR, has-miR-326 5′moR) were discovered in solid tumors.29 

In this study, we used massive sequencing of short RNAs with the aim to discover novel miRNAs in known hairpin precursors, characterize and quantify isomiRNAs, and identify moRNAs in SET2 cells, a JAK2V617F-mutated cell line. Furthermore, the expression relationships between sister miRNA pairs (miR/miR* of the same hairpin), between sister moRNAs and between moRNAs and the corresponding miRNAs were investigated. The effects of a JAK2 inhibitor on the modulation of miRNA expression in SET2 cells were also evaluated. This report shows an unforeseen complexity of small RNAs in this cellular system and represents the starting point for further functional studies and systematic analysis in primary MPN cells.

Methods

RNA isolation, short RNA library construction, and sequencing

Total RNA, prepared from exponentially growing SET2 cells with miRNeasy reagents (QIAGEN), was used to generate a small RNA library according to a standard Illumina procedure (supplemental Methods; see the Supplemental Materials link at the top of the article). The purified cDNA library was used for cluster generation on Illumina's Cluster Station and sequenced on Illumina GAIIx. Raw sequencing reads were obtained using Illumina's Pipeline Version 1.5 software after analysis of sequencing image by Pipeline Firecrest Module and base-calling by Pipeline Bustard Module. Library construction and sequencing were performed at LC Sciences (www.lcsciences.com)

Short RNA sequencing data processing

A detailed description of the computational pipeline used for data handling and analysis is reported in supplemental Methods, including a flow-chart outline of study procedures (supplemental Figure 1).

KEGG pathway enrichment

KEGG pathway enrichment for known miRNAs was obtained using DIANA-miRPath webtool30  based on TargetScan5 miRNA target predictions. mirPath is used to identify molecular pathways potentially altered by the expression of single or multiple microRNAs through enrichment analysis of multiple microRNA targets and comparison of each set of microRNA targets with all KEGG pathways. The combinatorial effect of coexpressed microRNAs is taken into account by simultaneous analysis of multiple microRNAs. Only pathways including at least 5 genes and displaying an enrichment P value of more than .001 were considered significant.

New miRNA target prediction

Gene targets of new miRNAs were predicted by TargetScan using as input new miRNA sequences and the set of 3′UTR sequences of all human transcripts and orthologs in 23 species (TargetScan database Version 5.2; for details, see supplemental Methods).

Analysis of miRNA expression in SET2 cells exposed to the JAK1/JAK2 inhibitor INC424

SET2 cells were incubated with increasing concentrations of INC424/ruxolitinib for 3 to 6 hours, and RNA was purified. Affymetrix Genechip miRNA array Version 2.0 was used to identify miRNA species showing a differential expression in treated cells compared with control cells (for a detailed description of procedure, see supplemental Methods).

Results

The Illumina GAIIx sequencing of the small RNA library from SET-2 cells produced 32 760 003 reads which, after extensive preprocessing and quality control, were reduced to 27 906 609 reads, 85% of sequenced reads. A first mapping phase aimed at discarding contaminations and repeats, yet tolerating for possibly unknown miRNA loci, produced 22 167 999 reads (68% of raw data). These were mapped to “extended hairpins” to identify and quantify known miRNAs and for discovery and characterization of novel isomiRs and other miRNA-associated expressed RNAs. In total, 1421 known hairpin precursors, corresponding to 1731 known mature miRNA sequences and including a number of miRNAs recently discovered by RNAseq, were finally considered (supplemental Table 1). A detailed description of sequencing data processing is provided in supplemental Methods. Raw sequencing data and processed files are available at GEO database (GSE30672).

Known miRNAs expressed in SET2 cells

A total of 652 known miRNAs were found expressed in SET2 cells, with expression levels ranging from 10 to 2 268 333 (mean, 29 830; median, 613); 300 and 124 miRNAs presented read counts of at least 103 and 104, respectively (supplemental Table 2). Twenty-one only highly expressed miRNAs accounted for 70% of known miRNA expression (Table 1) and are thus predicted to account for most of miRNA-mediated gene repression in SET2 cells. Supplemental Table 3 reports that KEGG pathways significantly enriched in genes predicted target of 21 most expressed miRNAs. These include pathways anticipated to have functional relevance for MPN-associated cellular abnormalities, such as the MAPK signaling pathway, TGF-β signaling pathway, mTOR signaling pathway, and Wnt signaling pathway. These findings produced meaningful insights about cell activities that are under the control of massively expressed miRNAs in SET2 cells. Supplemental Table 4 reports functional insights for the 21 most expressed miRNAs.

Table 1

Twenty-one miRNAs highly expressed account for 70% of total known miRNA expression in SET2 cells

miRNAExpression levelCumulative % of total known miRNA expression
hsa-miR-21 2268333 12 
hsa-miR-148a 1737364 21 
hsa-miR-146b-5p 1632959 29 
hsa-miR-101 1293074 36 
hsa-miR-142-3p 1037658 41 
hsa-miR-19b 1025647 46 
hsa-miR-378 561770 49 
hsa-miR-92a 538422 52 
hsa-miR-191 416870 54 
hsa-miR-425 391247 56 
hsa-miR-126 349719 58 
hsa-miR-17 314388 59 
hsa-miR-181b 310800 61 
hsa-miR-181a 301902 63 
hsa-miR-93 223220 64 
hsa-miR-126* 216687 65 
hsa-miR-30e 201602 66 
hsa-miR-99b 200438 67 
hsa-let-7f 196921 68 
hsa-miR-25 191693 69 
hsa-miR-20a 186436 70 
miRNAExpression levelCumulative % of total known miRNA expression
hsa-miR-21 2268333 12 
hsa-miR-148a 1737364 21 
hsa-miR-146b-5p 1632959 29 
hsa-miR-101 1293074 36 
hsa-miR-142-3p 1037658 41 
hsa-miR-19b 1025647 46 
hsa-miR-378 561770 49 
hsa-miR-92a 538422 52 
hsa-miR-191 416870 54 
hsa-miR-425 391247 56 
hsa-miR-126 349719 58 
hsa-miR-17 314388 59 
hsa-miR-181b 310800 61 
hsa-miR-181a 301902 63 
hsa-miR-93 223220 64 
hsa-miR-126* 216687 65 
hsa-miR-30e 201602 66 
hsa-miR-99b 200438 67 
hsa-let-7f 196921 68 
hsa-miR-25 191693 69 
hsa-miR-20a 186436 70 

miRNA sequence heterogeneity (isomiRs)

The majority of expressed mature miRNAs were not represented by a unique sequence corresponding to that annotated in miRBase. We evaluated the whole group of reads belonging to each miRNA, including the “classic” mature sequence annotated in miRBase (“exact” alignment) as well as those reads perfectly matching the precursor but overlapping the mature position by 3 nt (longer/shorter), those presenting 1 mismatch (1-Mismatch), and those presenting 2 mismatches at the 3′-end (2-3′-Mismatches; supplemental Methods).

Considering the whole set of variants presenting at least 10 reads each, 636 miRNAs were identified: 232 (36%) appeared “invariant,” whereas the remaining 404 (64%) represented a mixture of 2 to 6 sequence variants (supplemental Figure 2). To be conservative, we considered as “biologically meaningful” sequence variants (isomiRs) accounting for at least 10% of the total per miRNA read count, calculated over all variants (supplemental Figure 2). In this way, of 644 miRNAs represented, 224 miRNAs (35%) are “invariant” and 420 (65%) are associated with 2 to 6 isomiRs each. Nevertheless, it is interesting that some of miRNAs considered “invariant” are associated with an unique major sequence as well as to other sequences that have not been annotated because none of them satisfied the required number of reads.

We reasoned that the presence of isomiRs could be particularly relevant for miRNAs expressed at the highest level. In Figure 1A, we report the distribution of number of isomiRs per miRNA separately for the 25% most expressed miRNAs (MEm) and for the remaining weakly expressed miRNAs. MEms were associated with one to 4 isomiRs. The contribution of the major isomiR (ie, the isomiR with the highest read count among all isomiRs of the miRNA) to total count was highly variable, accounting for a fraction of total count, which was inversely related to the total number of isomiRs per miRNA (Figure 1B). Specifically, the contribution of the major isomiR averaged 52% in case of miRNAs with more than one isomiRs and 83% in “invariant” miRNAs. As known, 5′- and 3′-regions of mature miRNA sequence play a different role in target recognition: in case of canonical target sites, the seed region, whose pairing to the target is crucial, is included in the 5′ of mature miRNA. It is conceivable that isomiRs with different seed sequence recognize and regulate different targets. Thus, for each MEm, we considered the 5′- and 3′-half sequences separately and classified expressed isomiRs according to the observed difference between the isomiR sequence and the “classic” mature miRNA (mismatch or length difference) in the involved region. Of 819 isomiRs, 648 (79%) differed in the 3′-region, 64 in the 5′-region (8%), and 107 in both (13%). No “Mismatch” isomiRs were observed regarding the 5′-region. Supplemental Table 5 includes details about isomiR sequences and counts for the group of highly expressed miRNAs. Figure 2 reports the relative contribution of isomiRs for the 21 most expressed miRNAs in SET2 cells.

Figure 1

The majority of miRNAs expressed in SET2 cells are mixtures of different isomiRs. (A) Distribution of number of isomiRs per miRNA in the top 25% most expressed miRNAs and in the remaining more weakly expressed miRNAs. (B) Contribution of the major isomiR per miRNA to the total count, for classes of highly expressed miRNAs composed of 1 to 4 isomiRs.

Figure 1

The majority of miRNAs expressed in SET2 cells are mixtures of different isomiRs. (A) Distribution of number of isomiRs per miRNA in the top 25% most expressed miRNAs and in the remaining more weakly expressed miRNAs. (B) Contribution of the major isomiR per miRNA to the total count, for classes of highly expressed miRNAs composed of 1 to 4 isomiRs.

Figure 2

Relative isomiR contribution to the total read count of 21 most expressed miRNAs in SET2 cells. The total read count belonging to sequences not passing the stringent threshold for being considered as isomiRs is indicated in the white rectangle.

Figure 2

Relative isomiR contribution to the total read count of 21 most expressed miRNAs in SET2 cells. The total read count belonging to sequences not passing the stringent threshold for being considered as isomiRs is indicated in the white rectangle.

In the MEm set, 323 isomiRs were identified: 140 (43%) were coincident with the miRNA sequence reported in miRBase (classic isomiR), 154 (48%) represented sequence variants with longer or shorter 5′- and 3′-ends, and 29 (9%) showed 1-nt difference versus the classic isomiR. Interestingly, we found that the major isomiR was not coincident with the classic sequence annotated in miRBase in 53 of 161 highly expressed miRNAs (33%). Besides, the isomiR corresponding to the sequence annotated in miRBase was not detected in 21 (13%) of the miRNAs expressed in SET2 cells.

We then considered 113 highly expressed miRNAs associated with more than one isomiR to understand how much mismatch alignments contributed to sequence variability. We found that 87 of these (77%) did not include isomiRs with single nucleotide substitutions respective to the sequence annotated in miRBase, whereas only a minority, which may arise from SNPs, RNA editing, and sequencing errors, included one or 2 isomiRs aligning with one mismatch to the hairpin precursor (24 and 2 miRNAs, respectively). In terms of total expression, isomiRs with mismatch alignments contributed only to less than 3% of read count. Classic isomiRs and isomiRs showing 5′- and 3′-length variability almost equally contribute to 97% of total read counts of highly expressed miRNAs. Therefore, we concluded that the most important sequence variation is represented by 5′- and 3′-length variability, possibly occurring as a consequence of alternative, noncanonical, regulated processing of the precursor sequence.

Novel miRNAs expressed in SET2 cells were discovered in known hairpin precursors

For novel miRNA discovery, we operationally defined as “expressed RNA elements” those discrete hairpin regions covered by overlapping reads with a minimum count of 10 and a start position within 4 nt each from the following one. We found that a discrete number of regions located outside known mature miRNAs were expressed from detectable to high level; a number of expressed RNA elements able to pair with known miRNAs in the most probable duplex produced by Dicer processing of the hairpin structure were identified. Specifically, we considered 943 hairpins associated with only one annotated mature miRNA, whereas an additional 478 included 2 known sister mature miRNAs.

The analysis of known hairpin precursors associated with only one known miRNA produced a set of 78 novel miRNAs expressed in SET2 cells (supplemental Table 6). Plots in Figure 3 show, for the 4 most expressed new miRNAs (hsa-miR-1307*, hsa-miR-376a-2*, hsa-miR-382*, and hsa-miR-539*), the number of reads aligned per nucleotide position. This information was integrated with hairpin sequence and folding data to define new miRNA sequences. Table 2 reports name, sequence, and expression level of novel miRNAs presenting a read count of more than 500. Expression levels of new miRNAs ranged from 10 to 72 670 (mean, 1471; median, 63) and were significantly lower than known miRNAs (2-sample t test of mean equality, P = 6.521 × 10−06; supplemental Figure 3). Nevertheless, 11 new miRNAs (14%) showed an expression level higher than the median value observed for known miRNAs. In particular, hsa-miR-1307*, hsa-miR-376a-2*, and hsa-miR-382* resulted very highly expressed, at a level even more than 75% of known miRNAs.

Figure 3

Examples of new miRNAs discovered. Plots show the number of reads per nucleotide position mapping in extended hairpin loci for hsa-mir-1307, hsa-mir-376a-2, hsa-mir-382, and hsa-mir-539, expressing both known and new miRNAs. Known miRNA position respective to the extended hairpin is indicated in the plot. In the top part of panel, the sequence of mature miRNAs is shown (red indicates known miRNA; and pink, new miRNA).

Figure 3

Examples of new miRNAs discovered. Plots show the number of reads per nucleotide position mapping in extended hairpin loci for hsa-mir-1307, hsa-mir-376a-2, hsa-mir-382, and hsa-mir-539, expressing both known and new miRNAs. Known miRNA position respective to the extended hairpin is indicated in the plot. In the top part of panel, the sequence of mature miRNAs is shown (red indicates known miRNA; and pink, new miRNA).

Table 2

Name and predicted sequence of 15 new miRNAs with read count > 500

New miRNASequenceExpression
hsa-miR-1307* CTCGACCGGACCTCGACCGGCTCGT 72 670 
hsa-miR-376a-2* GGTAGATTTTCCTTCTATGGTTA 14 296 
hsa-miR-382* CGAATCATTCACGGACAACACTTTTT 8296 
hsa-miR-539* AATCATACAAGGACAATTTCTTTTTGA 3332 
hsa-miR-181b-1* CTCACTGAACAATGAATGCAACT 1542 
hsa-miR-561* ATCAAGGATCTTAAACTTTGCC 1315 
hsa-let-7c* CTGTACAACCTTCTAGCTTTCCT 1195 
hsa-miR-652* ACAACCCTAGGAGAGGGTGCCATTCA 982 
hsa-miR-301a* GCTCTGACTTTATTGCACTACT 880 
hsa-miR-487a* GTGGTTATCCCTGCTGTGTTCG 823 
hsa-miR-370* AAGCCAGGTCACGTCTCTGCAGTTACAC 624 
hsa-miR-412* TGGTCGACCAGTTGGAAAGTAAT 578 
hsa-miR-376c* GTGGATATTCCTTCTATGTTTAT 568 
hsa-miR-381* AAGCGAGGTTGCCCTTTGTATATTC 567 
hsa-miR-376b* GTGGATATTCCTTCTATGTTTA 532 
New miRNASequenceExpression
hsa-miR-1307* CTCGACCGGACCTCGACCGGCTCGT 72 670 
hsa-miR-376a-2* GGTAGATTTTCCTTCTATGGTTA 14 296 
hsa-miR-382* CGAATCATTCACGGACAACACTTTTT 8296 
hsa-miR-539* AATCATACAAGGACAATTTCTTTTTGA 3332 
hsa-miR-181b-1* CTCACTGAACAATGAATGCAACT 1542 
hsa-miR-561* ATCAAGGATCTTAAACTTTGCC 1315 
hsa-let-7c* CTGTACAACCTTCTAGCTTTCCT 1195 
hsa-miR-652* ACAACCCTAGGAGAGGGTGCCATTCA 982 
hsa-miR-301a* GCTCTGACTTTATTGCACTACT 880 
hsa-miR-487a* GTGGTTATCCCTGCTGTGTTCG 823 
hsa-miR-370* AAGCCAGGTCACGTCTCTGCAGTTACAC 624 
hsa-miR-412* TGGTCGACCAGTTGGAAAGTAAT 578 
hsa-miR-376c* GTGGATATTCCTTCTATGTTTAT 568 
hsa-miR-381* AAGCGAGGTTGCCCTTTGTATATTC 567 
hsa-miR-376b* GTGGATATTCCTTCTATGTTTA 532 

Seven new miRNAs associated with read counts > 103.

TargetScan custom target predictions (see “New miRNA target prediction”) for 15 novel miRNAs with a read count more than 500 are included in supplemental Table 7. Specifically, conserved and nonconserved target sites were predicted using TargetScan and filtered according to the context score. Only conserved target sites associated with top 25% scores and nonconserved sites included in top 5% scores were reported.

Sister miRNA expression prevalence

Considering known and new miRNAs expressed in SET2 cells as a whole, we found that both miRNA and miRNA* were expressed concurrently in 260 hairpins, corresponding to approximately one-half of those with at least one miRNA expressed (Tables 3 and 4). miRNA and miRNA* of the same hairpin, called a sister miRNA pair, have different sequences, thus targeting different sets of coding RNAs and contributing uniquely, but possibly in a coordinate way, to transcriptional regulation. When increasing thresholds of read count were applied to consider a miRNA as being “expressed,” the fraction of hairpins producing concurrently miRNA and miRNA* decreased (Tables 34), whereas the fraction of hairpins producing a meaningful quantity of both miRNAs remained considerable at all thresholds.

Table 3

Number of expressed miRNAs respective to number of known miRNAs, per hairpin

No. of known miRNAsHairpinsHairpins with 2 miRNAs expressed
Hairpins with 1 miRNA expressed
Hairpins with 0 miRNAs expressed
No.%No.%No.%
478 189 39.54 102 21.3 187 39.1 
943 71 7.53 175 18.6 697 73.9 
Hairpins 1421 260 18.30 277 19.5 884 62.2 
No. of known miRNAsHairpinsHairpins with 2 miRNAs expressed
Hairpins with 1 miRNA expressed
Hairpins with 0 miRNAs expressed
No.%No.%No.%
478 189 39.54 102 21.3 187 39.1 
943 71 7.53 175 18.6 697 73.9 
Hairpins 1421 260 18.30 277 19.5 884 62.2 
Table 4

Percentages of hairpins associated with concurrently expressed miRNA/miRNA* pairs, according to different thresholds on expression level

Expression of at leastHairpins with at least 1 miRNA expressedHairpins with 2 miRNAs expressed% of expressed hairpins with 2 miRNAs
10 522 250 47.9 
102 331 175 52.9 
103 229 94 41.0 
104 125 15 12.0 
Expression of at leastHairpins with at least 1 miRNA expressedHairpins with 2 miRNAs expressed% of expressed hairpins with 2 miRNAs
10 522 250 47.9 
102 331 175 52.9 
103 229 94 41.0 
104 125 15 12.0 

We observed no strand prevalence in expressed miRNAs (supplemental Figure 4). Of 277 hairpins expressing only one miRNA, 47% and 53% expressed only the 5′- or 3′-miRNA (130 and 147), respectively. Regarding 260 hairpins expressing both sister miRNAs, the most expressed miRNA was the 5′- or 3′-form in 135 and 125, respectively. A considerable fraction of concurrently expressed miRNA pairs showed comparable levels: 66 (25%) and 38 (14%) of concurrently expressed pairs were associated with an absolute log2(ratio) of expression values not higher that 2 and 1, respectively.

moRNA identification in known hairpin precursors

We discovered ERE also outside known and novel miRNAs. These were classified, as explained in supplemental Methods, as 5′-moRNAs, 3′-moRNAs, and expressed loops (Figure 4). In particular, we identified 58 moRNAs expressed from 56 hairpins, at moderate to high level (mean 127, median 52; supplemental Tables 7-8). The expression level of 12 and 4 moRNAs exceeded 100 and 500, respectively (Table 5). Three hairpins were associated each with 2 moRNAs expressed (supplemental Table 9, Figure 4); of these, hsa-mir-106b is an example of a “5-phased” precursor, with detected expression for 5 different regions. Figure 5 shows, for extended hairpins of the 4 most expressed new moRNAs, the plots of number of reads aligned against nucleotide position of extended hairpin sequences and predicted moRNA sequences.

Figure 4

Two miRNAs and 2 moRNAs may be produced by transcription and processing of a single miRNA locus. The plots show expression levels of miRNAs, moRNAs, and loops from 15 hairpins corresponding to most expressed moRNAs.

Figure 4

Two miRNAs and 2 moRNAs may be produced by transcription and processing of a single miRNA locus. The plots show expression levels of miRNAs, moRNAs, and loops from 15 hairpins corresponding to most expressed moRNAs.

Table 5

Name and predicted sequence of 16 new moRNAs with read count > 100

moRNASequenceExpression
hsa-5′-moR-103a-2* GAGCTGCGTCTTTGTGCTTTC 1098 
hsa-5′-moR-106b* CCGCTCCAGCCCTGCCGGGGC 918 
hsa-5′-moR-19b-1* TGTTACTGAACACTGTTCTATGGTT 700 
hsa-5′-moR-16-1* TAGCAATGTCAGCAGTGCCT 506 
hsa-5′-moR-20a TGATGTGACAGCTTCTGTAGCAC 418 
hsa-5′-moR-21 CATGGCTGTACCACCTTGTCGGG 345 
hsa-5′-moR-27a CTGTGCCTGGCCTGAGGAGC 306 
hsa-5′-moR-125a ACCATGTTGCCAGTCTCTAGG 300 
hsa-5′-moR-377 CCGTGCTGATGTTTGACCCTTGAGCA 273 
hsa-3′-moR-103a-2 AAGAACCAAGAATGGGCTGCC 206 
hsa-5′-moR-941-3 CACCCGGCTGTGTGCACATGTGC 155 
hsa-5′-moR-361 TTTTCCTGGGATTTGGGAGC 154 
hsa-5′-moR-92a-1 AACTCAAACCCCTTTCTACAC 125 
hsa-5′-moR-19a TTTGTTTGCAGTCCTCTGTT 119 
hsa-5′-moR-7-1 GGATGTTGGCCTAGTTCTGTG 106 
hsa-5′-moR-98 GGATTCTGCTCATGCCAGGG 102 
moRNASequenceExpression
hsa-5′-moR-103a-2* GAGCTGCGTCTTTGTGCTTTC 1098 
hsa-5′-moR-106b* CCGCTCCAGCCCTGCCGGGGC 918 
hsa-5′-moR-19b-1* TGTTACTGAACACTGTTCTATGGTT 700 
hsa-5′-moR-16-1* TAGCAATGTCAGCAGTGCCT 506 
hsa-5′-moR-20a TGATGTGACAGCTTCTGTAGCAC 418 
hsa-5′-moR-21 CATGGCTGTACCACCTTGTCGGG 345 
hsa-5′-moR-27a CTGTGCCTGGCCTGAGGAGC 306 
hsa-5′-moR-125a ACCATGTTGCCAGTCTCTAGG 300 
hsa-5′-moR-377 CCGTGCTGATGTTTGACCCTTGAGCA 273 
hsa-3′-moR-103a-2 AAGAACCAAGAATGGGCTGCC 206 
hsa-5′-moR-941-3 CACCCGGCTGTGTGCACATGTGC 155 
hsa-5′-moR-361 TTTTCCTGGGATTTGGGAGC 154 
hsa-5′-moR-92a-1 AACTCAAACCCCTTTCTACAC 125 
hsa-5′-moR-19a TTTGTTTGCAGTCCTCTGTT 119 
hsa-5′-moR-7-1 GGATGTTGGCCTAGTTCTGTG 106 
hsa-5′-moR-98 GGATTCTGCTCATGCCAGGG 102 
*

Four moRNAs associated with read counts > 500.

Figure 5

Selected examples of moRNA expression and sequence. Plots show the number of reads per nucleotide position mapping in extended hairpin loci expressing miRNAs and moRNAs.

Figure 5

Selected examples of moRNA expression and sequence. Plots show the number of reads per nucleotide position mapping in extended hairpin loci expressing miRNAs and moRNAs.

We observed a considerable prevalence of 5′-moRNAs: 95% of moRNAs derived from the 5′-arm of the hairpin precursor (n = 55 vs 5 for the 3′-moRNAs). The average expression level of 5′-moRNAs (129) was higher than 3′-moRNAs (89), although it did not reach the significance level. Of interest, the prevalence of 5′-moRNAs was not dependent on the type of the most expressed miRNA (supplemental Figure 5A).

Expressed miRNA length ranged from 16 to 27 nt, with an average of 21, at variance with moRNA sequences that were less variable in length, ranging from 18 to 25 nt (average, 20.2 nt).

To derive insights concerning moRNA biogenesis, we considered the position of moRNA sequences relatively to both the miRNA hairpin precursor (Figure 6A) and the mature miRNA from the same hairpin arm (Figure 6B). Regarding the region of the hairpins from which moRNAs derived, we noticed that 8 of 55 5′-moRNAs were included in the “classic” hairpin precursor sequence, whereas only one belonged to the region of 30 nt flanking the hairpin and 46 (84%) were partially overlapping the hairpin 5′-border. This indicates that moRNA sequence spanned the canonical Drosha cutting site and supports a role for noncanonical Drosha cleavage in moRNAs biogenesis. Moreover, we considered the positions of 50 couples of concurrently expressed 5′-moRNAs and 5′-miRNAs in extended hairpins. It is noteworthy that 5′-moRNAs and 5′-miRNAs were adjacent in 27 pairs (54%): the miRNA sequence started exactly after the moRNA end. In this case, the 2 elements might derive from a single cutting step of a common precursor sequence. Conversely, in only 1 pair (2%), the 2 sequences were disjointed and the miRNA sequence started at least 1 nt after the moRNA end. Finally, 22 moRNAs (44%) overlapped with the nearby miRNA sequence for 1 to 8 nt, indicating the presence of pairs of moRNAs and miRNAs that were probably produced alternatively from the processing of the same precursor.

Figure 6

Analysis of 5′-moRNA position provides clues about moRNA biogenesis. 5′-moRNA position relative to miRNA hairpin precursor (A) and to the mature miRNA expressed from the same hairpin arm (B). The number of expressed moRNA falling in each category is indicated close to the corresponding example.

Figure 6

Analysis of 5′-moRNA position provides clues about moRNA biogenesis. 5′-moRNA position relative to miRNA hairpin precursor (A) and to the mature miRNA expressed from the same hairpin arm (B). The number of expressed moRNA falling in each category is indicated close to the corresponding example.

miRNA and moRNA evolutionary conservation

Mean evolutionary sequence conservation in placental mammals of both detected miRNAs and moRNAs in SET2 cells varied widely from 0.0 to 1.0, representing extreme values of the conservation measure. Both distributions of conservation values are considerably negatively skewed. miRNA observed mean and median conservation values were 0.68 and 0.99, respectively. Supplemental Figure 6 shows the relationship between miRNA expression and average sequence conservation in placental mammals. Despite high conservation variation, highly expressed miRNAs do not include poorly conserved elements. Conservation levels of 5′- and 3′-miRNA groups were almost the same (P = .83). Besides, the paired comparison of 5′- and 3′-miRNAs of the same hairpin highlighted a slight but significant (P = .001) higher conservation of 3′-miRNAs. Mean and median moRNA conservation scores were 0.80 and 0.99, respectively, not significantly different from miRNAs (P = .165).

Changes in miRNA expression in SET2 cells treated with ruxolitinib

Information concerning expression of different miRNA species in SET2 cells might represent a starting point for studies aimed at identifying set of miRNas that lie in the JAK2 regulatory pathway. To this end, SET2 cells were treated with the JAK2 inhibitor INC424/ruxolitinib and subjected to miRNA quantification using microarrays. INC424/ruxolitinib caused a dramatic reduction of phophoSTAT5 and Pim-1 mRNA, know targets of activated JAK2 (Figure 7) supporting the effectiveness of inhibition. Using stringent selection criteria (supplemental Methods), we identified 9 microRNAs among those identified by RNASeq in untreated cells that showed significant dose-dependent modulated expression (hsa-let-7c, hsa-miR-1299, hsa-miR323-3p were up-regulated, while hsa-miR-330-5p, hsa-miR-548, hsa-miR-665, hsa-miR-935u, hsa-miR-381, hsa-miR-92a-1* were down-regulated; Figure 7). miRNA expression data are available at GEO (GSE33809).

Figure 7

Modulation of miRNAs by INC424/ruxolitinib treatment of SET2 cells. Cells were exposed to INC414/ruxolitinib at 160, 80, and 1600nM for 3 and 6 hours and then processed for miRNA expression analysis by microarray. The concentration of 160nM was selected as the IC50 concentration in a proliferation assay of SET2 cells (not shown). INC424/ruxolitinib caused dose-dependent reduction of the level of phosphorylated STAT5 (measured by FACS) and Pim-1 mRNA (measured by quantitative real-time PCR), indicating effective inhibition of JAK/STAT signaling (top panels). (Bottom panels) Examples of INC424/ruxolitinib-induced miRNA level modulation, as measured using miRNA array, are shown. The levels of hsa-miR-92a-1* and hsa-miR-330-5p dose-dependently decreased, whereas those of has-miR-let-7c increased compared with untreated control cells. Dose effects were statistically significant (P < .05).

Figure 7

Modulation of miRNAs by INC424/ruxolitinib treatment of SET2 cells. Cells were exposed to INC414/ruxolitinib at 160, 80, and 1600nM for 3 and 6 hours and then processed for miRNA expression analysis by microarray. The concentration of 160nM was selected as the IC50 concentration in a proliferation assay of SET2 cells (not shown). INC424/ruxolitinib caused dose-dependent reduction of the level of phosphorylated STAT5 (measured by FACS) and Pim-1 mRNA (measured by quantitative real-time PCR), indicating effective inhibition of JAK/STAT signaling (top panels). (Bottom panels) Examples of INC424/ruxolitinib-induced miRNA level modulation, as measured using miRNA array, are shown. The levels of hsa-miR-92a-1* and hsa-miR-330-5p dose-dependently decreased, whereas those of has-miR-let-7c increased compared with untreated control cells. Dose effects were statistically significant (P < .05).

Discussion

In this study, we characterized the fraction of short RNAs expressed in SET2 cells, a JAK2V617F-mutated cell line established from a patient with leukemic transformation of essential thrombocythemia, by exploiting deep sequencing data analyzed in integration with different levels of genomic information and metadata concerning known hairpin loci. We focused on known miRNA loci by considering 1421 known hairpin precursors and 1731 known mature miRNAs, including hundreds of new miRNAs discovered by massive sequencing approaches. Sequencing data were processed and filtered with very stringent criteria, based on sequence length and quality. Read mapping was obtained using a quality-based method first and with a new method based on comparative mapping of reads to the human genome and to hairpin precursors later to filter out contaminant RNAs from repetitive loci.

One accomplishment of this study is the production of the first catalog of known miRNAs in SET2 cells and the identification of most expressed ones, which are conceivably responsible of most miRNA-mediated gene repression. Considering predicted targets of highly expressed miRNAs, different signaling (and biochemical) pathways were identified, obtaining interesting clues about possible regulatory circuits that are deregulated in SET2 cells (supplemental Table 4). For instance, we noticed that several among the top 25% expressed miRNAs belonged to 2 clusters, the miR-106b-25 (located on chromosome 7q22.1, which includes miR106b, miR-93, and miR-25) and the miR-17-92 (located on chromosome 13q31.3, which includes miR-17, miR-19b, miR-20a, and miR-92a, miR19a, and miR18b); 6 of these miRNAs were among the 21 mostly expressed miRNAs. Recent evidence support the contention that uncontrolled expression of both or either these 2 clusters results in the loss of physiologic control of cell cycle arrest and apoptosis by TGF-β contributing to oncogenesis.31 

The expression of known mature miRNAs from sequencing data has been estimated with a method able to correct for multiple mapping issues. Read multiple mapping might arise from several reasons largely attributable to miRNA loci redundancy and similarity, to the existence of miRNA families, as well to reads-to-reference sequence mapping choices. The problem of mapping quality was considerably overlooked by different studies that recently exploited RNAseq for miRNA and isomiR discovery. It has been shown32  that short reads are prone to map to multiple genomic loci with an equal number of mismatches (or even without mismatches), in particular among multicopy miRNA precursors (identical mature miRNAs may be produced by different hairpins, transcribed from different loci) and miRNA families. Other studies considered only reads with unique mapping,33  although too restrictive mapping settings may affect quantitative estimation of multicopy miRNAs. Basically, we implemented a “multiple-mapping corrected read count” to quantify correctly expression levels relative to both hairpin and mature miRNAs, also including multilocus miRNAs and miRNA families.

Another main finding of the study regards mature miRNA sequence variation. We showed that mature miRNAs were not represented by a unique sequence (ie, the one annotated in miRBase); rather, they were found as mixtures of sequence variants, called isomiRs, that derive mainly from noncanonical processing of hairpin precursors. Only one-third of expressed miRNAs were associated with a unique isomiR corresponding to known mature miRNA sequence, with two-thirds being associated with different sequences. Surprisingly, in 33% of the miRNAs with at least 2 isomiRs, the most expressed isomiR was not identical to the miRBase-annotated mature miRNA sequence. Using a conservative approach by focusing only on highly expressed miRNAs and considering merely those isomiRs accounting for at least 10% of total miRNA expression, we showed that 60% of miRNAs are associated with 2 to 4 isomiRs, which differ from the most expressed isomiR only for 5′- and 3′-sequence length but align exactly to the hairpin precursor. This might imply that most isomiRs arise as a consequence of alternative processing of miRNA precursors. Interestingly, the sequence variability of isomiRs regards prevalently the 3′-region of mature miRNAs; however, 8% of isomiRs differ from the “classic” miRNA sequence in the 5′-region, possibly impacting on target recognition and/or regulatory activity and strength. Future evaluations of isomiR quality and proportion in different cell types will help to clarify if and how much this miRNA biogenesis feature is cell- and context-specifically regulated.

In the second part of the study, we exploited the discovery power of RNAseq to identify a consistent number of novel short RNAs expressed from miRNA loci in SET2 cells. It is noteworthy that RNA discovery results were produced using only read aligning exactly to hairpin loci to minimize artifacts. At variance with other studies, we mapped sequence reads to extended hairpin loci (ie, the hairpin precursor genomic regions plus 30 nt upstream and downstream). This allowed us to uncover novel short RNAs derived from extended hairpin loci transcription and canonical or noncanonical processing. More specifically, we discovered 78 novel miRNAs expressed from known hairpin precursors. These are all new miRNAs*, 11 of which are expressed in SET2 cells over the median value observed for the group of known miRNAs. These data confirm recent findings about the importance of miRNAs*.18  Considering known and new miRNAs together, we took into account expression behavior of pairs of miRNAs derived from 5 and 3′-strands of the same hairpin (miRNA/miRNA*, so-called “sister miRNAs”). Sister miRNAs have different sequences and target different sets of coding RNAs. We showed that for approximately one-half of hairpins both sister miRNAs are expressed concurrently; very likely, both contribute to target repression. Moreover, slightly less than one-fourth of concurrently detected sister miRNA pairs were expressed at the same or comparable level. No strand prevalence was observed. These results further support theories about the nondeterministic nature of strand selection bias, which appears to be regulated in a cell-, tissue-, and condition-specific way.19  In addition, our findings are compatible with the 2-step cleavage of hairpin RNA by Dicer16 ; bidirectional binding of processed dsRNAs by Dicer may indeed result in directional presentation of double-strand miRNA duplex to Argonaute, influencing the strand selection bias.

Interestingly, we found short expressed RNAs derived from regions of extended hairpins outside known and new miRNAs. These are members of a novel class of miRNA-related RNAs, called moRNAs,25  recently identified by massive short RNAseq. moRNAs are currently defined as approximately 20-nt-long RNAs that originate predominantly from the 5′-arm of pre-miRNAs, with a biogenesis linked to that of miRNAs but not necessarily interdependent. These studies also indicated that moRNAs are generally included in the miRNA hairpin precursor; and although the moRNA usually overlaps the miRNA position by a few nucleotides, sometimes it overhangs the miRNA hairpin. In the present study, we identified 58 moRNAs expressed from 56 hairpins at moderate to high level. In particular, hsa-5′-moR-103a-2, hsa-5′-moR-106b, hsa-5′-moR-19b-1, and hsa-5′-moR-16-1 were highly expressed in SET2 cells. It was also of interest that 5 of the 16 most expressed moRNAs (Table 5) were associated with the 17-92 (miR-19b-1, miR-20a, miR-92a, and miR-19a) and the 106-25 (miR-106b) cluster, as discussed in the second paragraph of this section for miRNAs. Of note, some of the highly expressed moRNAs were unique and did not perfectly overlap the most expressed miRNAs from the same cluster, suggesting that activation of the cluster might result in global up-regulation of the genes but also in a nonbalanced ratio between miRNAs and moRNAs (supplemental Table 10).

We observed a prevalence of moRNAs derived from the 5′-arm of hairpin precursors, similar to previous reports, and showed that the prevalence of 5′-moRNA is independent from miRNA strand selection bias. We found moRNAs longer than previously reported, with an average length of 20.2 nt and a maximum of 28 nt. A large majority (84%) of 5′-moRNAs overlapped the 5′-miRNA hairpin precursor boundary. Because moRNAs span the canonical Drosha cutting site, noncanonical Drosha cleavage might conceivably be responsible for moRNA production. Moreover, we observed that slightly less than one-half of moRNAs overlapped the nearby expressed mature miRNA of a few nucleotides; in these cases, moRNA and miRNA expression from the same precursor was mutually exclusive. Among the 9 most expressed moRNAs, 5 overlapped the corresponding miRNA. Interestingly, the expression ratio between 5′-moRNA and miRNA pairs was much variable: the hsa-miR-103a was expressed almost twice the hsa-5′-moR-103a-2, the most expressed moRNA, and the 2 expression values were of the same degree of magnitude (2809 vs 1098, respectively). A similar situation was observed for hsa-mir-19b-1 and hsa-mir-16-1. Conversely, other moRNAs that are produced from the same hairpin of miRNAs were expressed at definitely lower levels: expression of hsa-miR-21 was approximately 6000 times higher than the hsa-5′-moR-21. The observation that moRNA sequences are as conserved as miRNAs gives additional, although indirect, support to the hypothesis that moRNAs plays a defined biologic role. Evidence regarding possible functions of moRNAs is still fragmentary25 ; one hypothesis claims that moRNAs might guide RISC to complementary target mRNAs as miRNAs do.34  Nevertheless, the fact that moRNAs are nuclear enriched may support the alternative hypothesis that moRNAs intervene specifically in nuclear processes as other nuclear short and long RNAs do. Therefore, moRNAs can be viewed as a new class of regulators whose qualitative and/or quantitative abnormalities might impact on human diseases. In this view, the discovery of expressed moRNAs in SET2 cells is intriguing.

In conclusion, we have extensively characterized the profile of short RNAs expressed in SET2 cells and provided proof of concept that modulation of mature miRNA expression may be obtained via inhibition of JAK2. It must be stressed, however, that SET2 cells are not a physiologic model of MPN, so any novel information should be validated in primary MPN cells. It will be key to understand how the interplay of known and new miRNAs, isomiRs, and moRNAs contributes to the abnormal regulation of cell proliferation that characterizes MPN cells and whether and how drugs affect this complex system of regulators. Finally, results regarding isomiRs and moRNA expression fit well in the “RNA in pieces” phenomenon35 : many transcripts undergo posttranscriptional cleavage to release specific, functionally independent fragments, expanding the spectrum of regulatory RNAs produced by the human genome and from single loci.

This article contains a data supplement.

The publication costs of this article were defrayed in part by page charge payment. Therefore, and solely to indicate this fact, this article is hereby marked “advertisement” in accordance with 18 USC section 1734.

Acknowledgments

This work was supported by AIRC “Special Program Molecular Clinical Oncology 5 X 1000” to AGIMM (project 1005; a detailed description of the AGIMM project is available at http://www.progettoagimm.it). A.M.V. was supported by the AIRC (project 9034) and Ministero dell'Istruzione, dell'Università e della Ricerca (MIUR)–20084XBENM_002). S.B. was supported by Fondazione Cassa di Risparmio di Padova e Rovigo and the University of Padova. R.N. was supported by the Collegio Ghislieri Foundation (fellowship).

Authorship

Contribution: S.B. and A.M.V. designed research, analyzed data, interpreted results, and wrote the manuscript; A.B. designed and implemented the computational pipeline, analyzed data, and wrote the manuscript; S.B. implemented the computational pipeline, analyzed data, and wrote the manuscript; M.B. prepared target predictions and wrote the manuscript; P.G. performed research, interpreted results, and wrote the manuscript; F.B. and R.N. performed research; and R.M. performed and interpreted array expression experiments.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Alessandro M. Vannucchi, Section of Hematology, Department of Critical Care, University of Florence, Largo Brambilla 3, 50134 Florence, Italy; e-mail: amvannucchi@unifi.it.

References

References
1
Vannucchi
 
AM
Guglielmelli
 
P
Tefferi
 
A
Advances in understanding and management of myeloproliferative neoplasms.
CA Cancer J Clin
2009
, vol. 
59
 
3
(pg. 
171
-
191
)
2
Levine
 
RL
Pardanani
 
A
Tefferi
 
A
Gilliland
 
DG
Role of JAK2 in the pathogenesis and therapy of myeloproliferative disorders.
Nat Rev Cancer
2007
, vol. 
7
 
9
(pg. 
673
-
683
)
3
Tefferi
 
A
Novel mutations and their functional and clinical relevance in myeloproliferative neoplasms: JAK2, MPL, TET2, ASXL1, CBL, IDH and IKZF1.
Leukemia
2010
, vol. 
24
 
6
(pg. 
1128
-
1138
)
4
Vainchenker
 
W
Delhommeau
 
F
Constantinescu
 
SN
Bernard
 
OA
New mutations and pathogenesis of myeloproliferative neoplasms.
Blood
2011
, vol. 
118
 
7
(pg. 
1723
-
1735
)
5
Li
 
J
Kent
 
DG
Chen
 
E
Green
 
AR
Mouse models of myeloproliferative neoplasms: JAK of all grades.
Dis Model Mech
2011
, vol. 
4
 
3
(pg. 
311
-
317
)
6
James
 
C
Ugo
 
V
Le Couedic
 
JP
, et al. 
A unique clonal JAK2 mutation leading to constitutive signalling causes polycythaemia vera.
Nature
2005
, vol. 
434
 
7037
(pg. 
1144
-
1148
)
7
Vannucchi
 
AM
Guglielmelli
 
P
Rambaldi
 
A
Bogani
 
C
Barbui
 
T
Epigenetic therapy in myeloproliferative neoplasms: evidence and perspectives.
J Cell Mol Med
2009
, vol. 
13
 (pg. 
1437
-
1450
)
8
Tefferi
 
A
Abdel-Wahab
 
O
Cervantes
 
F
, et al. 
Mutations with epigenetic effects in myeloproliferative neoplasms and recent progress in treatment: proceedings from the 5th International Post-ASH Symposium.
Blood Cancer J
2011
, vol. 
1
 pg. 
e7
 
9
Bruchova
 
H
Yoon
 
D
Agarwal
 
AM
Mendell
 
J
Prchal
 
JT
Regulated expression of microRNAs in normal and polycythemia vera erythropoiesis.
Exp Hematol
2007
, vol. 
35
 
11
(pg. 
1657
-
1667
)
10
Bruchova
 
H
Merkerova
 
M
Prchal
 
JT
Aberrant expression of microRNA in polycythemia vera.
Haematologica
2008
, vol. 
93
 
7
(pg. 
1009
-
1016
)
11
Guglielmelli
 
P
Tozzi
 
L
Pancrazzi
 
A
, et al. 
MicroRNA expression profile in granulocytes from primary myelofibrosis patients.
Exp Hematol
2007
, vol. 
35
 
11
(pg. 
1708
-
1718
)
12
Lionetti
 
M
Biasiolo
 
M
Agnelli
 
L
, et al. 
Identification of microRNA expression patterns and definition of a microRNA/mRNA regulatory network in distinct molecular groups of multiple myeloma.
Blood
2009
, vol. 
114
 
25
(pg. 
e20
-
e26
)
13
Garzon
 
R
Marcucci
 
G
Croce
 
CM
Targeting microRNAs in cancer: rationale, strategies and challenges.
Nat Rev Drug Discov
2010
, vol. 
9
 
10
(pg. 
775
-
789
)
14
Guglielmelli
 
P
Tozzi
 
L
Bogani
 
C
, et al. 
Over-expression of microRNA-16-2 contributes to the abnormal erythropoiesis in polycythemia vera.
Blood
2011
, vol. 
117
 
25
(pg. 
6923
-
6927
)
15
Winter
 
J
Jung
 
S
Keller
 
S
Gregory
 
RI
Diederichs
 
S
Many roads to maturity: microRNA biogenesis pathways and their regulation.
Nat Cell Biol
2009
, vol. 
11
 
3
(pg. 
228
-
234
)
16
Burroughs
 
AM
Ando
 
Y
de Hoon
 
ML
, et al. 
Deep-sequencing of human Argonaute-associated small RNAs provides insight into miRNA sorting and reveals Argonaute association with RNA fragments of diverse origin.
RNA Biol
2011
, vol. 
8
 
1
(pg. 
158
-
177
)
17
Ro
 
S
Park
 
C
Jin
 
J
Sanders
 
KM
Yan
 
W
A PCR-based method for detection and quantification of small RNAs.
Biochem Biophys Res Commun
2006
, vol. 
351
 
3
(pg. 
756
-
763
)
18
Kuchenbauer
 
F
Mah
 
SM
Heuser
 
M
, et al. 
Comprehensive analysis of mammalian miRNA* species and their role in myeloid cells.
Blood
2011
, vol. 
118
 
12
(pg. 
3350
-
3358
)
19
Biasiolo
 
M
Sales
 
G
Lionetti
 
M
, et al. 
Impact of host genes and strand selection on miRNA and miRNA* expression.
PLoS One
2011
, vol. 
6
 
8
pg. 
e23854
 
20
Kozomara
 
A
Griffiths-Jones
 
S
miRBase: integrating microRNA annotation and deep-sequencing data.
Nucleic Acids Res
2011
, vol. 
39
 
suppl 1
(pg. 
D152
-
D157
)
21
Jima
 
DD
Zhang
 
J
Jacobs
 
C
, et al. 
Deep sequencing of the small RNA transcriptome of normal and malignant human B cells identifies hundreds of novel microRNAs.
Blood
2010
, vol. 
116
 
23
(pg. 
e118
-
e127
)
22
Vaz
 
C
Ahmad
 
H
Sharma
 
P
, et al. 
Analysis of microRNA transcriptome by deep sequencing of small RNA libraries of peripheral blood.
BMC Genomics
2010
, vol. 
11
 
1
pg. 
288
 
23
Martí
 
E
Pantano
 
L
Bañez-Coronel
 
M
, et al. 
A myriad of miRNA variants in control and Huntington's disease brain regions detected by massively parallel sequencing.
Nucleic Acids Res
2010
, vol. 
38
 
20
(pg. 
7219
-
7235
)
24
Fernandez-Valverde
 
SL
Taft
 
RJ
Mattick
 
JS
Dynamic isomiR regulation in Drosophila development.
RNA
2010
, vol. 
16
 
10
(pg. 
1881
-
1888
)
25
Bortoluzzi
 
S
Biasiolo
 
M
Bisognin
 
A
MicroRNA–offset RNAs (moRNAs): by-product spectators or functional players?
Trends Mol Med
2011
, vol. 
17
 
9
(pg. 
473
-
474
)
26
Shi
 
W
Hendrix
 
D
Levine
 
M
Haley
 
B
A distinct class of small RNAs arises from pre-miRNA-proximal regions in a simple chordate.
Nat Struct Mol Biol
2009
, vol. 
16
 
2
(pg. 
183
-
189
)
27
Langenberger
 
D
Bermudez-Santana
 
C
Hertel
 
J
Hoffmann
 
S
Khaitovich
 
P
Stadler
 
PF
Evidence for human microRNA-offset RNAs in small RNA sequencing data.
Bioinformatics
2009
, vol. 
25
 
18
(pg. 
2298
-
2301
)
28
Taft
 
RJ
Simons
 
C
Nahkuri
 
S
, et al. 
Nuclear-localized tiny RNAs are associated with transcription initiation and splice sites in metazoans.
Nat Struct Mol Biol
2010
, vol. 
17
 
8
(pg. 
1030
-
1034
)
29
Meiri
 
E
Levy
 
A
Benjamin
 
H
, et al. 
Discovery of microRNAs and other small RNAs in solid tumors.
Nucleic Acids Res
2010
, vol. 
38
 
18
(pg. 
6234
-
6246
)
30
Papadopoulos
 
GL
Alexiou
 
P
Maragkakis
 
M
Reczko
 
M
Hatzigeorgiou
 
AG
DIANA-mirPath: integrating human and mouse microRNAs in pathways.
Bioinformatics
2009
, vol. 
25
 
15
(pg. 
1991
-
1993
)
31
Petrocca
 
F
Vecchione
 
A
Croce
 
CM
Emerging role of miR-106b-25/miR-17-92 clusters in the control of transforming growth factor beta signaling.
Cancer Res
2008
, vol. 
68
 
20
(pg. 
8191
-
8194
)
32
Guo
 
L
Liang
 
T
Lu
 
Z
A comprehensive study of multiple mapping and feature selection for correction strategy in the analysis of small RNAs from SOLiD sequencing.
Biosystems
2011
, vol. 
104
 
2
(pg. 
87
-
93
)
33
Berezikov
 
E
Robine
 
N
Samsonova
 
A
, et al. 
Deep annotation of Drosophila melanogaster microRNAs yields insights into their processing, modification, and emergence.
Genome Res
2011
, vol. 
21
 
2
(pg. 
203
-
215
)
34
Umbach
 
JL
Strelow
 
LI
Wong
 
SW
Cullen
 
BR
Analysis of rhesus rhadinovirus microRNAs expressed in virus-induced tumors from infected rhesus macaques.
Virology
2010
, vol. 
405
 
2
(pg. 
592
-
599
)
35
Tuck
 
AC
Tollervey
 
D
RNA in pieces.
Trends Genet
2011
, vol. 
27
 
10
(pg. 
422
-
432
)

Author notes

*

S.B. and A.B. contributed equally to this study.