Gastropods were sampled twice a year from autumn 2021 to spring 2024 in anthropized green spaces invaded by O. nungara, across four regions of mainland France (Paris region, La Rochelle, Bordeaux region and Brittany). Sampling was carried out over a standardized 2-hour survey in all sheltering microhabitats within a 100 m² area at each site. Samples were stored in absolute ethanol at 4 °C prior to DNA extraction. In a preliminary identification step, samples were assigned to morphospecies using identification keys from Gargominy (2020) and from Mc Donnell et al. (2009). A total of 469 individuals were collected (ranging from 1 to 61 per morphospecies, depending on availability in the field and rarity of the species), and DNA barcoding was then used to assign specimens to species or genetic lineages. Total DNA was extracted from about 4 mm² of the mantle with the DNeasy Blood and Tissue kit (Qiagen, France) following the manufacturer instructions. The standard barcode fragment of 650 pb of the cytochrome oxidase I gene (COI) was amplified and sequenced using LCO1490/HCO2198 primer pair from Folmer et al. (1994). In total COI PCRs were performed on 400 individuals, yielding 349 successful COI barcodes. Using these barcodes, 38 species or genetic lineages of terrestrial snails and slugs were identified using BOLD Identification System (IDS) and BLAST from NCBI (BOLD dataset: http://dx.doi.org/10.5883/DS-MOLPWP1).

To develop and validate specific primers for the amplification and sequencing of DNA from terrestrial gastropods in predator gut contents, available primers and sequences for several genes were collected from the literature and the National Center for Biotechnology Information (NCBI) database. To select the most appropriate gene for this study, several criteria were applied: the gene had to be well documented, with sufficient sequences available to build a reference database, and suitable for targeting taxa within terrestrial gastropods, which encompass more than ten families. Although COI is commonly used, preliminary in silico and PCR tests did not reveal regions suitable for designing exclusive primers for gastropods. Moreover, previous studies have shown that primer bias in COI can lead to the underrepresentation or non-detection of some taxa, even within closely related groups (Deagle et al. 2014; Elbrecht and Leese 2015; Piñol et al. 2015). The mitochondrial 16S rRNA gene was therefore selected, as it combines conserved regions for primer design with sufficient variability for species-level discrimination. This choice is supported by several studies that have highlighted 16S as a suitable alternative marker for metabarcoding in other taxa (Clarke et al. 2014; Deagle et al. 2014; Elbrecht et al. 2016).

In order to build the mitochondrial 16S rDNA database, the primer pair 16sar-L (5’-CGCCTGTTTATCAAAAACAT-3’) (Palumbi 1996) and modified 16sbr-H (5’-CCGGTYTGAACTCAGATCAYGT-3’) (Palumbi 1996, modified in Bossuyt and Milinkovitch (2000)) was used to amplify a 550 bp 16S fragment (hereafter referred to as “large fragment”) from one to three samples per species/lineage when available, giving priority to maximizing genetic diversity (by covering distinct COI clades when applicable) and geographic spread when possible, in order to avoid redundancy and to better represent intraspecific diversity for a total of 110 individuals. Reactions were carried out in a volume of 25 µL containing 2.5 µL of DNA extract. Amplification mix contained 0.08 U of GoTaq® G2 DNA Polymerase (Promega, France), 1X Colorless GoTaq® Reaction Buffer, 2 mM of MgCl2, 0.5 mM of dNTPS, 0.2 µM of each primer and 0.25 mg/mL of bovine serum albumin. The following cycling conditions were used: an initial denaturing step at 94 °C for 3 min; 36 cycles of denaturation at 94 °C for 30 s; annealing at 51 °C for 30 s, extension at 72 °C for 1 min and final extension step at 72 °C for 10 min. PCR products were sent to Eurofins Genomics for sequencing using the 16sar-L primer. Sequences were manually corrected using Chromas 2.6.6 software (https://technelysium.com.au) and submitted to BOLD and GenBank (all references and data available in BOLD dataset: http://dx.doi.org/10.5883/DS-MOLPWP1).

To improve species coverage for the primer design, 244 sequences were retrieved from NCBI and added to the database. These sequences were selected to represent both species recorded at our study sites and those likely to occur in French habitats potentially colonized by O. nungara. This addition represents 20 new species and lineages, bringing the total to 58, which corresponds to 13% of the French terrestrial gastropod fauna, considering species that could potentially be encountered based on their geographical distribution. Sequences were then aligned using the MUSCLE algorithm in SeaView4 (Gouy et al. 2010).

215 earthworm 16S sequences from Barraux et al. (2024) and from our laboratory (not published), and the only available 16S sequence closely related to O. nungara from NCBI (AN: KP208777, listed as Obama sp.) were added to the gastropod alignment. Highly conserved regions surrounding variable regions were visually identified, and specific primer pairs (i.e. expected to amplify only gastropod DNA and not earthworm or O. nungara DNA) targeting ~100 bp variable sequences were manually designed. Standard primer design checks (e.g. annealing temperature balance, absence of strong secondary structures such as hairpins or primer–dimers, assessed with the Multiple Primer Analyzer web tool from Thermo Fisher Scientific) were performed before empirical validation.

To assess the overall coverage and specificity of the designed primer pairs for terrestrial gastropods, in silico PCR were first performed using the ecoPCR program (Ficetola et al. 2010). Each primer pair was tested against 18,210 mitochondrial reference sequences for Eukaryota (18,209 plus Obama sp.) retrieved from the NCBI Reference Sequence Database (RefSeq) (https://ftp.ncbi.nlm.nih.gov/genomes/refseq/mitochondrion/), and against 26,578 bacterial 16S rRNA sequences from the Bacterial 16S Ribosomal RNA RefSeq Targeted Loci Project (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA33175). Up to two mismatches per primer were allowed, with a single mismatch within the last three bases at the 3’ end considered a PCR failure.

PCR specificity cross-tests were also conducted on DNA extracts (see below CTAB extraction procedure) from three gastropod species, three earthworm species and eleven samples of O. nungara that were fed differently (with an empty stomach, previously fed with terrestrial gastropods or earthworms, or with gastropods and earthworms). The PCR reaction was carried out using the same reagent composition and volumes as previously described. The following cycling conditions were used: an initial denaturing step at 94 °C for 3 min; 40 cycles of denaturation at 94 °C for 30 s; annealing at variable temperature depending on primer pairs for 30 s, extension at 72 °C for 1 min and final extension step at 72 °C for 10 min.

To ensure broad coverage of terrestrial gastropod species, promising primers were tested in silico using the PRIMERMINER 0.21 R package (Elbrecht and Leese 2017). The penalty score, which quantifies mismatches between primer and target sequences—including mismatch type, position, and adjacency—was calculated for each primer. Sequences with a score below 120 were considered likely to be successfully amplified, whereas scores above 120 indicated low amplification probability, with 120 serving as the practical threshold for primer selection. PCR amplifications on 49 DNA samples from the 38 gastropod COI species/lineages were also conducted with the protocol described above.

Finally, a Neighbor-Joining tree was reconstructed with the software MEGA (Kimura 2-parameter model, 1000 Bootstrap replicates) based on the short DNA sequences (i.e. 16S metabarcoding region) targeted by the retained primers, to verify whether the different species or lineages occur in distinct clades.

To discuss the potential of using the primer pair in future studies studying terrestrial gastropods, its coverage in silico was also tested on available 16S sequences of terrestrial gastropod species from France (sequences available for 206 species out of 422 described, across 43 families) (list extracted from https://uicn.fr/liste-rouge-mollusques-continentaux/) and from Europe (sequences available for 698 species out of 2480 described, across 53 families and including sequences from this study) (list extracted from https://portals.iucn.org/library/node/48439) using the PRIMERMINER package.

Two types of mock community controls were set up (See Suppl. material 1). Three mock community controls (“POS-G1 to G3”) were each composed of DNA from five different gastropod species randomly introduced at the same concentration (21 ng/µL). Three other mock community controls (“NREP-G1 to G3”) were each composed of 27 species among the most abundant in anthropized sites. DNA volumes were diluted to the same concentration as above and randomly graduated from 10 µL to 1 µL depending on species (See Suppl. material 1). Both types of mock communities were designed to ensure that all introduced DNA from species were retrieved in proportional sequence counts, to assess the contamination level and to set a filtering threshold for clusters (i.e., sets of similar sequences grouped based on a predefined similarity threshold) (De Barba et al. 2014; Ficetola et al. 2015; Port et al. 2016; Thomas et al. 2016; Andruszkiewicz et al. 2017; Braukmann et al. 2019).

Six “controlled-feeding” controls were DNA from digestive contents of O. nungara individuals fed for three weeks with one prey per week. Three individuals (“UNI-G1 to G3”) received the same species, Cornu aspersum (non-differentiated genetic lineages, as they could not be recognized based on the external morphology). The other three individuals (“MIX-G1 to G3”) received a different species every week in the following order: Hygromia cinctella, Cornu aspersum and Arion hortensis L1. Consumption was checked after three days, and prey was always eaten. After three weeks, O. nungara individuals were killed and stored in absolute ethanol for further DNA extraction. Both types of “controlled-feeding” controls were used to test if all given prey would be detected by sequencing digestive contents after three weeks of experimentation and if there were potential amplification biases between DNA from different species. They were also used to assess contamination levels and to set a filtering threshold for clusters (De Barba et al. 2014; Port et al. 2016; Andruszkiewicz et al. 2017).

Three negative controls (“BLANK 1 to 3”), from three independent PCR amplifications without added DNA, were set up to check for possible contamination during PCR and to set the threshold of stochasticity for technical replicates. A pairwise distance matrix between negative replicates was created using the Bray-Curtis dissimilarity index (Bray and Curtis 1957) based on sequence frequency data. The distance value ranged from 0, indicating identical relative abundance profiles, to 1, representing completely distinct assemblages (no shared sequences). Since negative controls are not expected to have repeatable sequence frequencies, this allowed a threshold of stochasticity to be set for technical replicates, with samples exceeding a maximum value being discarded (De Barba et al. 2014; Corse et al. 2017).

The last three controls were DNA from digestive contents of three O. nungara individuals collected on the field during sampling campaigns (V6, V29 and V53 samples) and tested positive for gastropod DNA. Each of them was replicated three times (i.e. three independent PCR amplifications and sequencing per DNA to create technical replicates). Replicated field controls were used to apply parameters and thresholds (DNA quantity, distance threshold between technical replicates, threshold for clusters filtering) calibrated with other controls mentioned above, to ensure their reliability on real digestive contents from the field (De Barba et al. 2014; Pansu et al. 2015; Piñol et al. 2015; Beentjes et al. 2019; Mata et al. 2019).

Finally, 16 DNA extracts from O. nungara individuals (5 from Limeil-Brevannes, Paris region, 3 from Saint-Medard-en-Jalles, Bordeaux region, 3 from Cleguer, Brittany and 5 from Concarneau, Brittany) that tested positive for gastropod DNA were included in the experiment. By estimating the diversity of terrestrial gastropod species at these particularly rich collection sites, it was possible to experiment in vivo the retained primer pair and the filtering threshold for clusters calibrated from metabarcoding controls.

DNA from about 20 mg of gut content was extracted from the region corresponding to the opening of the mouth of O. nungara according to Roy et al. (2022). DNA extraction was performed with CTAB method adapted from Sambrook and Russell (2001). Samples were incubated overnight at 57 °C and DNA was extracted twice with an equal volume of chloroform: isoamyl alcohol (24:1).

DNA amplifications were then conducted on all 40 samples by using the validated primer pair and the corresponding PCR protocol (see Results). For each sample and replicate, a separate library was prepared with a target sequencing depth of 10,000 reads. Each library was assigned a unique index, following GenoScreen’s protocol for multiplexing. Sequencing was performed on the Illumina MiSeq platform (2 × 150 bp paired-end chemistry) with 15% PhiX.

The Migale platform (https://galaxy.migale.inrae.fr/) was used to trim, merge, filter and cluster sequences. Pair-end reads were first trimmed with FROGS SICKLE program, based on a Phred quality score Q < 30. Merging of pair-end reads followed by first filtering were realized with FROGS PRE-PROCESS program with the following parameters: maximum lengths of 127 bp for forward reads (R1) and 134 bp for reverse reads (R2) were selected for merging, based on reads length distribution reports from FASTQC and MULTIQC programs; >30 bp coverage zone and 97% nucleotide identity were required to assemble paired-end reads and merged amplicons not comprised between 119 and 135 bp were discarded.

Clustering was then performed using FROGS CLUSTERING SWARM program with d = 1 (a small local linking threshold, representing the maximum number of differences between sequences in each aggregation swarm step) and the fastidious option selected. Clustered sequences will be referred to and identified hereafter as Operational Taxonomic Units (OTUs) in this study. Sequences identified as chimera in all samples were discarded using FROGS REMOVE CHIMERA program. Finally, clusters were filtered using the FROGS CLUSTER FILTERS program: only singletons were discarded, and the OTU filtering threshold established from control analyses was applied to both the technical replicates and the sixteen samples from the field.

For the taxonomic assignment of OTUs, a double BLAST approach was used: first with the BLASTn program from NCBI and second, with the reference database. Only BLAST matches with at least 95% identity and 95% sequence coverage were used to assign a gastropod species or lineage to each OTU. Unassigned OTUs, including those under the threshold and those unknown from any reference database, were then discarded from the analysis to maintain clean assignments and real detections. To further validate the taxonomic assignments, a Neighbor-Joining phylogenetic tree was reconstructed including OTU sequences with the corresponding reference sequences. This tree was then inspected to confirm that OTUs assigned to the same nominal species were clustered together and that the reference sequences of each species were included within the same clade.

Sequencing results were finally analyzed in R 4.4.2 (R Core Team 2024) with the PHYLOSEQ 1.46.0 R package (McMurdie and Holmes 2013) for each type of metabarcoding controls and field data, including data transformation for qualitative interpretation and graphical visualizations. Following the methodology of De Barba et al. (2014), an approach based on the relative read abundance (RRA) (i.e. proportional summaries of count reads for each detected prey) was chosen at the individual scale to drive the analysis (Casey et al. 2019; Deagle et al. 2019; Ando et al. 2020).

Ninety-six 16S sequences were obtained from DNA extracts of gastropods, representing 38 species and lineages described by morphology and/or the COI gene from French sites invaded by O. nungara. Including sequences from NCBI, the final database contained 340 sequences from 58 species/lineages (1 to 40 sequences per species or lineage) (See Suppl. material 2) and the alignment spanned 489 bp, including gaps.

Six primers, ranging from 14 to 20 bp, were designed by visualizing the 16S alignment of sequences from terrestrial gastropods, earthworm and O. nungara (See Suppl. material 3). These primers were complemented by three primers specific to gastropods obtained from Taberlet et al. (2018), resulting in eight primer pairs targeting short DNA fragments (98–217 bp) (Fig. 1).

Figure 1.

Schema of 16S DNA sequence in terrestrial gastropods showing variable (white) regions, conserved (green) regions and positions of the primers designed and retrieved for this study, the lengths indicated include the primers; and table detailing names, sequence lengths, 5’-3’ sequences and references for each potential primer pair.

In silico PCR results indicated that none of the primer pairs amplified bacterial 16S sequences, with amplification restricted to eukaryotic taxa. Among them, three primer pairs (the two including the MVF primer and Gast01R-Gast01F) produced the highest number of amplicons. Gast01R-Gast01F amplified exclusively Animalia, while the MVF-containing pairs retrieved mostly Animalia but also some sequences from Fungi and a small number of sequences from Chromista, Plantae, Protozoa, and an unassigned Eukaryota taxon (NC_086778). The five remaining primer pairs amplified fewer sequences overall, all belonging to Animalia (Fig. 2A). Within Mollusca, the performance of primer pairs has been examined on Stylommatophora, an order that includes most of terrestrial gastropods. The reference dataset contained 45 sequences from 45 Stylommatophora species. In silico PCR showed that Gast01R-Gast01F could theoretically amplify all 45 species, while the other primer pairs amplified between 26 and 40 species, except for GastD-Gast01R-F, which retrieved only about 10 species. Apart from primer pairs involving MVF, all amplified exclusively terrestrial gastropods, with no co-amplification of Obama sp. or Lumbricidae (See Suppl. material 4). Based on these results, Gast01R-Gast01F, VRF-VRR, and VRF-MVR appeared to be the most promising primer pairs for broad coverage and specificity of terrestrial gastropods (Fig. 2B). Additional details on amplification across kingdoms, phyla, and other taxonomic levels are provided in Suppl. material 4.

Figure 2.

In silico evaluation of primer pairs across eukaryotic taxa. A Global amplification per primer pair. Each bar represents the number of species amplified per primer pair, colored by kingdom. B Amplification specifically within terrestrial gastropods (order Stylommatophora). Bars indicate the number of Stylommatophora species amplified per primer pair. Dashed vertical line indicates the maximum number of target species available (n = 45).

Specificity cross-tests showed that, contrary to in silico predictions, two primer pairs (VRF-VRR and MVF-VRR) amplified only the targeted gastropod taxa and O. nungara samples fed with gastropods. Two primer pairs (VRF-MVR and GastD-GastE) amplified some of the targeted gastropods but also produced unspecific amplification. Four other primer pairs showed broader cross-amplification: MVF-MVR co-amplified targeted gastropods together with O. nungara from empty or earthworm-fed stomachs and earthworm DNA; GastD-Gast01F co-amplified targeted gastropods and O. nungara from empty and earthworm-fed stomachs; and GastD-Gast01R-F amplified the targeted gastropods but also co-amplified O. nungara from empty stomachs and produced non-specific amplicons (Fig. 3). The primer pair Gast01R-Gast01F, although initially promising in terms of taxonomic coverage in silico (Fig. 2) and in the literature (Taberlet et al. 2018), showed limited resolution for some key genera when tested in our phylogenetic alignments (poor taxonomic resolution for Deroceras species, Theba and Cornu lineages). In contrast, VRF-VRR consistently yielded reproducible amplification across different DNA samples. Given its better performance and broader utility, VRF-VRR was selected for all subsequent analyses, and Gast01R-Gast01F was not pursued further.

Figure 3.

PCR specificity cross-tests results for the three taxonomic groups (O. nungara, earthworms, and gastropods) using different potential primer pairs for minibarcode design. Bars indicate the number of species and samples showing each amplification outcome. “Targeted amplification” includes the three gastropod species (A. parvipenis, C. aspersum L1, D. invadens), the three O. nungara fed with gastropod prey (g. prey) and the two fed with both gastropod and earthworm prey (both prey); “Co-amplification” includes non-targeted amplifications either with the three O. nungara with empty stomach (e.s.), the three O. nungara fed with earthworm prey (e. prey) or with earthworm species (A. icterica, L. terrestris, L. castaneus L1). Colors correspond to amplification categories as indicated in the legend, annealing temperatures are indicated for each primer pair.

In silico amplifications using PRIMERMINER 0.21 showed that DNA sequences from 30 out of the 38 lineages from the field could theoretically be amplified with the VRF-VRR primer pair (Table 1) based on the standard penalty score of 120. Twenty-six of these lineages had a penalty score of zero and were therefore perfectly matched to the primer pair, while four species had a penalty score between 36 and 51. Two lineages showed haplotype-dependent results (Arion hortensis L2: 0 to 152, and Theba pisana L1: 71 to 325), and PCR amplification was predicted to fail for each sequence of six species (Table 1). Except for the family Pomatiidae, represented only by Pomatias elegans, primers could amplify at least one species for each of the other eleven families (all species from this study for most of them). Otherwise, the results of the in silico amplifications extended to French and European species showed that the VRF-VRR primer pair could theoretically be able to amplify DNA from at least 178 species found in France across 37 families, and 607 species found in Europe across 47 families. All detailed results for France and Europe are available (See Suppl. materials 5, 6).

The results of PCR assays on gastropod DNA showed that 35 out of the 38 species and lineages tested showed the expected amplification band. Pomatias elegans was the only species that did not amplify, and two other lineages from a morphospecies (Theba pisana L1 and L2) showed unspecific results (Table 1).

Fifty-six out of 58 species and lineages were monophyletic in the NJ tree reconstructed based on the short fragment: they grouped into distinct clades or were distinct singletons (See Suppl. material 7). This is particularly relevant for the very close species of the genus Deroceras (Deroceras invadens and Deroceras panormitanum, and Deroceras agreste and Deroceras reticulatum) (Ventura et al. 2025) and the lineages of Arion hortensis (2 lineages never described) and Arion subfuscus (5 lineages) described in the literature (Pinceel et al. 2004, 2005) that are differentiated with the short 16S metabarcoding region. Furthermore, recently differentiated subspecies of Arion ater (Reise et al. 2020), Arion ater ater, Arion ater ruber and Arion ater rufus could also be distinguished. Finally, two species, Ambigolimax parvipenis and Ambigolimax valentianus, were indistinguishable, as they merged in the tree reconstructed based on the short 16S metabarcoding region.

High-throughput sequencing resulted in 714,602 paired-end reads for the 40 samples (metabarcoding controls and field samples). Details of reads and OTU counts through the bioinformatic process before analysis of the clustering filter threshold are shown in Table 2, and raw reads have been deposited in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB100783 (https://www.ebi.ac.uk/ena/browser/view/PRJEB100783). As expected, the three negative controls (“BLANK 1 to 3”) retained few sequences (5 to 8 per sample) after trimming and merging procedures. In these controls, the vast majority of reads were discarded because they corresponded to very short sequences (less than 25 bp). The mock communities (“POS-G1 to G3” and “NREP-G1 to G3”) and the three replicate samples (V6, V29 and V53, from digestive contents of O. nungara collected in the field) showed an average of 5,000 analyzable sequences. The “controlled-feeding” controls showed a considerable decrease after the trimming and merging of paired-end reads, with an average of 911 ± 498 sequences for “UNI-G1 to G3” and 1870 ± 540 for “MIX-G1 to G3” before cluster filtering. Finally, the sixteen field samples yielded an average of 17,000 analyzable sequences (Table 2).

In all controls, OTUs were either assigned to terrestrial gastropod species (based on more than 95% identity and 95% coverage with NCBI sequences and confirmed by alignment of OTU sequences with the reference database), or showed “no significant similarity” in both BLAST analysis and the database (i.e. likely artefactual sequences). The number of OTUs per nominal species or lineage varied from 1 to 6 OTUs across all samples. Some species showed high OTU richness, such as Cochlicella barbara (6 OTUs) and Cornu aspersum L1 and L2 (4 OTUs each), but phylogenetic reconstruction confirmed the clustering of OTUs assigned to the same nominal species with reference sequences (data not shown). As expected, the 16S region amplified by the VRF-VRR primer pair only failed to distinguish Ambigolimax parvipenis from Ambigolimax valentianus.

In POS-G1 and POS-G3 mocks, all five expected gastropod species were recovered, with more than 99% of sequences from these species on average, followed by unassigned clusters and artefactual sequences (Fig. 4A-C). DNA sequences from four of the five expected species from POS-G2 were also detected, representing 87% of the sequences. However, DNA sequences from the fifth species, Deroceras invadens, were not detected in POS-G2 or in any of the NREP mock communities, and instead sequences from Deroceras panormitanum were detected, suggesting a laboratory error (Fig. 4).

Figure 4.

Comparison of expected and observed sequencing results for both types of mock community: A–C (POS) and D–F (NREP). Results are displayed as relative read abundance (RRA, % of sequences) after applying a 0.1% threshold filter for unexpected and non-assigned OTUs (representing 99.9% of the target gastropod taxa) and clustering reads from the same expected species. For D–F (NREP), each point represents the observed abundance of a given taxon in one sample; the red dotted line indicates the theoretical 1:1 relationship between observed and expected abundances, and the gray dashed line shows the regression across all points, representing the observed trend. The color legend applies for both type of mock community.

The VRF-VRR primer pair amplified and discriminated 16S DNA sequences from 24 out of the 27 species and lineages included in the “NREP” mock communities (Fig. 4D-F). The 27 most abundant OTUs, representing on average more than 99% of the sequences in these samples, were successfully assigned to the 24 expected species. Indeed, Cornu aspersum L1 was represented by 4 OTUs and Pomatias elegans, Theba pisana L1 and Theba pisana L2 were not detected.

In both “POS” and “NREP” mock types, Cornu aspersum L1 showed the highest disparity between expected and observed sequence counts, as the results showed that it was overrepresented in the mock communities: +29% in POS-G1 and +10% on average in NREP mocks. Conversely, in the POS-G1 and POS-G2 mocks, Testacella haliotidea (-16% in POS-G1 and -3% in NREP mocks on average) and Hygromia cinctella (-18% in POS-G2 and -4% in NREP mocks on average) sequences were the most underrepresented in the observed communities (Fig. 4).

In the “controlled-feeding” controls UNI-G1 to G3 (fed only with Cornu aspersum), on average more than 99% of the sequences were assigned to Cornu aspersum L1 (one OTU for UNI-G1), to Cornu aspersum L2 (one OTU for UNI-G3) or to Cornu aspersum L1 and Cornu aspersum L2 (one OTU each for UNI-G2). In all three “controlled-feeding” controls MIX-G1 to G3, sequences from the last prey given (Arion hortensis L1, 1 OTU) were detected in up to 99% of the sequences in MIX-G3 and 54% on average in MIX-G1 and MIX-G2. MIX-G1 and MIX-G2 also showed sequences from the second prey (Cornu aspersum L1), whereas sequences from the first prey (Hygromia cinctella) were not detected in any sample. Again, sequences from the second prey Cornu aspersum L1 were highly abundant in the “MIX” samples, yielding 40% (MIX-G1) and 50% (MIX-G2) of the sequences compared to sequences from the last prey Arion hortensis L1, consumed and digested one week later (See Suppl. material 8).

The results obtained for controls were used to set the filtering threshold for OTUs. The mock communities showed that unexpected or artefactual OTUs (i.e. assigned to unexpected species or returned as “no significant similarity”) had an average occurrence of 0.1% across samples. On the other hand, unexpected OTUs from non-given prey for the “controlled-feeding” controls averaged 0.4%. These two potential thresholds were tested on PCR technical replicates from the field (V6, V29 and V53 samples). Of the 70 OTUs detected in these samples (see Suppl. material 8), 33 were assigned to 14 terrestrial gastropod species, representing on average 95.8% of the sequences. Two OTUs matched a chromosome 8 gene from Homo sapiens and Pan troglodytes, but accounted for only 0.016% of the sequences. Thirty-five OTUs, corresponding on average to 4.2% of the sequences, could not be assigned to any taxon, as no significant BLAST matches were retrieved. When applied, the 0.1% filter was not strict enough, leaving several unassigned OTUs and artefactual sequences in the samples. Conversely, the 0.4% filter allowed to discard each of these from the replicates – except for one, OTU 7 – while preserving all correctly identified OTUs. To ensure that any potential contamination was discarded and because this work was based on digestive contents, such as the “controlled-feeding” controls, the 0.4% filter was retained for further analyses.

Negative controls (BLANK 1–3) amplified very few sequences (5–8 reads in total, scattered across six species after discarding non-assigned OTUs), indicating minimal background noise and no evidence of systematic contamination. Bray-Curtis distances were calculated among these controls to provide an estimate of stochastic variation (d = 0.7), consistent with previous work using similar approaches (De Barba et al. 2014; Lopes et al. 2015; Zepeda-Mendoza et al. 2016). This threshold was used as a reference to evaluate technical replicates: any replicate with a Bray-Curtis distance greater than 0.7 would have been discarded. Replicates for field-collected digestive content samples (V6, V29 and V53) showed much lower distances (≤ 0.1) and were therefore all retained for merging. For the merging step, a conservative criterion was applied: only OTUs detected in all three technical replicates from the same sample were retained and their read counts averaged to generate consensus results (Alberdi et al. 2018; Beentjes et al. 2019) (See Suppl. material 9). Final retained OTUs were assigned to species: V6 sample, from La Rochelle, showed sequences from Boettgerilla pallens, Deroceras invadens and Oxychilus draparnaudi. V29, in Cleguer, Brittany, showed sequences from Arion distinctus, Boettgerilla pallens, Cornu aspersum L1 and the unassigned OTU 7. Finally, sequences from Ambigolimax parvipenis/valentianus, Cochlicella barbara, Cornu aspersum L1 and Deroceras invadens were recovered from the digestive contents of V53, sampled in Lampaul-Plouarzel, Brittany. A single unassigned OTU, but no possible artefactual sequences or intriguing assignments demonstrate the consistency of this filtering threshold for digestive contents from the field.

Gut contents from the 16 field samples showed 1 to 6 prey species per sample after applying the 0.4% OTU filtering threshold (See Suppl. material 8). Between 5 (in Saint-Medard-en-Jalles, region of Bordeaux) and 8 species (in Concarneau, Brittany) were recorded at the 4 sites studied (Fig. 5). 32 of the 43 (74%) detections in gut contents were species also observed at the sites during sampling phases. Deroceras invadens and Arion hortensis L1 were detected in all four sites and in the digestive contents of all sites except Cleguer, Brittany, for Arion hortensis L1. At the Cleguer site, two species not found in the field (Oxyloma elegans and Lauria cylindracea) were detected in the gut contents of two samples, and one OTU was assigned as Deroceras sp. due to a lack of taxonomic precision and the ambiguity between two potential species (Deroceras laeve and Deroceras lasithionense, the latter being endemic to the island of Crete in Greece). Finally, an OTU representing 2–4% of the sequences in two samples from Saint-Medard-en-Jalles, but below the assignment criteria, was named as “Unassigned OTU” (Fig. 5).

Figure 5.

Number of digestive contents in which each prey species was detected for the sixteen field samples, shown per sampling site. Site titles follow the nomenclature “Region, Site”. Each point corresponds to one prey species: red symbols indicate species also recorded during field sampling, whereas grey symbols indicate species not detected in the field but found in gut contents.

Universal primers detect a wide range of taxa but often lack specificity, leading to co-amplification of non-target organisms or contaminants. Their broad taxonomic range can also reduce resolution, sometimes limiting identification to family or genus level (Garwood et al. 2023). On the contrary, a prey-targeted approach maximizes the sequencing of gastropod DNA, allowing for a finer resolution of the predator’s diet. It also avoids the need for predator DNA masking with blocking primers, which, despite their effectiveness, introduce biases. Indeed, these primers can unintentionally inhibit the amplification of certain prey species or cause amplifications to fail altogether, distorting dietary analysis (Shehzad et al. 2012; Piñol et al. 2015; Robeson et al. 2018).

This study introduces a new primer pair specifically targeting short 16S DNA fragments of about 162 bp from terrestrial gastropods, representing a key advancement for the effective implementation of metabarcoding in trophic interaction and biodiversity studies. To date, only one other primer pair, developed by Taberlet et al. (2018) and already used in several studies (Mulero et al. 2021; Zhou et al. 2024), has been designed specifically for gastropods but it has not been validated in our cross-tests implying O. nungara, earthworms and gastropods DNA.

A key step in primer development is ensuring broad taxonomic coverage, i.e. the ability to amplify as many species or lineages as possible within the target taxa. Two complementary approaches, in silico amplifications based on sequence alignments and PCR amplifications on DNA from a variety of identified samples, were performed in this study and the results from both methods were only slightly divergent: among the 38 gastropod species and lineages tested, DNA amplification was expected for 32 referring to in silico tests, and PCR results showed successful amplification in 35 cases.

DNA from Pomatias elegans, Theba pisana L1, Theba pisana L2, Deroceras agreste, Deroceras reticulatum, Hygromia cinctella, and Trochulus hispidus was predicted in silico to fail amplification due to mismatches between their sequences and the retained primer pair, particularly at the 3’ end of the reverse primer. This region is critical for primer binding, and mismatches or gaps in this area can significantly reduce amplification efficiency (Elbrecht and Leese 2017). As expected, PCR amplification of DNA from Pomatias elegans and Theba pisana L2 were unsuccessful, and these species were not detected in mock community controls. For Theba pisana L1, in silico results varied depending on the sequence considered, but PCR amplification was ultimately unsuccessful and it was not detected in the mock communities.

However, despite similar mismatches, PCR amplification tests successfully amplified DNA from Deroceras agreste, Deroceras reticulatum, Hygromia cinctella, Trochulus hispidus, and all tested specimens of Arion hortensis L2, though in silico predictions for the latter showed some haplotype-dependent variability in amplification potential. These results underscore the need for experimental validation, as in silico predictions may not fully capture amplification potential. They also highlight the need to account for both interspecific and intraspecific sequence variability when designing primers for eDNA studies, ensuring their applicability across a wide range of targeted taxa (Elbrecht and Leese 2017; Taberlet et al. 2018; Bruce et al. 2021). Otherwise, the in silico analyses suggest that the VRF-VRR primer pair could be useful for the detection of terrestrial gastropod species from a wide range of habitats in France (i.e. not only anthropized areas) and in the rest of Europe in future studies using eDNA or digestive DNA. However, this still requires experimental and/or field confirmation.

Another important step in marker development is their resolution i.e. the potential of the amplified fragment to unambiguously distinguish taxa (species or genetic lineages, as described for the reference database). The phylogenetic reconstruction based on the short 16S fragment revealed that nearly all species and lineages considered could be clearly distinguished. These promising results confirm the 16S database’s value for future research and expansion. Indeed, some species still lack sufficient or even any available sequences in the literature, as is the case for Ambigolimax valentianus.

Maintaining a taxonomically verified and updated reference database is crucial for deriving ecological insights. Here, we emphasize that BLAST-based identification can be limited not by the tool itself, but by incomplete or erroneous reference databases and unresolved taxonomy, which can hinder accurate species or lineage assignment from environmental DNA. Indeed, there are two types of constraints on the identification of terrestrial gastropods from mitochondrial data: the existence of highly differentiated genetic lineages within nominal species, suggesting species complexes or cryptic species, and the existence of hybrids and introgressed forms, possibly leading to morpho-molecular incongruences.

Regarding the existence of genetic lineages within species, Pinceel et al. (2004, 2005) identified five distinct mitochondrial (mtDNA) groups (S1 to S5) within Arion subfuscus, which appears polyphyletic in phylogenetic reconstructions. Arion subfuscus samples from the French sites were assigned to group S1 using the updated reference database of these lineages. Pinceel et al. (2004, 2005) noted that the high divergence between mtDNA lineages may reflect cryptic diversity, although nuclear markers and morphological data do not provide conclusive evidence for distinct species. Similarly, genetic divergence has been observed within other species complexes, such as Trochulus hispidus (Kruckenhauser et al. 2014) which is polyphyletic and the subspecies of Theba pisana (Greve et al. 2010), underscoring the fact that multiple lineages can exist within a single species. These cases highlight the need to include all described lineages in reference databases when identifying taxa from environmental or digestive samples. Indeed, although their taxonomic status is not yet confirmed, their ecological role may be different.

Beyond cases of taxonomic identification challenges at the genetic lineage level, possible cases of misidentification could occur due to complex patterns of hybridization. They concern, for example, Arion vulgaris and Arion ater s.l. (i.e. Arion ater s.s. and the two morphotypes of A. rufus according to Reise et al. (2020). Arion vulgaris tends to replace local species in anthropized areas, producing hybrids with Arion ater s.l. and a continuum of morphological intermediates (Reise et al. 2020). In the present study, molecular analysis of the French samples with Arion rufus external morphology allowed them to be assigned to Reise et al.’s (2020) “br” lineage corresponding to the predominantly British subspecies Arion ater rufus (i.e. Rowson et al.’s (2014) Arion rufus and Chevallier’s (1972, 1974) Arion rufus rufus).

It should also be noted that studies of terrestrial gastropods have revealed remarkably high levels of both intra- and interspecific divergence in mtDNA, often exceeding 10% (Parmakelis et al. 2013). This complicates BLAST-based identification and reinforces the need for comprehensive databases and phylogenetic placement of OTUs. The 95% threshold arbitrarily chosen in this study seems rather strict in this sense, but limits misidentifications.

How to interpret high-throughput sequencing results has always been a concern in metabarcoding studies, especially when analyzing fragmented and/or partially digested DNA. Whether to rely on a quantitative or qualitative approach depends on several factors such as research aims, species-specific amplification biases, target taxa and dietary diversity, and should be considered when evaluating sequencing results (Deagle et al. 2010, 2013, 2019; De Barba et al. 2014; Alberdi et al. 2019; Casey et al. 2019; Lamb et al. 2019; Gold et al. 2023; Shelton et al. 2023).

Extracting DNA from digestive contents rather than feces implies that the sequencing results are highly dependent on the time elapsed since prey consumption, regardless of the species consumed. Consequently, the number of sequences per species does not necessarily reflect their abundance in the diet but rather indicates that the prey was consumed more recently – or conversely, earlier and more fully digested (Deagle et al. 2010, 2019). Analyzing metabarcoding data from a quantitative perspective (e.g., sequence counts or relative read abundance i.e. RRA) is also particularly challenging when considering the bias introduced by the preferential amplification of primer pairs (also referred to “species-specific bias” in Thomas et al. 2016; Yang et al. 2021; Gold et al. 2023). This bias results from several factors, including sequence motifs, GC content, fragment length, primer mismatches and DNA degradation, which affect amplification efficiency (Deiner et al. 2017; Braukmann et al. 2019). The relative abundance of detected taxa is biased because the observed sequence counts often fail to reflect the initial DNA proportions, even in controlled mock communities (Deiner et al. 2017; Shelton et al. 2023).

The results from mock communities obtained in this study are consistent with this trend. Indeed, Cornu aspersum L1 is systematically overrepresented in the samples. This suggests that this species may be preferentially amplified during the annealing step. In the MIX “controlled-feeding” controls, this species was consistently detected at a sequence number equivalent to that of Arion hortensis L1, which was consumed a week later. Interestingly, Arion hortensis L1 yielded sequencing results that were consistent with expectations in the mock communities, suggesting it is not negatively affected by amplification bias. This supports the idea that Cornu aspersum L1 DNA hybridizes more effectively and thus biases its true sequence counts in samples. Conversely, other species appear to show the opposite trend. For example, Hygromia cinctella was consistently underrepresented in mock communities, was completely absent from the “controlled-feeding” controls despite being the first prey eaten. It was also not detected in the field replicates or in the 16 field samples, despite being sampled in two of the four study sites. These results are consistent with an in silico binding score of 254, indicating that DNA from this species would theoretically not be expected to amplify efficiently with the VRF-VRR primers. However, PCR tests still showed amplification bands for this species, suggesting that although its DNA can be amplified, it is likely to be penalized during the annealing step due to mismatches. This likely results in low sequence counts, often insufficient to pass OTU filtering thresholds. The same thought process can be applied to other species such as Testacella haliotidea or Deroceras agreste, where the in silico score is either below 120 but not equal to 0, or above 120 and sequence counts in mock communities are significantly below expectations.

These results highlight a critical loss of diversity information as demonstrated by the “controlled-feeding” controls. The process is inherently multivariate, with the amplification of one taxon’s DNA influenced by the co-occurrence of others within the same samples (Shelton et al. 2023). In addition, species-specific biases remain, making it difficult to recover accurate relative abundance information from amplicon data (Yang et al. 2021), although several studies have attempted to do so in diet studies of different organisms (Port et al. 2016; Thomas et al. 2016; Ando et al. 2020). This issue is exacerbated in complex matrices, such as dietary or environmental samples, where DNA degradation, inhibitors, and preferential amplification further complicate accurate quantification (De Barba et al. 2014). Given that the quantitative aspect of metabarcoding data, particularly for digestive contents, presents numerous interpretative challenges, it was deliberately chosen to focus on the qualitative aspect, specifically on the detection (presence/absence) of prey through OTU clustering. This approach is not merely a way to circumvent the limitations of quantitative analyses but is actually better suited to the study objectives. Indeed, the goal was to gain a comprehensive understanding of the predator’s overall diet across multiple individuals and sampling sites, rather than quantifying the relative consumption of each prey item at the individual level. By emphasizing prey detection rather than abundance estimates, this qualitative framework allows for a broader and more ecologically relevant perspective on dietary diversity.

Another critical element of metabarcoding analysis is the selection of an appropriate OTU filtering threshold. This choice directly affects species detection and quantification by influencing how sequences are grouped and background noise is removed. In metabarcoding, contamination and artefacts can arise at multiple stages, including field sampling (e.g., cross-contamination between samples), DNA extraction (e.g., introduction of external DNA from reagents), PCR amplification (e.g., carry-over contamination or sequencing errors), and sequencing (e.g., through tag jumping). Low-abundance OTUs often correspond to such contaminants or errors, making the filtering strategy essential (Taberlet et al. 2018; Zinger et al. 2019). An adequate threshold reduces spurious sequences while preserving ecological signal, improving dietary diversity estimates, and helping to limit artefacts from rare events such as secondary predation.

Secondary predation, i.e., the detection of DNA from prey previously consumed by the predator’s prey, represents an additional source of potential bias (Da Silva et al. 2019; Deagle et al. 2019). Its impact depends on predator ecology: DNA from secondary ingestion is typically more degraded and occurs in lower quantities than DNA from directly consumed prey, reducing its likelihood of amplification (Hardy et al. 2017; Mata et al. 2021). For generalist or omnivorous predators, however, secondary consumption can broaden the apparent diet and introduce biases if not carefully interpreted (Bowser et al. 2013; Da Silva et al. 2019). In the case of O. nungara, secondary predation is considered unlikely. Carnivorous gastropods are rare among the species in our study area (e.g. Oxychilus spp.; Díaz et al. 2025). Other predatory species, such as Testacella spp., mainly feed on earthworms, which would not be amplified by our primers (Stokes and Hirst 1958). Therefore, any secondary predation signal is expected to be minimal relative to the main trophic signal.

In metazoan metabarcoding, a common filtering threshold of 1% of total read abundance is effective in reducing contaminants and rare erroneous sequences. Optimal thresholds, however, vary with target taxa, barcoding regions, sequencing depth, environment, and study goals (Deagle et al. 2019; Zinger et al. 2019; Ando et al. 2020). In this study, the availability of lab-fed individuals provided a robust reference for calibration, allowing thresholds to be refined with known diets—a situation not always possible in wild populations or environmental samples. In such contexts, threshold determination must rely more heavily on controls such as mock communities, positive and negative controls, and replicates to approximate artefact levels and avoid over-filtering ecologically relevant taxa.

Extensive controls—including positive and negative controls, mock communities, and technical replicates—were employed here to estimate artefact rates and refine thresholds. PCR replicates are generally recommended when feasible to increase reliability and allow stochastic amplification events to be identified, although they were not implemented in the present study (Deiner et al. 2017; Zinger et al. 2019), in order to balance the trade-off between sequencing costs and the large number of samples processed, while relying on other validation procedures. Positive controls with known DNA mixtures further guide contamination assessment and threshold selection (De Barba et al. 2014; Ando et al. 2020). Mock community analyses suggested a 0.1% threshold, whereas controlled-feeding data indicated 0.4% is more appropriate. Because DNA from mock communities is less degraded than gut content DNA, lower thresholds are reasonable to account for fewer artefacts. The 0.4% threshold performed well for digestive content samples: it removed unassigned OTUs while retaining species consistent with field observations. For instance, all species detected in samples V6, V29, and V53 matched species recorded at their respective collection sites.

The choice of OTU filtering thresholds primarily depends on the taxon and study objectives, as no universal value exists. In this study, a single threshold was applied across all samples, calibrated and validated with mock and controlled-feeding experiments to ensure that rare but ecologically relevant taxa were retained. The purpose of this threshold is not to provide a universal value, but to illustrate a rigorous approach to OTU filtering. Filtering thresholds are highly sensitive to the diversity of diets or ecosystems: in dietary metabarcoding, specialists with narrow diets may lose rare prey even under mild thresholds, exaggerating inter-individual variation, while generalists may appear artificially homogeneous if low-abundance taxa are removed (Deagle et al. 2019; Littleford-Colquhoun et al. 2022). Similarly, in environmental DNA studies, species-rich ecosystems risk losing rare taxa under strict filtering, whereas species-poor ecosystems may lose a large fraction of observed diversity even at moderate thresholds (Ando et al. 2020). Although precise thresholds depend on the system, lower thresholds may be necessary in highly diverse diets or ecosystems to retain rare taxa, whereas higher thresholds can be more appropriate in species-poor systems where low-abundance sequences are more likely to be artefacts. Optimal thresholds depend on sequencing depth, taxon abundance distributions, artefact prevalence, and study objectives; mock communities and controls are essential to guide appropriate selection (Drake et al. 2022). Importantly, for environmental samples where reference individuals or lab-fed systems are absent, empirical validation using controls, replicates, and mock communities becomes even more critical to balance artefact removal and biodiversity retention. Overall, these considerations highlight the importance of carefully calibrating filtering thresholds to the diversity and complexity of the system under study, as well as the specific ecological questions being addressed.

Diet estimation remains a major challenge in ecology. Here are presented the first results on the dietary diversity of the invasive species O. nungara, from specimens collected directly in their invaded habitats, three private gardens and a greenhouse. Prey diversity surveys were conducted at these sites, allowing for a comparison of gut content sequencing results with field observations. Species observed in the field were expected in gut contents, though not necessarily all, as O. nungara could exhibit selective feeding. If the predator is highly selective, or if certain prey are rare, the gut contents should contain fewer species than observed in the field. Poor sampling or incomplete coverage could lead to the reverse pattern. The results obtained for Saint-Medard-en-Jalles suggest that O. nungara might show selectivity, which could be influenced by factors such as prey accessibility (e.g., life cycles, defensive mechanisms) or an active feeding preference, as indicated by the observed disparity between field and gut contents species diversity. In contrast, at Concarneau, the positive correlation between species detected in the field and those found in the gut contents suggests a more opportunistic feeding behavior. At Cleguer, the detection of species in gut contents that are not recorded in field surveys indicates a potential gap in sampling effort. Indeed, field inventories focused on macrofauna, which can lead to an underrepresentation of micro-snails such as Lauria cylindracea and other small species. High-throughput sequencing of digestive DNA can therefore enhance diversity estimates, serving as a complementary tool to field surveys (Kartzinel et al. 2015; Jo et al. 2016; Amundsen and Sánchez‐Hernández 2019; De Sousa et al. 2019).

These findings highlight the importance of considering sampling effort and PCR biases when interpreting the predator’s feeding behavior. The contradictions also highlight the limited number of specimens and sites studied so far, preventing definitive ecological and trophic conclusions. However, O. nungara demonstrates a remarkable dietary diversity at this stage. Discarding the unidentified OTU and the taxon assigned to Deroceras sp., 16 gastropod species and lineages were detected in the digestive contents of only sixteen individuals, while a total of 28 species and lineages were recorded across sampled sites. Among these, there are widespread and common species across the sites, such as Cornu aspersum, Deroceras invadens, and Ambigolimax parvipenis, as well as rarer and more sporadic species like Arion hortensis L2, Arion intermedius, and Cochlicopa lubrica. A broader-scale analysis is now a thrilling perspective to determine whether O. nungara primarily adopts an opportunistic or selective feeding strategy, and what ecological impact it can have on its prey and their ecosystems.

Beyond this system, studies using environmental DNA (eDNA) to investigate terrestrial gastropods remain scarce, in contrast to aquatic mollusks where eDNA is increasingly applied to detect diverse species (Mulero et al. 2021; Zieritz et al. 2022; Vanderpool et al. 2024). This gap is particularly relevant given the growing presence of exotic predatory flatworms in regions with high gastropod endemism, such as Platydemus manokwari in Guadeloupe, Martinique, and Saint-Martin (Justine et al. 2021), and Bipalium vagum in Jamaica (Brown et al. 2022), with further expansions projected globally (Fourcade et al. 2022). Terrestrial gastropods are also consumed by a wide range of other taxa, including malacophagous beetles (Pseudoophonus rufipes, Pterostichus melanarius; Němec and Horsák 2024), snail-eating snakes (You et al. 2015; Poyarkov et al. 2022) and occasional snakes feeding on snails (Zamenis situla; Christopoulos and Kotselis 2023), rodents (Allen 2004), and omnivorous birds (Laporta-Ferreira and Salomão 2004; Morii et al. 2021). Studying these interactions with a species-level molecular assay can provide additional resolution and complement existing dietary studies.

The metabarcoding assay developed here, based on newly designed primers targeting gastropod DNA, provides a practical tool for both gut content analyses and environmental DNA studies. It can help investigate terrestrial gastropod consumption by different predator taxa, as well as explore gastropod presence in understudied habitats. While the primers were designed in the context of O. nungara predation, they could also be applied to other ecological situations involving terrestrial gastropods, offering additional opportunities for ecological research and conservation studies.

The authors have declared that no competing interests exist.

No ethical statement was reported.

No use of AI was reported.

No funding was reported.

MV and VR designed the study. Collection of the samples was assured by MV, NM, SN, LD and VR. Laboratory work was performed by MV, NM and VR. MV conducted the analyses. MV wrote the first draft of the manuscript. All authors contributed to the final version and approved the submitted manuscript.

Mathis Ventura https://orcid.org/0009-0001-6831-9654

Nicolas Mazuras https://orcid.org/0009-0004-0888-6659

Shanèze Noël https://orcid.org/0009-0001-5380-1085

Lise Dupont https://orcid.org/0000-0003-4824-2685

Virginie Roy https://orcid.org/0000-0003-1908-0074

All of the data that support the findings of this study are available in the main text or Supplementary Information.