Dynamics of tRNA fragments and their targets in aging mammalian brain

Small RNA analysis

We used small RNA sequencing libraries from brains of the rat Rattus norvegicus²⁶ publicly available from the European Nucleotide Archive (accession number, ERA365111). Using the sra-toolkit 2.8.0 (https://trace.ncbi.nlm.nih.gov/Traces/sra/sra.cgi?view=software) we converted the sra files to fastq format. We used fastx toolkit 0.0.13 (http://hannonlab.cshl.edu/fastx_toolkit/) to clip the adapter sequences and collapse identical reads. The reads of length above 16 nts were used for downstream analysis. We collapsed and mapped the reads to the rat genome (rn6, http://hgdownload.soe.ucsc.edu/goldenPath/rn6/bigZips/), and also to the union of rat tRNAs from two independent databases (http://gtrnadb.ucsc.edu and http://trnadb.bioinf.uni-leipzig.de, also including mitochondrial tRNA genes from the latter) using Bowtie (version 1.1.1, released on 10/1/2014, http://bowtie-bio.sourceforge.net/index.shtml). Bowtie parameters were set to output only perfect matches to tRNA sequences, including the post transcriptional CCA modification. Read counts in each experiment were normalized by the total number of reads detected, and averaged across three replicates for each of the three time points (ages of 6, 14 and 22 months). For further analysis, we selected only tRFs with read counts >0.1% of total reads in every replicate.

Seed sequence analysis

We generated 7-mer subsequences of tRFs by applying a 7-nt sliding window and shifting by one nt from the 5’ to the 3’ end. We then found the counts of exact matches for each of these subsequences to the 3’ UTR regions conserved in at least 15 species (always including human, mouse and rat) out of 23 (Table 1; alignments obtained from http://www.targetscan.org/). To estimate significance of the seed matches, we compared the observed match counts for each respective 7-mer in a tRF to (i) the expected number of matches by chance (estimated from 7-mer genomic frequency) and (ii) the average number of matches of all possible 7-mers with the same nucleotide composition in conserved 3’UTRs. Genes with exact matches of 7-mer and 7-mer_1a candidate seeds to the 3’UTR were considered potential targets.

RNA sequencing analysis

For target expression analysis, we downloaded files with pre-computed transcript expression levels for the rat cerebral cortex transcriptome (GEO data series; accession number, GSE34272²⁷). The expression levels in each experiment were normalized by the total number of reads detected and averaged across three replicates for each of the three time points (ages of 6, 12 and 28 months).

Statistical evaluation of downregulation levels for miRNA and tRF targets

For each set of predicted targets of a tRF or miRNA, we compared its ratio of down-regulated/up-regulated target transcripts from young to old rats with the distribution of such ratios calculated for 1,000 randomly selected transcript sets (from the same transcriptome) of the same size as the target set (different for each tRF and miRNA). This process was repeated for three different thresholds (up-regulated by >5% / downregulated by >5%; up-regulated by >10% / downregulated by >10%; and up-regulated by >20% / downregulated by >20%), and the statistical significance of the differences observed was obtained using two-tailed t-test in the R package (version 3.3.1; www.R-project.org).

Gene ontology enrichment analysis

The predicted targets for each tRF were used as input in order to perform GO enrichment analysis. Each set of targets was uploaded to PANTHER website (http://pantherdb.org/; version 11.1²⁸) and results were obtained using the website default parameters.

tRNA fragments in rat brain

We analyzed available datasets of nine different small-RNA libraries corresponding to three replicates for three distinct time points throughout a rat lifespan. These libraries were originally produced to study miRNA in the brains of young, middle-aged and old rats²⁶. We will refer to the results associated with these three time points (6, 14 and 22 months) as Y, M and O, respectively. After mapping short RNA reads from these libraries to the union of tRNA sequences obtained from two independent databases (http://gtrnadb.ucsc.edu; http://trna.bioinf.uni-leipzig.de), we observed that the vast majority of alignments localized preferentially to a 5'- or a 3'-end of a tRNA molecule. Only 1–7% of the reads among the nine sequencing experiments aligned elsewhere on the tRNA sequence. In the datasets we analyzed, a negligible number of reads aligned to 3' U tRFs, therefore we limited our focus to 5' and 3' tRFs, for which there was extensive evidence. The two dominant tRF classes appeared likely to be generated by different mechanisms of cleavage. For instance, there was a striking consistency regarding the cleavage site location in 3’ tRFs, compared to a wider distribution of those sites in 5' tRFs (Figure 2), supporting the notion that tRFs are not byproducts of random degradation, but have specific structure-dependent cleavage sites.

Figure 1. Typical endpoints of tRNA fragments.

Typical secondary structure representation for the PheGAA tRNA gene (from http://trnadb.bioinf.uni-leipzig.de). Blue arrows point to typical endpoints for a 5' tRF. Red arrows indicate the ends of the most frequent 3' tRF. The mature tRNA molecule also contains the post-transcriptional 3' CCA modification (as does the 3' tRF). A 3' U tRF would derive from the uracil-rich trailer sequence downstream of the end of the tRNA gene (not shown).

Figure 2. Differences in length distribution of tRFs.

Length distributions for 5' (blue) and 3' (red) tRFs. tRF length is shown on the x-axis and the frequency on the y-axis. Note the much broader variability for the 5’ tRFs.

Age-related patterns of tRF abundance

We then analyzed age-related abundance of 3' and 5’ tRFs in the brain. Interestingly, we observed a very common trend of an overall monotonous increase in the 3' tRF levels with age, Y < M < O (Figure 3; Table 2). In striking contrast, the 5’ tRFs displayed a much less consistent picture (Figure 4; Table 2), with several cases of monotonous increase or decrease with age, but mostly with a visibly different pattern of change M < Y < O (Figure 5). This difference, together with the cleavage site distributions (Figure 2), suggests that distinct processes are likely responsible for the generation of 3’ and 5’ tRFs, which may be relevant for their function.

Figure 3. Age-related change in 3' tRF levels.

Abundance of 3' tRFs in rat brains for 3 distinct time points (Y is shown in green, M in blue and O in red). An average of 3 replicates for each tRF is shown on the x-axis for each time point. Error bars indicate the range of read counts. The y-axis shows the normalized tRF abundance; the numbers on the x-axis correspond to tRNA genes listed in Table 2. Y, young (6 months); M, mid-age (14 months); O, old (22 months).

Figure 4. Age-related change in 5' tRF levels.

Abundance of 5' tRFs in rat brains for 3 distinct time points (Y is shown in green, M in blue and O in red). An average of 3 replicates for each tRF is shown on the x-axis for each time point. Error bars indicate the range of read counts. The y-axis shows the normalized tRF abundance; the numbers on the x axis correspond to tRNA genes listed in Table 2. Y, young (6 months); M, mid-age (14 months); O, old (22 months).

Figure 5. Total tRF abundance.

Change in total abundance levels with age for all 5' tRFs (blue) and 3' tRFs (red). Y, young (6 months); M, mid-age (14 months); O, old (22 months).

Computational prediction of tRF seeds and their targets

Given the significant levels of tRFs in rat brains and their dynamic changes with age, we aimed to investigate their possible effect on the brain transcriptome. Although the mechanism of tRF action is yet to be elucidated, there is recent evidence suggesting an animal miRNA-like pathway of action. Previous reports have detected tRNA fragments in the cytoplasmic fraction of various human cells, including B-cells and A549 cells^22,29, as well as mouse ES cells, plant cells, fission yeast cells and carcinoma cell lines, including HepG2, LNCap and LNCap-derived C4-2^13,30–33. It has been proposed that tRFs are likely to function similarly to a traditional miRNA-like mode, using perfect complementarity of the 5' seed sequence of the tRF (typically, positions 2–8 in miRNAs) to target a subsequence within a 3' UTR of a transcript²⁹. Contrary to the above, an alternative mode of action for tRFs has been suggested by a study²², which utilized luciferase reporter assays to demonstrate that a potential seed sequence resided in the 3' end of the tRNA fragment, ruling out a 5' and a middle segment seed binding.

A search for a near-perfect complementarity of tRF sequences against transcripts yielded very few results, both in the 12 Drosophila genomes²³ and in the present study, further suggesting a targeting mode similar to animal miRNAs. Assuming such an animal-like miRNA targeting mechanism for tRFs, we further investigated the targeting mechanism of tRFs and adjusted our computational pipeline, used previously to find targets in 12 Drosophila genomes²³, to perform the tRF seed search in mammalian genomes. This pipeline functions similarly to the approach used to identify such seed sequences for miRNAs³⁴. We used 7-nt sliding windows across the length of a tRF sequence and aligned them against conserved 3' UTR regions of 23 vertebrate species. The region was considered conserved if it was found in 15 genomes (always including rat, mouse and human) out of these 23 (Table 1). We took into consideration the following match types: 7-mer-m8 (full 7-mer match), 7-mer-1a (perfect match of the first 6 nts followed by an A in the 3' end of the targeted transcript) and 8-mer-1a (perfect 7-mer match followed by an A in the 3' end), which have been extensively confirmed for miRNAs in the past³⁵. Our results in finding seeds (Table 2) demonstrate that such conserved matches can be located both on the 5' end and on the 3' end of the tRF (Figure 6), concordant with the existing experimentally validated results for tRF targeting mechanisms^22,29. A similar arrangement of the seed regions on the 5' end and on the 3' end of the tRF has also been observed in Drosophila²³.

Figure 6. Candidate seed region locations for tRFs.

The numbers of exact sequence matches in the 3’ UTR regions are plotted against the starting position of a 7-mer. Expected number of matches in 3' UTRs is shown in yellow, average number of conserved matches for all other 7-mers with the same nucleotide composition as the given window is shown in blue, and the observed number of matches in the conserved regions of 23 vertebrates is shown in red. The letters on the top left corners of each plot correspond to individual tRFs: (A) ProTGG, (B) ValAAC, (C) PheGAA, (D) AlaTGC, (E) SerAGA and (F) SerGCT.

The success in finding seeds (Table 2) was overwhelmingly in favor of shorter 3’ tRFs (6 out of 24, 25%) compared to longer 5’ tRFs (2 out of 30, 6.67%). Given the number of differences between these two tRF types, we chose to focus on 3’ tRFs for the remainder of this paper. For our meta-analysis, detailed in the sections below, we combined experimental results performed in different labs, with different brain material, and at different ages (e.g., 22 months in small RNA-seq series is quite far from both 12 and 28 months in RNA-seq series). Given the small changes in gene expression and to avoid the effects of non-monotonous changes in many 5’tRFs, we limited our subsequent analysis to six 3’ tRFs (Figure 6), which showed monotonous changes in their levels from Y to M to O and clearly defined seed sequences.

Gene ontology enrichment analysis of conserved predicted targets

Following our seed region identification for tRFs, we focused on their predicted targets with conserved seed matches within their 3' UTR (Supplementary File 1). We explored potential functions of targets of six tRFs that showed clearly defined seed sequences (Figure 6). Gene Ontology (GO) enrichment analysis of conserved predicted targets of these tRFs revealed >150 significantly enriched GO terms for biological process. “Nervous system development” was found to be consistently enriched for all six tRFs, except ProTGG. Additional biological process GO terms, such as “central nervous system development”, “neurogenesis” and “axonogenesis”, were also enriched for multiple tRF targets (Supplementary File 2). Furthermore, the same tRF targets that showed an enrichment for nervous system functionality and development were also associated with significantly enriched neuron/axon-related cellular localization terms (Supplementary File 2). Overall, these results are in agreement with our previous work on D. melanogaster²³, where we have noted a similar enrichment for biological processes related to neuronal function and development for predicted targets of tRFs increasing with age from young to adult flies. However, in addition to these functions, ProTGG and other tRFs also appeared to target transcription and splicing regulators in rat brains (Supplementary File 1 and Supplementary File 2).

Expression patterns of predicted tRF targets

We compared our observations of tRF abundance changes with age to the measured expression levels of their targets. We compared the profiles of all mRNAs in the rat cerebral cortex transcriptome²⁷ with those predicted to be targeted by miRNAs (using Targetscan³⁴) and by tRFs (using perfect matching of the identified tRF seed sequence and a conserved target sequence located in the 3' UTR of a transcript). We calculated the ratios of down- to up-regulated transcripts for the whole rat cortex transcriptome and for the targets of six 3’ tRFs (in which seeds could be clearly seen, Figure 6) and five miRNAs (Table 3) that had >500 raw reads in the old age, and, similarly to 3' tRFs, showed a monotonous increase Y < M < O. We observed that both tRF and miRNA targets were significantly enriched for down-regulated transcripts at three different regulation thresholds (Table 3). Interestingly, the enrichment for down-regulation in the union set of all neuron-related tRF targets was also significant (p<0.05) for each of these three thresholds of regulation.

Comparing the distributions of de-regulation levels from young to old age for (i) all mRNAs detected in rat cortex, (ii) for miRNA-targeted mRNAs and (iii) for tRF targets (Figure 7), we observed a consistently higher proportion of down-regulated and lower proportion of up-regulated targets in both miRNA and tRF groups of targets compared to all mRNAs. Although these proportions for mRNAs and tRF targets were generally comparable, we noted a bimodal distribution for tRF targets, whereas such bimodality was much less pronounced for miRNA targets (Figure 7). Targets for both types of small RNAs show their most prominent peaks for low levels of down-regulation with age (these range from 0 to -5% and are possibly related to targeting relevant in other cellular contexts or false positives in target predictions). However, the proportion of tRF targets down-regulated in the range of 10.0–22.5%, and thus more likely to be relevant in the brain, is consistently higher compared to that of miRNA targets. Such (relatively low) level of change is not surprising, given that miRNAs are considered to be fine-tuning the transcriptional control by post-transcriptionally modulating the target transcript levels³⁶. The age-related decrease in the mRNA levels for tRF targets is generally more pronounced than that for miRNA targets.

Figure 7. Changes in target transcript levels with age.

Transcript level changes from young to old rat brains. Distributions of changes for all detectable mRNAs in rat brains (blue), miRNA-targeted (green) and 3' tRF-targeted transcripts (red) are shown using 2.5% bins.

Finding missed tRNA genes

In our effort to identify every possible tRF present in rat brains, we took into account a union of all annotated rat tRNAs from two databases (http://gtrnadb.ucsc.edu; http://trnadb.bioinf.uni-leipzig.de). Although the latter database is rather small compared to the former, we found that it contained a handful of rat tRNA genes (to which tRF fragments did map perfectly), which were missing from the UCSC database at the time of our first analysis. Upon subsequent checking, we found that most of the missing tRNA genes have been added correctly to the most recent update of the UCSC database (not including mitochondrial tRNAs). However, there is a tRNA gene (tdbD00000658-GluCTC; Table 2), which aligns perfectly to the rat genome (chr17:45,642,771-45,642,843 of rn6), and which is still absent in the latest version of the UCSC database. In our analysis, we detected tRFs from all nine sequencing libraries mapping to tdbD00000658-GluCTC sequence. Together with the fact that annotating tRNAs is not a typical priority in genome sequencing projects, our observations suggest that there are potentially other tRNA genes lacking annotation in the published genomes. However, such genes appear to be sources of detectable tRFs. Hence, analysis of tRFs can have an added value of revealing unannotated tRNA genes for multiple species.