RNA-sequencing reveals the expression profiles of tsRNAs and their potential carcinogenic role in cholangiocarcinoma.

BACKGROUND: Recently, the incidence of cholangiocarcinoma (CCA) has gradually increased. As CCA has a poor prognosis, the ideal survival rate is scarce for patients. The abnormal expressed tsRNAs may regulate the progression of a variety of tumors, and tsRNAs is expected to become a new diagnostic biomarker of cancer. However, the expression of tsRNAs is obscure and should be elucidated in CCA.
METHODS: High-throughput RNA sequencing technology (RNA-seq) was utilized to determine the overall expression profiles of tsRNAs in three pairs CCA and adjacent normal tissues and to screen the tsRNAs that were differentially expressed. The target genes of dysregulated tsRNAs were predicted and the biological effects and potential signaling pathways of these target genes were explored by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. Quantitative real-time polymerase chain reaction (qRT-PCR) was used to validate 11 differentially expressed tRFs with 12 pairs CCA and adjacent normal tissues.
RESULTS: High-throughput RNA-seq totally demonstrated 535 dysregulated tsRNAs, of which 241 tsRNAs were upregulated, such as tRF-21-YLKZKWE5D，tRF-16-9NF5W8B，tRF-27-78YLKZKWE52，tRF-19-RLXN48KP，tRF-33-IK9NJ4S2I7L7DV，tRF-19-F8DHXYIV, and 294 tsRNAs were downregulated (tRF-20-739P8WQ0, tRF-34-JJ6RRNLIK898HR, tRF-17-VL8RPY5, tRF-23-YP9LON4VDP, tRF-39-EH623K76IR3DR2I2, tRF-17-18YKISM, tRF-19-Q1Q89PJZ, etc.) in CCA compared with adjacent normal tissues (|log2 [fold change] | ≥ 1 and p value <0.05). GO and KEGG enrichment analyses indicated that the target genes of dysregulated tRFs (tRF-34-JJ6RRNLIK898HR, tRF-38-0668K87SERM492V, and tRF-39-0668K87SERM492E2) were mainly enriched in the Notch signaling pathway, Hippo signaling pathway, cAMP signaling pathway and in growth hormone synthesis, secretion and action, etc. qRT-PCR result showed that tRF-34-JJ6RRNLIK898HR/tRF-38-0668K87SERM492V/tRF-39-0668K87SERM492E2 was downregulated (p = 0.021), and tRF-20-LE2WMK81 was upregulated in CCA (p = 0.033).
CONCLUSION: Differentially expressed tRFs in CCA are enriched in many pathways associated with neoplasms, which may impact the tumor progression and have potential to be diagnostic biomarkers and therapeutic targets of CCA.

INTRODUCTION

The incidence rate of cholangiocarcinoma (CCA) is second in malignant tumor of the hepatobiliary system, with hepatocellular carcinoma as the most often. Recently, the global occurrence of CCA has gradually increased, and China is one of the countries with the highest incidence. , , With low 5‐year survival rates, the global cumulative mortality rate of CCA patients has increased markedly. , , , Hepatocellular carcinoma (HCC) and intrahepatic cholangiocarcinoma are mainly treated by surgical resection. However, many patients are diagnosed at an advanced stage, commonly resulting in recurrence after resection. There are many studies on the etiology of HCC but few on CCA. Comprehending the pathogenesis of CCA will be based on the exploration of the etiology of CCA. In addition, the study of the early diagnostic biomarkers of CCA is of great significance for improving the early diagnosis rate and survival rate of patients with CCA.

tRNA‐derived small RNAs (tsRNAs) is a kind of small noncoding RNA that is produced by specific enzymes that cleave‐specific sites of tRNA. According to the different sites of enzyme cleavage, tsRNAs mainly include two parts of tRNA‐derived fragments (tRFs) and tRNA halves (tiRNA). tRFs can be consisted of i‐tRF, tRF‐1, tRF‐2, tRF‐3, and tRF‐5. 5'‐tiRNA and 3'‐tiRNA belong to tiRNA. Recently, many features of tsRNAs have been explored, mainly in terms of the following two aspects. On the one hand, some tRFs can bind to the Argonaute protein, an important member of the RNA‐induced silencing complex, in a manner similar to miRNA or piRNA, to perform biological functions. Furthermore, tRFs can directly regulate protein translation at the post‐transcriptional level: 5'‐tiRNA (5'‐tiRNAAla and 5'‐tiRNACys) containing 5'‐terminal oligoguanine motifs can replace eIF4F (eIF4F is a eukaryotic transcription initiation factor) on the m7GTP of mRNA, thus inhibiting translation initiation and generating multiple mRNA protein complexes (mRNPs) to directly regulate protein translation. In addition, it has been reported that tsRNAs are critical to the initiation and development of tumors and have influenced on the carcinogenesis of tumors. Some researchers have demonstrated that there are four tsRNAs, including AS‐tDR‐000064, AS‐tDR‐000069, AS‐tDR‐000102, and AS‐tDR‐001391, can be upregulated in pancreatic carcinoma by using RNA‐seq. Further GO and KEGG enrichment analysis showed that these four tsRNAs differentially expressed in pancreatic cancer are enriched in diverse pathways associated with cancer, such as “the PI3K/Akt pathway,” “the Ras pathway,” and “axon guidance” and have the potential to become a new biomarker of pancreatic carcinoma. Moreover, one study indicated that tRF‐Leu‐CAG was remarkably upregulated in the serum, cells, and tissues of non‐small‐cell lung cancer, and tRF‐Leu‐CAG could contribute to cells proliferation and maintain cell cycle of lung cancer cells, which suggested that tRF‐Leu‐CAG made a difference in the initiation and development of lung cancer.

At present, it has been disclosed that tsRNAs can participate in the development of all kinds of tumors, but it is still unclear whether tsRNAs has an impact on the origin of CCA. Exploration of the changes in tsRNAs in CCA compared with normal control tissues could be of great significance in determining early diagnostic markers of CCA.

In this study, we collected CCA and adjacent normal tissues from three patients and detected the overall expression profiles of tsRNAs in the two groups by high‐throughput RNA sequencing (RNA‐seq) and further analysis and screening of differentially expressed tsRNAs. Subsequently, the target genes of differentially expressed tRFs were predicted. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were taken advantage of predicting the important biological roles and potential signal pathways of differentially expressed tsRNAs in CCA. Quantitative real‐time polymerase chain reaction (qRT‐PCR) was used to validate the differentially expressed tRFs with CCA and adjacent normal tissues. In general, the main intention of our research was to detect the occurrence and etiology of CCA, to identify early diagnostic biomarkers of CCA and to provide new molecular targets for the remedy of CCA.

MATERIALS AND METHODS

Patients and samples

Our research was implemented according to the Declaration of Helsinki and supported by the research ethics committee of the First Hospital of “China Medical University.” The written informed consents of our study were signed by all patients. We gathered CCA and adjacent normal tissues of 15 patients at the First Affiliated Hospital of “China medical university” between September and December 2020. The criteria of patient selection are as follows: (1) The patient had not received radiotherapy or chemotherapy before surgical resection; (2) patients are meeting the indications for surgical resection; (3) postoperative pathological diagnosis was cholangiocarcinoma, not hepatocellular carcinoma and lymphoid sarcoma. The patient information is listed in Table S1. After collection, the samples were frozen in RNAlater (Life Technologies) storage solution and stored at −80°C for RNA extraction. All specimens were confirmed to be CCA by histopathology. All patients were not treated with radiotherapy, chemotherapy, or targeted therapy before operation.

Preparation of RNA Library and RNA sequencing

Three pairs of tissues from CCA patients were used to acquire the extract of total RNA. High‐throughput sequencing technology was supplied by CloudSeq Biotech. The tRF&tiRNA sequencing library was generated by total RNA of each tissue. The steps of library preparation mainly included 3′‐adaptor ligation; 5′‐adaptor ligation; synthesis of cDNA; PCR amplification; recovery of ~150 bp PCR amplification fragment. The libraries were denatured as single‐stranded DNA molecules. The Illumina flow cells were used to capture single‐stranded DNA molecules and subsequently these molecules were amplified into clusters in situ. Finally, the sequence of 50 cycles was implemented by Illumina HiSeq sequencer following the supplier's directions.

Data analysis

The steps of generation raw data mainly included sequencing, analysis of image, base calling, and quality filtering, which were carried out by Illumina sequencer. Firstly, quality control was executed by Q30. Cutadapt software (v1.9.3) was used to trim the adaptor sequences and select the adaptor‐trimmed reads (≥16 nt). Then, the raw counts of each tRF&tiRNA (MINTbase v2.0) was calculated for all samples, considered as the amount of original expression of that tRF&tiRNA. EdgeR software (v3.16.5) was utilized to standardize the raw counts and differentially expressed tRF&tiRNAs between two groups were calculated by edgeR and filtered by fold change and p value. Raw data were uploaded to Gene Expression Omnibus (GEO), GSE147017.

