Genomic divergence between Spanish Littorina saxatilis ecotypes unravels limited admixture and extensive parallelism associated with population history

Abstract The rough periwinkle, Littorina saxatilis, is a model system for studying parallel ecological speciation in microparapatry. Phenotypically parallel wave‐adapted and crab‐adapted ecotypes that hybridize within the middle shore are replicated along the northwestern coast of Spain and have likely arisen from two separate glacial refugia. We tested whether greater geographic separation corresponding to reduced opportunity for contemporary or historical gene flow between parallel ecotypes resulted in less parallel genomic divergence. We sequenced double‐digested restriction‐associated DNA (ddRAD) libraries from individual snails from upper, mid, and low intertidal levels of three separate sites colonized from two separate refugia. Outlier analysis of 4256 SNP markers identified 34.4% sharing of divergent loci between two geographically close sites; however, these sites each shared only 9.9%–15.1% of their divergent loci with a third more‐distant site. STRUCTURE analysis revealed that genotypes from only three of 166 phenotypically intermediate mid‐shore individuals appeared to result from recent hybridization, suggesting that hybrids cannot be reliably identified using shell traits. Hierarchical AMOVA indicated that the primary source of genomic differentiation was geographic separation, but also revealed greater similarity of the same ecotype across the two geographically close sites than previously estimated with dominant markers. These results from a model system for ecological speciation suggest that genomic parallelism is affected by the opportunity for historical or contemporary gene flow between populations.

