TRMT6/61A-dependent base methylation of tRNA-derived fragments regulates gene-silencing activity and the unfolded protein response in bladder cancer

Abstract

RNA modifications are important regulatory elements of RNA functions. However, most genome-wide mapping of RNA modifications has focused on messenger RNAs and transfer RNAs, but such datasets have been lacking for small RNAs. Here we mapped N¹-methyladenosine (m¹A) in the cellular small RNA space. Benchmarked with synthetic m¹A RNAs, our workflow identified specific groups of m¹A-containing small RNAs, which are otherwise disproportionally under-represented. In particular, 22-nucleotides long 3' tRNA-fragments are highly enriched for TRMT6/61A-dependent m¹A located within the seed region. TRMT6/61A-dependent m¹A negatively affects gene silencing by tRF-3s. In urothelial carcinoma of the bladder, where TRMT6/61A is over-expressed, higher m¹A modification on tRFs is detected, correlated with a dysregulation of tRF targetome. Lastly, TRMT6/61A regulates tRF-3 targets involved in unfolded protein response. Together, our results reveal a mechanism of regulating gene expression via base modification of small RNA.

Fig. 1. Systematic mapping of m¹A sites in small RNA space.

a Chemical structure of m¹A (N¹-methyladenosine). b Overall m¹A mapping strategy for small RNAs by combining m¹A RIP and m¹A-induced mismatch analysis. Synthetic m¹A RNA is included as a positive control. c TGIRT identified as the optimal reverse transcriptase for m¹A-induced mismatch analysis. Briefly, synthetic m¹A-containing RNA at different m¹A stoichiometry was sequenced by three different reverse transcriptases (RTs): TGIRT (thermostable group II intron reverse transcriptase), PSII (ProtoScriptII, retrovirus RT) and RT-1306 (engineered HIV RT). For each RT, mismatch rate was calculated across all the reads that map to the synthetic RNA sequence (allowing 1 mismatch) as represented by the sequence logo. The mismatch rate at the known m¹A site is then plotted against known m¹A stoichiometry (R² derived from linear fitting forced to cross intercept at zero). d ProtoScriptII leads to under-representation of m¹A-containing RNAs. Cloning frequency is normalized to spike-ins. e m¹A RIP successfully enriches synthetic m¹A-containing RNA compared to input, when spiked in with total short RNAs. TGIRT captures enrichment better than ProtoScriptII. Data are based on four independent RIP experiments (two HEK293T and two U251). Boxplot center represents median, bounds represent 25 and 75%, and whiskers show the minimum or maximum no further than 1.5 * interquartile range from the bound. See also Supplementary Fig. 1 and Supplementary Table 1. Source data are provided as a Source Data file.

Fig. 2. Specific tRNA-derived fragments are highly enriched for m¹A.

a m¹A RNA immunoprecipitation (RIP) successfully enriches for tRNAs compared to IgG control RIP, by both Sybr gold staining and Northern blot against tRNA^Ala. b m¹A RIP enriches for reads that contain A->C/G/T mismatch in tRNA fragments. Bar graph represents RPM (reads per million) values for each group of small RNAs allowing no mismatch (black) and additional mapped reads when allowing 1-nt mismatch (red). Labeled percentage shows the percentage of A->C/G/T mismatch among all mismatches. c Major types of tRFs include tRF-3s (3′ end of mature tRNAs), tRF-5s (5′ end of mature tRNAs) and tRF-1s (trailer of precursor tRNAs). tRF-3s have two major isoforms, tRF-3b (22-nt) and tRF-3a (18-nt). d, e tRF-3s are significantly enriched by m¹A RIP. d Log2 fold enrichment (m¹A RIP/input) was plotted against small RNA abundance for different RNA groups. Only significantly (p.adj < 10⁻⁵ by DESeq2 Wald test with multiple hypothesis adjustment) enriched or depleted RNAs were colored. e tRF-3s (n = 114) and tRF-5s (n = 73), but not tRF-1s (n = 34) are significantly more enriched by m¹A RIP than miRs (n = 421) (p value by two-sided Wilcoxon test compared to miRs). Boxplot center represents median, bounds represent 25 and 75%, and whiskers show the minimum or maximum no further than 1.5 * interquartile range from the bound. f Bar graph shows 22-nt tRF-3b but not 18-nt tRF-3a are significantly enriched by m¹A RIP compared to input. Data in this figure are based on TGIRT-seq of two independent RIP experiments in HEK293T. See also Supplementary Fig. 2. Source data are provided as a Source Data file.

