Nanopore direct RNA sequencing maps the complexity of Arabidopsis mRNA processing and m6A modification

Abstract

Understanding genome organization and gene regulation requires insight into RNA transcription, processing and modification. We adapted nanopore direct RNA sequencing to examine RNA from a wild-type accession of the model plant Arabidopsis thaliana and a mutant defective in mRNA methylation (m⁶A). Here we show that m⁶A can be mapped in full-length mRNAs transcriptome-wide and reveal the combinatorial diversity of cap-associated transcription start sites, splicing events, poly(A) site choice and poly(A) tail length. Loss of m⁶A from 3' untranslated regions is associated with decreased relative transcript abundance and defective RNA 3' end formation. A functional consequence of disrupted m⁶A is a lengthening of the circadian period. We conclude that nanopore direct RNA sequencing can reveal the complexity of mRNA processing and modification in full-length single molecule reads. These findings can refine Arabidopsis genome annotation. Further, applying this approach to less well-studied species could transform our understanding of what their genomes encode.

Keywords: A. thaliana; RNA; chromosomes; computational biology; epitranscriptome; gene expression; nanopore; polyadenylation; systems biology.

Figure 1.. Diverse Arabidopsis RNAs are detected by nanopore DRS.

(A) Nanopore DRS 12.7 kb read alignment at AT1G48090, comprising 63 exons. Blue, nanopore DRS read alignments; black, Araport11 annotation. (B) Nanopore DRS read alignments at the snoRNA gene U3B (AT5G53902). Blue, sample of nanopore DRS read alignments; black, Araport11 annotation. (C) PIN4 (AT2G01420) antisense RNAs detected using nanopore DRS. Blue, Col-0 PIN4 sense Illumina RNAseq coverage and sample of nanopore PIN4 sense read alignments; orange, Col-0 PIN4 antisense Illumina RNAseq coverage and nanopore PIN4 antisense read alignments; green, hen2-2 mutant PIN4 sense Illumina RNAseq coverage; purple, hen2-2 mutant PIN4 antisense Illumina RNAseq coverage; black, PIN4 sense RNA isoforms found in Araport11; grey, PIN4 antisense differentially expressed regions detected in Illumina RNAseq with DERfinder.

Figure 1—figure supplement 1.. Properties of nanopore DRS sequencing data.

(A) Nanopore DRS identified a 12.8 kb transcript generated from the AT1G67120 gene that includes 58 exons. Blue, nanopore DRS isoform; black, Araport11 annotation. (B) Synthetic ERCC RNA spike-in mixes are detected in a quantitative manner. Absolute concentrations of spike-ins are plotted against counts per million (CPM) reads in log₁₀ scale. Blue, ERCC RNA spike-in mix 1; orange, ERCC RNA spike-in mix 2. (C) Overview of the sequencing and alignment characteristics of nanopore DRS data for ERCC RNA spike-ins. Left, distribution of the length fraction of each sequenced read that aligns to the ERCC RNA spike-in reference; centre, distribution of fraction of identity that matches between the sequence of the read and the ERCC RNA spike-in reference for the aligned portion of each read; right, distributions of the occurrence of insertions (black), substitutions (orange) and deletions (blue) as a proportion of the number of aligned bases in each read. (D) Substitution preference for each nucleotide (left to right: adenine [A], uracil [U], guanine [G], cytosine [C]). When substituted, G is replaced with A in more than 63% of its substitutions, while C is replaced by U in 73%. Conversely, U is rarely replaced with G (12%) and A is rarely substituted with C (16%). (E) Nucleotide representation within the ERCC RNA Spike-In reference sequences (black dots) compared with nucleotide representation within four categories from the nanopore DRS reads. Identity matches between the sequence of the read and the ERCC RNA spike-in reference (green crosses), insertions (blue pentagons), deletions (yellow stars) and substitutions (purple diamonds).G is under-represented and U is over-represented for all three categories of error (insertion, deletion and substitution) relative to the reference nucleotide distribution. C is over-represented in deletions and substitutions. A is over-represented in insertions and deletions and under-represented in substitutions. (F) Signals originating from the RH3 transcripts are susceptible to systematic over-splitting around exons 7–9 (highlighted using a purple dashed box), resulting in reads with apparently novel 5′ or 3′ positions. This appears only to occur at high frequency in datasets collected after May 2018 (Supplementary file 1) and may result from an update to the MinKNOW software. (G) PIN7 antisense RNAs detected using nanopore DRS. Blue, Col-0 PIN7 sense Illumina RNAseq coverage and nanopore PIN7 sense read alignments; orange, Col-0 PIN7 antisense Illumina RNAseq coverage and nanopore PIN7 antisense read alignments; green, hen2–two mutant PIN7 sense Illumina RNAseq coverage; purple, hen2–two mutant PIN7 antisense Illumina RNAseq coverage; black, PIN7 sense RNA isoforms found in Araport11 annotation; grey, PIN7 antisense differentially expressed regions detected with DERfinder.

