A miRNAs catalogue from third-stage larvae and extracellular vesicles of Anisakis pegreffii provides new clues for host-parasite interplay

Characterization of extracellular vesicles released from Anisakis pegreffii by nanoparticle tracking analysis

As a first step, to verify whether A. pegreffii third-stage larvae (L3) release extracellular vesicles (= EVs hereafter), we incubated L3 in cell culture medium and, as described in the method section, treated their supernatants with an exosome precipitation solution. Size distribution and concentration of vesicles in this precipitated fraction were estimated by Nanoparticle Tracking Analysis (NTA) using the Nanosight technology. The analysis revealed the presence of particles that in large majority have a size (mean 140.5 ± 0.08 nm, mode 107.4 ± 4.6 nm; Fig. 1) fully compatible with extracellular vesicles, which usually have a diameter between 30 and 150 nm. Moreover, the mean number of particles observed per ml in the sample was 2.49 × 10¹⁰(± 1.16 × 10⁹). Details of NTA settings and results are available in the Supplementary Fig. 1.

Finite track length adjustment (FTLA) Concentration/Size image for Nanoparticle Tracking Analysis (NTA) of extracellular vesicles secreted by third stage larvae of A. pegreffii.

Deep sequencing of small RNAs from Anisakis pegreffii L3 and EVs

Small RNA fractions (< 200 nt) were extracted from A. pegreffii L3 and EVs (three biological replicates per sample type) and used for the construction of 6 small RNA libraries suitable for Illumina high-throughput sequencing. A total of ~ 170 million raw reads (L3: ~ 100 million; EVs: ~ 70 million reads) were obtained as indicated in Table 1. After quality filtering, adapter trimming and size selection (≥ 14 nt), approximately 140 million reads were retained and mapped to the A. simplex genome (AS14). Reads mapping to AS14 (a total of ~ 55 corresponding to the 38.74% of the reads), excluding those representing ribosomal RNAs, were analysed for their size distribution (Fig. 2). In third stage larvae most reads were in the range of 20–24 nt in length, with two peaks at 21 nt and 23 nt. Although these lengths are fully compatible with the mean size of miRNAs, this bimodal distribution is not typical in small RNA-seq studies. For this reason we analyzed the ten most abundant sequences in the 21 nt and 23 nt peaks and found that they all represented miRNAs and accounted for 84% and 78% of the peak content, respectively. Among highly represented miRNAs were miR-100, miR-71, miR-9, miR-5358a and miR-5360. Extracellular vesicles enriched fractions showed a substantially different size distribution with no evident peaks and a higher frequency of reads 18–23 nt in length. Analysis of the most abundant sequences in this range also showed that miRNAs were highly represented, confirming the overall suitability of our RNA-seq libraries for the planned miRNA analysis.

Bar plots with the frequency and size distribution of reads 16–38 nt in length mapping to the A. simplex genome (AS14) and subtracted of those mapping to rRNAs.

A. pegreffii miRNAs prediction and counting

Given the absence of information on miRNAs from any Anisakis species, we used two complementary approaches to obtain a list of putative miRNAs from A. pegreffii. First, we retrieved from miRBase²³ all the available miRNAs from the pig roundworm Ascaris suum, which is the species with experimentally validated miRNAs evolutionary closer to Anisakis species. We then used these miRNAs to search for putative orthologues in the Anisakis simplex genome and, using this in silico approach, we obtained a list of 97 hypothetical A. simplex miRNA precursors and 115 mature miRNAs (dataset 1). In a second approach, as described in detail in the method section, reads from our samples were used to search the A. simplex genome using the miRNA prediction software miRDeep*. This way we obtained a list of 150 putative mature (and corresponding precursors) miRNAs from A. simplex (dataset 2). By merging these two datasets, we obtained a final non-redundant list of 206 putative Anisakis mature miRNAs. Alignments of reads from our L3 and EVs samples to this list, considering only miRNAs with reads in two replicates of at least one sample, provided evidence for the expression of 156 mature miRNAs. Among these, 67 (43%, from 46 miRNA precursors) were putative orthologues of A. suum and were renamed according to the current miRNAs nomenclature, as ape-miR- followed by the identification codes used for A. suum. The remaining putative novel 89 miRNAs (57%) were named as novelMiR followed by the identification number assigned by the miRDeep* software. Most of the miRNAs were observed in the L3 samples (126) while 30 were observed also in EVs samples and no miRNAs were reported exclusively in the EVs samples. Raw data from RNA-seq of A. pegreffii small-RNAs has been deposited to SRA under the Bioproject PRJNA786753. The list of 156 A. pegreffii miRNAs with related features is available as Supplementary Table 1.