directions in this research is in uncovering the extent of genetic parallelism and elucidating different processes that could produce varying degrees of this parallelism (Conte, Arnegard, Peichel, & Schluter, 2012;Elmer & Meyer, 2011;Nosil, 2012).
The extent of geographic and genetic separation between populations undergoing parallel adaptation may play a role in determining observed levels of shared genomic differentiation (Conte et al., 2012;Johannesson et al., 2010;Ord & Summers, 2015). Populations that show greater connectivity or a recent common origin are likely to share similar complements of standing genetic variation that may be included in the genetic architecture of ecologically adaptive traits (Barrett & Schluter, 2008). Connected or recently separated populations are also expected to share similar genetic backgrounds that determine constraints on adaptive de novo mutations (Conte et al., 2012). Additionally, migration between populations will dictate the likelihood that ecologically advantageous mutations will be transmitted and selected across sites (Johannesson et al., 2010;Ralph & Coop, 2010;Rosenblum, Parent, & Brandt, 2014).
Adoption of next-generation sequencing methods (NGS) in studies of ecological adaptation has allowed simultaneous identification of genome-wide distributions of locally selected loci and high-resolution estimates of population genetic structure (Bernatchez, 2016;Bonin, Ehrich, & Manel, 2007;Hohenlohe et al., 2010). Genome scans, which detect outlier regions of elevated genomic divergence that may reflect reduced gene flow due to adaptation (Nosil et al., 2009), can be conducted across replicate instances of ecological divergence to quantify genomic parallelism through identifying proportions of outlier sharing (Bernatchez, 2016;Lewontin & Krakauer, 1973;Westram, Panova, Galindo, & Butlin, 2016). However, shared ancestry or contemporary migration between populations may constrain the genetic independence of parallel phenotypic evolution. Observed levels of genomic parallelism can be contextualized by quantifying population genetic structure and phylogenetic distance between instances of parallel adaptation, to determine whether rates of genomic parallelism are elevated in more closely related populations (Lowe & Allendorf, 2010;Pritchard, Stephens, & Donnelly, 2000).
Next-generation sequencing studies in model species for ecological speciation have revealed substantial variation in genomic parallelism between systems, inferred from the proportion of shared genomic outlier regions exhibiting divergence beyond a neutral threshold. Investigation of shared parallel divergence in stick insect ecotypes (Timema cristinae) (Soria-Carrasco et al., 2014) has revealed modest (17%) sharing of genomic divergence between ecotype pairs across multiple populations, whereas comparisons of parallel genomic divergence in freshwater and marine stickleback (Gasterosteus aculeatus) ecotypes have uncovered greater (35%) genomic reuse (Jones et al., 2012). A meta-analysis of rates of parallel molecular evolution from a small sample of studies that passed inclusion criteria, conducted by Conte et al. (2012), identified substantial reuse of individual candidate genes (55%) or genomic regions (32%) associated with a selected trait. This rate of genomic parallelism declined with phylogenetic distance between compared groups, suggesting that the extent of evolutionary independence between lineages is important to how repeated molecular evolution occurs. Given the great heterogeneity of observed patterns, identification of the evolutionary and ecological contexts in which parallel genomic differentiation occurs across systems is a vital next step in understanding the genomics of parallel evolution.
Previous measures of population genetic structure of Spanish ecotypes using microsatellites, neutral allozymes, and AFLP loci have indicated greater divergence between separate sites than between ecotypes at different tidal levels on the same shore, providing support for the hypothesis that evolution of ecotypes has occurred in microparapatry with minimal spatial separation to prevent gene flow (Galindo, Morán, & Rolán-Alvarez, 2009;Rolán-Alvarez et al., 2004). Mitochondrial phylogenetic analysis of crab and wave ecotypes across the full Spanish geographic range has revealed mtDNA haplotype trees that cluster by geography rather than ecotype, indicating that neutral sequence divergence has accumulated between sites isolated by distance, rather than between ecotypes (Quesada et al., 2007). Comparison of different historical demographic hypotheses supports origin of northern and southern populations of Spanish ecotypes from two separate glacial refugia , consistent with previous phylogeographic studies (Doellman, Trussell, Grahame, & Vollmer, 2011;Panova et al., 2011).
Surprisingly, phylogenetic or geographic distance has not yet been shown to affect rates of genomic parallelism across L. saxatilis ecotype pairs. Transcriptomic comparison of parallel divergence across the European range of L. saxatilis ecotypes found only moderate genomic parallelism between Swedish, UK and Spanish populations (Westram et al., 2014) despite the expectation that UK and Swedish sites should show greater similarity due to colonization from shared glacial refugia (Panova et al., 2011). Analysis of RAD-derived SNP markers in Swedish ecotypes found at multiple sites within a single archipelago at short geographic distance (10 km) identified low-to-moderate (8%-28%) genomic parallelism between sites despite restriction of the spatial scale of comparisons (Ravinet et al., 2016). Thus, L. saxatilis ecotypes appear to defy the expectation that reduced separation corresponds to greater sharing of genomic divergence. A genome-wide study including sites at varying levels of geographic and evolutionary divergence has yet to be carried out to clarify the role of connectivity in dictating genomic parallelism.
Here, we study patterns of genomic divergence and parallelism in Spanish L. saxatilis wave and crab ecotypes and putative hybrids of intermediate phenotype in three geographically separated sites. We test the hypothesis that the degree of genomic parallelism associated with adaptation to replicated ecological conditions on the upper or lower shore is affected by the possibility of historical or contemporary gene flow between sites and will show correspondence with geographic separation between sites. We predict that samples from the same tidal level at geographically closer sites that shared a glacial refuge will exhibit less neutral genetic differentiation and will also demonstrate greater levels of shared genomic divergence than more distant sites. We used a recently developed, inexpensive doubledigest RAD sequencing protocol (Kess, Gross, Harper, & Boulding, 2016) to quantify population structure, phylogenetic separation, and parallel genomic divergence between crab and wave ecotype snails from three sites that have colonized intertidal zones from two glacial refugia. Consistent with our prediction, we find levels of genomic parallelism correspond to the degree of geographic separation.
Surprisingly, we identify greater genomic differentiation between tidal levels of the same site than has previously been revealed using neutral molecular markers, and find little evidence of admixture in phenotypically intermediate individuals from the mid-shore. This study contributes to our understanding of factors dictating genomic parallelism and provides new evidence of more extensive genomic isolation and progress toward potential microparapatric speciation in the model "crab" and "wave" ecotypes of L. saxatilis.

