A distinct class of small RNAs arises from pre-miRNA–proximal regions in a simple chordate

Distinct small RNAs encoded by miRNA loci

We prepared small RNA (∼16–26-nt) libraries from C. intestinalis at various developmental stages, including unfertilized eggs, early embryos, late embryos and adults. High-throughput sequencing of the resulting cDNAs was performed with an Illumina 1G Genome Analyzer. Combining earlier studies with a recently described miRNA-discovery algorithm, we defined 80 miRNA loci in the C. intestinalis genome^16–18. Detailed information regarding the encoded miRNAs and their potential target mRNAs is provided in Supplementary Tables 1–4 online and at the following website: http://flybuzz.berkeley.edu/cgi-bin/CionaMicroRNAs.cgi.

Half of these genes encode a single major product (the miRNA), along with a less abundant miRNA* sequence, as is typically seen in other organisms^19,20. For example, the C. intestinalis (Ci) miR-125 gene (ortholog of the prototypic lin-4 miRNA in Caenorhabditis elegans) encodes a predominant miRNA that is stably expressed at all developmental stages examined²¹ (Supplementary Fig. 1 online). Ci-miR-125 is most highly expressed in adults, and at the adult stage a single clone of miR-125* is also detected.

Unexpectedly, the remaining half of C. intestinalis miRNA loci encode previously uncharacterized small RNAs, in addition to conventional miRNA and miRNA* products. This new class of RNAs arises from sequences located adjacent to the predicted pre-miRNA stem-loop, and we hereafter refer to them as ‘moRs’, for miRNA-offset RNAs. Only small RNAs with 5′ monophosphates and free 3′ hydroxyl groups can be cloned by the method used in this study (see Methods), although they could contain modifications on the 2′ oxygen, as seen for Piwi-interacting RNAs (piRNAs) and some miRNA species²². Most moR sequences are 19–20 nt in length, whereas C. intestinalis miRNAs range in size between 19 nt and 22 nt (Supplementary Fig. 2a online). Overall, moRs are considerably less abundant than miRNAs, but just ∼50% less abundant than miRNA* sequences (1,552 total reads and 3,353 total reads, respectively) (Supplementary Table 4b). In general, moRs show greater 5′ heterogeneity than miRNA or miRNA* sequences (Supplementary Fig. 2b). However, several abundantly expressed moRs, such as 5′ moRs 124-1 and 219, contain a rigid 5′-terminal nucleotide identity and show developmental regulation, suggesting that particular moRs may be under selective pressure, as has been suggested for the 5′ ends of miRNAs²³.

It is possible that the C. intestinalis miRNA loci encoding moRs contain unique structural features, as compared to those that do not²⁴. Global comparisons of base-pairing probabilities across the extended pre-miRNA loci in C. intestinalis revealed only modest structural differences between the two classes of miRNA loci (Supplementary Fig. 3 online). Overall, C. intestinalis miRNA loci maintain a similar base-pairing probability trace as those seen in Drosophila melanogaster, suggesting that C. intestinalis miRNA genes lack an intrinsic, species-specific structure. Similarly, there is no obvious difference in the size of the loop sequences in pre-miRs that produce moRs and those that do not (∼13 nt and ∼15 nt, respectively; Supplementary Fig. 4a online). In addition, we analyzed sequence motifs for all small RNAs cloned in this study. Whereas C. intestinalis miRNAs retained the expected 5′-uracil bias, no obvious motifs were apparent in the moRs²⁵ (Supplementary Fig. 4b). Thus, it is currently unclear why they arise from particular miRNA loci.

Defining characteristics of moRs

The C. intestinalis miR-219 gene (Ci-miR-219) encodes a predicted 57-nt pre-miRNA hairpin that is processed to produce miR-219 and miR-219*. In addition, a 5′ moR product (5′-moR-219) arises from sequences located immediately adjacent to miR-219 (Fig. 1a,b). The predominant miR-219 and miR-219* sequences are each ∼21 nt in length, whereas the 5′-moR-219 sequence is 20 nt (Fig. 1b).

