Quantitative tRNA-sequencing uncovers metazoan tissue-specific tRNA regulation

DOI: https://doi.org/10.21203/rs.3.pex-1031/v1


Transfer RNAs (tRNA) are quintessential in deciphering the genetic code; disseminating nucleic acid triplets into correct amino acid identity. While this decoding function is clear, an emerging theme is that tRNA abundance and functionality can powerfully impact protein production rate, folding, activity, and messenger RNA stability.  Importantly, however, the expression pattern of tRNAs is obliquely known. Here we present Quantitative Mature tRNA sequencing (QuantM-tRNA seq), a technique to monitor tRNA abundance and sequence variants secondary to RNA modifications. With QuantM-tRNA seq we assess the tRNA transcriptome in mammalian tissues. We observe dramatic distinctions in isodecoder expression and known tRNA modifications between tissues. Remarkably, despite dramatic changes in tRNA isodecoder gene expression, the overall anticodon pool of each tRNA family is similar across tissues. These findings suggest that while anticodon pools appear to be buffered via an unknown mechanism, underlying transcriptomic and epitranscriptomic differences suggest a more complex tRNA regulatory landscape.



HEK293 T-Rex Flp-IN cells (ThermoFisher. cat# R78007)

Dulbecco’s modified essential media (Thermo Fisher)

fetal bovine serum (Thermo Fisher)

 penicillin and streptomycin (Thermo Fisher)

trypsin (Thermo Fisher)

Trizol (Thermo Fisher)

 Q5 2X master mix (NEB)

QiaQuick gel extraction kit (Qiagen)

HiScribe T7 High Yield RNA synthesis kit (NEB)

T4 polynucleotide kinase (NEB)

rtStarTM tRNA-optimized First-Strand cDNA Synthesis Kit, ArrayStar)

nanodrop spectrophotometer (ThermoFisher)

RNA ligase 2 (NEB) 

glycoblue (ThermoFisher)

Superscript IV (ThermoFisher)

SYBR gold (ThermoFisher)

CircLigase (Epicentre)

NEBnext Ultra II Q5 next-generation master mix (NEB)


Qubit Fluorometer (ThermoFisher)

Bioanalyzer (Agilent)

DNA HS bioanalyzer chip (Agilent)

NextSeq 550 (v2.5)


Cell culture and RNA isolation