| Sampling design and study sites
Samples were obtained from three sites across the known distribu- viously been used to identify individuals likely arising from recent hybridization between ecotypes (Galindo et al., 2013;Johannesson, Johannesson, & Rolán-Alvarez, 1993). By focusing on phenotypic intermediates, we intended to increase the chances of sampling admixed genotypes formed by hybridization between the two ecotypes. We froze all sampled individuals at −80°C and stored them in a −20°C freezer until DNA extraction.

| DNA extraction
We selected 32 individuals randomly from each sampled tidal level (upper, mid-shore, and lower) from each transect for DNA extraction. We separated body tissue from the shell using forceps and then further dissected individuals into foot and visceral mass sections.
Foot tissue was flash-frozen in liquid nitrogen and then ground in a mortar and pestle. Ground tissue was then extracted using E.Z.N.A Mollusc DNA kits (Omega Bio-Tek). DNA quality was visualized on 2% TAE agarose gel against a 100-bp DNA ladder (Invitrogen) and checked for concentration and purity on a NanoDrop spectrophotometer (Thermo Fisher, Inc.).

| Library construction and sequencing
Libraries were constructed using the protocol described in A total of 541 individuals produced PCR products that could be visualized on a 2% TAE agarose gel. From each dual-indexed sample, 10 uL of PCR product was pooled in a 1.5-mL tube with PCR products from other sampled individuals from the same tidal level and site, creating groups of twelve individuals per pool in a final volume of 120 μl. Each tube was concentrated in a SpeedVac for 15 min to a final volume of 60 μl and then size selected in the 300-to 600-bp size range using a Pippin Prep gel cassette system (Sage Biosciences) at the London Regional Genomics Centre in London, Ontario.
Samples were then combined in three libraries with a final volume ≥20 μl and a concentration ≤5 ng/μl, with each library representing a different sampled site (Burela, Corrubedo, Silleiro). Libraries were then sequenced using paired-end, 125-bp sequencing reagent kits on an Illumina HiSeq 2500, at the Clinical Genomics Centre (CGC) in Toronto, Ontario. Sequences were then demultiplexed at the CGC and attributed to a sequenced individual using Hiseq N7 and S5 indexes.

| Sequence quality analysis and control
Read number, base pairs sequenced, Q quality scores per lane and per sequenced individual, average and per base quality score, base content, and adapter content were assessed at the CGC using FastQC 0.11.0 (Andrews, 2014). We conducted all subsequent sequence quality checks and bioinformatic analyses using the Shared Hierarchical Academic Research Computing Network (SHARCNET) Linux cluster. Sequenced adapter content adjacent to restriction fragment inserts was identified at the end of several sequenced DNA fragments. We trimmed sequenced fragments for quality by removing regions of bases that fell below a Q score of 20 within a 10-bp sliding window and removed adapter and index material using Trimmomatic v 0.32 (Bolger, Lohse, & Usadel, 2014). A second round of quality analysis was carried out using the process_radtags module in Stacks 1.40 (Catchen, Amores, Hohenlohe, Cresko, & Postlethwait, 2011) to remove any sequenced fragments with uncalled bases or ambiguous restriction enzyme sequences at the beginning of forward or reverse reads. All reads were truncated to a length of 85 bp, and reads falling below this read length threshold were discarded.
The rescue_radtags option was enabled to retain reads with restriction enzyme cut sequences with individual base mismatch deviations from the expected cut site sequence. found that <5% of loci demonstrated genotypic frequencies that deviated from HWE. These loci were retained for outlier analysis as biological explanations for departure from HWE in a small subset of loci cannot be ruled out and may reflect genomic regions with distinct patterns of admixture rather than error in genotyping or SNP identification (Nosil, Parchman, Feder, & Gompert, 2012).

