Genetic and morphological divergence between Littorina fabalis ecotypes in Northern Europe

Low dispersal marine intertidal species facing strong divergent selective pressures associated with steep environmental gradients have a great potential to inform us about local adaptation and reproductive isolation. Among these, gastropods of the genus Littorina offer a unique system to study parallel phenotypic divergence resulting from adaptation to different habitats related with wave exposure. In this study, we focused on two Littorina fabalis ecotypes from Northern European shores and compared patterns of habitat‐related phenotypic and genetic divergence across three different geographic levels (local, regional and global). Geometric morphometric analyses revealed that individuals from habitats moderately exposed to waves usually present a larger shell size with a wider aperture than those from sheltered habitats. The phenotypic clustering of L. fabalis by habitat across most locations (mainly in terms of shell size) support an important role of ecology in morphological divergence. A genome scan based on amplified fragment length polymorphisms (AFLPs) revealed a heterogeneous pattern of differentiation across the genome between populations from the two different habitats, suggesting ecotype divergence in the presence of gene flow. The contrasting patterns of genetic structure between nonoutlier and outlier loci, and the decreased sharing of outlier loci with geographic distance among locations are compatible with parallel evolution of phenotypic divergence, with an important contribution of gene flow and/or ancestral variation. In the future, model‐based inference studies based on sequence data across the entire genome will help unravelling these evolutionary hypotheses, improving our knowledge about adaptation and its influence on diversification within the marine realm.


| INTRODUC TI ON
The marine rocky intertidal represents one of the most abrupt environmental gradients on Earth (Little & Kitching, 1996;Raffaelli & Hawkins, 1996;Tomanek & Helmuth, 2002). Different environmental conditions across the tidal range result in patterns of vertical and horizontal (along the shore) zonation both in terms of species diversity and in terms of intraspecific phenotypic variation (Connell, 1972;Little & Kitching, 1996;Raffaelli & Hawkins, 1996;Southward, 1957). Wave action is one of the most important physical selective agents across intertidal environments worldwide, shaping both axes of zonation (Le Pennec et al., 2017). Habitats exposed to, or sheltered from waves differ consistently in their biodiversity (Denny & Wethey, 2001;Helmuth & Denny, 2003). Biotic factors are also known to affect the intertidal community (e.g. presence of predators; Paine et al., 1994), which adds to the abiotic selective pressures (Seeley, 1986). Altogether, these environmental gradients make the intertidal a natural laboratory to study local adaptation and ecological speciation.
Taxa with low dispersal, where divergent selection can be strong enough to counteract gene flow among populations and promote divergence with gene flow, are particularly well suited for studies about local adaptation and ecological speciation (Sanford & Kelly, 2011;Smadja & Butlin, 2011). Among these, species with distinctive phenotypes associated with different microhabitats (i.e. ecotypes) in the intertidal can provide important information about how natural selection influences biological diversification (e.g. Coyer et al., 2011;Kess et al., 2018;Sanford & Kelly, 2011;Wilding et al., 2001). Instances of parallel evolution of ecotypes across similar environmental gradients in multiple locations across a species' distribution, that is parallel evolution, are viewed as support for a role of natural selection in driving divergence. This is because it is unlikely that ecotypes have repeatedly evolved in the same phenotypic direction if only purely stochastic processes were involved (Johannesson, 2001;Nosil, 2012;Schluter, 2000).
The characterization of the genetic variation underlying parallel evolution allows to distinguish if repeated events of phenotypic divergence tend to involve the same de novo mutations, different de novo mutations in the same or different genomic regions, ancestrally shared standing polymorphisms, the same alleles due to migration between populations or a combination of all the above (Elmer & Meyer, 2011;Faria et al., 2014;Johannesson et al., 2010;Nosil, 2012). Thus, the study of parallel phenotypic divergence across intertidal microhabitats can help us understand if adaptation usually involves the same genetic paths, as well as the relative contributions of ancestral polymorphism and/or gene flow.
Parallel ecotypic divergence in the intertidal zone has been particularly well documented in the rough periwinkle Littorina saxatilis (Gastropoda; Reid, 1996;Rolán-Alvarez et al., 2015). Two main ecotypes have been described: a large, thickshelled ecotype inhabiting sheltered microhabitats that faces intense crab predation (Crab ecotype); and a small, thin-shelled one facing heavy surf in exposed microhabitats (Wave ecotype) (reviewed in Johannesson et al., 2010). These ecotypes, found only tens of metres apart in Spain, the United Kingdom (UK) and Sweden, have diverged in parallel within each of these countries .
Amplified fragment length polymorphisms (AFLPs) allowed the initial identification of loci under divergent selection (hereafter "outliers") between ecotypes and provided the first insights on the proportion of outliers shared at different geographic scales, within and among countries (Galindo et al., 2009(Galindo et al., , 2013Grahame et al., 2006;Hollander et al., 2015;Wilding et al., 2001). Benefiting from the assembly of a reference genome and the construction of a genetic map for L. saxatilis, the heterogeneous genomic differentiation between ecotypes was recently confirmed with the identification of some genomic regions enriched for the presence of outliers, which tend to coincide with polymorphic inversions Westram et al., 2018). Moreover, whole-genome sequencing of pools of individuals from multiple populations of both ecotypes across the species' distribution range revealed that outlier sharing tends to be high even among distant populations, although it decreases with the geographic distance between populations .
A closely related species for which phenotypic variation is also found associated with an environmental cline in wave exposure is Littorina fabalis (Kemppainen et al., 2005;Reimchen, 1981;Tatarenkov & Johannesson, 1994, 1999. Individuals with large and thick shells (hereafter "large ecotype") are commonly found in moderately exposed shores, whereas individuals with smaller and thinner shells (hereafter "dwarf ecotype") predominate in sheltered habitats. Contrary to L. saxatilis that live on rocks, these two L. fabalis ecotypes dwell on brown macroalgae (Fucus spp. and Ascophyllum spp.), grazing on the epiphytes that grow on the algae fronds (Williams, 1990). The fronds are also thought to provide refuge against one of their main predators, the green crabs future, model-based inference studies based on sequence data across the entire genome will help unravelling these evolutionary hypotheses, improving our knowledge about adaptation and its influence on diversification within the marine realm.
Habitat-related variation in one allozyme locus (arginine kinase, Ark) was initially found in Swedish populations of L. fabalis, suggesting that this locus was under the influence of natural selection related to wave exposure and/or other associated factors (Tatarenkov & Johannesson, 1994). The differences in Ark allele frequencies between sheltered and moderately exposed habitats were also associated with variation at a random amplification of polymorphic DNA (RAPD) locus and with the size differences described above. This is true even for sites with intermediate exposure, suggesting some reproductive isolation between the ecotypes despite gene flow (Johannesson & Mikhailova, 2004;Tatarenkov & Johannesson, 1998).
Similar habitat-related phenotypic divergence has also been observed at least in the UK (Wales), France and Norway, where contrary to Sweden, L. fabalis also has to face high tidal amplitudes (Kemppainen et al., 2009(Kemppainen et al., , 2011Reimchen, 1981;Tatarenkov & Johannesson, 1999). Even if a northern refugium could have existed during the last glacial maximum (LGM) for L. fabalis , most of these shores were likely colonized after the LGM (Charbit et al., 2007). This, together with the absence of significant genetic differentiation between the ecotypes across this region using neutral markers , suggests a relatively recent (after the LGM) local establishment of habitat-related phenotypic divergence.
Contrary to neutral markers, Ark intron sequencing revealed highly significant divergence between individuals from sheltered and moderately exposed habitats (Kemppainen et al., 2011). Results show that one haplotype was almost fixed and shared across sheltered habitats of these different countries whereas wave-exposed habitats maintained similarly higher variation. This increases the possibility of "evolution in concert", where locally adapted alleles could arise once and subsequently spread to geographically distant populations inhabiting the same habitat by means of ecotype-specific selective sweeps (Johannesson et al., 2010;Kemppainen et al., 2011;Schluter & Conte, 2009). However, except for Ark and the putatively linked RAPD locus, whether the same genetic variation is involved in the evolution of these ecotypes is currently unknown.
Thus, studies that integrate both morphological data and a high number of nuclear markers from different locations are needed to evaluate the parallelism of phenotypic divergence in L. fabalis across its geographic range, as well as the genetic variation and processes involved.
In this study, we used shell geometric morphometrics and AFLPs to perform the first characterization of L. fabalis ecotypes across multiple Northern European populations (Norway, Sweden and the UK) with four main goals: (a) to understand if shell shape and size divergence evolved in the same direction among locations within the same country and across countries; (b) to identify loci involved in the differentiation between sheltered and exposed sites (i.e. outlier loci); (c) to quantify the degree of outlier sharing at two different geographic scales: within and among countries; and (d) to contrast population structure and relatedness using outlier (putatively adaptive) versus nonoutlier loci (putatively neutral) in order to gain insights about the evolutionary history of phenotypic divergence.
Although not fully conclusive, our results are compatible with phenotypic parallel divergence within L. fabalis associated with habitat and reveal a pattern of outlier sharing that suggests a relevant role of gene flow and/or retention of ancestral polymorphism in the evolution of ecotypes.

| Sampling
Littorina fabalis individuals from both sheltered and moderately exposed habitats were randomly collected with respect to their shell morphology between August and October 2012, the sample size in each population varied between 21 and 72 (average N = 39) (Table S1). A nested sampling design was implemented: samples were collected from each habitat within each location with replicates at two geographical scales: regional (two to three locations within each country, <50 km) and global (different countries, >1,000 km) ( Figure 1, Figure S1, Table S1). This allows the investigation of parallel evolution at these two contrasting scales.
Within each location, sites representing each type of habitat were preselected based on information from previous studies (Kemppainen et al., 2009;Tatarenkov & Johannesson, 1994), or on their orientation and topography as retrieved from the Google Earth Engine (Gorelick et al., 2017). This preclassification was further confirmed in situ, where sheltered habitats were distinguished from moderately exposed ones by the high abundance of Ascophyllum spp., a macroalga that is commonly used as an indicator of sheltered shores (Tatarenkov & Johannesson, 1998).
However, two putatively exposed sites did not adjust to these patterns: those in Anglesey South (UK) and in Ursholmen (Sweden) (Table S1). In Anglesey South, Ascophyllum spp. was only abundant in the upper part of the shore but individuals were collected from the lower part, where wave action is likely rather high. In Ursholmen, the moderately exposed site was chosen based on its orientation towards relatively open sea but Ascophyllum spp.
was highly abundant there. The site was still included in the study but re-classified as sheltered (Table S1). Since we were not able to sample a proper moderately exposed site in Ursholmen, this location was excluded from statistical tests on shell morphology. It is also important to emphasize that the distance between habitats within each location differs among countries. In the Scandinavian locations, snails are continuously distributed along the shore from sheltered to moderately exposed sites (i.e. parapatric), which are <500 m apart. In the UK, the two habitats within each sampled location are geographically isolated from each other, from 400 m to 8 km apart (Figure 1, Figure S1).
After collection, individuals were brought alive to the laboratory where they were frozen at −20°C and then stored in 95% ethanol.

| Sample processing and classification into species
Shells of adult individuals were photographed over graph paper (used for scale) in a standardized position following Carvajal-Rodríguez et al. (2005) and using a digital ICA video module fitted on a MZ12 dissection microscope (Leica) for subsequent geometric morphometric analyses. Shells were then broken, and individuals dissected and sexed under the same microscope. Since shell morphology is not completely diagnostic between the two closely related species of flat periwinkles (L. fabalis and L. obtusata), we could not guarantee that all collected individuals were indeed L. fabalis just based on their shells. Thus, we focused our analysis on males, which were classified into L. fabalis or L. obtusata based on the morphology of their genitalia, one of the most distinctive traits between the two species Reid, 1996).
In sites where the number of males was too low, females and immature individuals were also included. Nevertheless, all samples (including males) were later classified into species using the AFLP genotypes. To do so, 17-19 individuals from L. obtusata populations within each country were deliberately included in the genetic analysis as references.

| Shell geometric morphometrics
Geometric morphometric (GM) analyses were carried out following the methodology previously developed for flat periwinkles by Costa et al. (2020). This consisted in digitizing a total of 28 landmarks for each shell of adult individuals (males and females classified as L. fabalis), including 4 fixed and 24 sliding semilandmarks, using TPSDIG v1.40 (Rohlf, 2006).
A principal component analysis (PCA) of shape variables was conducted to assess overall patterns of variation. In addition, UPGMA dendrograms for size and shape were generated to evaluate clustering patterns between samples from each habitat and location, based on Euclidean distances (for size) and Procrustes distances (for shape) of population means. General linear models (GLMs) were then used to assess if shell shape and size (logCS) differed significantly between the two sampled habitats and to evaluate whether these differences varied among locations and countries. As such, GLMs included country, habitat and sampling location (nested within country) as main effects, as well as all interaction terms. Allometric effects were also investigated by performing a GLM on shape with logCS as a covariate, using the same factorial design. These analyses were carried out in the geomorph R package , using randomized residual permutation procedures (RRPP) of 1,000 permutations and Z-scores for significance assessment (Collyer & Adams, 2018, 2019. Deformation grids were used to visualize differences in shape between L. fabalis individuals from the two habitats within each location. Finally, a discriminant function analysis (DFA) based on shape was implemented using the R package MASS v7.3.50 to infer the probability of morphological assignment of individuals into moderately exposed or sheltered habitats. The DFA was constructed based on the 287 individuals from all locations under study except for Ursholmen, using a leave-one-out cross-validation procedure.
The resulting morphological posterior probability (PP) assignments were then compared with the genetic membership coefficients obtained in STRUCTURE (see below) for the 83 individuals for which both genetic and phenotypic data were available. All morphological data analyses where carried out in the R language for statistical computing (R Development Team, 2019).

| Genetic analyses
Genomic DNA was extracted using the CTAB-chloroform protocol described in Galindo et al. (2009) com/genev alab/AFLPT ools) that follows AFLPScore methodology (Whitlock et al., 2008). The final genotypes were obtained using this same package.

| Detection of outlier loci
Two different methodologies were used to identify outlier loci between L. fabalis exposed and sheltered habitats, BAYESCAN v.2.1 (Foll & Gaggiotti, 2008) and DFDIST (Beaumont & Nichols, 1996). In both cases, the analyses were carried out independently within each location. As previously mentioned, both Ursholmen sites were likely sheltered. Thus, we do not expect to detect outliers related to the level of wave-exposure or associated factors. Nevertheless, outliers between the two sites were still estimated as a control.
BAYESCAN calculates population-specific and locus-specific F ST coefficients and then estimates the posterior probabilities of two alternative models (including or excluding the effect of selection) for each locus using a reversible-jump Markov chain Monte Carlo (MCMC) approach. Ten pilot runs (10,000 iterations) were performed to tune the model parameters, followed by 400,000 iterations (100,000 as burn-in, 20 as thinning interval and 20,000 as sample size). Loci were identified as outliers when the posterior probability was higher than 0.75 (equivalent to a Bayes factor of 3 or greater) (see Foll & Gaggiotti, 2008

| Genetics substructure analysis
Genetic substructure analyses were performed with these two different sets of loci, all outliers (i.e. putatively under divergent selection) and nonoutliers (i.e. putatively neutral). In both cases, they represent the combination of loci detected in each locality. AFLP-SURV v.1.0 (Vekemans et al., 2002) was used to calculate population pairwise genetic differentiation (F ST ) and Nei's genetic distance (10,000 bootstrap) using Zhivotovsky's (1999) Bayesian approach.
Neighbour-joining (NJ) trees were constructed based on Nei's genetic distance using the NEIGHBOR routine implemented in the PHYLIP package (Felsenstein, 2013 regional, including all sites within a country; and (c) local, including both the sheltered and exposed site within a location. All analyses were performed with the outlier and nonoutlier data sets, considering only as outliers those detected within the respective hierarchical levels (i.e. for the global analysis using all detected outliers in all pairwise comparisons, for the regional analysis only those detected in the pairwise comparisons within the corresponding region/country, whereas for the local analysis using only the outliers detected in that specific location). The consistency of individuals' assignment based on outliers across the different hierarchical levels was assessed for each location using correlation tests (Pearson correlation coefficient (r)- Sokal & Rohlf, 1995), implemented in the R package Stats v.3.6.1.
For each K-value, we analysed five replicate runs of 500,000 iterations (100,000 as burn-in), assuming an admixture model, correlated allele frequencies and without prior population information. All analyses were carried out from K = 1 up to K = number of sites plus one, depending on the hierarchical level analysed. The method developed by Evanno et al. (2005) implemented in STRUCTURE HARVESTER was employed to determine the best K.

| Characterization of shell morphology and phenotypic divergence
Adult individuals from each habitat differed significantly in size, with individuals from moderately exposed habitats being consistently larger than those from sheltered habitats (Figure 2a, Table S2).
The same factors influencing size also influence shape variation, as reflected by significant effects of habitat, location and their interaction, but not of country (Table S2). When considering allometric effects, we found a significant influence of size on overall shape.
In addition, size seems to affect shape differently across locations (significant interaction between both effects) (Table S2). Taking this effect of size into account, the same effects of habitat, location and their interaction remained significant, suggesting that although size accounts for some of the observed shape variation, it is not the only factor influencing it.
The distribution of individuals in the morphospace of the two principal components of shape variation, which explained a total of 69.4% of variation (PC1: 42.1%; PC2: 27.3%), revealed that average shapes of individuals from the moderately exposed habitats tend to cluster in the lower-left quadrant, whereas the average shape of individuals from the sheltered habitats is found in the upper-right part, with the exception of Anglesey South (Figure 3a). This agrees with the GLM results that show habitat as the factor with the strongest effect (based on Z-scores) on shape variation (as well as on size), followed by location (Table S2). The UPGMA dendrogram based on shape revealed two main groups, one consisting on the moderately exposed habitat from Scandinavian locations except Syltøyna (as for size); and the second comprising the sheltered habitats from Scandinavian locations and Syltøyna as well (Figure 3b). Contrary to what was observed in terms of size, although populations within the UK clustered by habitat rather than by location, they were nested within the clade composed mostly by sheltered sites from other countries, suggesting also important geographic influence on shape, in accordance with GLM results ( Table 2).
The inspection of deformation grids depicting the mean shape for individuals from each habitat (Figure 3c) shows that individuals inhabiting the moderately exposed habitat generally exhibit a larger width of the aperture (relative to overall shell size), whereas individuals of the sheltered habitat tend to exhibit a smaller aperture.

| Detection of outlier loci and comparison among locations
A total of 746 AFLP loci were analysed in 299 individuals (average sample size per population of 21) ( Excluding Ursholmen, a total of 1.6%-5.6% of outliers (p > .99) and 6.3%-9.5% of outliers (p > .95) were detected across all locations using DFDIST (Table 1, Figure S2). Using the most conservative threshold (p > .99), the number of outliers per location ranged from 6 in Lökholmen (Sweden) to 29 in Hummelsund (Norway) with the F ST between habitats for these loci ranging from 0.0895 to 0.2922, respectively. This contrasts with lack of differentiation for nonoutlier loci (Table 1). In total, only three outliers (3.8%) are shared among all countries whereas the highest number of outliers are shared between locations within each country (Figure 4), with Norway presenting the highest values (between 8 and 14 outliers, 21% and 36.8%, respectively) followed by the UK (7 outliers, 20%) ( Table 2).
The number and proportion of shared outliers between locations of different countries were generally lower, ranging from 2 (6.9%) between Lökholmen (Sweden) and Syltøyna (Norway) to 7 (14.6%) between Anglesey North (UK) and Hummelsund (Norway; Table 2).
In total, 147 outliers (p > .95) were detected across locations; this is combining the outlier loci that were detected between habitats within each location. The NJ tree based on these 147 outliers shows that the populations group first by habitat and only then by country within each habitat (Figure 5a). The two Ursholmen populations cluster together within the "sheltered" clade confirming the absence of a moderately exposed site within this location in our dataset. Consequently, this site was excluded from the discussion on the main patterns revealed by this study, except when otherwise  Table S1 and Figure S1 between the two clusters, suggesting migration and interbreeding between individuals from the two habitats.

| Genetic substructure based on nonoutlier loci
In contrast, the NJ tree for nonoutlier loci (380 loci, combining the  Figure S4A). The genetic differentiation between locations is relatively higher in UK ( Figure S4B), followed by Sweden and Norway.
Within each country, the populations from sheltered habitats showed higher differentiation than the exposed ones, except in Norway where F ST between Syltøyna and Seløyna is close to zero in both sheltered and moderately exposed habitats ( Figure S4D). Pairwise F ST between populations living in different habitats within each location was zero, whereas this was generally not the case between populations from the F I G U R E 3 Results of the geometric morphometric (GM) analyses of shell shape. (a) Mean shape of individuals of each habitat and location based on the two first principal components of the PCA, where mean values are indicated by coloured symbols, whereas individuals are represented in grey; square and circle symbols represent the moderately exposed and sheltered habitats, respectively. (b) UPGMA dendrogram based on the Procrustes shape distances. (c) Deformation grids depicting the average shape of individuals of each habitat and location compared to the global mean. Mean shapes were magnified x2 to improve visualization. For Ursholmen, only samples of the sheltered habitat were included (see main text). Population codes are the same as in Table S1 and Figure S1  SLok 3 (9.7) 2 (9.1) 2 (6.9) 3 (10.7) 3 (8. same habitat among locations, apart from moderately exposed locations from Norway ( Figure S4 and Table S3).
Using a PP ≥ 0.90 criterion to classify individuals as one of two ecotypes, 49 individuals over a total of 83 analysed for both genetics and shape show concordant classification between the two types of information, with the majority (96%, N = 47) of these samples assigned to the cluster/form typical of the habitat where they were sampled (Table S4). The proportion of concordance was similar between sheltered and moderately exposed habitats (59% and 55%, respectively), and although this varies substantially between sites, we refrain to draw any conclusions given the small sample sizes per site (from 1 to 11, Table S4). Some individuals that are genetically classified as pure for the typical sheltered cluster show a typical shape of the exposed habitat (N = 9), and vice versa (N = 5).
However, most of these individuals (N = 9) belong to the genetic cluster typical of the habitat where they were sampled ( Figure S5, top left and bottom right). Finally, genetically admixed individuals (N = 14) present a wide range of shell shape, with the majority (N = 9) showing a morphology typical of moderately exposed habitats.

| D ISCUSS I ON
Studies of systems like L. fabalis, where different ecotypes coexist across a species´ distribution range, offer important insights about the genetic variation underlying traits involved in local adaptation and ecological speciation. Combining an AFLP genome scan for F I G U R E 5 Neighbour-Joining trees based on Nei´s genetic distance calculated with AFLP-SURV in two sets of loci: (a) outliers 95 (p > .95, DFDIST) (147 loci) and (b) nonoutliers (p < .80) (380 loci), in both cases combining the loci from all comparisons (see main text). Population codes are the same as in Table S1 and Figure S1. Dark squares represent populations from exposed sites, whereas white circles represent populations from sheltered sites. Only bootstrap support values above 70 are shown outlier loci with geometric morphometric analysis, we made the first attempt to detect genetic variation associated with divergent ecotypes in terms of wave-exposure and/or related ecological factors (environmental or biological) and to evaluate the level of sharing of these outlier loci between ecotypes across populations from the UK, Sweden and Norway.

| The role of natural selection on phenotypic divergence
In agreement with previous studies (Kemppainen et al., 2005(Kemppainen et al., , 2011Tatarenkov & Johannesson, 1998), mean shell size was larger in moderately exposed than in sheltered sites (Figure 2a). The remarkable F I G U R E 6 STRUCTURE plots (2 ≤ K ≤ 4) for: (a) three hierarchical levels for K = 2 (global, regional, and local; see Methods) using outlier loci (p > .95, DFDIST) (147 loci) and (b) nonoutliers (p < .80) (380 loci), in both cases combining the loci from all comparisons (see main text). Population codes are the same as in Table S1 and Figure S1. Each cluster is represented by a different grey shade. For the global analysis, shades are comparable across all countries and locations, whereas for the regional analysis shades are comparable only within countries. Shades are not comparable among locations in the local analysis

(b)
influence of habitat type in moulding size variation is evident in the size-based UPGMA dendrograms where moderately exposed sites (with a single exception) consistently cluster together, irrespective of their geographic origin (Figure 2b). Although shell shape variation was also influenced by geography, divergence between populations from the two contrasting habitats was significant across locations ( Figure 3a). Thus, despite local-specific phenotypic effects, we found consistent divergence in shell morphology between L. fabalis populations at the extremes of similar environmental transitions across localities from the same or different European regions. Although this could not be assessed in Ursholmen, observations from a moderately exposed site recently visited show that this pattern holds for this island too, at least for size (R. Faria, personal observation).
The patterns of phenotypic and genetic variation can be explained by both genetic and environmental factors. Although a plastic component cannot be excluded, the association between size/growth and genetic variation found by Tatarenkov and Johannesson (1998) in individuals originally from the same intermediately exposed habitat suggests an important heritable component. This is in agreement with variation in shell morphology being largely genetically inherited in other species within the genus, such as L. saxatilis (Conde-Padín et al., 2009;Conde-Padin et al., 2007;Galindo et al., 2019;Hollander & Butlin, 2010;Johannesson & Johannesson, 1996;Johannesson et al., 1997) and L. subrotundata (Boulding & Hay, 1993;Kyle & Boulding, 1998). Thus, shell divergence between habitats across multiple populations of L. fabalis strongly suggests local adaptation to different levels of wave exposure and/or related ecological factors. For example, a larger shell aperture is observed in exposed sites when compared with sheltered sites in L. fabalis (Figure 3c), which agrees with earlier findings in L. saxatilis (Johannesson et al., 2010). In contrast, L. fabalis from central and northern Europe are larger in moderately exposed habitats than in sheltered habitats (Reimchen, 1981;Tatarenkov & Johannesson, 1998), whereas the opposite is true in L. saxatilis (Johannesson et al., 2010).
However, both the sheltered and moderately exposed habitats of L. fabalis correspond to the "crab" sites of L. saxatilis, whereas the "wave" sites of L. saxatilis are far more exposed and crab free.
A previous ecological study of the two L. fabalis ecotypes suggested that morphological differences are due to a complex interaction between wave exposure and the algae canopy inhabited by the snails, which provides shelter against crabs living in the boulders below (Kemppainen et al., 2005). According to this hypothesis, individuals with a larger and thicker shell would be favoured against predators in moderately exposed habitats where dislodgment is more likely. The shape patterns observed here are consistent with this hypothesis as the large shell aperture characteristic of individuals from moderately exposed habitats can accommodate a larger foot and protect them against dislodgement. In this sense, the moderately exposed ecotype seems to be adjusted for both higher predation (large size) and avoiding dislodgement (large aperture relative to body size), in contrast to L. saxatilis where the exposed ecotype does not need protection against crab predation and is only adapted to fit into crevices and cling tightly to the rocks to avoid wave dislodgement (reviewed in Johannesson et al., 2010).
Strikingly, the moderately exposed population of Syltøyna (Norway) groups with the sheltered populations when considering both size and shape (Figures 2b and 3b). This suggests that selection related with wave action is weaker in this site, which is supported by our observations in the field. Independently of the causes, this highlights that besides general habitat-related differences there are also important location-specific effects on shell morphology (Table S2).

| Habitat-related genetic differentiation
As commonly found in other studies (reviewed by Ravinet et al., 2017), our results reveal heterogeneous differentiation between ecotypes, with a relatively small proportion of loci showing high levels of differentiation likely resisting the substantial gene flow that seems to be eroding differentiation to very low levels in the rest of the genome here assessed. However, an important caveat of genome scans is that the number (and proportion) of discovered outliers depends on how stringent is the cut-off used to detect them (Faria et al., 2014). Here, we used different methods, including BAYESCAN that is known to be conservative (Pérez-Figueroa et al., 2010), as well as different stringency levels to detect outliers. For instance, when considering the less stringent criterion and method (DFDIST, p > .95), we detected ~7% of outliers across locations, which is similar to that detected between other pairs of divergent ecotypes in a wide range of taxa (reviewed in Nosil et al., 2009), including the Crab and Wave ecotypes of L. saxatilis . Nonetheless, as with any genome scan, these outliers must be seen as a list of potential candidate loci influenced by divergent selection that need further confirmation by alternative approaches, as suggested by Ravinet et al. (2017).
Outlier loci identified by means of genome scans could also result from factors not related with local adaptation, such as background selection (Cruickshank & Hahn, 2014;Ravinet et al., 2017  this study), these patterns are compatible with differences in the strength of selection, which can result in geographic heterogeneity under migration-selection equilibria. Nevertheless, future studies of hybrid zones based on cline analysis of loci across the genome, like the one implemented by Westram et al. (2018), will be important to confirm the role of divergent selection in this system.

| Overlap between outliers across locations
A small number of outliers (3, p > .99) were shared among all countries but these estimates of outlier sharing need to be interpreted with caution. Although the genome size of L. fabalis is currently unknown, data from closely related species (L. obtusata and L. saxatilis) suggest a genome size of 1.  (Johannesson et al., 2010;Kemppainen et al., 2011). Alternatively, similar patterns could also have resulted from shared standing genetic variation inherited from a common ancestral population (Westram et al., 2016). In fact, accumulated evi-  (Berner & Roesti, 2017;Haenel et al., 2018;Roesti et al., 2012). This is particularly true within inverted regions, where recombination is heavily reduced in heterokaryotypes, as shown in multiple systems (Wellenreuther & Bernatchez, 2018), including L. saxatilis Morales et al., 2019). The presence of inversions in L. fabalis has been previously hypothesized based on strong association between snail size, Ark genotype, and a RAPD locus genotype (Johannesson & Mikhailova, 2004;Tatarenkov & Johannesson, 1999). Since the location of the AFLP loci (anonymous markers) in the L. fabalis genome is unknown, a contiguous reference genome together with a high-resolution genome scan and linkage maps will be needed to understand the genomic architecture of adaptation in this system and how it relates with the recombination landscape. In particular, the characterization of transects across the environmental gradient using targeted-capture or whole-genome sequencing, as those performed in L. saxatilis Westram et al., 2018), will be important to inform us about the role of inversions in L. fabalis diversification.
The contrast between the patterns of genetic structure based on outlier loci and on putatively neutral variation is noteworthy. The clustering of populations from both ecotypes by geography, when considering nonoutlier loci, suggests parallel evolution of L. fabalis ecotypes (Figures 4 and 5). However, gene flow between populations from different habitats is likely high enough to generate a pattern of parallel evolution, even if they had a single origin and came into secondary contact across multiple locations (Faria et al., 2014).
Thus, a modelling approach, possibly based on approximate Bayesian computation and using sequence data across the genome needs to be implemented in future studies to infer the demographic history of ecotypes and specifically test whether parallel evolution, as observed in L. saxatilis , is more likely than a single origin of ecotypes. Nevertheless, the fact that almost half of all outliers (79 of 147 for p > .95) are private to one location suggests that at least some components of divergent evolution are site-specific. Consequently, even if some of the genetic variation involved in ecotype differentiation had a single origin and was shared after secondary contact, each ecotype is likely to follow its own evolutionary trajectory in each location.

| Implications for the study of speciation
Studies of local adaptation in intertidal habitats of rocky shores are key to quantify the contribution of ecological speciation to marine biodiversity (Sanford & Kelly, 2011 (Saltin et al., 2013). However, while mate choice is known to play an important role in reproductive isolation in L. saxatilis (Perini et al., 2020), whether this partial reproductive barrier is likely to result in significant reproductive isolation remains to be evaluated in natural populations for L. fabalis. Nevertheless, the maintenance of habitat-related phenotypic differentiation in L. fabalis despite high gene flow likely involves multiple loci influencing different traits, including shell thickness, growth and shape, suggesting some degree of reproductive isolation generated by extrinsic (ecological) factors. Although more detailed studies are needed to confirm and extend some of the reported results, the heterogeneous patterns of differentiation here identified, where a relatively small proportion of loci resist the homogenizing effects of gene flow in comparison with lack of differentiation at most studied loci, is compatible with initial stages of ecological speciation (Nosil, 2012;Seehausen et al., 2014). The L. fabalis system thus comprises an interesting system from the marine environment where information from multiple instances of divergence across an environmental transition opens doors for further studies on how populations cope with environmental changes (Sgrò et al., 2011) and how different reproductive barriers accumulate during speciation.

ACK N OWLED G M ENTS
We would like to thank Carolina Pereira for her help in processing and extracting DNA of some samples. This study was supported by:

CO N FLI C T O F I NTE R E S T
The authors declare no conflicts of interest.