Figure 2.. Nanopore DRS reveals poly(A) tail length and maps 3′ cleavage and polyadenylation sites.

(A) Normalized histogram showing poly(A) tail length of RNAs encoded by different genomes in Arabidopsis. Blue, nuclear (n = 2,348,869 reads); orange, mitochondrial (n = 2,490 reads); green, chloroplast (n = 1,848 reads). (B) Distance between the RNA 3′ end positions in nanopore DRS read alignments and the nearest polyadenylation sites identified by Helicos data. (C) Nanopore DRS identified 3′ polyadenylation sites in RNAs transcribed from FPA (AT2G43410). Blue track, coverage of nanopore DRS read alignments; blue, collapsed read alignments representing putative transcript annotations detected by nanopore DRS; black, Araport11 annotation; pPAS, proximal polyadenylation site; dPAS, distal polyadenylation sites.

Figure 2—figure supplement 1.. 3′ end processing is revealed by nanopore DRS.

(A) Poly(A) tail length estimates for ERCC spike-in controls. Boxplots showing distribution of poly(A) length estimates for ERCC spike-in controls with more than 100 mapped reads for which tail length could be successfully estimated. Expected poly(A) tail lengths are shown as orange points. (B) The RNA poly(A) tail length negatively correlates with the gene expression level. Expression in log₂ scale of counts per million (CPM) obtained from nanopore DRS data is plotted against the median poly(A) tail length. ρ, Spearman’s correlation coefficient; black line, locally weighted scatterplot smoothing (LOWESS) regression fit. (C) Nanopore DRS identified known 3′ polyadenylation sites in RNAs transcribed from the IBM1 (AT3G07610) locus. Blue track, nanopore DRS coverage; blue, isoforms detected by nanopore DRS; black, Araport11 annotation; PAS, proximal polyadenylation site. (D) Nanopore DRS identified novel 3’ polyadenylation sites in RNAs transcribed from the PTM (AT5G35210) locus. Black (on the top), Helicos DRS 3’ coverage; blue track, nanopore DRS coverage; blue, isoforms detected by nanopore DRS; black (on the bottom), Araport11 annotation; green, five annotated transmembrane domain regions from Uniprot entry F4JYC8 (PTM_ARATH) that mapped to exons; pPAS, proximal polyadenylation site; dPAS, distal polyadenylation sites.

Figure 3.. Cap-dependent ligation of an adapter enables detection of authentic RNA 5′ ends.

