Transcriptome profiling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research
BMC Genomics volume 22, Article number: 144 (2021)
Host immune function can contribute to numerous ecological/evolutionary processes. Ecoimmunological studies, however, typically use one/few phenotypic immune assays and thus do not consider the complexity of the immune system. Therefore, “omics” resources that allow quantifying immune activity across multiple pathways are needed for ecoimmunological models. We applied short-read based RNAseq (Illumina NextSeq 500, PE-81) to characterise transcriptome profiles of Lymnaea stagnalis (Gastropoda), a multipurpose model snail species. We used a genetically diverse snail stock and exposed individuals to immune elicitors (injury, bacterial/trematode pathogens) and changes in environmental conditions that can alter immune activity (temperature, food availability).
Immune defence factors identified in the de novo assembly covered elements broadly described in other gastropods. For instance, pathogen-recognition receptors (PRR) and lectins activate Toll-like receptor (TLR) pathway and cytokines that regulate cellular and humoral defences. Surprisingly, only modest diversity of antimicrobial peptides and fibrinogen related proteins were detected when compared with other taxa. Additionally, multiple defence factors that may contribute to the phenotypic immune assays used to quantify antibacterial activity and phenoloxidase (PO)/melanisation-type reaction in this species were found. Experimental treatments revealed factors from non-self recognition (lectins) and signalling (TLR pathway, cytokines) to effectors (e.g., antibacterial proteins, PO enzymes) whose transcription depended on immune stimuli and environmental conditions, as well as components of snail physiology/metabolism that may drive these effects. Interestingly, the transcription of many factors (e.g., PRR, lectins, cytokines, PO enzymes, antibacterial proteins) showed high among-individual variation.
Our results indicate several uniform aspects of gastropod immunity, but also apparent differences between L. stagnalis and some previously examined taxa. Interestingly, in addition to immune defence factors that responded to immune elicitors and changes in environmental conditions, many factors showed high among-individual variation across experimental snails. We propose that such factors are highly important to be included in future ecoimmunological studies because they may be the key determinants of differences in parasite resistance among individuals both within and between natural snail populations.
Several fields of ecology and evolution increasingly recognise host immune function as an essential contributor to biological processes (see ). For instance, the immune system plays critical roles in life-history evolution [2, 3], sexual selection [4, 5], and responses/adaptations of organisms to environmental change [6,7,8]. This recognition has given rise to an interdisciplinary field of ecological immunology (or ecoimmunology; see ). Ecoimmunological studies, especially in invertebrates, typically measure the end products of one or few immunological cascades that are controlled by several genes (e.g., [10, 11]). Thus, the field relies on quantitative genetic theory, initially motivated by the assumed simplicity of innate-type invertebrate immune systems with non-specific recognition and killing mechanisms. Over the recent decades, however, comparative immunology with aid from genomics, has shown invertebrate immune systems to be complex and diversified across the phyla (e.g., [12,13,14,15,16]). In fruit flies, for instance, specific immune pathways respond towards Gram-positive and Gram-negative bacteria, as well as fungi [17,18,19].
The realised complexities of invertebrate immune systems provide new challenges and opportunities for studies focusing on ecological and evolutionary questions on immune function. This is because typical ecoimmunological studies that focus on one/few immunological mechanisms are incapable of describing the multivariate “immune phenotypes” (sensu ) of organisms. That, however, would be important because each immunological pathway can respond differently to various selective agents (e.g., parasite/pathogen species) and may be traded-off with different physiological, life-history, and other immune traits [10, 14, 21,22,23]. Additionally, the expression of immune traits, as well as the associated trade-offs, may depend on environmental conditions such as resource availability and temperature (e.g., [8, 24, 25]), which could be due to stress-related responses or altered metabolism. Therefore, a comprehensive characterisation of different immunological and physiological mechanisms, as feasible by genomics-type approaches, in ecologically relevant experiments utilising invertebrates is essential. This expansion has been done successfully in some insects such as bumblebee (e.g., [26, 27]) and red flour beetle [28, 29]. Considering the extensive diversity of invertebrate phyla, development and utilisation of such genomic data also in other taxonomic groups would further expand the potential of ecoimmunology.
Mollusca is the second-largest animal phylum after Arthropoda. In this phylum, Gastropoda represents the largest taxonomic class with 40,000 to 150,000 living species . Gastropods inhabit aquatic and terrestrial habitats, and incur disease from viruses and bacteria (e.g., [31,32,33]) but also from specialist flatworm parasites called digenetic trematodes [34, 35]. Trematodes have complex life cycles that typically involve a gastropod as an intermediate host. In snails, trematodes reproduce asexually to produce transmission stages that infect the next host in the life cycle (another intermediate host or a definitive host). Most of the attention on molecular immunology in gastropods has focused on Biomphalaria glabrata (Planorbidae, Hygrophila, Panpulmonata), a tropical freshwater snail that transmits the human blood fluke, Schistosoma mansoni [36,37,38]. Attention to other species has been comparatively limited although several gastropods, including pond snails of the family Lymnaeidae (Hygrophila, Panpulmonata), transmit medically and veterinary-relevant parasites (e.g., Fasciola hepatica, Diplostomum spp., Trichobilharzia spp., [39,40,41]) in temperate regions. In this family, diversity and function of circulating defence cells, haemocytes, have been investigated in the great pond snail, Lymnaea stagnalis [42,43,44,45,46,47,48,49]. Additionally, a draft genome of L. stagnalis is available .
Here, we applied short-read based RNAseq to characterise transcriptome profiles of L. stagnalis exposed to various immune elicitors (injury, bacterial and trematode pathogens) and changes in environmental conditions (temperature, food availability) using a genetically diverse laboratory population of snails. Lymnaea stagnalis is a model organism in multiple biological disciplines (reviewed in ), including ecological immunology (e.g., [8, 52,53,54,55,56]). Earlier ecoimmunological work, however, employs only a few phenotypic immune assays. Therefore, our data analyses focused on examining the suitability of a broad range of candidate immune genes for quantifying snail immune activity. Additionally, we investigated components of snail physiology/metabolism that may be related to immune responses. In a two-step process, we first annotated the transcriptome using previously characterised immune/physiology-related proteins/genes from other organisms. After that, we evaluated variation in the transcription of those candidates to detect transcripts that responded to experimental treatments and/or that showed high among-individual variation. We propose that the use of candidate genes exhibiting high variation in transcription allows a comprehensive examination of variation in polymorphic snail immune responses in future ecoimmunological studies.
Sequencing and assembly
Sequencing (Illumina NextSeq 500 platform) of all 48 libraries (generated from RNA samples extracted from homogenised whole-body tissues) yielded a total of 1.08 billion paired-end 81 nt long reads (PE-81). The raw reads representing a total of 175.2 Gb sequence data were deposited in the NCBI Sequence Read Archive (accession number PRJNA664475). The de novo assembly (trinity v2.0.6) of the reads from nine libraries (one per each experimental treatment; treatments: untreated, anaesthesia, wounding, injection with Escherichia coli, injection with Micrococcus lysodeikticus, injection with healthy snail tissue, injection with trematode-infected snail tissue, elevated temperature, food deprivation) yielded 264,746 transcripts with N50 of 1589 nt, mean length of 551 nt, and median length of 188 nt (146.0 Mb in total). These transcripts included 68,473 open reading frames (ORFs) that were longer than 100 aa (TransDecoder v3.0.1). From this initial full assembly, the removal of likely contaminant sequences that lacked a significant similarity hit with L. stagnalis genome (GenBank accession number GCA_900036025.1) and that were not previously characterised from this species, yielded a final reference transcriptome (139.1 Mb) consisting of 226,116 contigs with a mean transcript length of 615 nt (assembly available in https://doi.org/10.5281/zenodo.4044169). BUSCO analysis indicated high completeness of the reference transcriptome based on the detection of 98.8% (complete plus partial) of the set of metazoan universal single-copy orthologs [complete: 96.3% (duplicated: 28.2%), fragmented: 2.5%, missing: 1.2%, N = 978].
Sequences supporting species identification
The reference transcriptome included sequences that represent parts of the snail rDNA cassette, including the complete ITS1 and ITS2 regions (GenBank accession numbers MT864603 and MT864602, respectively). Lymnaea stagnalis was confirmed as the species identity of the experimental snails by the highest sequence similarities (BLASTN) of these sequences with GenBank entries from this species. Similarly, snails exposed to tissue extracts from trematode-infected gonads provided LSU-derived sequences with the highest similarity with Clinostomidae and Plagiorchiidae families of trematodes (GenBank accession numbers MT872505 and MT872506, respectively), which confirms the morphological identification of cercaria released by the donor snails.
The mining of the L. stagnalis reference transcriptome (226,116 contigs) yielded numerous factors that contribute to immune responses from non-self recognition to the elimination of pathogens (Fig. 1, Additional file 1). For accommodation of PRRs, L. stagnalis expressed seven variants (i.e., transcripts with unique ORFs) of long-form (no short forms) peptidoglycan recognition proteins (PGRPs), four variants of Gram-negative bacteria binding-proteins (GNBPs) and a repertoire of lectins. Representatives of the detected lectin families included two FREPs of the VIgLs (variable immunoglobulin and lectin-domain-containing molecules), galectins comprising of either one, two or four carbohydrate recognition domains, multiple Chi-lectins, as well as L- and M-type lectins. The latter two families may, however, function mostly in housekeeping roles.
Along with one TLR sequence encoding the canonical architecture of tandemly arranged leucine-rich repeats (LRRs), a transmembrane region and a C-terminal Toll-interleukin receptor (TIR) domain, the associated signalling NFκB pathway for immune activation was represented by multiple components such as adaptor proteins for intracellular signalling downstream of the receptor (MyD88, TRAF), transcription factor (Rel protein) NF-κB, and various regulators (e.g., SARM, IκB). Additionally, three categories of cytokines for intercellular communication (including immune responses and inflammation) were detected; one variant of macrophage migration inhibitory factor (MIF), three variants of interleukin 17 (IL-17) and 23 variants of tumor necrosis factor (TNF).
The search of the complete reference transcriptome for antimicrobial defences yielded a single family of macin antimicrobial peptides (AMPs, 6 variants), as well as several families of antimicrobial proteins, with an abundant representation of lipopolysaccharide-binding /bactericidal permeability-increasing proteins (LBP/BPI, 3 variants), L-amino acid oxidases (LAAOs, 3 variants), lysozymes (4 variants), and transcripts encoding for cytolytic β pore-forming toxins. These latter sequences were designated Lymnaea-lysins (3 variants) and stagnalysins (3 variants), in analogy to the orthologous Biomphalysins and glabralysins described from B. glabrata [57, 58]. Central components from a gene network facilitating the production and management of reactive oxygen species for oxidative killing mechanisms were detected, including the membrane-bound enzyme complex NADPH-oxidase (NOX, with Cytochrome B 245 as the main component, 3 variants) that produces superoxide anions, as well as dual oxidases (DUOX, 2 variants) that catalyse the synthesis of anion superoxide and hydrogen peroxide. Moreover, dual oxidase maturation factor (DUOXA, 1 variant) that activates DUOX was found, as well as additional proteins involved in the biosynthesis of other oxidative compounds. Furthermore, eight variants of laccase and two of tyrosinase that can contribute to PO/melanisation-type defence were recorded. The search also yielded transcripts encoding for antioxidant enzymes and proteins regulating oxidative damage. These transcripts included two types of superoxidase dismutase enzymes against ROS [SOD (2 variants) and MnSOD (1 variant)], catalase (CAT, 2 variants), glutathione peroxidase (3 variant), and glutathione reductase (1 variant).
Our analysis also revealed transcripts that encode for proteins involved in apoptosis (programmed cell death), a process that also functions in internal defence responses. Identified transcripts included FAIM1 that regulates the extrinsic signalling pathway leading to apoptosis. Representatives of the intrinsic signalling pathway included serine protease HTRA2 (2 variants) that contribute to the loss of the mitochondrial transmembrane potential, one apoptosis-inducing factor (AIF) that causes DNA fragmentation and chromatin condensation, and factors that regulate these processes (Baculovirus IAP Repeat domain-containing caspase inhibitors, Bcl-Bax, PARP). Finally, eight variants of caspases that destroy critical cellular proteins were detected.
The search for response factors related to environmental change yielded two families of heat shock proteins that are induced by a variety of stressors: HSP70 (5 variants) and HSP90 (1 variant). Along with the HSPs, a heat shock factor (HSF, 1 variant) that regulates the expression of HSPs was identified. Additionally, our analysis yielded transcripts linked to invertebrate metabolism. These included one protein phosphatase (PP1) that controls for cellular processes linked to metabolism, gene transcription and translation, cell movement and apoptosis, and eight variants of ubiquitin, a regulatory protein that marks proteins, for instance, for degradation and recycling by proteasome and affects protein activity. Furthermore, transcripts encoding for ferritin (3 variants) that participates in iron transport and storage in invertebrates and protects organisms from iron-induced oxidative stress, and alcohol dehydrogenase (ADH, 1 variant) that catalyses the interconversion between alcohols and aldehydes or ketones were identified. Lastly, the search yielded transcripts linked to the estrogen-related receptor (ERR, 2 variants) and retinoid acid receptor (RXR, 3 variants) that regulate, for instance, oxidative metabolism, energy homeostasis, and imposex development in invertebrates.
Transcriptomic responses to experimental treatments
Plots of principal component (PC) scores for the transcriptome-wide expression profiles of individual snails did not reveal any grouping based on experimental treatments or outlier libraries that deviated from others (Fig. 2, Additional file 2). Heatmaps for the transcription of the annotated immune factors (see Additional file 1) across experimental treatments indicated some immune-elicitor specific responses (Fig. 3). For instance, injection with lyophilised E. coli cells increased the transcription of TLR, as well as of some components of the TLR signalling pathway (IκB, NF-κB), when compared to the levels in snails injected with physiological saline (comparison 1 in Fig. 3, Additional file 3). Additionally, exposure to E. coli increased the transcription of effectors representing antibacterial defence [Lymnaea-lysins (4 out of 6 individuals); comparison 2 in Fig. 3, Additional file 4] and PO/melanisation-type reaction (laccase; comparison 3 in Fig. 3, Additional file 5). Exposure of snails to lyophilised M. lysodeikticus cells activated the transcription of the lectin FREP (3 out of 6 individuals; comparison 4 in Fig. 3, Additional file 6), the cytokine IL-17 (comparison 5 in Fig. 3, Additional file 7), one component of the TLR signalling pathway (IκB; comparison 6 in Fig. 3, Additional file 3), and the effector laccase (comparison 3 in Fig. 3, Additional file 5). The injection of soluble extracts from trematode-infected snail tissue increased the transcription of laccase (comparison 7 in Fig. 3, Additional file 5) when compared to the snails challenged with extracts from healthy snail tissue. Also, IκB was activated in some individuals exposed to extracts of both healthy and trematode-infected snail tissue (comparison 8 in Fig. 3, Additional file 3). Additionally, wounding led to upregulation in the transcription of DUOXA (Additional file 8) that contributes to ROS production. These responses conform to the notions of their importance for invertebrate immune function.
Befitting the innate type immunity of invertebrates, the transcription of other immune factors did not show systematic differences among immune activation treatments (Additional files 3-9). However, the transcription of many factors showed high among-individual variation within treatment groups. This was the case for PRRs GNBP (Additional file 6), galectin (Additional file 6) and TLR (Additional file 3), as well as for some components of the TLR signalling pathway (Myd88, TRAF, NF-κB p65; Additional file 3), cytokines MIF and TNF (Additional file 7), several components of ROS production (NOX, DUOXA, GST, nucleoredoxin, NOS, peroxidasin; Additional file 8), the effector tyrosinase (PO/melanisation-type reaction; Additional file 5), factors regulating apoptosis (HTRA2, AIF, Bcl-Bax, PARP; Additional file 9), and all examined antibacterial defence factors (Additional file 4). A few immune factors showed high among-individual variation that was limited to certain immune activation treatments: lectin FREP (Additional file 6), IκB (TLR pathway; Additional file 3), DUOX (ROS production; Additional file 8), laccase (PO/melanisation-type reaction; Additional file 5) and the cytokine IL-17 (Additional file 7). Interestingly, within certain components of the immune system, individual snails expressed distinct gene transcripts, which suggests different defence strategies against pathogens. For instance, snails challenged with E. coli showed high transcription in two to four of the examined six antibacterial defence factors and the factors with high levels were different among individuals (Additional file 4). Most strikingly, individuals with the highest transcription of Lymnaea lysins showed the lowest transcription of LBP/BPI and vice versa (comparison 9 in Fig. 3). Similarly, snails exposed to trematode-infected snail tissue extracts showed high transcription in only one to two of the examined non-self recognition components (GNBP, FREP, galectins, TLR), and the expressed factors were different among individuals (areas 10 in Fig. 3).
The examined environmental changes apart from immune challenge affected several components of the snail immune system. The elevated temperature increased the transcription of the cytokine MIF (comparison 11 in Fig. 3, Additional file 7) when compared to untreated snails. Food deprivation, on the other hand, downregulated the transcription of an M-type lectin (comparison 12 in Fig. 3, Additional file 6), and glutathione S-transferase (GST, comparison 13 in Fig. 3, Additional file 8), but enhanced the transcription of MIF (comparison 14 in Fig. 3, Additional file 7), LBP/BPI (comparison 15 in Fig. 3, Additional file 4) and lysozyme (comparison 16 in Fig. 3, Additional file 4). Additionally, some of the factors that contribute to stress-responses, antioxidation and metabolism were found to be affected by the experimental treatments (Fig. 4). First, heat shock proteins (both HSP70 and HSP90; comparison 1 in Fig. 4, Additional file 10) and glutathione reductase (comparison 1 in Fig. 4, Additional file 11) showed reduced transcription under food deprivation. Exposure of snails to immune elicitors (injection with lyophilised E. coli and M. lysodeikticus cells, injection with healthy snail gonad) increased the transcription of ADH (comparison 2 in Fig. 4, Additional file 12) when compared to the snails injected with physiological saline. Similarly to the immune defence genes, high among-individual variation was seen in some annotated factors with a potential role as antioxidant enzymes [SOD (all treatments), MnSOD (some treatments), glutathione reductase (most immune activation treatments); Additional file 11] or in snail metabolism [PP1 (most treatments), ADH (most immune activations), RXR, (most treatments); Additional file 12].
This study aimed to produce a resource for future ecoimmunological research on the freshwater snail L. stagnalis by recording snail transcriptomic responses to immune challenges and environmental changes, with a focus on know immune factors. Earlier work on snail ecoimmunology has relied on quantifying a narrow subset of phenotypic immune defence traits. Those traits, however, respond differently to immune elicitors  and environmental conditions (e.g., [8, 52, 55, 60]), and they contribute differently to snail fitness . Additionally, contrary to the earlier expectation of simple innate-immune system in invertebrates (i.e., non-specific defence mechanisms) the recent development in comparative immunology and genomics/transcriptomics has revealed invertebrate immune systems to be highly complex, diverse, and also specific against different parasite types (e.g., [12,13,14,15,16]). Together these findings call for ecological experiments that quantify snail “immune phenotypes” across a wide range of immunological mechanisms. Only such studies may evaluate the role of snail immune function as a whole in ecological and evolutionary processes.
In this study, the samples for transcriptome sequencing were extracted from whole-body tissues of individual snails to avoid any bias in the recovery of immune and stress-response factors that are expressed in specific tissues or cell types, even if this may lead to reduced sensitivity to detect rare transcripts (see e.g., [38, 61]). Focusing on specific tissues/organs could be more sensitive compared to the whole-body approach when comparing experimental treatments, but we chose to use a measure that describes an individual’s overall immune activity. The reference transcriptome assembled from nine individuals exposed to different experimental treatments may not fully represent the genome content of L. stagnalis because mRNA-based approaches only capture genes that are actively expressed. BUSCO analysis, however, indicated representation of 98.8% of the universally shared metazoan genes in the assembly. Therefore, the reference transcriptome can be assumed to be highly comprehensive. Additionally, the use of a genetically diverse lab stock may have led to the capture of allelic variants of genes that do not all occur within single snail individuals, yet exist in natural populations where they may contribute differently to organismal fitness. This is supported by the high proportion (28.2%) of core BUSCO genes that were found with more than one copy.
Along with the aim to support future ecoimmunological research, the characterisation of immune genes of L. stagnalis broadens comparative immunology of aquatic pulmonate snails (Hygrophila), previously available only for a few species from the families Lymnaeaidae, Physideae and Planorbidae [38, 62, 63]. The organisation of antimicrobial defences in L. stagnalis is in line with hygrophilid snails that all show considerable diversity of antimicrobial proteins (LBP/BPI, LAAO, lysozyme, lymnaealysin, cytolytic β pore-forming toxins), contrasted by a modest number of AMP genes (a single family of 5 macin-type transcripts in L. stagnalis). This is remarkable because of the great diversity and numbers of AMP genes and gene families that are usually recorded from other organisms, including bivalve molluscs (e.g., ). Perhaps hygrophylid gastropods employ novel categories of AMPs that remain to be characterised, but it appears that these snails have a unique approach to dealing with invading microbes.
The recovery of L. stagnalis FREPs is consistent with the distribution of VIgLs as likely multimeric immune recognition factors across Gastropoda . The low number of FREPs in L. stagnalis is similar to that of the physid Physella acuta (Schultz et al., 2017) but differs from the planorbid B. glabrata that generates unique repertoires of multiple FREPs at the individual level . Furthermore, the transcriptomic responses of L. stagnalis to the introduction of extracts of healthy gonad tissue from other individuals (used as a control for injection with trematode-infected gonad extracts) was remarkably similar to those elicited by wounding alone. The lack of an immune response to allotypic tissue extracts suggests the absence of (acute) allograft rejection in L. stagnalis, just as it lacks from closely related planorbid snails [67, 68]. Together these findings suggest that snail immune function has evolved in a lineage-specific manner. However, our study supports the view that the main aspects of innate immunity are conserved throughout animal evolution . For example, lectins and GNBP receptors are available to detect pathogens and signal through the TLR pathway and cytokines to activate defence responses that encompass cellular (e.g., ROS production) and humoral (antimicrobial proteins) branches of the immune system.
The obtained detailed information regarding the molecular immunology of L. stagnalis can now be integrated into future ecoimmunological investigations. Studies on L. stagnalis typically measure two phenotypic immune defence traits from snail haemolymph, namely antibacterial activity and PO-like activity (e.g., [8, 52,53,54,55,56, 60]). Antibacterial activity is quantified as a reduction in optical density of a solution in which lyophilised E. coli cells are mixed with snail haemolymph (see [52, 59]). The current study identified several different effector proteins that may contribute to antibacterial responses. These included AMPs, LBP/BPI, LAAOs, lysozymes and cytolytic β pore-forming toxins. Therefore, the phenotypic antibacterial activity assay is likely to estimate the combined activity of several components. The other commonly studied phenotypic immune defence trait, PO-like activity, measures an increase in optical density of a solution in which PO enzymes from snail haemolymph oxidise the substrate L-dopa (see [52, 59]). This assay may also measure the combined activity of different factors (see ). Our study revealed two categories of PO enzymes, laccase and tyrosinase, suggesting that they may be the most important contributors to the melanisation-type reaction in L. stagnalis. Further tests at the phenotypic level using enzyme-specific substrates can estimate their relative importance in snail immune function (see ).
Visual examination of the transcription patterns of the annotated immune defence factors suggested different responses both against Gram-negative (E. coli) and Gram-positive (M. lysodeikticus) bacteria. Exposure to E. coli activated TLR pathway and the transcription of two effectors, Lymnaea-lysins (antibacterial activity) and laccase (PO activity). Patterns after exposure to M. lysodeikticus suggested upregulation in one FREP lectin and cytokine IL-17. The heat map analysis also indicated trematode-infected snail tissue extracts to activate the transcription of laccase. The apparent pathogen-specificity of some immune components is in line with previous findings in other invertebrates (e.g., [14, 17,18,19, 21]). Two of the annotated factors (laccase, IκB), however, showed broader responses to various immune elicitors. IκB (regulator of TLR signalling), showed increased transcription to all immune challenges except wounding, although this effect was not seen in all individuals. This finding suggests that the TLR pathway generally controls immune responses in L. stagnalis, similar as in other organisms [71, 72].
Based on the heat maps, manipulating the environmental conditions altered the transcription of some immune factors. Elevated temperature (25 °C vs. 20 °C) activated cytokine MIF and antibacterial lysozymes. This result differs from previous observations of high temperature reducing antibacterial and PO-like activity of snail haemolymph at the phenotypic level (e.g., [8, 53, 55]). In those studies, however, adverse effects were observed after a one-week exposure to high temperature . In this study, the exposure lasted for two days. Together these findings suggest that only prolonged exposure to high temperature negatively impacts snail immune function. Instead, short-term exposure to high temperature could have positive effects, for example, by increasing general performance through faster metabolic rate. In fact, high temperature initially increases snail growth rate and reproductive output . Food deprivation had both positive (MIF, LBP/BPI, lysozymes) and negative (M-type lectins, glutathione) effects on the transcription of immune factors. Previously, food deprivation has been reported to reduce snail immune activity at the phenotypic level, although also this effect depends on the duration of reduced food availability .
The transcription patterns did not indicate activation of the heat shock proteins at elevated temperature (25 °C vs. 20 °C). This finding suggests that the high-temperature treatment may not have been highly stressful to snails. Although 25 °C is above the thermal optimum of adult L. stagnalis snails (21–23 °C depending on the trait), it was still not close to the critically high temperature of this species . Additionally, the exposure time was short compared to earlier studies (see the previous paragraph) and to a typical 8.4-day summer heatwave the snails are exposed to in nature . Interestingly, transcription of the heat shock proteins was reduced under food deprivation, which indicates their general role in stress responses and that the resource level of snails could affect their ability to tolerate exposure to high temperature. Such interactive effects should be examined in future studies. Furthermore, exposure of snails to both Gram-negative and Gram-positive bacteria increased the transcription of ADH. This suggests that immune challenge can alter the metabolic rate of snails. However, although immune activation is often expected to increase metabolic activity (reviewed in ), other annotated factors with a potential role in metabolism were not affected by immune challenge treatments.
Our study was conducted using snails that originated from a genetically diverse laboratory stock that was formed by initially combining individuals from different natural populations . Genetically variable snails were used to identify genes that show high among-individual variation in transcription. Such variation may arise, for example, from differences in the genetic background (expressed within and/or among populations) and physiological condition of snails. Identifying genes that vary in their transcription among individuals is very important for ecological and evolutionary studies that focus on understanding the sources and consequences of trait variation. Understanding among-individual variation can be critical because trait evolution may depend more strongly on variation in gene expression than differences in protein-coding sequences [75, 76]. Moreover, the suitability of considering transcription as a “trait”, and that it evolves in natural populations has been demonstrated in yeast (e.g., [77, 78]), fruit fly (e.g., [79, 80]), and fish (e.g., [81,82,83]). Among-individual variation in transcription is, however, easily overlooked in experiments that use genetically as homogeneous individuals as possible (e.g., one clone/inbred strain).
The heat map analysis suggested high variation in transcription among snail individuals in several components of the snail immune system from pathogen recognition (GNBP, lectins, TLR), to signalling/regulation (TLR pathway, cytokines, ROS production, apoptosis) and effectors (laccase, tyrosinase, antibacterial proteins). In several cases, such variation was seen even among untreated snails (e.g., PGRP, TLR, Myd88, MIF, LBP/BPI, LAAOs, lysozymes, NOX, DUOX, TNF receptor). Two potential explanations are proposed for the latter finding. First, control snails may have harboured obscure infections from opportunistic microbes or viruses that activate snail immune function , leading to variation in transcription profiles. Second, these components of the snail immune system are expressed at a constant level within each individual, but the actual level varies among individuals , perhaps due to differences among alleles that they carry. This would fit with innate type immunity, the expression levels reflecting the immediate capacity to respond to invading pathogens. Nevertheless, all the immune factors that showed high among-individual variation in transcription are potentially interesting targets for future ecoimmunological research. To our knowledge, such variation has been utilised in ecoimmunology only in bumblebee [26, 27, 85, 86]. In that species, condition dependence of immune defence and genetic specificity determining the outcome of a host-parasite interaction have been examined at the transcription level. Similar variation may underlie the observation that parasite exposure does not yield infection in all snails of a generally susceptible strain in the B. glabrata-S. mansoni interaction .
Despite high among-individual variation in the transcription of many components of the snail immune system, our data did not reveal any individuals with generally very low or very high immune activity across the annotated factors. However, individuals within some of the immune activation treatments showed investment in different components of defence. For example, each snail exposed to E. coli showed high transcription in two to four of the examined six antibacterial defence factors. Interestingly, which of these six factors showed high expression levels varied among individuals. Similarly, in snails exposed to trematode-infected snail tissue extracts, each individual typically showed high transcription in only one, but not the same of the examined non-self recognition components (GNBP, FREP, galectins, TLR). Yet other immune properties vary among individuals both within field- and laboratory-maintained L. stagnalis populations. Van der Knaap et al. [88, 89] described two genetically-defined types of L. stagnalis, with “type I” snails having a particular type of lectin-mediated immune recognition that was absent in the majority of snails, designated as type II. These observations suggest that snails may show different immune defence strategies that lead to diverse responses even under a similar immune challenge. Nevertheless, the data presented in this article does not allow estimating whether the observed among-individual variation in transcription arises from differences expressed within and/or among populations. Both sources of variation are important to consider, which calls for studies investigating variation in the transcription of different immunological mechanisms and their contribution to snail performance (e.g., survival, fecundity) at both within- and among-population scales. Such studies would, for example, allow examining evolution as a within-population process (e.g., natural selection) and an outcome such as differentiation in immune activity among populations that have experienced different pathogen histories.
Our results confirm uniform aspects of gastropod/molluscan immunity, but also demonstrate apparent differences between L. stagnalis and some previously examined taxa. For instance, PRRs and lectins are present to detect pathogens and signal through the TLR pathway and cytokines to activate both cellular and humoral defence mechanisms. However, contrary to the high number of AMPs in other organisms, including molluscs, only a modest number of macin type AMPs was identified in L. stagnalis. Similarly, the low number of FREPs is comparable to P. acuta but differs from B. glabrata that shows high FREP diversity. Our data also indicate that various defence factors can contribute to the phenotypic immunological assays (antibacterial activity and PO-like activity of haemolymph) commonly used in earlier ecoimmunological research. Additionally, immune activation treatments (bacteria, trematodes) and manipulations of environmental conditions (ambient temperature, resource availability) revealed factors that contribute to snail immune responses (TLR pathway, FREP, cytokines, antibacterial defence, PO activity) and dependence of immune activity on environmental variation (cytokines, antibacterial defence, ROS production), as well as how these responses may depend on general stress reactions (HSP). Interestingly, many examined immune factors (PRRs, lectins, TLR pathway, cytokines, ROS production, apoptosis, PO enzymes, antibacterial proteins) showed considerable among-individual variation in the genetically diverse experimental snails used in this study. In addition to the components responding to immune elicitors, these factors can be important determinants of fitness variation within and between natural snail populations and should be included in future ecoimmunological studies.
Study animals and experimental design
This study employed adult L. stagnalis snails (N = 48, shell length 24.6–30.5 mm) from an F2 generation of a laboratory stock population that was generated by interbreeding snails that originated from seven different natural populations in Switzerland (see ). Interbreeding was done to increase genetic variation among the experimental snails relative to natural populations that may show low genetic diversity . Higher genetic diversity among experimental snails increased the potential to identify immune genes that show variation in their transcription among individuals depending, for example, on their genotype and allelic variation.
Before the experiment, experimental snails were placed individually in perforated 0.2 L plastic cups and distributed in four water baths each containing 40 L of aged tap water at 20 °C, with twelve snails per water bath. Snails were fed with fresh lettuce ad libitum, and half of the water in each water bath was changed every second day. This set-up allowed water changes with minimal disturbance to the snails. Snails were acclimated to these conditions for 4 (treatments including an environmental change) to 6 days (treatments including an immune challenge) before exposure to experimental manipulations (see the next two paragraphs for the description of the treatments). The acclimation period differed between treatment types to allow simultaneous sampling of snail tissues for RNA extractions (i.e., duration of the treatments was not the same).
In the experiment, the snails were randomly assigned to nine experimental treatments (Table 1). In seven treatments, snails’ responses to immune elicitors, including the effects of handling that were necessary for immune challenge treatments, were tested. Snails in treatment one (untreated controls; N = 5) were maintained as during the acclimation period. Snails in treatment two (anaesthetised controls, N = 5) were anaesthetised by placing them into 2% diethyl ether for two and a half minutes. These snails were used to examine the effects of snail handling when compared to untreated controls (treatment 1 above) because anaesthesia was needed in following immune challenge treatments. In treatment three (wounding; N = 6), each anaesthetised snail was injected with 170 μL of snail saline (see ) in the head-foot using a syringe and needle (0.45 mm × 16 mm BD Microlane™ 3, New Jersey, USA). This treatment examined the effect of wounding in activating immune defence when compared to anaesthesia without exposure to immune elicitors (treatment 2 above). The sham injection also served as a control for the treatments involving bacterial injections (described in the next paragraph).
In treatments four and five (bacterial injections; N = 6 for each), anaesthetised snails were injected as above with 0.5 mg/mL lyophilised E. coli (Sigma-Aldrich, Steinheim, Germany; product number EC11303) or M. lysodeikticus (synonym Micrococcus luteus; Sigma-Aldrich, Steinheim, Germany; product number M3770) cells in snail saline, respectively. Both Gram-negative (E. coli) and Gram-positive (M. lysodeikticus) bacteria were used because lipopolysaccharides of Gram-negative bacteria and peptidoglycans of Gram-positive bacteria are recognised by different receptors that activate pathogen-specific immune responses [14, 21, 92,93,94].
In treatments six and seven (snail/trematode-tissue injections; N = 6 for each), anaesthetised snails were injected as above with soluble extracts from healthy or trematode-infected gonads (ovotestes) dissected from L. stagnalis snails collected from a pond in Winterthur, Switzerland (47° 28′N, 8° 43′E). In donor snails, parasite infections were evident from the shedding of strigeid and echinostome-type cercariae. Both healthy and infected gonads were disrupted in 1.5 mL Eppendorf tubes using plastic pestles, adding 1 mL of snail saline per 100 mg of tissue. Samples were vortexed and pelleted briefly (2000×g for 5 s). The supernatant was collected from three uninfected and three infected snails and the samples from each snail type were pooled for injections. The treatment using the extract from healthy gonads served as a control for injection with tissue extract from the trematode-infected snail. After injections, all experimental snails were placed back in their original cups with lettuce for 6 hours to allow their immune systems to respond to immune elicitors.
In addition to immune challenges, snails’ responses to environmental changes were tested in two treatments. Snails in treatment eight (elevated temperature; N = 4) were placed in water baths where the water temperature was raised to 25 °C during several hours and maintained at an elevated level for 2 days (see [55, 56]). After the transfer, additional snails were added to all water baths to keep the experimental snails continuously at a density of twelve snails per water bath (data were not collected from these extra snails). Snails in treatment nine (food deprivation; N = 4) were maintained without food for 2 days. Untreated snails (treatment 1 above) served as controls for these treatments. All snails survived the experimental treatments and displayed normal motility and activities during the experiment.
RNA extraction, library preparation and sequencing
At the end of the experimental treatments, each snail was removed from its shell, and whole-body soft tissues were submerged in self-made RNAlater solution (see protocol at https://www.protocols.io/view/RNAlater-Recipe-c56y9d) in which they were cut into small pieces. The cut tissues were stored in 10 mL of RNAlater at − 20 °C. Within the next 15 days, all samples were disrupted under liquid nitrogen using mortar and pestle. Total RNA was extracted from each homogenised tissue sample (71–102 mg per snail) using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions except for conducting RNA wash with 75% EtOH three times. Residual genomic DNA was removed from samples using Turbo DNA-free kit (Ambion, Austin, TX, USA) according to the manufacturer’s instructions. RNA quantity was measured using Qubit 2.0 Fluorometer (Invitrogen, Carlsbad, CA, USA), and sample quality and purity were verified using 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA; RNA Nano chips) and NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA), respectively. Complementary DNA (cDNA) libraries were constructed for each snail individual using TruSeq Stranded mRNA Sample Preparation Kit (Illumina, San Diego, CA, USA) according to the manufacturer’s instructions. All libraries were sequenced on Illumina NextSeq 500 platform using paired-end reads with 81 nt read length at the Genomics Facility Basel.
RNA-seq data handling and transcriptome assembly
Raw Illumina reads from all libraries were subjected to adaptor trimming and quality filtering. Reads with Illumina TruSeq adaptors were trimmed using Cutadapt v1.5 (min. Overlap 30 nt and max. Allowed error rate 0.05). Reads were then quality filtered using PRINSEQ-lite v0.20.4 [min. Read length: 50 nt, GC range: 10–90%, mean Q ≥ 5, trim base left and right with Q < 10, no ambiguous nucleotides allowed, trim poly A/T tail (min. 5) from both sequences, and remove duplicate reads]. Finally, Phix sequences with GenBank genome reference NC_001422.1, as well as SSU and LSU rRNA sequences with SILVA (Release 111) reference, were removed using BBDuk v2015.08.21. To generate a L. stagnalis de novo reference transcriptome, the cleaned reads from nine libraries (one randomly chosen library per experimental treatment) were combined and assembled with trinity v2.0.6 (parameters: min contig length: 100, min glue: 4, min kmer cov: 4, group pairs distance: 300, path reinforcement distance: 85, normalise reads, normalise max read cov: 50, and normalise by read set). Likely contaminant transcripts from other organisms were removed by mapping all transcripts on the L. stagnalis reference genome (GenBank accession number GCA_900036025.1) and removing those that lacked a significant similarity hit in the genome if they were not previously characterised from L. stagnalis. Completeness of the produced reference transcriptome was assessed by examining the detection of core BUSCO genes [95, 96] from metazoans (978 genes).
Annotation of immune-, stress- and metabolism-related factors
Amino acid sequences of previously characterised proteins/genes relevant for immune function, stress responses and metabolism primarily from molluscs (especially gastropods; e.g., [38, 97]), but also from other organisms, were collected from GenBank (Additional file 13) and used in BLAST (v2.6.0+) similarity searches to identify orthologs in the L. stagnalis reference transcriptome (226,116 contigs). Components of the immune system queried included non-self recognition through pathogen-recognition receptors (PRRs; all abbreviations are explained at the end of the article), Toll-like receptor (TLR) pathway, cytokines, antibacterial peptides and proteins, production of reactive oxygen species (ROS), phenoloxidase (PO)/melanisation-type defence, and apoptosis. Components of stress responses included cell protection and survival, oxidative stress and antioxidant enzymes. Targeted components linked to metabolism included regulatory proteins, proteins related to the transport of elements and compounds as well as a selection of receptors.
From the statistically best Blast hits (E-value ≤0.05), transcripts that were at least 60% of the length of the shortest reference sequence and showed minimum 40% similarity were used to identify ORFs (EMBOSS v220.127.116.11) that were between a start and a stop codon. ORFs that were at least 60% of the length of the shortest reference sequence (absolute minimum length: 30 amino acids) were examined for the presence of a signal peptide (SignalP v4.1) and the domain structure (SMART, PFAM, NCBI). For identification of fibrinogen-related proteins (FREPs), sequences that were computationally identified to contain fibrinogen-like (FBG) domains were manually inspected for upstream immunoglobulin superfamily (IgSF) domains because the diffuse sequence motifs of invertebrate IgSF domains frequently challenge automated detection.
Transcript expression level analysis
To examine variation in the expression levels of different transcripts among experimental snails, sequencing reads from each library were first indexed to transcripts in the reference transcriptome using Kallisto v.0.44.0 with 100 bootstraps . Sleuth was then used for downstream processing . The transcriptome-wide expression profiles in different libraries were visualised using their PC scores obtained from a principal component analysis (PCA, the first 5 PCs were used) using the internal normalisation in Sleuth. Very high variation in the expression profiles among snail individuals even within experimental treatments (see the Results section) violated the assumptions of formal statistical tests for differential expression of individual transcripts among treatments, thus negating their use (see ). Additionally, we are not aware of a statistical test that compares the amount of variation in transcription over biological replicates across transcripts within treatments. Therefore, both the within- and among-treatment variation in expression levels of individual transcripts that represented the annotated factors were examined visually using heatmaps that presented all experimental snails. Signal strength was calculated in units of transcripts per million (TPM). This measure calculates the proportion of counts per base for each transcript in the whole data set (multiplied by million).
Availability of data and materials
The raw sequencing reads are deposited in the NCBI Sequence Read Archive (accession number PRJNA664475). The reference transcriptome is available in https://doi.org/10.5281/zenodo.4044169.
Oxidase maturation factor
Fas apoptotic inhibitory molecule 1
Gram-negative bacteria binding-protein
Heat shock factor
70 kilodalton heat shock protein
90 kilodalton heat shock protein
Inhibitor of apoptosis proteins
Inhibitor of NF-κB kinase
Inhibitor of nuclear factor-κB
L-amino acid oxidase
Lipopolysaccharide-binding protein/bactericidal permeability-increasing protein
PS-induced TNF-α factor
Macrophage migration inhibitory factor
Macrophage expressed gene
Myeloid differentiation primary response protein 88
Nicotinamide adenine dinucleotide phosphate
NF- B essential modulator
- NF- B:
Nuclear factor B
Nitric oxide synthase
Open reading frame
Poly (ADP-ribose) polymerase
Principal component analysis
Peptidoglycan recognition protein
Retinoid acid receptor
Sterile alpha and armadillo motif-containing protein
- SOD and MnSOD:
TAK1 binding protein
Transforming growth factor β-activated kinase-1
Tumor necrosis factor
Transcripts per million
Tumor necrosis factor receptor-associated factor
Variable immunoglobulin and lectin-domain-containing molecules
Schmid-Hempel P. Evolutionary parasitology: the integrated study of infections, immunology, ecology, and genetics. New York: Oxford University Press; 2011.
Zuk M, Stoehr AM. Immune defense and host life history. Am Nat. 2002;160:S9–S22.
Nunn CL, Lindenfors P, Pursall ER, Rolff J. On sexual dimorphism in immune function. Philos T R Soc B. 2009;364(1513):61–9.
Hamilton WD, Zuk M. Heritable true fitness and bright birds: a role for parasites. Science. 1982;218(4570):384–7.
Rantala MJ, Koskimaki J, Taskinen J, Tynkkynen K, Suhonen J. Immunocompetence, developmental stability and wingspot size in the damselfly Calopteryx splendens L. Proc R Soc B. 2000;267(1460):2453–7.
Siva-Jothy MT, Thompson JJW. Short-term nutrient deprivation affects immune function. Physiol Entomol. 2002;27:206–12.
Wilson K, Thomas MB, Blanford S, Doggett M, Simpson SJ, Moore SL. Coping with crowds: density-dependent disease resistance in desert locusts. P Natl Acad Sci USA. 2002;99(8):5471–5.
Leicht K, Seppälä K, Seppälä O. Potential for adaptation to climate change: family-level variation in fitness-related traits and their responses to heat waves in a snail population. BMC Evol Biol. 2017;17:140.
Demas GE, Nelson RJ. Ecoimmunology. New York: Oxford University Press; 2012.
Cotter SC, Kruuk LEB, Wilson K. Costs of resistance: genetic correlations and potential trade-offs in an insect immune system. J Evol Biol. 2004;17:421–9.
Schwarzenbach GA, Hosken DJ, Ward PI. Sex and immunity in the yellow dung fly Scathophaga stercoraria. J Evol Biol. 2005;18:455–63.
Loker ES, Adema CM, Zhang SM, Kepler TB. Invertebrate immune systems - not homogeneous, not simple, not well understood. Immunol Rev. 2004;198:10–24.
Ghosh J, Lun CM, Majeske AJ, Sacchi S, Schrankel CS, Smith LC. Invertebrate immune diversity. Dev Comp Immunol. 2011;35(9):959–74.
Deleury E, Dubreuil G, Elangovan N, Wajnberg E, Reichhart JM, Gourbal B, et al. Specific versus non-specific immune responses in an invertebrate species evidenced by a comparative de novo sequencing study. PLoS One. 2012;7(3):e32512.
Dheilly NM, Adema C, Raftos DA, Gourbal B, Grunau C, Du Pasquier L. No more non-model species: the promise of next generation sequencing for comparative immunology. Dev Comp Immunol. 2014;45(1):56–66.
Doublet V, Poeschl Y, Gogol-Döring A, Alaux C, Annoscia D, Aurori C, et al. Unity in defence: honeybee workers exhibit conserved molecular responses to diverse pathogens. BMC Genomics. 2017;18:207.
Hoffmann JA, Reichhart JM. Drosophila innate immunity: an evolutionary perspective. Nat Immunol. 2002;3(2):121–6.
Royet J, Reichhart JM, Hoffmann JA. Sensing and signaling during infection in Drosophila. Curr Opin Immunol. 2005;17(1):11–7.
Hetru C, Hoffmann JA. NF-kappa B in the immune response of Drosophila. Csh Perspect Biol. 2009;1(6):a000232.
Pedersen AB, Babayan SA. Wild immunology. Mol Ecol. 2011;20(5):872–80.
Hanelt B, Lun CM, Adema CM. Comparative ORESTES-sampling of transcriptomes of immune-challenged Biomphalaria glabrata snails. J Invertebr Pathol. 2008;99(2):192–203.
Adema CM, Hanington PC, Lun CM, Rosenberg GH, Aragon AD, Stout BA, et al. Differential transcriptomic responses of Biomphalaria glabrata (Gastropoda, Mollusca) to bacteria and metazoan parasites, Schistosoma mansoni and Echinostoma paraensei (Digenea, Platyhelminthes). Mol Immunol. 2010;47(4):849–60.
Zhang LL, Li L, Guo XM, Litman GW, Dishaw LJ, Zhang GF. Massive expansion and functional divergence of innate immune genes in a protostome. Sci Rep-Uk. 2015;5:8693.
Feder D, Mello CB, Garcia ES, Azambuja P. Immune responses in Rhodnius prolixus: influence of nutrition and ecdysone. J Insect Physiol. 1997;43:513–9.
Lee KP, Cory JS, Wilson K, Raubenheimer D, Simpson SJ. Flexible diet choice offsets protein costs of pathogen resistance in a caterpillar. Proc R Soc B. 2006;273(1588):823–9.
Barribeau SM, Sadd B, du Plessis L, Schmid-Hempel P. Gene expression differences underlying genotype-by-genotype specificity in a host-parasite system. P Natl Acad Sci USA. 2014;111(9):3496–501.
Brunner FS, Schmid-Hempel P, Barribeau SM. Protein-poor diet reduces host-specific immune gene expression in Bombus terrestris. Proc R Soc B. 2014;281(1786):20140128.
Ferro K, Ferro D, Corrà F, Bakiu R, Santovito G, Kurtz J. Cu,Zn superoxide dismutase genes in Tribolium castaneum: evolution, molecular characterisation, and gene expression during immune priming. Front Immunol. 2017;8:1811.
Greenwood JM, Milutinović B, Peuß R, Behrens S, Esser D, Rosenstiel P, et al. Oral immune priming with Bacillus thuringiensis induces a shift in the gene expression of Tribolium castaneum larvae. BMC Genomics. 2017;18:329.
Lindberg DR, Ponder WF, Haszprunar G. The mollusca: relationships and patterns from their first half-billion years. In: Cracraft J, editor. Assemb Tree Life. Cary, North Carolina: Oxford University Press; 2004. p. 252–78.
Moore JD, Finley CA, Robbins TT, Friedman CS. Withering syndrome and restoration of southern California abalone populations. Cal Coop Ocean Fish. 2002;43:112–7.
Prince JS. Presumptive alphavirus in the gastropod mollusc, Aplyka californica. B Mar Sci. 2003;73(3):673–7.
Travers MA, Barbou A, Le Goïc N, Huchette S, Paillard C, Koken M. Construction of a stable GFP-tagged Vibrio harveyi strain for bacterial dynamics analysis of abalone infection. FEMS Microbiol Lett. 2008;289(1):34–40.
Cribb TH, Bray RA, Olson PD, Littlewood DTJ. Life cycle evolution in the Digenea: a new perspective from phylogeny. Adv Parasitol Vol 54. 2003;54:197–254.
Olson PD, Cribb TH, Tkach VV, Bray RA, Littlewood DTJ. Phylogeny and classification of the Digenea (Platyhelminthes : Trematoda). Int J Parasitol. 2003;33(7):733–55.
Bayne CJ. Successful parasitism of vector snail Biomphalaria glabrata by the human blood fluke (trematode) Schistosoma mansoni: a 2009 assessment. Mol Biochem Parasitol. 2009;165(1):8–18.
Coustau C, Gourbal B, Duval D, Yoshino TP, Adema CM, Mitta G. Advances in gastropod immunity from the study of the interaction between the snail Biomphalaria glabrata and its parasites: a review of research progress over the last decade. Fish Shellfish Immunol. 2015;46(1):5–16.
Adema CM, Hillier LW, Jones CS, Loker ES, Knight M, Minx P, et al. Whole genome analysis of a schistosomiasis-transmitting freshwater snail (vol 8, 15451, 2017). Nat Commun. 2017;8:15451.
Caron Y, Martens K, Lempereur L, Saegerman C, Losson B. New insight in lymnaeid snails (Mollusca, Gastropoda) as intermediate hosts of Fasciola hepatica (Trematoda, Digenea) in Belgium and Luxembourg. Parasite Vector. 2014;7:66.
Selbach C, Soldánová M, Georgieva S, Kostadinova A, Sures B. Integrative taxonomic approach to the cryptic diversity of Diplostomum spp. in lymnaeid snails from Europe with a focus on the 'Diplostomum mergi' species complex. Parasite Vector. 2015;8:300.
Horák P, Mikeš L, Lichtenbergová L, Skala V, Soldánová M, Brant SV. Avian schistosomes and outbreaks of cercarial dermatitis. Clin Microbiol Rev. 2015;28(1):165–90.
van der Knaap WPW, Sminia T, Kroese FGM, Dikkeboom R. Elimination of bacteria from the circulation of the pond snail Lymnaea stagnalis. Dev Comp Immunol. 1981;5(1):21–32.
Vanderknaap WPW, Adema CM, Sminia T. Invertebrate blood cells: morphological and functional aspects of the hemocytes in the pond snail Lymnaea stagnalis. Comp Haematol Int. 1993;3(1):20–6.
Adema CM, Vandeutekommulder EC, Vanderknaap WPW, Sminia T. Schistosomicidal activities of Lymnaea stagnalis hemocytes: the role of oxygen radicals. Parasitology. 1994;109:479–85.
Horák P, Deme R. Lectins and saccharides in Lymnaea stagnalis haemocyte recognition. Comp Haematol Int. 1998;8(4):210–8.
Lacchini AH, Davies AJ, Mackintosh D, Walker AJ. β-1,3-glucan modulates PKC signalling in Lymnaea stagnalis defence cells: a role for PKC in H2O2 production and downstream ERK activation. J Exp Biol. 2006;209(24):4829–40.
Wright B, Lacchini AH, Davies AJ, Walker AJ. Regulation of nitric oxide production in snail (Lymnaea stagnalis) defence cells: a role for PKC and ERK signalling pathways. Biol Cell. 2006;98(5):265–78.
Walker AJ, Lacchini AH, Sealey KL, Mackintosh D, Davies AJ. Spreading by snail (Lymnaea stagnalis) defence cells is regulated through integrated PKC, FAK and Src signalling. Cell Tissue Res. 2010;341(1):131–45.
Skála V, Walker AJ, Horák P. Snail defence responses to parasite infection: the Lymnaea stagnalis-Trichobilharzia szidati model. Dev Comp Immunol. 2020;102:103464.
Davison A, McDowell GS, Holden JM, Johnson HF, Koutsovoulos GD, Liu MM, et al. Formin is associated with left-right asymmetry in the pond snail and the frog. Curr Biol. 2016;26(5):654–60.
Fodor I, Hussein AAA, Benjamin PR, Koene JM, Pirger Z. The natural history of model organisms: the unlimited potential of the great pond snail, Lymnaea stagnalis. Elife. 2020;9:e56962.
Seppälä O, Jokela J. Maintenance of genetic variation in immune defense of a freshwater snail: role of environmental heterogeneity. Evolution. 2010;64:2397–407.
Seppälä O, Jokela J. Immune defence under extreme ambient temperature. Biol Lett. 2011;7:119–22.
Langeloh L, Behrmann-Godel J, Seppälä O. Natural selection on immune defense: a field experiment. Evolution. 2017;71(2):227–37.
Salo T, Stamm C, Burdon FJ, Räsänen K, Seppälä O. Resilience to heat waves in the aquatic snail Lymnaea stagnalis: additive and interactive effects with micropollutants. Freshw Biol. 2017;62(11):1831–46.
Salo T, Kropf T, Burdon FJ, Seppälä O. Diurnal variation around an optimum and near-critically high temperature does not alter the performance of an ectothermic aquatic grazer. Ecol Evol. 2019;9(20):11695–706.
Galinier R, Portela J, Moné Y, Allienne JF, Henri H, Delbecq S, et al. Biomphalysin, a new β pore-forming toxin involved in Biomphalaria glabrata immune defense against Schistosoma mansoni. PLoS Pathog. 2013;9(3):e1003216.
Lassalle D, Tetreau G, Pinaud S, Galinier R, Crickmore N, Gourbal B, et al. Glabralysins, potential new β-pore-forming toxin family members from the schistosomiasis vector snail Biomphalaria glabrata. Genes-Basel. 2020;11(1):65.
Seppälä O, Leicht K. Activation of the immune defence of the freshwater snail Lymnaea stagnalis by different immune elicitors. J Exp Biol. 2013;216:2902–7.
Leicht K, Jokela J, Seppälä O. An experimental heat wave changes immune defense and life history traits in a freshwater snail. Ecol Evol. 2013;3(15):4861–71.
Baron OL, Deleury E, Reichhart JM, Coustau C. The LBP/BPI multigenic family in invertebrates: evolutionary history and evidences of specialisation in mollusks. Dev Comp Immunol. 2016;57:20–30.
Schultz JH, Adema CM. Comparative immunogenomics of molluscs. Dev Comp Immunol. 2017;75:3–15.
Alba A, Tetreau G, Chaparro C, Sánchez J, Vázquez AA, Gourbal B. Natural resistance to Fasciola hepatica (Trematoda) in Pseudosuccinea columella snails: a review from literature and insights from comparative "omic" analyses. Dev Comp Immunol. 2019;101:103463.
Mitta G, Vandenbulcke F, Roch P. Original involvement of antimicrobial peptides in mussel innate immunity. FEBS Lett. 2000;486:185–90.
Gorbushin AM. Derivatives of the lectin complement pathway in Lophotrochozoa. Dev Comp Immunol. 2019;94:35–58.
Zhang SM, Adema CM, Kepler TB, Loker ES. Diversification of Ig superfamily genes in an invertebrate. Science. 2004;305(5681):251–4.
Sullivan JT, Brammer SR, Hargraves CD, Owens BS. Heterotopic heart-transplants in Biomphalaria glabrata (Mollusca, Pulmonata).: fate of xenografts from seven pulmonate genera. Invertebr Biol. 1995;114(2):151–60.
Sullivan JT, Galvan AG, Lares RR. Comparison of several types of allografts in Biomphalaria glabrata (Mollusca : Pulmonata). J Invertebr Pathol. 1998;71(1):1–8.
Buchmann K. Evolution of innate immunity: clues from invertebrates via fish to mammals. Front Immunol. 2014;5:459.
Le Clec'h W, Anderson TJC, Chevalier FD. Characterisation of hemolymph phenoloxidase activity in two Biomphalaria snail species and impact of Schistosoma mansoni infection. Parasite Vector. 2016;9:32.
Pila EA, Tarrabain M, Kabore AL, Hanington PC. A novel toll-like receptor (TLR) influences compatibility between the gastropod Biomphalaria glabrata, and the digenean trematode Schistosoma mansoni. PLoS Pathog. 2016;12(3):e1005513.
Nie L, Cai SY, Shao JZ, Chen J. Toll-like receptors, associated biological roles, and signaling networks in non-mammals. Front Immunol. 2018;9.
Meehl GA, Tebaldi C. More intense, more frequent, and longer lasting heat waves in the 21st century. Science. 2004;305:994–7.
Demas GE, Chefer V, Talan MI, Nelson RJ. Metabolic costs of mounting an antigen-stimulated immune response in adult and aged C57BL/6J mice. Am J Physiol. 1997;273(5):R1631–7.
Enard W, Khaitovich P, Klose J, Zöllner S, Heissig F, Giavalisco P, et al. Intra- and interspecific variation in primate gene expression patterns. Science. 2002;296(5566):340–3.
Gilad Y, Oshlack A, Smyth GK, Speed TP, White KP. Expression profiling in primates reveals a rapid evolution of human transcription factors. Nature. 2006;440(7081):242–5.
Ferea TL, Botstein D, Brown PO, Rosenzweig RF. Systematic changes in gene expression patterns following adaptive evolution in yeast. P Natl Acad Sci USA. 1999;96(17):9721–6.
Brem RB, Yvert G, Clinton R, Kruglyak L. Genetic dissection of transcriptional regulation in budding yeast. Science. 2002;296(5568):752–5.
Jin W, Riley RM, Wolfinger RD, White KP, Passador-Gurgel G, Gibson G. The contributions of sex, genotype and age to transcriptional variance in Drosophila melanogaster. Nat Genet. 2001;29(4):389–95.
Gibson G, Riley-Berger R, Harshman L, Kopp A, Vacha S, Nuzhdin S, et al. Extensive sex-specific nonadditivity of gene expression in Drosophila melanogaster. Genetics. 2004;167(4):1791–9.
Oleksiak MF, Roach JL, Crawford DL. Natural variation in cardiac metabolism and gene expression in Fundulus heteroclitus. Nat Genet. 2005;37(1):67–72.
Whitehead A, Crawford DL. Neutral and adaptive variation in gene expression. P Natl Acad Sci USA. 2006;103(14):5425–30.
Leder EH, McCairns RJS, Leinonen T, Cano JM, Viitaniemi HM, Nikinmaa M, et al. The evolution and adaptive potential of transcriptional variation in sticklebacks—signatures of selection and widespread heritability. Mol Biol Evol. 2015;32(3):674–89.
Bender RC, Goodall CP, Blouin MS, Bayne CJ. Variation in expression of Biomphalaria glabrata SOD1: a potential controlling factor in susceptibility/resistance to Schistosoma mansoni. Dev Comp Immunol. 2007;31(9):874–8.
Brunner FS, Schmid-Hempel P, Barribeau SM. Immune gene expression in Bombus terrestris: signatures of infection despite strong variation among populations, colonies, and sister workers. PLoS One. 2013;8(7):e68181.
Marxer M, Barribeau S, Schmid-Hempel P. Experimental evolution of a trypanosome parasite of bumblebees and its implications for infection success and host immune response. Evol Biol. 2016;43(2):160–70.
Théron A, Coustau C. Are Biomphalaria snails resistant to Schistosoma mansoni? J Helminthol. 2005;79(3):187–91.
van der Knaap WPW, Doderer A, Boerrigter-barendsen LH, Sminia T. Some properties of an agglutinin in the hemolymph of the pond snail Lymnaea stagnalis. Biol Bull. 1982;162(3):404–12.
van der Knaap WPW, Sminia T, Schutte R, Boerrigter-barendsen LH. Cytophilic receptors for foreignness and some factors which influence phagocytosis by invertebrate leukocytes: in vitro phagocytosis by amebocytes of the snail Lymnaea stagnalis. Immunology. 1983;48(2):377–83.
Kopp KC, Wolff K, Jokela J. Natural range expansion and human-assisted introduction leave different genetic signatures in a hermaphroditic freshwater snail. Evol Ecol. 2012;26(3):483–98.
Standen NB. Voltage-clamp studies of calcium inward current in an identified snail euron: comparison with sodium inward current. J Physiol-London. 1975;249(2):253–68.
Poltorak A, He XL, Smirnova I, Liu MY, Van Huffel C, Du X, et al. Defective LPS signaling in C3H/HeJ and C57BL/10ScCr mice: mutations in Tlr4 gene. Science. 1998;282(5396):2085–8.
Takehana A, Katsuyama T, Yano T, Oshima Y, Takada H, Aigaki T, et al. Overexpression of a pattern-recognition receptor, peptidoglycan-recognition protein-LE, activates imd/relish-mediated antibacterial defense and the prophenoloxidase cascade in Drosophila larvae. P Natl Acad Sci USA. 2002;99(21):13705–10.
Montminy SW, Khan N, McGrath S, Walkowicz MJ, Sharp F, Conlon JE, et al. Virulence factors of Yersinia pestis are overcome by a strong lipopolysaccharide response. Nat Immunol. 2006;7(10):1066–73.
Simão FA, Waterhouse RM, Ioannidis P, Kriventseva EV, Zdobnov EM. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs. Bioinformatics. 2015;31(19):3210–2.
Waterhouse RM, Seppey M, Simão FA, Manni M, Ioannidis P, Klioutchnikov G, et al. BUSCO applications from quality assessments to gene prediction and phylogenomics. Mol Biol Evol. 2018;35(3):543–8.
Gorbushin AM, Borisova EA. Lectin-like molecules in transcriptome of Littorina littorea hemocytes. Dev Comp Immunol. 2015;48(1):210–20.
Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7.
Pimentel H, Bray NL, Puente S, Melsted P, Pachter L. Differential analysis of RNA-seq incorporating quantification uncertainty. Nat Methods. 2017;14(7):687.
Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, et al. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? Rna. 2016;22(6):839–51.
We thank B. Hannah, S. Kobel, A. Minder, T. Torossi, and N. Zemp for help and guidance with the lab work and the analyses. We are grateful to M.-A. Coutellec, A. Dennis, B. Feldmeyer, J. Jokela, E. S. Loker, M. Neiman, and L. Shu for helpful discussions considering the experiment, and acknowledge J. Jokela and anonymous reviewers for commenting on the manuscript. Data produced and analysed in this paper were generated in collaboration with the Genetic Diversity Centre (GDC), ETH Zürich.
The study was funded by the Emil Aaltonen Foundation and the Swiss National Science Foundation (grant 31003A 169531) to OS. The funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.
Ethics approval and consent to participate
This study was carried out in accordance with the laws governing animal experimentation in Switzerland [see the animal welfare ordinance (Tierschutzverordnung): https://www.admin.ch/opc/de/classified-compilation/20080796/index.html#id-6], where work with snails does not require permissions.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Summary of the annotated transcripts linked to proteins with a role in immune defence, stress-responses, or metabolism. The total number of identified transcripts, the number of transcripts with unique open reading frames (ORFs), an example domain architecture using SMART and PFAM domains (when available in the annotation database), and transcripts representing unique ORFs for each examined factor are presented. Transcript ID the domain architecture refers to is in bold.
Principal component analysis (PCA) plots showing the variation in transcriptome-wide expression profiles of the experimental snails using the first five principal components (PCs) after internal normalisation in Sleuth. Additionally, the proportion of total variance each PC explained in the data is presented.
Expression levels of individual transcripts found to represent annotated factors related to TLR signalling pathway in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent the annotated antibacterial defence factors in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent annotated factors related to phenoloxidase/melanisation-type reaction in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using its dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent annotated factors related to non-self recognition in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent annotated cytokines in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent annotated factors related to the production of reactive oxygen species (ROS) in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent annotated factors related to apoptosis in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent the annotated stress-response factors in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent the annotated antioxidant enzymes in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
Expression levels of individual transcripts found to represent annotated factors related to metabolism in units of transcripts per million (TPM) for each experimental snail. Heatmap shows the variation for each factor using the dynamic range. Transcripts related to each factor are clustered according to their similarity.
GenBank accession numbers of the sequences of proteins/genes relevant for immune function, stress responses, antioxidation and metabolism that were used as references in BLAST similarity searches to identify orthologs in the L. stagnalis reference transcriptome.
About this article
Cite this article
Seppälä, O., Walser, JC., Cereghetti, T. et al. Transcriptome profiling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research. BMC Genomics 22, 144 (2021). https://doi.org/10.1186/s12864-021-07428-1
- Ecological immunology
- Great pond snail