| De novo assembly, SNP identification, and validation
We then used PGDSpider v 2.1.0.3 Lischer and Excoffier (2012) and VCFtools, and the genepopedit R package (Stanley, Jeffery, Wringe, DiBacco, & Bradbury, 2017) to convert Stacks output files to formats required by software used for outlier analysis, calculation of population structure, and phylogenetic comparison.

| Outlier analysis
Detection of SNPs potentially under selection was carried out through three outlier detection methods. We conducted F ST outlier analysis in ARLEQUIN v 3.5.2.1 (Excoffier & Lischer, 2010) implementing the standard FDIST method described by Beaumont and Nichols (1996). We also identified F ST outliers using BayeScan 2.1 (Foll & Gaggiotti, 2008) and the pcadapt R package (Luu, Bazin, & Blum, 2017). We conducted outlier analysis between upper and lower shore sampled individuals in each site separately using each method. For the FDIST analysis, we carried out 20,000 simulations of 100 demes to build the null distribution for outlier identification. Following outlier analysis, we applied Sequential Goodness of Fit (SGoF) correction at 5% significance to reduce the possibility of spurious positive detection of outlier loci (Carvajal-Rodríguez & de Uña-Alvarez, 2011).
BayeScan analysis was conducted with a sample size of 5000, running 20 pilot runs with 5,000 iterations each, followed by 10 thinning interval runs, each with 5,000 iterations, and burn-in of 100,000 iterations. Loci identified in BayeScan were considered significant using 5% false discovery rate (
Because demographic and evolutionary processes can generate statistical linkage, we categorized potentially physically linked loci as those exhibiting r 2 above 0.2 in all sites.

| Population structure and admixture
We used the Bayesian clustering program STRUCTURE v2.3.4 (Pritchard et al., 2000) to detect population structure and admixture between ecotypes from each sampled tidal level and site. We conducted separate analyses for the total set of 4,256 loci that passed quality filtering, 193 divergent loci, and 2,549 neutrally evolving loci that did not exhibit patterns of differentiation consistent with divergent or balancing selection in any comparison. STRUCTURE analyses were conducted using the parallelstructure R package (Besnier & Glover, 2013), with three replicates per K value for K = 1 to K = 9, with 100,000 burn-in steps and 500,000 MCMC sampling steps.
We visualized results in STRUCTURE HARVESTER (Earl & vonHoldt, 2012) to identify the value of K with highest ΔK (Evanno, Regnaut, & Goudet, 2005). To quantify admixture within and among ecotypes from separate sites, we investigated patterns of population structure of neutral loci for K = 6. Replicate runs for the K value chosen for each set of loci were summarized and graphically displayed using

| Genetic differentiation
We used hierarchical analysis of molecular variance (AMOVA) (Excoffier, Smouse, & Quattro, 1992) implemented in ARLEQUIN to identify the extent of population structure caused by separation between tidal levels (upper and lower shore) and sites. We included only upper and lower shore individuals in all AMOVA tests.

| Phylogenetic analysis
We calculated Cavalli-Sforza and Edwards' chord distances (1967) between all individuals separately for 2,000 randomly selected neutral loci and all 193 divergent outlier loci in POPULATIONS 1.2.33 (Langella, 1999

| Sequence quality analysis and control
We sequenced 541 individuals (173 Burela

| SNP identification
We identified 127,476 unique RAD loci across 477 sequenced individuals, and 300,874 SNP sites were variable within those loci.
Coverage per locus was high, with an average value of 20.55 reads and standard deviation of 6.50 reads per RAD locus across samples.
We obtained a set of 4,256 filtered SNPs present in a minimum of 70% of samples, which we used to assess signatures of selection and population genetic structure.

| Outlier analysis
We uncovered 240 unique loci demonstrating elevated patterns of divergence across all comparisons at 5% significance in FDIST, compared to only 59 unique outlier loci identified at 5% significance in BayeScan (Table 1). All significant loci identified by BayeScan were also identified using FDIST. Using pcadapt, we identified 275 loci that exhibited signals of potential selection and produced separation between ecotype clusters. Of these loci, 193 exhibited overlap with outliers identified using FDIST. A total of 1,098 loci at 5% significance exhibited low differentiation consistent with potential balancing selection occurring within sites using the FDIST method, whereas no loci under balancing selection were detected using BayeScan or pcadapt. BLASTX annotated seven outlier loci when the search was restricted to "Mollusca" and had an e-value cutoff of 10 −8 (Table 2).

| Sharing of outlier loci
Consistent with our prediction, we observed high levels of shared divergent outliers between the closer Corrubedo and Silleiro sites, but low levels of shared genomic divergence with the northern Burela site (

| Linkage disequilibrium
We observed a low level of potential physical linkage among SNPs

| Genetic differentiation
Quantification of genetic differentiation revealed strong roles for both site and tidal level in partitioning genetic variation. In AMOVA conducted with all SNPs, we observed greatest differentiation between site subgroups (F SC = 0.219), indicating that geography is the primary level of differentiation among sites or tidal levels (Table 4).
However, when Burela was excluded from analysis, the greatest differentiation was instead observed between the crab and wave ecotypes collected from the upper and lower tidal levels, respectively (F CT = 0.147). Among neutral loci, the greatest source of genetic differentiation was observed between geographic subgroups (F SC = 0.20). This pattern of differentiation was similar in comparisons of Corrubedo and Silleiro, but differences between the ecotypes from the upper and lower tidal level also contributed to a nearly equal degree to genetic divergence (Table 4) Information Table S1 and S2).

| Parallelism of divergent loci
High levels of shared divergence observed between samples from the same tidal levels in Silleiro and Corrubedo compared with Burela support a role for historical or contemporary separation in affecting levels of genomic parallelism. Genomic regions potentially under divergent selection were largely non-parallel between Burela and Silleiro or Burela and Corrubedo. In contrast, a large proportion (34.4%) of the total identified divergent outlier sites were shared between Silleiro and Corrubedo, even though the latter two sites were significantly differentiated for neutral loci (F SC = 0.126), and phylogenetically distinct. These results differ from those observed in L. saxatilis studies by Ravinet et al. (2016) and by Westram et al. (2014). In these studies, outlier loci sharing between the same ecotypes across the entire eastern Atlantic range did not differ strongly by geographic separation or level of genomic differentiation, and no increase in parallelism was observed even when spatial scale was restricted to ecotype comparisons between geographically close sites in Sweden.
The timing of adaptive divergence of ecotypes following coloni- also revealed a potential biogeographic barrier between northern and southern populations corresponding to variation in climatic and oceanographic features that may limit dispersal between sites in the present study (Piñeira, Quesada, Rolán-Alvarez, & Caballero, 2008). The presence of multiple refugia increases the recent evolutionary independence of northern and southern sites, such that genetic backgrounds of divergent traits may have continued to evolve through periods of isolation to be less suitable for the establishment of introduced alleles from distant sites (Conte et al., 2012). Even without genetic barriers to adaptive gene flow or variation in divergent traits between sites, many shared ecologically relevant traits in ecotypes may possess a polygenic inheritance or may arise from many overlapping pathways to produce similar phenotypes, increasing the probability that alternative alleles may be selected (Hoekstra & Nachman, 2003;McGee, Neches, & Seehausen, 2016). Crosses between individuals of similar ecotype from multiple sites may be useful in determining the genetic architecture of traits divergent between crab and wave ecotypes, and in revealing whether negative epistasis between alleles from different genetic backgrounds has played a role in limiting observed levels of genomic parallelism.
Our results show marked difference from other studies of parallel genomic divergence that emphasize a role for standing variation in driving local adaptation (Barrett & Schluter, 2008;Nelson & Cresko, 2018;Schluter & Conte, 2009). Interestingly, divergent loci exhibited clear separation between refugia in wave ecotypes, whereas shared phylogenetic grouping and population structure clustering were observed for these loci in crab ecotypes. These results suggest that shared variation has produced crab ecotypes, whereas distinct complements of adaptive alleles led to wave ecotype formation. The low overlap of outliers between glacial lineages in L. saxatilis suggests that divergent adaptation to intertidal conditions may be achieved without drawing on species-level ancestral variation. Although separate complements of standing variation that arose after separation of northern and southern sites may have enabled local ecotype formation, our results demonstrate that ancient variation may not be critical for rapid adaptation. These results are contrary to observation in a highly connected model system, marine and freshwater sticklebacks, in which ancient alleles exhibiting high differentiation are shared between replicate populations across a large range (Nelson & Cresko, 2018). This difference may be attributed to greater levels of connectivity between stickleback populations, allowing regular movement of adaptive variation to geographically distant sites, and ensuring a range-wide pool of genomic diversity (Schluter & Conte, 2009). In contrast to another model sys- Darwin's finches (Lamichhaney et al., 2018). Thus, the temporal context of parallel adaptation may also play a role in dictating genomic parallelism by increasing the timeframe during which these events may occur. These hypotheses may be tested in L. saxatilis ecotypes through characterization of molecular diversity within divergent genomic regions and adjacent sites using whole genome sequence.
These analyses, as have recently been conducted in stickleback marine and freshwater ecotypes (Nelson & Cresko, 2018), could allow differentiation between hypotheses positing that selection on ancestrally shared standing variation has produced shared regions of local divergence between ecotypes, or that shared genomic divergence has instead occurred through repeated selective sweeps of the same de novo mutation shared through infrequent migration. structure also exist for observations of regions of increased divergence and have been covered in recent discussions of the limitations of genome scans (see Bierne, Roze, & Welch, 2013;Bierne, Welch, Loire, Bonhomme, & David, 2011;Ravinet et al., 2017;Wolf & Ellegren, 2017). Greater confidence that identified divergent regions also contribute to ecologically relevant traits may be achieved through association of those regions with phenotypic and environmental variables, and through molecular characterization of those regions or linked loci.
Here, we annotated a very small proportion of the F ST outliers that exhibit functions that may be important in ecotype divergence and speciation (~4%) ( Table 2). One of the outliers, spondin-1-like protein, has also been detected between L. saxatilis ecotypes from populations in Britain . This type of protein has been detected as a candidate for shell formation genes in pearl oyster and may play a similar role in ecotype shell architecture (Kinoshita et al., 2011). Moreover, vertebrate SCO-spondin proteins have been regarded as similar to invertebrate gel-forming mucins (Lang et al., 2016) which could also be involved in shell formation and that have also been identified in a previous transcriptome scan in L. saxatilis (Galindo, Grahame, & Butlin, 2010). Two other annotated outliers, axonemal dynein heavy chain and G protein-coupled receptor, have been detected as outliers in Swedish L. saxatilis populations (Ravinet et al., 2016). Axonemal dynein heavy chain could be involved in flagellar movement in sperm and therefore linked to speciation (Carvalho, Lazzaro, & Clark, 2000). G protein-coupled receptors could also be related to adaptation like the case of MC1R in vertebrates (Hoekstra, Hirschmann, & Bundey, 2006). Characterization of the function of these outlier loci across multiple instances of ecotype formation is a clear objective for future follow-up studies.

| Population genetic structure and admixture
Population clustering with RAD-derived SNPs enabled separation of distinct site and tidal level groups with minimal admixture using neutral loci, although the presence of separate glacial refugia dominated the clustering pattern observed in STRUCTURE comparisons.
However, the overall pattern of divergence across measurements of population structure, individual phylogenetic separation, and differentiation indicates the presence of geographically distinct crab and wave ecotype pairs within each site.
We detected much less frequent hybridization between ecotypes than previously identified, as well as greater genetic separation between both tidal level and sites (see Galindo et al., 2013).
Greater precision achieved using a larger set of codominant SNP markers to detect differentiation and admixture could partially account for lower observed amounts of admixture not detected by AFLP markers. Similar improved resolution of population structure using SNPs has revealed greater levels of genomic differentiation than inferred from previous AFLP studies in phenotypically distinct cichlid fish species (Keller et al., 2013), as well as in delineating genetically distinct groups in Wilson's warblers (Cardelina pusilla) (Ruegg et al., 2014).

| Genetic differentiation and phylogenetic analysis
Consistent with observed infrequent admixture and clear genetic clustering of tidal levels, comparisons of sources of genomic differentiation revealed more extensive genetic differentiation between each tidal level than has been previously observed (Galindo et al., 2009;Rolán-Alvarez et al., 2004). Similarly, measurement of neutral F ST (0.17) in this study was substantially larger than values observed in comparisons between ecotypes using AFLP markers (F ST = 0.064; Galindo et al., 2009). Similar levels of divergence have been observed in comparisons between sympatric Lake Victoria cichlids (Keller et al., 2013), whereas slightly greater levels of overall divergence have been observed between sympatric Heliconius host races (Martin et al., 2013). This observed level of baseline genomic differentiation suggests that progression toward speciation between ecotypes in upper and lower tidal levels may be more extensive than previously believed, potentially to an early stage of genome-wide isolation (Feder, Egan, & Nosil, 2012;Seehausen et al., 2014). However, although current results indicate strong isolation and genomic divergence, sampling of RAD markers may not perfectly reflect genomic distribution of divergent sites and corresponding evolutionary history. Future comparison of whole genome sequence from individuals in either tidal level will help clarify the degree of genomic isolation and progress toward speciation between L. saxatilis ecotypes.
Our phylogenetic and genetic differentiation analyses of neutral loci suggest a large degree of genetic independence of ecotypes in each site. Repeated observation of mitochondrial phylogenies that cluster by sampling site rather than ecotype, including between southern sites (Quesada et al., 2007;Tirado et al., 2016), and support for independent ecotype origin between northern and southern sites using demographic modeling , together provide strong evidence for local emergence of ecotype pairs.
Although this scenario cannot be completely rejected between southern sites, the observation of separate neutral SNP phylogenies in the present study and the occurrence of high between-ecotype levels of divergence that nearly match levels of distance driven isolation provide evidence against ongoing gene flow between ecotypes eroding historical patterns of divergence. Additionally, recent simulations of divergence with gene flow suggest that unless high levels of between-ecotype migration occur, phylogenetic methods should be able to differentiate between local divergence and secondary contact (Pérez-Pereira, Quesada, & Caballero, 2017). Thus, our results provide additional support for local origin of ecotypes across sites, but future whole genome studies may more explicitly test this hypothesis through demographic reconstruction with coalescent models.

| CON CLUS IONS
In the present study, we identified substantial sharing of total divergent outlier loci between upper and lower shore ecotypes across two geographically close sites that was not observed at a third distant site colonized from a separate glacial refuge. These results support the hypothesis that ancestral or contemporary connectivity of populations may play a role in dictating genomic parallelism. These results also indicate that shared ancestral variation may not be necessary for ecological adaptation and divergence.
Additionally, we identified very little admixture between upper and lower shore ecotypes, despite the presence of phenotypically intermediate individuals found at the mid-shore. Lastly, we identified a greater role for differentiation between tidal levels contributing to genetic variation than was found in previous studies, potentially enabled by greater genomic sampling. Taken together, these results indicate greater progression of L. saxatilis ecotypes along a continuum of speciation toward genomic isolation than has been previously uncovered for this model system.

ACK N OWLED G M ENTS
Funding for snail collection, DNA library preparation, DNA sizeselection, snail culture assistance, and partial PhD and postdoctoral salary support for T.K. was from a NSERC Discovery Grant (RGPIN-2015-05468) to EGB. Support for T.K. was also provided through an Ontario