The tsRNAs (tRFdb-3013a/b) serve as novel biomarkers for colon adenocarcinomas

The identifications of tsRNAs within colon adenocarcinomas data

Firstly, we collected more than thousands of tsRNA signatures from tRFdb [9] and tRFexplorer databases, the tsRNA sequences have mapped within the specific regions of tRNA (transfer-RNA) as described previously [17]. As statistical analysis, these tsRNA fragments are mainly derived from 5′ end (5’-tRFs) and 3′ end of mature tRNA, as well as the 5′ leader (5′U-tRF) and 3′ trailer regions of primary tRNA-gene. With computational approaches, we determined two hundreds of tsRNAs and assessed their expression profiles within the sncRNA-seq data of TCGA-COAD, the summary of our identification pipeline was showed in Supplementary Figure 1A. After filtering, a total of 192 tsRNAs with available expression abundances were identified within COAD samples (Supplementary Table 3), of which about thirty-two percents were derived from the 3′ trailer regions (tRF-1s) of primary tRNA genes, and thirty percent tsRNAs from the 3′ end (3’-tRFs) of mature tRNAs (Supplementary Figure 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 (Supplementary Figure 1B).

Based on the differential expression analysis between paired colon adenocarcinoma samples and non-tumor controls, we discriminated a total of forty-two differentially expressed tsRNA-genes (DEGs), of which fifteen were down-regulated and twenty-seven were up-regulated (Figure 1A and Supplementary Table 4) in paired colon adenocarcinomas compared to non-tumor controls. Hierarchical cluster analysis and heatmap of differentially expressed tsRNA genes were performed and exhibited in Figure 1B. The expression patterns of tsRNAs were significantly distinguished between paired colon adenocarcinomas and non-tumor groups. As shown in the blue dotted line frame of Figure 1B, several tsRNAs (including ts-43, ts-44, tRFdb-3013a and tRFdb-3013b) caught our attention that derived from tRNA-His-GUG, which is a histidyl-transfer-RNA [18].

Significance of tRFdb-3013a and tRFdb-3013b within colon adenocarcinomas

As shown in Figure 2A, we have identified seven tsRNA fragments derived from tRNA-His-GUG or tRNA-His-GTG gene in the sncRNA-seq datasets of colon adenocarcinomas, there are two 3’-tRFs fragments (tRFdb-3013a and tRFdb-3013b) that can re-map to the cleavages of 3′ end of several different mature tRNA-His-GUG isoforms (Supplementary Figure 2A), and five tRF-1 fragments (ts-30, ts-43, ts-44, ts-1 and ts-88) that can re-map 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 patterns of above tsRNAs were analyzed in the samples of TCGA-COAD dataset. It showed that the expressions of four tsRNAs were down-regulated in colon adenocarcinoma samples compared to paired controls (Supplementary Figure 2B, all P <0.05, Supplementary Table 5), and compared with non-tumor controls, the expressions of tRFdb-3013a and tRFdb-3013b (known simply as tRFdb-3013a/b) were significantly down-regulated in colon adenocarcinoma and/or rectum adenocarcinoma samples with four stages (Figure 2B, 2C, all P <0.01). In our clinical tissue samples detections, it also demonstrated that the tRFdb-3013a/b expressions were down-regulated in colon adenocarcinoma compared to non-tumor controls (Supplementary Figure 3A, 3B).

Subsequently, we investigated the relationship between tsRNA expressions and overall survival of colon and/or rectum adenocarcinoma samples using KM (Kaplan–Meier) curve analysis [10]. Based on the median value of tsRNAs expression, primary tumor samples were divided into low-expressions and high-expressions group [19]. As shown in Figure 3B, the clinical survival of colon adenocarcinoma patients with low-expressions of tRFdb-3013a or tRFdb-3013b was worse than that of the high-expression patients (P =0.029 and P =0.026); As for colon adenocarcinoma plus rectum adenocarcinoma (COAD-READ), the clinical survival of patients with low-expressions of tRFdb-3013a or tRFdb-3013b were also significantly worse than that of the high-expression patients (P =0.018 and P =0.011, Figure 3A). What's more, the Cox regression analysis of overall survival within COAD-READ patients were performed [20]. As shown in Figure 4A, the forest plot described that univariate Cox regression analysis of patients showed that low tRFdb-3013a expression (P =0.012), increased age (P <0.001), high AJCC pathologic stage (P <0.01 for stage III and P <0.001 for stage IV), advanced metastasis pathologic spread (P <0.001 for M1) and lymph node pathologic spread (P <0.001 for N2) were the risk factors associated with prognosis. Subsequent multivariate Cox regression analysis results revealed that low tRFdb-3013a expression (HR=1.684, P = 0.006) were an independent prognosis factor for survival of patients (Figure 4B), in addition to the increased age, advanced metastasis and lymph node pathologic spread. Similar results were also obtained from the Cox regression analysis of tRFdb-3013b. These results suggested that down-regulated tRFdb-3013a/b was associated with poor survival prognosis of colon adenocarcinoma and rectum adenocarcinoma patients.

Nevertheless, does the survival prognosis role of tRFdb-3013a/b specific to colorectal adenocarcinomas, and how about other types of tumors? In order to understand the clinical significances of tRFdb-3013a and tRFdb-3013b within multiple types of tumors, Cox regression analysis with hazard ratio test based on tRFdb-3013a/b expressions were performed within many types of tumors on TCGA datasets. As shown in the left of Figure 5, the forest plot described that tRFdb-3013a and tRFdb-3013b down-regulation were the risk factors associated with prognosis for five types of tumors in TCGA datasets, including COAD, READ, HNSC, KIRC, and LGG (all P≤0.03). Additionally, the differential expression analysis with log2 fold-change test of tRFdb-3013a and tRFdb-3013b were performed between primary tumor samples and paired non-tumor controls. As expected, the bar graphs showed that tRFdb-3013a/b expressions were markedly down-regulated in primary tumor samples compared to non-tumor controls for other six types of tumors (P<0.01, |log(FC)| >1, Figure 5-right). Taken together, the results implicated that only within colon and rectum adenocarcinoma, the expressions of tRFdb-3013a and tRFdb-3013b were both down-regulated, and could be served as the independent risk biomarkers for prognosis.

The correlations of tRFdb-3013a/b expression with primary clinical pathological characteristics were analyzed using Chi-square test, and visualized via hierarchical-clustering heatmaps [21]. As shown in Supplementary Figure 4, tRFdb-3013a/b might tend to decreased in the patients with lymphatic or vascular invasion present. However, no significant relationship was found between tRFdb-3013a/b expressions and other pathological characteristics, such as grade, gender, MLH1 silencing, methylation subtype, colon polyps and MSI status, etc. (all P > 0.05). Above all, findings demonstrated 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 roles of tRFdb-3013a/b on colon adenocarcinomas. The mimic or inhibitor was used to over-expressed or inhibit tRFdb-3013a in cancer cells (Supplementary Figure 5A). Cell proliferations of colon adenocarcinomas were detected using CCK8 assay, it showed that the proliferations were significantly suppressed in cells transfected tRFdb-3013a mimic compared to negative control (NC; both P<0.05; Figure 6A). Transwell and scratch assays were performed to test the effect of tRFdb-3013a on invasion and migration abilities of colon adenocarcinoma cells. The results displayed that migratory rates of cells transfected with tRFdb-3013a mimic were reduced in comparison with negative controls (Figure 6B and Supplementary Figure 5B); the numbers of invasive cells with tRFdb-3013a mimic were significantly reduced, and that cell numbers with tRFdb-3013a inhibitor were increased than those of negative controls (Figure 6C). These suggested that tRFdb-3013a might be effective to regulate the cell proliferations, cell migration and invasion abilities of colon adenocarcinomas.

Moreover, the enrichment analyses were performed to provide a biological insight concerning the roles of tRFdb-3013a/b, mRNA expression profiling of TCGA-COAD data was used to find the 500-top related-genes of tRFdb-3013a/b via Spearman’s correlation method. As shown in Figure 7A, for gene ontology (GO) analysis [22], some tRFdb-3013a related-genes were enriched in the specific GO terms, such as extracellular matrix organization, collagen-containing extracellular matrix, and extracellular matrix structural constituent (EMILIN1, ADAT2, FBLN1). The dot-plots (Figure 7B) displayed that tRFdb-3013a related-genes were enriched in several KEGG pathway, including phagosome (C1R). In the GSEA analyses [23] (Figure 7C), a molecular signature (LMTK3 target genes) was found to associated with tRFdb-3013a, such as FBLN1 and ILVBL-AS1 were both the tRFdb-3013a related-genes. The enrichment analysis of tRFdb-3013b were also conducted and presented in Supplementary Materials (Supplementary Figure 6). These data suggested that tRFdb-3013a/b might play a critical role in tumor progression of colon adenocarcinomas.

The tRFdb-3013a/b target and regulate ST3GAL1 expression in colon adenocarcinomas

The researches have revealed that tRNA fragments could bind to argonaute (AGO) protein and lead to gene silencing by targeting the 3' untranslated region (3'UTR) of some mRNAs in a manner similar to miRNAs [5, 11]. Via using bioinformatics database and tools [24], many potential targets of tRFdb-3013a/b were predicted and selected out, the interaction relationship of tRFdb-3013a/b and their target genes were presented in Figure 8A.

