Freshwater snails play a key role in the transmission of schistosomiasis, a tropical parasitic disease affecting over 150 million people. Adaptation of these snails to local climatic conditions is a critical factor in determining how climate change and other environmental factors influence disease transmission dynamics, yet this potential adaptation has remained unexplored. Bulinus truncatus is the schistosome intermediate host snail with the widest geographic distribution and is therefore an important factor determining the maximum range of urogenital schistosomiasis. In this study, we assessed the local adaptation capacity of B. truncatus to temperature through an integrative approach encompassing phenotypic, ecophysiological, and genomic data. Ten snail populations from diverse thermal environments were collected in three countries, with eight populations reared in a common garden. The F2 generation (N= 2,304) was exposed to eight chronic temperature treatments (± 36 snails/population/temperature treatment) and various life-history traits were recorded for over 14 weeks. Subsequently, ecophysiological analyses were conducted on the ten last surviving snails per population. Genotyping the parental generation collected in the field using a genotyping-by-sequencing (GBS) approach, revealed 12,875 single nucleotide polymorphisms (SNPs), of which 4.91 % were potentially under selection. We observed a significant association between outlier SNPs, temperature, and precipitation. Thermal adaptations in life-history traits were evident, with lower survival rates at high temperatures of warm-origin snails compensated for by higher reproduction rates. Cold-origin snails, on the other hand, exhibited higher growth rates adapted to a shorter growing season. Ecophysiological adaptations included elevated sugar and haemoglobin contents in cold-adapted snails. In contrast, warm-adapted snails displayed increased protein levels but also more oxidative damage. Furthermore, heightened phenoloxidase levels indicated a more robust immune response in snails from parasite-rich regions. These morphological and physiological differences provide convincing evidence for a genetic basis of local adaptation. This in turn holds profound implications for the snail’s response to climate change, future schistosomiasis risk, and the effectiveness of schistosomiasis control measures.

All data were obtained through an extensive common garden experiment in which second generation snails from different origins were subjected to different temperature treatments. Life-history parameters such as growth, fecundity and survival were measured during 14 weeks. After the experiment, a subset of snails were analysed physiologically. The DNA from the parental generation snails coming from the field was extracted and single nucleotide polymorphisms (SNPs) were identified using a genotyping by sequencing approach.

Field Collection and experimental design

Bulinus truncatus snails were collected in the second half of 2021 from ten localities in three countries (France, Senegal and Zimbabwe) spanning a latitudinal and temperature gradient. About 100 snails were collected per locality. Additionally, a French (Corsican) strain that has been maintained in the lab since 2014 was included. Each locality was screened for B. truncatus and snails were collected from every habitat where they were found. We chose French sampling sites on the border between the Corsican montane broadleaf and mixed forests and the Tyrrhenian-Adriatic Sclerophyllous and mixed forests ecoregions in the hot summer Mediterranean climate as this area represents the cold limit of the distribution range of B. truncatus. This area is characterised by a high seasonality with average air temperatures ranging between 7 °C in winter and 27 °C in summer. Senegalese snails were collected in the Sahelian Acacia savanna ecoregion characterised by a hot desert climate which represents the warm limit of the distribution range with less strong seasonality and average air temperatures ranging from 30 to 39 °C. Finally, the snails from Zimbabwe originate from the Zambezian and Mopane woodlands (Triangle and Malilangwe populations, characterised by a hot, semi-arid steppe climate) and the Southern Miombo woodlands (Imire, characterised by a subtropical highland climate) ecoregions where temperatures show less seasonality and are more temperate with averages between 14 and 22.6 °C. The field-collected parental (P) snails were bred to the second (F2) generation in a common garden at 24 °C and a 12:12 h day-night regime, fed ad libitum with pesticide-free lettuce (dried at 55 °C for 12 h), and green macro-algae were added to enrich the water with oxygen. About 100 random egg masses from the P generation were transferred to start the F1 generation. As many eggs as possible were collected from the F1 generation to breed at least 300 F2 snails per locality, needed for the experiments.

