First of all, we collected more than thousands human tsRNA signatures data from tRFdb (http://genome.bioch.virginia.edu/trfdb/) [8] and tRFexplorer (https://trfexplorer.cloud/), all tsRNA sequences have mapped within the specific regions of tRNA (transfer-RNA) as described previously [15]. As statistical analysis, these tsRNA fragments are mainly derived from 5′ end (5’-tRFs) and 3′ end (3’-tRFs) of mature tRNA, as well as the 5′ leader (5′U-tRF) and 3′ trailer (tRF-1) regions of primary tRNA genes.
With computational approaches, we determined lots of tsRNAs and assessed their expression profile within the sncRNA-seq (small non-coding RNA sequence) data of TCGA-COAD dataset, the summary of our identification pipeline is shown in Fig. 1A. After filtering, a total of 192 tsRNAs with available expression abundance were identified within COAD samples (Table S1), of which about thirty-two percent derived from the 3′ trailer (tRF-1) regions of primary tRNA genes, and thirty percent tsRNAs from the 3′ end (3’-tRFs) of mature tRNA (Fig. 1C). As for the chromosomes, more than thirty-five percent of tsRNAs are originated from the human chr1 region, and nearly twenty percent originated from the chr6 region (Fig. 1B).
Base on the differential expression analysis, we discriminated a total of forty-two differentially expressed tsRNA genes (DEGs), of which fifteen were down-regulated and twenty-seven were up-regulated (Fig. 2A) between colon adenocarcinomas and non-tumor controls. Hierarchical cluster analysis and heat-map of differentially expressed tsRNA genes were performed and exhibited in Fig. 2B. The expression patterns of tsRNAs were significantly distinguished between colon cancers and non-tumor groups. As shown in the blue dotted line frame of Fig. 2B, we identified several tsRNAs (including ts-43, ts-44, tRFdb-3013a and tRFdb-3013b) that derived from tRNA-His-GUG, which is a histidyl-transfer-RNA [30].
Significance of tRFdb-3013a and tRFdb-3013b within colon adenocarcinomas
Among the identified tsRNA fragments, we are interested in the fragments derived from tRNA-His-GUG or tRNA-His-GTG gene. As shown in Fig. 3A, we have identified seven tsRNA fragments in the sncRNA-seq datasets of colon adenocarcinomas samples, there are two 3’-tRFs fragments (tRFdb-3013a and tRFdb-3013b) can re-map to the cleavage of 3′ end of several different mature tRNA-His-GUG (Fig. 3A and 4A), and five tRF-1 fragments (ts-30, ts-43, ts-44, ts-1 and ts-88) that can re-mapping to the down-stream regions (35 bases) of primary tRNA-His-GTG genes. To reveal the implication of these tsRNAs in colon adenocarcinomas, the expression pattens of above tsRNAs were analyzed in the samples of TCGA-COAD datasets. It showed that the expressions of four tsRNAs (ts-43, ts-44, tRFdb-3013a and tRFdb-3013b) were significantly decreased in colon adenocarcinoma samples compared to paired non-tumor controls (Fig. 3 all P < 0.05), however, no significant difference for other tsRNAs was observed between colon adenocarcinomas and controls (both P > 0.05).
Next, we investigated the relationship between tsRNAs expression and overall survival of colon adenocarcinoma samples using Kaplan–Meier curve analysis with a log-rank comparison. Base on the median value of tsRNAs expression, colon adenocarcinoma samples were divided into low expression group and high expression group. As shown in Fig. 4B, the clinical survival of colon adenocarcinoma patients with low tRFdb-3013a expression was significantly worse than that of the high expression patients (P = 0.029); the survival of patients with low tRFdb-3013b expression was also worse than that of high expression patients (P = 0.026). These results suggested that down-regulated tRFdb-3013a and tRFdb-3013b (known simply as tRFdb-3013a/b) were associated with poor survival prognosis of colon adenocarcinoma patients.
The primary clinical and molecular pathology characteristics of colon adenocarcinoma patients were listed in Table S4. The correlations of tRFdb-3013a/b expression with pathology characteristics were analyzed using Chi-square test, and that were visualized by Hierarchical clustering heatmap. As shown in Fig. 5A-5B, the expression of tRFdb-3013a and tRFdb-3013b was significantly down-regulated in colon adenocarcinoma samples with four stages (both P < 0.05) compared to non-tumor controls. Moreover, tRFdb-3013a/b may tend to down-regulated in the patients with lymphatic or vascular invasion present (Fig. 5A). However, no significant relationship was found between tRFdb-3013a/b expression and other pathology characteristics including gender, grade, MLH1 silencing, methylation subtype, colon polyps and MSI status, etc. (all P > 0.05). These findings implied that tRFdb-3013a and tRFdb-3013b may serve as novel biomarkers for diagnosis and prognosis of colon adenocarcinomas.
Biological insights of tRFdb-3013a/b on colon adenocarcinoma cells.
To further explore the effects of tRFdb-3013a/b on colon adenocarcinomas. The mimic or inhibitor of tRFdb-3013a/b was used to over-expressed or inhibit the tRFdb-3013a/b in cells. Cell proliferation of colon adenocarcinomas was determined by CCK-8 assays, it showed that cell proliferation was significantly suppressed in tRFdb-3013a mimic transfected cells compared to negative controls (both P < 0.05; Fig. 6A). Transwell and scratch assays were performed to test the effect of tRFdb-3013a on invasion and migration ability of colon adenocarcinoma cells. The results displayed that migratory rate of cells transfected with tRFdb-3013a mimic was significantly reduced in comparison with negative controls (Fig. 6B); the numbers of invasive cells with tRFdb-3013a mimic were significantly reduced and that cells with tRFdb-3013a inhibitor were increased compared with those of negative controls (Fig. 6C). These indicated that tRFdb-3013a may be effective to regulate the cell proliferation, cell migration and invasion of colon adenocarcinomas.
Moreover, some enrichment analyses were performed to provide biological insights concerning the role of tRFdb-3013a/b. TCGA-COAD mRNA expression profiling data were used to find some correlated genes of tRFdb-3013a/b via Spearman’s correlation method. As shown in Fig. 7A, for gene ontology (GO) analysis, some tRFdb-3013a-related genes were enriched in the specific GO terms, including extracellular matrix organization, collagen-containing extracellular matrix, and extracellular matrix structural constituent (EMILIN1, ADAT2, FBLN1). The dot plots (Fig. 7B) showed that tRFdb-3013a-related genes were enriched in several KEGG pathway, such as phagosome (C1R). In the GSEA analyses (Fig. 7C), a molecular signature (LMTK3 target genes) was found to associated with tRFdb-3013a, such as FBLN1 and ILVBL-AS1 are both the tRFdb-3013a correlated-genes. The enrichment analyses of tRFdb-3013b were also conducted and presented in supplemental material (Figure S1). These data suggested that tRFdb-3013a/b may play an important role in tumor progression of colon adenocarcinomas.
The tRFdb-3013a/b target and regulate ST3GAL1 expression in colon adenocarcinomas.
By using bioinformatics database and tools, many potential targets of tRFdb-3013a/b were predicted and screened out. The interaction relationship of tRFdb-3013a/b and their target genes were presented in Fig. 8A. Among these targets, we picked out a molecular with strong binding score, ST3GAL1, which is a type Ⅱ membrane protein with sialyl-transferase activity [31]. The positions in 3' untranslated region (3′UTR) of ST3GAL1 mRNA contained a putative binding site of tRFdb-3013a/b (Fig. 8B). Dual-luciferase reporter assay was performed to validate their direct binding affinity, it showed that relative luciferase activity was reduced in the co-transfection of pmirRB-ST3GAL1 wild-type reporter and tRFdb-3013a/b mimic (both P < 0.05, Fig. 9B-left), but no significant difference in that mutant reporter. Moreover, the analysis of TCGA-COAD and GSE39582 datasets showed that ST3GAL1 was highly expressed in colon adenocarcinomas compare to non-tumor controls (Fig. 9A), as was mentioned above, tRFdb-3013a and tRFdb-3013b were decreased in colon adenocarcinomas (Fig. 5B). In addition, the cell experiments demonstrated that expression levels of ST3GAL1 mRNA and protein were down-regulated in colon adenocarcinoma cells transfected with tRFdb-3013a/b mimic (both P < 0.01, Fig. 9B-right and 9C). These results suggested tRFdb-3013a and tRFdb-3013b might directly target ST3GAL1 and regulate ST3GAL1 expression in colon adenocarcinomas.