High-resolution quantitative profiling of tRNA abundance and modification status in eukaryotes by mim-tRNAseq

Abstract

Measurements of cellular tRNA abundance are hampered by pervasive blocks to cDNA synthesis at modified nucleosides and the extensive similarity among tRNA genes. We overcome these limitations with modification-induced misincorporation tRNA sequencing (mim-tRNAseq), which combines a workflow for full-length cDNA library construction from endogenously modified tRNA with a comprehensive and user-friendly computational analysis toolkit. Our method accurately captures tRNA abundance and modification status in yeast, fly, and human cells and is applicable to any organism with a known genome. We applied mim-tRNAseq to discover a dramatic heterogeneity of tRNA isodecoder pools among diverse human cell lines and a surprising interdependence of modifications at distinct sites within the same tRNA transcript.

Keywords: RNA modification; tRNA quantitation; transfer RNA; translation regulation.

Graphical abstract

Figure 1

An optimized workflow for full-length cDNA library construction from eukaryotic tRNA pools (A) Schematic of template-switching TGIRT reactions primed by an RNA/DNA duplex with a single-nucleotide 3′ overhang and a gel image of cDNA products from endogenously modified tRNA pools from S. cerevisiae (Sc), K562 cells (Hs), or a synthetic unmodified tRNA (Syn) at different reaction temperatures and salt concentration. Red, reaction conditions previously used for tRNA library construction; asterisks, premature stops to cDNA synthesis; hash, potential products from end-to-end linkage of tRNAs. (B) Schematic of the mim-tRNAseq library generation workflow. Top gel image: 3′ adapter ligation reactions with four barcoded adapters. Ligation efficiency was measured by normalizing input tRNA band intensity to that in reactions from which Rnl2trKQ was omitted. Bottom gel image: comparison of cDNA yield in short (1 h) or extended (16 h) primer-dependent TGIRT RT on a mix of adapter-ligated tRNA pools from S. cerevisiae and human K562 and HEK293T cells. See also Figure S1 and STAR methods.

Figure 2

The mim-tRNAseq computational pipeline: a comprehensive framework for tRNA sequencing data analysis (A) Bowtie and Bowtie 2 alignment strategies and mapping statistics for a tRNA library from HEK293T cells constructed with the mim-tRNAseq workflow (n = 1). (B) Outline of the mim-tRNAseq computational pipeline. (C) Alignment statistics of HEK293T data (as in A; n = 1) using the mim-tRNAseq pipeline. (D) Uniquely aligned read proportions for inosine 34 (I34)- and uridine 34 (U34)-containing Ser and Pro tRNA isoacceptors using the three alignment strategies on a HEK293T dataset generated as in Figure 1B. (E) Distribution of uniquely aligned reads among tRNA isotypes in published datasets and mim-tRNAseq from HEK293-derived cell lines (hydro-tRNAseq and QuantM-tRNAseq: HEK293 T-Rex Flp-IN; DM-tRNAseq control or AlkB-treated [+AlkB] and mim-tRNAseq library construction: HEK293T). Proportions were obtained from published counts per tRNA (“publ”) or after re-analysis of the datasets with the mim-tRNAseq pipeline (“new”). tRNA families that carry the same amino acid (isotypes) are sorted by the number of RT barriers annotated in MODOMICS (decreasing from top to bottom; grayscale, isotypes without MODOMICS annotation). See also Figure S1 and STAR methods.

Figure 3

mim-tRNAseq improves quantitative analysis of tRNA pools in cells from diverse eukaryotes (A) Alignment statistics for mim-tRNAseq datasets from the indicated cell types. Bars and labels indicate average values, dots show individual sample values (n = 2). (B) Metagene analysis of scaled sequence coverage across nuclear-encoded tRNA isotypes ordered per sample by differences between 3′ and 5′ coverage (decreasing order from top to bottom; n = 1). y axis values normalized to the second-to-last bin from the 3′ end. Each x axis bin represents 4% of tRNA length. Indicated are major known barriers to RT. (C) Boxplots of full-length fraction per tRNA transcript in datasets from (B) (center line and label: median; box limits: upper and lower quartiles; whiskers: 1.5 × interquartile range). (D) Correlation plots of unique tRNA gene copy number and corresponding proportion of uniquely aligned tRNA reads in single replicates (same samples as in B) from S. cerevisiae, S. pombe, and D. melanogaster BG3-c2 cells and hiPSCs. Solid blue lines: linear regression model; shaded gray: 95% confidence interval (CI). See also Figures S2 and S3.

