Native RNA sequencing in fission yeast reveals frequent alternative splicing isoforms

Native sequencing of poly(A)⁺ RNAs in S. pombe

We extracted total RNA from S. pombe cells growing at log-phase and subsequently performed poly(A)⁺ selection. Then we performed dRNA-seq of the polyadenylated RNA using an ONT Gridion instrument (Garalde et al. 2018). We obtained a total of 7,297,642 dRNA-seq reads from four sequencing runs. Each of these reads corresponds to a single native poly(A)⁺ RNA. The average read length was ∼650 nt (for more details, see Supplemental Table S1).

Nanopore reads are remarkably long compared with other short-read sequencing technologies, and they contain more errors, which need to be corrected (Amarasinghe et al. 2020). One commonly used approach to try to decrease the proportion of errors is to select reads that pass a certain quality score (typically Q ≥ 7). However, we found that eliminating reads with Q < 7 had nearly no effect on the error rate (Fig. 1A), and thus, we did not apply this filter. Instead, we performed a correction based on Illumina reads with the program fmlrc (Wang et al. 2018), taking advantage of a previous S. pombe Illumina RNA-seq experiment performed in the same growth conditions as here (Blevins et al. 2019). The error rate decreased to about half its original values, but some regions, such as the 3′-end of transcripts, still remained largely uncorrected. For this reason, we subsequently applied TranscriptClean, which uses the genome sequence as reference to correct ONT reads (Wyman and Mortazavi 2019). The final “clean” set had an average error rate of only 1.24%, which basically corresponded to short indels.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