(A) 5′ adapter RNA ligation reduces 3′ bias in nanopore DRS data at RCA (AT2G39730) from 0.92 to 0.03. Blue line, exonic read coverage at RCA for reads without (left) and with adapter (right); orange line, median coverage; orange shaded area, interquartile range (IQR). Change in 3′ bias can be measured using the IQR/median = quartile coefficient of variation (QCoV). (B) 5′ adapter RNA ligation reduces 3′ bias in nanopore DRS data. Histogram showing the QCoV in per base coverage for each gene in the Araport11 annotation, for reads with a 5′ adapter RNA (orange), and reads without a 5′ adapter RNA (blue). (C) Cap-dependent adapter ligation allows identification of authentic 5′ ends using nanopore DRS. The cumulative distribution function shows the distance to the nearest Transcription Start Site (TSS) identified from full-length transcripts cloned as part of the RIKEN RAFL project (left), or 5’ tag identified from nanoPARE data (right), for reads with a 5′ adapter RNA (orange) compared with reads without a 5′ adapter RNA (blue). (D) Cap-dependent adapter ligation enabled resolution of an additional 11 nt of sequence at the RNA 5′ end. Histogram showing the nucleotide shift in the largest peak of 5′ coverage for each gene in data obtained using protocols with vs without a 5′ adapter RNA ligation. (E) For RCA (AT2G39730), the 5′ end identified using cap-dependent 5′ adapter RNA ligation protocol was consistent with Illumina RNAseq and full-length cDNA start site data but differed from the 5′ ends in the Araport11 and AtRTD2 annotations. Grey, Illumina RNAseq coverage; blue, nanopore DRS 5′ end coverage generated without a cap-dependent ligation protocol; green/pink, nanopore DRS 5′ end coverage for read alignments generated using the cap-dependent ligation protocol with (green) and without (pink) 5′ adapter RNA; orange, 5’ coverage of nanoPARE data; blue, TSSs identified from full-length transcripts cloned as part of the RIKEN RAFL project; black, Araport11 and AtRTD2 annotations (with duplicated 5’ positions removed).

Figure 3—figure supplement 1.. Nanopore DRS with cap-dependent ligation of 5′ adapter RNA.

(A) Histogram showing the distribution of 5′ adapter RNA length in the nanopore raw current signal, as inferred from alignment of the mRNA sequence to the signal using nanopolish eventalign. The median signal length was 1441 points and 96% of adapter signals were 3000 points or less. (B) Out-of-bag receiver operator characteristic curve showing the performance of the trained convolutional neural network at detecting 5′ adapter RNA using 3000 points of signal. The curve was generated using five-fold cross validation. (C) Out-of-bag precision recall curve showing the performance of trained neural network, generated using five-fold cross validation. (D) Alternative transcription start sites were identified using nanopore DRS with cap-dependent ligation of a 5′ end adapter at the AT1G17050 and AT5G18650 genes. Orange, 5’ coverage for capped nanoPARE reads; blue track, nanopore DRS coverage with cap-dependent ligation of 5′ adapter RNA; blue, isoforms detected by nanopore DRS with cap-dependent ligation of 5′ adapter RNA; black, Araport11 annotation. (E) Reads mapping to ERCC RNA spike-ins lack approximately 11 nt of sequence at the 5′ end. Histogram showing the distance to the 5′ end for ERCC RNA spike-in reads (each spike-in is shown in a different colour; only those with > 1000 supporting reads are shown). (F) Reads mapping to in vitro transcribed mGFP lack approximately 11 nt of sequence at the 5′ end. Histogram showing the distance to the 5′ end for in vitro transcribed mGFP. (G) Araport11 annotation overestimates the length of 5′ UTRs. The cumulative distribution function shows the distance to the nearest TSS identified from full-length transcripts cloned as part of the RIKEN RAFL project (blue) and Araport11 annotation (orange). (H) Nanopore DRS detects miR170/miR171 cleavage products of Hairy Meristem 1 (HAM1, AT2G45160) transcripts. Orange, 5’ coverage from capped nanoPARE reads; purple, 5’ coverage from uncapped nanoPARE reads; blue, nanopore DRS 5’ coverage; grey, miRNA target site alignment is shown; black, Araport11 annotation.

Figure 4.. Nanopore DRS reveals the complexity of alternative splicing.

(A) Nanopore DRS identified the mutually exclusive alternative splicing of the second and third exons of FLOWERING LOCUS M (FLM, AT1G77080) and several novel isoforms. Black arrows indicate mutually exclusive exons. Orange, isoforms present in the AtRTD2 annotation but not identified using nanopore DRS; blue, isoforms common to both AtRTD2 and nanopore DRS; green, novel isoforms identified in nanopore DRS. (B) Comparison of splicing events identified in error-corrected and non-error-corrected nanopore DRS, Illumina RNA sequencing, and Araport11 and AtRTD2 annotations. Bar size represents the number of unique splicing events common to the set intersection highlighted using circles (see Supplementary file 5 for the exact values). GU/AG splicing events are shown on the top and non-GU/AG on the bottom of the plot: yellow, splicing events common to all five datasets; green, events common to both error-corrected and non-error-corrected nanopore DRS with support in orthogonal datasets; blue, events common to both nanopore DRS datasets without orthogonal support; orange, events found in uncorrected nanopore DRS (but not error corrected) with orthogonal support; pink, events found in error-corrected nanopore DRS (but not uncorrected) with orthogonal support. (C) Comparison of RNA isoforms (defined as sets of co-spliced introns) identified in error-corrected full-length nanopore DRS, Araport11 and AtRTD2 annotations. Bar size represents the number of splicing events common to a group highlighted using circles below (see Supplementary file 5 for the exact values): yellow, unique splicing patterns found in nanopore DRS and both reference annotations; blue, novel isoforms found only in nanopore DRS.

