The Liver Fluke Opisthorchis felineus Exosomal tRNA-Derived Small RNAs as Potential Mediators of Host Manipulation

TL;DR

This study shows that liver fluke exosomes contain specific tRNA fragments that may influence human cell behavior during infection.

Contribution

First characterization of tRNA-derived small RNAs in the liver fluke Opisthorchis felineus and their potential role in host manipulation.

Findings

Exosome-like vesicles from O. felineus are enriched with specific tRNA-derived small RNAs.

These tsRNAs target human genes involved in cell cycle, migration, and proliferation.

Target genes are significantly enriched in cholangiocytes treated with fluke exosomes.

Abstract

The role of extracellular vesicle non-coding RNAs in host–parasite interactions remains poorly understood, particularly for human liver flukes. Although tRNA-derived small RNAs (tsRNAs) are emerging as new regulatory molecules in parasite exosomes, they have not yet been characterized for the liver flukes. We performed small RNA sequencing to profile tsRNAs in the exosome-like vesicles derived from the liver fluke Opisthorchis felineus. Transcriptomic data from human cholangiocytes were analyzed to assess the enrichment of the predicted target genes among differentially expressed genes. We identified 247 functional tRNA genes in the O. felineus genome. Exosome-like vesicles were highly enriched for particular tsRNAs: derived from tRNA-Asp-GTC, tRNA-Ile-AAT, tRNA-Lys, tRNA-His, and tRNA-Tyr. This enrichment was independent of both genomic tRNA copy number and the amino acid composition of the trematode proteome. In silico prediction revealed that these tsRNAs target human genes involved in cell cycle, migration, and proliferation. Notably, these predicted target genes were significantly enriched among the differentially expressed genes in treated cholangiocytes. Our study provides the first evidence that O. felineus exosomes carry a specific repertoire of tsRNAs with the potential to regulate host gene networks. We propose that tsRNAs may contribute to host cell manipulation during O. felineus infection.

Figures

Full text

1. Introduction

The liver fluke Opisthorchis felineus is a trematode of the Opisthorchiidae family that infests the bile ducts of fish-eating mammals, including humans. Chronic infection with O. felineus is a known risk factor for biliary epithelial neoplasia [1,2]. Despite the epidemiological significance of this helminth, the molecular mechanisms underlying the neoplasia development remain poorly understood. Helminths secrete exosome-like vesicles (ELVs) containing non-coding RNAs, proteins, and other biomolecules that are taken up by host cells.

Among the regulatory molecules carried by exosome-like vesicles are tRNA-derived small RNAs (tsRNAs). This diverse class of small non-coding RNAs, generated from mature tRNAs or their precursors, includes tRNA halves (tHFs), tRNA-derived fragments (tRFs), and long tRNA-derived fragments (tLFs) [3]. tsRNAs are found in exosome-like vesicles across diverse species [4,5] and function as regulators in numerous biological and pathological contexts, acting as tumor suppressors, oncogenes, and modulators of protein synthesis, apoptosis, and neurodegeneration [6]. They are broadly categorized by their biogenesis and structure: tRNA halves (28–36 nucleotides) are produced in response to cellular stress (e.g., nutrient deprivation, oxidative stress) primarily through cleavage by angiogenin, yielding 5′- and 3′-halves [7,8]. In contrast, tRNA fragments result from Dicer-dependent or independent cleavage within the D- or T-loops [9], while long tRNA fragments are generated by cleavage adjacent to these loops [3].

Notably, tsRNAs have been identified in the exosome-like vesicles of numerous parasites and are hypothesized to mediate host gene regulation [10,11,12,13,14]. Their proposed mechanism mirrors that of microRNAs; they can bind to complementary sequences in target host mRNAs, leading to transcript degradation or translational repression, or they can interact with RNA-binding proteins to influence mRNA stability [15].

Given the critical gap in understanding how the liver flukes influences host cells, this study aimed to characterize the tsRNA repertoire within O. felineus exosome-like vesicles and investigate their potential regulatory functions. We combined small RNA sequencing of parasite vesicles with in silico prediction of tsRNA targets in the human genome and integrated these predictions with transcriptomic data from treated human cholangiocytes. Our findings provide first insights into mechanisms of host–parasite interaction, suggesting that vesicular tsRNAs may contribute to the pathogenic outcomes of liver fluke infection.

2. Materials and Methods

2.1. Animals

Two-month-old golden hamsters (Mesocricetus auratus) from the Animal Facility of the ICG SB RAS were used in the study. O. felineus metacercariae were isolated from freshwater fish (Leuciscus idus) caught in the Ob River near Novosibirsk. The hamsters were orally infected with 75 metacercariae each. Adult trematodes were removed from the hamsters’ bile ducts two months after the infection and washed 10 times with sterile saline [16]. The animal study protocol was approved by the Ethics Committee of the Institute of Cytology and Genetics, Siberian Branch of the Russian Academy of Sciences (protocol code 42 of 25 May 2018).

2.2. Isolation of Exosome-like Vesicles

Adult trematodes were maintained in the medium (RPMI 1640 medium, 1% glucose, 100 units/mL penicillin G, 100 μg/mL streptomycin and 0.25 μg/mL amphotericin B), then exosome-like vesicles were isolated from the medium as described previously [17]. Transmission electron microscopy of EVs derived from O. felineus is presented in Figure S1.

2.3. Small RNAs Extraction, Library Preparation and Sequencing

Small RNAs were extracted from exosome-like vesicles with the miRNeasy Serum/Plasma Kit (Qiagen, Shanghai, China). Libraries for sequencing were prepared using the MGIEasy Small RNA Library Prep Kit V2.0 (MGI Tech., Shenzhen, China).

The three DNA libraries were subjected to single-end (1 × 50) sequencing on the BGISEQ-500 platform (BGI Genomics, Hong Kong, China). Sequencing data quality was checked using FastQC (v.0.11.9). Basic sequencing statistics are presented in Table S1.

2.4. Prediction of the O. felineus tRNA

The amino acid composition of O. felineus proteins was determined based on the predicted trematode protein sequences (GCA_004794785.1).

For tRNA prediction, the O. felineus genome (GCA_004794785.1) was used, along with tRNAscan-SE 2.0 v.2.0.9 with the following parameters: Search Mode = Eukaryotic, Isotype-specific model scan = yes, Infernal first pass cutoff score = 10 [18]. Of the 4523 predicted sequences, pseudogenes (4276) were removed. The remaining sequences (247), identified as genes, were used for further analysis (Table S2).

2.5. Search for tsRNAs in Small Noncoding RNAs from O. felineus Exosome-like Vesicles

Small RNA libraries derived from O. felineus exosome-like vesicles were aligned to predicted tRNA sequences using STAR v.2.7.10b [19] (Table S3). Initial mapping revealed that many predicted tRNA sequences corresponded to different isoforms of the same tRNA, distinguished by a small number of nucleotide substitutions, insertions, or deletions. This sequence similarity caused a significant proportion of reads to map to multiple loci, precluding reliable quantitative analysis.