Fig. 3. m¹A on 22-nucleotides 3′ tRNA fragments is dependent on TRMT6/61A.

a Scheme of 3′ tRNA fragments (tRF-3s) that are derived from 3′ end of mature tRNAs. tRF-3b indicates 22-nt tRF-3s, and tRF-3a indicates 18-nt tRF-3s. b–e m¹A mismatch on 22-nt tRF-3bs are globally regulated by TRMT6/61A. b, c Mismatch rate is calculated for the A4 position on tRF-3bs, based on TGIRT-seq of two independent knock-down experiments in HEK293T. d Knock-down efficiency is confirmed by western blot and RT-qPCR. RT-qPCR data was represented as mean ± SD (p value based on two-tailed paired student’s t test from three independent knock-down experiments). e Coverage plot of all tRF reads mapped on the parental tRNA^Ala shows decrease of m¹A of tRF-3b but not tRF-3 abundance (TGIRT-seq in HEK293T, one replicate shown as example). f m¹A mismatch on A4 position of tRF-3b is also decreased by siTRMT6/61A in HeLa and U251, based on TGIRT-seq. g Northern blot confirms detection of both tRF-3a and tRF-3b isoforms, which are not altered by siTRMT6/61A in HEK293T. Full-length tRNA^Ala, tRNA^Tyr and U6 are probed as a loading control. See also Supplementary Fig. 3. Source data are provided as a Source Data file.

Fig. 4. m¹A attenuates gene silencing by tRF-3s.

a m¹A is located within the seed region of tRF-3b, which represents the guide RNA in RISC (RNA-induced silencing complex). m¹A may interfere with base pairing between tRF-3b and the target RNA, thus alleviate target repression. b Scheme of the dual luciferase reporter assay to measure tRF-3 gene-silencing activity. tRF-3b complementary target sequence is inserted in the 3′ UTR of Renilla luciferase (Rluc) gene. Firefly luciferase (Fluc) signal is used for normalization. Synthetic tRF-3b with unmodified A₄ or m¹A-modified A₄ is co-transfected with the reporter plasmid to measure relative effect of tRF-3b on the reporter activity. c–f Relative tRF-3 gene-silencing activity is measured by the dual luciferase reporter assay in HEK293T. Relative activity is calculated from Rluc signal divided by Fluc signal, and normalized to the empty site reporter; NT1 (non-targeting control1) is set as 1 in each replicate. g–j Relative tRF-3 gene-silencing activity is de-repressed after tRF-3 knock-down by LNAs (lock nucleic acids). Relative activity is calculated from Rluc signal divided by Fluc signal and normalized to the empty site reporter. Non-targeting LNA was used as control (Ctrl). All data are represented as mean ± SD; p value based on two-tailed paired student’s t test from three (four for d, e) independent experiments. Source data are provided as a Source Data file.

Fig. 5. Argonaute association of tRF-3b is not decreased by TRMT6/61A-dependent m¹A.

