Healthy, live, adult female Ascaris suum were collected from swine hosts from an abattoir in Marshalltown, Iowa, United States. These parasites were thoroughly washed multiple times in Ascaris Ringer’s Solution (ARS) [(13.14 mM NaCl, 9.67 mM CaCl₂, 7.83 mM MgCl₂, 12.09 mM Tris, 99.96 mM sodium acetate, 19.64 mM KCl) with gentamycin (100 μg/ml), ciprofloxacin hydrochloride (20 μg/ml), penicillin (10,000 units/ml), streptomycin (10,000 μg/ml), and amphotericin B (25 μg/ml) at pH 7.87 (all Sigma Aldrich, St Louis, MO)] and then incubated at 35°C. The following day, parasites were checked visually for signs of bacterial or fungal contamination, and worms were discarded if present. Worms tissue and EVs were collected on the second day of culture to allow for an overnight acclimation period and to limit negative impact on gene expression. Adult female parasites were chosen due to their large size, facilitating the ease of collection material for RNA extraction. To obtain tissue for RNA isolation, worms were cut along the ventral midline and the ovaries were gently removed for excision. Tissue was collected proximal to the bifurcation of the ovaries and rinsed with fresh ARS. Given the location of ovarian tissue collection, we acknowledge residual embryonated egg material may be present and circRNA data should be interpreted with that in mind. Body wall tissue was collected directly anterior to the genital aperture and musculature was scraped from the underlying cuticle with a single-edge razor blade. Approximately 200 mg of body wall muscle and ovary tissue samples were obtained in each sample isolation. Tissues were either used for immediate RNA extraction or stored at -80°C until use.

Individual worms were treated with 0.1 µM or 1 µM (final concentration) of ivermectin, diethylcarbamazine, or levamisole (all Sigma-Aldrich) for 24 h in 100 ml culture media in sterile 250 ml Erlenmeyer flasks. Drug concentrations were prepared from stock solutions dissolved in dimethyl sulfoxide (DMSO, Sigma-Aldrich). Conditioned media from drug treated and 0.1% DMSO vehicle control worms was collected after the 24-h time period and retained for downstream analysis. Body wall muscle and ovary-enriched tissue samples were collected from these parasites as described for immediate RNA extraction or storage at -80°C until use.

EVs were collected as previously described using differential ultracentrifugation (Zamanian et al., 2015; Harischandra et al., 2018; Loghry et al., 2020). Media was filtered through 0.2 µm PVDF vacuum filters (Sigma-Aldrich) and centrifuged at 120,000 x g for 90 min at 4°C. The supernatant was decanted, and pellets were filtered through a PVDF 0.2 µm syringe filter (GE Healthcare, Chicago, IL) and centrifuged further at 186,000 x g for 2 hours at 4°C. EV samples were then resuspended to 500 µl in dPBS (Thermo Fisher Scientific, Waltman MA) and stored at -80°C until use.

EV quantification and size determination were performed using nanoparticle tracking analysis (NTA; Nano-Sight LM10, Malvern Instruments, Malvern, United Kingdom). EV imaging was performed using transmission electron microscopy. A 2 µl aliquot of isolated EV preparation was placed onto a carbon film grid for 1 min. The drop was wicked to a thin film and 2 µl of uranyl acetate (2% w/v final concentration) was immediately applied for 30 s, wicked, and allowed to dry. Images were taken using a 200kV JEOL 2100 scanning and transmission electron microscope (Japan Electron Optics Laboratories, LLC, Peabody, MA) with a Gatan OneView camera (Gatan, Inc. Pleasanton, CA).

Total RNA was extracted from adult female A. suum body wall and ovary tissues using a two-step process. First, worm tissues were homogenized in TRIzol reagent (Life Technologies, Carlsbad, CA) according to manufacturer recommendations. Following RNA isolation in TRIzol, total RNA was further purified using the miRNeasy Mini kit (QIAGEN, Hilden, Germany), following the manufacturer’s protocol and assessed for purity and quantified using a NanoVue spectrophotometer (General Electric, Boston, MA). Similarly, total RNA was extracted from EV enriched samples isolated from conditioned media using the miRNeasy Micro kit (QIAGEN), again following the manufacturer’s protocol. Total RNA from EV isolation supernatants (i.e. EV-depleted media) was extracted using Zymo ZR urine RNA isolation kit (Zymo Research, Irvine, CA), following manufacturer’s instructions. All linear RNA was subsequently removed from each of these total RNA preparations by RNase R digestion (RNase R was provided by the Singh Laboratory, Iowa State University). 10 Units of RNase R was used for each reaction along with 2 µg total RNA. Reactions were incubated at 37°C for 45 min followed by heat inactivation at 65°C for 20 min circRNA was then stored at -80°C until use. circRNA samples were sent to LC Sciences (Houston TX, USA) for sequencing or transcribed into cDNA for qPCR validation.