Figure 4—figure supplement 1.. Patterns of splicing revealed using nanopore DRS.

(A) Nanopore DRS novel splicing events can be validated using RT-PCR. Blue, RNA isoforms found using nanopore DRS; orange, Sanger sequencing products obtained using RT-PCR; black, RNA Araport11 annotation. (B) Splice junction classification of unannotated GU/AG splice sites found in error-corrected nanopore DRS data that also have Illumina support. Counts are plotted in log₁₀ scale and the exact number of splice junctions in each set is indicated. (C) A novel exon skipping event at exon 3 of KH domain-containing protein AT5G56140, identified by nanopore DRS, leads to a frameshift resulting in a premature termination codon in exon 4 (orange asterisk). Blue track, nanopore DRS cleavage, blue, nanopore DRS read alignments; black, AtRTD2 annotation; orange, full novel open reading frame (ORF). (D) A selection of novel alternative acceptor site at XRCC (AT1G80420), identified in nanopore DRS, leads to a frameshift resulting in a premature termination codon in exon 6 (orange asterisk). Blue track, nanopore DRS cleavage, blue, nanopore DRS read alignments; black, AtRTD2 annotation; orange, full novel open reading frame (ORF).

Figure 5.. Differential error rate analysis identifies sites of VIR-dependent m⁶A modifications transcriptome-wide.

(A) Loss of VIR function is associated with reduced error rate in nanopore DRS. Histogram showing the log₂ fold change in the ratio of errors to reference matches at bases with a significant change in error profile in vir-1 mutant compared with the VIR-complemented line. Orange and green shaded regions indicate sites with increased and reduced errors in vir-1, respectively. (B) The motif at error rate sites matches the consensus m⁶A target sequence. The sequence logo is for the motif enriched at sites with reduced error rate in the vir-1 mutant. (C) Differential error rate sites are primarily found in 3′ UTRs. Bar plot showing the number of differential error rate sites per kb for different genic feature types of 48,149 protein coding transcript loci in the nuclear genome of the Araport11 reference. Upstream and downstream regions were defined as 200 nt regions upstream of the annotated transcription start sites and downstream of the annotated transcription termination sites, respectively. (D) Differential error rate sites and miCLIP peaks are similarly distributed within the 3′ UTR, without accumulation at the stop codon. Metagene plot centred on stop codons from 48,149 protein coding transcript loci, showing the frequency of nanopore DRS error sites (blue) and miCLIP peaks (orange). (E) The locations of differential error rate sites are in good agreement with the locations of sites identified by miCLIP. Histogram showing the distribution of distances to the nearest miCLIP peak for each site of reduced error. Most error sites (66%) are within 5 nt of a miCLIP peak. (F) Nanopore DRS differential error site analysis and miCLIP identify m⁶A sites in the 3′ UTR of PRPL34 RNA. Blue, miclip 5′ end coverage; purple, nanopore DRS differential error sites; black, RNA isoform from Araport11 annotation. Expanded region: blue, miCLIP coverage; purple, nanopore DRS differential error sites; black, motif scores using the m⁶A consensus position weight matrix (Figure 5B). A higher positive score denotes a higher likelihood of a match to the consensus sequence.

Figure 5—figure supplement 1.. Identification of VIR-dependent m⁶A transcriptome-wide.

