# What Is the Population Structure of Poecilogonic Species? Evidence From Large-Scale Genotyping in a Neogastropod Lineage (Conoidea: Raphitoma)

Giacomo Chiappa, Giulia Fassio, Maria Vittoria Modica, Nicolas Puillandre, Marco Oliverio

► To cite this version:

Giacomo Chiappa, Giulia Fassio, Maria Vittoria Modica, Nicolas Puillandre, Marco Oliverio. What Is the Population Structure of Poecilogonic Species? Evidence From Large-Scale Genotyping in a Neogastropod Lineage (Conoidea: Raphitoma). Molecular Ecology, In press, 10.1111/mec.70170. hal-05359932

HAL Id: hal-05359932

https://hal.science/hal-05359932v1

Submitted on 12 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

ORIGINAL ARTICLE

OPEN ACCESS

# What Is the Population Structure of Poecilogonic Species? Evidence From Large-Scale Genotyping in a Neogastropod Lineage (Conoidea: *Raphitoma*)

Giacomo Chiappa<fnref n="1" /> ID | Giulia Fassio<fnref n="1" /> | Maria Vittoria Modica² | Nicolas Puillandre³ ID | Marco Oliverio<fnref n="1" />

**Correspondence:** Marco Oliverio (marco.oliverio@uniroma1.it)

**Received:** 27 March 2025 | **Revised:** 16 September 2025 | **Accepted:** 23 October 2025

**Handling Editor:** Angus Davison

**Funding:** This study was supported by the Sapienza University of Rome SEED PNR 2021 grant (n° B89J21032850001) to G.F., the Sapienza 2020 INV-EVO/2 grant (n° RM12017276B846BB) to M.O., the Sapienza 2021 SEED-3840974 grant to M.O., and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 865101) to N.P. This study was also partially funded by the following awards granted to G.C.: Sapienza University of Rome: Avvio alla Ricerca 2024, Malacological Society of London—Early Career Research 2023 and Unitas Malacologica—Student Research Award 2024. The support by the Doctorate School in Evolutionary and Environmental Biology (Sapienza University of Rome) to G.C. is also acknowledged. The CORSICABENTHOS expeditions were supported by the Agence Française pour la Biodiversité and Collectivité Territoriale de Corse.

**Keywords:** dispersal | poecilogony | population genetic | *Raphitoma* | SLAF-seq

## ABSTRACT

Dispersal of marine benthic invertebrates is typically dependent on the developmental mode of their pelagic larvae, which can be prolonged and based on plankton feeding (planktotrophic), or short and rely on the nutrients from the egg (non-planktotrophic). The uncommon poecilogonic species can commit to both developmental modes, with remarkable implications for their population genetics, not yet fully investigated. In this study, we obtained reduced-representation genome-wide SNP data for three species of the neogastropod genus *Raphitoma* Bellardi, 1847 from the Mediterranean, including the putatively poecilogonic species *Raphitoma philberti* (Michaud, 1829) and *R. laviae* (R. A. Philippi, 1844). A total of 80 samples were sequenced to test the poecilogony and analyse the population genetics of *R. laviae*. Although a certain degree of segregation was highlighted between the planktotrophic and non-planktotrophic samples in both species, they were ultimately found within conspecific bounds, confirming poecilogony. A set of loci that split samples with different development was identified, suggesting that a genetic component may be involved in poecilogony in both species. The population genetics of *R. laviae* fit patterns of both long- and short-dispersal species: in Corsica, with only planktotrophic developers, no geographic structure was detected, whereas in Croatia, with only non-planktotrophic developers, the geographic structure was present among localities 15–40 km apart. Notably, the species delimitation based on genome-wide SNP data was contrasting with the one assessed in past studies, reiterating that a solid taxonomy (still not achieved) is paramount to correctly interpret the evolution of larval development in this group.

<footnote n="1">¹Department of Biology and Biotechnologies “Charles Darwin”, Sapienza University of Rome, Rome, Italy | ²Department of Biology and Evolution of Marine Organisms, Stazione Zoologica Anton Dohrn, Rome, Italy | ³Institut Systématique Evolution Biodiversité (ISYEB), Muséum National D'histoire Naturelle, CNRS, Sorbonne Université, EPHE, Université Des Antilles, Paris, CP, France</footnote>

<footnote>Giacomo Chiappa and Marco Oliverio submitting author.</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.</footnote>

<footnote>© 2025 The Author(s). *Molecular Ecology* published by John Wiley & Sons Ltd.</footnote>
In marine benthic invertebrates, where adults have lim-
ited (if any) vagility, larvae are usually entitled to dispersal
and are crucial for maintaining gene flow among popula-
tions (Cowen and Sponaugle 2009; Jablonski and Lutz 1983;
Modica et al. 2017; Travis and Dytham 1998; Wang and
Altermatt 2019; Wray 1995). A variety of reproductive strat-
egies and larval developments have evolved to maximise dis-
persal, including planktotrophic larvae, which can travel far
distances while feeding autonomously, and are usually pro-
duced in large numbers to cope with mortality during the pe-
lagic larval phase. Among the advantages of a planktotrophic
development is also the reduction of intraspecific competition.
In contrast, low-dispersal reproductive strategies sacrifice the
dispersal potential to favour higher local or short-distance re-
cruitment; invertebrates specialised in this approach rely on
lecithotrophic larvae, which are usually produced in lower
amounts and get their sustenance from the egg. In benthic
invertebrates, these contrasting reproductive strategies can
affect the speciation mode and rates (Jablonski and Lutz 1983;
Krug et al. 2015; Oliverio 1996).

Each species usually relies on a single strategy, with the rare exceptions of some polychaete annelids and sacoglossan gastropods (Krug 2007; Zakas and Wares 2012), where different larval developmental modes have been observed in the same species or population, a condition known as poecilogony (Giard 1905). Being able to choose between high- and low-dispersal developmental strategies allows combining the advantages of local recruitment and of gene flow, with profound implications for the ecology of marine species and their developmental evolution (Chia et al. 1996; Zakas and Wares 2012).

In a strict definition (Hoagland and Robertson 1988; Bouchet 1989), poecilogony was adopted to encompass only the larval phenotypic variation among the offspring of a single individual. The use of a restricted definition helps avoid the traps of mixing cryptic species with distinct developmental modes. However, such a restricted definition may mask the actual variety of larval phenotypic variation, hampering a comparative approach that may be crucial for understanding the evolutionary dynamics of developmental modes (Knott and McHugh 2012). At the same time, not all developmental variations are conveniently included within poecilogony: for instance, we agree that polymorphism in dispersal without other phenotypic differences among larvae than variation in the time to hatching (e.g., Clemens-Seely and Phillips 2011; Przeslawski 2011) should not be included in the definition of poecilogony (Krug 2009). Therefore, it is useful to return to something closer to the original Giard's (1905) broad definition of poecilogony, addressing for every case its geographic-ecological scale (variation in the offspring of a single individual; offspring of different individuals from the same population; offspring of individuals from different populations; offspring of the same individual or individuals from the same population, but in different seasons or in response to different environmental conditions). The adoption of a broader definition is therefore the most suitable for a truly comparative approach. However, it should be a priority for any study on poecilogony to take into careful consideration

the taxonomy of the group, and rule out potentially cryptic
species (Knott and McHugh 2012).

Current studies on poecilogonic species aim at understanding the biological mechanisms behind poecilogony, as well as the evolutionary processes that favour a mixed reproductive strategy against a purely planktotrophic or lecithotrophic one (Krug 2009). Regrettably, only very few species are known to exhibit poecilogony, and finding new cases implies the difficult task of observing and studying larval development. Shelled gastropods are optimal model organisms as they retain their embryonic/larval shell (protoconch) as adults, and the modality of the pelagic phase can be presumed directly from protoconch length and sculpture (Jablonski and Lutz 1983; Shuto 1974). The larval shell is, in fact, multi-spiral in planktotrophic larvae and paucispiral in non-planktotrophic larvae. Based on the presence of very similar or identical adult shells but different protoconchs, poecilogony has been suggested in several species of Gastropoda. However, almost all cases of poecilogony in non-sacoglossan gastropods have been confidently dismissed, recognising in each case either distinct taxa with different development, erroneous assignment of larvae to adults, or poecilogony assumed only on the egg morphology (Bouchet 1989; Fassio et al. 2022; Hoagland and Robertson 1988). Consequently, the type of protoconchs in shelled gastropods (and thus their developmental strategy) is still considered a species-level diagnostic character (Bouchet 1989; Oliverio 1996).

Nevertheless, a convincing poecilogonic species was subsequently discovered in the family Calyptraeidae (McDonald et al. 2014), and candidate cases were highlighted more recently in the neogastropod family Raphitomidae, the most diverse taxon of Conoidea (Bouchet et al. 2011; Kantor and Taylor 2002). Species delimitation within the genus *Raphitoma* was performed to test the taxonomic status of sibling species differing in larval development (Russini et al. 2020). Two species—preliminarily identified as *R. philberti* (Michaud, 1829) and *R. laviae* (R. A. Philippi, 1844)—included specimens exhibiting both multi-spiral and paucispiral protoconch, revealing two candidate cases of poecilogony. In a subsequent study, a third species (namely “*Raphitoma* sp. C”), composed of planktotrophic specimens only, was documented as a member of this species complex (Chiappa et al. 2025). The assessment by Russini et al. (2020) was based on a Sanger dataset of mitochondrial (*cox-1*, 16S rDNA, 12S rDNA) and nuclear (*ITS2*) markers. Although the claim of poecilogony seems well supported by this data, it is also argued that the analyses may have failed to split the multi-spiral and paucispiral samples within *R. philberti* and *R. laviae* because of a lack of resolution of the dataset (Russini et al. 2020).

Confirming poecilogony is of great importance from a taxonomic perspective, as it implies that larval development, and therefore protoconch morphology, may vary in some species (potentially leading to incorrect species delimitations and biased connectivity assessments: see Pante et al. 2015), and from an evolutionary perspective, as poecilogony could be the intermediate condition leading to speciation by loss of planktotrophy in at least one of the descending lineages (Oliverio 1996). Furthermore, the few sacoglossan species known so far as poecilogonic among gastropods are all herbivores, and it has been hypothesised that poecilogony could be an adaptation to trophic source fluctuations
(Krug 2009; Strathmann 1978b, 1978a). Nevertheless, raphitomids are predators (Chiappa et al. 2024), which would represent a novel, unknown ecological combination.

Using the candidate cases of poecilogony identified by Russini et al. (2020), we performed a reduced-representation sequencing (RRS) on samples of *R. philberti*, *R. laviae* and *Raphitoma* sp. C, to test the following hypotheses:

a. That one or two poecilogonic species are identified in the genus *Raphitoma* using genome-scale data, confirming the species delimitation pattern suggested by Sanger data (Russini et al. 2020), or, alternatively, that non-poecilogonic, discrete units are identified in the species complex, thus rejecting the claims of poecilogony;

b. Whether a poecilogonic species shows a global geographic pattern of genetic variation similar to that of a high-dispersal species (due to the presence of planktotrophic larvae), or the genetic structure features typical of a low-dispersal species (due to the presence of lecithotrophic larvae), or both.

# 2 | Material and Methods

## 2.1 | Dataset

The dataset included 80 specimens identified as *Raphitoma philberti*, *R. sp*. C and *R. laviae*, following genetic assessment (Chiappa et al. 2025; Fassio et al. 2019; Russini et al. 2020) sampled from three main regions in the Mediterranean (Corsica + Sardinia, Croatia and Greece). The type of larval development was assessed by observation of protoconch morphology: multi-spiral (two or more whorls) with cancellated sculpture for planktotrophic development; paucispiral with finer, orthogonal sculpture for non-planktotrophic development. For clarity, we will directly refer to specimens as planktotrophic or non-planktotrophic hereafter, instead of referring to their protoconch morphology. The vouchers are stored in the Department of Biology and Biotechnologies "Charles Darwin", Sapienza University of Rome (BAU), or in the Muséum National d'Histoire Naturelle, Paris (MNHN).

Samples were processed following a specific locus amplified fragment sequencing (SLAF-Seq): DNA was double digested using enzymes RsaI and HaeIII, fragments were range selected, barcoded, amplified via PCR and sequenced in 150 PE on a NovaSeq machine (BMKGene—Münster, Germany). For additional details on the SLAF-Seq protocol, see Sun et al. (2013). Following raw reads de-multiplexing and quality trimming, the SLAF tags were clustered and polymorphic tags were identified. For each tag, reads were mapped on the highest depth sequence using BWA-MEM v0.7.10 with -r 789 (Li and Durbin 2009) and SNP calling was performed by taking the output intersection of GATK v4.2.1.0 UnifiedGenotyper (McKenna et al. 2010) and samtools v1.9 mpileup (Li et al. 2009) commands, with default settings. The resulting dataset was filtered using vcftools v0.1.16 (Danecek et al. 2011), allowing a minimum allele frequency of 0.05, a minimum sequencing depth of 6x and only biallelic SNPs. Samples with more than 95% missing data were removed.

Bioinformatic analyses requiring heavy computing resources were performed on Terastat2 (Bompiani et al. 2020). Two types of datasets were finally produced:

i. A SNP dataset, including only allele information with less than 10% missing data. A second SNP dataset, including only a subset of samples with both larval development of *R. laviae*, was also assembled for population genetics (see results). A single random SNP per locus was selected for analyses requiring unlinked data.

ii. Three SLAF tag sets of ~100 high-quality tags each: one with tags with less than 10% missing data in a subset of 32 samples representing all lineages while maximising completeness (dataset SLAF90); two randomly resampled from tags with less than 25% missing data shared by all samples (datasets SLAF75.1 and SLAF75.2).

## 2.2 | Validation of Poecilogony

The SNP dataset was analysed to check if the species as identified previously with mitochondrial and nuclear markers (Chiappa et al. 2025; Russini et al. 2020) were consistently retrieved by a reduced genome dataset. Multiple approaches were adopted: a structure analysis performed with the R package LEA (Frichot and François 2015), estimating the number of ancestral populations (*K*) that minimises cross-entropy (testing 1 to 10*K*); a shared ancestry analysis computed using the R package fineSTRUCTURE (Lawson et al. 2012), ordering loci according to the Linkage Disequilibrium using the provided script for unmapped data; a Principal Component Analysis (PCA) performed on R with the prcomp command (R Core Team 2023) and plotted using the R package ggplot2 (Wickham 2016).

When a large number of loci is used, species delimitation methods are known to possibly “oversplit” (Sukumaran and Knowles 2017), as they may retrieve populations as distinct species regardless of gene flow (Leaché et al. 2019; Rannala and Yang 2020). In order to evaluate whether the populations identified from the SNP dataset analyses were conspecific, the SLAF tags datasets were analysed with a Bayesian approach using the species delimitation algorithm bpp 4.0 (Flouri et al. 2018; Yang 2015), by estimating the genealogical divergence index (gdi; Jackson et al. 2017), which gives the probability that samples in two groups can coalesce in a given time (Poelstra et al. 2021; Rannala and Yang 2020). Following Jorna et al. (2021), we performed multiple analyses on three SLAF tag subsets and took the consensus among them as the final results, allowing the retrieval of the most robust signal with approachable computational costs. Additional information on the calibration of this analysis is in Appendix S1-1. We used Jackson's et al. (2017) rule of thumb, where gdi < 0.2 indicates a single species, gdi > 0.7 indicates two species and values in-between represent an ambiguous grey zone. Nucleotide diversity ($\pi$) and divergence ($D_{xy}$), and net divergence ($D_a$) among the units retrieved from the structure analyses were also computed for each SLAF tag subset (Nei and Li 1979) using MEGA X (Kumar et al. 2018).

To retrieve loci possibly related to larval development, $F_{ST}$ per site was computed between a subsample (see Results) of
planktotrophic and non-planktotrophic specimens of *R. laviae* using Weir and Cockerham's (1984) formula with vcftools v0.1.16 (Danecek et al. 2011). Negative estimates were set to 0. Loci including sites with scores above the 95% quantile were selected and filtered from the SNP dataset. PCAs were then performed again to visualise the effects of this set of loci on the total variance.

## 2.3 | Divergence Rates and Time Estimation

To date the historical events within the group of interest, we performed a time-calibrated ML phylogenetic reconstruction of Raphitomidae using a *cox-1* + 16S rDNA dataset comprising 28 species of eight genera, after Almón et al. (2024).

Four constrained age intervals based on fossil data were used in the Raphitomidae analysis: (1) the origin of the family is assumed to have occurred 47.8–37.8 million years ago (Mya) based on the first known representative, the Middle Eocene *Pleurotomella polycolpa* (Cossmann, 1889) (Gougerot and Le Renard 1981; Lozouet 2017), a dating confirmed by molecular assessment (Abdelkrim et al. 2018); (2) the age of the mrca of European raphitomids (including the lineages corresponding to the modern *Leufroyia*, *Raphitoma* and *Cyrillia*) can be dated at 28–23 Mya, based on Upper Oligocene Raphitomidae with morphology more similar to *Leufroyia* than to *Raphitoma* (e.g., *R. subfragilis* Lozouet 2017); (3) the age of the mrca of the *Raphitoma*-like lineages (*Raphitoma* and *Cyrillia*) can be placed in the Middle Miocene (16.3–12.8 Mya) when species recalling modern forms – like the Badenian Paratethys *Raphitoma praehispida* (Boettger, 1906) – are first found; and (4) the split between the lineage leading to the Plio-Pleistocene *Raphitoma histrix* and then to the Recent *R. pseudohystrix* (with its oldest representative, the Tortonian *Raphitoma dellabellaorum* Landau et al. 2020), and the rest of the NE Atlantic and Mediterranean *Raphitoma*, is dated at 11.6–7.2 Mya (Landau et al. 2020). Additional information on this analysis is in Appendix S1-2.

## 2.4 | Population Genetics of a Poecilogonic Species

A SNP dataset including a subset of planktotrophic and non-planktotrophic samples in *R. laviae* (see Results) was analysed to explore the population genetics of a poecilogonic species, including a total of 47 specimens representing both larval developmental strategies. Instead of looking for the most likely number of populations in the dataset (e.g., the one with the lowest cross-entropy), we preliminarily delimited populations according to the three main sampling regions: Corsica, Croatia and Greece. Then, we looked further at the micro geographic scale by partitioning the dataset according to the distribution of the sampling sites within Corsica (with samples from the Northern, Western and Southern coasts), and within Croatia (with samples from Punta Križa and the region spanning ~25 km around Ravni Kotari). Population genetics statistics were computed using the “populations” module of STACKS v2.76 (Catchen et al. 2011, 2013). A structure analysis was performed using LEA (Frichot and François 2015) setting the initial number of ancestral populations to three to test for geographic structure in the sampled areas. Isolation by Distance (Wright 1943) was tested using

