A ligation-independent sequencing method reveals tRNA-derived RNAs with blocked 3' termini

Abstract

Despite the numerous sequencing methods available, the diversity in RNA size and chemical modification makes it difficult to capture all RNAs in a cell. We developed a method that combines quasi-random priming with template switching to construct sequencing libraries from RNA molecules of any length and with any type of 3' modifications, allowing for the sequencing of virtually all RNA species. Our ligation-independent detection of all types of RNA (LIDAR) is a simple, effective tool to identify and quantify all classes of coding and non-coding RNAs. With LIDAR, we comprehensively characterized the transcriptomes of mouse embryonic stem cells, neural progenitor cells, mouse tissues, and sperm. LIDAR detected a much larger variety of tRNA-derived RNAs (tDRs) compared with traditional ligation-dependent sequencing methods and uncovered tDRs with blocked 3' ends that had previously escaped detection. Therefore, LIDAR can capture all RNAs in a sample and uncover RNA species with potential regulatory functions.

Keywords: RNA-seq; ligation-independent; non-coding RNA; rTNA fragment; sequencing; small RNA; sperm RNA; tDR; tRF; template switch.

Figure 1.. A modified Smart-seq3 protocol captures small RNAs

(A) Schematic of LIDAR. The four key modifications compared to Smart-seq3 are indicated with numbers highlighted in yellow. (B) Gel electrophoresis of pre-amplification (step III) libraries constructed from total ESC RNA using the Smart-seq3 or LIDAR TSO, and random (R) or quasi-random (QR) RT primers. The black line on the side indicates productive libraries; the white arrowheads indicate adapter dimers. (C) Gel electrophoresis of pre-amplification (step III) libraries constructed from decreasing amounts of a synthetic 20 nts RNA. The black arrowhead indicates libraries with the 20 bp insert; the white arrowheads indicate adapter dimers. (D) Read coverage of a 20 nts RNA molecule cloned with LIDAR (blue) or a ligation-based protocol (gray). Read density is expressed as % of maximum coverage in each method. (E) Number of reads per million (RPM) sequenced mapping to 20 or 50 nts synthetic RNAs cloned by LIDAR or ligation. (F) Average (n = 3) biotype distribution of ncRNAs as % of mapped reads in LIDAR, ligation, or Takara SMARTer libraries from total ESC RNA. (G) Genome browser snapshot of average (n = 3) LIDAR and ligation-based read coverage (counts per million, CPM) on two miRNAs. (H) Average (n = 3) size distribution of reads mapping to miRNAs in LIDAR or ligation-based libraries from total ESC RNA. See also Fig. S1 and S2.

Figure 2.. LIDAR detects differential abundance of small and large RNAs

