Differential Expression Profiles and Functional Prediction of tRNA-Derived Small RNAs in Rats After Traumatic Spinal Cord Injury

Abstract

Spinal cord injury (SCI) is mostly caused by trauma. As the primary mechanical injury is unavoidable, a focus on the underlying molecular mechanisms of the SCI-induced secondary injury is necessary to develop promising treatments for patients with SCI. Transfer RNA-derived small RNA (tsRNA) is a novel class of short, non-coding RNA, possessing potential regulatory functions in various diseases. However, the functional roles of tsRNAs in traumatic SCI have not been determined yet. We used a combination of sequencing, quantitative reverse transcription-polymerase chain reaction (qRT-PCR), bioinformatics, and luciferase reporter assay to screen the expression profiles and identify the functional roles of tsRNAs after SCI. As a result, 297 differentially expressed tsRNAs were identified in rats' spinal cord 1 day after contusion. Of those, 155 tsRNAs were significantly differentially expressed: 91 were significantly up-regulated, whereas 64 were significantly down-regulated after SCI (fold change > 1.5; P < 0.05). Bioinformatics analyses revealed candidate tsRNAs (tiRNA-Gly-GCC-001, tRF-Gly-GCC-012, tRF-Gly-GCC-013, and tRF-Gly-GCC-016) that might play regulatory roles through the mitogen-activated protein kinase (MAPK) and neurotrophin signaling pathways by targeting brain-derived neurotrophic factor (BDNF). We validated the candidate tsRNAs and found opposite trends in the expression levels of the tsRNAs and BDNF after SCI. Finally, tiRNA-Gly-GCC-001 was identified to target BDNF using the luciferase reporter assay. In summary, we found an altered tsRNA expression pattern and predicted tiRNA-Gly-GCC-001 might be involved in the MAPK and neurotrophin pathways by targeting the BDNF, thus regulating the post-SCI pathophysiologic processes. This study provides novel insights for future investigations to explore the mechanisms and therapeutic targets for SCI.

Keywords: BDNF; bioinformatics; sequencing; spinal cord injury; tRNA-derived small RNA.

Figure 1

Study design illustration. SCI, spinal cord injury; FC, fold change; qRT-PCR, quantitative real-time polymerase chain reaction.

Figure 2

The expression level analysis. (A) Heat-map of the correlation coefficient from all samples. The color in the panel represents the correlation coefficient of the two samples. The blue represents the two samples with a high correlation coefficient, and the white represents the low similarity of the two samples. (B) Primary component analysis: X, Y and Z axes represent the three main factors that affected the expression level of the sample. The colored point represents the corresponding sample, and its location shows the main character of the sample. The space distance represents the similarity of the data size. (C) Venn diagram based on the number of commonly expressed and specifically expressed tsRNAs. This diagram shows the number of tsRNAs which were expressed in both groups and also indicated the number of the specifically expressed tsRNAs. (D) Venn diagram based on the number of tsRNAs known and detected. This diagram shows the number of tsRNAs detected in this project and collected in the tRFdb.

Figure 3

Pie chart of each tsRNA subtype. (A) Pie chart of the distribution of tsRNA subtypes in the sham group. (B) Pie chart of the distribution of tsRNA subtypes in the SCI group.

Figure 4

Stacked Bar Chart. (A) The number of tsRNA subtypes against tRNA isodecoders in the sham group. (B) The number of tsRNA subtypes against tRNA isodecoders in the SCI group. (C) The Frequency of Subtype against Length of the tsRNAs in the sham group. (D) The Frequency of Subtype against Length of the tsRNAs in the SCI group.

Figure 5

Hierarchical clustering and difference in the tsRNA expression between the two groups. (A) The hierarchical clustering heat-map for the tsRNAs. The color in the panel represents the relative expression level (log2-transformed). The colored bar top at the top panel showed the sample groups, and the colored bar at the right side of the panel indicated the divisions that were performed using K-means. (B) The scatter plots of differentially expressed tsRNAs. tsRNAs above the top line (red dots, up-regulation) or below the bottom line (green dots, down-regulation) indicate more than 1.5-fold change between the two compared groups. Gray dots indicate non-differentially expressed tsRNAs. (C) The volcano plots of differentially expressed tsRNA. The values of the X and Y axes in the volcano plot are log2 transformed fold change and −log10 transformed P-values between the two groups, respectively. Red/Green circles indicate statistically significant differentially expressed tsRNAs with a fold change of no less than 1.5 and a P-value ≤ 0.05 (Red: up-regulated; Green: down-regulated). Gray circles indicate non-differentially expressed tsRNA, with FC and/or q-value that are not meeting the cut off thresholds.

Figure 6

The tsRNA/mRNA network analysis. The network included the four candidate tsRNAs and their predicted target mRNAs (Nodes in red color are tsRNAs; nodes in light-blue color are mRNAs).

Figure 7

The general GO annotations for cellular component, molecular function, and biological processes of the target mRNAs regulated by the four candidate tsRNAs. (A) tiRNA-Gly-GCC-001. (B) tRF-Gly-GCC-012. (C) tRF-Gly-GCC-013. (D) tRF-Gly-GCC-016.

Figure 8

KEGG pathway analysis of the target mRNAs regulated by the four candidate tsRNAs. (A) The bar plot shows the top 10 enrichment score values of the significantly enriched pathway for tiRNA-Gly-GCC-001. (B) The bar plot shows the top 10 enrichment score values of the significantly enriched pathway for tRF-Gly-GCC-012. (C) The bar plot shows the top 10 enrichment score values of the significantly enriched pathway for tRF-Gly-GCC-013. (D) The bar plot shows the top 10 enrichment score values of the significantly enriched pathway for tRF-Gly-GCC-016.

Figure 9

Validation of the four selected tsRNAs using qRT-PCR. Compared with the sham control, tiRNA-Gly-GCC-001, tRF-Gly-GCC-012, tRF-Gly-GCC-013, tRF-Gly-GCC-016 were all significantly upregulated in the SCI group (A–D). The data were normalized using the mean ± standard error of the mean (SEM; n = 6 per group). *P < 0.05, **P < 0.01, ***P < 0.001.

Figure 10

The relationship of tiRNA-Gly-GCC-001 and brain-derived neurotrophic factor (BDNF) validated by PCR and luciferase reporter assay. (A) The BDNF expression in the SCI group was significantly downregulated compared with the sham control (n = 6 per group). *p < 0.05. (B) An opposite change of expression between tiRNA-Gly-GCC-001 and BDNF was identified after injury, indicating that BDNF might act as a target gene inhibited by tiRNA-Gly-GCC-001. (C) A schematic of the base pair regions of the BDNF 3′ UTR and the predicted binding site for tiRNA-Gly-GCC-001 is shown with the wild type and mutant seed sequences listed below. (D) The repressive effect of tiRNA-Gly-GCC-001 on the activity of the site of BDNF 3′UTR measured by luciferase reporter assay, indicating that tiRNA-Gly-GCC-001 directly targeted BDNF. ***p < 0.001.