geodetic pairwise oceanographic distances measured on QGis v. 3.12.2 and a pairwise nucleotide distance matrix computed using vcf2dis (available at https://github.com/BGI-shenzhen/VCF2Dis). Mantel tests (Mantel 1967) were performed with the Spearman model and 10,000 replicates with the R packages vegan v2.6.4 (Dixon 2003) and tidyverse v2.0.0 (Wickham et al. 2019).

# 3 | Results

## 3.1 | Dataset

The dataset composition (Table 1) included 13 specimens of *Raphitoma philberti* from Sardinia (*n* = 1), Corsica (*n* = 3), Greece (*n* = 2) and Croatia (*n* = 7), nine of *Raphitoma* sp. C. from Sardinia (*n* = 1), Corsica (*n* = 6) and Greece (*n* = 2), and 58 of *Raphitoma laviae* from Corsica (*n* = 35), Greece (*n* = 2) and Croatia (*n* = 21). The attribution of larval development by proto-conch observation revealed 7 non-planktotrophic specimens in *R. philberti*, including one from Corsica and 11 in *R. laviae*, all from Croatia.

A total of 1.57 billion reads were outputted after sequencing (19,165,278 reads per sample on average), and 287,633 polymorphic SLAF tags with x10.17 mean depth were retrieved in total, resulting in 1,329,954 total SNPs. Three samples of *R. philberti* with > 95% missing data were removed. The average miss rate of tags in the remaining samples was 70%, with values ranging from 60% to 84%. The SNP dataset comprised a total of 10,111 SNPs from 1789 loci. The high-quality SLAF tags datasets comprised 119 tags with > 90% completeness in a subset of 32 samples (dataset SLAF90), and 202 tags with > 75% completeness shared by all samples (split in SLAF75.1 and SLAF75.2).

## 3.2 | Validation of Poecilogony

The optimal number of ancestral populations (*K*) resulting from the LEA analysis (Figure 1, additional plots provided in Figure S2) was five, corresponding to *Raphitoma* sp. C, *R. philberti* and three populations within *R. laviae*, hereafter called L1 (with samples from all localities, all planktotrophic), L2 (with samples from Corsica, all planktotrophic) and L3 (with samples from Croatia, all non-planktotrophic). With *K*=3, the analysis did not retrieve the three species as predetermined, but lumped *R. philberti* with *R. sp*. C and split L1 from the rest of *R. laviae*. Samples of population L1 were the first to split with *K*=2. Three specimens of *R. philberti*, all planktotrophic, had mixed ancestry. In the fineSTRUCTURE analysis (Figure 1), high coancestry values were estimated within the ancestral populations from the LEA analyses, and a geographical structure was present in non-planktotrophic populations.

The variance explained by the first two principal components of the PCA was 15.27% and 11.72%, rapidly decreasing to 6.67%, 5.63%, and below 3% in the following PCs (Figure 2, additional plots provided in Figure S3). Three main groups were identified for each species by plotting the samples onto the first two axes, not corresponding to the three species, but lumping *R. philberti* and *R. laviae* and separating L1 from the
<table><caption>TABLE 1 | List of *Raphitoma* included in this study.</caption><thead><tr><th>Sample ID</th><th>Larval develop</th><th>Region</th><th>Locality</th><th>Pop</th><th>Total reads</th><th>SLAF number</th><th>Miss rate</th><th>NCBI Biosample</th></tr></thead><tbody><tr><td colspan="9"><em>Raphitoma philberti</em></td></tr><tr><td>BAU 1888**</td><td>non-pl</td><td>Italy</td><td>Taranto</td><td>—</td><td>43,580</td><td>4143</td><td>99.8</td><td>SAMN47597843</td></tr><tr><td>BAU 2241.2*</td><td>—</td><td>Croatia</td><td>Biograd</td><td>P</td><td>17,627,994</td><td>524,635</td><td>74.2</td><td>SAMN47597844</td></tr><tr><td>BAU 2248.2*</td><td>—</td><td>Croatia</td><td>Vrsi</td><td>P</td><td>33,829,538</td><td>546,128</td><td>76.1</td><td>SAMN47597845</td></tr><tr><td>BAU 2252.2</td><td>non-pl</td><td>Croatia</td><td>Zaton</td><td>P</td><td>19,806,532</td><td>501,471</td><td>76.8</td><td>SAMN47597846</td></tr><tr><td>BAU 2261.2</td><td>pl</td><td>Croatia</td><td>Biograd</td><td>P</td><td>21,060,522</td><td>376,808</td><td>84.3</td><td>SAMN47597847</td></tr><tr><td>BAU 3046*</td><td>pl</td><td>Greece</td><td>Vai Is.</td><td>P</td><td>19,121,126</td><td>566,470</td><td>71.8</td><td>SAMN47597848</td></tr><tr><td>BAU 3352.1</td><td>non-pl</td><td>Croatia</td><td>Punta Križa</td><td>P</td><td>8,678,402</td><td>402,672</td><td>79</td><td>SAMN47597849</td></tr><tr><td>BAU 3352.2**</td><td>non-pl</td><td>Croatia</td><td>Punta Križa</td><td>—</td><td>921,222</td><td>27,967</td><td>99</td><td>SAMN47597850</td></tr><tr><td>BAU 3352.3**</td><td>non-pl</td><td>Croatia</td><td>Punta Križa</td><td>—</td><td>7,319,556</td><td>104,765</td><td>95.1</td><td>SAMN47597851</td></tr><tr><td>BAU 4301</td><td>—</td><td>Italy</td><td>Cor-N</td><td>P</td><td>20,657,512</td><td>513,510</td><td>77.5</td><td>SAMN47597852</td></tr><tr><td>MNHN-1M-2019-6068*</td><td>non-pl</td><td>Corsica</td><td>Cor-N</td><td>P</td><td>23,640,440</td><td>528,642</td><td>76.4</td><td>SAMN47597853</td></tr><tr><td>MNHN-1M-2019-17259*</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>P</td><td>19,339,866</td><td>490,878</td><td>75.2</td><td>SAMN47597854</td></tr><tr><td>MNHN-1M-2019-18330*</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>P</td><td>17,869,208</td><td>511,477</td><td>75.5</td><td>SAMN47597855</td></tr><tr><td colspan="9"><em>Raphitoma sp. C</em></td></tr><tr><td>BAU 3545.1</td><td>pl</td><td>Greece</td><td>Leros Is.</td><td>C</td><td>21,833,116</td><td>414,178</td><td>82.3</td><td>SAMN47597856</td></tr><tr><td>BAU 3545.2*</td><td>pl</td><td>Greece</td><td>Leros Is.</td><td>C</td><td>25,183,848</td><td>507,083</td><td>76.4</td><td>SAMN47597857</td></tr><tr><td>BAU 4300*</td><td>—</td><td>Sardinia</td><td>Cor-S</td><td>C</td><td>30,721,574</td><td>480,241</td><td>81</td><td>SAMN47597858</td></tr><tr><td>MNHN-1M-2019-13292*</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>C</td><td>17,271,624</td><td>462,920</td><td>79.8</td><td>SAMN47597859</td></tr><tr><td>MNHN-1M-2019-13944*</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>C</td><td>22,629,494</td><td>514,508</td><td>78.4</td><td>SAMN47597860</td></tr><tr><td>MNHN-1M-2019-14607*</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>C</td><td>38,144,402</td><td>620,817</td><td>75</td><td>SAMN47597861</td></tr><tr><td>MNHN-1M-2019-17885</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>C</td><td>14,681,778</td><td>456,363</td><td>78.7</td><td>SAMN47597862</td></tr><tr><td>MNHN-1M-2019-18329</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>C</td><td>14,339,814</td><td>472,601</td><td>76.8</td><td>SAMN47597863</td></tr><tr><td>MNHN-1M-2019-18507*</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>C</td><td>13,864,880</td><td>467,896</td><td>75.8</td><td>SAMN47597864</td></tr><tr><td colspan="9"><em>Raphitoma laviae</em></td></tr><tr><td>BAU 2238.1*</td><td>non-pl</td><td>Croatia</td><td>Biograd</td><td>L3</td><td>17,862,178</td><td>703,998</td><td>64.5</td><td>SAMN47597865</td></tr></tbody></table>

(Continues)
<table><caption>TABLE 1 | (Continued)</caption><thead><tr><th>Sample ID</th><th>Larval develop</th><th>Region</th><th>Locality</th><th>Pop</th><th>Total reads</th><th>SLAF number</th><th>Miss rate</th><th>NCBI Biosample</th></tr></thead><tbody><tr><td>BAU 2238.3*</td><td>non-pl</td><td>Croatia</td><td>Biograd</td><td>L3</td><td>23,416,166</td><td>746,026</td><td>62.3</td><td>SAMN47597866</td></tr><tr><td>BAU 2241.1</td><td>non-pl</td><td>Croatia</td><td>Biograd</td><td>L3</td><td>25,428,916</td><td>593,366</td><td>73.8</td><td>SAMN47597867</td></tr><tr><td>BAU 2243.1*</td><td>pl</td><td>Croatia</td><td>Sukošan</td><td>L1</td><td>17,823,974</td><td>666,224</td><td>63</td><td>SAMN47597868</td></tr><tr><td>BAU 2243.2*</td><td>pl</td><td>Croatia</td><td>Sukošan</td><td>L1</td><td>19,820,810</td><td>643,919</td><td>65.6</td><td>SAMN47597869</td></tr><tr><td>BAU 2243.3</td><td>pl</td><td>Croatia</td><td>Sukošan</td><td>L1</td><td>20,512,588</td><td>619,748</td><td>68.4</td><td>SAMN47597870</td></tr><tr><td>BAU 2246.2</td><td>pl</td><td>Croatia</td><td>Zaton</td><td>L1</td><td>20,789,562</td><td>627,163</td><td>67.8</td><td>SAMN47597871</td></tr><tr><td>BAU 2246.3</td><td>pl</td><td>Croatia</td><td>Zaton</td><td>L1</td><td>19,563,282</td><td>564,078</td><td>72.4</td><td>SAMN47597872</td></tr><tr><td>BAU 2246.4</td><td>pl</td><td>Croatia</td><td>Zaton</td><td>L1</td><td>11,130,408</td><td>549,984</td><td>69</td><td>SAMN47597873</td></tr><tr><td>BAU 2247</td><td>non-pl</td><td>Croatia</td><td>Biograd</td><td>L3</td><td>15,149,460</td><td>637,950</td><td>68.2</td><td>SAMN47597874</td></tr><tr><td>BAU 2249.1*</td><td>non-pl</td><td>Croatia</td><td>Sukošan</td><td>L3</td><td>12,490,178</td><td>570,818</td><td>70.2</td><td>SAMN47597875</td></tr><tr><td>BAU 2249.2</td><td>non-pl</td><td>Croatia</td><td>Sukošan</td><td>L3</td><td>13,146,510</td><td>570,592</td><td>71.8</td><td>SAMN47597876</td></tr><tr><td>BAU 2249.3</td><td>non-pl</td><td>Croatia</td><td>Sukošan</td><td>L3</td><td>16,383,468</td><td>477,206</td><td>76.5</td><td>SAMN47597877</td></tr><tr><td>BAU 2251.2</td><td>pl</td><td>Croatia</td><td>Turanj</td><td>L1</td><td>16,126,654</td><td>538,295</td><td>72.7</td><td>SAMN47597878</td></tr><tr><td>BAU 2255*</td><td>non-pl</td><td>Croatia</td><td>Sabunike</td><td>L3</td><td>16,715,602</td><td>646,719</td><td>66.6</td><td>SAMN47597879</td></tr><tr><td>BAU 2267</td><td>non-pl</td><td>Croatia</td><td>Sabunike</td><td>L3</td><td>13,918,112</td><td>557,429</td><td>71</td><td>SAMN47597880</td></tr><tr><td>BAU 2268.2</td><td>non-pl</td><td>Croatia</td><td>Biograd</td><td>L3</td><td>12,450,926</td><td>573,088</td><td>72.6</td><td>SAMN47597881</td></tr><tr><td>BAU 2270.1</td><td>pl</td><td>Croatia</td><td>Biograd</td><td>L1</td><td>11,036,018</td><td>534,254</td><td>70.1</td><td>SAMN47597882</td></tr><tr><td>BAU 2270.2</td><td>pl</td><td>Croatia</td><td>Biograd</td><td>L1</td><td>11,598,640</td><td>526,275</td><td>71.5</td><td>SAMN47597883</td></tr><tr><td>BAU 2363.2</td><td>non-pl</td><td>Croatia</td><td>Biograd</td><td>L3</td><td>25,574,130</td><td>590,484</td><td>73</td><td>SAMN47597884</td></tr><tr><td>BAU 3358.2</td><td>pl</td><td>Croatia</td><td>Punta Križa</td><td>L1</td><td>17,136,184</td><td>553,114</td><td>71.2</td><td>SAMN47597885</td></tr><tr><td>BAU 3550*</td><td>pl</td><td>Greece</td><td>Leros Is.</td><td>*</td><td>25,231,780</td><td>674,351</td><td>67.6</td><td>SAMN47597886</td></tr><tr><td>BAU 3551</td><td>pl</td><td>Greece</td><td>Leros Is.</td><td>L1</td><td>16,148,146</td><td>539,193</td><td>72.2</td><td>SAMN47597887</td></tr><tr><td>MNHN-1M-2019-4014*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>25,825,530</td><td>763,427</td><td>60.8</td><td>SAMN47597888</td></tr><tr><td>MNHN-1M-2019-4020*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>13,659,394</td><td>587,962</td><td>69</td><td>SAMN47597889</td></tr><tr><td>MNHN-1M-2019-4066</td><td>—</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>13,766,474</td><td>583,750</td><td>69.2</td><td>SAMN47597890</td></tr><tr><td>MNHN-1M-2019-4068*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>15,591,622</td><td>678,562</td><td>62.4</td><td>SAMN47597891</td></tr></tbody></table>

(Continues)
<table><caption>TABLE 1 | (Continued)</caption><thead><tr><th>Sample ID</th><th>Larval develop</th><th>Region</th><th>Locality</th><th>Pop</th><th>Total reads</th><th>SLAF number</th><th>Miss rate</th><th>NCBI Biosample</th></tr></thead><tbody><tr><td>MNHN-IM-2019-4087</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>21,972,282</td><td>644,356</td><td>67.2</td><td>SAMN47597892</td></tr><tr><td>MNHN-IM-2019-4094*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>*</td><td>27,515,058</td><td>693,000</td><td>64.7</td><td>SAMN47597893</td></tr><tr><td>MNHN-IM-2019-5300</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L1</td><td>14,162,368</td><td>572,206</td><td>70.1</td><td>SAMN47597894</td></tr><tr><td>MNHN-IM-2019-5306</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L1</td><td>13,501,756</td><td>569,194</td><td>69.7</td><td>SAMN47597895</td></tr><tr><td>MNHN-IM-2019-5342</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>14,987,102</td><td>651,221</td><td>66.2</td><td>SAMN47597896</td></tr><tr><td>MNHN-IM-2019-6022</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>14,320,572</td><td>607,104</td><td>68.2</td><td>SAMN47597897</td></tr><tr><td>MNHN-IM-2019-6024</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>15,467,482</td><td>644,948</td><td>65.6</td><td>SAMN47597898</td></tr><tr><td>MNHN-IM-2019-6025*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L1</td><td>20,655,040</td><td>645,116</td><td>66.7</td><td>SAMN47597899</td></tr><tr><td>MNHN-IM-2019-6026</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>11,783,640</td><td>583,195</td><td>68.1</td><td>SAMN47597900</td></tr><tr><td>MNHN-IM-2019-6034*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>28,775,204</td><td>736,832</td><td>62.4</td><td>SAMN47597901</td></tr><tr><td>MNHN-IM-2019-6035*</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>23,458,906</td><td>692,349</td><td>64.4</td><td>SAMN47597902</td></tr><tr><td>MNHN-IM-2019-6067</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>13,916,090</td><td>607,161</td><td>69.1</td><td>SAMN47597903</td></tr><tr><td>MNHN-IM-2019-6069</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>14,397,358</td><td>613,299</td><td>69.2</td><td>SAMN47597904</td></tr><tr><td>MNHN-IM-2019-6090</td><td>pl</td><td>Corsica</td><td>Cor-N</td><td>L2</td><td>14,914,300</td><td>635,732</td><td>67.2</td><td>SAMN47597905</td></tr><tr><td>MNHN-IM-2019-12026</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>L2</td><td>11,158,584</td><td>537,028</td><td>70.4</td><td>SAMN47597906</td></tr><tr><td>MNHN-IM-2019-12932*</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>L2</td><td>18,311,870</td><td>637,215</td><td>66.4</td><td>SAMN47597907</td></tr><tr><td>MNHN-IM-2019-13703*</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>L2</td><td>36,880,462</td><td>820,390</td><td>59.6</td><td>SAMN47597908</td></tr><tr><td>MNHN-IM-2019-13959*</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>L1</td><td>19,585,334</td><td>635,153</td><td>66.9</td><td>SAMN47597909</td></tr><tr><td>MNHN-IM-2019-14537</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>L2</td><td>18,516,596</td><td>684,876</td><td>64.6</td><td>SAMN47597910</td></tr><tr><td>MNHN-IM-2019-14594</td><td>pl</td><td>Corsica</td><td>Cor-S</td><td>L2</td><td>18,330,742</td><td>667,113</td><td>66.9</td><td>SAMN47597911</td></tr><tr><td>MNHN-IM-2019-16533</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>18,990,116</td><td>672,372</td><td>65.8</td><td>SAMN47597912</td></tr><tr><td>MNHN-IM-2019-16949*</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L1</td><td>19,386,880</td><td>615,702</td><td>67.7</td><td>SAMN47597913</td></tr><tr><td>MNHN-IM-2019-17159</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>17,390,480</td><td>670,593</td><td>64.6</td><td>SAMN47597914</td></tr><tr><td>MNHN-IM-2019-17246*</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>40,912,682</td><td>765,119</td><td>61.8</td><td>SAMN47597915</td></tr><tr><td>MNHN-IM-2019-17406</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>19,034,446</td><td>679,430</td><td>65.5</td><td>SAMN47597916</td></tr><tr><td>MNHN-IM-2019-17663*</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>14,375,434</td><td>628,467</td><td>66.5</td><td>SAMN47597917</td></tr></tbody></table>

(Continues)
<table><caption>TABLE 1 | (Continued)</caption><thead><tr><th>Sample ID</th><th>Larval develop</th><th>Region</th><th>Locality</th><th>Pop</th><th>Total reads</th><th>SLAF number</th><th>Miss rate</th><th>NCBI Biosample</th></tr></thead><tbody><tr><td>MNHN-1M-2019-17690</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>12,192,298</td><td>589,108</td><td>67.9</td><td>SAMN47597918</td></tr><tr><td>MNHN-1M-2019-17888</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>17,024,372</td><td>655,702</td><td>66.4</td><td>SAMN47597919</td></tr><tr><td>MNHN-1M-2019-17942</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>10,882,632</td><td>562,118</td><td>70.1</td><td>SAMN47597920</td></tr><tr><td>MNHN-1M-2019-18144</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>14,249,610</td><td>632,210</td><td>66.8</td><td>SAMN47597921</td></tr><tr><td>MNHN-1M-2019-18145</td><td>pl</td><td>Corsica</td><td>Cor-W</td><td>L2</td><td>17,820,428</td><td>677,517</td><td>64.8</td><td>SAMN47597922</td></tr></tbody></table>

Note: Asterisks mark samples included in the SLAF90 reduced dataset (*) or removed from the analysis (**). “Pop” refers to the results of the structure analyses (see Figures 1 and 2). “Larval development” was inferred from protoconch morphology (pl: Planktotrophic; non-pl: Non-planktotrophic). “Miss rate” refers to the percentage of missing tags. For additional sampling information refer to Chiappa et al. (2024) Supporting Information.

rest of *R. laviae*. Groups including planktotrophic specimens were very condensed, despite including samples from distant localities, whereas those that included non-planktotrophic samples were more spread, often in accordance with their geographic distribution. The samples with different larval development were separated, with the exception of a single non-planktotrophic specimen of *R. philberti* from Croatia (BAU 3352.1). In the plot resulting from the third and fourth PCs, the non-planktotrophic samples were more evidently separated from the planktotrophic, with the planktotrophic specimens of *R. philberti* closer to the planktotrophic of *R. laviae* rather than the conspecific, non-planktotrophic samples, with the only exception of a planktotrophic specimen of *R. philberti* from Croatia (BAU 2261.2). This plot also showed a broader spread of non-planktotrophic samples compared to the others.

The estimation of the genealogical divergence index (gdi) on three independent high-quality datasets (Table 2) showed values (0.08–0.11) consistently below the 0.2 threshold when merging the samples belonging to units L2 and L3, indicating that the two planktotrophic and non-planktotrophic populations retrieved in *R. laviae*, although separated in the structure analysis and slightly apart in the PCA, belong to the same species. Including samples of the group L1, however, yielded values (0.19–0.26) not consistently below the threshold, thus not supporting conspecificity of L1 with L2 and L3. All remaining values, including those corresponding to the three preliminary species, ranged between 0.18 and 0.27. Nucleotide distance ($D_{xy}$) among units ranged from 0.0052 to 0.0021 (Table 3). In particular, $D_{xy}$ between L1 and the rest of *R. laviae* ranged from 0.0032 to 0.0041, whereas $D_{xy}$ between L2 and L3 was lower, ranging from 0.0021 to 0.0027. Similarly, net divergence ($D_{a}$) ranged from 0.0607 to 0.0004 (Table 3) and was higher between L1 and L2 or L3 (0.0028–0.0018) compared to the $D_{a}$ between L2 and L3 (0.0009–0.0004).

A total of 114 loci were identified with sites showing $F_{st}$ values above the 95% percentile threshold value (0.732) between the *R. laviae* planktotrophic samples in L2 and the non-planktotrophic samples in L3. The PCA plots (Figure 3) obtained using loci below the threshold were very similar to those obtained with the all-loci dataset. In *R. laviae*, one group of planktotrophic specimens (L2) was well separated from the other planktotrophic (L1), and was instead very close to a group of non-planktotrophic samples (L3). In the PCA using only the set of loci above the $F_{st}$ threshold (i.e., those with a divergence correlated with a difference in larval development in L2 v. L3), that same group became instead evidently closer to samples of the same development; furthermore, with these loci, the separation of non-planktotrophic and planktotrophic specimens also involved *R. philberti*, which became fully split, without the exceptions retrieved using the full dataset.

## 3.3 | Divergence Rates and Time Estimation

The likelihood-ratio tests rejected the strict molecular clock hypothesis ($p=0.0000$). For all major lineages, a clock-like behaviour was also rejected on most nodes, suggesting the need for methods to estimate the age of nodes that allowed for different rates of molecular evolution among lineages. Cross-validation analyses using a range of smoothing values from 0.1 to 10 yielded
<figure><img src="image_3.png" /><figcaption>**FIGURE 1** | (A) Ancestry proportion of samples to five ancestral populations (K). (B) Matrix of co-ancestry between samples. Colours correspond to ancestral populations retrieved in the structure analysis; all nodes had maximum support except the one crossed.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**FIGURE 2** | PCA plots with the first four principal components. (A) Proportion of variance explained by the first ten components. (B, C) Colour legend corresponding to groups retrieved in the structure analyses (Figure 1). (D, E) Colour legend corresponding to sampling localities (COR: Corsica; GR: Greece; KR: Croatia). Larval development inferred from protoconch morphology (*pl*: Planktotrophic; *non-pl*: Non-planktotrophic; *u*: Unknown).</figcaption></figure>

the best score for the lowest values. This indicated the need for a smoothing value that would allow for abrupt disparities in rates among close lineages. We thus used a value of 0.1 (however, we also used values up to 1, which produced similar results – data not shown).

Dating the Raphitomidae ML tree (Figure S1) yielded an age of 1.66 Mya (0.83–2.98) for the mrca of the three clades corresponding to *Raphitoma philberti*, *R. sp*. C and *R. laviae*. The mrca of

the sibling species *R. cordieri* (Payraveau, 1826) and *R. horrida* (Monterosato, 1884)—a planktotrophic and a non-planktotrophic species, respectively—was dated at 3.3 Mya (1.81–8.81).

## 3.4 | Population Genetics of a Poecilogonic Species

The patterns resulting from the SNP dataset analyses did not completely match those retrieved in the past from mitochondrial
<table><caption>**TABLE 2** | Gdi scores estimated from three SLAF tags datasets of 119 (SLAF90) and 100 (SLAF75) loci using three different priors for root age (τ) and genetic diversity (θ).</caption><thead><tr><th></th><th>τ=5θ</th><th>τ=0.5θ</th><th>τ=0.05θ</th></tr></thead><tbody><tr><td>SLAF90</td><td></td><td></td><td></td></tr><tr><td>L2L3</td><td>0.1088</td><td>0.1125</td><td>0.1082</td></tr><tr><td>L1L2L3</td><td>0.2326</td><td>0.2120</td><td>0.2613</td></tr><tr><td>PC</td><td>0.2383</td><td>0.2470</td><td>0.2390</td></tr><tr><td>SLAF75.1</td><td></td><td></td><td></td></tr><tr><td>L2L3</td><td>0.0852</td><td>0.081</td><td>0.078</td></tr><tr><td>L1L2L3</td><td>0.2203</td><td>0.1967</td><td>0.1963</td></tr><tr><td>PC</td><td>0.2328</td><td>0.2202</td><td>0.2309</td></tr><tr><td>SLAF75.2</td><td></td><td></td><td></td></tr><tr><td>L2L3</td><td>0.0971</td><td>0.0958</td><td>0.0957</td></tr><tr><td>L1L2L3</td><td>0.2453</td><td>0.2348</td><td>0.2336</td></tr><tr><td>PC</td><td>0.2731</td><td>0.2508</td><td>0.2711</td></tr></tbody></table>

Note: Letters refer to the units as derived from the structure analysis (see Figure 1). Values represent the probability that samples in two groups can coalesce in a given time. Red scores are below the threshold for conspecific populations (gdi < 0.2).

<table><caption>**TABLE 3** | Genetic distance estimates from three SLAF datasets of 119 (SLAF90) and 100 (SLAF75) loci among units as derived from the structure analysis (see Figure 1). Left: Nucleotide diversity within (π) and divergence between (Dxy) units. Right: Net divergence (Dn) between units.</caption><thead><tr><th>SLAF90</th><th>π</th><th>C</th><th>P</th><th>L1</th><th>L2</th><th>SLAF90</th><th>C</th><th>P</th><th>L1</th><th>L2</th></tr></thead><tbody><tr><td>C</td><td>0.00269</td><td></td><td></td><td></td><td></td><td>P</td><td>0.00154</td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.00271</td><td>0.00424</td><td></td><td></td><td></td><td>L1</td><td>0.00302</td><td>0.00247</td><td></td><td></td></tr><tr><td>L1</td><td>0.00099</td><td>0.00486</td><td>0.00432</td><td></td><td></td><td>L2</td><td>0.00188</td><td>0.00146</td><td>0.00186</td><td></td></tr><tr><td>L2</td><td>0.00220</td><td>0.00433</td><td>0.00392</td><td>0.00346</td><td></td><td>L3</td><td>0.00238</td><td>0.00192</td><td>0.00213</td><td>0.00042</td></tr><tr><td>L3</td><td>0.00124</td><td>0.00435</td><td>0.00390</td><td>0.00324</td><td>0.00214</td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>SLAF75.1</th><th>π</th><th>C</th><th>P</th><th>L1</th><th>L2</th><th>SLAF75.1</th><th>C</th><th>P</th><th>L1</th><th>L2</th></tr><tr><td>C</td><td>0.00227</td><td></td><td></td><td></td><td></td><td>P</td><td>0.00164</td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.00207</td><td>0.00381</td><td></td><td></td><td></td><td>L1</td><td>0.00270</td><td>0.00235</td><td></td><td></td></tr><tr><td>L1</td><td>0.00119</td><td>0.00443</td><td>0.00398</td><td></td><td></td><td>L2</td><td>0.00183</td><td>0.00176</td><td>0.00193</td><td></td></tr><tr><td>L2</td><td>0.00214</td><td>0.00403</td><td>0.00386</td><td>0.00360</td><td></td><td>L3</td><td>0.00241</td><td>0.00218</td><td>0.00217</td><td>0.00061</td></tr><tr><td>L3</td><td>0.00134</td><td>0.00422</td><td>0.00388</td><td>0.00344</td><td>0.00235</td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>SLAF75.2</th><th>π</th><th>C</th><th>P</th><th>L1</th><th>L2</th><th>SLAF75.2</th><th>C</th><th>P</th><th>L1</th><th>L2</th></tr><tr><td>C</td><td>0.00268</td><td></td><td></td><td></td><td></td><td>P</td><td>0.00199</td><td></td><td></td><td></td></tr><tr><td>P</td><td>0.00323</td><td>0.00495</td><td></td><td></td><td></td><td>L1</td><td>0.00320</td><td>0.00295</td><td></td><td></td></tr><tr><td>L1</td><td>0.00132</td><td>0.00520</td><td>0.00522</td><td></td><td></td><td>L2</td><td>0.00234</td><td>0.00207</td><td>0.00205</td><td></td></tr><tr><td>L2</td><td>0.00249</td><td>0.00492</td><td>0.00493</td><td>0.00395</td><td></td><td>L3</td><td>0.00304</td><td>0.00264</td><td>0.00278</td><td>0.00086</td></tr><tr><td>L3</td><td>0.00127</td><td>0.00501</td><td>0.00489</td><td>0.00407</td><td>0.00274</td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

and nuclear markers. However, at least one unit (L2 + L3) comprising enough planktotrophic and non-planktotrophic specimens was identified. We chose this subsample as our study unit to analyse the effect of poecilogony on population genetics

(Figure 4, additional plots in Figure S4). The samples spanned all areas of sampling, with 30 specimens from Corsica (planktotrophic), 1 from Greece (planktotrophic) and 11 from Croatia (non-planktotrophic).
<figure><img src="image_5.png" /><figcaption>**FIGURE 3** | Left: *F*<sub>st</sub> per site between planktotrophic (L2) and non-planktotrophic (L3) samples in *R. laviae*. Dotted line represents the 95% percentile above mean. Right: PCA plots obtained with loci below and above the threshold, including only *R. laviae* and the complete dataset. Colour legend corresponding to groups retrieved in the structure analyses (Figure 1). Bar plots show the variance explained by the first six PCs.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**FIGURE 4** | (A) Results of LEA population ancestry analysis of *R. laviae* (L2+L3) based on geographic localities. Specimens with non-planktotrophic development marked with ‘x’. (B) Sampling map. (C) *F*<sub>st</sub> statistics, with values below the 0.05 threshold highlighted in red and Mantel test correlation index with asterisks representing statistical significance (one: *P*-val < 0.01, two: *P*-val < 0.001, three: *P*-val < 0.0001). (D) Population genetic statistics (COR: Corsica; KR: Croatia).</figcaption></figure>

*F*<sub>st</sub> values were higher than 0.05 between Croatia and the other regions, and also within Croatia between localities 15 and 40 km apart. *F*<sub>st</sub> values between Corsica and Greece and within Corsica were < 0.02, suggesting a single metapopulation. Observed heterozygosity in Corsica (0.20 ± 0.05) and Croatia (0.15 ± 0.11) was congruent with the expected values. The population in Croatia had a lower percentage of private alleles (663 v. 2630) and polymorphic loci (48% v. 75%). The inbreeding coefficient (*F*<sub>IS</sub>) was

also lower (0.04 v. 0.07), but fell within the same confidence interval. The populations of this poecilogonic species retrieved from the structure analysis partially corresponded to those preliminarily assumed on a geographic basis. The non-planktotrophic samples from Croatia all shared the same ancestral population and Greece was not set apart from Corsica, where two different haplotypes were retrieved. Within Corsica, no geographic pattern was present and Isolation by Distance was not retrieved.
On the contrary, a clear structure emerged within Croatia, with a significant (*p* < 0.0001) correlation index of 0.77.

# 4 | Discussion

## 4.1 | Poecilogony in Raphitoma

Compared to the previous molecular assessments of this group (Chiappa et al. 2025; Russini et al. 2020), our large-scale genotyping yielded a molecular framework with a much higher resolution, retrieving differences among samples that did not emerge from the previous dataset of few Sanger genetic markers. With only a few exceptions in *R. philberti* (which were resolved when using a set of loci correlated with differences in development), samples of different larval development were segregated, and not mixed as appeared using predominantly mitochondrial markers. Nevertheless, the *gdi* and genetic distance analyses suggested that the distance retrieved is not enough to confidently separate them as distinct species, as occurred with other evaluations of poecilogony cases in shelled gastropods (Fassio et al. 2022). In fact, all *gdi* values between the two units of *R. laviae* with different development L2 and L3 were well below the threshold for conspecificity (Jackson et al. 2017), and the same result was verified from nucleotide diversity estimates, following Roux et al. (2016). Therefore, we were able to assume the specimens in the L2 and L3 groups as comprising a single poecilogonic species, which would result in a typical case of geographic poecilogony. One group, ranging from Corsica to Greece, shows the metapopulational features of planktotrophic developers: it has a wide range and no geographic structure. Although different haplotypes are retrieved in Corsica (Figure 4), these specimens may have originated from very far regions of the Mediterranean (e.g., Eastern v. Western, Alboran Sea v. Central Mediterranean). The population in Croatia, on the contrary, shows a distinctive geographic structure, even at a small geographic scale, which fits an intracapsular development with very limited dispersal abilities.

In Croatia, the lack of recruitment of planktotrophic specimens of this species may indicate that planktotrophic specimens from outside Croatia cannot reach this region anymore. The isolation of this restricted region may be the reason that promoted the loss of planktotrophy in this population, which combined with the prolonged absence of gene flow from outside may eventually prompt a speciation event (Krug 2011; Oliverio 1996). In this case, the barrier to gene flow between the planktotrophic and non-planktotrophic units in *R. laviae* may be merely of a geographic nature. An alternative explanation is that a drastically lower proportion of planktotrophic specimens can be found in the area at a given time, as their offspring would rapidly switch to a non-planktotrophic development. This would explain why they are not easily found sympatrically; however, in this case, one would expect not all non-planktotrophic specimens to be strictly related, since some of them should have originated from planktotrophic ancestors carrying haplotypes of distant populations. This did not emerge from our assessment, but we cannot completely exclude this hypothesis given the sample size of our Croatian dataset, and the absence of sampling from intermediate regions between those where our planktotrophic and non-planktotrophic populations reside. If there is a transition

from planktotrophic to non-planktotrophic development in the Adriatic (and/or the eastern Mediterranean: Oliverio 1996) favoured by the environmental conditions, with traces of gene flow from the planktotrophic towards the non-planktotrophic population, it may result from additional sampling in the Southern Adriatic, the Aegean Sea and the Levant basin.

In either case, if loss of planktotrophy were an irreversible process, strong enough to determine by itself a reproductive barrier between samples with different larval development (e.g., by involving asynchronous reproductive events), then poecilogony in those gastropods would be a very transient phenomenon (thus explaining its rarity), with the rapid emergence of reproductively isolated lineages and the following possible destinies: (a) both lineages persist as distinct species, (b) the non-planktotrophic lineage becomes extinct, as it is more vulnerable to changes in environmental conditions and should therefore be naturally more prone to extinction (especially if limited to a restricted range) and (c) the planktotrophic lineage becomes extinct due to counter-selection over a wide range, and is replaced by the non-planktotrophic one. At any short time frame, the perception of the researcher would be that poecilogony does not exist (Bouchet 1989).

Possible drivers of a switch to the non-planktotrophic strategy could also be related to trophic resources for planktotrophic larvae being insufficient or unpredictable (Krug 2009; Oyarzun and Strathmann 2011). Additionally, or alternatively, sea level lowering may have caused significant changes in hydrographical conditions, producing changes in circulation patterns and confinement of subbasins. Noteworthy, the estimated divergence time of the two species in the pair *R. cordieri*–*R. horrida* (~3.3 mya) correlates with the Late Pliocene Marine Isotope Stage M2, a globally recognisable cooling event, which has disturbed an otherwise relatively warm background climate state, and may have involved a global sea-level fall of between 40 m and 60 m (Dolan et al. 2015). Similarly, the differentiation within the *R. laviae*–*R. philberti* complex has started in the Early Pleistocene, i.e., the period following the onset of the glaciations and their southward extension (Thunnel and Douglas 1983, 1984). There is evidence that during glacial periods the Mediterranean underwent severe changes in the circulation patterns, often showing repeated strong confinement of the subsidiary basins (Med v. Atlantic, eastern v. western Med) (Guo et al. 2020; Kaboth et al. 2017), enhanced by sea-level lowering that caused reductions of the Sicily Channel up to three quarter width (Bethoux 1979, 1984). Under such conditions, fluctuations in the energy (food) input, restricted areas and higher predatory pressure may have easily been the main factors counterselecting the planktotrophic larvae (Strathmann 1978b, 1978a, 1985) in the Eastern Mediterranean.

Finally, the pattern emerged that *R. philberti* had many similarities to *R. laviae*, as non-planktotrophic specimens were mostly spread according to sampling locality, while the planktotrophic ones were not. Most non-planktotrophic specimens were also from Croatia, suggesting similar conditions promoting loss of planktotrophy in the genus *Raphitoma*, but two of them were from the Tyrrhenian Sea. Therefore, in this species, we observed the sympatric occurrence of specimens with different larval development, and we also retrieved one non-planktotrophic sample
genetically mixed with planktotrophic ones. Furthermore, when we used the same loci that in *R. laviae* correlated to different larval development, we observed a clear segregation of samples also in *R. philberti*. This result suggests not only that the switch to a non-planktotrophic strategy may have a genetic component, but also that it could involve the same genes in both species. Unfortunately, the reduced sampling size for this taxon prevented us from analysing in depth its population genetics and drawing evolutionary conclusions on how it is affected by poecilogony.

## 4.2 | Alternative Scenarios

Despite our focus being on the validation and study of poecilogony, our analyses also highlighted persisting uncertainties in the current species delimitation of the group. Our starting hypothesis of delimitation was based primarily on *cox-1* mitochondrial data, which is a robust marker for such aims (Miralles et al. 2024) with wide applications in gastropods (Fedosov and Puillandre 2012; Furfaro et al. 2023; Nocella et al. 2024), although in the most recent analyses (Chiappa et al. 2025) data did not show a clear barcoding gap in this group. Even using hundreds of loci, part of the species delimitation of this group remains in a grey zone (Jackson et al. 2017; Roux et al. 2016). The high co-ancestry values between *R. philberti* and *Raphitoma* sp. C, their distribution in the PCA and structure analyses and their relatively low genetic divergence, do not unequivocally support their separation as distinct species or their pooling into a single one. Conversely, *R. laviae* as now conceived could comprise multiple taxa (e.g., one planktotrophic species [L1] with wide distribution and one poecilogonic species [L2 + L3], which is the scenario explored in this study).

Evolutionary hypotheses on a group with such quickly changing taxonomy are inevitably burdened by a certain degree of uncertainty (Pante et al. 2015). This is even worse when attempting to validate poecilogony, as the species boundaries between and among planktotrophs and non-planktotrophs cannot easily be defined without performing cross-fertilisation experiments. This study explores a likely scenario, assuming in *R. laviae* two different species, one planktotrophic with a wide distribution and one poecilogonic, with different larval development in the two geographic localities explored. However, the taxonomy of this complex is yet to be fully resolved, which should be the aim of a species delimitation study done in an integrative approach. For example, another possible interpretation of the data is that all specimens of *R. laviae* included in the analysis were conspecific, including those in group L1. In this scenario, there is a single metapopulation in the Mediterranean, due to the high dispersal of planktotrophic specimens that may have maintained important gene flow over large distances, maintaining a condition of panmixia. Yet, the strong genetic divergence among groups would indicate a sort of isolation between the sympatric (often syntopic) subpopulations. This would be compatible with the high-dispersal capability of planktotrophic specimens. In fact, although being sampled together, these specimens may have originated in different, distant regions in the Eastern, Western and Central Mediterranean. Finally, in this scenario, sympatric specimens of differing larval morphology coexist in Croatia. Population genetic statistics for this scenario were computed as they may provide the basis for compelling interpretations of the

origin and effect of poecilogony in future studies, and are shown in Figures S5 and S6. Another possibility is that L2 and L3 constitute two different species, and there is no poecilogony in *R. laviae*. The low *gdi* and distance estimates may be explained by a recent speciation event. In this scenario, an allopatric speciation event occurred in the Adriatic Sea, with loss of planktotrophy in the new lineage.

In order to face these unsolved issues, as well as further pursue the study of the genetics underlying poecilogony, it will be necessary to take a step further and analyse the genomes of specimens of this complex. This would allow for the performance of association studies (Uffelmann et al. 2021) and draw evolutionary scenarios that better explain the reasons behind the contrasting levels of polymorphism detected by mitochondrial and nuclear data (Singhal et al. 2025). A second promising strategy to study poecilogony could be a comparative transcriptomic approach applied to eggs/embryos (Lesoway et al. 2016), and to the female gonads, to explore if and how the mother effect is implied in the offspring larval development.

### Author Contributions

Conceptualisation: G.C.; G.F.; M.O.; Data curation: G.C.; N.P.; Formal analysis: G.C.; Investigation: G.C.; G.F.; M.O.; Methodology: G.C.; G.F.; M.V.M.; M.O.; Software: G.C.; Visualisation: G.C.; G.F.; M.O.; Writing original draft: G.C.; Writing - review and editing: G.C.; G.F.; M.V.M.; N.P.; M.O.; Supervision: G.F.; M.V.M.; M.O.; Validation: G.C.; G.F.; M.O.; Resources: G.F.; N.P.; M.O.; Funding acquisition: G.C.; G.F.; N.P.; M.O.; Project administration: M.O.

### Acknowledgements

We are deeply grateful to three anonymous reviewers for the extensive work on this study. MNHN specimens were collected during scientific expeditions under the auspices of the Our Planet Reviewed program (CORSICABENTHOS 1-3). We thank Barbara Buge, Mélanie van Weddingen and Dario Zuccon for their help with the MNHN collection. We are thankful to Paolo Mariottini (Università Roma Tre) for support during sampling in Sardinia, to Paolo Franchini (Tuscia University of Viterbo) and to Paolo Colangelo (National Research Council, Rome) for their advice with analyses. Paul Zaharias (MNHN, Paris) and Yasunori Kano (University of Tokyo) are thanked for their valuable comments on preliminary results. We extend our gratitude to Marco Sollitto (University of Trieste) and Chiara Mancino (Sapienza University of Rome) for the help provided in several analyses. Open access publishing facilitated by Universita degli Studi di Roma La Sapienza, as part of the Wiley - CRUI-CARE agreement.

### Disclosure

Benefit-sharing: The material was collected during several expeditions operated under the regulations then in force in the countries in question and satisfies the conditions set by the Nagoya Protocol for access to genetic resources. All relevant data and results produced in this study have been shared on appropriate public databases as described above.

### Conflicts of Interest

The authors declare no conflicts of interest.

### Data Availability Statement

Raw sequence reads have been deposited in the SRA repository (www.ncbi.nlm.nih.gov/sra) with BioProject ID PRJNA1242061. The
Biosample accession numbers are reported in Table 1. The unfiltered variant call file (vcf) was uploaded to Zenodo with the following doi: 10.5281/zenodo.17086275.

## References

Abdelkrim, J., L. Aznar-Cormano, A. E. Fedosov, et al. 2018. "Exon-Capture-Based Phylogeny and Diversification of the Venomous Gastropods (Neogastropoda, Conoidea)." *Molecular Biology and Evolution* 35: 2355-2374.

Almón, B., J. Pérez-Dieste, J. E. Trigo, M. Oliverio, and G. Chiappa. 2024. "First Record of *Raphitoma ephesina* Pusateri, Giannuzzi Savelli & Stahlschmidt, 2017 (Gastropoda: Conoidea: Raphitomidae) From the Atlantic, With an Assessment of Its Systematics." *Cahiers de Biologie Marine* 65: 195-205. https://doi.org/10.21411/CBM.A.199E1833.

Bethoux, J. P. 1979. "Le Régime de la Méditerranée au Cours de Périodes Glaciaires." *Il Nuovo Cimento* 2, no. 2: 117-126.

Bethoux, J. P. 1984. "Paléo-Hydrologie de la Méditerranée au Cours Des Derniers 20 000 Ans." *Oceanologica Acta* 7, no. 1: 43-48.

Bompiani, E., U. F. Petrillo, G. J. Lasinio, and F. Palini. 2020. "High-Performance Computing With Terastat." Intl Conf on Dependable, Autonomic and Secure Computing, Intl Conf on Pervasive Intelligence and Computing, Intl Conf on Cloud and Big Data Computing, Intl Conf on Cyber Science and Technology Congress (DASC/PiCom/CBDCom/CyberSciTech), 499-506.

Bouchet, P. 1989. "A Review of Poecilogony in Gastropods." *Journal of Molluscan Studies* 55, no. 1: 67-78.

Bouchet, P., Y. I. Kantor, A. Sysoev, and N. Puillandre. 2011. "A New Operational Classification of the Conoidea (Gastropoda)." *Journal of Molluscan Studies* 77, no. 3: 273-308. https://doi.org/10.1093/mollus/eyr017.

Catchen, J. M., A. Amores, P. Hohenlohe, W. Cresko, and J. H. Postlethwait. 2011. "Stacks: Building and Genotyping Loci *De Novo* From Short-Read Sequences." *G3: Genes, Genomes, Genetics* 1, no. 3: 171-182. https://doi.org/10.1534/g3.111.000240.

Catchen, J. M., P. A. Hohenlohe, S. Bassham, A. Amores, and W. A. Cresko. 2013. "Stacks: An Analysis Tool Set for Population Genomics." *Molecular Ecology* 22, no. 11: 3124-3140. https://doi.org/10.1111/mec.12354.

Chia, F.-S., G. Gibson, and P.-Y. Qian. 1996. "Poecilogony as a Reproductive Strategy of Marine Invertebrates." *Oceanologica Acta* 19, no. 3-4: 203-208.

Chiappa, G., G. Fassio, M. V. Modica, et al. 2025. "Sexual Dimorphism in *Raphitoma purpurea* (Neogastropoda: Raphitomidae)." *Journal of Molluscan Studies* 91, no. 3: eyaf002. https://doi.org/10.1093/mollus/eyaf002.

Chiappa, G., G. Fassio, M. V. Modica, and M. Oliverio. 2024. "Potential Ancestral Conoidean Toxins in the Venom Cocktail of the Carnivorous Snail *Raphitoma purpurea* (Montagu, 1803) (Neogastropoda: Raphitomidae)." *Toxins* 16, no. 8: 348. https://doi.org/10.3390/toxin s16080348.

Clemens-Seely, K., and N. E. Phillips. 2011. "Effects of Temperature on Hatching Time and Hatchling Proportions in a Poecilogonous Population of *Haminoea zelandiae*." *Biological Bulletin* 221, no. 2: 189-196. https://doi.org/10.1086/BBLv221n2p189.

Cowen, R. K., and S. Sponaugle. 2009. "Larval Dispersal and Marine Population Connectivity." *Annual Review of Marine Science* 1, no. 1: 443-466. https://doi.org/10.1146/annurev.marine.010908.163757.

Danecek, P., A. Auton, G. Abecasis, et al. 2011. "The Variant Call Format and VCFtools." *Bioinformatics* 27, no. 15: 2156-2158. https://doi.org/10.1093/bioinformatics/btr330.

Dixon, P. 2003. "VEGAN, a Package of R Functions for Community Ecology." *Journal of Vegetation Science* 14, no. 6: 927-930. https://doi.org/10.1111/j.1654-1103.2003.tb02228.x.

Dolan, A. M., A. M. Haywood, S. J. Hunter, et al. 2015. "Modelling the Enigmatic Late Pliocene Glacial Event—Marine Isotope Stage M2." *Global and Planetary Change* 128: 47-60. https://doi.org/10.1016/j.gloplacha.2015.02.001.

Fassio, G., P. Bouchet, M. Oliverio, and E. E. Strong. 2022. "Re-Evaluating the Case for Poecilogony in the Gastropod *Planaxis sulcatus* (Cerithioidea, Planaxidae)." *BMC Ecology and Evolution* 22, no. 1: 13. https://doi.org/10.1186/s12862-022-01961-7.

Fassio, G., V. Russini, F. Pusateri, et al. 2019. "An Assessment of *Raphitoma* and Allied Genera (Neogastropoda: Raphitomidae)." *Journal of Molluscan Studies* 85, no. 4: 414-425. https://doi.org/10.1093/mollus/eyz022.

Fedosov, A. E., and N. Puillandre. 2012. "Phylogeny and Taxonomy of the *Kermia*-pseudodaphnella (Mollusca: Gastropoda: Raphitomidae) Genus Complex: A Remarkable Radiation via Diversification of Larval Development." *Systematics and Biodiversity* 10, no. 4: 447-477. https://doi.org/10.1080/14772000.2012.753137.

Flouri, T., X. Jiao, B. Rannala, and Z. Yang. 2018. "Species Tree Inference With BPP Using Genomic Sequences and the Multispecies Coalescent." *Molecular Biology and Evolution* 35, no. 10: 2585-2593. https://doi.org/10.1093/molbev/msy147.

Frichot, E., and O. François. 2015. "LEA: An R Package for Landscape and Ecological Association Studies." *Methods in Ecology and Evolution* 6, no. 8: 925-929. https://doi.org/10.1111/2041-210X.12382.

Furfaro, G., W. Renda, G. Nardi, and S. Giacobbe. 2023. "Integrative Taxonomy of the Bubble Snails (*Cephalaspidea*, *Heterobranchia*) Inhabiting a Promising Study Area: The Coastal Sicilian Faro Lake (Southern Italy)." *Water (Basel)* 15, no. 14: 2504. https://doi.org/10.3390/w15142504.

Giard, A. 1905. "La Poecilogonie." In *Compte-Rendu Des Séances Du Sixième Congrès International de Zoologie*, edited by M. Bedot, 617-646. Kündig & Fils.

Gougerot, L., and J. Le Renard. 1981. "Clefs de Détermination Des Petites Espèces de Gastéropodes de l'Eocène du Bassin Parisien. 15—Les Genres *Raphitoma* et *Mangelia*." *Cahiers Des Naturalistes* 36: 69-82.

Guo, Q., B. Li, A. H. L. Voelker, and J.-K. Kim. 2020. "Mediterranean Outflow Water Dynamics Across the Middle Pleistocene Transition Based on a 1.3 Million-Year Benthic Foraminiferal Record Off the Portuguese Margin." *Quaternary Science Reviews* 247: 106567. https://doi.org/10.1016/j.quascirev.2020.106567.

Hoagland, K. E., and R. Robertson. 1988. "An Assessment of Poecilogony in Marine Invertebrates: Phenomenon or Fantasy?" *Biological Bulletin* 174, no. 2: 109-125. https://doi.org/10.2307/1541778.

Jablonski, D., and R. A. Lutz. 1983. "Larval Ecology of Marine Benthic Invertebrates: Paleobiological Implications." *Biological Reviews* 58, no. 1: 21-89. https://doi.org/10.1111/j.1469-185X.1983.tb00380.x.

Jackson, N. D., B. C. Carstens, A. E. Morales, and B. C. O'Meara. 2017. "Species Delimitation With Gene Flow." *Systematic Biology* 66, no. 5: 799-812. https://doi.org/10.1093/sysbio/syw117.

Jorna, J., J. B. Linde, P. C. Searle, et al. 2021. "Species Boundaries in the Messy Middle—A Genome-Scale Validation of Species Delimitation in a Recently Diverged Lineage of Coastal Fog Desert Lichen Fungi." *Ecology and Evolution* 11, no. 24: 18615-18632. https://doi.org/10.1002/ece3.8467.

Kaboth, S., P. Grunert, and L. Lourens. 2017. "Mediterranean Outflow Water Variability During the Early Pleistocene." *Climate of the Past* 13, no. 8: 1023-1035. https://doi.org/10.5194/cp-13-1023-2017.

Kantor, Y. I., and J. D. Taylor. 2002. "Foregut Anatomy and Relationships of Raphitomine Gastropods (Gastropoda: Conoidea: Raphitominae)." *Bollettino Malacologico* 38, no. 4: 83-110.
Knott, K. E., and D. McHugh. 2012. "Introduction to Symposium: Poecilogony–A Window on Larval Evolutionary Transitions in Marine Invertebrates." *Integrative and Comparative Biology* 52, no. 1: 120–127.
https://doi.org/10.1093/icb/ics037.

Krug, P. J. 2007. "Poecilogony and Larval Ecology in the Gastropod Genus *Alderia*." *American Malacological Bulletin* 23, no. 1: 99–111.
https://doi.org/10.4003/0740-2783-23-1.99.

Krug, P. J. 2009. "Not My "Type": Larval Dispersal Dimorphisms and Bet-Hedging in Opisthobranch Life Histories." *Biological Bulletin* 216, no. 3: 355–372.
https://doi.org/10.1086/BBLv216n3p355.

Krug, P. J. 2011. "Patterns of Speciation in Marine Gastropods: A Review of the Phylogenetic Evidence for Localized Radiations in the Sea." *American Malacological Bulletin* 29, no. 1–2: 169–186.
https://doi.org/10.4003/006.029.0210.

Krug, P. J., J. E. Vendetti, R. A. Ellingson, et al. 2015. "Species Selection Favors Dispersive Life Histories in Sea Slugs, but Higher Per-Offspring Investment Drives Shifts to Short-Lived Larvae." *Systematic Biology* 64, no. 6: 983–999.
https://doi.org/10.1093/sysbio/syv046.

Kumar, S., G. Stecher, M. Li, C. Knyaz, and K. Tamura. 2018. "MEGA X: Molecular Evolutionary Genetics Analysis Across Computing Platforms." *Molecular Biology and Evolution* 35, no. 6: 1547–1549.
https://doi.org/10.1093/molbev/msy096.

Landau, B. M., F. V. Dingenen, and L. Ceulemans. 2020. "The Upper Miocene Gastropods of Northwestern France, 5. Conoidea." *Cainozoic Research* 20, no. 1: 3–107.

Lawson, D. J., G. Hellenthal, S. Myers, and D. Falush. 2012. "Inference of Population Structure Using Dense Haplotype Data." *PLoS Genetics* 8, no. 1: e1002453.
https://doi.org/10.1371/journal.pgen.1002453.

Leaché, A. D., T. Zhu, B. Rannala, and Z. Yang. 2019. "The Spectre of Too Many Species." *Systematic Biology* 68, no. 1: 168–181.
https://doi.org/10.1093/sysbio/syy051.

Lesoway, M. P., E. Abouheif, and R. Collin. 2016. "Comparative Transcriptomics of Alternative Developmental Phenotypes in a Marine Gastropod." *Journal of Experimental Zoology Part B: Molecular and Developmental Evolution* 326, no. 3: 151–167.
https://doi.org/10.1002/jez.b.22674.

Li, H., and R. Durbin. 2009. "Fast and Accurate Short Read Alignment With Burrows–Wheeler Transform." *Bioinformatics* 25, no. 14: 1754–1760.
https://doi.org/10.1093/bioinformatics/btp324.

Li, H., B. Handsaker, A. Wysoker, et al. 2009. "The Sequence Alignment/ Map Format and SAMtools." *Bioinformatics* 25, no. 16: 2078–2079.
https://doi.org/10.1093/bioinformatics/btp352.

Lozouet, P. 2017. "Les Conoidea de l'Oligocène Supérieur (Chattien) Du Bassin de l'Adour (Sud-Ouest de la France)." *Cossmanniana* 19: 3–180.

Mantel, N. 1967. "The Detection of Disease Clustering and a Generalized Regression Approach." *Cancer Research* 27, no. 2: 209–220.

McDonald, K. A., R. Collin, and M. P. Lesoway. 2014. "Poecilogony in the Caenogastropod *Calyitraea lichen* (Mollusca: Gastropoda)." *Invertebrate Biology* 133, no. 3: 213–220.
https://doi.org/10.1111/ivb.12057.

McKenna, A., M. Hanna, E. Banks, et al. 2010. "The Genome Analysis Toolkit: A MapReduce Framework for Analyzing Next-Generation DNA Sequencing Data." *Genome Research* 20, no. 9: 1297–1303.
https://doi.org/10.1101/gr.107524.110.

Miralles, A., N. Puillandre, and M. Vences. 2024. "DNA Barcoding in Species Delimitation: From Genetic Distances to Integrative Taxonomy." In *DNA Barcoding*, edited by R. DeSalle, vol. 2744, 77–104. Springer US.
https://doi.org/10.1007/978-1-0716-3581-0_4.

Modica, M. V., V. Russini, G. Fassio, and M. Oliverio. 2017. "Do Larval Types Affect Genetic Connectivity at Sea? Testing Hypothesis in Two Sibling Marine Gastropods With Contrasting Larval Development."

*Marine Environmental Research* 127: 92–101. https://doi.org/10.1016/j.marenvres.2017.04.001.

Nei, M., and W. H. Li. 1979. "Mathematical Model for Studying Genetic Variation in Terms of Restriction Endonucleases." *Proceedings of the National Academy of Sciences* 76, no. 10: 5269–5273. https://doi.org/10.1073/pnas.76.10.5269.

Nocella, E., G. Fassio, D. Zuccon, N. Puillandre, M. V. Modica, and M. Oliverio. 2024. "From Coral Reefs Into the Abyss: The Evolution of Corallivory in the Coralliophilinae (Neogastropoda, Muricidae)." *Coral Reefs* 43, no. 5: 1285–1302. https://doi.org/10.1007/s00338-024-02537-1.

Oliverio, M. 1996. "Life-Histories, Speciation, and Biodiversity in Mediterranean Prosobranch Gastropods." *Vie et Milieu / Life & Environment* 46, no. 2: 163–169.

Oyarzun, F. X., and R. R. Strathmann. 2011. "Plasticity of Hatching and the Duration of Planktonic Development in Marine Invertebrates." *Integrative and Comparative Biology* 51, no. 1: 81–90. https://doi.org/10.1093/icb/icr009.

Pante, E., N. Puillandre, A. Viricel, et al. 2015. "Species Are Hypotheses: Avoid Connectivity Assessments Based on Pillars of Sand." *Molecular Ecology* 24, no. 3: 525–544. https://doi.org/10.1111/mec.13048.

Poelstra, J. W., J. Salmona, G. P. Tiley, et al. 2021. "Cryptic Patterns of Speciation in Cryptic Primates: Microendemic Mouse Lemurs and the Multispecies Coalescent." *Systematic Biology* 70, no. 2: 203–218. https://doi.org/10.1093/sysbio/syaa053.

Przeslawski, R. 2011. "Notes on the Egg Capsule and Variable Embryonic Development of *Nerita melanotragus* (Gastropoda: Neritidae)." *Molluscan Research* 31, no. 3: 152. https://doi.org/10.11646/mr.31.3.2.

R Core Team. 2023. *R: A Language and Environment for Statistical Computing*. R Foundation for Statistical Computing. https://www.R-project.org/.

Rannala, B., and Z. Yang. 2020. "Species Delimitation." In *Phylogenetics in the Genomic Era*, edited by C. Scornavacca, F. Delsuc, and N. Galtier, 1–18. https://hal.inria.fr/PGE.

Roux, C., C. Fraïsse, J. Romiguier, Y. Anciaux, N. Galtier, and N. Bierne. 2016. "Shedding Light on the Grey Zone of Speciation Along a Continuum of Genomic Divergence." *PLoS Biology* 14, no. 12: e2000234. https://doi.org/10.1371/journal.pbio.2000234.

Russini, V., R. Giannuzzi-Savelli, F. Pusateri, et al. 2020. "Candidate Cases of Poecilogony in Neogastropoda: Implications for the Systematics of the Genus *Raphitoma* Bellardi, 1847." *Invertebrate Systematics* 34, no. 3: 293–318. https://doi.org/10.1071/IS19039.

Shuto, T. 1974. "Larval Ecology of Prosobranch Gastropods and Its Bearing on Biogeography and Paleontology." *Lethaia* 7, no. 3: 239–256. https://doi.org/10.1111/j.1502-3931.1974.tb00899.x.

Singhal, S., A. D. Leaché, M. K. Fujita, C. D. Cadena, and F. Zapata. 2025. "A Genomic Perspective on Species Delimitation." *Annual Review of Ecology, Evolution, and Systematics* 56: 467–489. https://doi.org/10.1146/annurev-ecolsys-102723-055311.

Strathmann, R. R. 1978a. "Progressive Vacating of Adaptive Types During the Phanerozoic." *Evolution* 32, no. 4: 907–914. https://doi.org/10.1111/j.1558-5646.1978.tb04643.x.

Strathmann, R. R. 1978b. "The Evolution and Loss of Feeding Larval Stages of Marine Invertebrates." *Evolution* 32, no. 4: 894. https://doi.org/10.2307/2407502.

Strathmann, R. R. 1985. "Feeding and Nonfeeding Larval Development and Life-History Evolution in Marine Invertebrates." *Annual Review of Ecology and Systematics* 16, no. 1: 339–361. https://doi.org/10.1146/annurev.es.16.110185.002011.

Sukumaran, J., and L. L. Knowles. 2017. "Multispecies Coalescent Delimits Structure, Not Species." *Proceedings of the National Academy of Sciences* 114, no. 7: 1607–1612. https://doi.org/10.1073/pnas.1607921114.
Sun, X., D. Liu, X. Zhang, et al. 2013. "SLAF-Seq: An Efficient Method of Large-Scale De Novo SNP Discovery and Genotyping Using High-Throughput Sequencing." *PLoS One* 8, no. 3: e58700. https://doi.org/10.1371/journal.pone.0058700.

Thunnel, R. C., and F. W. Douglas. 1983. "Paleotemperature and Paleosalinity History of the Eastern Mediterranean Sea During the Late Quaternary." *Palaeogeography, Palaeoclimatology, Palaeoecology* 44: 23–39.

Thunnel, R. C., and F. W. Douglas. 1984. "Anoxic Events in the Mediterranean Sea in Relation to the Evolution of the Late Neogene Climates." *Marine Geology* 59: 105–134.

Travis, J. M. J., and C. Dytham. 1998. "The Evolution of Dispersal in a Metapopulation: A Spatially Explicit, Individual-Based Model." *Proceedings of the Royal Society of London, Series B: Biological Sciences* 265, no. 1390: 17–23. https://doi.org/10.1098/rspb.1998.0258.

Uffelmann, E., Q. Q. Huang, N. S. Munung, et al. 2021. "Genome-Wide Association Studies." *Nature Reviews Methods Primers* 1, no. 1: 59. https://doi.org/10.1038/s43586-021-00056-9.

Wang, S., and F. Altermatt. 2019. "Metapopulations Revisited: The Area-Dependence of Dispersal Matters." *Ecology* 100, no. 9: e02792. https://doi.org/10.1002/ecy.2792.

Weir, B. S., and C. C. Cockerham. 1984. "Estimating F-Statistics for the Analysis of Population Structure." *Evolution* 38, no. 6: 1358. https://doi.org/10.2307/2408641.

Wickham, H. 2016. "ggplot2: Elegant Graphics for Data Analysis." In *Data Analysis*, edited by H. Wickham, 189–201. Springer International Publishing. https://doi.org/10.1007/978-3-319-24277-4_9.

Wickham, H., M. Averick, J. Bryan, et al. 2019. "Welcome to the Tidyverse." *Journal of Open Source Software* 4, no. 43: 1686. https://doi.org/10.21105/joss.01686.

Wray, G. A. 1995. "Evolution of Larvae and Development Modes." In *Ecology of Marine Invertebrate Larvae*, edited by L. McEdward, 413–448. CRC Press.

Wright, S. 1943. "Isolation by Distance." *Genetics* 28, no. 2: 114–138. https://doi.org/10.1093/genetics/28.2.114.

Yang, Z. 2015. "The BPP Program for Species Tree Estimation and Species Delimitation." *Current Zoology* 61, no. 5: 854–865. https://doi.org/10.1093/czoolo/61.5.854.

Zakas, C., and J. P. Wares. 2012. "Consequences of a Poecilogonous Life History for Genetic Structure in Coastal Populations of the Polychaete *Streblospio benedicti*." *Molecular Ecology* 21, no. 22: 5447–5460. https://doi.org/10.1111/mec.12040.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section. **Data S1:** mec70170-sup-0001-Supinfo1.docx.