Figure 4

mim-tRNAseq accurately captures differential tRNA expression and aminoacylation with single-transcript resolution (A) Differential expression analysis of unique tRNA transcripts in HEK293T and K562 relative to hiPSCs. Axes represent log-transformed normalized read counts from DESeq2, with significant down- and upregulation in hiPSCs indicated with closed orange and green triangles, respectively (false discovery rate [FDR]-adjusted one-sided Wald test p ≤ 0.01, n = 2). (B) Differential expression analysis as in (A) for counts per tRNA anticodon family. (C) Left panel: hierarchically clustered expression heatmap showing scaled z score of normalized unique transcript counts in HEK293T, K562, and hiPSCs (n = 2). Middle panels: differential expression for HEK293T and K562 relative to iPSCs (values, log₂ fold changes; bar plots, numbers of up- and downregulated genes in green and orange, respectively). Right panel: base mean normalized per tRNA transcript across all samples. (D) Northern blot analysis of tRNA-Arg-UCU-4 and tRNA-Gly-CCC-2 in HEK293T, K562, and hiPSCs (n = 2, matched samples to those used for mim-tRNAseq). Band intensities were quantified by densitometry and normalized to the mean value for HEK293T. (E) Relative abundance of tRNA-Arg-UCU-4 and tRNA-Gly-CCC-2 in HEK293T, K562, and hiPSCs measured by mim-tRNAseq (C) or northern blotting (D), normalized to the mean value for HEK293T (n = 2, matched samples). (F) tRNA charging analysis in wild-type and trm7Δ S. cerevisiae. Charged tRNA are represented by proportion of reads with 3′-CCA ends (light green, in %). Light green bars and tRNA-Phe-GAA labels, average charged tRNA fractions (% CCA; n = 3). See also Table S1 and Figure S3.

Figure 5

Near-complete modification readthrough in mim-tRNAseq datasets enables modification discovery and annotation (A) Average proportion of stops (red) and misincorporation rates (blue) per nucleotide for all tRNA unique transcripts (n = 2) in S. cerevisiae, S. pombe, and D. melanogaster BG3-c2 cells and hiPSCs. x axis, canonical tRNA position at major sites with known RT barriers. (B) RT readthrough per annotated modification aggregated for cytosolic and mitochondrial tRNA from the four species. (C) Boxplots of misincorporation signatures for annotated modified sites as in (B) (center line: median; box limits: upper and lower quartiles; whiskers: 1.5 × interquartile range). Signatures stratified by upstream context (rows) and modification type (columns); proportion per nucleotide scaled to total misincorporation at this site. (D) Boxplot of misincorporation signature at G37 of tRNA-Phe-GAA from WT and trm7Δ S. cerevisiae (n = 3). (E) Modified site discovery by mim-tRNAseq (“new”) compared with misincorporation-inducing modified sites previously annotated in MODOMICS (“annot.”). Labels indicate percentage of newly detected sites relative to annotated ones. See also Figure S4 and Table S2.

Figure 6

Misincorporation rates in mim-tRNAseq reflect modification stoichiometry (A) Global heatmap of average misincorporation proportions in S. cerevisiae per unique tRNA transcript with coverage above 2,000 reads (n = 2; top bar graph, mean misincorporation per position; right bar graph, number of sites per transcript with detectable misincorporation signatures in ≥10% of reads spanning that position). (B) Relative misincorporation proportions at G9 in samples from wild-type (WT) S. cerevisiae and trm10Δ (lacking m¹G9) or mixes thereof (filtered for clusters with ≥10% misincorporation in WT and scaled to WT proportion; solid blue line, linear regression model; shaded gray, 95% CI). (C) Analysis as in (B) but for misincorporation at G26 in samples from WT S. cerevisiae or trm1Δ (lacking m²₂G26). (D) Misincorporation proportions per canonical nucleotide position and identity (aggregated per species; e₂, second nucleotide of variable loop; n = 2). (E) Significant changes in misincorporation rates in trm1Δ relative to WT S. cerevisiae (FDR-adjusted chi-square p ≤ 0.01, log₂ fold change ≥ 0.5, n = 1). See also Figures S5 and S6 and Table S3.