dRNA sequencing (dRNA-seq) of S. pombe. (A) Error rate distribution of raw and clean reads. Error rate is the percentage of aligned positions that contain a mismatch or indel. (Raw) The original reads, (raw P7) original reads with quality score Q ≥ 7, and (clean) corrected reads used in this work. (B) Sequence length distribution of raw and clean reads. Number of raw reads, 7,097,130; number of clean reads, 5,054,233; median value raw reads, 500; and median value clean reads, 620. Differences in the distribution are significant by a Wilcoxon test (P-value < 2.2 × 10⁻¹⁶). (C) Correlation transcript abundance ONT dRNA versus Illumina. The reads were mapped to the PomBase transcriptome. In the case of dRNA reads, we simply divided the number of reads by the number of million reads (reads per million [RPM]). For Illumina reads, we calculated transcripts per million (TMP), normalizing by transcript length as well as number of million reads. We selected transcripts expressed in at least one of the two data sets; transcripts with multimapping Illumina reads were removed. Number of transcripts analyzed was 4999. (D) Estimated number of transcripts with at least one full-length read with respect to transcript length. The data are shown for different transcript length deciles: (71.0–633.2), (633.2–923.4), (923.4–1175.0), (1175.0–1395.0), (1395.0–1637.0), (1637.0–1911.2), (1911.2–2227.4), (2227.4–2695.0), (2695.0–3444.6), (3444.6–15,022.0]. Number of transcripts was 6453. (E) dRNA counts with respect to transcript length. Bins are the same as in D. (F) Poly(A) tail distribution in mRNAs and ncRNAs. Poly(A) tail is estimated as the mean of the poly(A) tail length of all the reads that map to each transcript. Differences are significant according to a Wilcoxon test (P-value < 10⁻⁵). (G) Relationship between poly(A) tail length and transcript abundance. For each transcript, the average poly(A) tail length of all the reads mapping to the transcript is taken. Only mRNAs are taken into account for this calculation (n = 4995). Genes related with reproduction (GO:0000003) and translation (GO:0006412) are highlighted. Highly expressed transcripts tend to have shorter poly(A) tails. The correlation is significant (Spearman's ρ = −0.376; P-value = 9.3 × 10⁻¹⁶⁸) (H) Distribution of poly(A) tail length with respect to transcript length. Bins are the same as in D. Poly(A) tail length is homogeneously distributed across different transcript length classes.

The reads were mapped to the PomBase gene annotations with minimap2 (Li 2018). The total number of mapped reads was 5,054,233. The longest mapped read was 13,899 nt long. The mapped reads had an average length of 756 nt and were significantly longer than the raw reads (Fig. 1B). We could see expression of the vast majority of the protein-coding transcripts (97.8%), as well as of a very large percentage of ncRNAs (87.8%) and smaller amounts of other RNA classes (Supplemental Fig. S1). In general, ncRNAs were expressed at much lower levels than mRNAs (average 70 dRNA reads vs. 1130 dRNA reads).

We inspected the correlation between dRNA and Illumina derived transcript abundances. Whereas each dRNA read corresponds to one native molecule, Illumina sequencing involves cDNA synthesis and PCR amplification, and the number of mapped reads needs to be normalized by length. In addition, we found that 13.41% of the Illumina reads were multimapping, increasing the uncertainty in the transcript abundance estimates. There was a high positive correlation between the abundance estimates obtained with the two technologies (Spearman's ρ = 0.849, P < 10⁻¹²) (Fig. 1C), after excluding transcripts with multimapping reads. Inclusion of the 1800 transcripts with Illumina multimapping reads caused an overestimation of transcript expression levels for some of the transcripts (Supplemental Fig. S2).

Nanopore mRNA sequencing starts from the 3′-end of the transcript and proceeds toward the 5′-end. Some of the mRNAs are sequenced to completion (full-length reads), whereas others are truncated at their 5′-end. We estimated the number of full-length reads by mapping the reads to the gene annotations and then comparing the length of each mapped read with the length of the corresponding annotated transcript. To be considered full length, the read had to be equal or longer than the annotated transcript, or in case it was shorter, the difference should be <50 nt. This accounted for the fact that the first 10–15 nt of the 5′ UTR are systematically missed with Nanopore and that the real 3′-end might also show some variation with respect to the annotated transcript. We estimated that the total number of full-length reads was 1,013,789 (20.06%). Perhaps more importantly, the vast majority of the transcripts with expression evidence had at least one full-length read (5165 out of 6453, 80.04%). As expected, the fraction of transcripts recovered as full-length reads decreased with transcript length, with the strongest effect being observed in transcripts >3.45 kb (Fig. 1D, last decile; Supplemental Fig. S3). We observed that this subset of very long transcripts also tended to be expressed at lower levels than transcripts of intermediate length (Fig. 1E; Supplemental Fig. S4). Transcripts in the first decile (<633 nt) were expressed at even lower values, but because of their short size, they were normally recovered as full-length reads.

Poly(A) tail length depends on expression level but not transcript length

Extended poly(A) tail lengths have been previously associated with increased transcript's stability and translatability (Dreyfus and Régnier 2002). We used nanopolish to measure poly(A) tail length directly from the dRNA data. The average poly(A) length was ∼50 nt, similar to that observed in humans (Workman et al. 2019). Poly(A) tails tended to be slightly shorter in mRNAs (median, 48.9 nt) than ncRNAs (median, 51 nt) (Fig. 1F). We found that poly(A) length and transcript abundance were negatively correlated (Spearman's ρ = −0.376 and P-value < 10⁻⁵) (Fig. 1G). The median number of counts for the top 10% transcripts with the shortest poly(A) tail (length <40) was 581.5, whereas the 10% of genes with the longest poly(A) tail (length >58) had a median of 131 counts. Gene Ontology (GO) term enrichment analysis indicated that genes with the shortest poly(A) tail were significantly enriched in translation-related functions, whereas those with the longest poly(A) tail were in sexual reproduction and meiosis-related functions (Supplemental Fig. S5). Consistently, these two groups also showed clear differences in their expression levels, with the translation genes being expressed at very high levels and the meiosis genes at much lower levels (Fig. 1G). No major differences in poly(A) length were observed in relation to transcript length (Fig. 1H).

Identification of hundreds of alternative transcript isoforms

We used StringTie2 to identify possible transcript isoforms supported by the dRNA reads (Kovaka et al. 2019). This program has the advantage that it does not require that the reads are full length, something which a priori cannot be determined for dRNA reads. StringTie2 yielded 5799 transcripts that showed a length distribution similar to that of annotated transcripts (Supplemental Table S2; Supplemental Fig. S6).

We identified a total of 332 alternative isoforms, in 262 different genes, that had an effect on the coding sequence. These events were novel and not annotated in PomBase. Because not all reads corresponded to full-length transcripts, a small proportion of the reads, 3.7%, mapped to different gene isoforms (7271 multimapping reads out of 189,281). The formation of alternative isoforms decreased the relative amount of the reference protein and, in some cases, could potentially lead to different protein products. We could distinguish between four types of events: intron retention (IR), intron inclusion (II), use of alternative splicing (AS) sites, and exon skipping (ES) (Fig. 2A). IR events were denoted by dRNA sequences in which the intron was not spliced out. In the case of II, a nonannotated intron was observed in a subset of the reads. AS sites implied the use of different splice site donor or acceptor signals in a subset of the mRNA molecules. Finally, ES was represented by sequences that lacked a complete exon. The most common event was IR, which represented ∼80% of all events, followed by II in ∼12% of the cases (Fig. 2B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Identification of alternative isoforms using dRNA-seq. (A) Alternative splicing (AS) isoform classes. We built a transcriptome using the dRNA reads with StringTie2. We identified 332 alternative isoforms in 263 different genes. The plot shows one example of each of the four main classes of alternative isoforms detected. Diagrams of the exons of the reference and the alternative isoform as shown, together with the dRNA coverage along the gene. (B) Number of different types of splicing isoforms. Intron retention (IR) represents ∼80% of the events. (C) Number of isoforms per gene. In most cases, only one alternative isoform was detected. The most extreme case corresponds to wtf19, with two annotated isoforms and nine additional alternative isoforms detected here. (D) Relative abundance of reference and alternative isoforms for each gene. Data are for the genes containing at least one alternative isoform. The abundance is computed using the number of mapped dRNA reads; the fraction is then calculated over all isoforms containing mapped reads. Number of reference isoforms is 263 (two for wtf19); number of new alternative isoforms detected here, 332; median fraction reference isoforms, 0.873; and median fraction alternative isoforms, 0.105. (E) Abundance of reference and alternative isoforms. Number of dRNA reads mapped to reference and alternative isoforms. Numbers of isoforms as in D. (F) Skipped exons tend to be smaller than the complete set of exons in the reference annotations. Median length reference exons is 170; median length skipped exons, 63. P-value = 0.00752 Wilcoxon test.

In general, the number of alternative isoforms observed was one or at most two, although in some cases, a larger number of isoforms could be observed (Fig. 2C). The latter cases corresponded to genes of the killer meiotic drive system, a rapidly evolving family of parasitic and antidote genes (Eickbush et al. 2019). The maximum number of alternative isoforms recovered by StringTie2 was nine in wtf19. These isoforms apparently originated from different types of exon/intron inclusion and exclusion events, as well as by the use of alternative splice sites. The capacity of this gene to generate many alternative isoforms might be important for the ongoing arms race that characterizes the gene family.

To quantify the isoform expression levels, we mapped the dRNA reads to the transcriptome and used only uniquely mapped reads. This allowed us to unambiguously distinguish between the reference and alternative isoform transcripts. As expected, alternative isoforms were, in general, found at lower frequencies than the annotated isoform (Fig. 2D). Nevertheless, some nonannotated isoforms were found at very high frequencies, and there was a clear overlap between the expression levels of alternative isoforms and already annotated transcripts (Fig. 2E). As many as 92 alternative isoforms had a frequency >20%; 172 cases, >10%. For example, retention of the first intron in rpl22, a gene encoding 60S ribosomal protein 22, showed a frequency of 30% (12,003 dRNA reads vs. 28,910 for the reference mRNA). In gdt2, coding for a Golgi calcium ion transporter, the IR isoform was supported by 40% of the dRNA reads (565 vs. 858). In the case of elo1, encoding an enzyme involved in fatty acid elongation, an isoform in which the third intron was included represented 42% of all transcripts (201 dRNA reads vs. 276 for the reference mRNA). An extreme case was etp1, a gene involved in the adaptation to high concentrations of ethanol (Snowdon et al. 2009). In this case, the transcript containing the intron was the predominant one (85% of the dRNA reads, 182 out of 212). These examples were further validated by RT-PCR (Supplemental Fig. S7; Supplemental Table S7).

We observed a moderate but significant tendency for the first intron to be retained. In genes with two introns and for which only one of the introns was retained, we found 64 cases in which the first intron was retained and 36 in which the second intron was retained (P-value = 0.007 compared with 50/50, proportion test). We identified nine ES events, mostly affecting very small exons (Fig. 2F). One example was sir2, encoding a histone deacetylase. An isoform in which the fourth exon was skipped represented 17.2% of the transcripts (37 dRNA reads vs. 178 for the reference isoform).

Virtual translation of the sequences of the alternative isoforms indicated that, except in 10 cases, they resulted in proteins that were shorter than the annotated one (Fig. 3A). We sought evidence of protein translation using previously published Ribo-seq data (Duncan and Mata 2017). We focused on IR isoforms, which are the easiest to analyze, because we do not expect to have Ribo-seq reads mapping to the intron except if the intron is retained and translated. In 18 cases, we found a minimum of five Ribo-seq reads supporting the alternative protein. One remarkable example was the translation of an ORF coding for a protein of only 13 amino acids (aa) in the IR isoform of rpl22 (Fig. 3B). The 13-aa protein was supported by 311 isoform-specific Ribo-seq reads. The number of Ribo-seq reads that map to a sequence can be used as a proxy of translation level, because each mapped Ribo-seq read potentially corresponds to a translating ribosome (Ingolia et al. 2009). Examination of the Ribo-seq coverage in the isoform-specific intronic region indicated that, although the 13-aa isoform was translated at levels estimated to be around 1/10 of the canonical 117-aa-long protein, it was still among the top 10% most expressed proteins in the cell.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

Proteome expansion by AS. (A) Difference in size between the putative alternative protein and the reference protein. Data are shown for the four main classes of AS events. In the vast majority of cases, the alternative protein would be smaller than the canonical protein. (B) Alternative 13-aa protein isoform in the rpl22 gene. The diagram shows the putative coding sequence in the alternative IR isoform. A stop codon in frame in the intronic region results in the translation of a 13-aa protein. We obtained Ribo-seq support for the 13-aa alternative protein. We also estimated the Ribo-seq coverage of the Rpl22 canonical protein (SPAC11E3.15.1) and the 13-aa alternative isoform (STRG.1957.2) using isoform-specific coding sequences. The values are compared with those for the coding sequences of all transcripts with five or more Ribo-seq reads mapped to the P-site (n = 5669). The gray area covers 90% of cases. (C) Alternative 38-aa protein isoform in the uap2 gene. IR in uap2 generates a shorter coding sequence, encoding a putative 38 aa protein. (D) Alternative 361-aa protein isoform in the etp1 gene. IR in etp1 results in a protein that is 15 aa longer than the reference one.

Another example of IR supported by ribosome profiling data was uap2, encoding the U2 snRNP-associated protein Uap2. About one-fourth of the transcripts corresponded to retention of the first intron, resulting in a putative protein of 38 aa instead of the standard 367-aa protein (Fig. 3C). Other genes displaying similar patterns were mal3, encoding a microtubule protein (des Georges et al. 2008); rpb4, encoding a RNA polymerase II complex subunit (Sakurai et al. 1999); and not11, encoding a CCR4-NOT complex subunit involved in the shortening of the poly(A) length and initiation of cytoplasmic mRNA decay (Ukleja et al. 2016). In mal3, IR was found at a frequency of 38% and resulted in a putative protein of 26 aa; in rpb4, 26.8% and a protein of 20 aa; and in not11, 60% and a protein of 55 aa. A very different case was etp1; the intron contained no stop codon in frame, and for this reason, the resulting protein was predicted to be 15 aa longer than the reference one (Fig. 3D).

In other genes, different isoforms were generated by the use of AS sites, ES, or II. One example was pat10 (SPAC18B11.08c), a gene that encodes an endoplasmic reticulum protein that is part of a chaperone complex involved in the biogenesis of proteins with multiple transmembrane domains (Chitwood and Hegde 2020). The reference transcript is composed of five exons and encodes a protein that is 95 aa long. The dRNA reads provided direct evidence of an alternative isoform arising from a downstream alternative splice site in intron 3 and skipping of exon 4. The alternative isoform had a frequency of 27% (161 dRNA reads alternative isoform vs. 431 for the reference) and resulted in a putative protein of 74 aa (Supplemental Fig. S8).

IR is associated with extended poly(A) tails

We next investigated poly(A) tail length with respect to AS events. First, for each type of event, we compared the poly(A) length of the dRNA reads in the alternative isoform and the reference isoform. Collectively, we could observe significant differences for IR events and AS site events (Fig. 4A). In both cases, poly(A) tails tended to be longer in the alternative isoform. However, when we compared the differences in poly(A) length for each gene, taking the average value for all reads mapping to the same isoform, a consistent difference was only observed for IR events (Fig. 4B).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Poly(A) length in alternative transcript isoforms. (A) Distribution of poly(A) length by isoform type. We computed poly(A) length for all dRNA reads with nanopolish. Only reads with the label PASS were considered. (n) Number of reads with poly(A) length information; (m) median poly(A) length. (IR) Intron retention; (II) intron inclusion; (AS) alternative splicing site; and (ES) exon skipping. Significant differences were identified for isoform retention and AS site events compared with the corresponding reference isoforms. (***) P-value < 10⁻⁵, Wilcoxon test. (B) Difference in the average poly(A) length between the alternative and reference isoforms. Longer poly(A) lengths were consistently observed for IR isoforms. (C) Negative correlation between average poly(A) length and expression level for reference and alternative transcript isoforms. The data are only for genes in which we detected alternative isoforms. Reference refers to the annotated isoform. Spearman's ρ = −0.41, P = 3.14 × 10⁻²³. (D) Examples poly(A) length differences between the reference and the IR isoforms. We computed the poly(A) length for the dRNA reads that correspond to each of the isoforms. The first example corresponds to RNA polymerase II subunit 4 (rpb4, SPBC337.14), with ∼27% of the reads corresponding to the IR isoform. The second example corresponds to a nucleolar RNA-binding protein also implicated in mRNA processing (tam10, SPBC14C8.19), with ∼18% of the reads mapping to the IR isoform. In both cases poly(A) length showed a significant tendency to be longer in the IR isoform. (***) P-value < 10⁻⁵; (ns) nonsignificant, Wilcoxon test.

In a recent study in S. cerevisiae, the investigators concluded that the poly(A) tail of newly transcribed transcripts is 50 adenosines long on average and is shortened in the cytoplasm to 40 adenosines on average (Tudek et al. 2021). We thus considered the possibility that longer poly(A) tails could be indicative of not fully processed transcripts still retained in the nucleus. However, isoforms with translation evidence, and thus presumably located in the cytoplasm, also had longer poly(A) tails than the reference isoforms (Supplemental Fig. S9). Thus, the data fit quite well the previously observed negative correlation between expression level and poly(A) length (Fig. 4C). Another possible explanation was that only a fraction of the molecules was being translated, whereas the rest was retained in the nucleus, resulting in overall longer poly(A) tails. Because the latter possibility cannot be tested with the current data, the question remains open.

In general, IR isoforms tended to be less abundant than the reference isoform and also tended to have longer poly(A) tails, as shown in the examples in Figure 4D. When we examined cases in which the alternative and reference isoforms had relatively similar abundances, the results varied depending on the gene. In some cases, such as not11, the poly(A) tail length of the alternative and reference isoform was not significantly different. In other cases, including mal3, rps13, and vps38, the differences were significant, although relatively small (Supplemental Fig. S10). In contrast, slm3 and red1 showed very strong and significant differences in poly(A) tail length between the alternative and reference isoforms (approximately 80 vs. 50 nt, respectively) (Supplemental Fig. S11).

Discovery of new transcribed loci

The reconstruction of the transcriptome using the dRNA reads also resulted in the discovery of 214 completely novel transcripts, whose coordinates on the genome did not show any overlap to annotated genes on the same strand. Mapping the Illumina reads to these transcripts confirmed the expression of the majority of them (168 out of 214). We found that about three-fourths of them, 158 (74%), overlapped other genes on the opposite orientation and were classified as antisense. The remaining 56 transcripts were located in regions with no other annotated features and were classified as intergenic.

The novel transcripts tended to be shorter than the annotated ones, especially the intergenic ones (Fig. 5A; Supplemental Table S2). They also tended to be expressed at lower levels than annotated transcripts, although in this case, there were no significant differences between antisense and intergenic transcripts. Because novel transcripts are lowly expressed, their detection might largely depend on the sequencing coverage. To explore this, we generated saturation curves by subsampling the number of original mapped dRNA reads. Whereas the number of known genes that could be detected reached a plateau at approximately 1.5 million reads, the number of novel transcripts showed an approximately linear relationship with the number of sequencing reads (Supplemental Fig. S12). We also investigated the fitness associated with the novel transcripts by using data from a previous saturating transposon mutagenesis experiment (Grech et al. 2019). We found that the level of constraints in the set of novel transcripts showed no significant differences to that of annotated ncRNAs and was clearly weaker than in coding sequences (Supplemental Fig. S13).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

New transcripts and peptides. (A) Transcript length and gene expression for novel antisense and intergenic transcripts. Novel antisense transcripts tend to be longer than intergenic transcripts, and both classes of transcripts are shorter than already characterized ones (Wilcoxon test, P-value < 10⁻⁵). Gene expression levels are lower in novel transcripts compared with not novel ones (Wilcoxon test, P-value < 10⁻⁵). Novel antisense and intergenic transcripts show a similar expression level distribution. (TPM) Transcripts per million, quantification provided by StringTie. (B) Prediction of translated ORFs. The plot shows the RibORF score versus the number of Ribo-seq mapped reads for different classes of transcripts. ORFs with at least 10 mapped reads and a RibORF score higher than 0.7 were selected as translated. The score is based on the Ribo-seq read 3-nt periodicity and homogeneity. Nanopore novel indicates transcripts that did not overlap any annotated transcript, 12 of these transcripts contained ORFs with translation signatures. PomBase noncoding indicates annotated ncRNAs that also had translation signatures. (C) Novel translated ORFs are shorter than annotated ones. Comparison of aa length for ORFs with evidence of translation in novel transcripts and annotated coding sequences. Differences are statistically significant (Wilcoxon test, P-value < 10⁻⁵). (D) New protein identified in S. pombe. The protein labeled as STRG.2094.1 showed significant similarity to other bacteria proteins (BLASTP e-values between 10⁻² and 10⁻⁷) and, in eukaryotes, only to a protein from the fungus Macrophomina phaseolina (e-value <10⁻⁹). The blue box represents the homologous region; lines at the side represent additional protein sequence.

We next used the ribosome profiling data to investigate the translation patterns in the novel transcripts. We predicted putatively translated ORFs using RibORF (Ji et al. 2015). This program produces a score based on 3-nt periodicity and homogeneity of the Ribo-seq reads along the ORF. In previous studies, we established that a RibORF score >0.7 was associated with significant translation activity (Blevins et al. 2021; Moro et al. 2021). As expected, the vast majority of the annotated coding sequences in mRNAs were classified correctly by the program (4514 out of 4560, 98.99%) (Fig. 5B). In addition, 16% of the annotated ncRNAs also contained ORFs with evidence of translation (Supplemental Table S3; Supplemental Fig. S14). These findings are in line with the translation signatures observed in a large fraction of the long ncRNAs in other biological systems (Ruiz-Orera et al. 2014; Ji et al. 2015; Chen et al. 2020).

Among the newly discovered transcripts we identified 12 cases with evidence of translation: eight antisense and four intergenic. One of the intergenic transcripts contained two putatively translated ORFs. The encoded proteins were small, with a median length of 44.5 aa for antisense transcripts and 25 aa for intergenic transcripts compared with 393 aa for canonical ORFs (Fig. 5C). The orientation of four of these novel transcripts, as well as the distance to the nearest transcription start site (<400 bp), suggested divergent transcription from a bidirectional promoter (Supplemental Table S4). One of the newly identified proteins showed significant homology with several prokaryotic proteins as well as to an uncharacterized protein from the fungus Macrophomina phaseolina (Fig. 5D). Given the sparse species distribution, it seems likely that this protein has originated by horizontal gene transfer, probably from bacteria. The rest of genes did not have homology with any other annotated protein or to a set of novel translated ORFs recently discovered in S. cerevisiae (Blevins et al. 2021). Therefore, these genes might have originated de novo in the S. pombe lineage.

Finally, we also used the dRNA-seq data to annotate 5′ and 3′ UTRs, focusing on those that were not yet annotated in PomBase (266 mRNA without a 5′ UTR and 337 mRNAs without a 3′ UTR). Nanopore data are expected to be very accurate for the 3′-end but less so for the 5′-end, as the first 10–15 nt of the mRNA are normally not recovered. Using the dRNA-based transcriptome, we annotated 105 5′ UTRs and 75 3′ UTRs that were previously missing. The median size of these sequences was 217 and 256, respectively, which was comparable to the length of the annotated ones (median 168 for 5′ UTR and 259 for 3′ UTR, respectively). Thus, dRNA provides an effective way to annotate 3′ UTRs and, to some extent, also 5′ UTRs.