Identification of tRNA-Derived Fragments Expression Profile in Breast Cancer Tissues.

INTRODUCTION

Breast cancer has now become the most common cancer all over the world as a major health problem [1]. It has both a higher incidence and mortality rate in women malignant tumours [2]. Effective treatment of breast cancer depends on early diagnosis and therapy [3]. However, it is difficult to achieve this goal because the diagnosis of breast cancer is lack of sensitive and effective biomarkers [4]. In recent years, studies have shown that tRNA-derived fragments are related to various pathological conditions, even pathogenic factors, such as respiratory syncytial virus infection [5],

cancer [6-8], neurodegenerative diseases [9], inherited metabolic diseases [10], etc. The tRNA-derived fragments are roughly divided into two categories: tRFs (transfer RNA-Derived Fragments) and transfer RNA-Derived Stress-induced RNAs (or tRNA halves) with specific molecular size, nucleotide composition, physiological function and biogenesis [11]. Several studies have revealed that the expression of tRNA-derived fragments is linked with proliferation, especially in tumor cells [12], but the generation of tRNA-derived fragments may only be produced because of the rampant transcription of cancer cells [13]. In addition, increased studies also show that the expression level of tRNA halves increases under stress conditions, although stress induces tRNA halves accumulation, the concentration of full-length tRNAs does not change [14]. The composition and quantity of tRNA-derived fragments are highly dependent on cell type and disease state [15]. Telonis, et al. independently validated several tRFs from the database in clinical samples from different breast cancer subtypes which demonstrated that the heterogeneity and stability of tRNA-derived fragments make them possible to be biomarkers for cancer diagnosis [15, 16]. The recent increase in rapid and inexpensive RNA sequencing may allow us to accurately identify tRNA-derived fragments from sRNA(small RNA) deep-sequencing data and evaluate their expression in multiple cancers [17, 18]. In particular, the Mintbase and the Mintmap software based on this database can help us better identify the tRNA-derived fragments subfamily and unknown tRNA-derived fragments [19, 20]. With the help of abundant data and information, those differentially expressed tRNA-derived fragments would be further studied.

Therefore, in our study, we explored the expression level of tRNA-derived fragments in the initial diagnosis of breast cancer. We constructed a database of tRNA-derived fragments expressed in clinical breast cancer samples and identified the most differentially expressed tRNA-derived fragments. We conducted qPCR (Quantitative Real-time PCR) analysis in the clinical samples of two paired groups, to verify the difference of expression of the selected tRNA-derived fragments, and explored the biological pathway analysis of GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes biological pathway).

METHODS

Tissue Specimens

In this study, six patients with breast cancer were selected for gene screening. A further ten cases were reported then and finally, there were sixteen cases of breast cancer tissue used for independent cohort validation. The tissue samples were resected in Jiangsu Cancer Hospital (Nanjing, China), which were also paired with normal tissues adjacent to cancerous tissues. All specimens were confirmed by histology and no treatment was administered to patients prior to diagnosis. The samples were frozen freshly in liquid nitrogen and stored in -80 degrees celsius until the specimens were used. The study was approved and supervised by the Clinical Research Ethics Committee of Nanjing Medical University. Written informed consent was obtained from all the patients.

Pretreatment of tRNA-derived Fragments

