Characterization and discovery of novel miRNAs and moRNAs in JAK2V617F-mutated SET2 cells

The JAK2V617F mutation, which occurs in most patients with polycythemia vera and more than or equal to 60% of essential thrombocythemia and primary myelofibrosis,¹  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.⁵  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.⁶  Recent information highlighted that deranged epigenetic gene regulation^7,8  and abnormal expression of microRNAs^9-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,¹²  miRNAs are considered prominent tumor markers, oncogenes (oncomiRNAs), relevant targets for therapy, or therapeutic agents depending on the specific context.¹³  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.¹⁴ 

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.¹⁵  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),¹⁶  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,²⁰  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.²⁵ 

moRNAs were first reported in the simple chordate ascidian Ciona intestinalis²⁶  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.²⁶  However, shortly thereafter, Langenberger et al reported moRNAs, derived from 78 miRNA loci, as being weakly expressed in human prefrontal cortex,²⁷  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.²⁸  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.²⁹ 

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.

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)

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 for known miRNAs was obtained using DIANA-miRPath webtool³⁰  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.

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).

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).

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).

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 10³ and 10⁴, 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.

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.

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.

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⁻⁰⁶; 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.

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.

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 3–4), whereas the fraction of hairpins producing a meaningful quantity of both miRNAs remained considerable at all thresholds.

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 log₂(ratio) of expression values not higher that 2 and 1, respectively.

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.

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.

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).

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).

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.³¹ 

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 shown³²  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,³³  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*.¹⁸  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.¹⁹  In addition, our findings are compatible with the 2-step cleavage of hairpin RNA by Dicer¹⁶ ; 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,²⁵  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 fragmentary²⁵ ; one hypothesis claims that moRNAs might guide RISC to complementary target mRNAs as miRNAs do.³⁴  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” phenomenon³⁵ : 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 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).

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.