Skip to main content

Transcriptome profiling of Lymnaea stagnalis (Gastropoda) for ecoimmunological research



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 [1]). 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 [9]). 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 [20]) 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 [30]. 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 [50].

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 [51]), 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 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.

Transcriptome annotation

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.

Fig. 1
figure 1

Summary of the identified immune defence factors in Lymnaea stagnalis reference transcriptome. Factors are organised across different immunological mechanisms/pathways and steps of the immune response (i.e., non-self recognition, signalling/regulation, effectors). Numbers in brackets show how many transcripts with unique open reading frames (ORFs) were detected from those factors for which determining the completeness of ORFs was possible

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.

Fig. 2
figure 2

Principal component analysis (PCA) plot showing variation in transcriptome-wide expression profiles of the experimental snails. The first two principal components (PCs) after internal normalisation in Sleuth are used. PCA plots for the first five PCs, as well as the proportion of total variance each of them explained in the data, are presented in Additional file 2

Fig. 3
figure 3

Expression levels of selected transcripts that represent different components of the immune system. Heatmap shows transcripts that deviated in their transcription between certain experimental treatments and their specific controls and components of the immune system in which individuals expressed distinct transcripts. Examined immunological pathways/mechanisms included non-self recognition, Toll-like receptor (TLR) signaling pathway, cytokines, antibacterial defence, production of reactive oxygen species (ROS), and phenoloxidase (PO)/melanisation-type reaction. Transcripts within each pathway/mechanism are clustered according to their similarity. Heatmap shows the variation for each factor among all experimental snails (each column represents one snail) using its dynamic range in units of transcripts per million (TPM). Red (injections with bacteria), blue (injections with snail/trematode tissue extracts), purple (injections with bacteria and snail/trematode tissue extracts) and black (exposures to environmental change) rectangles (dashed line), arrows and numbers refer to the specific results mentioned in the text

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].

Fig. 4
figure 4

Expression levels of selected transcripts that represent antioxidation (ANOX), stress responses and metabolism. Heatmap shows transcripts that deviated in their transcription between certain experimental treatments and their specific controls. Transcripts within each pathway/mechanism are clustered according to their similarity. Heatmap shows the variation for each factor among all experimental snails (each column represents one snail) using the dynamic range in units of transcripts per million (TPM). Purple (immune challenge) and black (exposures to environmental change) rectangles (dashed line), arrows and numbers refer to the specific results mentioned in the text


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 [59] and environmental conditions (e.g., [8, 52, 55, 60]), and they contribute differently to snail fitness [54]. 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., [64]). 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 [65]. 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 [66]. 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 [69]. 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 [70]). 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 [70]).

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 [60]. 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 [60]. 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 [52].

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 [56]. 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 [73]. 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 [74]), 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 [54]. 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 [59], 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 [84], 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 [87].

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 [54]). Interbreeding was done to increase genetic variation among the experimental snails relative to natural populations that may show low genetic diversity [90]. 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 [91]) 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).

Table 1 Summary of the experimental design

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 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 v6.6.0.0) 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 [98]. Sleuth was then used for downstream processing [99]. 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 [100]). 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



Alcohol dehydrogenase


Apoptosis-inducing factor


Antimicrobial peptide




Dual oxidase


Oxidase maturation factor


Estrogen-related receptor


Fas apoptotic inhibitory molecule 1


Fibrinogen-like domain


Fibrinogen-related protein


Gram-negative bacteria binding-protein


Glutathione S-transferase


Heat shock factor


70 kilodalton heat shock protein


90 kilodalton heat shock protein


Serine protease


Inhibitor of apoptosis proteins


Immunoglobulin superfamily


Inhibitor of NF-κB kinase


Interleukin 17


Inhibitor of nuclear factor-κB


L-amino acid oxidase


Lipopolysaccharide-binding protein/bactericidal permeability-increasing protein




PS-induced TNF-α factor


Leucine-rich repeat


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


Nuclear factor-κB


Nitric oxide synthase




Open reading frame


Poly (ADP-ribose) polymerase


Principal component


Principal component analysis


Peptidoglycan recognition protein




Protein phosphatase


Pathogen-recognition receptor


Retinoid acid receptor


Sterile alpha and armadillo motif-containing protein

SOD and MnSOD:

Superoxidase dismutases


TAK1 binding protein


Transforming growth factor β-activated kinase-1


Toll-interleukin receptor


Toll-like receptor


Tumor necrosis factor


Transcripts per million


Tumor necrosis factor receptor-associated factor


TRAF-like protein