Total RNA was qualified by agarose gel electrophoresis and quantified using NanoDrop ND-1000 (Shanghai, China). Total RNA was first pretreated as follows to remove some RNA modifications that interfere with small RNA-seq library construction by rtStar™ tRF&tiRNA Pretreatment Kit (Cat# AS-FS-005, Arraystar, MD, USA). The rtStar™ tRF&tiRNA Pretreatment Kit is designed to remove the modifications that interfere with small RNA cDNA library construction of qPCR. These modifications include terminal modifications that block adaptor ligation to the RNA ends and internal methylations that hinder reverse transcription for cDNA synthesis. As described in the manufacturer's instructions, 3’-aminoacyl (charged) deacylation to 3’-OH for 3’ adaptor ligation, 3’-cP (2’,3’-cyclic phosphate) removal to 3’-OH for 3’ adaptor ligation, 5’-OH (hydroxyl group) phosphorylation to 5’-P for 5’-adaptor ligation, m1A and m3C were demethylated for efficient reverse transcription.

Sequencing Library Preparation

The pretreated total RNA was used to prepare the sequencing library by rtStar™ First-Strand cDNA Synthesis Kit (3’ and 5’ adaptors) (Cat#: AS-FS-003, Arraystar, MD, USA) in the following steps: 1) 3'-adapter ligation; 2) 5'-adapter ligation; 3) cDNA (complementary DNA) synthesis; 4) PCR amplification; 5) size selection of 135~160bp PCR amplified fragments (corresponding to 15~40nt small RNA size range). The prepared tRF&tiRNA-seq libraries were finally quantified using BioAnalyzer2100(Agilent, California, USA) by Agilent DNA 1000 chip kit (Agilent, part # 5067-1504) to quantified.

Sequencing and Expression Analysis

The libraries were denatured with 0.1M NaOH to generate single-stranded DNA molecules and diluted to a loading volume of 1.3ml and loading concentration of 1.8pM. Diluted libraries were loaded onto reagent cartridge and forwarded to sequencing run on the Illumina NextSeq 500 system using NextSeq 500 V2 kit (#FC-404-2005, Illumina, California, USA), according to the manufacturer’s instructions. Raw data files in FASTQ format were generated by the Illumina sequencer which were examined by FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) software (v0.11.3). Subsequently, the 3’ adapter sequence was trimmed by Cutadapt (https://doi.org/10.14806/ej.17.1. 200) from the clean reads and the reads with lengths shorter than 15 nt were discarded. As the 5’-adaptor was also used as the sequencing primer site, the 5’-adaptor sequence is not present in the sequencing reads. The trimmed reads were recorded in the FASTA format and then were aligned to mature-tRNA and pre-tRNA sequences from GtRNAdb (Genomic tRNA Database) using NovoAlign software (v2.07.11).

Different Expression of tRNA-derived Fragments for Validation

The tRNA-derived fragments expression levels were measured and normalized as tag counts per million of total aligned tRNA reads (TPM, Trans Per Million). When comparing two groups of profile differences, we used the normalized tag number of tRNAs annotated in GtRNAdb, including the tag count of each sample. The “fold change” (i.e. the ratio of the group averages) between the groups for each tRNA-derived fragment was computed. tRNA-derived fragments having fold changes ≥2 and with FDR (false discovery rate) modified p-value<0.05 were selected as the significantly, differentially expressed tRNA-derived fragments. Image analysis and base calling were performed using Solexa pipeline v1.8 (Off-Line Base Caller software, v1.8). TPM values of all tRNA-derived fragments were plotted in the Scatter plot. The Volcano plot provides a visualization method to perform a quick visual identification of the tRNA-derived fragments displaying large-magnitude changes which are also statistically significant. Outputs analysis could be filtered by deleting the duplicate sequence and ranking the differentially expressed genes according to fold change and p-value by using Microsoft Excel’s Data/Sort & Filter functionalities. Studies have shown that using GtRNAdb database alone will generate incorrect result [15, 20]. Hence, complete matching sequences were selected, and each sequence was checked using the MINTbase v2.0 framework (https://cm.jefferson.edu/MINTbase/) to make sure that each sequence has a MINTbase unique ID. The repeated annotation of the same tRNA-derived fragments was deleted. The tRNA-derived fragments with the potential of coding protein and low expression were removed. By cluster analysis, the eligible indexes were obtained and the Heatmap was drawn. The Scatter plot, Volcano plot and Heatmap were all performed in R environment for statistical computing and graphics.

Primer Design and Quantitative Real-time PCR

Six of the high different degrees and the qualified requirements for signal values were selected. Total RNA was extracted from 16 clinical samples with TRIZOL (Invitrogen life technologies). cDNA was synthesized with the rtStar™ First-Strand cDNA Synthesis Kit (3’ and 5’ adaptor) (Cat# AS-FS-003, Arraystar, MD, USA). The primers were designed for the tRNA-derived fragments particularly by Primer 5.0 and were synthesized with 2X PCR master mix (Arraystar, MD, USA). U6 was utilized as an internal control. Quantitative Real-time PCR was performed on ViiA 7 Real-time PCR System (Applied Biosystems, MA, USA). The relative expression level of each tRNA-derived fragments was calculated with the 2-△Ct method.

GO and KEGG Biological Pathway Analysis

The Gene Ontology project provides a controlled vocabulary to describe gene and gene product attributes in any organism (http://www.geneontology.org). The ontology covers three domains: Biological Process, Cellular Component and Molecular Function. Fisher’s exact test is used to detect whether there is more overlap between the DE list and the GO annotation list than would be expected by chance. The p-value denotes the significance of GO terms enrichment in the DE genes. The lower the p-value, the more significant the GO Term. Pathway analysis is functional analysis of mapping genes to KEGG pathways (http://www.genome.jp/kegg/). We regarded the -log10 (p) as the enrichment score that indicated the significance of correlation.

Statistical Analysis

Statistical analyses were performed with SPSS 23.0 (SPSS Inc, Chicago, IL, USA) and GraphPad prism 7.0 (Graphpad Software Inc, CA, USA). Wilcoxon test was used to examine significant differences of expression of tRNA-derived fragments in cancerous and adjacent normal tissues. The association between tRNA-derived fragments expression and clinicopathological characteristics of breast cancer was analyzed by Pearson’s chi-squared test and Fisher’s exact test. The value of p <0.05 was declared statistically significant.

RESULTS

Expression Profiles of tRNA-derived Fragments

The expression levels of tRNA-derived fragments in six tissue specimens of breast cancer patients and adjacent normal tissues were analyzed. The expression signatures of tRNA-derived fragments were reviewed using scatter plot, volcano plot and hierarchical clustering analyses. The scatter plots revealed that tRNA-derived fragments were differentially expressed between breast cancer and adjacent normal tissues (Fig. ). TPM values of all tRNA-derived fragments showed a more than 2.0fold change (Default fold change value is 2.0). In the scatter plot, 496 of the 1910 total genes were up-regulated and 639 down-regulated expressed in the cancer tissue compared to the normal tissue adjacent to the cancerous tissues. The Volcano plot provides a visualization method to perform quick visual identification of the tRNA-derived fragments displaying large-magnitude changes. The plot is constructed by plotting -log10 (p) on the y-axis, and tRF&tiRNA expression log2-scaled fold change (FC ≥ 2 and p<0.05) between the two experimental groups on the x-axis (Fig. ). In the difference results list, each detected tRF&tiRNA had an ID number. In comparison to the genome, one molecule might be compared to the precursor of several tRNAs. The duplicate, mismatched and imperfect sequences were deleted, and we finally found 30 significantly differentially expressed genes. A total of 17 tRNA-derived fragments were up-regulated, and 13 tRNA-derived fragments were down-regulated (Table ). Heatmap of gene expression data obtained from the control group and the experimental group, includes the genes that were selected by the Coefficient of Variation (CV) based on TPM counts (Fig. ). Each row represents one gene and each column represents one sample. The color in the panel represents the relative expression level (log2-transformed). The sequencing raw data of all tRNA-derived fragments genes have been stored in GEO (Gene Expression Omnibus). The assigned GEO accession numbers are GSE123967 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE123967).

Validation of Quantitative Real-time PCR

The 6 selected tRNA-derived fragments were up-regulated tRF-32-XSXMSL73VL4YK, tRF-28-PSQP4PW3FJD0, tRF-33-PSQP4PW3FJI0V, tRF-31-PSQP4PW3FJI0B and down-regulated tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, respectively. Table enumerated the primer sequences and the optimal annealing temperature. U6 was used for tRNA-derived fragments template normalization. The results demonstrated that tRF-32-XSXMSL73VL4YK, tRF-32-Q99P9P9NH57SJ and tRF-17-79MP9PP were significantly expressed in 16 breast cancer tissues, compared with adjacent samples (Fig. ). The up and down-regulated expression patterns were similar to the results which were observed in the sequencing analyzed in heatmap as said above (Fig. ). Expression of tRNA-derived fragments was quantified by qPCR in breast cancer tissues and adjacent samples in which expressions were normalized by U6 RNA in each sample (Fig. ). Compared to the adjacent tissues, relative expression of tRF-32-XSXMSL73VL4YK, tRF-32-Q99P9P9NH57SJ and tRF-17-79MP9PP was differential significantly (p<0.0001).

tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXMSL73VL4YK in MINTbase and Correlations with Clinical Characteristics

We checked the information of tRF-32-Q99P9P 9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXMSL73 VL4YK

in version 2.0 of MINTbase and also listed the number and location of all genomic loci (Table ). From MINTbase, we have found that tRF-32-XSXMSL73VL4YK has 9 genomic locations and tRF-32-Q99P9P9NH57SJ has only one. Both of them have exclusive tRNA space. tRF-17-79MP9PP has 18 genomic locations but it has no exclusive tRNA space. Fig. also showed the column row links to their parental tRNA isodecoders and pointed out their anticodons. The TCGA module of MINTbase provides the information separately for each of 32 TCGA cancer types (https://doi.org/10.1101/135517). tRF-32-XSXMSL 73VL4YK and tRF-32-Q99P9P9NH57SJ have no reported database in TCGA currently but tRF-17-79MP9PP has. The correlation of tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXMSL 73VL4YK expression levels with clinicopathological factors of breast cancer patients was also analyzed (Table ). tRF-32-Q99P9P9NH57SJ expression was significantly correlated with TNM stage (p = 0.012) and lymph node metastasis (p = 0.013), but there were no obvious changes between tRF-17-79MP9PP (p = 0.614 and 0.091 respectively) and tRF-32-XSXMSL73VL4YK (p = 0.317 and 0.302 respectively). Furthermore, we examined the relationship of tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP and tRF-32-XSXMSL 73VL4YK between breast cancer subtypes and expression levels, but none of them were statistically significant (p =0.375,0.375 and 0.522 respectively).

GO and KEGG Pathway Analysis of tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXMSL 73VL4YK

For subsequent research, we predicted the target gene of tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXM SL73VL4YK through GO and KEGG pathways. The top ten significantly results of GO and KEGG pathways enrichment analysis are presented in Fig. . The GO analysis covered three domains: BP (Biological Process), CC (Cellular Component) and MF (Molecular Function). The p-value denoted the significance of GO terms enrichment in the differential genes (p<0.05) and the enrichment score was regarded as the -log(p). In the BP analysis, the most enriched tRF-32-Q99P9P9NH57SJ was in the regulation of nucleobase-containing compound metabolic process (GO: 0019219). In the CC, the analysis showed that enrichment mainly occurred at intracellular (GO:0005622). For MF, the analysis revealed that the most significant enrichment was transcription regulator activity (GO: 0140110). In the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, tRF-32-Q99P9P9NH57SJ was meaningfully enriched in the regulation of lipolysis in adipocytes (hsa04923) (Fig. ). The most significant results in BP, CC, MF and KEGG of tRF-17-79MP9PP (Fig. ) were respectively Golgi vesicle transport (GO:0048193), protein complex (GO:0043234), DNA binding transcription factor activity (GO:0003700) and AMPK signaling pathway (hsa04152). For tRF-32-XSXMSL 73VL4YK (Fig. ), the positive regulation of myelination (GO:0031643), Golgi apparatus (GO:0005794), transcriptional activator activity, RNA polymerase II transcription regulatory region sequence-specific binding (GO:0001228) and drug metabolism-other enzymes (hsa00983) were the most significantly enriched term, separately.

DISCUSSION

In the early studies, it was confirmed that tRNA-derived fragments were not a random by-product of tRNA degeneration or synthesis, but rich and novel classes of short RNA with an exact sequence structure, specific expression patterns and biological roles [21]. With the rapid development of sequencing technology in recent years, the regulatory mechanisms of tRNA-derived fragments have been further studied. In addition to the MINTbase and MINTmap [19, 20] mentioned in the introduction section, Ling-Ling Zheng et al. developed an integrated web-based computing system to identify tRNA-derived fragments from sRNA deep-sequen- cing data and evaluate their expression in multiple cancers [22]. Goodarzi et al. investigated tRNA-derived fragments suppressed growth under serum-starvation, cancer cell invasion, and metastasis by breast cancer cells. Tumor suppressive pathway by attenuating the induction of tRFs might generalize to other tRNA fragments [23]. Sobala et al. presented that 5'-tRFs could inhibit the process of protein translation without the need for complementary target sites in the mRNA [24]. Olvedy, et al. have established some tRFs as novel candidate biomarkers for the early detection of recurrent aggressive in prostate cancer [25].

Telonis, et al. have screened the tRNA-derived fragments from the data of the BRCA repository of The Cancer Genome Atlas (TCGA) at the National Institutes of Health (NIH). They have reported that the specific type of i-tRFs were identified first and verified in tested breast tumor and adjacent normal samples [15]. Our study was based on the extraction of fresh breast cancer tissues, sequencing to obtain a list of different expressed tRNA-derived fragments. Ebhardt, et al. also found that there were diversities in the expression of tRNA-derived fragments in different tissues and they were heavily decorated by RNA modifications that interfere with small RNA-seq library construction [26]. Hence, before the experiment started, we preprocessed total RNA samples using the specialized kits in accordance with the manufacturer’s instructions to ensure the endpoints. Then we sequenced tRNA-derived fragments expression profiles on Illumina NextSeq 500 system and finally identified 30 significantly expressed tRNA-derived fragments. A total of 17 tRNA-derived fragments were up-regulated, and 13 tRNA-derived fragments were down-regulated. After comparing with the Mintbase, we found that the 30 tRNA-derived fragments contained various types.

We did not select i-tRF for validation because of its special type, which required special analysis [15, 20, 27]. We selected 6 tRNA-derived fragments from other types by using qPCR in 16 paired breast cancer samples. The results suggested that tRF-32-XSXMSL73VL4YK was significantly up-regulated, tRF-32-Q99P9P9NH57SJ and tRF-17-79MP9 PP were significantly down-regulated, which were consistent with Dhahbi, et al. They found that there were significant changes in up- and down-regulated of specific 5′-tRNA halves in the serum samples of breast cancer [28].

When we investigated the relationship between tRNA-derived fragments expression and clinical characteristics, there are connections in tumour stage and lymph node metastasis of tRF-32-Q99P9P9NH57SJ with no obvious correlation in tRF-17-79MP9PP and tRF-32-XSXMSL73VL4YK. However, significant literatures has emphasized that the different abundance and expression differences of tRNA-derived fragments depend on tissue type and tissue state. Both their studies were conducted under the subtype of three-negative breast cancer [29, 30]. We are currently doing an elementary exploration of the different expressions in tRNA-derived fragments between normal tissue and clinical cancer tissue specimens. Therefore, there was no special subdivision of subtype in breast cancer when the specimens were selected. The subtypes of triple-negative and non-triple negative breast cancer were both tested in the tissues we used to sequence and validate sequencing results. No more sequencing of tissue specimens was possible due to funding constraints. More information is required to be collected and explored the implications. In the following studies, we will expand the sample size for a certain subtype of breast cancer, and conduct further analysis.

In GO and KEGG pathway analyses, for subsequent research reference, we conducted target gene prediction analysis. We have found the biological functions and potential pathways of tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXMSL73VL4YK. The top ten results of GO enrichment analysis were presented in results showed above. The results of the pathway analysis needed to be supported by the later protein functional experiments, which have laid a foundation for the follow-up mechanism research. Further studies are needed to investigate the molecular mechanism and biological function of tRNA-derived fragments.

The above studies have suggested that these genes were likely to be biomarkers in the diagnosis of breast cancer. The tRNA-derived fragments are highly enriched in biological tissues and fluids [31, 32]. Although the current screening of humoral biomarkers is mainly concentrated on miRNA, it opens up a broad prospect for us to develop less invasive and sensitive biomarkers for tRNA-derived fragments. They have high equivalent and stability in a variety of body fluids, extensive involvement in pathological processes [33-35], differently expressed in solid tumours and hematological malignancies [36], and also have a strong resolution between cancer patients and normal controls [37].

CONCLUSION

Taken together, the detection of tRNA-derived fragments expression profile provides a new diagnosis process for breast cancer. We forecast that tRNA-derived fragments would become promising biomarkers for cancer diagnosis.

Table 1

Thirty significantly differential expressed tRNA-derived fragments.

MINTbase Unique ID (sequence derived)	Sequence (5’-3’)	tDF Type	Anticodon(s)	Fold Change	p value	FDR
Up-regulated Genes
tRF-32-XSXMSL73VL4YK	TGCCGTGATCGTATAGTGGTTAGTACTCTGCG	5‘-tRF	HisGTG (n)	6.5781	0.0226	0.17708
tRF-33-PSQP4PW3FJI0W	GCCCGGCTAGCTCAGTCGGTAGAGCATGAGACT	5‘-half	LysCTT (n)	3.7331	0.0001	0.01100
tRF-33-PSQP4PW3FJI0V	GCCCGGCTAGCTCAGTCGGTAGAGCATGAGACC	5‘-half	LysCTT (n)	3.6596	0.0001	0.01100
tRF-18-8R6546D2	TCCCCGGCACCTCCACCA	3'-tRF	AlaAGC (n)AlaTGC (n)	3.3616	0.0364	0.19936
tRF-18-8R6Q46D2	TCCCCGGCATCTCCACCA	3'-tRF	AlaTGC (n)AlaCGC (n)	3.3406	0.0248	0.18437
tRF-35-PW5SVP9N15WV7W	GCCGTGATCGTATAGTGGTTAGTACTCTGCGTTGT	5‘-half	HisGTG (n)	3.281	0.0371	0.19936
tRF-34-6SXMSL73VL4YHE	GGCCGTGATCGTATAGTGGTTAGTACTCTGCGTT	5'-half	HisGTG (n)	3.1888	0.0044	0.08534
tRF-32-PSQP4PW3FJI01	GCCCGGCTAGCTCAGTCGGTAGAGCATGAGAC	5‘-half	LysCTT (n)	2.9049	0.0022	0.05282
tRF-28-PSQP4PW3FJD0	GCCCGGCTAGCTCAGTCGGTAGAGCATG	5'-tRF	LysCTT (n)	2.9021	0.0086	0.09614
tRF-33-P4R8YP9LON4VDP	GCATGGGTGGTTCAGTGGTAGAATTCTCGCCTG	5‘-half	GlyGCC (n)	2.8576	0.009	0.09614
tRF-32-P4R8YP9LON4V3	GCATGGGTGGTTCAGTGGTAGAATTCTCGCCT	5‘-half	GlyGCC (n)	2.6216	0.0175	0.14423
tRF-34-PNR8YP9LON4VHM	GCATTGGTGGTTCAGTGGTAGAATTCTCGCCTGC	5‘-half	GlyGCC (n)	2.6003	0.0167	0.14278
tRF-32-PW5SVP9N15WVN	GCCGTGATCGTATAGTGGTTAGTACTCTGCGT	5'-half	HisGTG (n)	5.5853	0.0062	0.08971
tRF-31-PSQP4PW3FJI0B	GCCCGGCTAGCTCAGTCGGTAGAGCATGAGA	5'-tRF	LysCTT (n)	4.89	0.0008	0.05255
tRF-31-PW5SVP9N15WV0	GCCGTGATCGTATAGTGGTTAGTACTCTGCG	5'-tRF	HisGTG (n)	4.6602	0.0142	0.12892
tRF-33-PW5SVP9N15WV0E	GCCGTGATCGTATAGTGGTTAGTACTCTGCGTT	5‘-half	HisGTG (n)	6.0763	0.0047	0.08533
tRF-34-XSXMSL73VL4YHE	TGCCGTGATCGTATAGTGGTTAGTACTCTGCGTT	5'-half	HisGTG (n)	8.6323	0.0017	0.05281
Down-regulated Genes
tRF-31-PIR8YP9LON4VD	GCACTGGTGGTTCAGTGGTAGAATTCTCGCC	5‘-half	ValCAC (n)GlyCCC (n)	-2.0206	0.0337	0.19601
tRF-30-PNR8YP9LON4V	GCATTGGTGGTTCAGTGGTAGAATTCTCGC	5'-tRF	GlyGCC (n)GlyCCC (n)	-2.2326	0.0017	0.05281
tRF-31-PER8YP9LON4VD	GCAATGGTGGTTCAGTGGTAGAATTCTCGCC	5‘-half	GluTTC (n)	-2.3354	0.0118	0.11362
tRF-32-Q99P9P9NH57SJ	GCTTCTGTAGTGTAGTGGTTATCACGTTCGCC	5‘-half	ValCAC (n)	-2.6476	0.0189	0.14836
tRF-20-40KK5Y93	CTAAGCCAGGGATTGTGGGT	i-tRF	ArgCCT (n)	-3.3527	0.0378	0.19936
tRF-17-79MP9PP	GTTTCCGTAGTGTAGTG	5'-tRF	ValCAC (n)ValAAC (n)	-4.8984	0.0276	0.18437
tRF-36-QKF1R3WE8RO86JD	GCGGGAGACCGGGGTTCGATTCCCCGACGGGGAGCC	3'-tRF	AspGTC (n)	-6.1712	0.0382	0.19936
tRF-24-XSXMSL732Z	TGCCGTGATCGTATAGTGGTTAGT	5'-tRF	HisGTG (n)	-9.1749	0.016	0.14278
tRF-18-BS68BFD2	AACCGGGCGGAAACACCA	3'-tRF	ValCAC (n)ValAAC (n)	-114.2616	0.0314	0.19318
tRF-24-6SXMSL732Z	GGCCGTGATCGTATAGTGGTTAGT	5'-tRF	HisGTG (n)	-127.0568	0.0321	0.19318
tRF-30-34HWH3RXSINH	CCGTGATCGTATAGTGGTTAGTACTCTGCG	i-tRF	HisGTG (n)	-137.1639	0.0327	0.19318
tRF-24-FSXMSL732Z	AGCCGTGATCGTATAGTGGTTAGT	5'-tRF	HisGTG (n)	-229.9317	0.029	0.18908
tRF-23-XSXMSL730H	TGCCGTGATCGTATAGTGGTTAG	5'-tRF	HisGTG (n)	-570.2959	0.0179	0.14556

Table 2

The primer sequences and the optimal annealing temperature in quantitative real-time PCR.

Genes	Two-way Primer Sequence	Optimal Annealing Temperature (°C)	Length(bp)
U6	F:5’GCTTCGGCAGCACATATACTAAAAT3’R:5’CGCTTCACGAATTTGCGTGTCAT3’	60	89
tRF-31-PSQP4PW3FJI0B	F:5’ AGTTCTACAGTCCGACGATCGC 3'R:5’ GATCTTCTCATGCTCTACCGACTG 3’	60	56
tRF-32-XSXMSL73VL4YK	F:5’ ATCTGCCGTGATCGTATAGTGGTT 3'R:5’ ACGTGTGCTCTTCCGATCTCG 3’	60	54
tRF-28-PSQP4PW3FJD0	F:5’ CGGCTAGCTCAGTCGGTAGA 3'R:5’ CGTGTGCTCTTCCGATCTCAT 3’	60	43
tRF-17-79MP9PP	F:5’ TCTACAGTCCGACGATCGTTTC 3'R:5’ TGCTCTTCCGATCTCACTACACTA 3’	60	48
tRF-33-PSQP4PW3FJI0V	F:5’ ATCGCCCGGCTAGCTCAGT 3'R:5’ TGTGCTCTTCCGATCTGGTCTC 3’	60	52
tRF-32-Q99P9P9NH57SJ	F:5’ TTCTACAGTCCGACGATCGCT 3'R:5’ GGCGAACGTGATAACCACTACA 3’	60	50

Table 3

tRF-32-Q99P9P9NH57SJ, tRF-17-79MP9PP, tRF-32-XSXMSL73VL4YK in MINTbase.

MINTbase Unique ID (sequence derived)		Genomic Locations		Fragment Length		MINTbase Alternative IDs (GRCh37 assembly-derived)	Anticodon(s)		Exclusive to tRNA space?
Up-regulated Genes
tRF-32-XSXMSL73VL4YK		9		32		trna111_HisGTG_1_-_147774845_147774916@-1T.31.32,trna118_HisGTG_1_-_145396881_145396952@-1T.31.32,trna16_HisGTG_1_+_146544773_146544844@-1T.31.32,trna21_HisGTG_1_+_147753471_147753542@-1T.31.32,trna1_HisGTG_15_+_45493349_45493420@-1T.31.32,trna8_HisGTG_15_-_45492611_45492682@-1T.31.32,trna9_HisGTG_15_-_45490804_45490875@-1T.31.32,trna33_HisGTG_6_+_27125906_27125977@-1T.31.32,trna7_HisGTG_9_-_14433938_14434009@-1T.31.32	HisGTG (n)		yes
Down-regulated Genes
tRF-32-Q99P9P9NH57SJ	1		32		trna152_ValCAC_6_-_27248049_27248121@1.32.32			ValCAC (n)		yes
tRF-17-79MP9PP	18		17		trna37_ValAAC_6_+_27203288_27203360@1.17.17trna136_ValAAC_6_-_27648885_27648957@1.17.17trna139_ValAAC_6_-_27618707_27618779@1.17.17trna15_ValAAC_5_-_180615416_180615488@1.17.17trna2_ValAAC_3_+_169490018_169490090@1.17.17trna12_ValAAC_5_-_180645270_180645342@1.17.17trna4_ValAAC_5_+_180591154_180591226@1.17.17trna5_ValAAC_5_+_180596610_180596682@1.17.17trna132_ValAAC_6_-_27721179_27721251@1.17.17trna133_ValCAC_6_-_27696327_27696399@1.17.17trna90_ValCAC_1_-_149684088_149684161@1.17.17trna98_ValCAC_1_-_149298555_149298627@1.17.17trna85_ValCAC_1_-_161369490_161369562@1.17.17trna10_ValCAC_5_-_180649395_180649467@1.17.17trna18_ValCAC_5_-_180529253_180529325@1.17.17trna2_ValCAC_5_+_180524070_180524142@1.17.17trna6_ValCAC_5_+_180600650_180600722@1.17.17trna9_ValCAC_6_+_26538282_26538354@1.17.17			ValCAC (n)ValAAC (n)		no

Table 4

Correlation between significantly expressed tRNA-derived fragments and their clinicopathological characteristics.

Characteristics	N	tRF-32-Q99P9P9NH57SJ Low Expression (≤ median)	tRF-32-Q99P9P9NH57SJ High Expression (> median)	p	tRF-17-79MP9PP Low Expression (≤ Median)	tRF-17-79MP9PP High Expression (> Median)	p	tRF-32-XSXMSL73VL4YK Low Expression (≤ Median)	tRF-32-XSXMSL73VL4YK High Expression (> Median)	p
Number	16	7	9	-	7	9	-	8	8	-
Age (years)	-	-	-	-	-	-	-	-	-	-
≤ 50	7	2	5	0.280	1	6	0.036*	3	4	0.614
>50	9	5	4	-	6	3	-	5	4	-
Tumor stage	-	-	-	-	-	-	-	-	-	-
I-II	8	1	7	0.012*	3	5	0.614	5	3	0.317
III- IV	8	6	2	-	4	4	-	3	5	-
Lymph node metastasis	-	-	-	-	-	-	-	-	-	-
No	10	2	8	0.013*	6	4	0.091	4	6	0.302
Yes	6	5	1	-	1	5	-	4	2	-
Subtype	-	-	-	-	-	-	-	-	-	-
Triple negative	3	2	1	0.375	2	1	0.375	1	2	0.522
Non-tri-negative	13	5	8	-	5	8	-	7	6	-

* Statistical significance (p<0.05).