Like most miRNAs, nearly all 5′-moR-219 clones maintain an invariant 5′ end (223 of 232 total reads)²³. Each of the three small RNAs observed at the miR-219 locus showed developmental regulation (Fig. 1a). Only miR-219 was detected in unfertilized eggs and adults (Fig. 1a), whereas both miR-219* and 5′-moR-219 were seen during embryogenesis.

In some cases, two distinct moRs are produced from a single miRNA gene, in addition to miRNA and miRNA* sequences (Fig. 2). The Ci-miR-124 locus encodes a pri-miRNA containing two tandem, but slightly different, ∼58-nt pre-miRNAs (Fig. 2b). The resulting miRNAs, miR-124-1 and miR-124-2, are identical, and the sequence shows peak expression, as evidenced by increased read counts (see miR-133 example below), in advanced-stage embryos (Fig. 2a). Both pre-miRNAs produce 5′ and 3′ moRs during embryogenesis (Fig. 2a). We observed the 3′ moR from the pre–miR-124-2 hairpin in both early embryos and late embryos, but the 3′ moR from the pre-miR-124-1 hairpin was detected only in early embryos. Moreover, 5′-moR-124 RNAs are considerably more abundant than the 3′-moR-124 RNAs, a result that is typical of the moRs and reminiscent of the processing of miRNA and miRNA* sequences, as well as processing of pri-miRNA 5′ and 3′ arms by Drosha^26–28.

Notably, alignment of coincident 5′ and 3′ moR sequences from numerous miRNA loci suggests that they arise from RNAse III processing (∼21-nt duplexed RNAs with ∼2-nt 3′ overhangs)²⁹ (Fig. 2c and Supplementary Fig. 5 online).

Despite the high prevalence of moRs associated with miRNA loci in C. intestinalis, we found that, overall, moR sequences are poorly conserved as compared to miRNAs between C. intestinalis and a related ascidian species, Ciona savignyi, and moRs are even less conserved than miRNA* sequences (Supplementary Fig. 4c). However, it has been noted that well-conserved small RNAs are expressed at higher levels than those lacking conservation^19,30. This is true for most miRNAs when comparing C. intestinalis to C. savignyi (Supplementary Fig. 4c). Similarly, abundant moRs are also better conserved than those found at low copy number. Nonetheless, the general lack of conservation raises the possibility that moRs may represent unstable processing intermediates during the biogenesis of miRNAs. Such intermediates might be produced through a generic RNA-degradation mechanism that leaves behind spurious and variably sized small RNAs. However, as with miRNAs, the high copy number and near uniformity of clones at each locus suggests that moRs are produced mainly as ∼20-nt RNAs. To further address this point, we used northern assays to directly examine the expression and size distribution of miRNA and moRs in C. intestinalis embryos (see below).

Direct detection of moRs as discrete small RNAs

Vertebrate miR-133 genes are often part of a bicistronic pri-miRNA that also contains miR-1, and the two miRNAs work together to promote mesodermal fates³¹. A similar genomic linkage is seen in C. intestinalis, and previous studies have shown that the primary transcript containing miR-1 and miR-133 is selectively expressed in developing tail muscles during C. intestinalis embryogenesis³². The C. intestinalis miR-133 locus encodes separate miRNA, miRNA* and 5′ moR products (Fig. 3). miR-133 reads steadily increase during embryogenesis and reach peak levels in adults (Fig. 3a). We found that the 5′-moR-133 RNA is most abundant in late embryos and is present at an equal or higher read count than miR-133 and miR-133* at all embryonic stages examined.

The levels of miR-133 and 5′-moR-133 detected in northern assays are in agreement with the sequencing frequencies obtained from the cDNA libraries (Fig. 3c). There is a progressive increase in the steady-state levels of miR-133 in unfertilized eggs, early embryos, late-stage embryos and adults (Fig. 3c, above). Similarly, the predicted 5′-moR-133 RNA was detected as a stable product (appearing as a doublet of ∼19–20-nt species in adults), with peak levels seen in late embryos. There was no indication of a smear or ‘ladder’ of higher- or lower-molecular-weight products, as would be expected if moRs represented incompletely degraded hairpin sequences or cleaved pri-miRNA transcripts. Moreover, ectopic expression of Ci-miR-133 directed by a Ci-Brachyury enhancer in the developing C. intestinalis notochord—the primitive chordate backbone—resulted in increased accumulation of both 5′-moR-133 and miR-133, indicating that expression of a discrete moR is correlated with that of the host miRNA transcript³³ (Fig. 3d.)