(A) Gene track showing the alignment of nanopore DRS reads from unmodified synthetic MALAT1 RNAs (MALAT1 control) and MALAT1 RNAs with m⁶A modification (MALAT1 m⁶A). Alignments are coloured by the basecall at each position (A = green, C = blue, G = yellow, U = orange). Grey, coverage tracks (insertions to the reference are not included in coverage counts); black, test statistic derived from the G-test of aligned basecall profiles between the two alignments; square brackets, modified GGm⁶ACU k-mer. (B) vir-1 shows reduced levels of m⁶A compared with Col-0 and restored m⁶A levels in the VIR-complemented line (VIRc). The ratio of m⁶A/A obtained using LC-MS. Blue, Col-0; orange, vir-1; green, VIRc. (C) Frequency of m⁶A motifs detected at vir-1 reduced error sites, as detected by FIMO using the motif detected de novo by MEME and an FDR threshold of 0.1. (D) Bar plot showing the number of miCLIP peaks per kb of different genic feature types in the Araport11 reference. Upstream and downstream regions were defined as 200 nt regions upstream of the annotated transcription start sites and downstream of the annotated transcription termination sites, respectively. (E) Subsampling analysis of MALAT1. Different coloured lines show detection rate at varying subsample sizes for unmethylated sample vs. mixed unmethylated and methylated sample. Fractions in legend denote the fraction of reads in mixed sample which originate from methylated data. (F) Subsampling analysis of a selection of four methylation sites in three genes, in Arabidopsis. Different coloured lines show detection rate at varying subsample sizes for each methylation site.

Figure 6.. Reduction in m⁶A RNA modification leads to disruption of the circadian clock.

(A) Genes with differential error rate sites have lower detectable RNA levels. Histogram showing the log₂ fold change in protein coding gene expression based on counts from nanopore DRS reads in the vir-1 mutant compared to the VIR-complemented line. Blue, genes with differential error rate sites (n = 5083 genes); orange, genes without differential error rate sites (n = 14,675 genes). (B) The circadian period is lengthened in the vir-1 mutant. Mean delayed fluorescence measurements in arbitrary units: blue, Col-0 (n = 61 technical reps); green, the VIR-complemented line (VIRc; n = 29 technical reps); orange, vir-1 (n = 24 technical reps). Shaded areas show bootstrapped 95% confidence intervals for the mean. (C) Boxplot showing the period lengths calculated from delayed fluorescence measurements shown in (b): blue, Col-0; green, VIRc; orange, vir-1 (orange). (D) Poly(A) tail length is altered in the vir-1 mutant. Histogram showing the poly(A) tail length distribution of CAB1 (AT1G29930): blue, Col-0 (n = 40,841 reads); green, VIRc (n = 65,810 reads); orange, vir-1 (n = 68,068 reads). Arrows indicate phased peaks of poly(A) length at approximately 20, 40 and 60 nt. vir-1 distribution is significantly different from both Col-0 and VIRc, according to the Kolmogorov–Smirnov test (p < 1 × 10⁻¹⁶, p < 1 × 10⁻¹⁶ respectively).

Figure 6—figure supplement 1.. Changes in the gene expression of circadian clock components in the vir-1 mutant.

(A) Histogram showing log₂ fold changes in gene expression based on Illumina RNAseq data for the vir-1 mutant and VIRc. Blue, genes with differential error rate sites (n = 5,095 genes); orange, genes without differential error rate sites (n = 14,686 genes). (B) Expression of core circadian clock components is perturbed in the vir-1 mutant. Boxplots showing normalized gene expression measured using Illumina RNAseq in log₂counts per million (CPM): green, the VIR-complemented line (VIRc); orange, the vir-1 mutant. Asterisks, significant expression changes (using an FDR threshold of 0.05); orange labelled genes, genes with 3′ UTR m⁶A detectable by nanopore DRS and miCLIP. Each scatter point represents a single biological replicate. UBIQUITIN LIGASE 21 (UBC21; AT5G25760) was used as a control. (C) Expression of CCA1, encoding a regulator of the circadian rhythm in Arabidopsis, is increased in the vir-1 mutant. Boxplot showing the gene expression change from Col-0 (measured by RT-qPCR) for VIRc (green) and vir-1 (orange). Three technical replicates of three biological replicates were conducted. Each scatter point represents the comparison of a technical replicate of treatment (VIRc or vir-1) against control (Col-0). The expression change in vir-1 is significant (p = 0.02; marked by asterisk). UBIQUITIN LIGASE 21 (UBC21; AT5G25760) was used as a control. (D) Expression of the Flowering Locus C (FLC) gene is decreased in the vir-1 mutant. Boxplot showing gene expression change from Col-0 (measured by RT-qPCR) for VIRc (green) and vir-1 (orange). Three technical replicates of three biological replicates were conducted. Each scatter point represents the comparison of a technical replicate of treatment (VIRc or vir-1) against control (Col-0). Expression change in vir-1 is significant (p = 2.4 × 10⁻¹⁴; marked by asterisk). The following source data are available for Figure 6—figure supplement 1.