Four miRNAs (ape-miR-100a-5p, ape-miR-1-3p, ape-miR-71-5p and ape-miR-9-5p) were by far the most abundant in L3, with > 100,000 mean CPM. The remaining miRNAs could be classified in four additional categories according to CPM values: (i) > 10,000 n = 8; (ii) > 1,000 n = 15; (iii) > 100 n = 33; (iv) ≥ 10 n = 63). Interestingly, ape-miR-100a-5p was also the most abundant miRNAs in EVs (629,841 mean CPM, 8720 mean counts). The 20 most abundant miRNAs identified in L3 are listed in Table 2. The predicted secondary structures for the first three most abundant miRNAs are shown in Fig. 3.

Secondary structures of the three most abundant miRNAs observed in the infective third stage larvae of A. pegreffii with the first nucleotide of matures sequence indicated in red, according to RNA fold.

The mature miRNAs correlation and cluster analyses confirmed the overall good quality of replicates, showing a higher homogeneity in the larval sample than in extracellular vesicles sample, which maintain their clustering trend despite the higher variation, probably due to sample normalization effect on triplicates (Supplementary Fig. 2). The expression heatmap from matures miRNAs highlighted groups with specific profile signatures, of which a few corresponded to enriched miRNAs in EVs samples (Fig. 4).

Heatmap and hierarchical clustering of mature miRNAs expression profiles of the third stage larvae of A. pegreffii (L) and of its released extracellular vesicles (EX). Each line corresponds to the mean-centered log2-transformed CPM, colored according to upregulation (yellow) and downregulation (violet). Upregulated miRNAs in EVs sample are indicated with a black dot.

Sample-specific miRNA enrichment was evaluated by pairwise comparisons between the two samples (L3 vs EVs): fold change (FC) and false discovery rates (FDR) were calculated to provide statistical validation (Supplementary Table 2). Using as threshold parameters |log2(FC)|> 1 and FDR < 0.05, we found that 38 miRNAs were differentially expressed, with 26 upregulated in L3 and 12 in EVs as shown in the Volcano plot (Fig. 5). Number of differentially expressed miRNAs applying progressively more stringent statistical thresholds are shown in Table 3.

Volcano plot with the differential abundance of miRNAs in the pairwise comparisons between third-stage larvae (L3) and in the extracellular vesicles-enriched fraction (EVs) of in Anisakis pegreffii. The log2 fold change (FC) versus the negative log10 of false discovery rate (FDR) as calculated by the Fisher’s exact test are reported. Vertical dotted lines mark logFC = 2, horizontal dashed lines mark FDR threshold equal to 0.05.

Stem and Loop RT-PCR validation of miRNAs

With the aim to confirm the presence of miRNAs in the samples studied, a list of ten miRNAs has been selected for their experimental validation in Stem and Loop RT-PCR. The list has been elaborated based on abundance category for L3 (three representative of the first categories of abundance, as > 100,000 mean CPM; > 10,000 mean CPM and > 1,000 mean CPM), on the novelty and on their differential expression (upregulated in EVs). All the ten selected miRNAs were successfully amplified and confirmed using Real Time Stem and Loop PCR (novel-miR-19, miR-7-5p, novel-miR-65, novel-miR-27, miR-72-5p, novel-miR-184, miR-100a-1-5p, lin-4-5p, miR-1-3p, novel-miR-131).

The following average of Ct were obtained for the larval sample: 17.5 for miR-100a-5p; 18 for miR-1-3p; 19.1 for novel-miR-131; 21.75 for miR-72-5p; 28.6 for novel-miR-19, 29.6 for miR-7-5p; 33.1 for novel-miR-184; 35.2 for lin-4-5p; 36.5 for novel-miR-65 and novel-miR-27.