DNA was extracted from the whole body (from the P generation collected in the field) using the E.Z.N.A.® mollusc DNA kit according to the manufacturer’s protocol (Omega Bio-tek). To verify the morphological identification from the field, a subsample of five snails from each population was genotyped through DNA barcoding using the protocol described in Maes et al. (2022) and identification of the specimens relied on a BLAST search against the NCBI database GenBank. After DNA barcoding, eight populations from three countries identified as B. truncatus (see results section) were selected for the life history experiment while all populations identified as B. truncatus were included in the genetic analyses. The Imire population from Zimbabwe was identified as Bulinus tropicus (Krauss, 1848) and excluded from the experiment and genetic analyses.

The F2 snails were individually placed in 100 ml plastic cups filled with aged tap water at four to eight weeks old. The experiment was started in 12 different batches, each batch containing three snails per locality for each temperature treatment (3 snails x 8 localities x 8 temperatures x 12 batches = 2304 snails in total, 36 snails/population/temperature). The temperature was increased/decreased by 2 °C/day starting from 24 °C until the snails reached their experimental temperature (4, 8, 12, 18, 24, 28, 32 or 36 °C) at the start of the experiment (time t=0). A constant temperature was permanently monitored using Hobo onset data loggers (tidbit v2 Temp logger). The snails were fed ad libitum with dried lettuce and green algae, and the water was refreshed weekly or when the oxygen level dropped below 5 mg/ml. The lettuce was changed daily in the high-temperature treatments (24, 28, 32 and 36 °C) to prevent its decomposition and associated oxygen depletion. The snails were kept at the experimental temperatures for 14 weeks.

Response variables

Each week the number of egg masses per snail was counted, dead snails were counted and removed, and shell height (apex to bottom of aperture) was measured under a stereomicroscope (Ceti Steddy-B) with a built-in ruler. The ten last surviving snails per treatment and locality combination were taken out of their shell, the body was weighed with an accuracy of 0.01 mg (Mettler Toledo AB135-S, Columbus, Ohio, USA) and stored at -80 °C for ecophysiological analyses. Egg masses were collected in weeks one to seven and in week 12 and stored in 98 % ethanol to later quantify the number of eggs per clutch under a stereomicroscope.

The critical thermal maximum and minimum (CT_max and CT_min, respectively; Bartnicki et al., 2021; Morgan et al., 2018) were measured on a separate subset of F2 snails that were not included in the chronic temperature experiment. For measuring the CT_max, the snails were individually put in 12 ml plastic tubes filled with aged tap water and placed in a heating block (Thermo Fisher Scientific). The water temperature increased by 0.1 °C/min (following Johansson & Laurila, 2017). For measuring the CT_min, the snails were also individually put in 12 ml plastic tubes filled with aged tap water and placed in a water bath. The water temperature was decreased by 0.1 °C/min using a cooler. Before both runs, the size of each snail was measured under a stereomicroscope with a built-in ruler. Snails were considered fainted when they showed no response to tactile stimuli (i.e. the snails did neither retreat in their shell nor retreat their tentacles). After the test, snails were allowed to recover and only snails that fully recovered within 10 min were considered for the analysis (5 out of 120 snails died in the CT_max experiment, none in the CT_min experiment).

For the preferred temperature experiment, ten aluminium U-profiles (lanes) of 2 m in length were attached to create ten replicas and filled with aged tap water to establish a temperature gradient to assess the preferred temperature (T_p) (cf. Johansson & Laurila, 2017). One side of the profiles was placed on a cooling plate, and the other side on a heating plate. The temperatures of both plates were adjusted to generate a stable temperature gradient of 24 °C in the middle of each lane and an increase/decrease of 1 °C/10 cm in the direction of the heater/cooler (min: 14 °C max: 34 °C). Ten snails that were not included in the chronic temperature experiment were randomly selected from each locality and randomly assigned to one of the ten lanes. The snails were placed in the middle at 24 °C and the temperature at the position of each snail was measured every 15 min for 4 h. The first 2 h were considered as acclimation time while the preferred temperature was defined as the average temperature during the last two hours. We assumed that there were no differences in movement speed between the populations (Dillon et al., 2012).