Figure 7.. RNA methylation is closely linked to poly(A) site selection.

(A) Global shift to proximal 3′ end usage observed in the vir-1 mutant compared with the VIRc line. Histogram showing the distance in nucleotides between the most reduced and most increased 3′ end positions for genes in which the 3′ end profile is altered in vir-1 (detected with the Kolmogorov–Smirnov test, FDR < 0.05). A threshold of 13 nt was used to detect changes in the 3′ end position. (B) A shift to the use of more proximal poly(A) sites is observed in m⁶A containing transcripts in the vir-1 mutant. Histogram showing distance from the error site (n = 17,491 error sites) to upstream and downstream 3′ ends: orange, vir-1; green, VIRc. (C) Poly(A) site choice is phased with m⁶A at PRPL34 (AT1G29070). PRPL34 has four distinct poly(A) site clusters identified from VIRc reads. Analysis of the error profiles of these clusters at the closest upstream m⁶A motif indicates that proximally polyadenylated isoforms are less likely to be methylated. Green, VIRc 3’ coverage and poly(A) clusters; purple, nanopore differential error site tested; black, Araport11 annotation. Zoomed view shows stacked barplots of VIRc error profiles for the tested nanopore differential error site. Pairs of clusters with significantly different profiles (G-test, FDR < 0.05) are indicated with asterisks. (D) PRPL34, which is methylated in at least two positions in the 3’UTR (Figure 5F), also displays an increase in alternative polyadenylation at a proximal poly(A) site in the vir-1 mutant. The proximal poly(A) site is approximately 17 nt downstream of the closest methylated position. Purple, nanopore differential error sites; green, VIRc coverage; orange, vir-1 coverage; black, Araport11 annotation. (E) RNA methylation sites overlap with canonical and degenerate poly(A) site motifs. Line plot showing the enrichment of poly(A) site motifs (with an edit distance of one) around m⁶A motifs detected de novo under nanopore differential error sites: blue, the canonical PAS motif AAUAAA; grey, other similar motifs (excluding AAAAAA); orange, the UAUAAA motif; green, the AACAAA motif; yellow, the AAGAAA motif; black, the m⁶A motif (includes all sequences shown in Figure 5—figure supplement 1E). (F) Readthough events and chimeric RNAs are detected in vir-1. Green, nanopore DRS and Illumina RNAseq data for the VIRc line; orange, nanopore DRS and Illumina RNAseq data for the vir-1 mutant; grey and white, differentially expressed regions between vir-1 and VIRc detected using Illumina RNAseq data with DERfinder that are upregulated and downregulated, respectively; black, Araport11 annotation. Intergenic readthrough regions are highlighted by the purple dashed rectangle. (G) ALG3-GH3.9 chimeric RNAs transcribed in the vir-1 mutant. Chimeric readthrough RNAs were validated using RT-PCR. RT-PCR products were separated on agarose gel and stained with SYBR Safe. RT+ and RT-, cDNA produced in the presence and absence of reverse transcriptase, respectively; ladder, DNA size ladder; Water, no template control. TUBULIN was used as a control.

Figure 7—figure supplement 1.. Changes in RNA 3′ end formation in the vir-1 mutant.