CircRNA quality was assessed with a Bioanalyzer 2100 and RNA 6000 Nano LabChip Kit (Agilent, CA, USA), allowing a minimum RNA integrity number (RIN) of 7 (Schroeder et al., 2006) before fragmentation using NEBNext^® Magnesium RNA Fragmentation Module (NEB, Ipswich, MA) into short fragments using divalent cations under high temperature. The cleaved RNA fragments were then reverse-transcribed to create the cDNAs using SuperScript™ II Reverse Transcriptase (Thermo Fisher Scientific), which were next used to synthesize U-labeled second-stranded DNAs with E. coli DNA polymerase I (NEB), RNase H (NEB) and dUTP Solution (Thermo Fisher). An A-base was added to the blunt ends of each strand, preparing them for ligation to the indexed adapters. Each adapter contains a T-base overhang for ligating the adapter to the A-tailed fragmented DNA. Single- or dual-index adapters were ligated to the fragments, and size selection (300-600bp) performed with AMPureXP beads (Beckman Coulter, Brea, CA). After the heat-labile UDG enzyme (NEB) treatment of the U-labeled second-stranded DNAs, the ligated products were amplified with PCR by the following conditions: initial denaturation at 95°C for 3 min; eight cycles of denaturation at 98°C for 15 s, annealing at 60°C for 15 s, and extension at 72°C for 30 s; and then final extension at 72°C for 5 min. The average insert size for the final cDNA library was 300 ± 50 bp. Finally, 2 × 150bp paired-end sequencing (PE150) was performed on an Illumina Novaseq™ 6000 (Illumina) following the vendor’s recommended protocol.

Cutadapt (Martin, 2011) and custom perl scripts were used to remove adaptors, low quality bases and undetermined bases, followed by quality assessments with FastQC (Andrews, 2010). Bowtie2 (Langmead & Salzberg, 2012) and Tophat2 (Kim et al., 2013) were used to map reads to the genome of A. suum (Wang J. et al., 2017) (Accession number: PRJNA62057; AGO1), with remaining unmapped reads remapped to the genome using Tophat-fusion (Kim & Salzberg, 2011). Dual de novo assemblies of circular RNAs were performed with CIRCExplorer (Zhang et al., 2014; Zhang et al., 2016), one with Bowtie2 and Tophat2 mapped reads and another with Tophat-fusion back-spliced reads. Since samples were not prepared simultaneously, the distinct circRNA assemblies were concatenated and filtered for uniqueness with duplicated sequences being removed.

Sequenced RNA reads (SRR) (SRR15295818-SRR15295823) were aligned to circRNAs from both body wall and ovary-enriched tissues, and to the A. suum genome (PRJNA62057), to reduce bias by avoiding creation of mapping reads that are artifacts from flawed methodology (Wang L. et al., 2017) using Hisat 2.2.0 (Kim et al., 2019). Samtools 1.10 (Li et al., 2009) was used to convert sam alignments to bam alignment files to map RNA sequencing reads to the A. suum genome and the three circRNA samples. Mapping statistics were assessed using Picard 2.17.0 (Institute, 2019). Read counts were taken using featureCounts from the Subread package 1.6.0 (Liao et al., 2014). Differential expression was assessed using DESEQ2 1.20.0 (Love et al., 2014), with both unique and multiple mapping reads considered in separate comparisons. Differentially expressed circRNAs were subjected to GO and KEGG enrichment analyses using clusterProfiler (Yu et al., 2012) and Ontologizer (Bauer et al., 2008). The number of exonic, intergenic, and intronic circRNAs were compiled from sequencing data and visualized using GraphPad Prism version 9.2.0 (GraphPad Software Inc., San Diego, CA). Notes and scripts used to produce expression analysis are available at https://github.com/ISUgenomics/Kimber. Raw data and circRNA sequences can be viewed using bio-project number PRJNA750737 with SRA numbers, SRR15295818 - SRR15295823.

Validation that circRNAs identified using circRNA-seq are expressed in Ascaris tissue or EV enriched samples was performed using quantitative real-time PCR (qRT-PCR). RNase R treated RNA (generated as described above) was reversed-transcribed to cDNA with random hexamers using Invitrogen SuperScript III First Strand Synthesis System (Thermo Fisher Scientific) according to the manufacturer’s instructions. A total of 10 sets of divergent primers targeting distinct circRNAs were designed to back-splice junction sites (Supplementary Table S2) and qRT-PCR was performed using Power Up Sybr Green Master Mix according to the manufacturer’s protocol (Thermo Fisher Scientific). Conditions for qRT-PCR were as follows: 2 min at 50°C, 2 min at 95°C, then 40 cycles of 15 s at 95°C, 15 s at 55–60°C, and 1 min at 72°C. CircRNA abundance was quantified by extrapolating qRT-PCR C_T values from exogenous (spiked in) Homo sapiens actin alpha cardiac muscle one RNA (GenBank accession number: NM_005159). This exogenous spike in RNA was generated using MEGAScript T7 transcription kit (Thermo Fisher Scientific) following manufacturer’s protocols corresponding to 517 bp–1,165 bp of the published transcript (primer sequences are detailed in Supplementary Table S2). The exogenous spike in RNA standard curve analysis was created using Graph Pad Prism version 9.2.0 (GraphPad Software Inc.) using a nonlinear second-degree polynomial, least squares fit. circRNA expression X values were interpolated to quantify the concentration of circRNA in each tissue type and EVs.

Two complementary programs were used to predict miRNAs and circRNA interactions: miRanda (https://github.com/hobywan/miranda) and Targetscan (Agarwal et al., 2018) using the identified Ascaris circRNA sequences and miRNA datasets for human and A. suum downloaded from miRbase.org (Release 22.1) (Kozomara and Griffiths-Jones, 2011). Within these programs, a higher miRanda free energy (200–140) and lower Targetscan score (-0.13–0) were used to assign interaction confidence. The number of Ascaris and human miRNA interactions for each Ascaris circRNA was totaled and visualized using GraphPad Prism (GraphPad Software Inc.). Similarly, Ascaris and human miRNAs with high numbers of predicted Ascaris circRNA binding partners were also collated and visualized using GraphPad Prism (GraphPad Software Inc.).

Drug treated tissues and EV circRNA expression levels were calculated from qRT-PCR C_T values using 2^−ΔΔCq (Livak & Schmittgen, 2001). Following fold change analysis, the data was log (2) transformed using GraphPad Prism (GraphPad Software Inc.). To compare circRNA concentrations between treatments and samples, a two-way ANOVA with multiple comparisons (GraphPad Software Inc.) was used with a p-value less than 0.05 being considered significant. Each N number represented a new and biologically distinct batch of worms, with an individual no-treatment control for each batch.

To examine the presence of circRNA in A. suum tissues, a total of six independently prepared samples (three ovary-enriched, three body wall) were used for circRNA sequencing. After the removal of redundant or duplicated circRNAs, we identified 1,982 circRNAs in body wall tissue and 1,978 circRNAs in ovary-enriched tissue, for a total of 1,997 unique and distinct circRNAs (Figure 2A). There were a significant number of circRNAs shared between the two tissue types (1,963) and only 34 circRNAs were identified through circRNA-seq analysis as having tissue-specific expression: 15 circRNAs were identified only in the ovary-enriched samples and 19 were body wall specific. Tissue specific circRNAs are listed in Supplementary Table S1 If Ascaris circRNAs are functional, such a general spatial distribution pattern suggests either that the majority of A. suum circRNAs are linked to the regulation of transcriptional processes that are broadly conserved across different cell types, or that any functional specificity within tissues is driven by a more restricted temporal or spatial expression of interacting partners rather than the circRNAs themselves. Raw data and circRNA sequences can be viewed using bio-project number PRJNA750737 with SRA numbers, SRR15295818 - SRR15295823.

Of the total 1,997 circRNAs identified in ovary-enriched and body wall tissues, 1,178 (59%) were exonic, 779 (39%) were intronic and 40 (2%) were intergenic (Figure 2B). Whilst there is a lack of data distinguishing the functional relevance of exonic versus intronic circRNAs, there is the potential for exonic circRNAs to be translated into proteins (Legnini et al., 2017), and those translated proteins could have important biological roles.

The number of exons per circRNA was calculated and on average, circRNAs contained approximately three exons with 83% of circRNAs composed of multiple exons (two or more). The number of circRNAs per A. suum gene can be viewed in Figure 2C. Interestingly, 37% (752) of circRNAs seem to be derived from chromosome AgR001, while one specific gene locus, AgR002_g270, had the most derived circRNAs, a remarkable 15% of the total (303) (Figure 2C). AgR002_g270 does not have a known or annotated function, or any identified ortholog or paralogs associated with it, but blast analysis (https://blast.ncbi.nlm.nih.gov) of the AgR002_g270 coding sequencing returns ribosomal proteins from various nematode species, including A. lumbricoides, Ascaridia galli, Toxocara cati, and Baylisascaris procyonics (U.S. National Library of Medicine, 2022).

CircRNA have the potential to encode proteins (Legnini et al., 2017) so understanding the function of parental genes from which Ascaris circRNA derive could provide valuable information about circRNA function. GO and KEGG annotation analyses were conducted to predict possible functions of parental genes (Figure 3). Significant GO and KEGG terms were calculated by hypergeometric equation (Erdélyi et al., 1953) and terms with p-values less than 0.05 were defined as significant. Significant GO terms were divided into three groups, biological process, cellular component, and molecular function. 49 GO terms were involved in the biological process category for all three replicates (duplicates were excluded). The most enriched GO terms in the biological process category included “regulation of transcription; DNA templated” (GO:0006355), “transcription, DNA templated” (GO:0006351); and “phosphorylation” (GO:0016310) (Figures 3A–C). For cellular component, 41 individual GO terms were significantly enriched. They mainly consisted of “nucleus” (GO:0005634) and “cytoplasm” (GO:0005737). Molecular function had one enriched GO term, “nucleotide binding” (GO:0000166) (Figures 3A–C).

KEGG pathway analysis was also carried out to determine further significant pathways of circRNA parental genes and identify enriched pathways. There was a total of three different KEGG comparison groups, one for each of the ovary and body wall samples that were submitted from the same adult female worm. Of the three different comparisons that were analyzed for differential expression, one sample comparison did not have any significant differentially expressed KEGG terms and is not included in this analysis. In the two other comparisons, the most enriched KEGG pathways for differential expression between ovary-enriched and body wall tissues include PPAR signaling pathway, MAPK signaling pathway, and retrograde endocannabinoid signaling (Figures 3D,E). 10 of 28 differentially expressed KEGG pathways were involved in “signaling”, suggesting that circRNA parental genes are involved with signaling, signal transduction pathways and other important cellular processes that may support worm viability. If A. suum circRNAs are translated and yield functional proteins, this GO and KEGG term analysis points to possible functions that circRNAs could be performing in the worm, based on the parental gene.

We used qRT-PCR to confirm and validate the expression of six individual circRNAs identified in A. suum samples using circRNA-seq. The six initial circRNAs were selected due their high-count numbers from sequencing data in both ovary and body wall tissues. Divergent primers spanning back-splice junction sites were designed for each circRNA and can be viewed in Supplementary Table S2. Using an RNA exogenous spike in approach allowed us to calculate circRNA concentration levels using a standard curve, and as expected we did not observe any statistical significance in the abundance of individual circRNAs between the two tissue types (N = 5) (Figure 4A). This data validated the circRNA-seq approach as a means to broadly describe circRNA expression but our subsequent qRT-PCR analyses of differentially expressed circRNAs underscored the importance of verifying the circRNA-seq data with secondary methods (Figures 4B,C). The circRNA-seq datasets identified AgE14_g005_t01:37667–38355 as specifically expressed in the ovary-enriched samples and AgR024_g060_t03: 1256654–1259251 in the body wall samples. A full list of tissue specific circRNA expression for both ovary and body wall can be viewed in Supplementary Table S1. Although our qRT-PCR data confirmed that AgE14_g005_t01:37667–38355 was indeed localized to ovary tissue, AgR024_g060_t03: 1256654–1259251 was found to be expressed in both tissue types (Figure 4B). This result could be due to the ability of qRT-PCR to amplify partially degraded transcripts and has been observed in other species (Westholm et al., 2014; Zhao et al., 2019). Alternative explanations could be contamination of the ovary-enriched sample with body wall circRNA during dissection or from cross-contamination of RNA samples. It underscores the need for increased depth of circRNA sequencing married with additional validation measures to confirm spatial localization of circRNAs.

Further validation of select circRNA expression was performed, specifically, of AgR015: 277523–408464 and AgR001_g15_t01: 3040550–3048317. These circRNAs were prioritized due to their large size: AgR015: 277523–408464 was 130,941 nt long and AgR001_g15_t01: 3040550–3048317 was 7,767 nt long. The full length of these circRNAs was calculated from next generation sequencing data. We were able to confirm that these two RNA molecules are circRNAs through qPCR validation (Figure 4D) by creating primers specific to back-splice junction sites. This approach demonstrated both these RNAs form a circular structure and are not spurious background artiffacts or RNA molecules residual from RNase R digestion, even though they are larger in size than typical circRNAs (approximately 500–600 nt) (Ding et al., 2018).

circRNAs have been found to be secreted from mammalian parental cells in extracellular vesicles (EVs) (Lasda et al., 2016). Many species of nematodes are known to secrete EVs (Buck et al., 2014; Zamanian et al., 2015; Tzelos et al., 2016; Tritten et al., 2017; Eichenberger et al., 2018a; Eichenberger et al., 2018b; Harischandra et al., 2018; Shears et al., 2018; Hansen et al., 2019) but the presence of circRNAs in those vesicles has not been explored. We hypothesized that A. suum EVs would contain circRNAs. To test this hypothesis, we used qRT-PCR to quantify the abundance of select circRNA transcripts in A. suum EVs using the spike-in approach as previously described in section 3.3. We first isolated EVs from conditioned media, imaged the vesicles using TEM and performed nanoparticle tracking analysis (NTA) on the isolated samples to confirm EV morphology, size profile and concentration (Figures 5A,B). We tested these EVs for the presence of the same six circRNAs that were also tested for in tissue qRT-PCR validation (Figure 4A). Of these, two circRNAs (AgR007_g109_t01: 1908852–1911650 and AgB08X_g209_t01:3249044–3252560) could be consistently amplified from EV RNA samples (N = 3) (Figure 5C). These data indicate that parasitic nematode circRNAs are secreted into the host milieu in EVs. We also looked for the presence of circRNAs in A. suum supernatants generated through EV isolation, representing non-EV mediated mechanisms of circRNA secretion. This approach yielded insufficient amounts of total RNA to conduct RT-qPCR analysis suggesting EVs may represent the primary mechanism of circRNA secretion from this worm.

Given previous data published by our laboratory on the inhibitory effect of ivermectin (IVM) on parasitic nematode EV secretion, we examined whether IVM would inhibit circRNA secretion via EVs. Parasites were cultured in the presence or absence of 0.1 µM (Figure 5D) or 1 µM IVM (Supplementary Figure S1) to model a therapeutically relevant dose. After 24 h, parasite media was collected, total EV RNA extracted and used in RT-qPCR. When worms were treated with 0.1 µM or 1 µM IVM we did not observe any decrease in AgR007_g109_t01: 1908852-1911650 or AgB08X_g209_t01:3249044-3252560 abundance in isolated EVs (Figure 5D, Supplemental Figure 3, N = 3). This observation was perhaps surprising, given the strong and consistent evidence for an inhibitory effect of IVM on EV secretion in parasitic nematodes, including Ascaris (Harischandra et al., 2018; Loghry et al., 2020). This may point to other non-EV mediated routes of circRNA release from these worms. Other anthelmintic drugs are reported to have sporadic inhibitory effects on EV secretion by some life stages of filarial parasitic nematodes (Loghry et al., 2020). Therefore, we also looked at the effect of diethylcarbamazine (DEC) and levamisole (LEV) treatment on circRNA expression in A. suum EVs (Figure 5D, Supplementary Figure S1). Consistent with the IVM data, we did not see any inhibition in EV circRNA abundance, collectively indicating inhibition of circRNA secretion via EVs is not clearly associated with the mode of action of anthelmintic drugs.

To fully evaluate the effect of anthelmintic drug treatment on circRNA expression, we lastly examined whether IVM, DEC, or LEV had any effect on circRNA expression in Ascaris tissues using the same cohort of six prioritized circRNAs. Treatment of worms with 0.1 µM (Figure 6) or 1 µM IVM (Supplementary Figure S2) did not alter expression of any tested circRNA in ovary tissues (N = 4). Similarly, five of the six circRNAs in body wall tissue were unaffected by IVM treatment although AgR002_g269_t04:4591332–4600485 was downregulated by 95% compared to control (p = 0.0071, N = 4) in body wall tissue at 0.1 µM (Figure 6), but not at 1 µM (Supplementary Figure S2). Consistent with the results from IVM treated tissues, we did not observe any effect of DEC or LEV (Supplementary Figure S3) on circRNA expression in A. suum ovary or body tissues. Collectively our data do not support the hypothesis that anthelmintic drug mechanism of action involves a direct impingement of normal circRNA expression or secretion.

A well-established functional role for circRNAs is to regulate gene expression by binding miRNAs, effectively acting as miRNA sponges. circRNAs can contain numerous binding sites for individual or multiple miRNAs (Capel et al., 1993; Li F. et al., 2015; Zheng et al., 2016) For instance, murine CDR1as (ciRs-7) has 63 conserved binding sites for the miRNA mir-7 (Hansen T.B. et al., 2013), while circHIPK3 can sponge nine different human miRNAs (Zheng et al., 2016). Here, we wanted to probe potential miRNA interactions with the A. suum circRNA dataset to support the hypothesis that Ascaris circRNAs can act as miRNA sponges. We used the miRanda algorithm (Enright et al., 2003) to predict interactions between Ascaris circRNAs and endogenous A. suum miRNAs.

Approximately 10% of A. suum circRNAs (202 out of 1,997) were predicted to interact with A. suum miRNAs (Figure 7A). The number of miRNA interactions per circRNA varied, with AgR015:277522-408464 found to have highest number of distinct miRNA interactions 174) at 19 discrete binding sites on the circRNA molecule, illustrating that different miRNAs can bind to the same sites on circRNAs. Interestingly, for both human and Ascaris miRNAs, exonic circRNAs exhibited a higher number of predicted miRNA interactions than intronic circRNAs (Figure 7). The highest number of interactions between exonic circRNAs and Ascaris miRNAs was 174 (AgR015: 277523-408464), while 135 (AgB16X_g049_t01:1133282-1216863) interactions was observed to be the highest for intronic circRNAs. Interactions between exonic and intronic circRNAs with Ascaris and host miRNAs is summarized in Table 1.

To determine if the number of miRNA interactions per circRNA was simply a reflection of circRNA length, we normalized the miRNA interaction number to length of each circRNA (Figure 8A). There was not a strong linear correlation between length of A. suum circRNA length and A. suum miRNA interaction number (R ² = 0.613), suggesting that the number of A. suum miRNAs with which an A. suum circRNA interacts is not strongly correlated with the length of circRNAs and by extension, that some A. suum circRNAs are explicitly enriched in miRNA interaction sites. For example, AgR030:203455–226217 (14 KB) and AgR030_g080_t11: 1244134–1246599 (2 KB) are both enriched by a high number of miRNA interactions (Figure 8A).

In addition to examining the number of miRNA interactions for each A. suum circRNA, we also wanted to probe the interaction from the opposite direction and determine if any A. suum were specifically enriched in binding to A. suum circRNAs. In total, we observed that 180 distinct A. suum miRNAs were predicted to interact with A. suum circRNAs (Figure 8C). Two worm miRNAs had over 200 circRNA interactions, asu-miR-1822-5p and asu-miR-5348-5p (Figure 8C), but there are no known functions or phenotypes associated with these two miRNAs so postulating some functional relevance to this miRNA sponging activity is not possible. asu-miR-1822-5p did have the highest amount of circRNA interactions of all of Ascaris miRNAs at 394 total predicted interactions.

The secretion of A. suum circRNAs into the host environment in EVs and the potential for delivery of those circRNAs to host tissues seeded the possibility that secreted circRNAs could be acting as sponges for host (human) miRNAs. Therefore, we also used miRanda to predict interactions between Ascaris circRNAs and host (human) miRNAs. There were more predicted interactions between Ascaris circRNAs and human miRNAs, with a total of 398 distinct circRNAs interacting with host miRNAs (Figure 7B). AgR015:277522–408464 was predicted to have the highest number of miRNA interactions (2,308). This disparity may reflect the greater number of annotated miRNAs in the human genome compared to that of Ascaris. AgR015: 277523–408464 (exonic) and AgB16X_g049_t01:1133282–1216863 (intronic) were predicted to have the highest number of interactions across all miRNA-circRNA interactions, suggesting that these circRNAs contain a high number of binding sites. There is no supporting data in the circRNA literature, however, to support a functional difference between miRNA sponging by exonic versus intronic circRNAs.

There was a total of 577 circRNAs that were predicted to interact with A. suum miRNAs and 645 circRNAs that were predicted to interact with human miRNAs. While the total number of circRNA interactions for both A. suum and human miRNA were similar, the number of A. suum miRNAs interactions per individual circRNAs was almost 10 times fewer than the number of human miRNA interactions per circRNA (Figure 7), suggesting that A. suum circRNAs could have greater binding affinity for human miRNAs. 79% of circRNA–A. suum miRNA interactions had under 10 interactions per circRNA, while only 41% of circRNA–human miRNA interactions were under 10 interactions per circRNA. This suggests that while there is not a large difference in the number of parasite or human miRNAs that are predicted to bind Ascaris circRNAs, the number of miRNAs sponging to each circRNA is discrepant and A. suum circRNAs tend to have more binding sites for human miRNAs than A. suum miRNAs.

To determine if host miRNA binding was correlated with circRNA length, we also normalized circRNA length to the number of predicted host miRNA interactions (Figure 8B). In contrast to endogenous Ascaris circRNA-Ascaris miRNA interactions, there does seem to be a stronger linear correlation in this relationship (R ²: 0.932), suggesting that there is no specific enrichment of human miRNA binding to A. suum circRNAs.

Again, we observed a larger number of circRNA interactions occurring with human miRNAs as compared to A. suum miRNAs. Human miRNAs had a total of 2,414 predicted circRNA interactions with 16 miRNAs having over 200 circRNA interactions (Figure 8D). Relative to the Ascaris miRNAs that were predicted to frequently bind to Ascaris circRNAs, more functional information is known about these human miRNAs that are predicted to bind to Ascaris circRNAs, among them hsa-miR-4668-5p, which had the highest number of predicted circRNA interactions at 291. hsa-miR-4668-5p has been implicated in regulating TGF-beta signaling (Bhardwaj et al., 2020). TGF-beta is a cytokine involved in the proliferation, differentiation and function of lymphocytes, macrophages, and dendritic cells (Kubiczkova et al., 2012) and dysregulation of TGF-beta through circRNA sponging could potentially be another strategy these parasitic worms use to manipulate and modulate the host immune response. hsa-miR-6780b-5p had the second highest number of predicted circRNA interactions at 289. This miRNA has been linked to insulin resistance in hepatocellular carcinoma cells (HepG2) (Li et al., 2020) and whilst not directly related to immune system function, there could be other processes that could be altered as a result of its sequestering by circRNA sponging that could be advantageous to parasite infection. The recognition that a human immunomodulatory miRNA, hsa-miR-4668-5p, strongly interacts with secreted parasite circRNAs prompted us to examine whether any of the other human miRNAs that were predicted to bind to Ascaris circRNAs had known immunoregulatory functions. The percentage of miRNAs with documented immunoregulatory roles as a function of the total number of miRNAs predicted to bind parasite circRNAs was calculated (Table 2). Although this did not present a strong argument that immunomodulatory miRNAs are explicitly enriched for parasite circRNA binding, some interesting miRNAs were noted, including hsa-let-7, which had the highest percentage (0.5461) and is known to influence T-cell activation and mediates cytokine expression (Gilles and Slack, 2018).

To determine if circRNA could form a stable bond to worm and human miRNAs, we examined at the miRanda free energy score. A lower miRanda free energy score is associated with a more thermodynamically stable bond, and therefore, a stronger bond between circRNA and miRNA. We found that circRNAs form a more stable bond with human miRNAs as compared to Ascaris miRNAs based on these free-energy scores (Figure 8E). The free energy scores associated with human miRNAs are almost two times lower than Ascaris miRNA free energy scores suggesting the secreted parasite circRNAs could be forming tight bonds to host miRNAs and this could lead to changes in host gene expression through secretion of A. suum circRNAs in EVs.