Among these targets, we picked out a molecular with strong binding score, ST3GAL1 gene, which encodes a type II membrane protein with sialyl-transferase activity [25]. The 3′UTR position of ST3GAL1 mRNA contained a putative binding site of tRFdb-3013a/b (Figure 8B and Supplementary Figure 7). Dual-luciferase reporter assays were used to validate their direct targeting affinity, it showed that relative luciferase activity was decreased in the co-transfection of pmirRB-ST3GAL1 wild-type reporter and tRFdb-3013a/b mimic (both P < 0.05, Figure 9B-left), but no significant difference in the mutant reporter. Moreover, the analysis of TCGA-COAD and GSE39582 data showed that ST3GAL1 was highly expressed in colon adenocarcinomas compare to non-tumor controls (Figure 9A, both P <0.01), as mentioned above, tRFdb-3013a and tRFdb-3013b were conversely decreased in colon adenocarcinomas (Figure 2B). In addition, the 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, Figure 9B-right, 9C). These findings indicated tRFdb-3013a and tRFdb-3013b might directly target ST3GAL1 and regulate ST3GAL1 expressions in colon adenocarcinoma cells.

Bioinformatic analysis

The RNA expression profiles of TCGA datasets were downloaded and processed using TCGAbiolinks package [8]. The differential expression analysis between paired colon adenocarcinoma samples and non-tumor controls was performed with the limma package [33]. An absolute log2 fold-change (|log₂FC|) more than one and P-value less than 0.05 were set as cut-off point. Hierarchical clustering of differentially expressed tsRNA genes was performed based on expression values to verify the difference between paired colon adenocarcinomas and non-tumors [19]. Visualization of the tsRNAs including volcano plots and heatmaps was performed with the ggplot2 [34] and pheatmap packages of R, respectively. Survival analysis with Kaplan-Meier (KM) curves plus the log-rank test was performed via the R survival and survminer packages. Hierarchical cluster analysis was performed on tsRNAs expression and primary pathology characteristics, its corresponding heatmap was visualized by pheatmap package (https://cran.r-project.org/web/packages/pheatmap/). Spearman’s correlation was calculated for each pair of tsRNA and mRNA, its scatter diagrams were plotted by ggplot2 and ggpubr package in R software [34]. Some enrichment analyses were performed to provide biology views that regard tRFdb-3013a/b-related genes, the molecular signatures database (MSigDB) was used as gene set for gene set enrichment analysis (GSEA) [23]. Gene ontology, KEGG pathway, and GSEA analysis were visualized by clusterProfiler and enrichplot packages [22, 35]. The computational binding interactions between tsRNA and targets were predicted by tRFtarget databases [24]. The interaction networks of top-50 targets and tRFdb-3013a/b were performed and visualized via Cytoscape (v3.7.2).

Cell culture and transfections

Human colon adenocarcinoma cell lines SW480 and Caco-2 were purchased from the Cell Bank of Chinese Academy of Sciences (Shanghai, China). The cells were cultured in Leibovitz's L-15 medium or Eagle's Minimum Essential Medium (MEM; Servicebio, Wuhan, China) supplemented with 10% or 20% fetal bovine serum (FBS; Servicebio) at 37° C in a humidified incubator with 5 % CO₂. The cell lines at 50–70% confluences were transfected with micromolar either mimic for tRFdb-3013a/b (tRF-3013a/b mimic), or inhibitor for tRFdb-3013a (tRF-3013a inhibitor), or negative controls (NC mimic, NC inhibitor) using the Lipofectamine RNAimax reagent (Invitrogen, Thermo Fisher Scientific, USA). The mimic and inhibitor for tRFdb-3013a/b and negative controls were purchased from RiboBio BioTech (Guangzhou, China).

RNA extraction and qRT-PCR analysis

Total RNAs were extracted from cultured cells using the RNAex reagent (AGbio, Changsha, China). One microgram total RNAs of each sample were reversely transcribed into cDNA under standard conditions by using PrimeScript RT reagent Kit with gDNA Eraser (Takara, Shiga, Japan). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed with the SYBR^® Premix DimerEraser™ (Takara) on the CFX connect real-time system (Bio-Rad, USA). The Bulge-Loop miRNA qRT-PCR Stater Kit (RiboBio, China) with specific stem-loop RT primers was used to quantify the expression of tsRNAs, following the manufacturer's protocol. RNU6 (U6) was used as the internal control for tsRNAs template normalization, and ACTB (β-actin) for mRNA template normalization. Relative quantification of gene expression was calculated by the comparative cycle-threshold (CT) method [19]. Stem-loop primers for tRFdb-3013a/b, and RNU6 were designed and synthesized by RiboBio BioTech (Guangzhou, China). The general primers for ST3GAL1 and ACTB were designed and synthesized by Sangon Biotech (Shanghai, China), and the primers sequences were listed in Supplementary Table 2.

Protein extraction and Western blot

Total proteins were extracted from cultured cells using RIPA buffer (Beyotime, Shanghai, China) and the protein concentrations were determined using a BCA protein assay kit (Servicebio). The specific primary antibody (anti-ST3GAL1 antibody diluted at 1:1000, Abcam, USA; anti-β-actin antibody diluted at 1:1500, Sigma-Aldrich, USA). The β-actin served as an endogenous protein for normalization. The protein bands were visualized using ChemiScope S6 imaging system and processing software (Clinx, Shanghai, China).

Cell proliferation, migration and invasion assays

Cell proliferation assays were performed with Cell Counting Kit-8 (CCK-8; Selleck, USA), the absorbance (A) in each well was measured at 450 nm with multimode Microplate Reader (Varioskan LUX, Thermo Fisher Scientific). Cell migration and invasion abilities were tested via the scratch and transwell assays, the images of scratch gap and invasion cells were recorded and photographed using Nikon eclipse Ti-S microscope (Tokyo, Japan).

Luciferase reporter assay

The wild-type or mutant sequences of ST3GAL1 3’UTR containing putative tRFdb-3013a and tRFdb-3013b binding sites were subcloned into pmir-RB-Report vector (RiboBio, China). SW480 cells were co-transfected with pmir-RB-Report vector plus or without tRFdb-3013a/b. The luciferase activity was measured using the Dual-luciferase Reporter Assay System (Promega, USA).

Statistical analysis

The statistical analyses were performed with the R software (https://www.r-project.org/), GraphPad Prism version 8.0 (GraphPad Software, Inc., USA) was used for graphing and analysis. Data were exhibited as means ± standard deviation (SD). The Chi-square test was used to compare the categorical variables. Regarding the numerical variables, statistical significance of differences between two groups was assessed using two-sided Student’s t test; and comparisons of multiple groups were made by one-way analysis of variance. All experiments were performed in triplicate and P values less than 0.05 were considered a statistically significant difference.

Data availability statement

The present study used secondary data which are available in the public domain. The dataset has no identifiable information of the survey participants. All data and materials used during the current study are available from the corresponding author on reasonable request. Some limitations about tsRNA identifications were elaborated and described in Supplementary Materials.