(A) Splicing is moderately disrupted in the vir-1 mutant. Results of differential exon usage analysis with DEXseq are shown for contiguous regions (‘exon chunks’), which occur in the same sets of transcripts in the Araport11 reference. Regions were classified as a 5′ or 3′ variation if they were bounded by the TSS of one or more transcripts. Orange, features with increased usage; blue, features with reduced usage. Boxplots show the distribution in absolute log₂ fold change for each feature set. (B) DaPars analysis of Illumina RNAseq data indicates global 3’UTR shortening in the vir-1 mutant. Histogram showing the change in the percentage of distal poly(A) site usage in vir-1 compared to VIRc. Only genes which are significantly differentially polyadenylated (FDR < 0.05) are shown. There is a clear directional effect towards proximal poly(A) site use. (C) vir-1 mutant exhibits increased readthrough of an intronic proximal poly(A) site in intron 5 of the Symplekin domain encoding gene TANG1 (AT1G27595). Green track, nanopore DRS coverage from VIRc; green, nanopore DRS isoforms from VIRc; orange track, nanopore DRS reads from the vir-1 mutant; orange, nanopore DRS isoforms from vir-1; black, Araport11 annotation; dashed black box highlights the site of proximal polyadenylation.

Author response image 1.. Low amounts of input RNA increase the number of Enolase II reads and short noise reads.

Read length distribution of mapped (blue) and unmapped (orange) reads originating from a standard DRS experiment (left panel) and an adapter ligation experiment (right panel).

Author response image 2.. Noise reads have low quality scores.

Boxplot showing distribution of median per-read quality scores for mapped and unmapped reads from an adapter ligation DRS experiment.

Author response image 3.. Differential error rate analysis of S. cerevisiae SK1 WT vs ime4.

(a) Histogram showing the log₂ fold change in the ratio of errors to reference matches at bases with a significant change in error profile in ime4 mutant compared with wild type S. cerevisiae SK1. Orange and green shading indicates sites with increased and reduced errors in ime4, respectively. (b) Sequence logo showing the motif enriched at sites with reduced error rate in ime4. (c) Frequency of m⁶A motifs detected at ime4 reduced error sites, as detected by FIMO using the motif detected de novode novo by MEME and an FDR threshold of 0.1.

Author response image 4.. Smoothed distribution of m⁶A sites detected by different approaches around stop codons in S.cerevisiae SK1.

Author response image 5.. Antisense Illumina RNA-seq reads at RCA are likely to be spurious.

(a) Scatter plot showing correlation of sense and antisense read abundances derived from Illumina RNAseq data mapping to ERCC spike-in controls, shown in blue. Each point represents a single spike in transcript detected in a single RNAseq replicate (both Col-0 and hen2-2). Since ERCC spike-ins should not have antisense RNAs, this indicates the level of spurious antisense signal in the dataset. RCA, shown in the orange, does not display antisense RNA levels greater than the ERCC spike-ins, indicating that there are not likely to be genuine antisense RNAs at this locus. FLC, which is known to have antisense RNAs, is shown in green as a positive control. (b) Strip plot showing the sense/antisense read ratio at RCA does not change in the hen2-2 mutant (T test p=0.88). Many genuine antisense RNAs, e.g. those at FLC, are unstable and are therefore upregulated in the hen2-2 mutant relative to sense transcripts.

Author response image 6.. Comparison of differential error rate analysis in MALAT1 data using older (3.1.2) and newer (3.2.4) Guppy releases.

Gene track showing the change in the per- reference base test statistic attained when comparing methylated and unmethylated synthetic MALAT1 fragments basecalled using Guppy 3.1.2, which uses recurrent neural network basecalling, and Guppy 3.2.4, which uses recurrent neural network basecalling with a flip-flop architecture.

Author response image 7.. Prevalence of concatemerised reads in nanopore DRS data.

Barplot showing the percentage of read alignments in nanopore DRS data where softclipped regions (i.e. unmapped sequence at 5’ and 3’ ends of alignments) which is mappable elsewhere in the TAIR10 genome, possibly indicating undersplitting of read signal by the MinKNOW software. Adapter ligated datasets do not show increased levels, suggesting that RNA ligase I does not cause a detectable increase in the number of concatemers/chimeric reads.

Author response image 8.. Nanopore DRS captures downstream products of miRNA cleavage at transcripts encoded by TAS genes.

Orange, 5’ coverage from capped nanoPARE reads; purple, 5’ coverage from uncapped nanoPARE reads; blue, nanopore DRS 5’ coverage; black miRNA target site alignment; black, Araport11 annotation.