The following average of Ct were obtained for the extracellular vesicles: 26 for novel-miR-131; 26.9 for novel-miR-19; 27.5 for miR-1-3p; 28.3 for miR-100a-1-5p; 29.6 for miR-72-5p; 31.8 for novel-miR-184; 34.4 for miR-7-5p; 35 for novel-miR-27 and lin-4-5p; 36.8 for novel-miR-65. A barplot of mean Ct values and SE obtained for L3 and EVs assays is available in Supplementary Fig. 3.

Target analysis and seed conservation in parasitic helminths

The complementarity between the 3′-UTR of the mRNA target and the miRNA seed region (nucleotides 2–8) is a crucial feature of miRNA-mRNA interaction. Identity of the entire miRNA or of the seed region may be indicative of evolutionary conservation of its function.

Comparison of A. pegreffii miRNAs to those from other helminths as Nematoda (A. suum, Trichuris muris, B. malayi, Nippostrongylus brasiliensis, Heligmosomoides polygyrus, Haemonchus contortus), Trematoda (Fasciola hepatica, Schistosoma mansoni, Schistosoma japonicum) and Cestoda (Dibothriocephalus dendriticus, Mesocestoides corti, Taenia crassiceps, Taenia asiatica, Echinococcus multilocularis) showed different levels of conservation. A complete seed conservation was reported between miRNAs from A. pegreffii with those of parasitic nematodes A. suum as shown in Table 4 (85% of the most abundant in L3 and 46% of EVs enriched) followed by B. malayi (75% of the most abundant in L3 and 31% of EVs enriched). Additionally, complete seed conservation and miRNAs homology were observed with other helminths as trematodes and cestodes as well as with human miRNAs, suggesting a potential functional role conserved across phylogenetically distant taxa.

The most abundant miRNA in both larvae and extracellular vesicles is ape-miR-100a-5p that shows complete homology to miR-100 from other parasitic helminths and humans. It belongs to the miR-10 family, which besides miR-100 also includes miR-51, miR-57 and miR-99 (Rfam RF00104). The predicted putative gene target of ape-miR-100a-5p is TRIB2, a crucial gene for regulation of apoptosis and thymocyte cellular proliferation²⁴, and its dysregulation was found associated with tumors, including colorectal cancer²⁵.

Another abundant miRNA is lin-4-5p: originally identified in C. elegans in relation to developmental timing²⁶, it shows complete seed identity with miR-125 of several parasitic helminths and humans. The dendrogram obtained from the comparison of the available orthologues showed similarities between ape-lin-4 in A. suum and B. malayi (Fig. 6). miR-125 seems to be involved in fundamental aspects of parasitic infection: F. hepatica miR-125b (fhe-miR-125b) was the most abundant miRNA trafficked by EVs observed in peritoneal macrophages during the early phase of infection. Given the homology with the mammalian miRNA hsa-miR-125b, it has been suggested the hijacking of the miRNA machinery as a parasitic strategy to control host innate cell function²⁷. Ape-lin-4-5p showed several interesting potential gene targets: ARID3B involved in crucial cellular processes as transcriptional regulation²⁸, STARD13 involved in cell proliferation and tumor suppression and in particular, miR-125b induces metastasis by targeting STARD13 in in-vitro breast cancer cells²⁹. Additional gene targets are BMF, related to apoptosis activation³⁰ and FREM1, a gene related to IL-1 inflammatory response³¹.

Alignment of miRNA-125 and lin-4 from several helminths and human is reported in (a) (sja: S. japonicum; sma: S. mansoni; emu: E. multilocularis; fhe: F. hepatica, hco: H. contortus; hpo: H. polygyrus; asu: A. suum; ape: A. pegreffii present study; bma: B. malawi; hsa: H. sapiens; str: S. ratti). Complete identity is highlighted by an asterisk, while dashes indicate nucleotide deletions. The dendrogram of aligned miRNAs is shown in (b) and hairpin structure of ape-lin4 with the first nucleotide of matures sequence indicated with a red arrow is available in (c), according to RNA fold.

Among miRNAs selectively packaged into extracellular vesicles, novel-miR-19 showed no match with other helminths but a match with human miRNAs was observed (Table 4): its putative gene target is FAM83C, which is involved in regulation of MAPK signaling in cancer cells³².