Mapping species birth across the recombination landscapes of marine snails

Mapping species birth across the recombination landscapes of marine snails | Authorea try { document.documentElement.classList.add('js'); } catch (e) { } var _gaq = _gaq || []; _gaq.push(['_setAccount', 'G-8VDV14Y67G']); _gaq.push(['_trackPageview']); (function() { var ga = document.createElement('script'); ga.type = 'text/javascript'; ga.async = true; ga.src = ('https:' == document.location.protocol ? 'https://ssl' : 'http://www') + '.google-analytics.com/ga.js'; var s = document.getElementsByTagName('script')[0]; s.parentNode.insertBefore(ga, s); })(); Skip to main content Preprints Collections Wiley Open Research IET Open Research Ecological Society of Japan All Collections About About Authorea FAQs Contact Us Quick Search anywhere Search for preprint articles, keywords, etc. Search Search ADVANCED SEARCH SCROLL Molecular Ecology This is a preprint and has not been peer reviewed. Data may be preliminary. 22 January 2025 V1 Latest version Share on Mapping species birth across the recombination landscapes of marine snails Authors : Emily Giles 0000-0002-3297-3013 [email protected] , Romuald Laso-Jadart 0000-0001-8410-1121 , Daniel Ortiz-Barrientos 0000-0002-7493-416X , Paulina Carimán , Erwan Delrieu-Trottin , Marie-Laure Guillemin , Stefano Mona , Xavier Pochon , and Pablo Saenz-Agudelo 0000-0001-8197-2861 Authors Info & Affiliations https://doi.org/10.22541/au.173753320.03570913/v1 Published Molecular Ecology Version of record Peer review timeline 413 views 268 downloads Contents Abstract INTRODUCTION METHODS RESULTS DISCUSSION FIGURE LEGENDS ACKNOWLEDGEMENTS DATA ACCESSIBILITY BENEFITS GENERATED AUTHOR CONTRIBUTIONS Supplementary Material References Information & Authors Metrics & Citations View Options References Figures Tables Media Share Abstract Understanding the drivers of heterogenous genomic divergence is essential for uncovering the mechanisms that generate and constrain biodiversity. The extent to which adaptation and speciation are facilitated by reorganization of the recombination landscape remains untested in many systems. Marine ecosystems, with their dynamic and fluid habitats, offer a compelling context to investigate genomic divergence. In this study we mapped genomic divergence and selection across recombination landscapes of parapatric marine snail sister species that we show have recently undergone secondary contact. Regions of reduced recombination were enriched for genes exhibiting signatures of negative selection, whereas regions of high recombination were associated with genes under putative positive selection. Notably, the recombination landscape of the population in parapatry of one species (Scurria viridula) differs markedly from that of the other population within this same species, highlighting the role of introgression in reshaping recombination landscapes. In the other species (Scurria zebrina), conservation of the recombination landscape and divergent selection among populations suggest trapping of beneficial allele combinations in regions of low recombination maintain the identity of this species. Among species, signals of divergence with gene flow consistently cluster within specific genomic regions characterized by high recombination rate variation among the populations of S. viridula. These results challenge theoretical expectations of recombination evolution by showing that the causes of genomic divergence can be population- specific. This study demonstrates that recombination landscapes are key modulators of genomic divergence, with contemporary evolutionary shifts that could enable populations to adapt to distinct environments. Our findings provide new insights into the interplay between recombination, selection, and gene flow during speciation, underscoring the complexity of evolutionary trajectories in marine systems. INTRODUCTION Empirical evidence has shown that genomic divergence is heterogenous (Kautt et al. 2020; Plessy et al. 2024; Roux et al. 2016; Wolf and Ellegren 2017; Ravinet et al. 2017). This finding has raised questions about what patterns of genomic divergence can tell us about speciation, particularly the evolution of reproductive isolation. Some students of speciation have suggested that genomic divergence primarily helps quantify reductions in gene flow over time (Westram et al. 2022). In contrast, others argue that understanding these patterns and their causes is important for understanding the mechanisms that generate biodiversity (Moyle 2022). Heterogeneous patterns of genomic divergence can result from multiple, non-mutually exclusive mechanisms. These genomic signatures are detectable and improve our understanding of evolution (Via 2012; Strasburg et al. 2012; Nachman et al. 2012; Fan et al. 2012). Divergent selection can promote and maintain differentiation at specific loci in a neutral and undifferentiated genomic background (Nosil 2009; Feder et al. 2012). To detect such loci, we would expect advantageous alleles to be in genes clustered throughout the genome and displaying high relative (F ST ) and absolute divergence (D XY ). These regions of divergence are typically impervious to gene flow and are expected to differ between parapatric populations and sister species (Nosil et al. 2009). Models of continuous divergence with gene flow are unfortunately often indistinguishable from a model of secondary contact where gene flow erodes past genetic differentiation in regions without barrier loci (Barton and Bengtsson 1986, Nosil et al. 2009). In both scenarios, regions of high relative and absolute divergence are expected. In cases of secondary contact, previously isolated populations come back into contact, allowing gene flow to resume. Importantly, adaptive alleles do not necessarily have to be located within the regions of high genomic divergence (Cruickshank and Hahn 2014) or confined to areas of low recombination (Cruickshank and Hahn 2014). However, regions with reduced recombination can facilitate the maintenance of divergence by preserving advantageous gene combinations (Cruickshank and Hahn 2014; Burri et al. 2015). To differentiate between divergent selection with ongoing gene flow and secondary contact eroding past divergence, we must estimate gene flow through time (ex: Rougemont et al. 2023). This temporal analysis helps determine whether gene flow has been continuous, supporting ongoing adaptation despite interbreeding, or if there was a distinct period of isolation followed by recent mixing of genetic material, which would indicate that gene flow is now reducing previously established genetic differences. Linked selection further complicates the interpretation of genomic divergence patterns. Genomic divergence can result from the hitchhiking of neutral alleles that are linked to beneficial mutations during selective sweeps (Maynard Smith and Haigh 1974, Nosil et al. 2009). When a beneficial allele spreads through a population, nearby alleles can increase in frequency, reducing genetic variation in that region. Moreover, negative selection against deleterious alleles (background selection) can similarly reduce variation around those loci (Hill and Robertson 1966; Charlesworth et al. 1993). Consequently, both positively and negatively selected loci tend to cluster in regions of reduced recombination. This clustering can mimic patterns expected from divergent selection alone (Cruickshank and Hahn 2014), making it challenging to identify the drivers of genomic divergence. To accurately interpret these patterns, it is necessary to account for the effects of linked selection by incorporating recombination rates and distinguishing between direct selection on specific loci and indirect effects caused by linkage. In marine systems, evidence of divergence in the presence of gene flow is abundant. For example, episodic gene flow between freshwater and oceanic three-spine stickleback populations allows for adaptive divergence in the genomes of individuals residing in both habitats (Marques et al. 2016) and genomic divergence detected in expanding ranges of clownfish populations is facilitated not just by drift paired with low Ne, but by adaptive selection aided by moderate migration (Clark et al. 2021). A recent review modeling divergence across 66 marine species, including 28 species of Chordata, 14 species of Mollusca, and 8 species of Cnidaria (De Jode et al. 2023), hypothesizes that divergence with gene flow is widespread in marine taxa. Conversely, a study in intertidal limpets found limited evidence for divergence with gene flow at early-stages of the speciation continuum (Carimán Soto et al. in review ). This supports the idea that divergence with gene flow might not be more common in marine compared to terrestrial organisms (Roux et al. 2016). Thus, periods of speciation in allopatry followed by secondary contact in later stages of speciation might be more prevalent in the marine realm than previously thought (De Jode et al. 2023). In-depth studies of genomic divergence mechanisms in marine organisms will help clarify whether speciation with continuous gene flow is common in these ecosystems. The intertidal of the southeastern Pacific (SEP) possesses distinct characteristics that influence species’ range distributions and population dynamics across its expensive span of over 3,000 km (Broitman et al. 2001). Three well-recognized biogeographic provinces have been defined: the Peruvian province (4°S to 30°S), the Magellanic province (42°S to 56°S), and a transition zone in between these provinces from 30°S to 42°S (Camus 2001). Population structure of macroalgae (Fraser et al. 2010, Montecinos et al. 2012, Tellier et al. 2009), invertebrates (Ewers-Saucedo et al. 2016, Haye et al. 2014, Quesada-Calderón et al. 2021, Peluso et al. 2023, Saenz-Agudelo et al. 2022, Sánchez et al. 2011, Zakas et al. 2009), and fish (Gonzalez-Wevar et al. 2015) align with these biogeographic provinces. Despite strong evidence for biogeographic barriers influencing gene flow in this region, other studies report high levels of panmixia (Chilean abalone, Cárdenas et al. 2009; squat lobster, Haye et al. 2010, jumbo squid, Ibáñez et al. 2011). This suggests that specific life history traits determine the degree to which long-range gene flow is possible in the SEP. Marine snails in the SEP are an interesting model for exploring genomic divergence and speciation. Unlike most gastropods, members of the subclass Patellogastropoda, including the genus Scurria endemic to the SEP, reproduce via broadcast spawning and form larvae that remain pelagic for up to 10 days (Nakano and Sasaki 2011). This planktonic release of gametes may facilitate gene flow across broad geographic scales and enable colonization of new habitats (Burgess et al. 2015). Despite this, comparisons of the genomes of three Scurria species has shown that the genomes of this group are structurally stable, notably in content of coding material. However, functional divergence in duplicated and unique genes suggests that ecological factors drive the diversification of this group (Giles et al. in press ). Furthermore, all Scurria are intertidal grazers (Camus et al. 2008) but have vastly different ecologies and physiologies. For example, the sister species Scurria viridula and Scurria zebrina aggregate in shelters or other protected microhabitats but they do not form aggregations at local scales (Aguilera et al. 2013) and likely compete where their ranges overlap (Aguilera et al. 2020). Additionally, the thermal breadth of S. zebrina is several degrees higher than that of S. viridula (Broitmain et al. 2018), and S. viridula has been shown to evade elevated temperatures while S. zebrina clamps down to the substrate (Broitman et al. 2018). Data suggest that in the past several decades, the range of S. viridula has expanded southwards such that the distributions of these two species currently overlap in a narrow transition zone from 30 to 32°S (Aguilera et al. 2018) where these species potentially hybridize (Saenz-Agudelo et al. 2022). Here we used a comparative population genomics framework to uncover genomic footprints of divergent selection, secondary contact, and recombination rate variation, and we explore how these have shaped speciation of Patellogastropod limpets. We performed whole-genome resequencing of the sister snail species, S. viridula and S. zebrina , to characterize the genomic architecture of divergence among populations within species, among parapatric populations across species, and among allopatric populations across species. Additionally, we estimated recombination rates, searched for signals of selection and inferred the demographic history of the two sister species. We hypothesized that genomic divergence among populations within species would show signatures of divergent selection given potentially high levels of gene flow. In contrast, we predicted that divergence among species would indicate the influence of secondary contact following the southward range expansion of Scurria viridula . Our results provide evidence of the diversity of mechanisms shaping divergence even among closely related species. METHODS Sample collection and Sequencing Muscle tissue from a total of 92 samples ( S. viridula , and S. zebrina ) were collected from sites indicated in Supplemental Table 1, Figure 1A. The sample sites were selected to represent a) the majority of the range distributions of these species (12- 33°S S. viridula ; 31- 42°S S. zebrina ), b) areas where the ranges of the two species overlap (30- 35°S), and c) across common biogeographic breaks that have been shown to structure populations within species (22-25°S and 25-29°S S. viridula ; 31-34°S S. zebrina ; Saenz-Agudelo et al. 2022). With this sampling scheme, we aimed to encompass potential divergence expected for these sister species that have diverged approx. 15mya (Giles et al. in press ). DNA was extracted using the GenJet Genomic DNA purification kit (Thermo Scientific, Waltham, MA, USA) following the manufacturer’s protocol. Quality of DNA extractions were accessed using a Nanodrop and 1.5% agarose gels. Whole genome sequencing was performed using DNBseq at Beijing Genome Institute (BGI). Adapters were then removed, and sequences were quality filtered using standard parameters in SOAPnuke (Chen et al. 2018) allowing less than or equal to 0.1% ambiguous bases, discarding trimmed reads shorter than 10bp, removing reads with quality less than 0.4, and allowing for 25% mismatches to adapter sequences for adapter removal. Alignment and Variant calling The quality filtered reads were aligned to the Scurria scurra reference genome (Giles et al. in press ) using bwa-mem2 v. 2.2.1 (Li and Durbin 2009). Quality of the alignments was assessed using SAMtools v. 1.6 flagstat (Danecek et al. 2021). Then, BCFtools v. 1.9 (Danecek et al. 2021) was used to generate genotype likelihoods and call genotypes, ignoring indels. These were then quality filtered using standard procedures retaining loci with mapping quality > 30, base quality > 20, depth > 4, maximum missing data per loci = 0.5, maximum missing data per sample = 0.1. Separate files were created for variant (minor allele count >1), invariant (minor allele frequency = 0), and combined variant invariant sites. Admixture A Principal Component Analysis (PCA) was performed to summarize genotypic variation across sampled individuals. The PCA was generated using approximately 2 million SNPs filtered for linkage disequilibrium. An admixture analysis was also performed to evaluate the number and composition of genetic groups within the dataset for use in further analyses and test for potential introgression among species. Admixture was assessed by considering the first two eigenvectors of PCA performed in EIGENSOFT v (Patterson et al. 2006; Price et al. 2006). Then, ancestral admixture of each individual was estimated using maximum likelihood in ADMIXTURE v.1.3 (Alexander et al. 2009) and evaluating K = 1 to 9. The Ks with the lowest five-fold cross validation errors (CV) were compared to the PCA to determine genetic groups. Here, K populations were selected to minimize the CV while exploring the greatest range of potential populations. We believe that this more relaxed definition of genetic groups allowed us to account for a wider range of potential divergence scenarios. Introgression We used three analyses to quantify and time potential gene flow between populations of S. viridula and S. zebrina . First, we used TreeMix v.1.12 (Pickrell and Pritchard 2012) to explore whether a phylogeny incorporating all individuals of each population followed a simple bifurcating tree; we explored relationships among K genetic populations with potentially 0-5 admixture events. For this, two samples of Scurria scurra (Carimán Soto et al. in review ) were incorporated in the vcf following the same quality control protocol as mentioned above. These samples were used as the outgroup. Second, we used an ABBA-BABA test to identify recently derived alleles between populations experiencing migration (as per the TreeMix results), and we determined where these alleles occur in the genome. For this, excess of shared derived alleles (f d ) was measured in 20kb windows across all chromosomes for the population comparisons that fit the topology (((P1, P2), P3), O) and showed consistent migration from the TreeMix results to warrant testing of allele sharing between P3 and P2 (here P1 = 4sz, P2 = 3sz, P3 = 2sv, see Results for details). Two samples of Scurria scurra (from Carimán Soto et al. in review ) were used as the outgroup. Then, windows containing significantly high f d values (>0.125; Van Belleghem et al. 2021) that coincided with D XY outliers were identified as genomic regions containing loci having shared derived alleles indicative of recent introgression. Third, to better characterize the timing of gene flow, while accounting for variation in effective population size, demographic modelling was conducted via fastsimcoal2 (Excoffier and Foll 2011). Four models were considered to evaluate gene flow among parapatric populations of the sister species S. viridula and S. zebrina (2sv and 3sz) while considering gene flow between populations within S. viridula (1sv and 2sv) and S. zebrina (3sz and 4sz): a model of no migration, a model of recent migration postdating the population splits of each species, a model of ancient migration predating the population splits, and a model of continuous migration. In all models, an ancestral population of size NANC_S splits at TS into two species Sz and Sv of sizes NANC_Sz nd NANC_Sv respectively. Each species then split at times TSv and TSz into four populations of sizes N1Sv, N2Sv, N3Sz and N4Sz respectively. In recent times, we estimated one migration rate for each species (mig_Sv and mig_Sz). We also designed, depending on the model : i) recent migration rates between 2Sv and 3Sz (RecmigSvSz and RecmigSzSv for the recent and the continuous migration models), ii) migration rates between Sv and 3Sz or between 2Sv and Sz (IntmigSvSz and IntmigSzSv for the continuous migration model), iii) ancient migration rates between Sz and Sv (for the ancient and the continuous migrations models). For each pair of populations, a 2D-Site Frequency Spectrum (2D-SFS) was computed on a dataset filtered for SNPs with missing data and a heterozygous rate > 0.8, and 100 runs of each model were conducted. The likelihood of each model was evaluated by comparing the Akaike information criterion (AIC) scores. Parametric bootstrap computation for the best model was run on 100 generated datasets, selecting the best run among 30 simulations for each dataset. Modelling was conducted using a mutation rate of 0.1x10-8 per base pair per generation. This mutation rate has been used for land snails (Chueca et al., 2021). Additionally, the Nacella mutation rate has been estimated at 1% per million years for cytochrome oxidase (González-Wevar et al. 2016). Furthermore, we used the generation time of one year for the demographic analyses based on observed time to maturity for Patella ranging from 1-2 years (Guallart et al. 2020). Divergence Population genetic statistics were estimated to measure diversity within and divergence among the populations. Nucleotide diversity, 𝜋 (Nei and Li 1979), was assessed for each population, and pairwise F ST (Weir and Cockerham 1984) and D XY (Korunes and Samuk 2021) were measured among populations within species and among sister species. These statistics were averaged across nonoverlapping windows (VCFtools v. 0.1.16, Danecek et al. 2011) for all nine chromosomes (Giles et al. in press ). Nonoverlapping windows of a set length (in bp) were used to facilitate comparisons across population pairs that could differ in position and size if windows were defined by number of SNPs. Window size was chosen to minimize the sampling error of estimates while reducing the number of tests performed. We used a spline analysis (GenWin, Beissinger et al. 2015) based on F ST calculations of sites (VCFtools v. 0.1.16, Danecek et al. 2011) along the first 50Mb of chromosome 1 to determine means and upper bounds of window sizes (Supplemental Figure 4); here one within species population (1sv_2sv) and one among species (1sv_3sz) comparison were evaluated. Approximate mean window sizes (20Kb) were used for statistic calculations, and the approximate upper bound (250Kb) was used for visualizations to reduce within group variation in order to graphically explore differences among groups. Windows (20Kb) with high (outlier) F ST and D XY values were identified based on 99 th and 95 th quantiles, respectively. These are hereafter referred to as windows of divergence with gene flow (Cruickshank and Hahn 2014; Irwin et al. 2016, 2018). Abundance of windows of divergence with gene flow were quantified for each population comparison and chromosome. Windows with F ST and D XY outliers that were separated by no more than 80Kb were pooled to estimate average size of islands of divergence for each population. Additionally, positions in genomic windows of divergence with gene flow were classified as genic or intergenic and lists of genes containing divergence windows were retrieved using bedtools (Quinlan & Hall, et al. 2010). Selection The footprints of selection on protein coding genes were assessed for each population individually. First, degenotate (gitHub harvardinformatics) was used to count synonymous and non-synonymous mutations in all protein-coding genes for each population. For this, counts were derived for each individual compared to the reference S. scurra, and then these were summed for the population. Custom scripts were then used to calculated pN/pS ratios by considering degeneracy (GitHub emgiles). Selection on genes was then determined using the standard approach (positive selection: pN/pS > 1; negative selection: pN/pS 5; negative selection: pN/pS < 0.01). This two-pronged approach was considered to get a broad scale picture of the landscape of selection across the genome, and to reduce the signal to noise ratio when determining potential functions under selection, where the mutations of genes with moderately inflated or deflated levels of pN/pS could be lost in time before becoming substitutions. Then, the de-novo genome annotation of S. scurra was used with bedtools to extract genes present in each genomic window of divergence with gene flow for each population comparison. These genes were compared with those under selection to identify population-specific differences in the number of genes experiencing positive and negative selection within genomic regions of divergence under gene flow. Enrichment of functional ontologies of genes in windows of divergence were determined for each population individually and then compared among populations to determine functions shared among populations. The top 30 most significantly enriched nodes were retained using topGO (Alexa and Rahnenfuhrer 2023). Recombination Variation in recombination along the nine chromosomes was assessed individually for each population down sampling to a maximum of 25 individuals per population to facilitate the use of precalculated likelihood tables. First, the variant call files (vcfs) were phased using Beagle (v4.1, Browning and Browning 2007). Following this, the program interval from LDHat (McVean and Auton 2007) was used to estimate recombination among pairs of SNPs spanning sliding windows. Here 100Kb windows with 10Kb overlaps were sampled every 2,000 iterations of 1,000,000 total iterations and with a block penalty of 5. Theta per site was set as 0.001. Results are shown as average ⍴= 4Ne*r in 250Kb windows in order to compare to F ST and D xy averages. We then estimated conservation of the recombination landscape by estimating correlation coefficients for the linear relationship between ⍴ of pairwise populations. We explored the contribution of recombination to diversity by estimating correlation coefficients for the relationship between ⍴ and 𝜋 within populations. The effect of recombination on divergence with gene flow was assessed by plotting windows of genomic divergence along the recombination landscapes and by comparing boxplots of ⍴ of windows of divergence with gene flow to that of non-windows for each population. Finally, the contribution of linked selection to divergence was determined by plotting genes under selection along the recombination landscapes to determine if windows of divergence occurring in regions of decreased recombination coincide with clusters of positively and negatively selected genes. RESULTS Data Summary Approximately 4GB of quality filtered sequencing data were retained for each of the 92 samples, representing at least 10X coverage of the reference genome. All sequences had greater than 90% mapping to the reference genome, and missing genotypes were similar across samples (Supplemental Figure 1). For one sample, viridula_dnbseq, we had 30X coverage which is reflected in its low % of missing data (Supplemental Figure 1). The filtered dataset contained over 24 million variant loci; the dataset filtered for linkage disequilibrium included approximately 2 million SNPs. Admixture The PCA indicated that over 20% of genetic variation was explained by the first three components when all samples were included in the analysis (Figure 1B). Notably, PC1 discriminated well between the two species and suggested the presence of some potentially admixed S. viridula individuals (viridula_ZA_10, viridula_HU_01; Figure 1B). Furthermore, PC2 and PC3 discriminated populations within species. The results of the PCA were consistent with the ADMIXTURE analysis (Figure 1C) where K=2 had the lowest cross-validation error clearly separating the two species. K=3 to K=5 had CV values differing by less than 0.04 (Supplemental Figure 2) and separated populations within S. zebrina first (K=3) followed by separation of populations within S. viridula (K=4). The ADMIXTURE results consistently showed admixture in the viridula_ZA_10 and viridula_HU_01 individuals, suggesting hybridization. Given the results of PCA and ADMIXTURE together K=4 populations were considered for subsequent analyses. We refer to these populations (1sv, 2sv, 3sz, 4sz) in all subsequent analyses. Thus forward, comparisons among populations within species include 1sv_2sv for S. viridula and 3sz_4sz for S. zebrina , among parapatric populations across species is 2sv_3sz, and among allopatric populations across species includes 1sv_3sz, 1sv_4sz, and 2sv_4sz. Introgression The TreeMix results showed that the log likelihood of migration events for K=4 populations reached an asymptote at 3 migration events (Figure 2A). Migration among parapatric populations across species (2sv_3sz) was consistently recovered and further explored using the ABBA-BABA test f d (Martin et al. 2015). Genomic windows on chromosomes 4, 5, 7 contained loci with shared derived alleles (f d > 0.125, and D XY outliers) indicative of recent introgression (Supplemental Table 2). There were nine genes in these windows, and five of these showed signs of putative positive selection in both 2sv and 3sz (Supplemental Table 2; Note, here the standard metric of pN/pS >1 was used to indicate putative positive selection). All of the nine genes were annotated as hypothetical proteins without full GO id assignment. This being said, two of these nine genes had homologs to known genes in other organisms. These included: g22040 which is homologous to genes involved in immune (Cd164L2) and nervous system functions (Kcng4); and gene g6498 which is homologous to genes involved in oxidative stress response (alt5 and GSTM1) and transcription regulation (Nr2c2). We constructed four demographic models to evaluate the timing and magnitude of gene flow among the parapatric populations across species (2sv and 3sz) while considering gene flow between populations within species (1sv_2sv and 3sz_4sz). These included a model of no migration, a model of recent migration postdating the population splits of each species, a model of ancient migration predating the population splits, and a model of continuous migration (Supplemental Figure 3). The most likely model based on AIC scores included recent migration between 2sv and 3sz (Figure 2B), with Nm ranging from 0.13 to 0.16 between the two species. Population sizes were very high for each population (from respectively), and intra-species migration rates revealed high gene flow within species, from Nm = 3 to Nm = 853 within Sz and Sv, respectively. Full results of models can be found in Supplemental Table 3. Bootstrap computations offered unreliable parameter distributions and were not used for confidence interval computation. Genetic diversity and divergence Across most positions along all chromosomes, nucleotide diversity was similar among populations though slightly lower in S. viridula (1sv mean 𝜋=0.015, 2sv mean 𝜋=0.015, 3sz mean 𝜋=0.018, and 4sz mean 𝜋=0.019) (Figure 3, Table 1). Across all nine chromosomes, mean F ST among populations across species was moderate (0.489 to 0.492). F ST averaged across all positions was extremely low among populations of S. viridula (0.007) and was low among S. zebrina populations (0.027). Across all genomic positions and group comparisons, D XY clearly showed two levels of divergence with comparisons among populations across species ranging from 0.046 – 0.047 and comparisons among populations within species ranging from 0.015 to 0.020 (Figure 3, Table 1). A total of 17,210 20Kb genomic windows were analyzed for all comparisons, yielding variable amounts of windows of divergence with gene flow among populations within and across species. Very few windows of divergence with gene flow were recovered among S. viridula populations, while the most windows of divergence with gene flow in the entire dataset were found among S. zebrina populations (Figure 4). Windows of divergence with gene flow were found in different positions and chromosomes among populations within the species. Positions of windows of divergence were very consistent among allopatric populations across species, and many, though, not all of these were shared with the across species parapatric comparison (Figure 5). Among populations within species, windows of divergence with gene flow were found in chromosomes 1, 2, 3, 5, and 8 for S. viridula and in chromosomes 1, 2, 3, 4, 5, 7, and 8 for S. zebrina . Windows of divergence with gene flow were found in chromosomes 1, 4, 5, 6, 7, 8 and 9 when comparing populations among species (Figure 5). Among populations within species, larger than average (>20Kb) genomic windows of divergence with gene flow were found on chromosomes 1, 3, and 5 ( S. viridula ) and chromosomes 1, 3, 4, and 5 ( S. zebrina ; Figure 6). Genomic windows of divergence with gene flow were consistently greater than 20Kb on chromosomes 4, 6, and 9 for comparisons among populations across species (Figure 6). Genomic windows of divergence with gene flow contained a total of 574 genes across all comparisons (Figure 6). Genomic windows of divergence with gene flow containing high amounts of genes (n > 50) were found among populations of S. zebrina on chromosomes 4 and 5, and among populations across species on chromosomes 6 and 9. Selection Across the genomes, many genes were found to be subject to putative selection when considering standard (positive pN/pS >1; negative pN/pS 5; negative pN/pS < 0.01; Supplemental Figure 5, Table 2 on the diagonal), and differences among populations were minimal. The number of genes under selection that were located in genomic windows of divergence with gene flow were few compared to the overall number of genes showing signs of selection across the genomes (bottom of Figure 7; Table 2). Populations of S. zebrina had the most genes in genomic windows of divergence with gene flow subject to selection (Table 2). We highlight the functions of some of the genes that differed in selection among populations within S. zebrina : protein processing such as protein synthesis (positive selection in 3sz; chrm 5; unnamed gene), dephosphorylation (negative selection in 3sz; chrm 6; g16992; SSU72), and cellular interaction with the environment (positive selection in 3sz; chrm 9; g9073; STBD1, Ttc) (Supplemental Figure 6). In contrast, among S. viridula populations, genes differing in selection included those involved in RNA splicing (positive selection in 2sv; chrm 3; g25406), rRNA acetylation (negative selection in 2sv; chrm 5; g4106; tmcA, At3g57940), and cell cycle regulation (negative selection in 1sv; chrm 5; g5093; CDC14, IGAX, GATB) (Supplemental Figure 6). Differences among species in genes under selection included those involved in DNA repair, splicing, and genome stability (chrm 5; g4103, g4104; HERC2, CHD8, ddx4) which were under positive selection in S. viridula but neutral in S. zebrina . Also among species a gene involved in nucleic acid synthesis and cell division (chrm 4; g22121; Spta1, purB, tmc1) was under negative selection in S. zebrina but not in S. viridula (Supplemental Figure 6). Recombination Recombination landscapes were well conserved among populations of S. zebrina (R 2 all_chromosomes 0.86, Table 3, Figure 7 top panel). Conservation of recombination landscapes was moderate and similar between the S. viridula populations (R 2 all_chromosomes = 0.69) and across species comparisons (R 2 all_chromosomes ranged from 0.58 to 0.70). The relationship between genetic diversity (𝜋) and population scaled recombination (⍴) was consistently strongest for one of the S. viridula populations (1sv), while it was weakest for almost all chromosomes for the other S. viridula population (2sv) (Table 4). Interestingly, this relationship between genetic diversity (𝜋) and population scaled recombination (⍴) was lower for the S. viridula and S. zebrina populations in parapatry. These results suggest that the contributions of recombination to diversity could vary even among populations within species, and could be shaped by hybridization. The contribution of recombination to divergence with gene flow differed among parapatric populations of the two species (3sv and 4sz). Specifically, windows of genomic divergence with gene flow were found in regions of increased recombination rate for the S. viridula population (2sv) that is parapatric with the S. zebrina distribution. Conversely, windows of divergence were found in regions of decreased recombination for the S. zebrina population (3sz) in the parapatric zone (Figure 8). Additionally, windows of divergence with gene flow among species were often located in regions of high variation in recombination among populations. Finally, genomic regions of reduced recombination were characterized as harboring positively and negatively selected genes while regions of increased recombination rate harboured positively selected genes, and this was consistent across all populations (Figure 7 bottom panel). DISCUSSION This study supports the growing body of literature indicating that the landscape of divergence is heterogenous during speciation. We challenge theoretical expectations by showing that recombination landscapes are key modulators of genomic divergence, with contemporary evolutionary shifts that could enable speciation and population divergence when closely related species come into contact. We find that as expected, regions of reduced recombination harbor large swaths of genes showing signs of negative selection while regions of high recombination are associated with genes under putative positive selection. Windows of genomic divergence sometimes coincided with areas of reduced recombination partially supporting the ”speciation island” model (Noor and Bennett 2009; Rieseberg 2001; Burri et al. 2015). However, this pattern was inconsistent across all chromosomes and comparisons. The high conservation of recombination landscapes within S. zebrina aligns with models of recombination rate conservation over short evolutionary timescales, where recombination hotspots are known to be ephemeral, but broad measures of recombination rate variation are thought to be conserved (Dumont and Payseur 2008, 2011; Stevison et al. 2016; Liu et al. 2018). In contrast, the variable relationship between population scaled recombination and genetic diversity in S. viridula populations (1sv and 2sv) indicates that the interplay between recombination and diversity may be population specific. This supports recent work proposing that recombination can be under selection during speciation (Ortiz-Barrientos et al. 2016). These findings highlight the need for more comprehensive theoretical frameworks that can account for the diverse recombination rate evolution patterns observed in natural systems undergoing speciation. Signals of divergence with gene flow were high and shared positions among populations across species. Demographic modelling, despite inconclusive confidence intervals, showed that the best model included low recent gene flow between species though the true underlying model is likely more complex. This, together with a test for allele sharing suggest that historically S. viridula and S. zebrina have diverged in allopatry, yet secondary contact has led to hybridization among parapatric 2sv and 3sz populations. Ecological and physiological studies of these sister species of Scurria have also hypothesized that they have only come into recent contact due to the southward range shift of S. viridula brought on with global climate change (Broitman et al. 2018). Our results support this idea as the timing of secondary contact was found to post date that of the species split and could be very recent. Interestingly, hybrids were only recovered in S. viridula populations, suggesting that viable introgression might be unidirectional and potentially related to missegregation during meiosis or aneuploidy (Butlin 2005) as has been found in other systems as a result of chromosomal inversions (Fishman et al. 2013; Lee et al. 2017; Hooper et al. 2019) and fusions (Yoshida et al. 2023) and anti-recombination (Bozdag et al. 2021). Recombination heterogeneity is shared among closely related species and populations of Pacific oyster (Teng et al. 2023); yet recent work shows that chromosomal inversion polymorphisms maintain species barriers in other marine snails (Le Moan et al. 2024) and a bivalve (Hollenbeck et al. 2022). A lack of fine-scale recombination maps in mollusks limits our understanding of the degree to which these inversions and other genomic features modulate recombination, divergence, and ultimately speciation in this diverse phylum. In terms of functional divergence among S. viridula and S. zebrina , genes found in windows of divergence with gene flow and involved in DNA repair, splicing, and genome stability were under selection in S. viridula but neutral in S. zebrina . Additionally, other genes involved in nucleic acid synthesis and cell division were under negative selection in S. zebrina but neutral in S. viridula . Interestingly, the functions of genes in windows of divergence subject to selection in S. viridula matched functions of duplicated and unique genes in this species (Giles et al. in press), suggesting that evolution of transcription regulation has greatly affected divergence in S. viridula . Indeed, in other mollusks, evolution of splicing (Huang et al. 2016) and DNA repair (Regan et al. 2021) have been linked to stress adaptation; which could be the case for the range expanding population of S. viridula as it moves southward (Saenz et al. in review; Broitman et al. 2018). Differences in patterns of genomic divergence found among populations within species suggest that incipient divergence is species-specific. This has also been found for other species of Scurria (Carimán Soto et al. in review ). Here, many windows of genomic divergence with gene flow were recovered for populations of one species ( S. zebrina ) while few windows of genomic divergence were recovered for the other species ( S. viridula ); this is in line with that previously found for these species using RADseq (Saenz et al. 2022). Additionally, recombination landscapes were conserved among S. zebrina populations, showing correlations among populations higher than that found for populations of D. melanogaster (Adrion et al. 2020). Additionally, windows of genomic divergence within this species contained more genes than did windows within the other species. Together, this suggests that the contributions of recombination to diversity are consistent within this species, while divergent selection in regions of high recombination is driving divergence among populations. Specifically, genes involved in protein synthesis, dephosphorylation, and cellular interaction with the environment showed signs of selection in the S. zebrina population in parapatry with S. viridula . In other mollusks (Prado-Alvarez et al. 2009; Stanovova et al. 2022) and in humans (Tang et al. 2023) genes with these functions have been shown to play key roles in immune response and could potentially be important to early divergence in S. zebrina . Here we found that divergence within S. viridula is very low, yet conservation of the recombination landscape among populations was only moderate. These results suggest that variance in recombination rates among Mollusks (Stapley et al. 2017) has likely been underestimated. Similar to that seen in only a few hundred generations of diverging populations of mice (Wang et al. 2017), we show that fine-scale recombination patterns can vary among incipient populations within species. Furthermore, contributions of recombination to diversity differed among populations of S. viridula . Specifically, the relationship between nucleotide diversity (𝜋) and population scaled recombination (⍴) was highest for the non-introgressed population of S. viridula (1sv) while it was lowest for the population in parapatry and hybridizing with S. zebrina (2sv). Additionally, windows of divergence in the introgressed S. viridula population were associated with high recombination. Together, these results suggest that reduced recombination is a strong determinant of diversity in the non-introgressed population of S. viridula (1sv), while elevated recombination offers novel allele combinations to selection during range expansion and introgression. Indeed, more genes involved in RNA splicing, RNA acetylation, and cell cycle regulation were subject to selection in the introgressed population of S. viridula (2sv) than in the non-introgressed population (1sv). While RNA splicing and acetylation are tightly linked to environmental temperature (Biamonti and Caceres 2009, Sas-Chen et al. 2020) and can cause a number of phenotypic changes (Reddy 2007), these functions are also important during hybridization and can influence nucleolar dominance (Zhou et al. 2011; Oliveira et al. 2006; Scascitelli et al. 2010). It is yet unknown how selection in these genes might affect gene expression or other biological processes across S. viridula populations, but these results suggest that they are important to incipient divergence in this species. We have used robust measures to assess factors contributing to genomic divergence in true limpets, yet some limitations of the current work should be discussed to promote future work in this system. Future work should be directed towards understanding the mechanisms promoting variation in fine-scale recombination patterns. Chromosome level assemblies for S. zebrina and S. viridula would allow for assessment of how large-scale rearrangements and fusions (Rudman et al. 2018) affect recombination and divergence, however the draft genomes of S. viridula and S. zebrina (Giles et al. in press ) could be used in the meantime to investigate the effects of small-scale features such as inversions. The ever-decreasing cost of sequencing technologies has spurned a massive revolution in what can be done with genomic scale resources. Pan-genomes now available for a number of taxa have shown that the idea of a “fixed” species genome is inaccurate (Vernikos et al. 2015). Future work exploring hemizygosity in this group, would be interesting to see if the high gene presence/absence variation found in bivalves (Cacino et al. 2021) also holds for patellogastropods and influences recombination rate variation across short-time scales. It should also be noted that the amount of genes in Scurria assemblies that have functional annotations is only around 50% and extremely lacking for genes with recently derived alleles. The lack of functional annotation greatly restricts our interpretation of the effects of a gene on organism function. In this genome-rich age, the next hurdle to overcome will be full and accurate annotations for taxa spanning the tree of life; genome-wide association and expression studies can be useful for filling these gaps (Giuffra et al. 2019; Kellis et al. 2014). Lastly, characterization of recombination patterns in other species of Scurria spanning clades will help determine the extent to which recombination rate variation contributes to later stages of speciation in this group. Here, among and within closely related species, we show that genomic divergence is shaped by secondary contact which is related to modification of recombination landscapes and few but potentially meaningful differences in genes under selection. Recombination patterns in Scurria limpets likely influence the adaptation and evolution of these species. In S. zebrina , the recombination landscape appears stable and conserved, suggesting that advantageous gene combinations are maintained across the species’ range. This preservation may facilitate adaptation to similar environmental pressures. In contrast, in Scurria viridula , limited recombination in some genomic regions could preserve species-specific traits by reducing mixing of adaptive alleles. Consequently, parapatric populations of Scurria species may maintain significant ecological and genetic divergence despite coexisting in the same environments; though the extent to which this is the case across clades of Scurria remains to be tested. Additionally, different populations of S. viridula show varying relationships between recombination and genetic diversity, suggesting population-specific adaptation to distinct environmental conditions. In regions of low recombination, genetic incompatibilities may accumulate more rapidly, potentially slowing adaptation unless structural variations, such as chromosomal inversions, capture blocks of preexisting adaptive alleles—a phenomenon warranting further investigation. Conversely, areas with high recombination might facilitate rapid adaptation but risk disrupting beneficial gene combinations (Ortiz-Barrientos et al. 2016). Overall, a balance between these mechanisms likely contributes to the heterogeneous genomic divergence observed among Scurria species. Future studies targeting other species of this genus will elucidate the extent to evolved recombination plays a main role in speciation within Patellogastropods. Certainly! Apologies for the previous omissions. Below is the complete LaTeX document that includes all the requested sections, arguments, code snippets, and proofs, organized logically into a single cohesive document. “‘latex FIGURE LEGENDS Figure 1 (A) Historical distribution ranges of S. viridula (purple vertical line) and S. zebrina (blue vertical line) and sample sites (points, purple S. viridula and blue S. zebrina ). Sample site details can be found in Supplemental Table 1. PCA of genetic variation including the two species (B) with admixed S. viridula individuals labeled next to points. Admixture plot (C) for k means clustering of genetic variation K=2 to K=5 is shown from top to bottom. Genetic groups (K=4) are shown at the bottom of the admixture plot and are organized by species and by sample site (north – south). Figure 2 (A) Evolutionary trees highlighting genetic exchange among populations. One to four migration events were considered for the K=4 populations; S. scurra was used as an outgroup to root the tree. Log-likelihood of migration events is given at the bottom right. (B) Demographic modeling of changes in migration and lineage splits for two S. viridula (1sv and 2sv) and two S. zebrina populations (3sz and 4sz). Four models were considered: a model of no migration, a model of recent migration postdating the population splits of each species, a model of ancient migration predating the population splits, and a model of continuous migration. Fit of model is given with AIC scores. Number of migrations per generation (Nm), effective population sizes (Ne) and time of population (Tsv and Tsz) and species (TS) splits were inferred. Figure 3 Genetic diversity (lowest panels) within groups across chromosomes 1-9. Colors of lines correspond to K=4 groups shown in Figure 1 (purple S. viridula , blue S. zebrina ). Absolute (middle panels) and relative (top panels) genetic divergence among groups; purple shows comparisons among S. viridula populations, blue shows comparisons among S. zebrina populations, and green and yellow show comparisons among the sister species S. viridula and S. zebrina . Figure 4 Number of genomic windows (20Kb) with high relative (F ST ) and absolute (D XY ) divergence based on 99 and 95% quantiles, respectively. Number of windows is shown for population comparisons. Figure 5 Distribution of windows of genomic divergence with gene flow across chromosomes (ss1-9) for population comparisons. Figure 6 Size of windows of divergence with gene flow across chromosomes for group comparisons. Windows in adjacent 80Kb regions were pooled, and bars represent average size of windows in bp. Black points represent the number of genes found at least partially within all windows of divergence with gene flow for given chromosomes and group comparisons. Figure 7 Mapping divergence and selection (bottom of plots) across population-scaled recombination landscapes (purple and blue lines on top of plots) for nine chromosomes and four populations. Populations are as follows: 1sv (light purple line), 2sv (purple line), 3sz (bright blue line), 4sz (light blue line). Location of genomic windows of divergence with gene flow in relation to recombination landscapes are indicated as points below the recombination landscapes. Colors of points are as follows: among populations within S. viridula – purple points; among populations within S. zebrina – blue points, among parapatric populations across species – yellow points, among allopatric populations across species – green points. Footprints of selection in genes are ordered sequentially along chromosomes, shown in the bottom panel. Four selection landscapes are presented for each of the four populations (top – 1sv, second to top – 2sv, second to bottom – 3sz, bottom – 4sz). Here genes are represented by vertical lines which are blue if found to show standard signs of putative positive selection (pN/pS > 1) or red if found to show standard signs of putative negative selection (pN/pS < 1). Functions of genes can be found in Supplemental Figure 6. Figure 8 Comparisons of mean population-scaled recombination within (windows) and outside (non-windows) of windows of genomic divergence with gene flow for four populations of Scurria limpets: 2sv (light purple), 3sv (purple), 4sz (bright blue), 6sz (light blue). Window only includes comparisons involving the focal population. The number of genomic windows considered is given (n), along with p-values and means for two sample t-tests. Certainly! Apologies for the previous omissions. Below is the complete LaTeX document that includes all the requested sections, arguments, code snippets, and proofs, organized logically into a single cohesive document. “‘latex ACKNOWLEDGEMENTS The work herein presented is supported by the FONDECYT grant No. 1190710 and the scholarship BECA CONICYT-PCHA/Doctorado Nacional Chile/2019 No. 2190878. The authors thank members of the Saenz Lab for assistance in the field and processing samples. DATA ACCESSIBILITY Raw sequence reads for S. viridula are being deposited in NCBI (BioProject: forthcoming) Raw sequence reads for S. zebrina have been deposited in NCBI (BioProject: PRJNA1212230) Scripts for analyses (GitHub emgiles) BENEFITS GENERATED A research collaboration was developed with scientists from the countries providing genetic samples, all collaborators are included as co-authors, the results of research have been shared with the provider communities and the broader scientific community. More broadly, our group is committed to international scientific partnerships, as well as institutional capacity building. AUTHOR CONTRIBUTIONS PS-A, ECG, and DO-B conceived the project. PC and PS-A conducted sampling and data generation. ECG, RL-J, and PC contributed to data analysis. ECG wrote the manuscript with edits from all authors particularly MLG, DO-B, and PS-A. Supplementary Material File (figure2.pdf) Download 81.67 KB File (table 1.docx) Download 15.19 KB File (table 2.docx) Download 15.71 KB References 1. Adrion JR, Galloway JG, Kern AD. (2020) Predicting the landscape of recombination using deep learning. Molecular Biology and Evolution , 37 (6):1790–1808. Alexa A, Rahnenfuhrer J. (2023) topGO: Enrichment analysis for gene ontology. R package version 2.54.0. Available at: https://bioconductor.org/packages/topGO. Barton N, Bengtsson BO. (1986) The barrier to genetic exchange between hybridising populations. Heredity , 57 (3):357–376. Beissinger TM, Rosa GJM, Kaeppler SM, Gianola D, de Leon N. (2015) Defining window-boundaries for genomic analyses using smoothing spline techniques. Genetics Selection Evolution , 47 (1). Bozdag GO, Ono J, Denton JA, Karakoc E, Hunter N, Leu JY, Greig D. (2021) Breaking a species barrier by enabling hybrid recombination. Current Biology , 31 (4):180–181. Broitman BR, Aguilera MA, Lagos NA, Lardies MA. (2018) Phenotypic plasticity at the edge: Contrasting population-level responses at the overlap of the leading and rear edges of the geographical distribution of two Scurria limpets. Journal of Biogeography . Broitman BR, Navarrete SA, Smith F, Gaines SD. (2001) Geographic variation of southeastern Pacific intertidal communities. Marine Ecology Progress Series , 224 :21–34. Browning SR, Browning BL. (2007) Rapid and accurate haplotype phasing and missing-data inference for whole-genome association studies by use of localised haplotype clustering. The American Journal of Human Genetics , 81 (5):1084–1097. Burri R, Nater A, Kawakami T, Mugal CF, Olason PI, Smeds L, Suh A, Dutoit L, Bures S, Garamszegi LZ, Hogner S. (2015) Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers. Genome Research , 25 (11):1656–1665. Butlin RK. (2005) Recombination and speciation. Molecular Ecology , 14 :2621–2635. Carimán Soto PJ, Guillemin ML, Giles EC, Narváez Guíñez GF, Suescún VA, Sáenz-Agudelo P. in review . Exploring evolutionary mechanisms of genomic divergence in Scurria limpets. https://doi.org/10.21203/rs.3.rs-4682442/v1 Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li S, Li Y, Ye J, Yu C, Li Z, Zhang X. (2018) SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience , 7 (1):gix120. Calcino AD, Kenny NJ, Gerdol M. (2021) Single individual structural variant detection uncovers widespread hemizygosity in molluscs. Philosophical Transactions of the Royal Society B , 376 :20200153. Camus PA. (2001) Biogeografía marina de Chile continental. Revista Chilena de Historia Natural , 74 :587–617. Cárdenas L, Castilla JC, Viard F. (2009) A phylogeographical analysis across three biogeographical provinces of the south-eastern Pacific: the case of the marine gastropod Concholepas concholepas . Journal of Biogeography , 36 (5):969–981. Charlesworth B. (1993) Directional selection and the evolution of sex and recombination. Genetics Research , 61 (3):205–224. Cheng JY, Mailund T, Nielsen R. (2017) Fast admixture analysis and population tree estimation for SNP and NGS data. Bioinformatics , 33 (14):2148–2155. Clark RD, Aardema ML, Andolfatto P, Barber PH, Hattori A, Hoey JA, Montes HR Jr, Pinsky ML. (2021) Genomic signatures of spatially divergent selection at clownfish range margins. Proceedings of the Royal Society B , 288 (1952):20210407. Cruickshank TE, Hahn MW. (2014) Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. Molecular Ecology , 23 (13):3133–3157. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Handsaker RE, Lunter G, Marth GT, Sherry ST, McVean G. (2011) The variant call format and VCFtools. Bioinformatics , 27 (15):2156–2158. De Jode A, Le Moan A, Johannesson K, Faria R, Stankowski S, Westram AM, Butlin RK, Rafajlovic M, Fraisse C. (2022) Ten years of demographic modelling of divergence and speciation in the sea. Evolutionary Applications , 15 (3):383–402. Dumont BL, Payseur BA. (2008) Evolution of the genomic rate of recombination in mammals. Evolution , 62 (2):276–294. Dumont BL, Payseur BA. (2011) Evolution of the genomic recombination rate in murid rodents. Genetics , 187 (3):643–657. Ewers-Saucedo C, Pringle JM, Sepúlveda HH, Byers JE, Navarrete SA, Wares JP. (2016) The oceanic concordance of phylogeography and biogeography: a case study in Notochthamalus . Ecology and Evolution , 6 (13):4403–4420. https://doi.org/10.1002/ece3.2205 Excoffier L, Foll M. (2011) Fasimcoal: a continuous-time coalescent simulator of genomic diversity under arbitrarily complex evolutionary scenarios. Bioinformatics , 27 (9):1332–1334. Fan S, Elmer KR, Meyer A. (2012) Genomics of adaptation and speciation in cichlid fishes: recent analyses in African and Neotropical lineages. Philosophical Transactions of the Royal Society B , 367 (1587):385–394. Feder JL, Egan SP, Nosil P. (2012) The genomics of speciation with gene flow. Trends in Genetics , 28 (7):342–350. Fishman L, Stathos A, Beardsley PM, Williams CF, Hill JP. (2013) Chromosomal rearrangements and the genetics of reproductive barriers in Mimulus (monkey flowers). Evolution , 67 (9):2547–2560. Fraser CI, Thiel M, Spencer HG, Waters JM. (2010) Contemporary habitat discontinuity and historic glacial ice drive genetic divergence in Chilean kelp. BMC Evolutionary Biology , 10 :203. https://doi.org/10.1186/1471-2148-10-203 Giuffra E, Tuggle CK, FAANG Consortium. (2019) Functional Annotation of Animal Genomes (FAANG): current achievements and roadmap. Annual Review of Animal Biosciences , 7 (1):65–88. Giles EC, Gonzalez V, Carimán P, Leiva C, Suescún AV, Lemer S, Guillemin M-L, Ortiz-Barrientos D, Saenz-Agudelo P. in press Comparative genomics points to ecological drivers of genomic divergence among intertidal limpets. Molecular Ecology Resources. González-Wevar CA, Salinas P, Hüne M, Segovia NI, Vargas-Chacoff L, Astorga M, Cañete JI, Poulin E. (2015) Phylogeography in Galaxias maculatus (Jenyns, 1848) along two biogeographical provinces in the Chilean coast. PLoS One , 10 (7):e0131289. https://doi.org/10.1371/journal.pone.0131289 Guindon S, Dufayard JF, Lefort V, Anisimova M, Hordijk W, Gascuel O. (2010) New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Systematic Biology , 59 (3):307–321. Haye PA, Salinas P, Acuña E, Poulin E. (2010) Heterochronic phenotypic plasticity with lack of genetic differentiation in the southeastern Pacific squat lobster Pleuroncodes monodon . Evolution & Development , 12 (6):628–633 Haye PA, Segovia NI, Muñoz-Herrera NC, Gálvez FE, Martínez A, Meynard A, Pardo-Gandarillas MC, Poulin E, Faugeron S. (2014) Phylogeographic structure in benthic marine invertebrates of the southeast Pacific coast of Chile with differing dispersal potential. PLoS One , 9 (2):e88613. https://doi.org/10.1371/journal.pone.0088613 Hill WG, Robertson A. (1966) The effect of linkage on limits to artificial selection. Genetics Research , 8 (3):269–294. Hoang DT, Chernomor O, von Haeseler A, Minh BQ, Vinh LS. (2018) UFBoot2: Improving the ultrafast bootstrap approximation. Molecular Biology and Evolution , 35 (2):518–522. Hollenbeck CM, Portnoy DS, Garcia de la Serrana D, Magnesen T, Matejusova I, Johnston IA. (2022) Temperature-associated selection linked to putative chromosomal inversions in king scallop ( Pecten maximus ). Proceedings of the Royal Society B , 289 :20221573. Hooper DM, Griffith SC, Price TD. (2019) Sex chromosome inversions enforce reproductive isolation across an avian hybrid zone. Molecular Ecology , 28 (6):1246–1262. Huang B, Zhang L, Tang X, et al. (2016) Genome-wide analysis of alternative splicing provides insights into stress adaptation of the Pacific oyster. Marine Biotechnology , 18 :598–609. Ibáñez CM, Cubillos LA, Tafur R, Argüelles J, Yamashiro C, Poulin E. (2011) Genetic diversity and demographic history of Dosidicus gigas (Cephalopoda: Ommastrephidae) in the Humboldt Current System. Marine Ecology Progress Series , 431 :163–171. Irwin DE, Alcaide M, Delmore KE, Irwin JH, Owens GL. (2016) Recurrent selection explains parallel evolution of genomic regions of high relative but low absolute differentiation in a ring species. Molecular Ecology , 25 (18):4488–4507. Irwin DE, Milá B, Toews DP, Brelsford A, Kenyon HL, Porter AN, et al. (2018) A comparison of genomic islands of differentiation across three young avian species pairs. Molecular Ecology , 27 (23):4839–4855. Kalyaanamoorthy S, Minh BQ, Wong TKF, von Haeseler A, Jermiin LS. (2017) ModelFinder: Fast model selection for accurate phylogenetic estimates. Nature Methods , 14 :587–589. Kautt AF, Kratochwil CF, Nater A, et al. (2020) Contrasting signatures of genomic divergence during sympatric speciation. Nature , 588 :106–111. Kellis M, Wold B, Snyder MP, Bernstein BE, Kundaje A, et al. (2014) Defining functional DNA elements in the human genome. PNAS , 111 :6131–6138. Korunes KL, Samuk K. (2021) pixy: unbiased estimation of nucleotide diversity and divergence in the presence of missing data. Molecular Ecology Resources , 21 (4):1359–1368. Lawson DJ, van Dorp L, Falush D. (2018) A tutorial on how not to over-interpret STRUCTURE and ADMIXTURE bar plots. Nature Communications , 9 :3258. Lee C-R, Wang B, Mojica JP, et al. (2017) Young inversion with multiple linked QTLs under selection in a hybrid zone. Nature Ecology & Evolution , 1 :0119. Le Moan A, Stankowski S, Rafajlovic M, Ortega-Martinez O, Faria R, Butlin RK, Johannesson K. (2024) Coupling of twelve putative chromosomal inversions maintains a strong barrier to gene flow between snail ecotypes. Molecular Ecology , in press . Li H, Durbin R. (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics , 25 (14):1754–1760. Liu H, Maclean CJ, Zhang J. (2019) Evolution of the yeast recombination landscape. Molecular Biology and Evolution , 36 (2):412–422. Marques DA, Jones FC, Di Palma F, Kingsley DM, Reimchen TE. (2016) Experimental evidence for rapid genomic adaptation to a new niche in an adaptive radiation. Nature Ecology & Evolution , 2 (7):1128–1138. Martin SH, Davey JW, Jiggins CD. (2015) Evaluating the use of ABBA-BABA statistics to locate introgressed loci. Molecular Biology and Evolution , 32 (1):244–257. Maynard Smith J, Haigh J. (1974) The hitch-hiking effect of a favorable gene. Genetical Research , 23 :23–35. McVean G, Auton A. (2007) LDhat 2.1: A package for the population genetic analysis of recombination. Department of Statistics, Oxford, UK. Montecinos A, Broitman BR, Faugeron S, Haye PA, Tellier F, Guillemin M-L. (2012) Species replacement along a linear coastal habitat: Phylogeography and speciation in the red alga Mazzaella laminarioides along the southeast Pacific. BMC Evolutionary Biology , 12 (1):97. https://doi.org/10.1186/1471-2148-12-97 Moyle LC. (2022) Forty-two, and other precise answers to difficult questions: a commentary on Westram et al., 2022. Journal of Evolutionary Biology , 35 (9):1183–1187. Nachman MW, Payseur BA. (2012) Recombination rate variation and speciation: theoretical predictions and empirical results from rabbits and mice. Philosophical Transactions of the Royal Society B , 367 :409–421. Nei M, Li WH. (1979) Mathematical model for studying genetic variation in terms of restriction endonucleases. Proceedings of the National Academy of Sciences , 76 (10):5269–5273. Noor MAF, Bennett SM. (2009) Islands of speciation or mirages in the desert? Examining the role of restricted recombination in maintaining species. Heredity , 103 :439–444. Nosil P, Funk DJ, Ortiz-Barrientos D. (2009) Divergent selection and heterogeneous genomic divergence. Molecular Ecology , 18 :375–402. Oliveira CI, Bicudo HE, Itoyama MM. (2006) New evidence for nucleolar dominance in hybrids of Drosophila arizonae and Drosophila mulleri . Genetics and Molecular Research , 5 :632–637. Ortiz-Barrientos D, Engelstadter J, Rieseberg LH. (2016) Recombination rate variation and the origin of species. Trends in Ecology and Evolution , 31 (3):226–236. Patterson N, Price AL, Reich D. (2006) Population structure and eigen analysis. PLoS Genetics , 2 (12):e190. Peluso L, Broitman BR, Lardies MA, Nespolo RF, Saenz-Agudelo P. (2023) Comparative population genetics of congeneric limpets across a biogeographic transition zone reveals common patterns of genetic structure and demographic history. Molecular Ecology , 32 (14):3812–3825. Pickrell JK, Pritchard JK. (2012) Inference of population splits and mixtures from genome-wide allele frequency data. PLoS Genetics . Plessy C, Mansfield MJ, Bliznina A, Masunaga A, West C, Tan Y, Liu AW, Grašič J, Del Río Pisula MS, Sánchez-Serna G, Fabrega-Torrus M, Ferrández-Roldán A, Roncalli V, Navratilova P, Thompson EM, Onuma T, Nishida H, Cañestro C, Luscombe NM. (2024) Extreme genome scrambling in marine planktonic Oikopleura dioica cryptic species. Genome Research , 34 (3):426–440. https://doi.org/10.1101/gr.278295.123 Prado-Alvarez M, Gestal C, Novoa B, Figueras A. (2009) Differentially expressed genes of the carpet shell clam Ruditapes decussatus against Perkinsus olseni . Fish and Shellfish Immunology , 26 (1):72–83. Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, Reich D. (2006) Principal components analysis corrects for stratification in genome-wide association studies. Nature Genetics , 38 (8):904–909. Quesada-Calderón S, Giles EC, Morales-González S, Saenz-Agudelo P. (2021) Pinpointing genetic breaks in the southeastern Pacific: Phylogeography and genetic structure in a commercially important tunicate. Journal of Biogeography . https://doi.org/10.1111/jbi.14227. Quinlan AR, Hall IM. (2010) BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics , 26 (6):841–842. Ravinet M, Faria R, Butlin RK, Galindo J, Bierne N, Rafajlovic M, Noor MAF, Mehlig B, Westram AM. (2017) Interpreting the genomic landscape of speciation: a road map for finding barriers to gene flow. Journal of Evolutionary Biology , 30 (8):1450–1477. Reddy ASN. (2007) Alternative splicing of pre-messenger RNAs in plants in the genomic era. Annual Review of Plant Biology , 58 :267–294. Regan T, Stevens L, Peñaloza C, Houston RD, Robledo D, Bean TB. (2021) Ancestral physical stress and later immune gene family expansions shaped bivalve mollusc evolution. Genome Biology and Evolution , 13 (8):evab177. https://doi.org/10.1093/gbe/evab177. Rieseberg LH. (2021) Chromosomal rearrangements and speciation. Trends in Ecology and Evolution , 16 :351–358. Roux C, Fraïsse C, Romiguier J, Anciaux Y, Caltier N, Bierne N. (2016) Shedding light on the grey zone of speciation along a continuum of genomic divergence. PLoS Biology , 14 (12):e2000234. Rudman SM, Barbour MA, Csillery K, Gienapp P, Guillaume F, Hairston NG, Hendry AP, Lasky JR, Rafajlovic M, Rasanen K, Schmidt PS, Seehausen O, Therkildsen NO, Turcotte MM, Levine JM. (2017) What genomic data can reveal about eco-evolutionary dynamics. Nature Ecology and Evolution , 2 :9–15. http://doi.org/10.1038/s41559-017-0385-2. Saco A, Novoa B, Greco S, Gerdol M, Figueras A. (2023) Bivalves present the largest and most diversified repertoire of Toll-like receptors in the animal kingdom, suggesting broad-spectrum pathogen recognition in marine waters. Molecular Biology and Evolution , 40 (6):msad133. Saenz-Agudelo P, Peluso L, Nespolo R, Broitman B, Haye P, Lardies M. (2022) Population genomic analyses reveal hybridization and marked differences in genetic structure and demographic history of Scurria limpet sister species with parapatric distributions across the southeastern Pacific. Ecology and Evolution , 12 (5): e8888 Saenz-Agudelo P, Nespolo R, Broitman B, Caballero P, Lardies M. (unpublished) Phenotypic and genetic diversity in Scurria limpets’ populations is associated with biogeographic zones along the South-Eastern Pacific. Samuk K, Owens GL, Delmore KE, Miller SE, Rennison DJ, Schluter D. (2017) Gene flow and selection interact to promote adaptive divergence in regions of low recombination. Molecular Ecology , 26 (17):4378–4390. Sánchez R, Sepúlveda RD, Brante A, Cárdenas L. (2011) Spatial pattern of genetic and morphological diversity in the direct developer Acanthina monodon (Gastropoda: Mollusca). Marine Ecology Progress Series , 434 :121–131. https://doi.org/10.3354/meps09184. Sas-Chen A, Thomas JM, Matzov D. (2020) Dynamic RNA acetylation revealed by quantitative cross-evolutionary mapping. Nature , 583 :638–643. Scascitelli M, Cognet M, Adams KL. (2010) An interspecific plant hybrid shows novel changes in parental splice forms of genes for splicing factors. Genetics , 184 :975–983. Smith JM, Haigh J. (1974) The hitch-hiking effect of a favourable gene. Genetics Research , 23 (1):23–35. Stanovova MV, Gazizova GR, Gorbushin AM. (2023) Transcriptomic profiling of immune-associated molecules in the coelomocytes of lugworm Arenicola marina (Linnaeus, 1758). Journal of Experimental Zoology Part B: Molecular and Developmental Evolution , 340 (1):34–55. Stapley J, Feulner PG, Johnston SE, Santure AW, Smadja CM. (2017) Variation in recombination frequency and distribution across eukaryotes: patterns and processes. Philosophical Transactions of the Royal Society B: Biological Sciences , 372 (1736): 20160455. Stevison LS, Woerner AE, Kidd JM, Kelley JL, Veeramah KR, McManus KF, Great Ape Genome Project, Bustamente CD, Hammer MF, Wall JD. (2016) The Time Scale of Recombination Rate Evolution in Great Ape. Molecular Biology and Evolution , 33 (4):928–945. Strausfeld NJ, Hirth F. (2016) Introduction to ‘Homology and convergence in nervous system evolution’. Philosophical Transactions of the Royal Society B , 371 (1685): 20150034. Strasburg JL, Sherman NA, Wright KM, Moyle LC, Willis JH, Rieseberg LH. (2012) What can patterns of differentiation across plant genomes tell us about adaptation and speciation? Phil. Trans. R. Soc. B. , 367 :364–373. Tang Q, Liu M, Zhao H, Chen L. (2023) Glycogen-binding protein STBD1: Molecule and role in pathophysiology. J Cell Physiol. , 238 (9):2010–2025. Tedesco M, Giantess F, Lazarevic G, Giansanti V, Romano D, Monzani S, Catalano I. (2022) Chromatin Velocity reveals epigenetic dynamics by single-cell profiling of heterochromatin and euchromatin. Nature Biotechnology , 40 (2):235–244. Tellier F, Meynard AP, Correa JA, Faugeron S, Valero M. (2009) Phylogeographic analyses of the 30°S south-east Pacific biogeographic transition zone establish the occurrence of a sharp genetic discontinuity in the kelp Lessonia nigrescens : Vicariance or parapatry? Molecular Phylogenetics and Evolution , 53 (3):679–693. https://doi.org/10.1016/j.ympev.2009.07.030. Teng W, Fu H, Li Z, Zhang Q, Xu C, Yu H, Kong L, Liu S, Li Q. (2023) Parallel evolution in Crassostrea oysters along the latitudinal gradient is associated with variation in multiple genes involved in adipogenesis. Molecular Ecology , 32 (19):5276–5287. Van Belleghem SM, Ruggieri AA, Concha C, Livraghi L, Hebberecht L, Rivera ES, Ogilvie JG, Hanly JJ, Warren IA, Planas S, Ortiz-Ruiz Y. (2023) High level of novelty under the hood of convergent evolution. Science , 379 (6636):1043–1049. Via S. (2012) Divergence hitchhiking and the spread of genomic isolation during ecological speciation with gene flow. Phil. Trans. R. Soc. B. , 367 :451–460. Wang J, Street NR, Scofield DG, Ingvarsson PK. (2016) Variation in linked selection and recombination drive genomic divergence during allopatric speciation of European and American Aspens. Molecular Biology and Evolution , 33 (7):1754–1767. Wang RJ, Gray MM, Parmenter MD, Broman KW, Payseur BA. (2017) Recombination rate variation in mice from an isolated island. Molecular Ecology , 26 (2):457–470. Weir BS, Cockerham CC. (1984) Estimating F-statistics for the analysis of population structure. Evolution , [pages]:1358–1370. Wersebe MJ, Sherman RE, Jeyasingh PD, Weider LJ. (2022) The roles of recombination and selection in shaping genomic divergence in an incipient ecological species complex. Molecular Ecology , 32 (6):1478–1496. Westram AM, Stankowski S, Surendranadh P, Barton N. (2022) What is reproductive isolation? Journal of Evolutionary Biology , 35 (9):1143–1164. Wolf JB, Ellegren H. (2017) Making sense of genomic islands of differentiation in light of speciation. Nature Reviews Genetics , 18 (2):87–100. Yoshida K, Rodelsperger C, Roseler W, Riebesell M, Sun S, Kikuchi T, Sommer RJ. (2023) Chromosome fusions repatterned recombination rate and facilitated reproductive isolation during Pristionchus nematode speciation. Nature Ecology & Evolution , 7 :424–439. Zakas C, Binford J, Navarrete SA, Wares JP. (2009) Restricted gene flow in Chilean barnacles reflects an oceanographic and biogeographic transition zone. Marine Ecology Progress Series , 394 (2):165–177. https://doi.org/10.3354/meps08265. Zhou R, Moshgabadi N, Adams KL. (2011) Extensive changes to alternative splicing patterns following allopolyploidy in natural and resynthesized polyploids. PNAS , 108 (38):16122–16127. Crossref Google Scholar Information & Authors Information Version history V1 Version 1 22 January 2025 Peer review timeline Published Molecular Ecology Version of Record 21 Sep 2025 Published Copyright This work is licensed under a Non Exclusive No Reuse License. Collection Molecular Ecology Keywords demography limpets mollusk recombination selection speciation Authors Affiliations Emily Giles 0000-0002-3297-3013 [email protected] Universidad Austral de Chile View all articles by this author Romuald Laso-Jadart 0000-0001-8410-1121 Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum National d'Histoire Naturelle, EPHE-PSL, Université PSL, CNRS, SU, UA View all articles by this author Daniel Ortiz-Barrientos 0000-0002-7493-416X The University of Queenland View all articles by this author Paulina Carimán El Instituto de Ciencias Ambientales y Evolutivas View all articles by this author Erwan Delrieu-Trottin Université Montpellier View all articles by this author Marie-Laure Guillemin Universidad Austral de Chile View all articles by this author Stefano Mona Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum National d'Histoire Naturelle, EPHE-PSL, Université PSL, CNRS, SU, UA View all articles by this author Xavier Pochon Cawthron Institute View all articles by this author Pablo Saenz-Agudelo 0000-0001-8197-2861 Universidad Austral de Chile View all articles by this author Metrics & Citations Metrics Article Usage 413 views 268 downloads .FvxKWukQNSOunydq8rnd { width: 100px; } Citations Download citation Emily Giles, Romuald Laso-Jadart, Daniel Ortiz-Barrientos, et al. Mapping species birth across the recombination landscapes of marine snails. Authorea . 22 January 2025. DOI: https://doi.org/10.22541/au.173753320.03570913/v1 If you have the appropriate software installed, you can download article citation data to the citation manager of your choice. Simply select your manager software from the list below and click Download. For more information or tips please see 'Downloading to a citation manager' in the Help menu . Format Please select one from the list RIS (ProCite, Reference Manager) EndNote BibTex Medlars RefWorks Direct import Tips for downloading citations document.getElementById('citMgrHelpLink').addEventListener('click', function() { popupHelp(this.href); return false; }); $(".js__slcInclude").on("change", function(e){ if ($(this).val() == 'refworks') $('#direct').prop("checked", false); $('#direct').prop("disabled", ($(this).val() == 'refworks')); }); View Options View options PDF View PDF Figures Tables Media Share Share Share article link Copy Link Copied! Copying failed. Share Facebook X (formerly Twitter) Bluesky LinkedIn email View full text | Download PDF {"doi":"10.22541/au.173753320.03570913/v1","type":"Article"} Now Reading: Share Figures Tables Close figure viewer Back to article Figure title goes here Change zoom level Go to figure location within the article Download figure Toggle share panel Toggle share panel Share Toggle information panel Toggle information panel Go to previous graphic Go to next graphic Go to previous table Go to next table All figures All tables View all material View all material xrefBack.goTo xrefBack.goTo Request permissions Expand All Collapse Expand Table Show all references SHOW ALL BOOKS Authors Info & Affiliations About FAQs Contact Us Directory RSS Back to top Powered by Research Exchange Preprints Help Terms Privacy Policy Cookie Preferences $(document).ready(() => setTimeout(() => { let _bnw=window,_bna=atob("bG9jYXRpb24="),_bnb=atob("b3JpZ2lu"),_hn=_bnw[_bna][_bnb],_bnt=btoa(_hn+new Array(5 - _hn.length % 4).join(" ")); $.get("/resource/lodash?t="+_bnt); },4000)); (function(){function c(){var b=a.contentDocument||a.contentWindow.document;if(b){var d=b.createElement('script');d.innerHTML="window.__CF$cv$params={r:'9ff4c3653caa300f',t:'MTc3OTM3ODkwMw=='};var a=document.createElement('script');a.src='/cdn-cgi/challenge-platform/scripts/jsd/main.js';document.getElementsByTagName('head')[0].appendChild(a);";b.getElementsByTagName('head')[0].appendChild(d)}}if(document.body){var a=document.createElement('iframe');a.height=1;a.width=1;a.style.position='absolute';a.style.top=0;a.style.left=0;a.style.border='none';a.style.visibility='hidden';document.body.appendChild(a);if('loading'!==document.readyState)c();else if(window.addEventListener)document.addEventListener('DOMContentLoaded',c);else{var e=document.onreadystatechange||function(){};document.onreadystatechange=function(b){e(b);'loading'!==document.readyState&&(document.onreadystatechange=e,c())}}}})();