Variable immunoglobulin and lectin-domain-containing molecules


  1. Schmid-Hempel P. Evolutionary parasitology: the integrated study of infections, immunology, ecology, and genetics. New York: Oxford University Press; 2011.

    Google Scholar 

  2. Zuk M, Stoehr AM. Immune defense and host life history. Am Nat. 2002;160:S9–S22.

    Article  PubMed  Google Scholar 

  3. Nunn CL, Lindenfors P, Pursall ER, Rolff J. On sexual dimorphism in immune function. Philos T R Soc B. 2009;364(1513):61–9.

    Article  Google Scholar 

  4. Hamilton WD, Zuk M. Heritable true fitness and bright birds: a role for parasites. Science. 1982;218(4570):384–7.

    Article  CAS  PubMed  Google Scholar 

  5. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Siva-Jothy MT, Thompson JJW. Short-term nutrient deprivation affects immune function. Physiol Entomol. 2002;27:206–12.

    Article  Google Scholar 

  7. 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.

    Article  CAS  Google Scholar 

  8. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Demas GE, Nelson RJ. Ecoimmunology. New York: Oxford University Press; 2012.

    Google Scholar 

  10. 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.

    Article  CAS  PubMed  Google Scholar 

  11. Schwarzenbach GA, Hosken DJ, Ward PI. Sex and immunity in the yellow dung fly Scathophaga stercoraria. J Evol Biol. 2005;18:455–63.

    Article  CAS  PubMed  Google Scholar 

  12. Loker ES, Adema CM, Zhang SM, Kepler TB. Invertebrate immune systems - not homogeneous, not simple, not well understood. Immunol Rev. 2004;198:10–24.

    Article  PubMed  PubMed Central  Google Scholar 

  13. Ghosh J, Lun CM, Majeske AJ, Sacchi S, Schrankel CS, Smith LC. Invertebrate immune diversity. Dev Comp Immunol. 2011;35(9):959–74.

    Article  CAS  PubMed  Google Scholar 

  14. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. 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.

  17. Hoffmann JA, Reichhart JM. Drosophila innate immunity: an evolutionary perspective. Nat Immunol. 2002;3(2):121–6.

    Article  CAS  PubMed  Google Scholar 

  18. Royet J, Reichhart JM, Hoffmann JA. Sensing and signaling during infection in Drosophila. Curr Opin Immunol. 2005;17(1):11–7.

    Article  CAS  PubMed  Google Scholar 

  19. Hetru C, Hoffmann JA. NF-kappa B in the immune response of Drosophila. Csh Perspect Biol. 2009;1(6):a000232.

  20. Pedersen AB, Babayan SA. Wild immunology. Mol Ecol. 2011;20(5):872–80.

    Article  CAS  PubMed  Google Scholar 

  21. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. 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.

    Article  CAS  PubMed  Google Scholar 

  23. 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.

  24. 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.

    Article  CAS  PubMed  Google Scholar 

  25. 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.

    Article  CAS  PubMed  Google Scholar 

  26. 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.

    Article  CAS  Google Scholar 

  27. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  28. 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.

  29. 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.

  30. 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.

  31. 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.

    Google Scholar 

  32. Prince JS. Presumptive alphavirus in the gastropod mollusc, Aplyka californica. B Mar Sci. 2003;73(3):673–7.

    Google Scholar 

  33. 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.

    Article  CAS  PubMed  Google Scholar 

  34. 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.

    Article  Google Scholar 

  35. 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.

    Article  CAS  PubMed  Google Scholar 

  36. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  37. 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.

    Article  CAS  PubMed  Google Scholar 

  38. 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.

  39. 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.

  40. 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.

  41. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  42. 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.

    Article  PubMed  Google Scholar 

  43. 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.

    Article  Google Scholar 

  44. Adema CM, Vandeutekommulder EC, Vanderknaap WPW, Sminia T. Schistosomicidal activities of Lymnaea stagnalis hemocytes: the role of oxygen radicals. Parasitology. 1994;109:479–85.

    Article  PubMed  Google Scholar 

  45. Horák P, Deme R. Lectins and saccharides in Lymnaea stagnalis haemocyte recognition. Comp Haematol Int. 1998;8(4):210–8.

    Article  Google Scholar 

  46. 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.

    Article  CAS  PubMed  Google Scholar 

  47. 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.

    Article  CAS  PubMed  Google Scholar 

  48. 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.

    Article  CAS  PubMed  Google Scholar 

  49. 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.

  50. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. 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.

  52. 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.

    PubMed  Google Scholar 

  53. Seppälä O, Jokela J. Immune defence under extreme ambient temperature. Biol Lett. 2011;7:119–22.

    Article  PubMed  Google Scholar 

  54. Langeloh L, Behrmann-Godel J, Seppälä O. Natural selection on immune defense: a field experiment. Evolution. 2017;71(2):227–37.

    Article  PubMed  Google Scholar 

  55. 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.

    Google Scholar 

  56. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  57. 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.

  58. 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.

  59. 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.

    Article  PubMed  CAS  Google Scholar 

  60. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  61. 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.

    Article  CAS  PubMed  Google Scholar 

  62. Schultz JH, Adema CM. Comparative immunogenomics of molluscs. Dev Comp Immunol. 2017;75:3–15.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. 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.

  64. Mitta G, Vandenbulcke F, Roch P. Original involvement of antimicrobial peptides in mussel innate immunity. FEBS Lett. 2000;486:185–90.

    Article  CAS  PubMed  Google Scholar 

  65. Gorbushin AM. Derivatives of the lectin complement pathway in Lophotrochozoa. Dev Comp Immunol. 2019;94:35–58.

    Article  CAS  PubMed  Google Scholar 

  66. Zhang SM, Adema CM, Kepler TB, Loker ES. Diversification of Ig superfamily genes in an invertebrate. Science. 2004;305(5681):251–4.

    Article  CAS  PubMed  Google Scholar 

  67. 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.

    Article  Google Scholar 

  68. 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.

    Article  CAS  PubMed  Google Scholar 

  69. Buchmann K. Evolution of innate immunity: clues from invertebrates via fish to mammals. Front Immunol. 2014;5:459.

  70. 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.

    Article  CAS  Google Scholar 

  71. 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.

  72. Nie L, Cai SY, Shao JZ, Chen J. Toll-like receptors, associated biological roles, and signaling networks in non-mammals. Front Immunol. 2018;9.

  73. Meehl GA, Tebaldi C. More intense, more frequent, and longer lasting heat waves in the 21st century. Science. 2004;305:994–7.

    Article  CAS  PubMed  Google Scholar 

  74. 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.

    CAS  PubMed  Google Scholar 

  75. 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.

    Article  CAS  PubMed  Google Scholar 

  76. 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.

    Article  CAS  PubMed  Google Scholar 

  77. 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.

    Article  CAS  Google Scholar 

  78. Brem RB, Yvert G, Clinton R, Kruglyak L. Genetic dissection of transcriptional regulation in budding yeast. Science. 2002;296(5568):752–5.

    Article  CAS  PubMed  Google Scholar 

  79. 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.

    Article  CAS  PubMed  Google Scholar 

  80. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  81. Oleksiak MF, Roach JL, Crawford DL. Natural variation in cardiac metabolism and gene expression in Fundulus heteroclitus. Nat Genet. 2005;37(1):67–72.

    Article  CAS  PubMed  Google Scholar 

  82. Whitehead A, Crawford DL. Neutral and adaptive variation in gene expression. P Natl Acad Sci USA. 2006;103(14):5425–30.

    Article  CAS  Google Scholar 

  83. 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.

    Article  CAS  PubMed  Google Scholar 

  84. 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.

    Article  CAS  PubMed  Google Scholar 

  85. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  86. 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.

    Article  Google Scholar 

  87. Théron A, Coustau C. Are Biomphalaria snails resistant to Schistosoma mansoni? J Helminthol. 2005;79(3):187–91.

    Article  PubMed  Google Scholar 

  88. 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.

    Article  Google Scholar 

  89. 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.

    PubMed  PubMed Central  Google Scholar 

  90. 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.

    Article  Google Scholar 

  91. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  92. 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.

    Article  CAS  PubMed  Google Scholar 

  93. 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.

    Article  CAS  Google Scholar 

  94. 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.

    Article  CAS  PubMed  Google Scholar 

  95. 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.

    Article  PubMed  CAS  Google Scholar 

  96. 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.

    Article  CAS  PubMed  Google Scholar 

  97. Gorbushin AM, Borisova EA. Lectin-like molecules in transcriptome of Littorina littorea hemocytes. Dev Comp Immunol. 2015;48(1):210–20.

    Article  CAS  PubMed  Google Scholar 

  98. Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nat Biotechnol. 2016;34(5):525–7.

    Article  CAS  PubMed  Google Scholar 

  99. Pimentel H, Bray NL, Puente S, Melsted P, Pachter L. Differential analysis of RNA-seq incorporating quantification uncertainty. Nat Methods. 2017;14(7):687.

    Article  CAS  PubMed  Google Scholar 

  100. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

Download references


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.

Author information

Authors and Affiliations



OS, JCW and KS conceptualised and designed the study. OS and KS conducted the experiment. JCW assembled the tanscriptome, and OS, TC, TS and CMA annotated it. OS and CMA examined the variation in transcription and wrote the first draft of the manuscript. JCW, TC, KS and TS provided comments on the manuscript. All authors have read and approved the manuscript.

Corresponding author

Correspondence to Otto Seppälä.

Ethics declarations

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):], where work with snails does not require permissions.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1.

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.

Additional file 2.

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.

Additional file 3.

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.

Additional file 4.

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.

Additional file 5.

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.

Additional file 6.

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.

Additional file 7.

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.

Additional file 8.

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.

Additional file 9.

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.

Additional file 10.

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.

Additional file 11.

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.

Additional file 12.

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.

Additional file 13

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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: