Hybridization patterns between two marine snails, Littorina fabalis and L. obtusata

Abstract Characterizing the patterns of hybridization between closely related species is crucial to understand the role of gene flow in speciation. In particular, systems comprising multiple contacts between sister species offer an outstanding opportunity to investigate how reproductive isolation varies with environmental conditions, demography and geographic contexts of divergence. The flat periwinkles, Littorina obtusata and L. fabalis (Gastropoda), are two intertidal sister species with marked ecological differences compatible with late stages of speciation. Although hybridization between the two was previously suggested, its extent across the Atlantic shores of Europe remained largely unknown. Here, we combined genetic (microsatellites and mtDNA) and morphological data (shell and male genital morphology) from multiple populations of flat periwinkles in north‐western Iberia to assess the extent of current and past hybridization between L. obtusata and L. fabalis under two contrasting geographic settings of divergence (sympatry and allopatry). Hybridization signatures based on both mtDNA and microsatellites were stronger in sympatric sites, although evidence for recent extensive admixture was found in a single location. Misidentification of individuals into species based on shell morphology was higher in sympatric than in allopatric sites. However, despite hybridization, species distinctiveness based on this phenotypic trait together with male genital morphology remained relatively high. The observed variation in the extent of hybridization among locations provides a rare opportunity for future studies on the consequences of different levels of gene flow for reinforcement, thus informing about the mechanisms underlying the completion of speciation.


| INTRODUC TI ON
It is widely accepted that the number of traits contributing to reproductive isolation generally increases as speciation progresses (Seehausen et al., 2014;Smadja & Butlin, 2011). However, how traits under different evolutionary forces such as sexual selection (Ritchie, 2007), natural divergent selection (Nosil, 2012), and selection against maladaptive hybridization (Butlin, 1987;Hollander et al., 2018) interact with each other toward completing speciation is still largely unknown.
Distinguishing whether sister species are or are not completely reproductively isolated is a key step to identify traits involved in speciation. Traits differentiating species that are fully reproductively isolated could have evolved after speciation was complete and have not necessarily contributed to reduce gene flow between them (Butlin, 1987;Coyne & Orr, 2004;Nosil & Schluter, 2011). For example, a recent study on mangrove snails, Littoraria cingulata and L. filosa, found reproductive character displacement in assortative mating, but no gene flow between the two species, suggesting that displacement was caused by reproductive interference after speciation was complete (Hollander et al., 2018). By contrast, sister species that hybridize can be an important source of knowledge about the buildup of traits that act as barriers during the progress of speciation (Abbott et al., 2013).
Although divergence in the presence of gene flow is a widely accepted mechanism (Pinho & Hey, 2010), our view of speciation remains oversimplified. Even when models of divergence take into account multiple rates of gene flow across the genome and through time (Roux et al., 2014;Sousa, Carneiro, Ferrand, & Hey, 2013), inferences about the genomic architecture of speciation often assume that there is one landscape of divergence across the genome that is consistent in all contacts between a pair of species. However, divergence between two evolutionary units often involves multiple geographical replicates (e.g., sticklebacks, Jones et al., 2012;rough periwinkles, Butlin et al., 2014) or a wide distribution with multiple opportunities for hybridization in contact zones that are relatively independent from each other (e.g., house mouse, Smadja, Catalan, & Ganem, 2004; fire-bellied toads, Szymura & Barton, 1986). Thus, species interactions are likely to be context-dependent both in space and time, resulting in different rates of hybridization with locality-specific evolutionary consequences (Harrison & Larson, 2016). Such heterogeneity is particularly relevant to understand how different barrier traits accumulate and work in concert to strengthen reproductive isolation. The mechanisms completing speciation, such as reinforcement, depend not only on the opportunity for hybridization but also on its costs and benefits, which may vary idiosyncratically among contact zones. Studies of multiple replicates where sister species contact and have the opportunity to hybridize, but also comprising different environmental conditions, demographic scenarios, and geographic contexts of divergence, are thus needed to gain a more comprehensive understanding of how hybridization and barrier traits vary across the species' range.
These two sibling species present a largely overlapping distribution along the European Atlantic shores (Reid, 1996). However, at a local scale, pockets of allopatric populations can be found, especially in north-western Iberia (Sotelo et al., 2020). Moreover, L. obtusata is usually found in more sheltered habitats than L. fabalis. These include areas of Galician bays ("Rías") in Spain, as well as sites in northern Portugal that are somewhat protected from waves. Littorina fabalis, on the other hand, occupies more exposed areas with stronger wave action. In Iberian locations where the two species co-occur (hereafter referred to as sympatric for simplicity), they tend to show some level of vertical zonation, with L. obtusata occupying the mid to upper part of the shore, while L. fabalis tends to inhabit the lower part.
Three L. fabalis ecotypes were previously identified in this region, facing different wave exposure regimes and dwelling in different macroalgae/seagrass. The Mastocarpus Exposed (ME) ecotype is usually found in exposed sites on Mastocarpus spp.; the Zostera Sheltered (ZS) ecotype is found in a single sheltered region associated with Zostera spp.; and the Fucus Intermediate (FI) ecotype is commonly found in Fucus spp. (Carvalho et al., 2016;Rolán & Templado, 1987). Some variation in shell morphology associated with wave exposure has also been described within L. obtusata (Reid, 1996). However, in contrast to L. fabalis, there is no association between this variation and macroalgae species. Since the phenotypic differences and distribution of these variants in Iberia have not been characterized in a systematic manner, the ecotype terminology is generally not used for L. obtusata.
Contact between L. obtusata and all different L. fabalis ecotypes has been observed during fieldwork for this and previous studies (Carvalho et al., 2016) but, given the distribution of both species, those involving the FI ecotype are the most common. How this diversity in terms of geographic context of divergence and local environmental conditions influences the prevalence of hybridization between L. obtusata and L. fabalis remains unclear.
Flat periwinkles also exhibit high intraspecific shell polymorphism in color patterns, as well as in size and shape (Rolán-Alvarez et al., 2015). Although L. fabalis tends to be smaller and have a different shell shape (typically a rounder shell with a wider aperture) when compared to L. obtusata (typically a more elongated shell and smaller aperture relative to size), male genital morphology is the most reliable trait for distinguishing sister species (Reid, 1996). A comparative analysis across Littorininae revealed greater male genital shape divergence between sympatric/parapatric sister species when compared with allopatric pairs (Hollander, Smadja, Butlin, & Reid, 2013), with flat periwinkles standing out as a strong candidate for prezygotic isolation to have evolved as a consequence of hybridization. However, similar patterns could have resulted from reproductive interference to reduce direct costs associated with interspecific mating after reproductive isolation was complete (Hollander et al., 2018). Thus, the comparisons between populations with different levels of hybridization is a prerequisite for further tests of reinforcement or other processes leading to completion of reproductive isolation in this system.
Here we have analyzed multiple north-western Iberian populations representing different geographical contexts (allopatric and sympatric) between flat periwinkles. We analyzed genetic data (microsatellites and mtDNA) together with shell and male genital morphology of snails from 27 Iberian sites in order to (a) characterize the extent of hybridization between the two sister species; (b) evaluate differences in hybridization frequency and dynamics across distinct geographic settings; and (c) assess the influence of hybridization on the phenotypic differences (shell and male genitalia morphology) between species across sites.          (sympatric), as well as sites where only one species was found (locally allopatric) ( Figure 1, Table 1). For the analyses requiring reference allopatric populations, these were chosen based on several field surveys performed in these locations over the years, where we only found one species. Since fully grown individuals were needed for morphological analyses, sampling efforts were directed toward adults. Otherwise, individuals were randomly collected in terms of shell shape and size to avoid biasing our sampling toward either of the species or potential hybrids. The ecotype of L. fabalis individuals was recorded based on the algae/seagrass where they were found. Individuals were collected at the lowest tides (<0.75 m) and were brought alive to the laboratory where they were frozen at −20°C.

| Analysis of penis morphology
After carefully removing the soft tissue from the shell, individuals were inspected for the presence of male genitalia. Males (N = 818) were initially preclassified into species based on their penis "appearance" (i.e., visual classification). This consisted of comparing the length of the filament with respect to total penis length. A relative length of 10%-25% was considered typical of L. obtusata and of 30%-60% typical of L. fabalis. Individuals were classified as intermediate when the proportion was 25%-30% or as unknown when they differed substantially from the typical proportions for the two species (Carvalho et al., 2016;Reid, 1996). The penis was then dissected at the base of its insertion and preserved separately from the body, both in 96% ethanol.
A discriminant function analysis (DFA) based on linear measurements of the penis was performed to evaluate the accuracy of this visual classification, using the mass r package v7.3.49 (Venables & Ripley, 2002). For those individuals where we were able to retain an intact penis (N = 278), seven morphological features were measured (Appendix S1) under a stereomicroscope (

| Geometric morphometric analyses of the shell
To evaluate variation in shell morphology between L. fabalis and L. obtusata, identify putative hybrids, and examine the effect of geographical context (i.e., allopatry vs. sympatry) on phenotypic variation, the shells Landmark coordinates were superimposed using Generalized Procrustes Analysis (Rohlf & Slice, 1990), to standardize the scale, location, and rotation, and to optimize the position of semilandmarks by minimizing bending energy (Mitteroecker & Gunz, 2009). and Morás (L. obtusata), because most sampled individuals were juveniles. Nevertheless, samples from these five locations were included in the genetic analyses.

| Shell variation across species and ecotypes
To characterize variation in shell morphology across species and ecotypes of L. fabalis without the influence of ongoing hybridization, we first conducted a set of analyses on individuals from allopatric populations (L. obtusata, N = 108; L. fabalis, N = 174; Figure 1, Table 1). We used a principal component analysis (PCA) of shape variables to explore general patterns of shape variation. To test whether L. obtusata and the three ecotypes of L. fabalis differed in shell size and shape, we examined general linear models (GLM) for logCS and the multivariate set of shape variables in two sequential analyses: the first with species (grouping all L. fabalis ecotypes) as the main factor and location of collection as a factor nested within species; and the second within L.
fabalis, with ecotype as the main factor and location also as a nested factor. To account for possible allometric effects of size on shape, we also examined a GLM for shape using size as a covariate, with the same main factors as before, and location as a nested factor. The significance of different terms was evaluated using residual randomization in permutation procedures based on Procrustes distances and consisting of 1,000 permutations, as implemented in the function procD.lm of geomorph R package (Adams, Collyer, & Kaliontzopoulou, 2018), and Zscores were used for significance testing, as recommended by Adams and Collyer (2016). Shape differences between groups (species and ecotypes of L. fabalis) were visualized using deformation grids.

| Hybridization effects on shell morphology
In order to identify shell morphological features that differ the most between the two species and examine whether genetic hybrids exhibited intermediate shell morphology, we performed a DFA on shape using the r package mass v7.3.50. The DFA was constructed using individuals from reference allopatric populations of both species (N = 282), including all L. fabalis ecotypes, aiming to capture as much morphological variation as possible. The Pp of assignment of each individual to each species based on shell shape was estimated using a leave-one-out cross-validation procedure. Based on this discriminant function, we then inferred the morphological Pp of each individual from the remaining populations (from sympatric sites) and for which a genetic assignment was available (N = 339), in order to compare the morphological predicted shapes with the genetic membership coefficients obtained with structure (global analysis, total N = 556; see below).

| Species differences in shell morphology across different geographic contexts
To examine how the contact between species may influence shell shape and size irrespective of hybridization, we compared the overall differences between the two species in distinct geographic contexts (allopatry vs. sympatry). As we were interested in examining morphological variation of the genetically pure individuals of each species, samples for which no genetic information was available and those genetically classified as hybrids in the global structure analysis were removed from the dataset. This resulted in a new subset F I G U R E 2 Geometric morphometric (GM) analysis of the shell. (a) Landmarks (LM) digitized in specimens of Littorina obtusata (left) and Littorina fabalis (FI ecotype; right), placed in the standard position. LM1-LM4 represent fixed LMs, whereas all the other points represent semilandmarks. (b-e) Results for the characterization of shape and size for the allopatric reference populations. (b) Mean log-centroid size (CS) ordered by species and ecotypes (vertical bars denote 95% confidence intervals). Population codes are described in Table 1. (c) Plot of the two first principal components (PC1 and PC2) of shape variation and deformation grids at maximum and minimum PC1 values compared to the global mean. (d, e) Deformation grids depicting the mean shapes for each species and ecotype, respectively. Mean shapes were magnified 2× to enhance visualization of 540 individuals: 117 L. fabalis and 100 L. obtusata from allopatric populations and 183 L. fabalis and 140 L. obtusata from sympatric populations. A GLM analysis was implemented to examine how shell morphology responds to species sympatry. We sequentially examined shell size, shape, and shape while taking size variation into account and fitted GLMs that included geographical context and species, as well as their interaction. As in previous analyses, we also included sampling location as a nested factor to account for local variation among populations.

| Molecular methods
The genetic characterization of the individuals collected during this study was based on two types of markers: microsatellites and mtDNA. While microsatellites were shown to distinguish the two species and identify recent-generation hybrids (Carvalho et al., 2015(Carvalho et al., , 2016, mtDNA demonstrated introgressive hybridization between flat periwinkles (Kemppainen et al., 2009;Sotelo et al., 2020).

| DNA extraction
Genomic DNA was extracted from head-foot tissue using a modified version of the standard high-salt protocol (Sambrook, Maniatis, & Fritsch, 1989), by replacing the lysis buffer by cetyltrimethyl ammonium bromide (Winnepenninckx, Backeljau, & Wachter, 1993). The final DNA concentration was standardized across samples (~5 to 10 ng/μl).  tetra-nucleotides repeat motifs) were selected based on their informativeness for species discrimination, as well as on their genotyping reliability after combining them into two multiplexes (Appendix S2).

| Microsatellite genotyping
Amplification reactions were performed as in Carvalho et al. (2016) and run on a ABI 3730 sequencer (Applied Biosystems) at STABVIDA.

| Mitochondrial DNA amplification and sequencing
A fragment of the mitochondrial gene cytochrome b (Cyt-b) was amplified using the cytbF-cytbR primer pair (Panova et al., 2011).

TA B L E 3
General linear model analysis (GLMs) of morphological differences in terms of size and shape, as well as shape accounting for the influence of size (CS) results, newhybrids analyses were only performed for the global dataset without a priori information on allele frequencies or admixture proportions. We used a "Jeffrey's-like" prior, which considers that some alleles may be rare or absent in the different populations and so more accurately determines the assignment of hybrid individuals to their respective categories (Anderson, 2003). To assess the consistency of estimates, three replicates of 1,000,000 MCMC iterations, after a burn-in of 100,000, were performed. A procedure to estimate the threshold of Pp (TPp) for the classification of individuals as pure or hybrid, similar to the one applied for TQ (structure), was implemented (Appendix S2).

| Classification of males based on genital morphology
The males analyzed both in terms of visual appearance and using the DFA (based on seven morphological features, with a Pp > .99) were classified into species with high concordance (94.3%, excluding 68 individuals from allopatric sites used as reference; Table 2).
The differences correspond to individuals classified as intermediate using one approach and as pure of one species using the other approach and with a single exception were all observed in sympatric sites. The discriminant axis contrasted the length of the filament to the length of the gland row, the number of rows, and the number of glands (Appendix S1). Since the number of samples visually classified was much higher than those available for the DFA, for which only intact penis could be used, the former was used in subsequent analyses.

| Variation in shell morphology
The GM analysis provided 52 shape variables plus CS. The two species differed significantly (Z = 1.511, p = .009) in terms of size, with L. obtusata being larger than L. fabalis (Figure 2b, Table 3), despite significant variation among locations (Z = 6.460, p = .001). Within L. fabalis, significant variation in CS was observed among sampling locations (Z = 3.736, p = .001), but no significant differences were observed among ecotypes (Table 3). The PCA of shape variables retrieved two components that cumulatively explained over 75% of total shape variation ( Figure 2c). The two sister species occupy different parts of this morphospace despite some overlap, while the distinction among L. fabalis ecotypes was less evident, as expected.
Accordingly, GLMs revealed significant shape differences both Note: N GM, total number of individuals analyzed; GM ≥ .90, number of individuals with a posterior probability equal to or higher than .90 to either species; Intermediate shape, number of individuals not classified to species (posterior probability below .90 to both of them). a Allopatric sites were analyzed separately from the sympatric ones (see Section 2).

TA B L E 4 Classification of individuals
for each location based on a geometric morphometric (GM) analysis of the shell between species and among ecotypes of L. fabalis, despite significant local variation (Table 3). Differences in shape among locations and ecotypes of L. fabalis were also significant irrespective of size variation (Z = 2.039, p = .016; and Z = 2.604, p = .006, respectively), but this was not the case between species (Table 3). Indeed, location was the factor with the strongest effect on shape variation, fol-

| Morphological variation in allopatric versus sympatric populations
The DFA of allopatric populations rendered a high percentage  Table 4).
After excluding hybrids, and despite significant local variation as represented by sampling location, species and geographic context (i.e., allopatric vs. sympatric) interacted significantly in their effect on size (Z = 1.572, p = .005), but not on shape (Figure 4, Table 5). Species and geographic context also did not interact significantly when taking size variation into account (Table 5). The examination of CS variation across locations indicated that the two species tend to be more similar in size in sympatry, as L. obtusata becomes slightly smaller and L. fabalis slightly larger than in allopatry ( Figure 4).

| Classification of individuals using genetics and assessment of hybridization
The initial structure analysis revealed a clear distinction between the two species ( Figure 5

L. obtusata -Sympatric
Our simulations show that under the current settings both structure and newhybrids erroneously assigned more hybrid genotypes as pure than vice versa, suggesting that our estimates of the number of hybrids are conservative.
Concerning the real dataset, 1,012 from the 1,059 analyzed individuals were unambiguously assigned to genotype classes, whereas the remaining individuals (4% of the total) could not be assigned with confidence to any of the six classes (Pp ≥ .80). Among the confidently assigned individuals, only 32 were classified as hybrids by newhybrids, distributed over three locations ( Two main clades of mtDNA haplotypes were identified (Figure 6).
From a total of 762 sequenced individuals that were also analyzed for microsatellites, 423 carried haplotypes from clade I (typical of L. fabalis) and 339 presented haplotypes from clade II (typical of L. obtusata) ( Table 7). The total proportion of L. obtusata-typical haplotypes present in L. fabalis individuals was slightly higher than the other way around (18.5% and 16.2%, respectively). However, the proportion of atypical haplotypes was very heterogeneous among sites (   Table 8).

| Comparison of classification based on genetics, male genitalia, and shell morphology
None of the individuals assigned to one species by the structure analysis was classified into the other species by either approach based on the genital morphology (N = 565, 337 L. fabalis and 228 L. obtusata). The converse was also true: genital classification to parental groups was 100% congruent with structure-based classifica-

| D ISCUSS I ON
Characterizing the patterns of hybridization between closely related species pairs is an important step toward understanding the role of gene flow in speciation (Abbott et al., 2013). Introgressive hybridization between flat periwinkles was previously suggested based on mtDNA analyses (Kemppainen et al., 2009), together with the detection of hybrids between the two species in Cabo do Mundo using microsatellites (Carvalho et al., 2016). However,

| Morphological and genetic differences and hybrid identification
The concordance found between genetic data and adult males' genital shape confirms that the latter is a reliable trait to distinguish the two species. On the other hand, the landmark-based GM analysis of the shell was less reliable, especially when comparing individuals from sympatric sites. Furthermore, neither shell shape nor male genitalia traits allowed for an accurate identification of hybrids between species.
The microsatellite panel was powerful not only to classify individuals correctly into species (>97.5% of the simulated genotypes to each parental group) but also to identify hybrids (from 85% to 97.8% using newhybrids and structure, respectively). Although there was a lower percentage of hybrid genotypes correctly assigned to a specific hybrid class (78%), this was based on 11 loci. The use of additional markers will certainly increase the statistical power to distinguish among different hybrid categories.
The two main mtDNA clades support two divergent lineages that largely correspond to the two species of flat periwinkles. However, shared haplotypes between the two species were observed in substantial proportions, as previously shown by Kemppainen et al. (2009). Although shared haplotypes between species could result F I G U R E 5 Membership coefficient (Q ranging from 0 to 1) of all genotyped individuals to the clusters identified by the initial structure analysis for k = 2. Each vertical bar represents an individual, where typical Littorina fabalis ancestry is represented in yellow and typical Littorina obtusata ancestry in blue. Vertical bars with both yellow and blue represent individuals with admixed ancestry. On top, location names follow Table 1 and "*" indicates allopatric locations. Results are consistent across replicates

| Extent and patterns of hybridization
The number of locations where hybrids were detected using microsatellites varied between three and seven (depending on the approach), with the majority of hybrids found to be BCO, followed by BCF and F2s. Although no signatures of hybridization were found in TA B L E 6 Number of hybrids between Littorina fabalis and Littorina obtusata detected across locations using structure (left) and newhybrids (right) The proportion of sites with signatures of mtDNA introgression is much higher than those with evidence for early generation hybrids using microsatellites. This suggests that hybridization could have been more frequent in the past, prior to the evolution of multiple reproductive barriers, including in currently allopatric sites. However, because mtDNA is unlinked to selected nuclear loci, we cannot exclude that mtDNA introgression can persist and spread over long into L. obtusata (Sotelo et al., 2020). Asymmetric mtDNA introgression into L. fabalis is in line with previous work showing that males of both species prefer larger females (Saltin, 2013), which suggests that most successful interspecific crosses in Northern European sites preferentially involve a L. obtusata female and a L. fabalis male. However, different densities of one species relative to the other (or of both species) could also interfere with the choosiness in the field and partially explain differences in the direction of introgression among sites (Carvalho et al., 2016).

| Comparison between allopatric and sympatric sites
The number of hybrids and introgressed individuals was clearly higher in sympatric than allopatric locations. This is true whether we conditions in the sympatric geographic context. However, these results need to be interpreted with caution. Differences between locations, the main factor contributing to shell size variation in this system, could at least partially be explained by individuals' age. Even though only adult individuals were analyzed, size is known to increase with age. Thus, we cannot exclude that some of the observed patterns in terms of size are heavily influenced by the age at which individuals were collected, which could not be measured.

| Hybridization implications for flat periwinkles' diversification
The relatively high variation among sites in terms of hybridization levels supports the suggestion that the different contacts between flat periwinkle species at local geographical scales can have different evolutionary outcomes. These outcomes depend on factors such as environmental conditions, strength of selection, population densities, and genetic background, among others. Consequently, the degree of reproductive isolation, as well as the reproductive barriers that have accumulated between L. fabalis and L. obtusata, may well differ among locations but the causes will be hard to disentangle.
Distinct hybridization outcomes among replicates can result from genetic (Borge et al., 2005), behavioral (Bettles et al., 2005), and/ or environmental differences (Muhlfeld et al., 2014;Renaut et al., 2012). Environmental change, in particular, has been suggested as the cause of temporal changes in hybridization rates between two trout species, with a recent increase related with changes in precipitation and water temperature (Muhlfeld et al., 2014). At least some of the populations analyzed in this study face marked environmental fluctuations, where exposure to strong wave action may lead to substantial demographic oscillations (Carvalho et al., 2016;Reid, 1996), contributing to temporal instability in potential for hybridization.
Furthermore, the recent decline of fucoid macroalgae inhabited by flat periwinkles at the southern limit of their distribution (Nicastro et al., 2013), including the Iberian Peninsula, could contribute to the increase of hybridization in this region.
The fact that hybrids were relatively rare in most of the sites analyzed here suggests that substantial reproductive barriers exist between flat periwinkle species, at least in most sites where they coexist. The difference in penis morphology between flat periwinkles makes this trait an important candidate for playing a role in prezygotic reproductive isolation between these species (Hollander et al., 2013;Reid, 1996). Furthermore, understanding the causes for the high number of hybrids observed in one site and considering whether reinforcement could have evolved in some of these sites are two important issues that need to be addressed in future studies.
The possibility of capitalizing on multiple local contacts between L.
fabalis and L. obtusata, with idiosyncratic dynamics, opens exciting research avenues in this system. Performing mating experiments between individuals from sites with varying rates of gene flow Abbreviations: N, total number of hybrids according to each software sequenced for mtDNA; Clade I, number of hybrids with mtDNA haplotype from clade I (typical of L. fabalis); Clade II, number of hybrids with mtDNA haplotype from clade II (typical of L. obtusata). a It was not possible to clearly distinguish between ME and FI.
TA B L E 8 Mitochondrial DNA clades of the hybrids identified with structure in each location and Littorina fabalis ecotype involved, and the same relatively to the hybrids identified with NEWHYRBIDS for each class between species, phenotypic analyses of traits involved in mating, and the characterization of environmental conditions and species densities among sites are important steps forward to understand the late stages of speciation in this system.
Although the process of divergence between two lineages (forms, ecotypes, or species) tends to be oversimplified as if there was a single possible outcome, the availability of multiple natural replicates suggests that different evolutionary routes are possible.
Thus, systems like the one presented here are crucial to understand the different paths that can be taken in the late stages of speciation and to determine which paths are more likely to lead to complete speciation.

ACK N OWLED G M ENTS
We would like to thank David Reid for initial discussions and training on analyses of shell and genitalia morphology. This study was supported by European Regional Development Fund (FCOMP-01-0124- Note: Shown is: N GEN and GM data, the number of individuals analyzed for both genetics and shell morphology; GEN and GM L. fabalis N (%), number (and percentage) of concordant assignment using genetics and shell morphology to L. fabalis; GEN and GM L. obtusata N (%), number (and percentage) of concordant assignment using genetics and shell morphology to L. obtusata; Intermediate shape N (%), number (and percentage) of individuals genetically pure with intermediate shell shape; GEN and GM Mismatch N (%), number (and percentage) of individuals genetically classified into one species with shell shape typical from the other species; GEN Hybrids, number of individuals genetically hybrid classified as pure from each species based on shell morphology. a For the GM analysis, allopatric sites were analyzed separately from the sympatric ones (see Section 2).