(A) Schematic of the differentiation protocol. (B) Transcri[ptome changes between NPCs and ESCs measured by mRNA-seq (left) or LIDAR (right) on total RNA. Dark grey, genes with significant changes (adjusted p < 0.05, n = 3). Red, NPC markers. Blue, ESC markers. (C) Z-score-converted RNA abundance, calculated as transcripts per million (TPM), of differentially expressed protein-coding genes between NPCs (N) and ESCs (E), as measured by LIDAR or mRNA-seq. Individual replicates are shown. The average rolling mean of the log₂(fold-change) is shown on the right. (D) Differentially expressed protein-coding genes detected in LIDAR (blue) or mRNA-seq (grey). The p value is from a hypergeometric test. (E) Correlation of average (n = 3) log₂(fold-changes) of differentially expressed protein-coding genes detected by LIDAR and mRNA-seq. r_s, Spearman’s rank correlation coefficient. (F–H) Same as in (C–E) but for miRNAs. See also Fig. S2.

Figure 3.. LIDAR uncovers long 3’ tDRs not efficiently captured by ligation

(A) Types of tDRs. (B) Average (n = 3) tDR distribution as % of reads mapping to any tDR in LIDAR and ligation-based libraries from total, < 200 nts, and < 50 nts ESC RNA. (C) Average (n = 3, ± SEM) read end (left) or beginning (right) position frequency, expressed as % of all reads mapping to the corresponding tDR type, in LIDAR (blue) or ligation-based (grey) libraries starting from < 200 nt or < 50 nts ESC RNA. (D) Left, genome browser snapshots showing reads mapping to two different tRNAs in LIDAR or ligation-based libraries starting from < 200 nt ESC RNA. Position of various loops, as in (A), are indicated on top. The 3’-terminal CCA is depicted as a blue box. Right, northern blots probing for the corresponding tDR (black arrow). The two lanes are from two independent biological replicates. See also Fig. S3.

Figure 4.. Tissue-specific 3’ tDR expression detected by LIDAR

(A) Z-score-converted abundance (TPM) of miRNAs differentially expressed in brain (B) vs. testis (T), as measured by LIDAR or ligation-based sequencing. Individual replicates are shown. (B) Overlap of differentially expressed miRNAs in brain and testis detected by LIDAR or ligation. The p value is from a hypergeometric test. (C) Reads mapping to piRNAs in brain or testis expressed as % of aligned reads. Bars indicate mean ± SEM (n = 3). (D) Average (n = 3) % of reads mapping to different classes of tDRs from brain and testis RNA. (E) Differential abundance of 3’ tDRs in brain vs. testis according to LIDAR (n = 3). Differentially expressed 3’ tDRs (adjusted p < 0.05) are in dark gray. Among these, the 3’ tDRs with 10-fold higher levels in LIDAR compared to ligation-based libraries are in blue. The 3’ tDRs validated by northern blot are highlighted. (F) Average (n = 3) read density at two example tRNA loci in brain or testis. (G) Left: northern blot probing for Arg-TCG-1-1 3’ tDR (black arrow) in brain and testis. Individual replicates are shown. Right: density plot showing size distribution (x axis) and abundance (CPM, y axis) of reads for Arg-TCG-1-1 tRNA in brain and testis. (H) Same as (G) but for Gly-ACC-1-1 3’ tDR. See also Fig. S4.

Figure 5.. Base modifications on 3’ tDRs are revealed by LIDAR

(A–B) Average (n = 3) misincorporation rate (as % detected) for every canonical position (column) in every long 3’ tDR (A) or short 3’-tRF (B) iso-acceptor (rows) in LIDAR or ligation-based libraries starting from < 200 nts ESC RNA. No coverage is shown in gray. Position of tRNA loops (as in Fig. 3A) are on top. (C) Average (n = 3) % of reads mapping to different tDRs from < 50 nts ESC (grown in serum-free medium) RNA in LIDAR or PANDORA-seq libraries. (D) Detected 5’ tDRs (top) and 3’ tDRs (bottom) in LIDAR or PANDORA-seq libraries used in (C). (E) As in (A–B) but for all 3’ tDRs comparing LIDAR with PANDORA-seq. PANDORA-seq data was obtained from Shi et al.. See also Fig. S5.

Figure 6.. LIDAR detects full-length tRNAs and tDRs with blocked 3’ ends

(A) Top: schematic of RNAs with possible blocked 3’ ends (cP: 2’−3’ cyclic phosphate; A.A.: amino acid; ?: unknown). Bottom: only RNAs (red) with free 3’ ends can be cloned via ligation, whereas LIDAR captures RNA with free or blocked 3’ ends. (B) Ratio of reads mapping to 3’ biotinylated (blocked) vs. 3’OH synthetic 20 nts RNA in LIDAR and ligation libraries. (C) Average (n = 3) size distribution of all reads mapping to tRNAs as % of all tRNA reads, in LIDAR or ligation-based libraries from ESC RNA < 200 nts. (D) Average (n = 3) coverage and representation of ESC tRNA anticodons (color-coded) in LIDAR or ligation libraries from ESC RNA < 200 nts. v, collapsed variable loop position. (E) Scheme for the enrichment of RNA with blocked 3’ ends. Input RNAs with 3’OH, 3’P, and 2’−3’cP (cP) were end-repaired with T4 PNK, polyadenylated with polyA polymerase (PAP), and depleted with oligo-dT beads (“B”). RNA with blocked 3’ end cannot be polyadenylated and remain in the flow-through (“FT”). (F) Urea PAGE of total and < 200 nts ESC RNA enriched for blocked transcripts (“FT”) as in (E). 100% input (In) and polyA bound (“B”) are shown. (G) Z-score-converted abundance (TPM) of differentially expressed (adjusted p < 0.05) snoRNAs, miRNAs, and tRNAs between LIDAR libraries from total or “blocked” ESC RNA. Z-score-converted abundance in ligation-based libraries from unfractionated (all) ESC RNA is also shown. (H) Differential abundance, calculated over all reads mapping to tDRs, of 3’-tRF (left) and long 3’ tDRs (right) isoacceptors between LIDAR libraries from all < 200 nts RNA (n = 3) vs. blocked (n = 2) RNA from ESCs. Comparisons with adjusted p value < 0.05 are colored according to their anticodon. (I) Average (n > 2) read density of LIDAR from all or blocked ESC RNA < 200 nts and ligation-based libraries on two tRNAs. See also Fig. S6.

Figure 7.. Sperm contains full-length tRNAs and long 3’ tDRs

(A) Average (n = 2) biotype distribution of ncRNAs, as % of mapped reads, in LIDAR from sperm RNA before or after enriching for blocked RNA and in ligation-based libraries. (B) Average (n = 2) size distribution of reads mapping to cytosolic or mitochondrial tRNAs, expressed as % of reads of all tRNA reads. Libraries obtained with LIDAR from all or blocked sperm RNA are shown, as well as ligation-based libraries. (C) Average (n = 2) tDR distribution, expressed as % of reads mapping to any tDR, in LIDAR from all sperm RNA, blocked RNA, or ligation-based libraries. (D) Average (n = 2) read density of LIDAR (from all or blocked sperm RNA) or ligation-based libraries on three example tRNAs. (E) Differential 3’ tDR (any length) abundance in LIDAR libraries on blocked (n = 2) vs. all (n = 3) RNA from sperm. Isoacceptors for which the comparison has an adjusted p value < 0.05 are colored according to their anticodon. See also Fig. S7 and S8.