Drosophila pri-miRNAs produce moRs in the Ciona tadpole

The preceding analysis suggests that moRs arise from an intrinsic property of the C. intestinalis small RNA–biogenesis machinery (see Discussion). To test this possibility, the miR-309 miRNA cluster (also known as ‘8-miR’) from D. melanogaster was selectively expressed in C. intestinalis^34,35 (Fig. 4). We reasoned that the pri-miR-309 transcript would be more likely to produce detectable moRs when expressed in the C. intestinalis tadpole because it seems to produce such products, albeit rarely, in D. melanogaster (Fig. 4a).

We separately placed the entire miR-309 cluster under the control of three different tissue-specific enhancers from C. intestinalis that direct expression in the notochord, epidermis and mesenchyme, respectively^33,36. All three transgenes were coelectroporated into fertilized eggs, and the embryos were allowed to develop to the tailbud stage (after neurogenesis). Total RNA was extracted from these embryos and subjected to high-throughput sequencing or used for northern assays.

Drosophila melanogaster moRs are produced at high steady-state levels in C. intestinalis, and here we focused on the miR-3 and miR-5 genes within the miR-309 cluster. We detected only four 3′-moR-3 RNA reads in the D. melanogster embryo, whereas in C. intestinalis we observed nearly 2,000 copies (Fig. 4a,b). There is also a marked increase in the levels of the 5′-moR-5 RNA produced in C. intestinalis as compared with those in D. melanogaster. Nearly all copies of this moR RNA contain homogenous 5′ and 3′ termini (1,616 of 1,629 cloned copies are identical; Fig. 4b). In contrast, miR-3 was cloned at high frequency in D. melanogaster and C. intestinalis. Using northern assays, we identified similar levels of miR-3 in C. intestinalis and D. melanogaster embryos, a result that is consistent with the similar number of reads detected by sequence analysis. However, using a specific 5′-moR-3 hybridization probe, we detected a discrete band, without any obvious intermediate products, only in C. intestinalis embryos ectopically expressing the miR-309 cluster (Fig. 4c ).

There is no obvious correlation between the efficiency of moR biogenesis and the size of the loop sequence in the pre-miRNAs or conservation of other features. For example, the pre-miRNAs encoding miR-3 and miR-5 contain loops of 13 nt and 18 nt, respectively, but nonetheless produce similar yields of moRs. These experiments clearly demonstrate that the stable expression of moRs is an intrinsic feature of the C. intestinalis small RNA–processing machinery.

Small RNA cloning and detection

We collected adult C. intestinalis animals from Half Moon Bay, California, and maintained them in an artificial seawater tank. We carried out fertilization, dechorionation and electroporations as previously described³³. Total RNA was extracted from unfertilized eggs, cleavage stage, tadpole-stage embryos and adults using the miRVana miRNA Isolation Kit (Ambion). Small RNA cloning was carried out as previously described⁴¹. Basically, from ∼30 μg of total RNA, only 17–25-nt RNAs were size selected via 15% denaturing PAGE. The 3′‘modban-1’ adaptor (IDT) was ligated to the RNAs from this fraction with RNA ligase (Ambion) in ATP-free reaction buffer⁴¹, and appropriately ligated RNAs were size selected via 15% denaturing PAGE. The modified RNAs were subsequently ligated to a 5′ linker (Solexa linker) in the presence of RNA ligase and in reaction buffer with ATP. The resulting RNA library was reverse transcribed to a cDNA library with SuperScript II (Invitrogen). cDNA was amplified using Illumina sequencing–specific primers, and the resulting libraries were sequenced on an Illumina 1G Genome Analyzer. In parallel, small RNAs were extracted using TRIZOL (Invitrogen), cloned and sequenced, as above, from staged, 2–4-hour-old D. melanogaster Toll^10b embryos³⁴. Northern blotting assays were performed as described previously⁴².