a Construction of stable FH-Ago2 cells, detected by anti-pan-Ago western blot (*indicates non-specific band detected by the antibody). The pan-Ago band detects both endogenous Ago proteins and tagged FH-Ago2 protein. b Top 100 ranking small RNAs (gray dots: microRNAs, red dots: tRFs) detected in Ago2 RIP (RNA immunoprecipitation) are plotted (Y-axis: read counts on log10 scale). c Ago2 RIP was performed in HEK293T siCtrl and siTRMT6/61A. d TGIRT-seq of the Ago2-bound and input fractions identified lower m¹A mismatch on tRF-3b. The color-shaded mismatch rate was calculated for all tRF-3b excluding the one from tRNA^iMet. e Volcano plot of the Ago2-bound small RNA changes upon siTRMT6/61A (differential analysis and p value by DESeq2 Wald test). Data in this figure are based on TGIRT-seq of three independent Ago2 RIP experiments in HEK293T stable Flag-HA-Ago2 cells. See also Supplementary Fig. 4 and Supplementary Table 3. Source data are provided as a Source Data file.

Fig. 6. m¹A-dependent changes in global gene silencing by tRF-3b.

a Workflow to identify potential tRF-3 target RNAs that could be regulated by m¹A status. b Ago-binding tRF-3b seeds sequences are identified from Ago2-RIP. c m¹A-dependent tRF-3b seed sequences are identified from Log2 Fold Change of the m¹A mismatch rate for the Ago2-bound fractions upon TRMT6/61A knockdown, as indicated by the color. d tRF-3 targets are globally repressed compared to the non-targets upon siTRMT61A in HEK293T. Distribution of expression changes (X-axis: Log2 fold change from RNA-seq) is visualized by CDF (Cumulative Distribution Function) plot. P value by one-sided Kolmogorov–Smirnov test to compare overall distribution between each target type versus non-targets. Significantly repressed tRF-3 targets (7/8-mers, DESeq2 adjusted p value <0.01, Log2FoldChange <−0.5) from (e) RNA-seq are selected for further validation by (f) RT-qPCR. g Dual-luciferase assay with tRF-3 target site was performed to measure the effect after siTRMT6/61A. MBTPS1 and CREB3L2 3′ UTR sequences were inserted after Renilla luciferase gene. h tRF-3004b mimic over-expression represses tRF-3 targets by RT-qPCR. i Dual-luciferase assay was performed to measure the effect after tRF-3004b mimic over-expression. tRF-3004b target sites from endogenous MBTPS1 and CREB3L2 3′ UTR were cloned as 4X tandem repeats. Data are based on b, c TGIRT-seq of three independent knock-down and Ago2 RIP experiments in Flag-HA-Ago2 HEK293T; d, e RNA-seq of two independent knock-down experiments; f n = 3, (g) n = 3, (h) n = 4, (i) n = 3 biologically independent experiments in HEK293T (mean ± SD). RT-qPCR: relative expression is normalized to siCtrl (f) or NT1 (non-targeting RNA control1 (h), and compared with ACTB expression. Dual luciferase assay: Renilla luciferase read is normalized to Firefly luciferase read and further normalized to siCtrl (g) or NT1 (i). The significance was based on student’s t test (two-tailed unpaired, *p < 0.05, **p < 0.01, N.S. = p > 0.05). Exact p values from left to right: f 9.6e−6, 9.1e−4, 0.0050, 2.2e−5, 0.018, 0.0095, 3e−4, 2.6e−5, and 0.042; g 0.49, 0.03, 0.0042; h 0.0097, 3.5e−6, 0.026, 0.0068, 0.003, 0.22, and 0.0001545; i 0.0011, 0.00063, 5.7e−5, and 4.1e−6. See also Supplementary Tables 4–5 and Supplementary Fig. 5. Source data are provided as a Source Data file.

Fig. 7. Alteration of tRF-3 m¹A levels and tRF-3 targets in bladder cancer.

a TRMT6 expression in TCGA BLCA (n = 404) compared to the normal (n = 28, TCGA and GTEX). Expression (Y-axis) is derived from RNA-seq on Log2 scale. Data visualized by GEPIA2 (p value by one-way ANOVA). Boxplot center represents median, bounds represent 25 and 75%, and whiskers show the minimum or maximum no further than 1.5 * interquartile range from the bound. b, c TRMT6 and TRMT61A protein expression by Western blots in BLCA tumor versus normal samples, with β-actin as loading control. Relative protein expression is normalized to β-actin; p values based on two-tailed unpaired student’s t test. d, e m¹A mismatch measured by TGIRT-seq in BLCA tumor samples and the paired normal (n = 5). d Overall tRF-3b normalized mismatch index (y-axis, index defined by mismatch% over non-mismatch%) is plot against TRMT6 expression (x-axis, log2 scale from RNA-seq). Significantly higher tRF-3b mismatch index is observed in tumors (p value by paired student’s t test). e tRF-3b A4-mismatch% is compared between tumor and paired normal for each patient (red color in heatmap suggests higher mismatch% in tumor, p value by unpaired two-sided Wilcoxon test). f, g Global de-repression of tRF-3b targets observed in BLCA tumor compared to the paired normal (f, n = 5) or in high-TRMT6 expressed (upper quartile) BLCA tumor (n = 102) compared to the low-TRMT6 expressed (bottom quartile) tumors (n = 102) in TCGA (g). CDF (Cumulative Distribution Function) p value is calculated by one-sided Kolmogorov–Smirnov test to compare overall distribution between tRF-3 targets versus non-targets. h Correlation between TRMT6/61A expression and tRF-3 targets in TCGA BLCA patients (n = 404). Pearson correlation was based on mRNA expression levels in TPM (transcripts per million) by GEPIA2 (visualized on log2 scale). Correlation coefficient summarized by color (red = positive correlation, blue = negative correlation), p value from paired correlation test, *p < 0.05, **p < 0.01. Exact p values from up to down: e 0.0075, 0.0075, 0.0075, 0.025, 0.12, 0.66, 0.66, 0.66, 0.0075, 0.65, 0.66, 0.12, and 0.66; h TRMT6 – 9.1e−8, 7.8e−9, 0.019, 0.78, 8,4e−9, 1.2e−6, 0, 0, 0.005, 4.5e−8; TRMT61A – 9.1e−8, 0.81, 0.08, 0.37, 0.68, 0.94, 0.29, 0.96, 0.17, and 0.47. See also Supplementary Figs. 6, 7 and Supplementary Table 6. Source data are provided as a Source Data file.

Fig. 8. tRF-3b modification by TRMT6/61A is important for maintaining the unfolded protein response.

a GSEA from RNA-seq in BLCA and HEK293T siTRMT61A data reveals Unfolded Protein Response (UPR) as a candidate pathway regulated by TRMT6/61A. b Gene Set Enrichment Analysis shows UPR as positively enriched pathway in BLCA tumor compared to the paired normal samples (n = 5, x-axis: ranking based on differential analysis from upregulated genes to downregulated genes). c TRMT6/61A regulates UPR response by UPR reporter assay in HEK293T and the BLCA tumor cell line T24. Three tandem repeats of ATF6 response element ERSE2 is inserted in the Firefly luciferase gene promoter. UPR response is measured by Firefly luciferase signals divided by Renilla luciferase signals (on co-transfected plasmid) and normalized to siCtrl basal level. Tunicamycin (Tm) is added to trigger UPR response. d m¹A status on tRF-3 regulates UPR reporter output. tRF-3004b synthetic mimic and UPR reporter plasmid (same as above) are co-transfected into HEK293T. NT1 and NT2 are non-targeting control RNAs. e, f tRF-3004 targets have increased levels after tRF-3004 knock-down and tunicamycin treatment in HEK293T cells. Data are based on c, d n = 3 independent experiments, e, f n = 2 independent experiments. c, d Data represent mean ± SD, p value based on student’s t test (two-tailed paired sample, *p < 0.05, **p < 0.01). Exact p values from left to right: c 0.00019, 0.00409, 0.00088, and 0.00062; d 0.00013, 0.00066. See also Supplementary Fig. 7. Source data are provided as a Source Data file.