To resolve this, we selected a single, representative sequence for each tRNA isoacceptor molecule. This procedure yielded a non-redundant set of 50 tRNA sequences. The small RNAs were then re-mapped to this refined reference set using STAR. Read counts mapping to each tRNA sequence were quantified using samtools coverage v.1.22, (https://www.htslib.org/doc/samtools-coverage.html [accessed on 30 January 2026]). The results are presented in Table S1.

The dominant types of tRNA-derived small RNAs (tsRNAs) were classified based on their length and cleavage position relative to the canonical tRNA structure using MEGA v.11.0.13 (https://www.megasoftware.net/ [accessed on 30 January 2026]). Data visualization was performed using LibreOffice Calc v.24.2.7.2 (https://www.libreoffice.org/discover/calc/ [accessed on 30 January 2026]), and figures were assembled in Inkscape v.1.4.2 (https://inkscape.org/ru/ [accessed on 30 January 2026]).

2.6. Functional Analysis of tsRNAs from O. felineus Exosome-like Vesicles

To identify potential target genes of dominant tsRNAs in the human genome, we performed an in silico analysis using tRFtarget 2.0 [20]. It utilizes the RNAhybrid and IntaRNA algorithms to predict binding sites and then assess where their predictions overlap to find reliable binding sites based on the consensus of the two predictions. We selected for this analysis tsRNAs that represented more than 2% of all mapped reads from the trematode exosome-like vesicles. The target sequences comprised the 3′-untranslated regions (3′-UTRs) of human genes involved in key pathways, including cell cycle (MsigDB: M7963), cell migration (MsigDB: M2001), proliferation (GO:0008283), and the negative regulation of these processes (cell cycle: GO:M12829; migration: GO:0030336; proliferation: GO:M6566; apoptotic process: GO:0043066). 3′-UTR sequences were retrieved from the Ensembl Biomart database (Ensembl Genes 115, GRCh38.p14). The list of genes is presented in Table S4. An interaction between a tsRNA and a gene was considered significant if the predicted binding energy was below −15 kcal/mol for both algorithms.

2.7. Transcriptome Data Analysis

To validate the biological relevance of the target gene prediction, we analyzed the transcriptome data from human H69 cholangiocytes treated with O. felineus excretory-secretory product (NCBI BioProject (PRJNA1203356: ESP treated samples—SAMN45991479, SAMN45991480, SAMN45991481; untreated—SAMN45991473, SAMN45991474, SAMN45991475). Quality of the raw data was assessed using FastQC v.0.11.9. The reads were mapped to the human GRCh38 genome using STAR. Analysis of differentially expressed genes (DEGs) was conducted using the DESeq2 v.1.42.0 R package [21] (Wald test with a threshold of 0.05), and genes with an adjusted p-value (P_adj_ < 0.05) were defined as DEGs. The data on the transcriptome analysis are presented in Table S5. PCA plot showing transcriptome differences between untreated and ESP-treated H69 cholangiocytes is presented in Figure S2. A heatmap for DEGs in H69 cholangiocytes genes (P_adj_ < 0.05) is presented in Figure S3 (color corresponds to per-gene z-score).

The statistical significance of the enrichment of differentially expressed genes among the predicted tsRNA targets was assessed using a Pearson’s chi-square test with Yates’s correction (stats v.4.3.3 R package). The list of target genes for transcription factor c-FOS was performed from database Harmonizome 3.0 (https://maayanlab.cloud/Harmonizome/ [accessed on 30 January 2026]).

Visualization of the relationships between tsRNAs and their target genes was created using the SankeyMATIC web service (https://sankeymatic.com/ [accessed on 30 January 2026]).

3. Results

3.1. Identification of tRNA Genes and Description of the Amino Acid Composition of Proteins Based on the O. felineus Genome Data

We identified 4523 putative tRNA genes within the O. felineus genome. Based on primary sequence and secondary structure analysis, 4276 of these were classified as pseudogenes [20], leaving 247 confirmed tRNA genes. For these functional tRNAs, the anticodon and cognate amino acid were determined. The distribution of genes for each isoacceptor tRNA is shown in Figure 1.

Analysis of 21,413 O. felineus protein sequences revealed their amino acid composition (Figure 2). The most prevalent amino acid residues were leucine (9.87%) and serine (9.45%), in contrast to the least prevalent, methionine (2.02%) and tryptophan (1.13%).

3.2. Small RNA Libraries from O. felineus Exosome-like Vesicles Are Enriched in Fragments and Halves of tRNAs of the Aspartic Acid, Isoleucine, Lysine, Histidine, and Tyrosine

Initial mapping of small RNA libraries from O. felineus exosome-like vesicles to the predicted tRNA genes revealed that reads mapped redundantly to multiple genomic copies of the same tRNA. This indicates that the trematode genome contains numerous near-identical tRNA gene duplicates, differing by minor polymorphisms. Due to this multi-mapping, the number of uniquely mapped reads was insufficient for reliable quantification, as they primarily represented these polymorphic sites.

To enable a quantitative analysis, we curated a non-redundant reference set by selecting the most abundant tRNA sequence for each anticodon variant, based on the initial read mappings. The small RNA libraries were then realigned to this curated set. This analysis revealed that the dominant tsRNAs were derived from tRNA-Asp (38.66% of all mapped reads), tRNA-Ile (20.56%), tRNA-Lys (6.93%), tRNA-His (6.18%), and tRNA-Tyr (5.15%) (Figure 2).

The pronounced abundance of these specific tsRNAs prompted us to investigate its basis. We hypothesized that their enrichment in exosome-like vesicles could be driven by: (i) the frequency of their cognate amino acids in the O. felineus proteome, (ii) the genomic copy number of the corresponding tRNA genes, or (iii) their intrinsic molecular stability. To test these hypotheses, we compared the relative abundance of tsRNAs with the amino acid composition of the trematode proteome and the number of tRNA genes for each amino acid. We found no correlation, ruling out the first two factors as primary drivers of the observed tsRNA enrichment in EVs (Figure 2).

3.3. Most tsRNAs from O. felineus Exosome-like Vesicles Are 5′-tRF and 5′-tHF

To characterize the diversity of tsRNA types, we analyzed the most abundant tsRNAs derived from the dominant isoacceptor tRNAs: tRNA-Asp-GTC, tRNA-His-GTG, tRNA-Ile-AAT, tRNA-Ile-TAT, tRNA-Lys-CTT, tRNA-Lys-TTT, and tRNA-Tyr-GTA. The most prevalent fragments from these tRNAs, which collectively accounted for >60% of all mapped reads, were classified (Figure 3A).

The 5′-tHF type was the most abundant, representing 39.41% of classified fragments and originating primarily from tRNA-Asp-GTC, tRNA-His-GTG, and tRNA-Ile-AAT. The 5′-tRF type accounted for approximately 19% of fragments, derived from tRNA-Ile-AAT, tRNA-Lys-CTT, tRNA-Lys-TTT, and tRNA-Tyr-GTA. In contrast, 3′-end fragments were less common, with 3′-tHF (1.36%) originating from tRNA-Asp-GTC and 3′-tRF (0.48%) and 5′-tLF (0.29%) being the least represented types (Figure 3B).

3.4. Functional Characteristics of Major tsRNAs

To predict the regulatory functions of O. felineus tsRNAs in human cells, we selected tsRNAs that were highly enriched in exosome-like vesicles, each constituting at least 2% of all mapped reads (Figure 4A). Since tsRNAs do not have generally accepted names, for convenience we designated them according to the scheme “Amino acid”–“Anticodon”–“fragment number by abundance.” Thus, Asp-GTC-1 represents the most abundant fragment of aspartic acid tRNA with the GTC anticodon. To search for tsRNA target genes, we analyzed the 3′UTRs of genes involved in the regulation of the cell cycle, migration, proliferation, and their negative regulation, as well as those participating in the negative regulation of apoptosis, since these processes were previously clearly observed on human cholangiocytes treated with O. felineus ESP [22].

The number of predicted target genes for each cellular pathway is presented in Table 1. The highest number of target genes was identified for tsRNA 5′-tHF-Asp-GTC-3 (206 genes), and the lowest for tsRNA 5′-tHF-His-GTG-1 (34 genes).

Given that tRNA fragments likely interact with target mRNAs through a mechanism analogous to microRNA-mediated RNA interference [23], we investigated whether their predicted target genes were differentially expressed in human cholangiocytes treated with the O. felineus excretory-secretory product (which contains exosome-like vesicles [24]). We observed a significant enrichment of differentially expressed genes among the predicted tsRNA targets (Pearson’s Chi-squared test with Yates’ continuity correction = 26.573, df = 1, p-value = 2.538 × 10^−7^). Specifically, 38 out of 1301 differentially expressed genes were identified as tsRNA targets (Figure 4B).

Furthermore, we found that some genes were targeted by multiple tsRNAs simultaneously. The EREG, NF1, and S1PR2 genes were targeted by five different tRNA fragments. In contrast, 19 target genes were unique to a single tsRNA. For the tsRNA Asp-GTC-3, nine genes (IL11, SLC35F6, NDRG1, RASSF5, TGFB1l1, TMEM127, ULK1, APPL1, and ITGB4) fell into this category (Table S6).

4. Discussion

Small tRNA-derived RNAs (tsRNAs) are an emerging class of regulatory noncoding RNAs studied as disease biomarkers and signaling molecules. Some studies revealed that tsRNAs are known components of parasitic exosomes [3,14,25], their role in host–parasite interactions remains largely unexplored. Our study provides the prediction of tRNA genes in the trematode Opisthorchis felineus and an analysis of tsRNAs derived from parasitic exosome-like vesicles.

Our analysis revealed that O. felineus genome contains multiple copies of tRNA genes, often organized in clusters [26], which is in line with many other eukaryotic species. However, the genomic architecture of these loci requires further investigation.

Strikingly, the abundance of tsRNAs in exosome-like vesicles did not reflect the genomic copy number of their corresponding tRNA genes nor the amino acid composition of the O. felineus proteome. Instead, we observed an enrichment of tsRNAs derived from aspartic acid, isoleucine, histidine, tyrosine, and lysine tRNAs.

Unfortunately, no data are available on the selective packaging of tsRNAs into exosome-like vesicles. In contrast, for microRNAs the roles for some proteins for such selective packaging were revealed. In particular, some RNA-binding proteins such as Y-box-binding protein 1 (YBX1), heterogeneous nuclear ribonucleoprotein A2/B1 (hnRNP A2/B1), and synaptotagmin-binding cytoplasmic RNA-interacting protein (SYNCRIP) are involved in their selective loading into vesicles. The interaction of these proteins with RNA requires the presence of specific short sequence motifs in the RNA molecule [27]. Various post-translational modifications, such as SUMOylation or phosphorylation, also are important in the selective shuttle of microRNAs into EVs. Notably, the endoplasmic reticulum is also implicated in the formation of RNA-enriched exosomes [28].

This enrichment suggests a functional, non-stochastic role for these particular tsRNAs, probably involved in RNA interference or other regulatory processes within the host. This hypothesis is supported by some data on regulatory functions of tsRNAs in other systems. For instance, human tsRNA-Asp-GTC induce inflammatory response via the tRF-Asp-GTC/galectin-3/TLR4/NF-κB pathway in vascular cells [29], while another tsRNA suppresses tumor growth by displacing oncogenic transcripts from YBX1 under hypoxia stress [30]. Notably, tsRNA-Asp-GTC is also highly abundant in Entamoeba histolytica exosomes, is capable of binding to Argonaute proteins, and is stress-inducible [31], in line with our findings of Asp-GTC enrichment. Furthermore, a tsRNA-Ile-AAT has been shown to downregulate SERPINE1 expression in human keratinocytes [32]. These parallels strengthen the possibility that the abundant O. felineus tsRNAs are functionally active in host cell manipulation.

The enrichment of particular tsRNAs in O. felineus exosome-like vesicles suggests a targeted packaging mechanism beyond passive accumulation. While actors such as intrinsic tsRNA stability or stress-induced fragmentation could contribute, a comparison with other trematodes argues against a universal explanation. For instance, the most abundant tsRNAs in exosomes from another liver fluke Fasciola hepatica (5′-tHF-Gly-GCC, 5′-tLF-Gln-CTG, and 5′-tHF-Cys-GCA) [3] are distinct from those we identified in O. felineus. This species-specific tsRNA signature indicates that abundance is not dictated solely by general stability or stress responses but likely reflects functional specialization unique to each parasite. Nevertheless, the predominance of 5′ tsRNAs may be associated with the greater stability of these fragments compared with 3′ tsRNAs. In silico analysis has shown that 5′ halves derived from tRNA-Gly and tRNA-Glu can form homo- and heterodimers, whereas spontaneous dimerization of 5′ halves from tRNA-Gly has been observed in vitro. The formation of such structures protects 5′ tsRNAs from RNase degradation [33]. Moreover, some 5′ tsRNAs (stress-induced 5′-tsRNA-Ala and 5′-tsRNA-Cys) are capable of forming G-quadruplex (G4) structures [34], which are known to be highly stable.

This hypothesis about functional activity of tsRNAs is supported by our findings in human cholangiocytes. Treatment with the O. felineus excretory-secretory product (containing exosome-like vesicles) is known to alter the cell proliferation and migration rates [22]. We found that the predicted target genes of the tsRNAs are significantly overrepresented among the differentially expressed genes in these treated human cells. This suggests that exosomal tsRNAs are capable of modulation of these host cellular processes during infection. The observed upregulation of some target genes could be explained by non-canonical mechanisms of gene regulation, as previously documented for microRNAs [35].

Among the predicted target genes of tsRNAs, the NF1 gene deserves special attention. It is a target for five of the most abundant OF tsRNAs, and its expression level in human cholangiocytes was significantly decreased after ESP treatment. After the cellular uptake as a cargo of exosome-like vesicles via clathrin-dependent endocytosis [17], tsRNAs may bind to the NF1 mRNA, leading to the silencing of the tumor suppressor NF1 protein output, thereby activating the Ras signaling pathway [36]. The Ras pathway strongly activates the c-FOS gene, a key player in cell differentiation and proliferation. Notably, we observed c-FOS among DEGs in human cholangiocytes treated with ESP (Padj < 0.05; Log_2_FC > 0.6). c-FOS transcription factor regulates over 13,000 target genes according to the Harmonizome 3.0 database. Among the 1301 DEGs in cholangiocytes (Padj < 0.05), the mRNA expression level of 82 c-FOS known target genes was upregulated (Log_2_FC > 0.6) (Table S5). These genes are primarily associated with apoptosis and cell proliferation pathways.

Finally, we acknowledge that the host cell phenotype is likely the result of a synergistic effect from the complex mixture of biomolecules in the excretory-secretory product and exosomes [18]. The contribution of other regulatory molecules, such as microRNAs, which are known to be functional in the exosomes of other trematodes like F. hepatica, O. viverrini, and Clonorchis sinensis [3,37,38,39,40,41], requires further investigation. Furthermore, given evidence that exosomes mediate inter-species communication [29,42], the potential autocrine or paracrine effects of tsRNAs on the trematodes themselves represent an interesting new way for future research.

5. Conclusions

Overall, the mechanisms by which the liver flukes contribute to the development of biliary neoplasia, including the role of parasite-derived exosomes, remain poorly defined. This study provides the first characterization of tRNA-derived small RNAs (tsRNAs) in O. felineus exosome-like vesicles. We also identified their potential target genes in the human genome. These tsRNAs may function as host interfering regulatory molecules capable of modulating host gene expression. Unfortunately, in silico analysis alone does not allow us to establish a definitive causal relationship between the presence of specific tsRNAs in exosome-like vesicles and the altered expression of their predicted target genes in recipient cells. Another limitation is that the pathway “Negative regulation of cell migration” has only 15 genes, which may limit statistical robustness.

Nevertheless, to move from prediction to established mechanism, the functionality of these tsRNAs needs further validation. Future research is needed to employ molecular and cellular approaches to confirm the activity, such as (i) assessing cell proliferation and migration rates of human cells treated with tsRNA mimics or inhibitors; and (ii) validating target interactions using reporter constructs containing the 3′-uncoded regions of mRNA target genes. Thus, further experiments are needed to elucidate the contribution of trematode tsRNAs to the host–parasite interaction and their potential role in biological carcinogenesis.