We cloned the D. melanogaster miR-309 cluster by amplifying the locus from yw genomic DNA using pfuUltra High Fidelity polymerase (Stratagene) and the TOPO TA cloning system (Invitrogen).

Ci-Brachyury, Ci-FoxF and Ci-Twist enhancers were used to drive transgene expression in the C. intestinalis notochord, epidermis and mesenchyme, respectively^33,36.

Primers used for amplification of the Ci-Twist enhancer were Ci-Twist-F (forward), 5′-ACCACAGCTTCTATTATATA-3′, and Ci-Twist-R (reverse), 5′-CATCGTGTGTTGATTGATTT-3′.

Probe sequences for the Ci-miR-133 northern assay were Ci-miR-133 (5′-CAGCTGGTTGAAGGGGACCAAA-3′), Ci-5′-moR-133 (5′-GACCGACACCCGCAATGTTT-3′) and Ci-U6 (5′-GTCATCCTTGCGCAGGGGCCATGCTA ATCTTCTCTGTATCGTTCC-3′).

The C. intestinalis miR-133 amplification primers were Ci-miR-133-F (forward), 5′-CGTTTTATACGGTTATATACAGG-3′, and Ci-miR-133-R (reverse), 5′-TATTTCCGACTACTGAGCG-3′.

The Drosophila miR-309 cluster amplification primers were Dme-8miR-F (forward), 5′-TGCAGACAAATGACGAATTGA-3′, and Dme-8miR-R (reverse), 5′-CCGACCCTTTCAGGTAACAA-3′.

The probe sequences for the Drosophila miR-3 northern assay were Dme-miR-3, 5′-TGAGACACACTTTGCCCAGTGAT-3′ and Dme-5′-moR-3, 5′-CAGGATCGGGACCTTAGGTG-3′.

Data analysis

The standard Illumina pipeline (GAPipeline-0.3.0) was used to extract sequenced reads. Nucleotide positions 1 to 26 were aligned to the C. intestinalis (JGI version 1.0) or D. melanogaster (version 4.3) genomes using ELAND, and for the calculation of position-specific error rates^18,43. Supplementary Figure 6a online shows the average error rate, defined as the estimated probability of a base call being incorrect as a function of nucleotide position for each of the four lanes (libraries) studied. The error rate model for the Illumina pipeline was calibrated on the basis of uniquely aligned reads to the genome, and then applied to all reads. The average error rate (averaged over all reads) rises sharply beyond the twenty-first base, consistent with an assumption that the reads should be dominated by miRNA sequences of roughly 21 nt, as subsequent unaligned bases of the 3′ adaptor would be scored as low quality. Reads were trimmed so as to optimize the total nucleotide quality in a dynamic programming approach that produced trimmed reads such that the maximum acceptable error rate over the trimmed sequence is less than 10% (Q_PHRED = 10), the total quality of the read is optimized globally over all start and stop positions, and the resulting length is greater than or equal to 17 nt⁴⁴.

The trimming procedure can be described formally as follows. An optimal trimming can be achieved by defining a penalty P associated with making an incorrect base call at a given nucleotide n. Using the position-specific error probability, e_n, one can define an expected score for a given nucleotide as s_n = 1 ⋅ (1 – e_n)+P ⋅ e_n = 1 – (P +1) ⋅ e_n. The total expected score for a trimming of nucleotide sequence to start at position and end at position j is then given by:

One is then free to choose the penalty, such that the expected score is zero when the error rate is the maximum tolerated, so $P=\frac{1}{e_{\max }}-1$ and any error rate greater than e_max will produce a negative contribution to the score. A dynamic programming search then globally optimizes S(i, j) over all start and stop positions⁴⁴. Further details of the data analysis rationale and methodology are available in the Supplementary Methods online. A meta-analysis of the distribution for all processed reads across a miRNA locus is presented in Supplementary Figure 7 online.