Ecophysiological analyses

The last surviving individuals per population and temperature treatment were weighed and stored at -80 °C until further analyses. Whole-body snail samples were homogenised, 15 X diluted in phosphate-buffered saline (PBS) and centrifuged at 4 °C for 10 min at 13,000 g. A minimum of 8 mg snail wet weight was required to run all ecophysiological analyses. If this weight was not reached, we pooled two snails from the same locality and temperature treatment to obtain the minimum required weight.

Various ecophysiological parameters were examined to gain a comprehensive understanding of the snail’s overall well-being and to identify nuanced signs of local adaptation that might not be evident through life history traits alone. Assessing the snails’ general condition involved determining fat, total sugar, and protein contents (for the full protocols, see Appendix S1). From these variables, the cellular energy allocation (CEA) was calculated as outlined in Gomes et al. (2015). The CEA integrates the energy available (Ea) and energy consumption (Ec) of an organism. Additionally, the haemoglobin content was measured to gauge the oxygen-binding capacity. Haemoglobin is absent in most freshwater snails, except for pulmonate species like B. truncatus (Lieb et al., 2006). Furthermore, the malondialdehyde (MDA) levels were quantified to assess oxidative damage to lipids (Miyamoto et al., 2012). Phenoloxidase (PO), instrumental in the immune response to various parasite species in snails (Le Clec’h et al., 2016), was also measured to evaluate the strength of the snails’ immune system. Lastly, the snails’ metabolic rate was quantified at the cellular level by assessing the activity of the electron transport system (ETS) (De Coen & Janssen, 2003).

SNP genotyping

DNA from the parental (P) generation snails was used to build a paired-end GBS reduced representation library (RRL) according to Elshire et al. (2011) and sequenced on an Illumina NovaSeq 6000 platform (Genomics Core, KU Leuven). High molecular weight DNA was digested with restriction enzyme NsiI; Illumina sequencing primers P1 and P2 and adapters containing the barcodes were ligated to the resulting fragments.

The NovaSeq run produced 2.01 x 10⁹ paired-end raw 101 and 92 bp reads (forward and reverse reads, resp.). Reads were demultiplexed using process_radtags from Stacks 2.5 (Catchen et al., 2013) and the sequence quality was checked using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc). FastX trimmer (http://hannonlab.cshl.edu/fastx_toolkit/) was used to trim the first 7 bp of the forward reads and the first 6 bp of the reverse reads. Trimmed forward and reverse reads were mapped against the B. truncatus reference genome (Young et al., 2022) using Bowtie2 (Langmead et al., 2019) and SNPs were identified using GATK UnifiedGenotyper 3.7. The resulting vcf file containing all variants was further filtered to only include biallelic SNPs, a read depth between 20 and 100, a minimum genotyping quality of 20, a minimum allelic depth of 6, a maximum of 40 % missing data per locus, a maximum of 70 % missing data per individual and a minor allele frequency of 0.02. Since B. truncatus is presumed to be autotetraploid (Young et al., 2022), a ploidy test was carried out using nQuire (Weib et al., 2018). NQuire uses next-generation sequencing data to distinguish between different ploidy levels based on the frequency distributions at variant sites where only two bases are segregating. SNP identification was carried out with both diploid (n=2) and tetraploid (n=4) settings in UnifiedGenotyper and a PCA was carried out on both filtered vcf files to check if the two datasets gave similar outcomes. Since both datasets gave the same output, all downstream analyses were carried out using the diploid dataset since more analysis tools are available for diploid data.