Introduction
Pulmonate freshwater snails (Gastropoda: Mollusca) play significant roles in the transmission of parasitic helminths, most notably flukes (Trematoda: Platyhelminthes) of global medical and veterinary importance (Lu et al., Reference Lu, Gu, Limpanont, Song, Wu, Okanurak and Lv2018; Skala et al., Reference Skala, Walker and Horak2020). They serve as intermediate hosts supporting the asexual multiplication phase of the parasites’ life cycle. Given the significance of the diseases they transmit and their importance as a model in biological research, Biomphalaria, Bulinus and Lymnaea species are among the more studied freshwater snails. Biomphalaria and Bulinus snails transmit species of Schistosoma blood flukes (S. mansoni and S. haematobium, respectively) that are responsible for schistosomiasis, a public-health significant disease that causes serious intestinal, urogenital, hepatosplenic and neurological morbidities in over 250 million people across 78 countries (Buonfrate et al., Reference Buonfrate, Ferrari, Adegnika, Russell Stothard and Gobbi2025). On the other hand, L. stagnalis naturally transmits Trichobilharzia szidati, a zoonotic avian blood fluke that causes cercarial dermatitis (swimmer’s itch) in humans (Skala et al., Reference Skala, Walker and Horak2020). Outside their parasitological significance, these snails (especially L. stagnalis) are great models in translational neuroscience and ecotoxicology (Amorim et al., Reference Amorim, Abreu, Rodrigues, Peixoto, Pinheiro, Saraiva, Carvalho, Guimaraes and Oliva-teles2019; Rivi et al., Reference Rivi, Benatti, Colliva, Radighieri, Brunello, Tascedda and Blom2020; Souza-Silva et al., Reference Souza-Silva, De Souza, CaJ, Dos Santos Lima, Mol and Silveira2023).
While functional genomics elucidates the associations between the genetic makeup and phenotypic traits of an organism, RNA interference (RNAi) has been a powerful tool used in validating these associations. It silences the expression of specific genes to study their functions through observable loss-of-function mutations (Fire et al., Reference Fire, Xu, Montgomery, Kostas, Driver and Mello1998). Although RNAi is a natural gene regulatory system that is evolutionarily conserved among most eukaryotes, its mechanisms can be triggered artificially to suppress endogenous expression of target genes through in vivo or in vitro delivery of synthetic long double-stranded RNAs (dsRNAs, typically ≥100 nucleotides in length), short interfering RNAs (siRNAs, ∼22 nt) or microRNAs (miRNAs, ∼22 nt) that are complementary to the target gene sequence. When introduced into the cell, long dsRNAs are first processed into siRNA duplexes by Dicer (DCR) enzymes. The small silencing RNAs (referred to here as siRNAs and miRNAs) within the cell are then loaded into the RNA interference silencing complex (RISC), where they guide an Argonaute (AGO) protein to recognize and bind target mRNA. This binding induces cleavage via the siRNA pathway or translational repression via the miRNA pathway to directly target mRNA. Both pathways ultimately cause degradation (knockdown) of the target mRNA, preventing protein expression of the target sequence (Wilson and Doudna, Reference Wilson and Doudna2013).
RNAi has permitted the study of gene functions in many genetically intractable systems, including intermediate snail hosts of helminths, especially Biomphalaria and Lymnaea species (Korneev et al., Reference Korneev, Kemenes, Straub, Staras, Korneeva, Kemenes, Benjamin and O’shea2002; Jiang et al., Reference Jiang, Loker and Zhang2006; Smith et al., Reference Smith, Yadav, Fagunloye, Pels, Horton, Alsultan, Borns, Cousin, Dixon, Mann, Lee, Brindley, El-sayed, Bridger and Knight2021; Xu et al., Reference Xu, Y-w-q, Habib, Li, Yuan, Ke, Jiang, Dong and Zhao2023). Methodologically advanced applications of RNAi also demonstrate high potential as a biocontrol tool for insect pests and aquaculture diseases (Ghosh et al., Reference Ghosh, Hunter, Park and Gundersen-rindal2018; Charoonnart et al., Reference Charoonnart, Taunt, Yang, Webb, Robinson, Saksmerprome and Purton2023). The efficiency of RNAi, however, varies between organisms due to inter- and intra-specific differences in the essential effector proteins. These instances are evident among various groups of organisms, including protozoans (Baum et al., Reference Baum, Papenfuss, Mair, Janse, Vlachou, Waters, Cowman, Crabb and De Koning-ward2009), nematodes (Nuez and Felix, Reference Nuez and Felix2012), arachnids (Nganso et al., Reference Nganso, Sela and Soroker2020) and insects (Arraes et al., Reference Arraes, Martins-de-sa, Noriega Vasquez, Melo, Faheem, De Macedo, Morgante, Barbosa, Togawa, Moreira, Danchin and Grossi-de-sa2021). Therefore, to exploit the full potential of RNAi technology, species-specific profiling and functional screening of the pathway effector compositions are critical. Currently, in disease-transmitting freshwater snails, description of the RNAi effectors is fragmentary (Queiroz et al., Reference Queiroz, Silva, Jeremias, Baba, Caldeira, Coelho and Gomes2017; Huang et al., Reference Huang, Yoshitake, Asaduzzaman, Kinoshita, Watabe and Asakawa2021; Formaggioni et al., Reference Formaggioni, Cavalli, Hamada, Sakamoto, Plazzi and Passamonti2024). While the available studies provided some fundamental insights into snail RNAi biology, the compositional diversity of snail RNAi effector components has not been captured. Consequently, the limitations of snail RNAi techniques as tools for experimentation or disease control are not yet understood.
To provide basic insight into RNAi biology that will guide future functional genomics in freshwater snails that transmit flukes, this paper aims to determine the diversity of RNAi effector components and extrapolate the consequences of those components to RNAi competence in selected species of Biomphalaria (B. glabrata, B. pfeifferi), Bulinus (Bu. truncatus) and Lymnaea (L. stagnalis) snails, using bioinformatics, computational and gene expression approaches. While we identified core RNAi components that primarily mediate gene silencing in these snails, the absence of other important effectors known to enhance RNAi efficiency provides new mechanistic understanding of the snail RNAi pathways.
Materials and methods
Curation of query sequences and reciprocal BLAST searches
Known RNAi effectors from the well-studied Caenorhabditis elegans, Drosophila melanogaster and Homo sapiens models were identified from the literature in April 2024 and grouped based on their functional roles in the RNAi pathway. Other model organisms belonging to Order Rhabditida, Digenea, Coleoptera, Hemiptera, Orthoptera and Decapoda were included for further validation of the effector orthologues. The literature searches were conducted using Google Scholar and 3 scholarly databases (PubMed, Web of Science and Scopus), using general search terms such as ‘RNAi effectors Drosophila’, ‘RNAi effectors C. elegans’, ‘RNAi export C. elegans’ and more specific terms such as ‘sid-1 human’ with Boolean operators (AND, OR) where appropriate. Additional full-text articles were retrieved by manual review of the reference lists of relevant articles.
Query protein sequences of the identified effectors were retrieved from online resources, including WormBase for C. elegans and other selected nematodes (http://www.wormbase.org/) (Sternberg et al., Reference Sternberg, Auken, Wang, Wright, Yook, Zarowiecki, Arnaboldi, Becerra, Brown, Cain, Chan, Chen, Cho, Davis, Diamantakis, Dyer, Grigoriadis, Grove, Harris, Howe, Kishore, Lee, Longden, Luypaert, Müller, Nuin, Quinton-tulloch, Raciti, Schedl, Schindelman and Stein2024), WormBase ParaSite for platyhelminths (https://parasite.wormbase.org/index.html) (Howe et al., Reference Howe, Bolt, Shafie, Kersey and Berriman2017), FlyBase for D. melanogaster (https://flybase.org/) (Öztürk-çolak et al., Reference Öztürk-çolak, Marygold, Antonazzo, Attrill, Goutte-gattat, Jenkins, Matthews, Millburn, Dos Santos and Tabone2024), UniProt for H. sapiens (https://www.uniprot.org/) (The Uniprot Consortium, 2024) and NCBI for other sequences (https://www.ncbi.nlm.nih.gov/protein) (Goldfarb et al., Reference Goldfarb, Kodali, Pujar, Brover, Robbertse, Farrell, Oh, Astashyn, Ermolaeva, Haddad, Hlavina, Hoffman, Jackson, Joardar, Kristensen, Masterson, Mcgarvey, Mcveigh, Mozes, Murphy, Schafer, Souvorov, Spurrier, Strope, Sun, Vatsan, Wallin, Webb, Brister, Hatcher, Kimchi, Klimke, Marchler-bauer, Pruitt, Thibaud-nissen and Murphy2025). The protein sequences were used as search strings to identify orthologous RNAi proteins in the snails (Figure 1). Accessions of all query proteins used are provided in Supplementary Material 1. For the 4 snail species, genome assembly completeness was analysed using BUSCO (Benchmarking Universal Single-Copy Orthologs) v5.8.0 (Simao et al., Reference Simao, Waterhouse, Ioannidis, Kriventseva and Zdobnov2015) on Galaxy (https://usegalaxy.eu/). For B. glabrata, B. pfeifferi and Bu. truncatus, with annotated genome assemblies available through the NCBI GenBank database, reciprocal Basic Logic Alignment Search Tool BLASTp searches were performed against the NCBI non-redundant (nr) protein database v. 2.15.0, with default parameters. In cases where a protein orthologue was confirmed in only one of the snail species, the snail sequence was used as a query for other snail species. Since a structurally annotated L. stagnalis genome was not publicly available at the time of this research, a transcriptome was generated in-house by us using previously published data from other authors. The adapter-trimmed, high-quality RNA sequencing (RNA-Seq) reads (SRA Bioproject PRJNA420639) were mapped to an L. stagnalis genome (https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_900036025.1/), using HiSAT2 v2.1.0 (Kim et al., Reference Kim, Langmead and Salzberg2015), with default settings. Samtools v1.9 was used to convert SAM (Sequence Alignment Map) files to BAM (Binary Alignment Map) files (Li et al., Reference Li, Handsaker, Wysoker, Fennell, Ruan, Homer, Marth, Abecasis and Durbin2009). New Gene Transfer Format files were merged into a single matrix, and reads that uniquely mapped to a gene were counted using StringTie v1.3.6 (Pertea et al., Reference Pertea, Pertea, Antonescu, Chang, Mendell and Salzberg2015). Using ggffread, the resulting sequences were extracted into a FASTA file, which was used as a BLAST search database. The RNAi protein queries were searched using tBLASTn against the L. stagnalis transcriptome. BLAST searches were performed via the Northern Ireland High Performance Computing (NI-HPC) facility (https://www.ni-hpc.ac.uk/). Hits were converted to protein sequences using the ExPASy Translate tool (Artimo et al., Reference Artimo, Jonnalagedda, Arnold, Baratin, Csardi, De Castro, Duvaud, Flegel, Fortier, Gasteiger, Grosdidier, Hernandez, Ioannidis, Kuznetsov, Liechti, Moretti, Mostaguir, Redaschi, Rossier, Xenarios and Stockinger2012), followed by reciprocal BLASTp of the longest open reading frames on NCBI (Figure 1).
Methodological overview of the identification and gene expression analysis of RNAi effector components in snail vectors of parasitic flukes.
Figure 1 Long description
The diagram illustrates the RNA interference pathway analysis in snails, detailing various components and processes. Protein databases such as WormBase, UniProtKB, FlyBase and NCBI are listed. Main model organisms include C. elegans, D. melanogaster and H. sapiens. Other models include nematodes, coleopterans, insects, crustaceans and molluscs. Curated RNAi pathway protein queries are shown leading to FASTA files. Unannotated snail genomes like Lymnaea stagnalis and annotated snail genomes like Biomphalaria glabrata and Bulinus truncatus are depicted. Transcriptome, genome mapping, BLASTp and reciprocal BLASTp processes are indicated. Snail orthologues are identified, leading to FASTA files. Multiple sequence alignment tools such as Clustal Omega, Muscle 5, WebLogo and Jalview are used. Domain structure analysis involves InterPro, Phobius, HMMER, Pfam, SMART, AlphaFold 3, Biopython and ChimeraX. Phylogenetic tree construction uses MEGA 11 and FigTree. Data analysis and visualization are performed using R packages, qPCR, PCR, gel electrophoresis, cDNA and RNA isolation with TRIzol. Live Biomphalaria glabrata snails are shown at the bottom right.
Conserved domain/motif annotation and 3-dimensional structure prediction
For further identification of the snail orthologues and to resolve any ambiguity after applying the above parameters, conserved functional domains of selected snail orthologues were assessed by searching the sequences against multiple protein family databases including Pfam, NCBIfam, SMART, CDD, Prosite and Panther, using hmmscan algorithm of the hidden Markov models (Finn et al., Reference Finn, Clements, Arndt, Miller, Wheeler, Schreiber, Bateman and Eddy2015) and InterProScan (Blum et al., Reference Blum, Andreeva, Florentino, Chuguransky, Grego, Hobbs, Pinto, Orr, Paysan-lafosse, Ponamareva, Salazar, Bordin, Bork, Bridge, Colwell, Gough, Haft, Letunic, Llinares-lópez, Marchler-bauer, Meng-papaxanthos, Mi, Natale, Orengo, Pandurangan, Piovesan, Rivoire, Sigrist, Thanki, Thibaud-nissen, Thomas, Tosatto, Wu and Bateman2025), adopting default settings in all cases.
To predict effector function in dsRNA import, SID-1 (systemic RNAi-defective 1)-like orthologues were characterized in silico from cholesterol uptake protein 1 (CHUP-1), a SID-1 paralogue in C. elegans (Whangbo et al., Reference Whangbo, Weisman, Chae and Hunter2017). This was done by first assessing the presence or absence of the 2 mirroring transmembrane (TM) cholesterol-recognition motifs; the CRAC (Cholesterol-Recognition Amino acid Consensus) and CARC (inverse CRAC) motifs, characterized by L/V-X(1–5)-Y-X(1–5)-R/K and K/R-X(1–5)-Y/F-X(1–5)-L/V, respectively (Fantini et al., Reference Fantini, Di Scala, Evans, Williamson and Barrantes2016), followed by scrutiny of the N-terminal extracellular microdomains (Whangbo et al., Reference Whangbo, Weisman, Chae and Hunter2017). Topologies of the TM motifs were predicted using DeepTMHMM v1.0.42 (https://dtu.biolib.com/DeepTMHMM/) (Hallgren et al., Reference Hallgren, Tsirigos, Pedersen, Almagro Armenteros, Marcatili, Nielsen, Krogh and Winther2022) and Phobius (https://phobius.sbc.su.se/) (Kall et al., Reference Kall, Krogh and Sonnhammer2007) and were visualized using the Protter web server v1.0 (https://protter.ethz.ch/start/) (Omasits et al., Reference Omasits, Ahrens, Muller and Wollscheid2014).
For DCRs, multidomain analysis was performed by assessing the 2 RNase III catalytic domains and the tandemly arrayed C-terminal dsRBD (dsRNA-binding domain), the PAZ (Piwi–Argonaute–Zwille) and Platform domains, the N-terminal helicase domain and the Domain of Unknown Function DUF283 (Dicer dimerization domain) (Ciechanowska et al., Reference Ciechanowska, Pokornowska and Kurzynska-kokorniak2021; Zapletal et al., Reference Zapletal, Kubicek, Svoboda and Stefl2023). The 2 functionally distinct DCR proteins in D. melanogaster (DmDCR-1 and DmDCR-2) were included as models, and conservation of the pre-miRNA 5′-end binding pockets within the Platform–PAZ region (Park et al., Reference Park, Heo, Tian, Simanshu, Chang, Jee, Patel and Kim2011; Lee et al., Reference Lee, Lee, Kim, Kim and Roh2023) was analysed in the snail orthologues. For AGO orthologues, the catalytic triad DDX or tetrad DEDX (where X is either Asp (D) or His (H) amino acid) of the Piwi domain (Nakanishi et al., Reference Nakanishi, Weinberg, Bartel and Patel2012) was considered for predicting the slicer activity. Staufen (STAU-1/2) proteins were assessed by their characteristic dsRBD repeats (Yoon et al., Reference Yoon, Mogilicherla, Gurusamy, Chen, Chereddy and Palli2018), and the presence/absence of an N-terminal signal peptide (Koo and Palli, Reference Koo and Palli2024) was predicted using the SignalP 6.0 model (Teufel et al., Reference Teufel, Almagro Armenteros, Johansen, Gislason, Pihl, Tsirigos, Winther, Brunak, Von Heijne and Nielsen2022).
AlphaFold3 Server (Abramson et al., Reference Abramson, Adler, Dunger, Evans, Green, Pritzel, Ronneberger, Willmore, Ballard, Bambrick, Bodenstein, Evans, Hung, O’neill, Reiman, Tunyasuvunakool, Wu, Žemgulytė, Arvaniti, Beattie, Bertolli, Bridgland, Cherepanov, Congreve, Cowen-rivers, Cowie, Figurnov, Fuchs, Gladman, Jain, Khan, Low, Perlin, Potapenko, Savy, Singh, Stecula, Thillaisundaram, Tong, Yakneen, Zhong, Zielinski, Žídek, Bapst, Kohli, Jaderberg, Hassabis and Jumper2024) was used to predict 3-dimensional protein fold structures. AlphaFold files were converted to Protein Data Bank files with predicted local distance difference test (pLDDT) confidence scores assigned to the B-factor field using Biopython (v1.85) (Cock et al., Reference Cock, Antao, Chang, Chapman, Cox, Dalke, Friedberg, Hamelryck, Kauff, Wilczynski and De Hoon2009), before further analysis and visualization using ChimeraX v1.9 (Meng et al., Reference Meng, Goddard, Pettersen, Couch, Pearson, Morris and Ferrin2023) Figure 1. Structural regions with very low prediction confidence (pLDDT < 50) were trimmed out.
Multiple sequence alignments and phylogenetic analysis
Alignments of multiple protein sequences were performed using Clustal Omega v1.2.4 (Sievers and Higgins, Reference Sievers and Higgins2018) and visualized in Jalview v2.11.4.1 (Procter et al., Reference Procter, Carstairs, Soares, Mourao, Ofoegbu, Barton, Lui, Menard, Sherstnev, Roldan-martinez, Duce, Martin and Barton2021). WebLogo v3.7 (Crooks et al., Reference Crooks, Hon, Chandonia and Brenner2004) was used to generate graphical representations of patterns within some multiple sequence alignments (Figure 1). Phylogenetic trees were constructed using MEGA11 software v11.0.13 (Tamura et al., Reference Tamura, Stecher and Kumar2021). Multiple full-length or partial protein sequences were subjected to the plug-in MUSCLE sequence alignment programme, and evolutionary relationships were inferred using maximum likelihood (ML) clustering with 1000 non-parametric bootstrap replicates. The non-bilaterian Nematostella vectensis (sea anemone) and Amphimedon queenslandica (sponge) were used as outgroups. Best-fitting substitution models for ML clustering were selected based on the lowest Bayesian information criterion and corrected Akaike information criterion scores. To account for sequence divergence and composition across alignment matrices, an 80% site coverage cut-off was applied to manage gaps uniformly during both the best-fitting substitution model selection and the final ML tree construction. Bootstrap consensus trees were visualized and edited with FigTree v1.4.4 (http://tree.bio.ed.ac.uk/).
Protein–protein interaction (PPI) network and Gene Ontology (GO) analysis
Putative B. glabrata RNAi effector protein sequences were mapped to the STRING v12.0 database (Szklarczyk et al., Reference Szklarczyk, Nastou, Koutrouli, Kirsch, Mehryary, Hachilif, Hu, Peluso, Huang, Fang, Doncheva, Pyysalo, Bork, Jensen and Von Mering2025) to generate interaction networks and GO terms. PPI networks were built with 70% confidence and clustered using the k-means algorithm. Clustered networks were loaded into Cytoscape v3.10.3 (Shannon et al., Reference Shannon, Markiel, Ozier, Baliga, Wang, Ramage, Amin, Schwikowski and Ideker2003) for further annotation, analysis and visualization. GO enrichments for biological process (BP), cellular component (CC) and molecular function (MF) were obtained with enrichment false discovery rate (FDR) ≤ 0.05 and P-values corrected for multiple testing for each category using the Benjamini–Hochberg procedure (Benjamini and Hochberg, Reference Benjamini and Hochberg1995).
Gene expression analysis
Snail husbandry
Biomphalaria glabrata snails (uninfected NMRI strains, obtained from the Wellcome Sanger Institute, Hinxton, UK, and maintained for over 4 years) were used in this study. The snails were separated from the main aquarium and reared in aerated, transparent 1 L plastic containers filled with diluted 1× artificial river water (10× ARW contained 5 L H2O, 2.78 g CaCl2, 6.15 g MgSO4·7H2O, 2.1 g NaHCO3, 215 mg K2SO4 and 250 µL FeCl3·6H2O, at pH 7.0) and kept at 25 ℃, under 12:12 light/dark photoperiod. Snail density was kept at a maximum of 10 adults per 1 L ARW. The snails were fed tropical fish food flakes (Pets at Home, UK). Both the rearing water and the food were changed 3 times weekly.
Experimental design
Gene transcripts of 6 putative RNAi effectors in B. glabrata (Ago2, Dcr1, Eri1, Sid1L, Chc1 and Stau) were selected for mRNA expression evaluation, based on their primary relevance in the RNAi machinery. Expression levels of key effectors such as Dicers and Argonautes may vary significantly with age (Sanz-ros et al., Reference Sanz-ros, Mas-bargues, Romero-garcia, Huete-acevedo, Dromant and Borras2023). To examine whether snail size influences expression levels of 3 major genes, Ago2, Dcr1 and Eri1, whole snail tissues of juveniles and adult snails, which were further categorized into 4 subgroups based on shell size (diameter), were analysed. The first juvenile group (J1) was 3–4 mm in size, the second juvenile group (J2) was 6–7 mm in size, the first adult group (A1) was 10–11 mm in size and the second adult group (A2) was 13 mm in size. Each group had 5 biological replicates. In another assay, gene expression profiles of Sid1L, Chc1 and Stau were analysed in various isolated tissues rather than the whole soft body. Five adult snails (13 mm) were sacrificed by shell crushing and removal and were anatomically dissected to isolate the head-foot (HF) region, ovotestis (OV), hepatopancreas (digestive gland) (HP/DG) and albumen gland (AG) from each snail. The remaining tissue, here termed ‘trunk’ T (consisting of the mantle lobe, mantle cavity, pseudobranch, intestine, stomach, heart and kidney), was also assessed for gene expression (Figure 2A).
Gene expression analysis of selected putative RNAi genes in Biomphalaria glabrata. (A) Anatomy of the whole Biomphalaria glabrata soft body. HF, head-foot region; ML, mantle lobe; H, heart; AG, albumen gland; S, stomach; HP/DG, hepatopancreas/digestive gland; Ov, ovotestis; trunk (T, a multi-organ section containing the mantle lobe, lung, pseudobranch, intestine, stomach, heart and kidney). (B) Agarose gel image of real-time PCR products. (C) Relative mRNA expression level of Argonaute 2, Ago2, (D) Dicer 1, Dcr1, and (E) enhancer of RNA interference, Eri1, in B. glabrata juvenile and adult snails. Comparisons of relative mRNA expression of (F) clathrin heavy chain 1, Chc1, (G) systemic RNA interference defective protein 1-like, Sid1L, (H) Staufen, Stau1/2, and (I) Chc1 vs Sid1L, in B. glabrata HF, OV, AG, HP/DG and T tissue lysates. *P ≤ 0.05, **P ≤ 0.01 and ***P ≤ 0.001.
Figure 2 Long description
A snail outline diagram with internal labels HF, ML, H, AG, S, Hp Dg, Ov and T. A scale label reads 2 millimeter. A vertical label at the right reads Head foot. The image B showing a gel band strip with multiple lanes. Lane labels read M, HF, AG, Ov, Hp Dg, T, Ago2, Dcr1, Eri1, Sid1L and Stau. The image C showing a bar chart titled Ago2. The x-axis label is Size group, with categories J1, J2, A1 and A2. The y-axis label is Relative expression, ranging from 0 to 0.05. Four bars are shown with overlaid points and error bars. The bar for A2 is the tallest and the bar for J2 is the shortest. The image D showing a bar chart titled Dcr1. The x-axis label is Size group, with categories J1, J2, A1 and A2. The y-axis label is Relative expression, ranging from 0 to 0.08. Four bars are shown with overlaid points and error bars. The bar for J1 is the tallest and the bars for A1 and A2 are shorter than J1. The image E showing a bar chart titled Eri1. The x-axis label is Size group, with categories J1, J2, A1 and A2. The y-axis label is Relative expression, ranging from 0 to 0.004. Four bars are shown with overlaid points and error bars. The bar for J1 is the tallest and the bars for A1 and A2 are near the lower end of the axis. The image F showing a box plot titled Chc1. The x-axis categories are AG, HF, Hp Dg, Ov and T. The y-axis label is Relative expression, ranging from 0 to 0.30. Five box plots are shown with overlaid points. The Ov category has the highest box plot values and AG has the lowest. The image G showing a box plot titled Sid1L. The x-axis categories are AG, HF, Hp Dg, Ov and T. The y-axis label is Relative expression, ranging from 0 to 0.20. Five box plots are shown with overlaid points. The Ov category has the highest box plot values and AG, HF, Hp Dg and T are clustered near the lower end of the axis. The image H showing a box plot titled Stau1 2. The x-axis categories are AG, HF, Hp Dg, Ov and T. The y-axis label is Relative expression, ranging from 0 to 1.00. Five box plots are shown with overlaid points. The Hp Dg category has the highest box plot values and Ov is near the lower end of the axis. The image I showing a grouped bar chart titled Chc1 vs Sid1L. The x-axis categories are AG, HF, Hp Dg, Ov and T. The y-axis label is Relative expression, ranging from 0 to 0.20. Two bars per category are shown, with a legend at right labeled Gene and entries Chc1 and Sid1L. Across AG, HF, Hp Dg, Ov and T, the Chc1 bars are higher than the Sid1L bars, with the largest Chc1 bar at Ov.
RNA isolation, cDNA synthesis and quantitative PCR
Total RNA was extracted from snail tissues using TRIzol (Invitrogen), followed by DNase treatment using DNA-free™ kit (Invitrogen). RNA samples were quantified spectrophotometrically before and after DNase treatment, using a DeNovix DS-11 FX instrument. First-strand cDNA was synthesized by reverse transcription (RT) of 1.8 µg of each DNase-treated RNA sample using High-Capacity RNA-to-cDNA™ Kit (Applied Biosystems). The RT-PCR conditions were 37 ℃ for 1 h, 95 ℃ for 5 min and 4 ℃ until storage. For quantitative PCR (qPCR), oligonucleotide primers (Metabion) were designed using Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) and Primer3Plus (https://www.primer3plus.com/). Supplementary Material 1 contains the list of all primer sequences used. All primer sets were tested on endpoint PCR with 1.5 µL of cDNA template and electrophoresed on a 2% agarose gel in 1× Tris-acetate-EDTA buffer to confirm amplicon size and evaluate amplification specificity. Gene expression levels in snail cDNA samples were subsequently assessed by real-time qPCR. Each qPCR reaction was set up in 3 technical replicates, each having a total reaction volume of 10 µL, containing 1.5 µL of cDNA, 5.0 µL of SensiFAST SYBR No-ROX mix 2× (Meridian Bioscience), 0.5 µL of the F primer, 0.5 µL of the R primer and 2.5 µL of molecular-grade water. Negative no-template controls were included for each sample by replacing cDNA with molecular-grade water. qPCR was run on a Qiagen Rotor-Gene Q instrument by 40 cycles of 95 ℃ for 10 s, 60 ℃ for 15 s and 72 ℃ for 20 s, after an initial hold for 10 min at 95 ℃. Thermal cycling results were processed on Rotor-Gene Q series software v2.3.1. Biomphalaria glabrata β-Actin gene was selected as an endogenous control for normalization, as it has shown stable expression in B. glabrata of different ages and sizes (Xu et al., Reference Xu, Y-w-q, Habib, Li, Yuan, Ke, Jiang, Dong and Zhao2023). Relative gene expression was calculated using Pfaffl’s E −∆Ct method (Pfaffl, Reference Pfaffl2001).
qPCR data analysis and visualization
All qPCR data were statistically analysed, compared among groups and visualized using R software v4.3.2 (R Core Team, 2025). To ensure that replicate-level variation did not bias downstream statistical analyses, we screened all biological replicates for outliers using a median absolute deviation (MAD) method (Leys et al., Reference Leys, Ley, Klein, Bernard and Licata2013), applied within each group. Replicates deviating more than 3 times the MAD from the group median were flagged as potential outliers and removed. In each group, at least 3 replicates were inliers and retained for statistical analysis. Outlier detection was performed in R using base functions and the dplyr package v1.1.4, prior to normality and variance testing. Data that assumed normal distribution (Shapiro–Wilk’s P > 0.05) were subjected to 1-way ANOVA for multiple comparisons, followed by Tukey’s HSD test. Otherwise, the non-parametric Kruskal–Wallis test was performed for multiple comparisons among groups, and the Wilcoxon rank-sum (Mann–Whitney U) test was performed for comparison between 2 independent variables. Non-parametric multiple comparisons were followed up by Dunn’s post hoc test with Bonferroni correction. In all cases, P-values ≤ 0.05 were considered statistically significant.
Results
Comparative BUSCO analyses showed uneven genome completeness between snail species
Based on BUSCO estimates, B. glabrata, B. pfeifferi and Bu. truncatus genomes showed high completeness, although higher gene duplication was found in Bu. truncatus (Table 1). Under the same parameters, the L. stagnalis genome showed much lower BUSCO scores for completeness. Relative to other genomes, many more pan-eukaryotic genes may be fragmented or missing in the L. stagnalis genome (Table 1).
Comparative genome completeness estimates for the snail species based on BUSCO analysis
Table 1 Long description
Genome completeness is reported for four snail species using BUSCO categories, with columns for Species, GenBank, C percent, S percent, D percent, F percent, M percent, and n. Biomphalaria glabrata, GenBank, GCA_947242115.1, has C 96.9 percent, S 95.7 percent, D 1.2 percent, F 2.4 percent, M 0.8 percent, n 255. Biomphalaria pfeifferi, GenBank, GCA_030265305.1, has C 96.1 percent, S 95.7 percent, D 0.4 percent, F 2.7 percent, M 1.2 percent, n 255. Bulinus truncatus, GenBank, GCA_021962125.1, has C 96.5 percent, S 67.5 percent, D 29.0 percent, F 2.7 percent, M 0.8 percent, n 255. Lymnaea stagnalis, GenBank, GCA_900036025.1, has C 73.3 percent, S 71.8 percent, D 1.6 percent, F 21.2 percent, M 5.5 percent, n 255. Across species, C values are near 96 percent for the two Biomphalaria species and Bulinus truncatus, while Lymnaea stagnalis has lower C and higher F and M.
BUSCO (Benchmarking Universal Single-Copy Orthologs) categories were reported as: C (Complete), subdivided into S (Complete and single-copy) and D (Complete and duplicated); F (Fragmented) and M (Missing). n indicates the total number of BUSCO groups searched. All assemblies were evaluated under the same parameters (eukaryota_odb10 lineage, euk_genome_met mode and metaeuk gene predictor) to ensure comparability.
Snail vector genomes encode RNAi pathway components
A total of 115 RNAi effector proteins were collated, comprising 47 uptake and transport proteins, 11 DCR and DCR-associated proteins, 4 AGOs (excluding worm-specific AGOs), 11 RISC proteins, 7 RNAi amplification effectors, 15 RNAi inhibitors and 20 nuclear RNAi effectors (Table 2). Out of these, 83 orthologues were found in L. stagnalis, 76 in B. glabrata, 75 in Bu. truncatus and 74 in B. pfeifferi (Table 2, Figure 3). Each snail species possesses at least 1 orthologue in each functional group of the RNAi effectors (Table 2). Nevertheless, all the snails examined lack orthologues of the RNA-dependent RNA polymerases (RdRPs) (RRF-1 and EGO-1) and associated proteins (RDE-10 and RDE-11) (Figure 3), which are key players in secondary RNAi amplification. Similarly, core components responsible for heritable or transgenerational RNAi in C. elegans (the nuclear RNAi defective NRDE-1, NRDE-2, NRDE-3 and NRDE-4) were absent from the snails (Figure 3). While these absences may limit RNAi potency, other preserved effector components are sufficient for successful RNAi.
A compendium of predicted RNAi effector components in Biomphalaria (B. gla, B. pfe), Bulinus (B. tru) and Lymnaea (L. sta) snails. Query sequences used were retrieved from C. elegans (c), D. melanogaster (d) or H. sapiens (h). In Biomphalaria spp., WHT-1 and WHT-3 protein orthologues are encoded by a single gene. Human TRBP-2 and PACT shared a single orthologous protein in B. glabrata and Bu. truncatus. In L. stagnalis, Drh-1 and Drh-3 returned separate sequences. Full accessions and sequences are available in Supplementary Material 1.
Figure 3 Long description
RNAi effector proteins across four snail species: B. glabrata, B. pfeifferi, Bu. truncatus and L. stagnalis. Columns are organized into categories: RNAi uptake and transport (Classical proteins), RNAi uptake and transport (Vesicle-trafficking proteins), Dicer proteins, RNAi inhibitor proteins, RNAi amplification proteins, Argonaute proteins, RISC proteins and Nuclear RNAi effector proteins. Each protein row has four markers indicating presence (filled circle) or absence (empty circle) for each species. Key observations include broad conservation in RNAi uptake and transport proteins, while RNAi amplification proteins show notable absences. Dicer and RISC proteins are mostly present across species, with slight variations. The table highlights differences in nuclear RNAi effector proteins, where L. stagnalis shows more presence compared to others. This structured summary provides an overview of protein distribution without listing each protein individually.
Functional grouping and diversity of snail RNAi effector proteins
Table 2 Long description
Counts are listed by effector protein category, with columns for No. queries and four species: B. glabrata, B. pfeifferi, Bu. truncatus, and L. stagnalis. Dicer and DCR-associated has 11 queries, with species counts of 8, 6, 8, and 9. Argonautes (AGOs) has 4 queries, with 1 in each species. RISC has 11 queries, with 9 in each species. RNAi amplification has 7 queries, with 1 in B. glabrata, 1 in B. pfeifferi, 1 in Bu. truncatus, and 3 in L. stagnalis. Nuclear RNAi has 20 queries, with 11, 11, 10, and 12. RNAi inhibition has 15 queries, with 6, 6, 6, and 9. Total is 115 queries, with species totals of 76 for B. glabrata, 74 for B. pfeifferi, 75 for Bu. truncatus, and 83 for L. stagnalis.
RNAi uptake and spreading may be mediated by a vesicle-trafficking pathway, rather than the TM SID-dependent pathway
Reciprocal BLASTp screens indicated that snail SID-1-like sequences displayed higher sequence similarities with C. elegans CHUP-1 and human SID-1-like orthologue (HsSIDT-1) (Figure 4A). Although B. truncatus SID-1 shared slightly higher per cent sequence identity with the C. elegans SID-1 (CeSID-1) orthologue (Figure 4A), the best match (higher bit score of 169 with 81% query cover and a lower E-value of 2e − 47) was found between the B. truncatus SID-1-like orthologue and C. elegans CHUP-1. Further TM sequence analysis identified 6 putative cholesterol-binding CRAC motifs in each of the 2 Biomphalaria species, and 3 CRAC motifs in Bu. truncatus and L. stagnalis (Figure 4B; Table 3). Each snail sequence also possessed multiple mirrored TM CARC motifs (Table 3). Conserved CRAC and CARC motifs were also found across other molluscan SID-1-like proteins, including slugs, limpets, bivalves and cephalopods (Table 3). Interestingly, single CRAC (TM1) and CARC (TM9) motifs were found within the CeSID-1 sequence (Figure 4C; Table 3), suggesting that the structural framework for cholesterol binding is conserved in both SID-1 and CHUP-1 proteins. However, the characteristic N-terminal microdomains of CeSID-1 were not observed in any of the molluscan or non-molluscan SID-1-like orthologues examined (Supplementary Material 1). Similar to SID-1, orthologues of other C. elegans dsRNA-transporting SID proteins (SID-2 and SID-5) appeared absent from snail genomes (Figure 3). These results suggest that the canonical SID-dependent RNAi transport pathway, well characterized in C. elegans, may be absent from gastropod lineages.
Characterization of snail SID-1-like proteins. (A) Heatmap matrix for reciprocal BLASTp showing per cent identity between orthologues. (B) Predicted topology of Bu. truncatus SID-1 showing the TM and non-TM CRAC/CARC motifs. (C) Predicted topology of C. elegans SID-1 secondary structure showing single TM CRAC and CARC motifs. The first 19 amino acid residues at the N-terminus of C. elegans SID-1 protein (M1–N19) were predicted to constitute a signal peptide domain. In figures (B) and (C), numeric annotations in orange colour indicate the number of TM domains. CRAC, cholesterol-recognition amino acid consensus; CARC, inverse CRAC motif; IPML, inner plasma membrane layer; OPML, outer plasma membrane layer; PM, plasma membrane; TM, transmembrane.
Figure 4 Long description
Three parts are shown. The first part is a heatmap matrix labeled A, showing percent identity based on reciprocal BlastP between various SID-1-like proteins. The matrix includes rows and columns labeled with different protein names such as Cele SID-1, Cele CHUP-1, Hsap SIDT-1 and others, with numerical values indicating percent identity. The second part, labeled B, depicts the predicted topology of SID-1 with transmembrane and non-transmembrane motifs. It includes TM CARC motifs in blue, TM CRAC motifs in red, non-TM CRAC motifs in green and non-TM CARC motifs in black. The diagram shows the extracellular, intracellular and plasma membrane layers, with numeric annotations indicating transmembrane domains. The third part, labeled C, illustrates a similar topology with single TM CRAC and CARC motifs, also showing extracellular, intracellular and plasma membrane layers with numeric annotations for transmembrane domains.
Amino acid residues of the putative transmembrane CRAC and CARC motifs in selected SID/CHUP-1-like sequences
Table 3 Long description
Protein rows list residue ranges for two columns, transmembrane CRAC motifs and transmembrane CARC motifs. CeleSID 1 has CRAC L674 to K681 and CARC K308 to L320. CeleCHUP 1 has CRAC L440 to R448 and CARC K265 to V274, K370 to L378, and K636 to L643. BglaSID 1 lists six CRAC ranges from L528 to R534 through L786 to K798, and three CARC ranges K505 to L514, K667 to V677, and K738 to L750; BpfeSID 1 lists six CRAC ranges from L477 to R483 through L735 to K745, and three CARC ranges K454 to L463, K616 to V626, and K687 to L699. BtruSID 1 lists three CRAC ranges V55 to R65, L134 to K144, and V247 to R258, and two CARC ranges K15 to V25 and R86 to L98; LstaSID 1 lists three CRAC ranges L61 to R70, V117 to R127, and L196 to K206, and two CARC ranges K77 to V87 and R148 to L160. AcalSID 1 lists five CRAC ranges from L579 to K585 through L837 to K847, and three CARC ranges K556 to L565, K718 to V728, and R789 to V800. CvirSID 1 lists four CRAC ranges L421 to K427, L468 to R476, L543 to R552, and L560 to R569, and four CARC ranges K273 to V279, R398 to L407, K559 to V570, and R630 to L642. LgigSID 1 lists three CRAC ranges L56 to K62, L173 to R182, and L309 to K318, and two CARC ranges K33 to L42 and R260 to L272; MgalSID 1 lists two CRAC ranges L488 to K494 and L610 to R619, and three CARC ranges K265 to V271, K465 to L474, and K626 to V637. OvulSID 1 lists four CRAC ranges L517 to K523, L639 to R648, V695 to R705, and L774 to K784, and two CARC ranges K494 to L503 and K729 to L738; HsapSIDT 1 lists three CRAC ranges L584 to R593, L643 to R652, and L722 to K732, and four CARC ranges K302 to L311, K540 to L552, R677 to L686, and R707 to L714.
CRAC, cholesterol-recognition amino acid consensus; CARC, inverse CRAC motif; K, lysine; L, leucine; R, arginine; V, valine; Cele, Caenorhabditis elegans (nematode); Bgla, Biomphalaria glabrata (gastropod), Bpfe, Biomphalaria pfeifferi (gastropod); Btru, Bulinus truncatus (gastropod); Lsta, Lymnaea stagnalis (gastropod); Acal, Aplysia californica (gastropod slug); Cvir, Crassostrea virginica (bivalve); Lgig, Lottia gigantea (limpet); Mgal, Mytilus galloprovincialis (bivalve); Ovul, Octopus vulgaris (cephalopod); Hsap, Homo sapiens (mammal). Supplementary Material 1 contains the multiple sequence alignments of the transmembrane CRAC/CARC motifs and the accession numbers of all protein sequences.
Protein components of an alternative, SID-independent, RNAi-importing clathrin-mediated endocytosis (CME) pathway were expectedly preserved in the intermediate host snails (Figure 3). These were reflected in clathrin heavy chain 1 (CHC-1) and associated proteins, such as clathrin adapter protein AP-50/AP-2, clathrin interactor CLINT-1 (the mammalian ENTH [epsin N-terminal homology]-domain-containing orthologue of C. elegans RSD-3), dynamin-1 (DYN-1), vacuolar ATPase VHA-16, VhaSFD (a VHA-16 subunit), RAB-7 (a Ras-like GTPase), ARF72A/ARL-1 (an ADP-ribosylation factor), a non-receptor tyrosine kinase ACK (the mammalian orthologue of C. elegans SID-3), scavenger receptors (Sr-CI and Eater) and protein chaperones (HSC70-3 and HSC70-4). Snail orthologues corresponding to vesicle proteins that facilitate cellular export of silencing RNAs, REXD-1, TBC-3, EHBP-1 and RAB-10 were also present (Figure 3). These results suggest that the import of exogenous silencing RNAs and their systemic spreading in vector snails may be mediated majorly via the clathrin-dependent vesicle-trafficking pathway, given the absence of SID-1-like sequences in our datasets.
Snail Dicer (DCR) and Argonaute (AGO) proteins contain functional domains and amino acid residues for miRNA/siRNA biogenesis and processing
Sequences similar to DCR and AGO proteins from humans, D. melanogaster and C. elegans were identified in the snails (Figure 3). Single orthologues of DCR proteins were found in all snails, except Bu. Truncatus, which returned 2 putative isoforms (Figure 5A, Supplementary Material 1). Two putative AGO isoforms were found in B. glabrata (Figure 6A, Supplementary Material 1).
Comparative domain structures of Dicer in H. sapiens, C. elegans, D. melanogaster and the studied snail species. (A) Cartoon representation of the domains with their amino acid sequence ranges. The ‘head’ consists of the Platform–PAZ domain, the ‘base’ is made up of the helicase domain, while the ribonuclease domain forms the ‘core’. InterProScan and Pfam accessions: DExH, DEAD box helicase domain (IPR011545/PF00270); PBD, Dicer partner-binding domain (IPR048513/PF20930); Helicase C, helicase conserved C-terminal domain (IPR001650/PF00271); Dimer, Dicer dimerization domain (IPR005034/PF03368); PAZ (Piwi–Argonaute–Zwille) domain (IPR003100/PF02170); RNase III, Ribonuclease III (IPR000999/PF00636); dsRBD, dsRNA-binding domain (IPR014720/PF20932). AlphaFold 3-dimensional (3D) protein folding architectures, showing the domains in same colour with cartoon representation, for H. sapiens (model), B. glabrata, Bu. truncatus and L. stagnalis Dicers are presented on either side. (B) Hypothetical motifs for 5′-end recognition pockets in the Platform–PAZ region of the studied snail species. BglaDCR (R822, R824, R855, R1030 and R1037); BpfeDCR (R763, R765, R796, R971 and R978); BtruDCR-a (R804, R806, R837 and R945); BtruDCR-b (R819, R821, R852, R1028 and R1035); and LstaDCR (R837, R839, R870, R1049 and R1056). Conserved tyrosine (Y) residues of the 3′-end pockets (Lee et al., Reference Lee, Lee, Kim, Kim and Roh2023) within the captured peptides are highlighted in dashed lines.
Figure 5 Long description
The content consists of two parts labeled A and B. In part A, the comparative domain structures of Dicer proteins from various species are illustrated. Each species is represented with a cartoon model of the protein structure on the left and right sides, showing the base, head and core domains. The species include HsapDCR-1, DmelDCR-1, DmelDCR-2, CeleDCR-1, LstaDCR-1, BglaDCR, BpfeDCR, BtruDCR-a and BtruDCR-b. The domains are color-coded and labeled as helicase, Dicer partner-binding, helicase conserved C-terminal, Dicer dimerisation, Piwi-Argonaute-Zwille, Ribonuclease III and dsRNA-binding. Amino acid sequence ranges are provided for each domain. Part B shows a detailed alignment of amino acid sequences for the Dicer proteins, highlighting conserved motifs and residues. Specific amino acids are marked with arrows and numbers, indicating positions such as R822, R824, R855, R1030 and others. Below the sequence alignment, motifs are illustrated with sequences showing conserved tyrosine residues and other amino acids, emphasizing the recognition pockets in the Platform to PAZ region. The motifs are represented with letters and numbers, such as R763, R765, R796, R971 and R978, along with graphical representations of the amino acid sequences. The diagram provides a detailed comparison of the structural and sequence features of Dicer proteins across different species.
Argonaute (AGO) proteins in intermediate host snails. (A) Multiple alignments of the Piwi domain sequences highlighting conserved catalytic residues. (B) Maximum likelihood phylogeny reconstructed from the PAZ/Mid/Piwi regions. Sequences represent freshwater snails Biomphalaria glabrata (BglaAGO), Biomphalaria pfeifferi (BpfeAGO), Bulinus truncatus (BtruAGO) and Lymnaea stagnalis (LstaAGO); D. melanogaster (DmelAGO), C. elegans (CeleAGO), H. sapiens (HsapAGO) and 2 selected outgroups (Amphimedon queenslandica and Nematostella vectensis). Outgroups were used to root the tree, enabling polarization of bilaterian AGO diversification relative to early-branching metazoans. Alignment positions were filtered using partial deletion at 80% site coverage, and phylogenetic inference was performed under the LG + G model. Branch labels represent bootstrap support values ranging between 0 and 1. Values closer to 1 denote higher confidence in the branching relationship. This analysis involved 17 amino acid sequences, and there were a total of 582 positions in the final dataset. Scale bar indicates the expected number of amino acid substitutions per site.
Figure 6 Long description
Two panels labeled A and B. Panel A displays a multiple sequence alignment of AGO proteins, highlighting conserved catalytic residues. The sequences are aligned vertically, with positions marked at 715, 754, 797 and 920. Specific residues are highlighted in different colors: blue, red, green and yellow, indicating conservation across species. The sequences include CeleRDE-1, CeleALG-1, CeleALG-2, CeleALG-4, DmelAGO-1, DmelAGO-2, HsapAGO-1, HsapAGO-2, BpfeAGO, BtruAGO, BglaAGO-a, BglaAGO-b, LstaAGO, AguaAGO and NvecAGO. Panel B presents a phylogenetic tree with three clades labeled Clade I, Clade II and Clade III. Clade I is colored red, Clade II is blue and Clade III is cyan. Bootstrap support values are shown at nodes: 0.91 for Clade I, 0.87 for Clade II and 0.99 for Clade III. The tree is rooted with outgroups, indicated by an arrow pointing left. A scale bar at the bottom represents the expected number of amino-acid substitutions per site, marked as 2.0.
The functional divergence of DCR enzymes is well established in D. melanogaster, where DmDCR-1 and DmDCR-2 specialize in miRNA and siRNA biogenesis, respectively, through distinct substrate preferences and ATP requirements. Guided by these known structural features, we examined the architecture of snail DCR orthologues. Our analysis revealed that the overall domain organization of the snail sequences closely resembles that of HsDCR (mi/siRNA), CeDCR (mi/siRNA) and DmDCR-2 (siRNA only). Notably, in B. glabrata and B. pfeifferi, the DExH (HEL1) lobe was absent from the helicase N-terminus, a configuration reminiscent of ATP-independent DmDCR-1 (miRNA only) (Figure 5A). Further inspection of the Platform–PAZ region, however, demonstrated perfect conservation of both pre-siRNA 3′-end and pre-miRNA 5′-end recognition motifs across all snail orthologues (Figure 5B).
At the amino sequence level, all snail AGO sequences shared higher similarities with D. melanogaster AGO-1/2 isoforms. Among the 4 human AGO isoforms (HsAGO 1-4), the snail orthologues also showed higher reciprocal similarities with HsAGO-2. When the Piwi domains were further examined, all snail AGO sequences encoded the conserved DEDH catalytic motif (Figure 6A), supporting their potential to mediate siRNA-guided mRNA cleavage activity. ML phylogeny constructed using human, D. melanogaster, C. elegans and snail sequences showed 3 distinct clades, grouping snail AGOs, D. melanogaster AGO-1 (DmAGO-1) and C. elegans ALG-1 and ALG-2 within the same clade (Clade II) (Figure 6B).
Putative RNAi inhibitor proteins, including Staufen, are well conserved in the snail vectors
Six orthologues of known RNAi inhibitors were encoded in Biomphalaria and Bulinus genomes. In L. stagnalis, 9 putative RNAi inhibitors were respectively present (Table 2). Out of these, 6 putative RNAi inhibitors, adenosine deaminases ADAR-1 and ADAR-2, enhanced RNAi protein ERI-1, exoribonucleases XRN-1 and XRN-2 and Staufen STAU-1/2, were found across the snail datasets (Figure 3). Further domain characterization of the STAU proteins revealed that an N-terminal signal peptide is preserved only in coleopteran STAU orthologues (StauC) (Leptinotarsa decemlineata M1-G21, Pr = 0.53; Tribolium castaneum M1-G16, Pr = 0.79; Diabrotica virgifera virgifera M1-S19, Pr = 0.98). Conversely, an RBD4 domain found within other STAU sequences was absent (Figure 7A). Phylogenetic ML analysis similarly showed distinct branching of StauC sequences (bootstrap value = 1) from other STAU sequences, which were also clustered by species (Figure 7B).
Characterization of Staufen (Stau) proteins in intermediate host snails of helminths. (A) Domain profile of Staufen proteins among human (Hsap), coleopteran insects (Dvir, Ldec, Tcas), Caenorhabditis elegans (Cele), Drosophila melanogaster (Dmel) and snail vectors (Bgla, Bpfe, Btru, Lsta). Domain annotations and coordinates were predicted based on multiple sequence alignments and InterProScan searches and were plotted to scale. (B) A maximum likelihood cladogram was constructed based on full-length sequences of Stau proteins from all species. Tree was inferred using LG + G model, and all positions with less than 80% site coverage were eliminated. Bootstrap values lower than 50% were hidden from the tree. This analysis involved 14 amino acid sequences, and there were a total of 437 positions in the final dataset. Bgla, Biomphalaria glabrata; Bpfe, Biomphalaria pfeifferi; Btru, Bulinus truncatus; Dvir, Diabrotica virgifera virgifera; Hsap, Homo sapiens; Ldec, Leptinotarsa decemlineata; Lsta, Lymnaea stagnalis; Tcas, Tribolium castaneum. Hsap_Stau-1 (577aa isoform) and Hsap_Stau-2 (570aa isoform) were used.
Figure 7 Long description
The diagram consists of two parts labeled A and B. Part A displays domain profiles of Staufen proteins across various species, including CelecStau-1, DmelStau-1, HsapStau-1, HsapStau-2, LdecStauC, DvirStauC, TcasStauC, BglaStau-a, BglaStau-b, BpfeStau-a, BpfeStau-b, BtruStau and Lsta_Stau. Each protein is represented by a horizontal line with colored boxes indicating different domains such as RBD1, RBD2, RBD3, RBD4 and others. The scale at the bottom ranges from 0 to 1000. Part B is a circular phylogenetic tree showing evolutionary relationships among Staufen proteins from different species, with branches labeled by species names like Celec, Dmel, Hsap, Ldec, Dvir, Tcas, Bgla, Bpfe, Btru and Lsta. The tree scale is marked as 0.05, indicating branch length proportional to evolutionary distance.
PPI network and GO enrichments supported functional snail orthologues
A PPI network constructed using B. glabrata protein sequences resulted in 15 clusters (Figure 8A) and showed an enrichment P-value < 1.0e − 16. Clathrin-1 (CHC-1) showed biological interactions with AP-2, DYN-1 and CLINT-1. Interactions between DCR-1 and TRBP-2, DRSH-1 and DGCR8, RAN and exportins (XPOs), as well as AGO-2 and deadenylation complexes and decapping enzymes (e.g. TNRC6 and DDX6, respectively) were also represented (Figure 8A). GO enrichments also showed strong involvement of the snail protein orthologues in the BP of gene silencing by RNA (signal = 4.29, FDR = 9.83e − 17), with marked MFs in dsRNA binding (signal = 3.29, FDR = 3.55e − 12) and ribonuclease activity (signal = 1.82, FDR = 1.31e − 07). In the CCs, putative snail proteins associated with the RISC complex (signal = 3.14, FDR = 9.90e − 11) and the processing bodies (P-bodies; signal = 2.53, FDR = 8.51e − 10) showed the strongest signals (Figure 8B, Supplementary Material 2). Six protein orthologues were predicted to constitute the RISC complex (DCR-1, AGO-2b, TSN-1, TRBP-2, XPO-5 and DCP-2), while 8 P-body protein orthologues (DDX6, AGO-2b, TNRC6, PAN-2, PAN-3, CCR4-NOT, XRN-1 and DCP-2) were predicted (Figure 8B, Supplementary Material 2). Figure 9 summarizes the hypothetical RNAi mechanism in freshwater snails based on the findings from this study.
Protein–protein interaction (PPI) and Gene Ontology (GO) analyses. (A) PPI network of putative RNAi complements in Biomphalaria glabrata. Only connected nodes are shown. Straight unbroken edges between nodes indicate direct functional interactions, while dashed edges indicate indirect functional associations. (B) Gene Ontology (GO) enrichment plot showing predicted biological process (BP), cellular component (CC) and molecular function (MF) terms associated with the identified B. glabrata protein sequences. False discovery rate (FDR) values show the significance level of enrichment. Signal is the weighted harmonic mean between the observed-versus-expected ratio and −log (FDR), which ensures more intuitive ordering of enriched terms. Gene count is the number of protein-coding genes in each enrichment term. For each GO category, the top 10 terms ranked by their signals were shown. The list of genes under each GO term is available in Supplementary Material 2.
Figure 8 Long description
A showing a network diagram made of labeled circular nodes connected by lines. Node labels include Chc1, Ap2, Clint1, Dyn1, Dcr1, Trbp2, Drsh1, Dgcr8, Ran, Xpo1, Xpo2, Xpo5, Ago2, Tsn1, Ddx6, Tnrc6, Pan2, Pan3, Ccr4, Not, Xrn1 and Dcp2. Some connections are drawn with solid lines and some with dashed lines. B showing a dot plot with a categorical y-axis listing gene ontology terms and an x-axis labeled, Signal. No units are shown for Signal. The x-axis tick labels are 0, 1, 2, 3 and 4. The y-axis terms, from top to bottom, are RNA binding; Ribonuclease activity; pre microRNA binding; nuclease activity; Exoribonuclease activity producing 5 phosphate ends; Exonuclease activity; Double stranded RNA binding; Catalytic activity acting on RNA; ATPase coupled transmembrane transporter activity; 5 3 exoribonuclease activity; RNA catabolic process; pre microRNA processing; Post transcriptional gene silencing; Nuclease containing exonuclease activity; Negative regulation of translation; Negative regulation of gene expression; RNA catabolic process; microRNA processing; RNA mediated gene silencing; double stranded RNA binding; RNA induced silencing complex; RISC loading complex; RISC complex; Ribonucleoprotein granule; P body; Endoribonuclease; Deadenylase; Cytoplasmic ribonucleoprotein granule; Clathrin coated vesicle. A legend labeled, minus log subscript 10 open parenthesis F D R close parenthesis, shows tick labels 0, 5, 10 and 15. A legend labeled, G O Category, shows three marker shapes. A legend labeled, Gene Count, shows marker sizes labeled 10, 20 and 30. Each listed term has one plotted marker. Markers appear at Signal values from near 0 up to about 4. Markers at the highest Signal values are near 4 and markers at the lowest Signal values are near 0. Marker sizes vary according to the Gene Count legend and marker shapes vary according to the G O Category legend.
Hypothetical model for mechanistic RNAi pathways in a single snail cell. In the miRNA pathway (traced with dotted arrows), primary miRNA (pri-miRNA) is first transcribed from genomic DNA by RNA polymerase II. Drosha associates with DGCR8 (referred to as Pasha in flies and nematodes) and cleaves pri-miRNA to generate precursor miRNA (pre-miRNA), which is exported from the nucleus, primarily by Exportin (XPO)-5, into the cytoplasm where it is diced into smaller ∼22-nt miRNA duplexes by the Dicer protein. The miRNA duplex is incorporated into the RNA-induced silencing complex (RISC) (now known as miRISC), where it is unwound into 2 separate strands by Argonaute 2 (AGO-2). The guide (antisense) strand, through AGO-2, locates and pairs insufficiently with a cognate cytosolic mRNA sequence, while AGO-2 mediates translational repression, and hence decay, of the cognate mRNA by recruiting TNRC6, a GW182 protein. The siRNA pathway (traced with bold arrows) is triggered by the introduction of exogenous long dsRNA into the cell. The long dsRNA substrate is diced into ∼22-nt siRNA duplexes. Short ∼22-nt siRNAs introduced directly into the cell do not undergo Dicer cleavage. The siRNA duplex is loaded into the RISC, now called siRISC, and follows the same unwinding and pairing process. The guide siRNA strand, however, forms a full complementary pairing with the cognate mRNA sequence, culminating in degradative cleavage of the mRNA sequence. P-bodies are cytoplasmic ribonucleoprotein granules that act as sites for mRNA translational repression or degradation. Negative RNAi regulators (RNAi inhibitors) may interfere with the overall process at specific stages.
Figure 9 Long description
The diagram illustrates RNA interference pathways within a snail cell, highlighting miRNA and siRNA processes. In the nucleus, primary miRNA is transcribed from genomic DNA by RNA polymerase II. Drosha and DGCR8 cleave pri-miRNA to form precursor miRNA, which is exported to the cytosol by Exportin-5. In the cytosol, Dicer processes pre-miRNA into miRNA duplexes, which are incorporated into the miRISC complex. Argonaute 2 unwinds the miRNA duplex, allowing the guide strand to pair with mRNA, leading to translational repression and mRNA decay via TNRC6. The siRNA pathway begins with exogenous long dsRNA, diced into siRNA duplexes by Dicer and loaded into the siRISC complex. The guide siRNA strand pairs fully with mRNA, resulting in mRNA cleavage. P-bodies are sites for mRNA repression or degradation. The diagram includes various proteins and complexes such as ADAR dimer, Stau and Xrn1 and indicates RNAi spread and import processes. Inhibitors and various protein interactions are depicted, including receptors, dynamin, clathrin and ABC transporters.
Expression levels of Ago2, Dcr1 and Eri1 genes were independent of snail size
Agarose gel analysis showed the expected amplicon size for each specifically amplified primer set (Figure 2B). Based on snail size, the expression of Dcr1, Ago2 and Eri1 suggested differences in the average relative mRNA expression levels, but these differences were not significant (P > 0.05) between groups (Figure 2C–E).
Enriched expression of Chc1, Sid1-like and Stau1/2 genes were tissue-specific
Tissue-specific expressions of Chc1, Sid1L and Stau were assessed to determine their localization and extrapolate their roles in the snail RNAi. Across the isolated tissues, expression of the endocytic pathway Chc1 gene was found to be significantly different (Kruskal–Wallis X 2 = 13.863, df = 4, P = 0.008). The expression was notably higher in the HP than in other tissues, especially the albumen gland (P = 0.011) and the trunk (P = 0.03) (Figure 2F). Sid1L expression also varied significantly among the isolated tissues (Kruskal–Wallis X 2 = 17.361, df = 4, P = 0.002), with the highest expression in the trunk, especially in comparison with the HP (P = 0.0009) and the HF region (P = 0.039) (Figure 2G). Out of the 4 B. glabrata STAU orthologues, 1 transcript (XM_056016125.1) examined for expression analysis showed higher expression in the OV than other tissues (Kruskal–Wallis X 2 = 11.784, df = 4, P = 0.019), most significantly relative to the trunk region (P = 0.013) (Figure 2H). Except in the trunk (P = 0.841), the expression levels of Chc1 were comparatively higher than Sid1L in all tissues examined (P = 0.008 in AG and HP; P = 0.038 in HF and OV) (Figure 2I).
Discussion
In reverse-genetics studies, CRISPR (Clustered Regularly Interspaced Short Palindromic Repeats)-aided knockout approaches are gaining immense traction but remain technically challenging in many malacology laboratories, although recent major breakthroughs have been made (Abe and Kuroda, Reference Abe and Kuroda2019; Batzel et al., Reference Batzel, Wang, Bock, Chen, Neal, Lopez-anido, Lee, Tjeerdema, Ignatoff, Patil, Ramirez, Lesoway, Hamdoun and Lyons2025; Oonuma et al., Reference Oonuma, Kuroda, Uchida and Zhang2025). Until these approaches become widely established, RNAi remains the main reverse-genetics tool in molecular snail research and deserves full consideration. In this study, we delineate the diversity of RNAi effectors in 4 freshwater snails to provide insight into the scope and constraints of their functional RNAi.
The reduced BUSCO completeness for the L. stagnalis genome suggests that the genome does not contain as many pan-eukaryotic genes as other snail genomes examined in this study. In the original study, where the genome assembly was assessed under a mollusc-specific lineage (mollusca_odb10), BUSCO completeness scored 91.8% (Koene et al., Reference Koene, Jackson, Nakadera, Cerveau, Madoui, Noel, Jamilloux, Poulain, Labadie, Da Silva, Davison, Feng, Adema, Klopp, Aury, Wincker and Coutellec2024). This indicates that the BUSCO outcome was strongly influenced by gene prediction settings, including the lineage choice. The high level of gene duplication observed with Bu. truncatus genome could be linked to the tetraploid nature of the snail (Young et al., Reference Young, Stroehlein, Wang, Korhonen, Mentink-kane, Stothard, Rollinson and Gasser2022). Although RNAi effectors are considered pan-eukaryotic, differences in the genome completeness of the snails did not impact the recovery of RNAi effector proteins. Interestingly, the number of RNAi effector proteins identified in L. stagnalis was slightly higher than in the Biomphalaria and Bulinus species examined. Two factors may have contributed to this pattern. First, searches for RNAi effector proteins were performed using the L. stagnalis transcriptome, which may capture genes that are fragmented or ‘missing’ but not biologically absent in the genome assemblies. Second, phylogenetic divergence between lymnaeid (Lymnaea) and planorbid (Biomphalaria and Bulinus) snails (Saadi et al., Reference Saadi, Davison and Wade2020) might have caused family-specific loss or expansion of RNAi gene sets. Altogether, the results explain the highly conservative nature of the core RNAi effectors, even when genome completeness is less robust.
While we found that the snails lack classical orthologues of the dsRNA-importing SID-1 protein, they may instead employ the CME pathway for dsRNA internalization, as has been found in other organisms including insects (Xiao et al., Reference Xiao, Gao, Xu, Liang, Li, Yao and Zhu2015), mites (Zhou et al., Reference Zhou, Wan, Jian, Guo, Zhang, Shi, Yang, Li, Liu and Ding2023) and fungi (Wytinck et al., Reference Wytinck, Sullivan, Biggar, Crisostomo, Pelka, Belmonte and Whyard2020). However, this alternative pathway may be slower and may require higher concentrations of dsRNA to initiate a systemic RNAi response (Feinberg and Hunter, Reference Feinberg and Hunter2003; Kobayashi et al., Reference Kobayashi, Tsukioka, Komoto, Uchino, Sezutsu, Tamura, Kusakabe and Tomita2012). Importantly, CME is only one of the several core endocytic trafficking pathways in eukaryotes. It is highly selective in its endocytosed cargo, repelling many other substances that are imported via the non-CME pathways (Kumari et al., Reference Kumari, Mg and Mayor2010).
Snail SID1-like orthologues showed closer evolutionary proximities to CeCHUP-1 (also known as TAG-130) and HsSIDT-1, both of which are more exclusively involved in cholesterol uptake, rather than dsRNA uptake (Valdes et al., Reference Valdes, Athie, Salinas, Navarro and Vaca2012; Mendez-Acevedo et al., Reference Mendez-Acevedo, Valdes, Asanov and Vaca2017). In C. elegans, SID-1 is ubiquitously expressed in non-neuronal tissues, especially in the cell periphery of externally exposed tissues where it readily imports exogenous dsRNAs (Winston et al., Reference Winston, Molodowitch and Hunter2002). Conversely, the paralogous CeCHUP-1 is densely localized in the nematode’s intestine where most dietary cholesterol uptake occurs (Valdes et al., Reference Valdes, Athie, Salinas, Navarro and Vaca2012). In this study, the Sid1L gene also showed localized and high expression in the trunk section, which included the snail’s intestine and other organs listed previously. In B. glabrata, however, the digestive gland–gonad (DGG) complex, encompassing HP and OV, forms the main site for cholesterol uptake (Shetty et al., Reference Shetty, Fried and Sherma1990). In these organs, we observed no increased expression of the Sid1L gene, suggesting that while the snail SID1L/CHUP1L orthologues are structurally competent for cholesterol uptake, their main activities are limited to certain organs and tissues outside the DGG complex. A major limitation here is that we were unable to achieve organ-specific resolution of the trunk area where Sid1L was significantly overexpressed. A similar study in Litopenaeus vannamei shrimp found high expression of the Sid1L gene in the gills, while CME mainly facilitated uptake of ingested dsRNA into the hepatopancreatic cells (Maruekawong et al., Reference Maruekawong, Namlamoon and Attasart2022). This provides the most proximate insight, as gills in shrimp are analogous to the pseudobranch in Biomphalaria, which partly constitutes the snail trunk tissue. Moreover, we similarly found in the snail HP (and OV) a significant enrichment of Chc1, the main effector of the CME pathway. Previously, fluorescence microscopy showed HP and OV as the major sites for RNAi uptake in B. glabrata (Knight et al., Reference Knight, Miller, Liu, Scaria, Woodle and Ittiprasert2011). Findings from the present study, reinforced by the observations from previous works, further support the clathrin-mediated cellular uptake of silencing RNAs in B. glabrata, mainly through the HP and OV.
Interestingly, more recent works showed that the cholesterol-trafficking HsSIDT-1 protein facilitates dsRNA uptake via its phospholipase D activity (Sun et al., Reference Sun, Xu, Yang, Hou, Luo, Zhang, Shan, Huang, Yao, Chen, Li and Zhou2024) at low pH (Zheng et al., Reference Zheng, Yang, Guo, Qi, Lu, Xiao, Gao, Liu, Yang, Zhou, Nguyen, Zhu, Sun, Zhang and Ji2024), although knockdown of clathrin in the same human embryonic kidney cell line (HEK293T) most markedly attenuated dsRNA uptake (Sun et al., Reference Sun, Xu, Yang, Hou, Luo, Zhang, Shan, Huang, Yao, Chen, Li and Zhou2024). Mammalian studies further revealed that CHUP1-like proteins may specialize in importing only cholesterol-conjugated dsRNAs (Brück et al., Reference Brück, Pascolo, Fuchs, Kellerer, Glocova, Geisel, Dengler, Yazdi, Röcken and Ghoreschi2015; Mendez-Acevedo et al., Reference Mendez-Acevedo, Valdes, Asanov and Vaca2017), a method that has also significantly enhanced RNAi delivery in B. glabrata (Portet et al., Reference Portet, Galinier, Lassalle, Faille, Gourbal and Duval2021). Together, this suggests that intermediate snail hosts import external silencing RNAs potentially via both CME- and CHUP1L/SID1L-dependent pathways, as similarly noted in Leptinotarsa decemlineata beetle (Cappelle et al., Reference Cappelle, De Oliveira, Van Eynde, Christiaens and Smagghe2016) and L. vannamei shrimp (Maruekawong et al., Reference Maruekawong, Namlamoon and Attasart2022), in a likely tissue-specific manner.
Dicer (DCR) as well as AGO proteins vary in their domain, microdomain and amino acid structures. These architectural variations commonly underscore evolutionary specialization of their isoforms towards distinct small RNA pathways (Park et al., Reference Park, Heo, Tian, Simanshu, Chang, Jee, Patel and Kim2011; Nakanishi et al., Reference Nakanishi, Weinberg, Bartel and Patel2012). While the Platform–PAZ region of DCR enzymes consistently conserves its 3′-end binding pockets required for pre-siRNA processing, the 5′-end binding pockets that are critical for pre-miRNA processing are not conserved in siRNA-restricted DCRs (e.g. in DmDCR-2) (Park et al., Reference Park, Heo, Tian, Simanshu, Chang, Jee, Patel and Kim2011; Lee et al., Reference Lee, Lee, Kim, Kim and Roh2023). Results from this study, however, indicate that each snail species examined retained the structural requirements in their DCR to mediate both miRNA and siRNA biogenesis. This is potentially in an ATP-independent manner for Biomphalaria species, as DExH (HEL1) moiety of the helicase domain that facilitates ATP hydrolysis was not found. Such DCRs may be inefficient in cleaving blunt-ended siRNA precursors (dsRNAs without 2-nt 3′ terminal overhangs) (Donelick et al., Reference Donelick, Talide, Bellet, Aruscavage, Lauret, Aguiar, Marques, Meignin and Bass2020).
Prior to this study, phylogenetic analysis suggested an evolutionary loss of siRNA-processing AGO proteins in molluscs (Formaggioni et al., Reference Formaggioni, Cavalli, Hamada, Sakamoto, Plazzi and Passamonti2024). Our phylogenetic inference similarly suggested that snail AGOs evolved within the miRNA-specific AGO clade. Although phylogenetic analyses did not recover clear cleavage-competent siRNA AGOs in the snails, multiple sequence alignments revealing full conservation of the Piwi catalytic residues, coupled with many successful reports of dsRNA/siRNA-based silencing in snails (Portet et al., Reference Portet, Galinier, Lassalle, Faille, Gourbal and Duval2021; Xu et al., Reference Xu, Y-w-q, Habib, Li, Yuan, Ke, Jiang, Dong and Zhao2023), show otherwise. These findings suggest that AGO-2-mediated mRNA cleavage in the snails may be evolutionarily divergent, yet mechanistically canonical.
Across the groups of organisms examined, STAU protein exhibits functional plasticity in RNA metabolism. STAU-1 and STAU-2 are negative regulators of RNAi (Ren et al., Reference Ren, Veksler-lublinsky, Morrissey and Ambros2016; Ehses et al., Reference Ehses, Schlegel, Schroger, Schieweck, Derdak, Bilban, Bauer, Harner and Kiebler2022), except for coleoptera-specific StauC proteins that promote dsRNA-to-siRNA processing (Yoon et al., Reference Yoon, Mogilicherla, Gurusamy, Chen, Chereddy and Palli2018). Unlike the putatively cytosolic snail STAU orthologues, N-terminal signal peptides found in StauC sequences indicate their typical localization in the endoplasmic reticulum, which suggests an ability to process dsRNA prior to secretion for inter- or intra-cellular transport (Koo and Palli, Reference Koo and Palli2024). Conversely, the fourth dsRBD4 found only in non-coleopteran STAU sequences may be involved in the RNAi inhibition process. Many metazoans, including C. elegans (Legendre et al., Reference Legendre, Campbell, Kroll-conner, Anderson, Kimble and Wickens2013), D. melanogaster (Micklem et al., Reference Micklem, Adams, Grunert and St Johnston2000), Danio rerio (Ramasamy et al., Reference Ramasamy, Wang, Quach and Sampath2006), Xenopus laevis (Allison et al., Reference Allison, Czaplinski, Git, Adegbenro, Stennard, Houliston and Standart2004; Bilogan and Horb, Reference Bilogan and Horb2012) and Mus musculus (Cao et al., Reference Cao, Du, Chen, Wang, Zhang, Liu, Liu, Weng, Liang and Ma2016), exhibit high and stable expression of STAU-1 and STAU-2 in the oocytes, ovaries and embryos. We similarly found the highest expression of Stau in the snail OV. Considering both the function and the pattern of expression, it may be suggested that STAU-1/2 shields key developmental transcripts from the cell’s own RNAi machinery.
The adenosine deaminases acting on RNAs (ADARs), found in all snail species examined, are also RNAi inhibitors that may play contrasting roles in RNAi processing. These proteins, typically containing one or multiple dsRBDs, an A-deaminase domain, and one or more Z-binding domains, can edit dsRNAs by converting adenosine (A) residues into inosines (I). This A-to-I editing antagonizes RNAi process by altering sequence regions that are required for small RNA precursor binding and/or cleavage by Drosha and Dicer proteins (Nishikura, Reference Nishikura2010). While the RNA editor role of ADAR is facilitated by its ADAR/ADAR homodimerization, ADAR/DCR heterodimerization may rather define its role as an RNAi promoter (Ota et al., Reference Ota, Sakurai, Gupta, Valente, Wulff, Ariyoshi, Iizasa, Davuluri and Nishikura2013).
The generation of secondary endogenous siRNAs (endo-siRNAs), following the introduction of initial RNAi triggers, drives a more robust RNAi response in some organisms. RdRPs, majorly RRF-1 and EGO-1, and their upstream proteins, RDE-10 and RDE-11, facilitate the production of endo-siRNAs (Zhang and Ruvkun, Reference Zhang and Ruvkun2012), none of which were present in the snails. The phylogenetic distribution of RdRPs among metazoans is sporadic. They appear to have been lost in numerous lineages, including gastropods (Pinzon et al., Reference Pinzon, Bertrand, Subirana, Busseau, Escriva and Seitz2019). However, an RdRP-independent RNAi amplification mechanism mediated by endogenous reverse transcriptases has been suggested in Drosophila (Tassetto et al., Reference Tassetto, Kunitomi and Andino2017) and may represent an alternative mechanism for endo-siRNA production in many other organisms.
Core components responsible for transgenerational RNAi in C. elegans are NRDE-1, NRDE-2, NRDE-3 and NRDE-4. NRDE-3, a secondary AGO, binds and shuttles secondary siRNAs from the cytoplasm into the nucleus, where it associates with other NRDE proteins (NRDE-2, -4 and -1, sequentially) to mediate transcriptional silencing of nascent mRNA through histone methylation. This way, the silencing effect is heritable by the progeny (Burton et al., Reference Burton, Burkhart and Kennedy2011). Our in silico analysis recovered no orthologues of the core nuclear RNAi effectors in the snails examined, and no evidence for heritable RNAi has been published from the phylum Mollusca. Substantiating all findings from the in silico data, the significantly low P-value of the PPI network indicates that the putative protein orthologues found in the snails could interact to form functional RNAi units. The PPI networks were consistent with the common functional grouping of the RNAi effectors, as well as the general mechanistic RNAi pathways.
Argonaute-2 (Ago-2) and Dcr-1 are 2 major dsRNA-processing genes, while Eri1 is a major RNAi-inhibiting gene. Our qPCR analyses of these genes between snail developmental stages displayed stochastic variation among size-matched replicates. This indicated that RNAi efficiency is not influenced primarily by snail size. Shell size itself has imperfect correlation with snail’s age. Age-related differences, along with microenvironmental heterogeneity such as individual responses to stress, could contribute to temporal fluctuations in mRNA expression, even under uniform rearing conditions. Importantly, these observations remain consistent with existing experimental evidence. In B. glabrata, for example, RNAi-mediated knockdown has been successfully achieved across snails of different sizes, ranging from juveniles to adults (Portet et al., Reference Portet, Galinier, Lassalle, Faille, Gourbal and Duval2021; Xu et al., Reference Xu, Y-w-q, Habib, Li, Yuan, Ke, Jiang, Dong and Zhao2023).
Overall, our study has contributed new mechanistic insights into the RNAi biology of freshwater snails. Although our literature search may not have captured all known RNAi effectors, and new RNAi effectors may have been published after our search was conducted, this study presents the most comprehensive diversity of RNAi effectors in freshwater snails. We showed that the diversity of RNAi effectors between our target snail species is similar, suggesting similar RNAi mechanisms and degree of susceptibility to exogenous gene silencing across the species. Despite lineage-specific absence of some effector components, our findings underscore the wide potential for RNAi use in trematode intermediate host snails as an experimental tool. These approaches can be used for functional genomics or as a future RNAi-based disease control strategy. Further work, including spatial and functional validation of these putative RNAi effectors, will strengthen the mechanistic framework established here. The development of RNAi approaches that can be manipulated more reliably and efficiently in parasite-transmitting snails will substantially enhance our understanding of these important organisms.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0031182026102819.
Data availability statement
All codes and data files used in data analysis are available on GitHub at https://github.com/DamiF25/Snail_RNAi_effectors.
Acknowledgements
We are grateful to Dr Gabriel Rinaldi for kindly providing the initial snail populations.
Author contributions
P.M., D.F. and C.L. conceived and designed the study. D.F., C.L. and P.M. performed bioinformatic analyses. D.F. performed gene expression analysis. D.F. collated all data for analysis and visualization. D.F. wrote the initial draft. C.L., D.F., D.W., G.G. and P.M. revised the initial draft and approved its final form.
Financial support
D.F. and C.L. received Postgraduate Research Studentships from the Northern Ireland Department for the Economy (DfE, NI) and the Northern Ireland Department of Agriculture, Environment and Rural Affairs (DAERA, NI), respectively. P.M. received an Agility+ award from Queen’s University Belfast. Funders had no role in the study design, analysis and interpretation of the results.
Competing interests
The authors declare that there are no conflicts of interest.
Ethical standards
Not applicable.
Open access