1. HEK293 T-Rex Flp-IN cells (ThermoFisher. cat# R78007) were cultured at 37°C with 5% CO2 in complete Dulbecco’s modified essential media (Thermo Fisher) supplemented with 10% fetal bovine serum (Thermo Fisher) and 1% penicillin and streptomycin (Thermo Fisher).

2. All passaging was performed with trypsin (Thermo Fisher) according to manufacturer’s suggested conditions. Following passaging, cells were plated at 25% confluency in 10 cm tissue culture treated dishes and cultured for 48-72 hours until ~90% confluent.

3. At 90% confluency, media was aspirated and 1mL of ice cold Trizol (Thermo Fisher) was added to the culture and mixed on ice for 30 seconds to assure a homogeneous solution.

4. Samples were stored in Trizol at -80°C until further processing. RNA was isolated according to manufacturer’s protocol with two 75% ethanol washes following isopropanol precipitation.

5. Samples were resuspended in distilled H2O and stored at -80°C until library preparation.


In vitro transcribed tRNA spike-in preparation

1. Five mature tRNA sequences derived from E. coli tRNA (gtRNAdb v8) were purchased as gBlocks (Integrated DNA Technologies; Supplementary Table 1) with a Hepatitis delta virus (HDV) ribozyme sequence at the 3’ end of the sequence to generate a precise CCA 3’ end.

2. gBlocks were amplified using the Q5 2X master mix (NEB) and the T7 F and HDV R primers in Supplementary Table 1 according to manufacturer’s suggested conditions. Design of the transcripts was in accordance with previously published protocols36.

3. Amplification products of the appropriate length were purified from a native agarose gels stained with 0.05% EtBr using the QiaQuick gel extraction kit (Qiagen).

4. Up to 1 µg of double stranded template DNA was added to the HiScribe T7 High Yield RNA synthesis kit (NEB) and transcribed according to manufacturer’s suggested conditions. The HDV cleavage reaction occurs spontaneously in the reaction conditions required for in vitro transcription.

5. The cleaved tRNA product was then purified from a denaturing polyacrylamide gel using the crush and soak method.

6. Removal of the 2’, 3’-cyclo-phosphate group on the 3’ end of the purified tRNA product was performed by T4 polynucleotide kinase (NEB).

7. The repaired tRNA product was quantified using a Qubit fluorometer and individual tRNA species were mixed to cover approximately 3.5 orders of magnitude. Appropriate mixing of the spike-in mix was assessed with a Qubit Fluorometer (ThermoFisher) and Bioanalyzer (Agilent). 17 ng of control tRNAs were spiked into 1 µg of total RNA from HEK293 total RNA treated with demethylase.


QuantM-tRNA seq library preparation

1. Total RNA samples were quantified using a nanodrop spectrophotometer (ThermoFisher) prior to library preparation and RNA integrity was checked on a 1.2% denaturing formaldehyde agarose gel.

2. To remove 3’ conjugated amino acids, total RNA was deacylated in deacylation buffer (final concentration 20mM Tris-HCL pH=9.0) at a final concentration of 1 µg/µL.

3. Where indicated, deacylated total RNA was treated with demethylase (rtStarTM tRNA-optimized First-Strand cDNA Synthesis Kit, ArrayStar) and cleaned up per the manufacturer’s instructions.

4. 1 µg of deacylated total RNA from each sample was subject to library preparation. 10 pmol of the 3’ and 10 pmol of the 5’ single-stranded adapter mix (2.5 pmol of each adapter 5’-TGrGrA-3’, 5’-TGrGrT-3’, 5’-TGrGrG-3’, 5’-TGrGrC-3’; Supplementary Table 1) were added to a 200 µL thin-walled amplification tube and denatured at 95°C for 2 minutes.

5. Then annealing buffer was added to a final concentration of 5 mM Tris-HCl (pH 8.0), 0.5mM ethylenediaminetetraacetic acid (EDTA), and 10 mM MgCl2 and incubated at 37°C for 15 minutes to hybridize the annealed double-stranded adapter to tRNA.

6. The ligation reaction was catalyzed by 5 U/ µL of RNA ligase 2 (NEB) with manufacturer’s suggested conditions at 37°C for 60 minutes then 4°C at 60 minutes.

7. All reactions were ethanol precipitated with glycoblue (ThermoFisher) followed by two 75% ethanol washes, then suspended in 10 µL of dH2O.

8. Following ligation, synthesis of cDNA began with hybridization of the reverse transcriptase primer to the ligated total RNA with a final concentration of 0.5 pmol/µL (10 pmol total).

9. The samples were incubated at 70°C for 2 minutes and temperature was reduced to 37°C by 0.1°C/sec. Synthesis of cDNA was achieved using Superscript IV at 55°C for 60 minutes.

10. To remove DNA-RNA dimers following cDNA synthesis, RNA was hydrolyzed with a final concentration of 0.1 N NaOH in dH2O at 98°C for 20 minutes.

11. All reactions were ethanol precipitated with glycoblue (ThermoFisher) followed by two 75% ethanol washes, then suspended in 12 µL of dH2O.

12. cDNA libraries were separated using 7M urea 6% denaturing polyacrylamide gels. Gels were stained with 1X SYBR gold (ThermoFisher) in 1X TBE for 15 minutes and regions representing tRNA derived cDNAs were excised on a UV light box.

13. Gel slices were sheared through the bottom of a 0.5 mL tube nested in a 1.7 mL tube by centrifugation then suspended in 400 µL of DNA elution buffer (300 mM NaCl, 10 mM Tris (pH=8.0, 1 mM EDTA), incubated on dry ice for 30 minutes, and allowed to incubate at room temperature overnight on a standing rotator.

14. cDNA was isopropanol precipitated with glycoblue followed by two 75% ethanol washes then was resuspended in 12 µL of dH2O.

15. Circularization of cDNA libraries was performed with CircLigase (Epicentre) at 0.5U/µL using manufacturer’s suggested conditions at 60°C for 1 hour. The reaction was terminated with incubation at 80°C for 20 minutes.

16. All reactions were ethanol precipitated with glycoblue followed by two 75% ethanol washes, then suspended in 12.5 µL of dH2O.

17. cDNA libraries were amplified using the NEBnext Ultra II Q5 next-generation master mix (NEB) with manufacturer’s suggested conditions. HEK293 libraries were amplified for 7 cycles.

18. Amplified libraries were gel purified from 2% agarose gels stained with 0.05 mg/ml ethidium bromide. Regions of interest (100-250bp) were excised on a UV light box and purified using the Qiaquick gel extraction kit (Qiagen) taking care to dissolve gel slices at room temperature and using all optional steps.

19. All libraries were ethanol precipitated with glycoblue (thermo) with two 75% ethanol washes, then suspended in 10 µL of dH2O before submitting for sequencing.

20. Library concentration was assessed using a Qubit (ThermoFisher), quality was assessed on a DNA HS bioanalyzer chip (Agilent), and library multiplexing directed by qPCR.

21. Sequencing was performed as single-end reads for 110 cycles on a NextSeq 550 (v2.5).


Read quality control and alignment

         Reads were first processed to remove 5’ adapter sequences using cutadapt --cut 2 then cutadapt -g TCCAACTGGATACTGGN -e 0.2 followed by cutadapt –a CCAGTATCCAGTTGGAATT -e 0.2 to remove 3’ CCA and adapter sequences. Custom human or mouse tRNA references were generated by collapsing identical tRNA sequences from gtRNAdb Release 18 hg38 or mm10 high confidence mature tRNA fasta files. Mapping of human or mouse reads with the corresponding reference was done with bowtie2 using the parameters: --quiet --min-score G,1,8 --local -D 20 -R 3 -N 1 -L 10 -i S,1,0.5. Isodecoder-level read count tables for further analyses were produced by counting reads with MAPQ > 10 over reference tRNAs using the Rsubread package’s featureCounts function in R. Anticodon-level read count tables were then created by summing reads from all isodecoders with the same anticodon. In addition to raw read count tables, tables of both isodecoder and anticodon-level reads per million mapped read (RPM) values were generated by dividing raw read counts * 1,000,000 by the number of reads mapped.


Differential expression analysis

         The raw read count tables at both the anticodon-level and isodecoder-level across all 7 mouse tissues described in the previous section were next used to perform differential tRNA expression analysis. The likelihood ratio test was applied to these tables using DESeq2 in R as detailed in https://hbctraining.github.io/DGE_workshop/lessons/08_DGE_LRT.html (Command: DESeq(raw_count_table, test="LRT", reduced = ~ 1) using default settings and p-value adjustment (Benjamin-Hochberg correction). Downstream data visualization and plotting were performed using ggplot2, gplots (heatmap.2), ggrepel, and ggforce in custom R scripts.


Variant analysis

         In order to analyze variants in tRNA sequencing reads, a custom Python script was used to generate variant counts at each position in every tRNA across all 7 mouse tissues. In brief, bam files were read into the script and the CIGAR string and MD tags for each read were used to tabulate each mutation, insertion, or deletion across every ribonucleotide base of all tRNA in the mouse reference. In addition, 5’ ends of reads internal to tRNA were used to infer sites of RT stalling or fall-off. These 4 types of variants were summed at each position of each tRNA for a total variant count, and then a read coverage at each position was also calculated.

         To identify significantly changed sequencing variants across tissues, we performed DEXseq analysis on the raw variant counts table in R. DEXseq was originally devised to identify alternative processing events in mRNA, but we reasoned that co-transcriptional splicing is similar in principle to RT-mediated misincorporation/stalling at RNA modifications. To ensure robust detection of variants that change across tissues, we added two additional filtering steps. First, for a given tRNA base, we required that variant percentage be > 1% on average in every tissue. Next, we only accepted base-level variants that changed variant percentage at least 1.5-fold across tissues. Downstream data visualization and plotting were performed using ggplot2, gplots, ggrepel, and ggforce in custom R scripts as well as matplotlib in custom python scripts.


Time Taken

In vitro transcription and mixing of e. coli spike-in control

6 days

Library preparation

4-10 days depending upon number of libraries

Data analysis

variable (3-14 days depending upon needs and hardware)

Anticipated Results

Please see 'procedure' for anticipated results.


1. Nirenberg, M.W. & Matthaei, J.H. The dependence of cell-free protein synthesis in E. coli upon naturally occurring or synthetic polyribonucleotides. Proc Natl Acad Sci U S A 47, 1588-602 (1961).

2. Phizicky, E.M. & Hopper, A.K. tRNA biology charges to the front. Genes Dev 24, 1832-60 (2010).

3. Schaffer, A.E., Pinkard, O. & Coller, J.M. tRNA Metabolism and Neurodevelopmental Disorders. Annu Rev Genomics Hum Genet 20, 359-387 (2019).

4. Pang, Y.L., Poruri, K. & Martinis, S.A. tRNA synthetase: tRNA aminoacylation and beyond. Wiley Interdiscip Rev RNA 5, 461-80 (2014).

5. Krutyholowa, R., Zakrzewski, K. & Glatt, S. Charging the code - tRNA modification complexes. Curr Opin Struct Biol 55, 138-146 (2019).

6. Pan, T. Modifications and functional genomics of human transfer RNA. Cell Res 28, 395-404 (2018).

7. Chan, P.P. & Lowe, T.M. GtRNAdb 2.0: an expanded database of transfer RNA genes identified in complete and draft genomes. Nucleic Acids Res 44, D184-9 (2016).

8. Agris, P.F. et al. Celebrating wobble decoding: Half a century and still much is new. RNA Biol 15, 537-553 (2018).

9. Crick, F.H. Codon--anticodon pairing: the wobble hypothesis. J Mol Biol 19, 548-55 (1966).

10. Bornelov, S., Selmi, T., Flad, S., Dietmann, S. & Frye, M. Codon usage optimization in pluripotent embryonic stem cells. Genome Biol 20, 119 (2019).

11. Graczyk, D., Ciesla, M. & Boguta, M. Regulation of tRNA synthesis by the general transcription factors of RNA polymerase III - TFIIIB and TFIIIC, and by the MAF1 protein. Biochim Biophys Acta Gene Regul Mech 1861, 320-329 (2018).

12. Geslain, R. & Pan, T. Functional analysis of human tRNA isodecoders. J Mol Biol 396, 821-31 (2010).

13. Hanson, G. & Coller, J. Codon optimality, bias and usage in translation and mRNA decay. Nat Rev Mol Cell Biol 19, 20-30 (2018).

14. Gingold, H. et al. A dual program for translation regulation in cellular proliferation and differentiation. Cell 158, 1281-1292 (2014).

15. Dittmar, K.A., Goodenbour, J.M. & Pan, T. Tissue-specific differences in human transfer RNA expression. PLoS Genet 2, e221 (2006).

16. Abbott, J.A., Francklyn, C.S. & Robey-Bond, S.M. Transfer RNA and human disease. Front Genet 5, 158 (2014).

17. Goodarzi, H. et al. Modulated Expression of Specific tRNAs Drives Gene Expression and Cancer Progression. Cell 165, 1416-1427 (2016).

18. Santos, M., Fidalgo, A., Varanda, A.S., Oliveira, C. & Santos, M.A.S. tRNA Deregulation and Its Consequences in Cancer. Trends Mol Med 25, 853-865 (2019).

19. Zhang, Z. et al. Global analysis of tRNA and translation factor expression reveals a dynamic landscape of translational regulation in human cancers. Commun Biol 1, 234 (2018).

20. Kapur, M., Monaghan, C.E. & Ackerman, S.L. Regulation of mRNA Translation in Neurons-A Matter of Life and Death. Neuron 96, 616-637 (2017).

21. Fujishima, K. & Kanai, A. tRNA gene diversity in the three domains of life. Front Genet 5, 142 (2014).

22. Goodenbour, J.M. & Pan, T. Diversity of tRNA genes in eukaryotes. Nucleic Acids Res 34, 6137-46 (2006).

23. Zheng, G. et al. Efficient and quantitative high-throughput tRNA sequencing. Nat Methods 12, 835-837 (2015).

24. Gogakos, T. et al. Characterizing Expression and Processing of Precursor and Mature Human tRNAs by Hydro-tRNAseq and PAR-CLIP. Cell Rep 20, 1463-1475 (2017).

25. Shigematsu, M. et al. YAMAT-seq: an efficient method for high-throughput sequencing of mature transfer RNAs. Nucleic Acids Res 45, e70 (2017).

26. Xu, H., Yao, J., Wu, D.C. & Lambowitz, A.M. Improved TGIRT-seq methods for comprehensive transcriptome profiling with decreased adapter dimer formation and bias correction. Sci Rep 9, 7953 (2019).

27. Kurschat, W.C., Muller, J., Wombacher, R. & Helm, M. Optimizing splinted ligation of highly structured small RNAs. RNA 11, 1909-14 (2005).

28. Kietrys, A.M., Velema, W.A. & Kool, E.T. Fingerprints of Modified RNA Bases from Deep Sequencing Profiles. J Am Chem Soc 139, 17074-17081 (2017).

29. Potapov, V. et al. Base modifications affecting RNA polymerase and reverse transcriptase fidelity. Nucleic Acids Res 46, 5753-5763 (2018).

30. Chou, C.C., Chen, C.H., Lee, T.T. & Peck, K. Optimization of probe length and the number of probes per gene for optimal microarray analysis of gene expression. Nucleic Acids Res 32, e99 (2004).

31. Dittmar, K.A., Mobley, E.M., Radek, A.J. & Pan, T. Exploring the regulation of tRNA distribution on the genomic scale. J Mol Biol 337, 31-47 (2004).

32. Liu, H., Bebu, I. & Li, X. Microarray probes and probe sets. Front Biosci (Elite Ed) 2, 325-38 (2010).

33. Ishimura, R. et al. RNA function. Ribosome stalling induced by mutation of a CNS-specific tRNA causes neurodegeneration. Science 345, 455-9 (2014).

34. Liu, F. et al. ALKBH1-Mediated tRNA Demethylation Regulates Translation. Cell 167, 1897 (2016).

35. Clark, W.C., Evans, M.E., Dominissini, D., Zheng, G. & Pan, T. tRNA base methylation identification and quantification via high-throughput sequencing. RNA 22, 1771-1784 (2016).

36. Schurer, H., Lang, K., Schuster, J. & Morl, M. A universal method to produce in vitro transcripts with homogeneous 3' ends. Nucleic Acids Res 30, e56 (2002).

37. Lee, C.P., Mandal, N., Dyson, M.R. & RajBhandary, U.L. The discriminator base influences tRNA structure at the end of the acceptor stem and possibly its interaction with proteins. Proc Natl Acad Sci U S A 90, 7149-52 (1993).

38. Limmer, S., Hofmann, H.P., Ott, G. & Sprinzl, M. The 3'-terminal end (NCCA) of tRNA determines the structure and stability of the aminoacyl acceptor stem. Proc Natl Acad Sci U S A 90, 6199-202 (1993).

39. Wende, S., Bonin, S., Gotze, O., Betat, H. & Morl, M. The identity of the discriminator base has an impact on CCA addition. Nucleic Acids Res 43, 5617-29 (2015).


We thank members of the Coller lab, Daniel Fischer, and Drs. Harvey Lodish, Peter Eimon, and Hsin Chen for important scientific discussions related to this work. We also thank Dr. Ashleigh Schaffer for helpful discussion. This research was supported by the Genomics Core Facility of the CWRU School of Medicine's Genetics and Genome Sciences Department.