Target gene prediction of differentially expressed tRFs

Similar to miRNAs, tsRNAs also can bind to Argonaute proteins. , It is reported that tRFs, similar to miRNAs, can target mRNA. , , Furthermore, some studies suggested that the seed sequences of tRFs have spectra analogous to those of miRNAs. , , In view of these findings, the target genes of the screened tRFs were predicted by miRanda (http://www.microrna.org/microrna/home.do; August 2010 release), which mainly depend on the binding sites existing in the 3'‐untranslated region (3'‐UTR). The tRFs that targeted genes more than 90 was included in subsequent analysis. Cytoscape 3.7.2 was used to visualize. All analyses were performed by RStudio 3.6.1.

Enrichment analysis of tRF's target genes

Subsequently, the clusterProfiler R package was used to perform GO (biological process, molecular function and cellular component) and KEGG enrichment analysis for the target genes of tRFs. The results were visualized by Cytoscape 3.7.2. The cutoff criteria is p.adjustment values <0.05 and q values <0.05.

qRT‐PCR validation

Differentially expressed tRFs were selected to perform qRT‐PCR validation. Total RNA was extracted from 12 pairs CCA and adjacent normal samples by TRIZOL (TIANGEN BIOTECH). Total RNA was pretreated by rtStar™ tRF&tiRNA Pretreatment Kit (Cat# AS‐FS‐005). cDNA was synthesized with the rtStar™ First‐Strand cDNA Synthesis Kit (3' and 5' adaptor) (Cat#AS‐FS‐003‐02; Arraystar). The primers were designed for Predesigned Human tRNA Primer Sets V2.0. U6 was utilized as an internal control. Quantitative Real‐time PCR was performed by Arraystar SYBR Green Real‐Time qPCR Master Mix (Cat#AS‐MR‐006‐5). The relative expression level of each tRNA‐derived fragments was calculated with 2−ΔCt. SPSSv24.0 (IBM; SPSS) and GraphPad Prism V8.0 (GraphPad software) were used to perform statistical analysis. p < 0.05 was statistically significant.

RESULTS

Expression profiles of tsRNAs in CCA

The results indicated that RNA‐seq totally detected 20,102 tsRNAs in the two groups, of which 9616 were upregulated, and 10,486 were downregulated. Difference in tsRNAs expression between CCA and adjacent normal tissues were displayed by scatter plot (Figure 1). The dysregulated tsRNAs were selected under the condition of |log2 (fold change) | ≥ 1 and p value <0.05. We found out 535 differentially expressed tsRNAs in this study. We discovered that the upregulated and downregulated of tsRNAs were 241 (Table S2) and 294 (Table S3) respectively in the CCA tissues. The volcano map showed a significant difference in tsRNAs in the CCA (Figure 2). The results of hierarchical clustering demonstrated different expression of tsRNAs between CCA and adjacent normal tissues (Figure 3). Tables 1 and 2 displayed the top 30 tsRNAs that were significantly upregulated and downregulated, respectively.

FIGURE 1

Scatter plot of the differentially expressed tRFs&tiRNA in CCA and adjacent normal tissues, with red dots indicating upregulated and green dots indicating downregulated. The default multiple change threshold is 2.0. CCA, cholangiocarcinoma; C, CCA tumor tissues; N, adjacent normal tissues; tiRNA, tRNA halves; tRFs, tRNA‐derived fragments.

FIGURE 2

Volcano plot of the differentially expressed tRFs&tiRNA. Red rectangle represents differential expression of tRFs&tiRNA in patients with CCA and adjacent normal tissues (|log2 (fold change) | ≥ 1; p < 0.05). CCA, cholangiocarcinoma; C, CCA tissues; N, adjacent normal tissues; tiRNA, tRNA halves; tRFs, tRNA‐derived fragments.

FIGURE 3

Hierarchical clustering of the differentially expressed tRFs in CCA and adjacent normal tissues. The color scale showed the expression values. The green represented relatively lower expression and red represented relatively higher expression. One tissue sample was represented by each column, and a single tRFs was represented by each row. CCA, cholangiocarcinoma; C, CCA tissues; N, adjacent normal tissues; tRFs, transfer RNA‐derived fragments.

TABLE 1

Top 30 upregulated tRFs in CCA versus adjacent normal tissues

MINTbase Unique ID (sequence derived)	Fragment sequence	logFC	p Value	FDR	Regulation	#Type	MINTbase Alternative IDs (GRCh37 assembly‐derived)	Chromosome
tRF‐21‐YLKZKWE5D	TTATACCCTTCCCGTACTACC	5.180160643	0.000127872	0.991038238	Up	3'‐tRF	trnaMT_MetCAT_MT_ + _4402_4469@50.70.21	MT
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna42_GlnCTG_6_ + _27263212_27263283@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna146_GlnCTG_6_‐_27515531_27515602@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna7_GlnCTG_15_‐_66161400_66161471@16.31.16	15
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna3_GlnCTG_17_ + _8023070_8023141@16.31.16	17
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna1_GlnCTG_6_ + _18836402_18836473@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna49_GlnCTG_6_ + _27487308_27487379@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna99_GlnCTG_6_‐_28909378_28909449@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna130_GlnTTG_6_‐_27763640_27763711@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.000127872	0.991038238	Up	i‐tRF	trna173_GlnTTG_6_‐_26311975_26312046@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.002412437	0.991038238	Up	i‐tRF	trna174_GlnTTG_6_‐_26311424_26311495@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.002412437	0.991038238	Up	i‐tRF	trna64_GlnTTG_6_ + _28557156_28557227@16.31.16	6
tRF‐16‐9NF5W8B	TGGTTAGCACTCTGGA	4.765299139	0.006323658	0.991038238	Up	i‐tRF	trna16_GlnTTG_17_ + _47269890_47269961@16.31.16	17
tRF‐27‐78YLKZKWE52	GTTGGTTATACCCTTCCCGTACTACCA	5.690385949	0.006323658	0.991038238	Up	3'‐tRF	trnaMT_MetCAT_MT_ + _4402_4469@45.71.27	MT
tRF‐19‐RLXN48KP	GGGTATGATTCTCGGTTTG	5.160191852	0.006323658	0.991038238	Up	i‐tRF	trna12_ProTGG_11_‐_75946869_75946940@16.34.19	11
tRF‐33‐IK9NJ4S2I7L7DV	ATGGGTGGTTCAGTGGTAGAATTCTCGCCTGCC	4.952823031	0.006323658	0.991038238	Up	i‐tRF	trna35_GlyGCC_1_ + _161413094_161413164@3.35.33	1
tRF‐33‐IK9NJ4S2I7L7DV	ATGGGTGGTTCAGTGGTAGAATTCTCGCCTGCC	4.952823031	0.006323658	0.991038238	Up	i‐tRF	trna37_GlyGCC_1_ + _161420467_161420537@3.35.33	1
tRF‐33‐IK9NJ4S2I7L7DV	ATGGGTGGTTCAGTGGTAGAATTCTCGCCTGCC	4.952823031	0.006323658	0.991038238	Up	i‐tRF	trna39_GlyGCC_1_ + _161427898_161427968@3.35.33	1
tRF‐33‐IK9NJ4S2I7L7DV	ATGGGTGGTTCAGTGGTAGAATTCTCGCCTGCC	4.952823031	0.006323658	0.991038238	Up	i‐tRF	trna41_GlyGCC_1_ + _161435258_161435328@3.35.33	1
tRF‐33‐IK9NJ4S2I7L7DV	ATGGGTGGTTCAGTGGTAGAATTCTCGCCTGCC	4.952823031	0.006323658	0.991038238	Up	i‐tRF	trna2_GlyGCC_21_‐_18827107_18827177@3.35.33	21
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006323658	0.991038238	Up	i‐tRF	trna61_MetCAT_6_ + _27745664_27745735@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna32_MetCAT_1_ + _153643726_153643797@14.32.19	1
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna20_MetCAT_17_‐_80452597_80452668@14.32.19	17
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna129_MetCAT_6_‐_27870271_27870342@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna142_MetCAT_6_‐_27560600_27560671@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna150_MetCAT_6_‐_27300764_27300835@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna169_MetCAT_6_‐_26330529_26330600@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.006984158	0.991038238	Up	i‐tRF	trna171_MetCAT_6_‐_26313352_26313423@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.007054099	0.991038238	Up	i‐tRF	trna2_MetCAT_6_ + _26286754_26286825@14.32.19	6
tRF‐19‐F8DHXYIV	AGCGGAAGCGTGCTGGGCC	6.121699305	0.007054099	0.991038238	Up	i‐tRF	trna151_ThrCGT_6_‐_27271568_27271639@14.32.19	6

TABLE 2

Top 30 downregulated tRFs in CCA versus adjacent normal tissues

MINTbase Unique ID (sequence derived)	Fragment sequence	logFC	p Value	FDR	Regulation	Type	MINTbase Alternative IDs (GRCh37 assembly‐derived)	Chromosome
tRF‐20‐739P8WQ0	GTGGTTGTAGTCCGTGCGAG	−6.165104314	5.10975E‐05	0.991038238	Down	i‐tRF	trnaMT_GluTTC_MT_‐_14674_14742@46.65.20	MT
tRF‐20‐739P8WQ0	GTGGTTGTAGTCCGTGCGAG	−6.165104314	5.10975E‐05	0.991038238	Down	i‐tRF	trnalookalike8_GluTTC_5_‐_93905172_93905240@46.65.20	5
tRF‐34‐JJ6RRNLIK898HR	CAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAAC	−4.753211723	5.10975E‐05	0.991038238	Down	3'‐tRF	trna3_GluTTC_13_‐_45492062_45492133@40.73.34	13
tRF‐34‐JJ6RRNLIK898HR	CAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAAC	−4.753211723	5.10975E‐05	0.991038238	Down	3'‐tRF	trna11_GluTTC_15_‐_26327381_26327452@40.73.34	15
tRF‐17‐VL8RPY5	TAGTATCCCCGCCTGTC	−6.201249994	0.000468303	0.991038238	Down	i‐tRF	trna5_AspGTC_12_ + _98897281_98897352@20.36.17	12
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000468303	0.991038238	Down	i‐tRF	trna18_GlyGCC_16_ + _70822597_70822667@11.33.23	16
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna25_GlyGCC_16_‐_70812114_70812184@11.33.23	16
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna68_GlyGCC_1_‐_161493637_161493707@11.33.23	1
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna19_GlyGCC_16_ + _70823410_70823480@11.33.23	16
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna24_GlyGCC_16_‐_70812942_70813012@11.33.23	16
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna5_GlyGCC_17_ + _8029064_8029134@11.33.23	17
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna19_GlyGCC_2_‐_157257659_157257729@11.33.23	2
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna128_GlyGCC_6_‐_27870686_27870756@11.33.23	6
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna35_GlyGCC_1_ + _161413094_161413164@11.33.23	1
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna37_GlyGCC_1_ + _161420467_161420537@11.33.23	1
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna39_GlyGCC_1_ + _161427898_161427968@11.33.23	1
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna41_GlyGCC_1_ + _161435258_161435328@11.33.23	1
tRF‐23‐YP9LON4VDP	TTCAGTGGTAGAATTCTCGCCTG	−7.673641325	0.000959116	0.991038238	Down	i‐tRF	trna2_GlyGCC_21_‐_18827107_18827177@11.33.23	21
tRF‐39‐EH623K76IR3DR2I2	ACGCGGGAGACCGGGGTTCAATTCCCCGACGGGGAGCCA	−5.821002219	0.000959116	0.991038238	Down	3'‐tRF	trna5_AspGTC_12_ + _98897281_98897352@37.75.39	12
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna25_GluCTC_2_‐_71273488_71273560@14.30.17	2
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna59_GluCTC_1_ + _249168447_249168518@14.30.17	1
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna116_GluCTC_1_‐_145399233_145399304@14.30.17	1
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna71_GluCTC_1_‐_161439189_161439260@14.30.17	1
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna74_GluCTC_1_‐_161431809_161431880@14.30.17	1
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna77_GluCTC_1_‐_161424398_161424469@14.30.17	1
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna80_GluCTC_1_‐_161417018_161417089@14.30.17	1
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna77_GluCTC_6_ + _28949976_28950047@14.30.17	6
tRF‐17‐18YKISM	AGTGGTTAGGATTCGGC	−7.30441932	0.001194796	0.991038238	Down	i‐tRF	trna87_GluCTC_6_‐_126101393_126101464@14.30.17	6
tRF‐19‐Q1Q89PJZ	GCGCCGCTGGTGTAGTGGT	−7.095420921	0.001589483	0.991038238	Down	5'‐tRF	trna34_GlyCCC_16_‐_686736_686806@1.19.19	16
tRF‐19‐Q1Q89PJZ	GCGCCGCTGGTGTAGTGGT	−7.095420921	0.001589483	0.991038238	Down	5'‐tRF	trna27_GlyCCC_2_‐_70476123_70476193@1.19.19	2

Identification target genes of differentially expressed tRFs

Based on miRanda, the target genes of differentially expressed tRFs were predicted. The tRFs with target genes more than 90 was included in subsequent analysis. Finally, the result demonstrated that only three differentially expressed tRFs met the inclusion criteria, including tRF‐38‐0668K87SERM492V, tRF‐39‐0668K87SERM492E2, and tRF‐34‐JJ6RRNLIK898HR. The differentially expressed information of these tRFs are listed in Table 3. The target genes with tRFs are visualization in Figure 4.

TABLE 3

tRFs enriched in GO and KEGG

MINTbase Unique ID (sequence derived)	Fragment sequence	logFC	p Value	FDR	Regulation	Type	MINTbase Alternative IDs (GRCh37 assembly‐derived)	Chromosome
tRF‐34‐JJ6RRNLIK898HR	CAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAAC	−4.753211723	0.045042521	0.991038238	Down	3'‐tRF	trna3_GluTTC_13_‐_45492062_45492133@40.73.34	13
tRF‐34‐JJ6RRNLIK898HR	CAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAAC	−4.753211723	0.045042521	0.991038238	Down	3'‐tRF	trna11_GluTTC_15_‐_26327381_26327452@40.73.34	15
tRF‐39‐0668K87SERM492E2	ACCCAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAACCA	−6.392275579	0.020148165	0.991038238	Down	3'‐tRF	trna3_GluTTC_13_‐_45492062_45492133@37.75.39	13
tRF‐39‐0668K87SERM492E2	ACCCAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAACCA	−6.392275579	0.020148165	0.991038238	Down	3'‐tRF	trna11_GluTTC_15_‐_26327381_26327452@37.75.39	15
tRF‐38‐0668K87SERM492V	ACCCAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAACC	−6.48157606	0.016704858	0.991038238	Down	3'‐tRF	trna3_GluTTC_13_‐_45492062_45492133@37.74.38	13
tRF‐38‐0668K87SERM492V	ACCCAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAACC	−6.48157606	0.016704858	0.991038238	Down	3'‐tRF	trna11_GluTTC_15_‐_26327381_26327452@37.74.38	15

FIGURE 4

Target genes of differentially expressed tRFs. The orange color represents tRFs and green color represents target genes. (A) tRF‐34‐JJ6RRNLIK898HR. (B) tRF‐38‐0668K87SERM492V. (C) tRF‐39‐0668K87SERM492E2.

Enrichment analyses revealing the correlation of tRFs with CCA

Subsequently, the target genes of three tRFs were performed GO and KEGG analysis by clusterProfiler R package. tRF‐38‐0668K87SERM492V enrichment analysis showed that it was related to formulate cellular transport complex, regulate different voltage‐gated channels activity, and participate in Hippo signaling pathway and Notch signaling pathway (Figure 5A). tRF‐39‐0668K87SERM492E2 was also associated with different voltage‐gated channels activity, Hippo signaling pathway, and Notch signaling pathway (Figure 5B). tRF‐34‐JJ6RRNLIK898HR was mainly involved in epithelial cell differentiation and participated in cAMP signaling pathway and growth hormone synthesis, and secretion and action (Figure 5C). The relationships between target genes and pathways are visualized in Figure 6.

FIGURE 5

Bulb map of GO and KEGG analysis of dysregulated tRFs. Enrichment degree of target genes was showed by GeneRatio. The name of enrichment pathways was showed in Y axis. The area of each node showed the number of the enriched target genes of differentially expressed tRFs (p < 0.05). (A) tRF‐38‐0668K87SERM492V. (B) tRF‐39‐0668K87SERM492E2. (C) tRF‐34‐JJ6RRNLIK898HR. GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes.

FIGURE 6

Relationship among differentially expressed tRFs, their target genes and related signaling pathways. Blue color represents differentially expressed tRFs, orange color represents signaling and green color represents target genes of tRFs.

qRT‐PCR validation of differentially expressed tRFs

A total of 11 differentially expressed tRFs were selected to perform qRT‐PCR validation by 12 pairs CCA tissues. The result showed that tRF‐34‐JJ6RRNLIK898HR/tRF‐38‐0668K87SERM492V/tRF‐39‐0668K87SERM492E2 was downregulated (p = 0.021, Figure 7A), and tRF‐20‐LE2WMK81 was upregulated in CCA (p = 0.033, Figure 7B), which were consistent with our RNA‐seq results. Other validated molecules have no statistical significance, including tRF‐19‐R118LOJX (Figure 8A), tRF‐19‐BR29N3E2 (Figure 8B), tRF‐21‐RK9P4P9L0 (Figure 8C), tRF‐16‐9NF5W8B (Figure 8D), tRF‐16‐KWEKK1B (Figure 8E), tRF‐17‐K5KKOV2/tRF‐18‐P6KP6HD2 (Figure 8F), tRF‐31‐6978WPRLXN4VE/tRF‐22‐8B8SOUPR2 (Figure 8G), and tRF‐18‐1SS2P4X/tRF‐17‐YEKPRSP/tRF‐17‐ML5FX23/tRF‐16‐V0J8O9E (Figure 8H).

FIGURE 7

Relative expression of tRFs detected by qRT‐PCR in CCA tissues (p < 0.05). Measurement data were expressed as mean ± SEM. (A) tRF‐34‐JJ6RRNLIK898HR/tRF‐38‐0668K87SERM492V/tRF‐39‐0668K87SERM492E2. (B) tRF‐20‐LE2WMK81. CCA, cholangiocarcinoma; qRT‐PCR, quantitative real‐time polymerase chain reaction.

FIGURE 8

Relative expression of tRFs detected by qRT‐PCR in CCA tissues(p > 0.05). Measurement data were expressed as mean ± SEM. (A) tRF‐19‐R118LOJX. (B) tRF‐19‐BR29N3E2. (C) tRF‐21‐RK9P4P9L0. (D) tRF‐16‐9NF5W8B. (E) tRF‐16‐KWEKK1B. (F) tRF‐17‐K5KKOV2/tRF‐18‐P6KP6HD2. (G) tRF‐31‐6978WPRLXN4VE/tRF‐22‐8B8SOUPR2. (H) tRF‐18‐1SS2P4X/tRF‐17‐YEKPRSP/tRF‐17‐ML5FX23/tRF‐16‐V0J8O9E CCA, cholangiocarcinoma. qRT‐PCR, quantitative real‐time polymerase chain reaction.

DISCUSSION

Noncoding RNAs(ncRNAs) are consisted of microRNAs (miRNAs), long ncRNAs (lncRNAs), circular RNAs (circRNAs), and tsRNAs, which are small RNAs that are not capable of encoding proteins. , , Recently, with the application of high‐throughput sequencing technology and biological information analyses in ncRNAs, researchers have revealed many new research theories, mainly including the theory that ncRNAs such as tsRNAs (including tRFs and tiRNAs) are of great importance in the initiation and development of tumors and can be used as molecular markers for tumor diagnosis or molecular targets for treatment.

Researches have proved that circRNAs and lncRNAs can be considered molecular markers of tumor diagnosis or therapeutic targets. , Similarly, tsRNAs can also act as molecular markers of disease diagnosis or molecular targets of treatments. , However, there is little information about differentially expressed tRFs in CCA. Therefore, in our study, we explored the differentially expressed tRFs in CCA by RNA‐seq and further analyzed the potential function and biological role of tRFs. In addition, we employed qRT‐PCR identify the differentially expressed tRFs with CCA tissues.

In this study, we initially detected 241 upregulated and 294 downregulated tsRNAs of CCA tissues by RNA‐seq (|log2 (fold change) | ≥ 1 and p value <0.05). It is reported that the CD5+ diffuse large B‐cell lymphoma, chronic lymphocytic leukemia, lung cancer, and breast cancer were found to have differentially expressed tRFs and could have an influence on tumor progression. , , We found the differentially expressed tRFs in CCA tissues for the first time. Analogous to results of differentially expressed tRFs before, we speculated that differentially expressed tRFs may have an impact on the progression of tumor and have potential to be diagnostic biomarkers of CCA.

Subsequently, the target genes of differentially expressed tRFs were predicted and were carried out GO and KEGG enrichment analysis. The luciferase activity analysis indicated that tRFGluTTC (tRF‐34‐JJ6RRNLIK898HR, tRF‐38‐0668K87SERM492V, and tRF‐39‐0668K87SERM492E2) could directly target Kruppel‐like factor family (KIF) including KLF9, KLF11, and KLF12 to regulate adipogenesis. In this study, we found that KIF16 was one of the targets molecular of tRFGluTTC by bioinformatics methods. Combined with previous studies, we predicted that tRFGluTTC can regulate CCA progression by directly target KIF16. The result revealed that three tRFs were downregulated, and their target genes were enriched in cancer‐related pathways (p < 0.05, adjusted p < 0.05). There are many similarities in the three differentially expressed tRFs screened. For instance, although the UniqueIDs are different in MINTbase, they share the same sequence (CAGGCGGCCCGGGTTCGACTCCCGGTGTGGGAAC), and similar regions of chromosomes 13 and 15 can be cut to generate these sequence fragments. These are interesting phenomena that we found in this research, but the specific mechanism remains to be further studied in the future. The target genes of dysregulated tRFs in CCA were mainly enriched in the Notch signaling pathway, Hippo signaling pathway, and cAMP signaling pathway and in growth hormone synthesis, secretion and action, etc. It has been reported that when the cAMP signaling pathway was activated, it may suppress the migration of breast cancer cells. The Hippo signaling pathway are relevant to the procession of hepatocellular carcinoma, breast carcinoma, and gastric carcinoma, influencing the proliferation, differentiation, and migration of tumor cells. The Notch signaling pathway is significant to the development of diversified cancers and can regulate the growth, survival, apoptosis, invasion, and migration of all kinds of tumor cells, such as pancreatic carcinoma and liver cancer cells. , , , In addition, research has identified the relationship between tRFs and the Notch signaling pathway and shown that tRF/miR‐1280 can bind with the 3'‐UTR of JAG2 mRNA to decrease the synthesis of JAG2 (a ligand of Notch signaling pathway), thereby inactivating the Notch signaling pathway and reducing the proliferation and metastasis of CRC. tRF/miR‐1280 can also depress the Notch signaling pathway, thus reducing the expression level of the CD133+ stem cell phenotype in CRC cells, decreasing tumor motility and migration ability, and hindering the formation of a tumor metastasis‐favorable microenvironment. Furthermore, the Notch signaling pathway is great of significance in the carcinogenesis of CCA. Some research claimed that the level of cyclin E protein is positive regulated by the Notch signaling pathway in CCA. Some scholars have indicated that when the Notch signaling pathway is activated, it can prompt hepatic progenitor cells and differentiate hepatocytes to turn into biliary lineage cells, thus promoting the carcinogenesis of CCA. The Notch signaling pathway may be a main driving factor of human CCA occurrence and a perspective molecular target of therapy. Our study found that the target genes of dysregulated tRFs in CCA could enrich in different cancer‐related pathway. According to the study above, we speculated that the dysregulated tRFs in CCA may result in the progression of CCA by regulating these cancer‐related enriched pathways. However, its effect on the biological behavior of tumor cells and its specific mechanism needs to be further explored.

qRT‐PCR was utilized to identify randomly selected 11 differentially expressed tRFs with 12 pairs CCA tissues. The result illustrated that only tRF‐34‐JJ6RRNLIK898HR/tRF‐38‐0668K87SERM492V/tRF‐39‐0668K87SERM492E2 and tRF‐20‐LE2WMK81 were consistent with RNA‐seq data and have statistical significance (p < 0.05). The result hinted again that these tRFs may play important role in CCA progression and is more likely to be diagnostic and prognostic biomarker as well as therapeutic targets.

In summary, we used high‐throughput RNA‐seq to determine differentially expressed tsRNAs in CCA and adjacent normal tissue. Subsequent GO and KEGG analyses of dysregulated tsRNAs in the two groups indicated that the target genes of tsRNAs were enriched in multiple cancer‐related pathways. qRT‐PCR also validated the differentially expressed tRF‐34‐JJ6RRNLIK898HR/tRF‐38‐0668K87SERM492V/tRF‐39‐0668K87SERM492E2 and tRF‐20‐LE2WMK81. Our study suggests that tsRNAs may make a difference in the origin of CCA, which may contribute us to discover the etiological mechanism about CCA, and tsRNAs may become a new molecular marker of diagnosis and a molecular target of therapy in CCA. This study built a solid foundation for the further pathogenic mechanism research of tRFs in CCA.

AUTHOR CONTRIBUTIONS

Ben‐gang Wang conceived and designed this study. Yan‐ke Li and Li‐rong Yan preformed the experiment. Ang Wang was responsible for the data analysis and performed data interpretation. Yan‐ke Li and Li‐rong Yan wrote the article. Li‐yue Jiang, Qian Xu and Ben‐gang Wang revised the article.

FUNDING INFORMATION

This work is supported partly by grants from the Natural Science Foundation of Liaoning Province in China (Ref No. 2022‐MS‐196, Ref No. 2022‐YGJC‐59).

CONFLICT OF INTEREST

All authors disclose no conflicts of interest that might bias their work.

Appendix S1

Click here for additional data file.