Exploring evolutionary mechanisms of genomic divergence in Scurria limpets

doi:10.21203/rs.3.rs-4682442/v1

Download PDF

Article

Exploring evolutionary mechanisms of genomic divergence in Scurria limpets

https://doi.org/10.21203/rs.3.rs-4682442/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Decades of research in population genetics have revealed that genetic divergence between populations and species is not uniformly distributed throughout the genome, but rather exhibits a high degree of heterogeneity. This phenomenon has sparked debate regarding the evolutionary mechanisms responsible for this variability. This study delves into the patterns of genomic divergence in populations of three marine limpet species, Scurria scurra, Scurria araucana, and Scurria ceciliana, which inhabit both sides of two biogeographic breaks (30-34°S and 41-43°S). To investigate this, genomic divergence was analyzed using estimates of F_ST and D_XY with a 20 kb sliding window approach across the entire genomes of 56 individuals. Highly divergent genomic regions ranging from 20 to 280 kilobases in size were identified, primarily associated with the allopatric divergence model in all three species, although significant signals of divergence with gene flow were observed in some regions. Sympatric populations of S. scurra and S. araucana, spanning the 30-34°S latitudinal range, exhibited heterogeneous divergence patterns. In S. ceciliana, the southernmost species with populations across the 41-43°S latitudinal range, the allopatric divergence and divergence with gene flow processes seem to have deeply affected the genomic architecture of divergence. This pattern was not expected for this species with an evolutionary history of populations being isolated in glacial refugia. Genetic divergence appears to be species-specific, with few shared regions, and genes relevant to lipid metabolism, response to oxidative stress, and mitochondrial functions were identified within regions of divergence. These findings highlight the complexity of genomic diversification in marine limpets.

Biological sciences/Genetics/Genomics/Comparative genomics

Biological sciences/Evolution/Evolutionary genetics

Multiple evolutionary forces, such as genetic drift, natural selection, mutation, and migration, contribute to genomic divergence among populations (Nosil et al. 2009). Understanding how these forces interact and shape genomic divergence remains one of the most active fields in evolutionary genomics. Genetic drift stochastically alters allele frequencies, eventually leading to the random fixation or loss of genetic variants, with its intensity and speed depending solely on population size. Consequently, populations diverging solely under genetic drift should exhibit a relatively homogeneous pattern of divergence across the genome (Wu 2001; Gavrilets and Vose 2005). In contrast, populations diverging under natural selection should show heterogeneous patterns of divergence across the genome, where regions directly under selection diverge more rapidly than neutral zones. Under this scenario, the magnitude of divergence and the number of rapidly diverging sites depend on the strength of selection and the number of sites or genes under selection. Generally, the prevalence of genomic heterogeneity has been associated with the predominance of divergent selection. Genomic regions diverging among populations are known as "genomic islands of divergence" (Nosil et al. 2009; Delmore et al. 2018; Van Belleghem et al. 2021; Zhou et al. 2024), and they have been linked to the formation of phenotypic differences (Wang et al. 2020) and barriers to gene flow (Poelstra et al. 2014; Wang and Coop 2022). However, they are not always associated with reproductive isolation (Clarkson et al. 2014) and may result from a hitchhiking effect around selected regions (Feder et al. 2012) or variation in local recombination rate (Cruickshank and Hahn 2014; Burri et al. 2015; Ortiz-Barrientos et al. 2016). Certainly, expanding the number of taxa for which the architecture of genomic divergence is studied can help us determine whether or not certain mechanisms are more general than others.

The observed diversity of patterns of genomic divergence among species has led to a rich debate on the underlying causes of genomic divergence. At least two explanatory models have been proposed to describe the heterogeneous nature of genomic divergence based on the comparisons between patterns of absolute (D_XY) and relative (F_ST) measures of genomic differentiation between groups (Irwin et al. 2016; 2018). These models invoke natural selection as the main mechanism underlying these genomic regions of high differentiation, but the models can be divided into two main classifications: Divergence with gene flow and Divergence in allopatry (Figure 1). The "Divergence with gene flow" model posits that genomic regions with high F_ST values should also exhibit high D_XY values. This is because reproductive isolation caused by barrier loci in these regions impedes gene flow between populations within these areas (Cruickshank and Hahn 2014). Meanwhile, gene flow between populations may be possible in other genome regions, keeping F_ST and D_XY values low. The second model, "Divergence in allopatry" posits that genomic regions with high F_ST values should not exhibit high D_XY values; rather, the increase in relative differentiation in areas with loci under selection is solely due to a decrease in variation within the group (i.e., π). While additional models of allopatric genetic divergence have been proposed (e.g. "Recurrent selection" and "Geographic sweep before selective differentiation."), testing these two main conceptual models is a good first approach to better understand the evolutionary drivers of genomic divergence.

Studying how genomic divergence accumulates in marine organisms, leading to the formation of new species in widely connected environments, is crucial to the evolutionary debate regarding the importance of physical isolation during speciation (Faria et al. 2021). In these environments, dispersal (or at least its potential) is greater than in terrestrial habitats, where physical barriers abound (Kinlan et al. 2003; Burgess et al. 2016). Therefore, divergence due to selection in the presence of gene flow is expected to be predominant in the marine realm and is supported by a recent review modeling divergence across 66 marine species including 28 species of Chordata, 14 species of Mollusca, and 8 species of Cnidaria (Jode et al. 2023). This being said, the extent to which these results are validated by empirical investigations across the marine tree of life remains to be seen.

The Southeastern Pacific (SEP) provides an interesting scenario for studying organismal adaptation to heterogeneous environmental conditions. The SEP has two biogeographic transitions occurring between 30 and 34°S and 41 and 43°S (Camus 2001; Thiel et al. 2007). The transitional zone of 30-34°S marks a significant change in coastal upwelling regimes (Hormazábal et al. 2004). North of this zone, relatively weak but persistent coastal upwelling conditions are experienced (Navarrete et al. 2005). In contrast, to the south, the coastal upwelling is more intense but highly variable (seasonal variation). This shift in coastal upwelling regimes is reflected in species’ distributions. A widely congruent phylogeographic discontinuity has been detected in various marine organisms, including gastropods, tunicates, and seaweeds (Tellier et al. 2009; Sánchez et al. 2011; Brante et al. 2012; González et al. 2012; Lara et al. 2019; Quesada-Calderón et al. 2021). The southern transitional zone between 41°S and 43°S is located near the northern extent of the ice cap from the Last Glacial Maximum (LGM) and is associated with the oceanic barrier formed by the division of the West Wind Drift into the Humboldt Current to the north and the Cape Horn Current to the south (Strub et al. 2019). Characterized by high precipitation and strong salinity gradients (Dávila et al. 2002; Garreaud et al. 2013; Kilian et al. 2013), this transitional zone marks the northern limit for some cold-water seaweed species and appear as a genetic discontinuity in invertebrates with limited dispersal (Montecinos et al. 2012; Guillemin et al. 2016; Pardo‐Gandarillas et al. 2018; Fernández Iriarte et al. 2020; Montecinos et al. 2021). Here, we investigate the population differentiation of three limpet species of the monophyletic genus Scurria (S. scurra, S. araucana, and S. ceciliana) across two biogeographic breaks at 30°S - 34°S, and 41°S - 43°S. A recent study describing the population genetic structure of four of these species (Peluso et al. 2023) - S. scurra, S. araucana, S. ceciliana, and S. variabilis - has shown that S. scurra and S. araucana are comprised of three populations whose genetic breaks coincide with the division of the transition zone at 30 and 34°S, while S. ceciliana presents a genetic break separating populations located north and south of the transition zone at 41 and 43°S. Additionally, atypical loci were found in all three species, suggesting diversification processes driven by selection. These results suggest that biogeographic transition zones may have acted as barriers to gene flow, driving population diversification.

This study builds upon these previous results and, using whole-genome resequencing data, aims to describe the distribution of genetic divergence along the genome of these species, and tests if putative regions under selection follow one of two main models of selection, with gene flow or in isolation. We test this using the models proposed in the literature that compare relative and absolute divergence rates (Irwin et al. 2018). We compare two populations for each of the three species. Two species, S. scurra and S. araucana, were sampled across the same biogeographic break (30°S - 34°S), associated with environmental heterogeneity. The third pair crossing the southernmost biogeographic break (41°S - 43°S) originated most likely by vicariance, possibly driven by the last glaciation (Sanmartín 2012; Peluso et al. 2023). We identify regions of relative and absolute divergence for each pair and compare their distribution, sizes, and the genes that they harbour. By doing so, we give a novel perspective regarding the most likely scenarios under which these limpets are evolving and the genes that appear to be associated with this process.

Sample Preparation and Whole Genome Resequencing

Adult individuals of S. scurra, S. araucana, and S. ceciliana were collected between August 2018 and February 2020 at six sites ranging from 25°S to 53°S (Figure 2). All samples were collected under permit R. Ex. N 2036, 2019 from the Subsecretaría de Pesca y Acuicultura (SUBPESCA) of Chile. Two populations per species and ten individuals per population were sampled. For each species, the two populations were located north and south of the previously reported genetic breaks (i.e. located at 30 and 34°S, and 41 and 43°S in Chile, Figure 2; Peluso et al. 2023). For S. scurra and S. araucana, individuals from Temblador (29° 28' 16.2444" S) and Concepción (36° 31' 26.4864" S) were sampled. For S. ceciliana, individuals from Chiloé (41° 50' 48.12" S) and Magallanes (53° 37' 17.6988" S) were sampled.

Genomic DNA was extracted from foot and mantle tissues of limpets using the GeneJET Genomic DNA Purification Kit (Thermo Scientific). DNA concentration and quality were determined using the Qubit dsDNA BR Fluorometer (Thermo Scientific) and 1% (w/v) agarose gel electrophoresis. Fifty-one genomic DNA samples met library preparation quality requirements and where sent to the Beijing Genomics Institute (BGI) for resequencing using the DNBseq platform. Sequencing was performed at a depth of 10x with 150 bp paired-end reads. Details of sequenced samples can be found in Table S1 in the supplementary material.

Read Filtering Before Alignment

After sequencing, raw reads were filtered using the SOAPnuke software developed by BGI (Chen et al. 2018). Filtering steps included adapter removal (if the sequence read matched 25.0% or more of the adapter sequence, allowing a maximum of two base discrepancies), read length filtering (if the sequencing read length was less than 150 bp, the entire read was discarded), N removal (if the N content in the sequencing read represented 0.1% or more of the entire read, the entire read was discarded), and low-quality data filtering (if more than 40.0% of bases in a sequencing read had a Phred score below 10, the entire read was discarded). Clean reads were obtained with a Phred+33 quality readout.

Alignment and Variant Calling

The S. scurra genome (393 Mb), which consists of 9 putative chromosomes and approximately 1225 small contigs larger than 1 kb (Giles et al. in review), was used as a reference genome to map sequence reads. Filtered reads from each re-sequenced sample were mapped to the reference genome using the 'bwa mem' sequence alignment algorithm in Burrows-Wheeler Aligner (BWA version 0.7.17-r1188; Li and Durbin, 2009) using default settings (Li 2013). The generated SAM files were converted to BAM and sorted using SAMtools view and SAMtools sort, respectively, with SAMtools v.1.9 (Li et al. 2009). Subsequently, BCFtools v1.9 (Danecek et al. 2021) was used to compile BAM files into Variant Call Format (VCF) using the 'bcftool mpileup' function. Indels were excluded, and minimum mapping qualities for using an alignment of 30 (-min-MQ 30) and minimum base qualities of 20 (--min-BQ 20) were set. To determine the genomic positions with polymorphisms, the 'bcftools call' function was used. The generated VCF file was filtered using VCFtools (version 0.1.16; parameters: --max-alleles 2 --min-alleles 2 –minDP 4 --minQ 30 --mac 3 -- max-missing 0.8; Danecek et al. 2011).

Population genomic analyses

Genetic structure was explored using a principal component analysis (PCA) with the 'glPca()' command from the 'adegenet' package in R (R Core Team, 2022), employing a subset of 100,000 randomly selected SNPs. This data subset was used to reduce the analysis runtime. Subsequently, for each pairwise population comparison, VCFTools (Danecek et al. 2011) was used to perform a comprehensive genome scan with non-overlapping 20 kb windows. Afterward, for each window, we calculated relative divergence (F_ST, Weir and Cockerham) and absolute divergence (D_XY, Korunes and Samuk 2021).

Definition of outlier values

We classified highly divergent genomic regions into two categories according to the distribution of F_ST and D_XY estimates (for more details see Figure 1 that presents the patterns of expected F_ST and D_XY predicted under four evolutionary models of genomic divergence; Cruickshank and Hahn, 2014; Irwin et al. 2016; 2018). Genomic regions of high genomic divergence were defined as those with F_ST values above the 99^th percentile. For D_XY, the 95^th percentile of the D_XY distribution was used to define regions with low (below the 95_th percentile) and high (above the 95^th percentile) D_XY values. Based on these definitions, genomic regions with atypical values were identified and classified as follows:

Divergence with gene flow (DGF): 1% highest F_ST (P ≥ 99) and 5% highest D_XY (P ≥ 95)

Divergence in Allopatry (DA): 1% highest F_ST (P ≥ 99) and below 95% D_XY (P ≤ 0.95)

Identification of genes in genomic regions with outliers

All 20 kb adjacent windows considered atypical (see above) were grouped together to form a single genomic region. The size of each genomic region in kilobase pairs was calculated. Once identified, BED files were generated with one line assigned to each genomic region. The de-novo genome annotations of S. scurra were used with BEDtools (Quinlan and Hall, et al. 2010) to extract genes present in each genomic region. Subsequently, the Degenotate tool (https://github.com/harvardinformatics/degenotate) was used to identify synonymous and non-synonymous mutations in protein-coding genes. For this analysis, the VCF file from the samples was used as input, and the 6 individuals of S. ceciliana from the Magallanes population were defined as external groups. Both synonymous and non-synonymous mutations were counted for each gene, and potential mutation ratios were calculated (pN/pS). The relationship between non-synonymous mutations and the dataset of genes representative of the DGF and DA models was examined.

Statistics analysis

The Kruskal-Wallis test was used to determine if there were significant differences in the mapping rate between species. This non-parametric test is suitable for comparing medians among multiple groups when the data are not normal. Subsequently, a Dunn's post hoc test was performed to identify which pairs of species showed specific differences in the mapping rate. The Benjamini-Hochberg correction was applied to adjust p-values and control the Type I error. Furthermore, the size of genomic windows with outlier values (DGF and DA) and the number of genes in genomic windows with outlier values (DGF and DA) between species were compared using Kruskal-Wallis test of the "stats" package in R. The assumptions of normality and homoscedasticity for all tests were assessed using the shapiro.test() and leveneTest() functions from the "stats" and "car" packages in R, respectively. Independence assumption was evaluated graphically using model residuals.

An average of 4.82 GB of clean data was generated for each sample. The average mapping rate was 97.076% (0.40 SD) for S. scurra, 96.69% (0.26 SD) for S. araucana, and 96.105% (0.48 SD) for S. ceciliana with respect to the reference genome of S. scurra (Table S1). Dunn's post hoc test revealed specific comparisons with significant differences. Between S. araucana and S. ceciliana, a significant difference was observed with an adjusted p-value of 0.0179. Similarly, between S. araucana and S. scurra, the difference was also significant with an adjusted p-value of 0.0290. The comparison between S. ceciliana and S. scurra showed the most notable difference, with an adjusted p-value of 0.00000264 (Figure S1). A significant number of SNPs were analyzed for each species and population (Table S2). For S. araucana, 6,888,985 SNPs were analyzed in the Temblador population (n=10) and 3,801,687 SNPs in the Concepción population (n=5). For S. scurra, 5,756,211 SNPs were analyzed in the Temblador population (n=10) and 6,027,035 SNPs in the Concepción population (n=10). For S. ceciliana, 5,939,690 SNPs were analyzed in the Chiloé population (n=10) and 4,471,003 SNPs in the Magallanes population (n=6)

The principal component analysis (PCA) revealed that the first two components (PC1 and PC2) explained 51.89% of the total variability in the data (Figure 3). These components clearly discriminated three groups: S. scurra (Temblador and Concepción), S. araucana (Temblador and Concepción), and S. ceciliana (Chiloé and Magallanes). For S. araucana, there were five samples from Concepción that clustered entirely with S. ceciliana and were excluded from further analyses. Another five individuals, which were slightly separated from those in the Temblador population and slightly associated with the samples of S. ceciliana, were included in the data set (Figure 3).

Distribution of F_ST and D_XY Values

F_ST values (0-0.504) and D_XY values (0-0.061) varied substantially among the three studied Scurria species (Figure 4). The Kruskal-Wallis analysis revealed significant differences between the species regarding the F_STindex (Kruskal-Wallis χ² = 6697.9, df = 2, p < 0.05) and the D_XY index (Kruskal-Wallis χ² = 10.44, df = 2, p < 0.05). When comparing populations on either side of the first transition zone (30°S and 34°S), S. scurra had the highest average F_ST = 0.074, followed by S. araucana = 0.047; while when comparing populations of S. ceciliana on either side of the second transition zone (41°S and 43°S), the average F_ST was lower =0.009 (Figure S2-A). Chromosomes with F_ST outliers varied among species (Table 1 and Figure S2-B): ch02, ch03, ch04, ch05, ch06, ch07, ch08, and ch09 in S. scurra; ch01, ch02, ch03, ch04, ch05, ch06, ch07, ch08, and ch09 in S. araucana; and ch02, ch03, ch04, ch05, ch06, ch08, ch09 and in S. ceciliana. Average D_XYvalues were similar across species. Unlike the observed trend in average F_ST values, the species that displayed the lowest D_XY values was S. scurra (0.014), followed by S. araucana and S. ceciliana, which had the same average (0.017) (Figure S2-A). Although all chromosomes had D_XY outliers, some had significantly more outliers, and this varied among species (Table 1 and Figure S2-B): ch01, ch02, ch04, and ch08 in S. scurra; ch01, and ch02 in S. araucana; and ch01, and ch02 in S. ceciliana. In all species, chromosome 1 (ch01) showed the highest average D_XY.

Divergence patterns revealed through F_ST and D_XY outliers

In all species, the predominant model was DA, i.e. regions with high F_ST (P ≥ 0.99) and low D_XY (P ≤ 0.95) were the most abundant (Table 1 and Figure 4). The DGF model in S. scurra presented 45 regions with outlier values distributed across chromosomes ch02 (1), ch03 (1), ch04 (18), ch05 (3), ch06 (5), ch07 (13), ch08 (3), and ch09 (1); in S. araucana it was represented by 5 regions with outlier values distributed across chromosomes ch02 (3), ch04 (1), and ch08 (1); in S. ceciliana, the 72 regions with identified outlier values were distributed across chromosomes ch02 (8), ch05 (1), ch08 (61), and ch09 (2). For the DA model, in S. scurra 114 regions with outlier values were observed distributed across chromosomes ch02 (12), ch03 (9), ch04 (15), ch05 (23), ch06 (32), ch07 (12), ch08 (1), and ch09 (10); in S. araucana, 154 regions with outlier values were observed distributed across chromosomes ch01 (6), ch02 (39), ch03 (14), ch04 (27), ch05 (7), ch06 (4), ch07 (36), ch08 (7), and ch09 (14); while in S. ceciliana, the models of DGF (72) and DA (87) were almost equally represented. The 87 regions with outlier values associated with DA in S. ceciliana were distributed across chromosomes ch02 (1), ch03 (1), ch04 (1), ch05 (2), ch06 (1), ch08 (31), and ch09 (50).

Twelve genomic regions with outliers were shared between the species S. scurra and S. araucana. The first genomic region (DGF model) is located on chromosome 4 and contains a single gene (g22915; Table S3). The other 11 genomic regions are associated with the DA model and are located on chromosomes 4 and 7. They contain 3 genes (g24601, g24602, and g24603) and 8 genes (g16034, g16035, g16036, g16037, g16038, g16039, g1604, and g16072), respectively. All the genes shared between these species presented a high rate of non-synonymous mutations (pN/pS > 1; Table S3). The gene associated with the DGF model (g22915) is involved in lipid metabolic processes. Regarding the genes associated with the DA model, some exhibit unique functions while others share related functions. Among the genes with unique functions is g24603, which is involved in signaling and plasma membrane processes. The genes with related functions include: g16035 and g16038, involved in catalytic and hydrolase activity on proteins; g16036 and g16037, involved in metabolic processes of superoxide, copper ion binding, and oxidoreductase activity, as well as responding to oxidative stress and toxic substances.

S. scurra presented nine unique genes with a ratio pN/pS > 1 under the DFG model, which can be grouped into various functional categories (Table S3). In terms of membrane-related functions, several genes were identified: g1122 (pN/pS = 1.15), g15741 (pN/pS = 2.34), and g15747 (pN/pS = 2.29), related to the membrane; g19258 (pN/pS = 2.18), involved in transmembrane transport, symporter activity, and anion transport; and g9109 (pN/pS = 1.2), which participates in the G-protein coupled receptor signaling pathway, amine receptor activity coupled to G-protein, and the plasma membrane. On the other hand, gene g15576 (pN/pS = 1.48) is related to DNA binding, transcription initiation at RNA polymerase III promoter, the transcription factor TFIIIC complex, cellular anatomical entity, ncRNA transcription, and cellular macromolecule biosynthetic process. Regarding signaling and programmed cell death, gene g15578 (pN/pS = 2.01) is involved in these functions. Additionally, the gene g3038 (pN/pS = Inf) is associated with oxidoreductase activity, and the gene g3039 (pN/pS = 2.35), whose product, for example, a protein or a metabolite, is localized or has functions in the extracellular space.

Compared to S. scurra, S. araucana presented only two unique genes; however, these genes exhibit greater functional complexity and a wider diversity of specific cellular localizations (Table S3). The functions of these genes span various functional categories. In terms of ATP hydrolysis activity and transport, the genes g15029 (pN/pS = 1.01) and g15031 (pN/pS = 1.16) participate in ATP hydrolysis activity, transmembrane transport of xenobiotics, ABC transporter activity, long-chain fatty acid transport, and phospholipid transport. Additionally, they are involved in the transport and detoxification of xenobiotics across the plasma membrane and the blood-brain barrier. The blood-brain barrier protects the brain from potentially harmful substances present in the blood, but it also regulates the passage of nutrients and other elements necessary for brain function. In this way, they may participate in the response to stress and external compounds such as estradiol, oxidative stress, and lipopolysaccharide. Principio del formularioFinal del formulario

Regarding cellular and metabolic processes, these genes are associated with cobalamin metabolism, bile acid and bile salt transport, regulation of biological quality, nucleotide transport, and response to internal stimuli such as peptides and steroid hormones. Finally, in terms of cellular localization, these genes are associated with various cellular structures, such as the lateral and basolateral plasma membrane, the Golgi apparatus, cytoplasmic vesicles, the endosome, the perinuclear region of the cytoplasm, and fungal-type vacuoles.

For its part, S. ceciliana presented twelve unique genes, which are grouped into various functional categories (Table S3). In terms of biological regulation and stress response, gene g12028 (pN/pS = 1.11) is involved in biological regulation, transition metal ion homeostasis, superoxide metabolic processes, antioxidant activity, and cellular response to oxidative stress and toxic compounds. Regarding DNA binding and transcription regulation, gene g3630 (pN/pS = 1.21) participates in DNA binding and oxidoreductase activity, while gene g3855 (pN/pS = 1.18) is related to the nucleus, nucleosome binding, cellular component assembly, and chromatin organization. Additionally, gene g7048 (pN/pS = 1.03) is implicated in DNA repair and ubiquitin-mediated protein modification. Genes related to mitochondrial activity and energy processes include gene g3850 (pN/pS = 2), associated with the mitochondria, and gene g3892 (pN/pS = 1.76), involved in ATP hydrolysis activity, DNA replication, and double-strand break repair through homologous recombination. In terms of enzymatic and metabolic activity, genes g3859 (pN/pS = 1.54) and g3887 (pN/pS = 1.55) are related to transferase activity, phosphate metabolic process regulation, and endoplasmic reticulum stress response. Functionally, gene g3909 (pN/pS = 1.68) is associated with zinc ion binding and carbonate dehydratase activity in the extracellular region and plasma membrane. Additionally, some genes have structural and reproductive functions, such as gene g3833 (pN/pS = 1), which is a structural constituent of the ribosome, and gene g9865 (pN/pS = 1.15), which is involved in the reproduction of multicellular organisms. Finally, gene g9866 (pN/pS = 2.28) participates in biological quality regulation, retinal homeostasis, and chaperone-mediated protein complex assembly.

Size of genomic regions and number of genes

In all three species, 52.83% of the regions of genomic divergence were found in coding regions (Table 2 and Figure S3). Both S. scurra and S. araucana had the same number of genomic regions with coding regions (92). S. scurra had 23 coding regions associated with the DGF model and 69 associated with the DA model. On the other hand, S. araucana had 3 coding regions associated with the DGF model and 89 associated with the DA model. Finally, S. ceciliana had 68 coding regions, with 27 associated with the DGF model and 41 associated with the DA model. The Kruskal-Wallis test revealed significant differences between the studied species regarding the number of genes contained within the genomic regions (χ² = 21.576, gl = 2, p < 0.05). On the other hand, no significant differences were observed in the size of the genomic regions (χ² = 4.8446, df = 2, p = 0.08872). The distribution of the number of genes and the size of the genomic regions by model is illustrated in Figure 5. For the DGF model, the genomic regions in S. ceciliana contained the highest number of genes (136), followed by S. scurra (71) and S. araucana (10). The largest genomic regions also corresponded to S. ceciliana (45000 bp), followed by S. scurra (31515 bp) and S. araucana (26666 bp). For the DA model, S. araucana had the highest number of genes (225) compared to S. scurra (155) and S. ceciliana (148). Similar to the DGF model, the largest genomic regions were found in S. ceciliana (34800 bp), followed by S. araucana (29333 bp) and S. scurra (25618 bp).

In this study, we investigated genomic divergence in three closely related species of marine limpets inhabiting populations on both sides of two biogeographic breaks (30-34°S and 41-43°S). Highly divergent genomic regions ranging from 20 to 280 kb were detected in all three species, providing valuable insights into speciation in marine environments. Overall, the results indicate that allopatry has been the primary driver of diversification in this group of organisms. However, signals of divergence with gene flow were also detected in some regions of the genome. Genetic divergence appears to be mostly species-specific, with few shared regions among taxa. Interestingly, for S. scurra and S. araucana the distribution of F_ST and D_XY in the regions of high genomic divergence suggests that population divergence in these species occurred mostly in allopatry. Surprisingly, among population genomic divergence pattern in S. ceciliana was less homogeneous than expected under a purely neutral scenario and the patterns of genomic divergence suggest that selection played an important role in the differentiation among populations of this species in the south of Chile. The data suggest that selection acted in conjunction with both gene flow and the absence of gene flow, with both models contributing similarly to the process of population differentiation. The following sections discuss the levels of genomic heterogeneity found among the studied populations and the contribution of the two scenarios of evolution, divergence with gene flow and allopatric divergence, as drivers of the observed patterns.

The magnitude of genomic heterogeneity (variance) in F_ST and D_XY appears to be moderately correlated with the divergence times between genetic groups previously estimated in these species (Peluso et al. 2023). Species with longer divergence times between genetic groups (S. scurra: 84,000 years and S. araucana: 170,000 years) exhibit higher variation and higher F_ST values across their genome, than the species with a shorter divergence time between genetic groups (S. ceciliana: 56,000 years). Results obtained in S. ceciliana are consistent with recent population growth that would have reduced genetic differentiation between its populations (Peluso et al. 2023). However, Scurria scurra displays a more heterogeneous pattern and a higher average F_ST value than S. araucana, despite its shorter divergence time. We can only speculate about the causes of this discrepancy in heterogeneity, but one potential cause could be linked to the strong dependency of S. scurra on brown algae of the genus Lessonia and Durvillea. It is possible that population sizes of these macroalgae have modulated population sizes of S. scurra in the past. Smaller population sizes mean larger genetic drift which in turn can explain larger variance in genomic divergence. Another interesting observation is that the population structure of S. scurra appears to coincide with the genetic breaks observed in the cryptic species complex Lessonia berteroana and Lessonia spicata (Peluso et al. 2023). L. berteroana extends from Peru (17°S) to central Chile (30°S), while L. spicata ranges from central-southern Chile (29°S) to Puerto Montt (41°S), with their contact zone at 30°S. Additionally, south of 30°S, S. scurra is also associated with the recently described species Durvillaea incurvata (Fraser et al. 2020). These macroalgae differ in their morphology, chemical composition, and distribution (Nardelli et al. 2023). If local adaptation is linked to the association of these limpets to these different macroalgae species, and this is concomitant to adaptation to different environmental conditions on both sides of the discontinuity, these together could also result in larger than expected heterogeneity of genomic divergence compared to S. araucana (that lives on rocky substrates).

A general expected outcome in S. scurra and S. araucana was to identify adaptive genes, both shared between the species and unique to each, associated with environmental changes, such as sea surface temperature (SST) and salinity, characteristic of the transition zone between 30-34°S. Under this approach, we identified several regions of high genomic divergence shared between the sister species S. scurra and S. araucana. These regions harbored some interesting genes. The functions of gene g22915, shared between S. scurra and S. araucana, located on chromosome 4 include lipid metabolic process, extracellular region, and hydrolase activity. Previous studies on the intertidal limpet Cellana toreuma (Han et al. 2013) have shown that changes in sea surface temperature (SST) induce the expression of genes related to energy metabolism, including anabolic and catabolic processes to provide energy and cope with stress. In this regard, SST variation has been identified as a key factor in shaping disparities in thermal physiology among S. araucana populations at sites located between 27°S, 30°S, and 32°S off the Chilean coast (Rodríguez-Romero et al. 2022). These physiological differences correlate with the heterogeneity of the thermal marine landscape along the biogeographical break of 30-34°S. Thus, it is suggested that the conserved region on chromosome 4 could be involved in regulating adaptations to SST changes or in protecting essential protein functions. The high rate of non-synonymous mutations observed in this gene is a strong indication of positive selection at work. This suggests that beneficial mutations in this region have been favored by natural selection, enhancing the species' fitness to cope with changes in sea surface temperature (SST). The fact that both species show similar patterns of genetic divergence for this genomic region suggests local adaptation to environmental differences across this biogeographic barrier might be part of the mechanism behind population differentiation within both species of Scurra.

Most of the genomic divergence regions recovered here were unique to each species. This suggests that divergence between populations within species is mainly driven by species-specific evolutionary processes. The functions of the unique regions under the DGF and DA models were diverse in each species. In this study, we focused on the DGF model and those unique regions that exhibited non-synonymous mutations. Principio del formularioFinal del formulario

Among the genes identified in the unique regions of S. scurra, the gene g19258 stood out due to its high rate of non-synonymous substitutions (pN/pS = 2.18). Functional analysis of the gene g19258 revealed that it regulates processes related to transmembrane transport, symporter activity, and anion transport. Transmembrane transport is a relevant functional category for habitat transitions in mollusks (Aristide and Fernández, 2023), as it is associated with nutrient absorption, detoxification, osmoregulation, and shell formation, all of which are key physiological functions linked to adaptation to new environments. For example, in herbivores, transmembrane transporters are essential for protection against plant toxins and lipid metabolism (Sorensen and Dearing, 2006). The gene g15029 in S. araucana, similar to S. scurra, is closely linked to membrane transporters (pN/pS = 1.01). This gene plays a crucial role in the transfer of solutes across cell membranes in the digestive system, specifically in the transport of ABC-type bile acids. Additionally, g15029 regulates cellular signaling processes, including the response to oxidative stress. Oxidative stress in aquatic organisms can be induced by changes in salinity, as noted by Lesser (2006) and Bal et al. (2021). This type of stress occurs due to an altered but insufficient redox capacity that organisms experience when facing salinity fluctuations (Rivera-Ingraham and Lignot, 2017). It is important to highlight those similar tissues in different animals, as well as different tissues within the same animal, can respond differently and adopt various adaptation mechanisms to environmental salinity changes (Paital and Chainy, 2010). Consequently, the demand and supply of energy, along with redox regulation responses, vary among animals (Paital and Chainy, 2012; 2013; 2014). Similar redox regulation responses have been observed in mollusks (Rivera-Ingraham et al. 2016), annelids (Freitas et al. 2016), arthropods (Delgado-Gaytán et al. 2020), echinoderms (Castellano et al. 2016), and fish (Paul and Rashid, 2017). This suggests that the gene g15029 may confer oxidative stress responses associated with abiotic factors, such as salinity, on both sides of the biogeographic barrier of 30-34°S.

Principio del formulario

Final del formulario

This study had anticipated that the species S. ceciliana, whose distribution extends across the area covered by the Last Glacial Maximum experienced an effect of glacial vicariance. This biogeographic process involves the fragmentation of populations due to the emergence of geographical barriers such as glaciers, which restricts gene flow between glacial refugia. In this context, it was expected to observe significant genomic divergence among the isolated populations, resulting from the accumulation of genetic differences due to drift over time. Glacial vicariance has been identified as a driver of speciation in the diversification of snails of the genus Laevilitorina in the Southern Ocean during the Quaternary geological processes (González-Wevar et al. 2023). It is likely that small snail populations, with low dispersal capacity, became isolated in refugia during these glacial events. This would have generated substantial genetic drift, promoting the diversification of the genus Laevilitorina through the accumulation of mutations and adaptations to new local environments, like what has been observed in other marine organisms in the region, such as the bivalve Aequiyoldia eightsii (González-Wevar et al. 2019), the gastropods Nacella (González-Wevar et al. 2011; 2017), and Tritoniella (Schächinger et al. 2022). The genetic divergence observed in S. ceciliana suggests an incipient differentiation with a very low F_ST, indicating that the populations have not yet fully differentiated. The presence of regions with high F_ST suggests differential selection in different parts of the genome. Although allopatry is the predominant model for this divergence, there is also significant evidence of gene flow between populations. This suggests that despite initial geographic isolation, populations have come into secondary contact, allowing gene flow and hindering complete differentiation. It is worth noting that, among the three species studied, S. ceciliana exhibited the highest number of genomic regions indicative of divergence with gene flow (72). Twelve unique regions with a pN/pS>1 were identified in S. ceciliana. A gene that exhibited high divergence in S. ceciliana was g3850 (pN/pS=2), which regulates the mitochondrion, a key cellular component. This gene, which plays an essential role in mitochondrial function, showed notable differentiation between the Chiloé and Magallanes populations of S. ceciliana. The observed high divergence could be related to variations in the efficiency of cellular respiration and energy production, which are fundamental aspects for the survival and biological performance of individuals in their respective habitats. The mitochondria of intertidal mollusks, such as clams, oysters, and mussels, demonstrate remarkable tolerance to changes in the osmolarity and ionic strength of their environment (Sokolova, 2018; Sokolov and Sokolova, 2018). These findings suggest that S. ceciliana has developed specific genomic adaptations to cope with environmental challenges. The presence of genes related to mitochondrial function in regions of divergence with gene flow indicates that these adaptations may have been positively selected in response to stressful environmental conditions.

It is crucial to note that speciation is a complex process that cannot be fully captured by a simple model. Future studies focusing on specific recombination estimates for these species would allow the identification of highly divergent genomic regions associated with low-recombination regions (Booker et al. 2020) and, thus, with mechanisms different from those evaluated in this study. Nevertheless, our findings provide a comprehensive insight into the mechanisms of divergence in non-model marine organisms and highlights the importance of biogeographical barriers and natural selection in species evolution. Our research presents molecular evidence confirming that divergence begins with small areas of differentiation in the early stages, spreading across the genome as the process progresses. Additionally, we introduce a novel method for identifying signatures of selection in biogeographical diversification scenarios. We identify signatures of selection in the same genomic regions in the two species living in sympatry (S. scurra and S. araucana), suggesting that the same selective pressure has acted on the same regions of the genome in both species, driving population differentiation. Overall, unexpected findings in the species S. scurra (high heterogeneity) and S. ceciliana (incipient divergence) highlight the complexity of divergence in non-model organisms and the need to continue exploring the underlying mechanisms. Future research on these species could shed light on the evolutionary dynamics in intertidal environments and contribute to a broader understanding of speciation in marine environments.

Acknowledgements

We acknowledge financial support from the Agencia Nacional de Investigación y Desarrollo ANID (grant: FONDECYT 1190710, and the Millennium Nucleus NUTME NCN2023_004, to PS-A). We are grateful to all Saenz, Navarrete and Chango lab members who assisted in the field and laboratory.

Author Contributions

Conceptualisation: PC, PS-A, MLG, Designed research: PC, PS-A, Performed research: PC, MLG, ECG, GN, VAS, PS-A, contributed new reagents or analytical tools: PC, MLG, ECG, GN, VAS, PS-A. Analysed the data: PC, ECG. Writing original draft: PC, Writing reviewing and editing manuscript: PC, MLG, ECG, GN, VAS, PS-A

Competing interests

The authors declare no competing interests.

Andrew RL, Rieseberg LH (2013) Divergence is focused on few genomic regions early in speciation: incipient speciation of sunflower ecotypes. Evolution 67(9):2468-82. https://doi.org/10.1111/evo.12106
Aristide L, Fernández R (2023) Genomic Insights into Mollusk Terrestrialization: Parallel and Convergent Gene Family Expansions as Key Facilitators in Out-of-the-Sea Transitions. Genome Biol Evol 15(10):evad176. https://doi.org/10.1093/gbe/evad176
Bal A, Panda F, Pati SG, Das K, Agrawal PK, Paital B (2021) Modulation of physiological oxidative stress and antioxidant status by abiotic factors especially salinity in aquatic organisms. Comp Biochem Physiol Part C Toxicol Pharmacol 241:108971. https://doi.org/10.1016/j.cbpc.2020.108971
Booker TR, Yeaman S, Whitlock MC (2020) Variation in recombination rate affects detection of outliers in genome scans under neutrality. Mol Ecol 29(22):4274-9. https://doi.org/10.1111/mec.15501
Brante A, Fernandez M, Viard F (2012) Phylogeography and biogeography concordance in the marine gastropod Crepipatella dilatata (Calyptraeidae) along the southeastern Pacific coast. J Hered 103(5):630-7. https://doi.org/10.1093/jhered/ess030
Burgess SC, Baskett ML, Grosberg RK, Morgan SG, Strathmann RR (2016) When is dispersal for dispersal? Unifying marine and terrestrial perspectives. Biol Rev 91(3):867-82. https://doi.org/10.1111/brv.12198
Burri R, Nater A, Kawakami T, Mugal CF, Olason PI, Smeds L, et al. (2015) Linked selection and recombination rate variation drive the evolution of the genomic landscape of differentiation across the speciation continuum of Ficedula flycatchers. Genome Res 25(11):1656-65. https://doi.org/10.1101/gr.196485.115.
Camus PA. Biogeografía marina de Chile continental (2001) Rev Chil Hist Nat 74(3):587-617. http://dx.doi.org/10.4067/S0716-078X2001000300008
Castellano GC, Santos IA, Freire CA (2018) Maintenance of ionic gradients and tissue hydration in the intertidal sea cucumber Holothuria grisea under hypo-and hyper-salinity challenges. J Mar Biol Assoc U K 98(2):325-32. https://doi.org/10.1017/S0025315416001314
Chen Y, Chen Y, Shi C, Huang Z, Zhang Y, Li S, et al. (2018) SOAPnuke: a MapReduce acceleration-supported software for integrated quality control and preprocessing of high-throughput sequencing data. Gigascience 7(1):gix120. https://doi.org/10.1093/gigascience/gix120
Clarkson CS, Weetman D, Essandoh J, Yawson AE, Maslen G, Manske M, et al. (2014) Adaptive introgression between Anopheles sibling species eliminates a major genomic island but not reproductive isolation. Nat Commun 5(1):4248. https://doi.org/10.1186/1756-3305-7-345
Cruickshank TE, Hahn MW (2014) Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. Mol Ecol 23(13):3133-57. https://doi.org/10.1111/mec.12796
Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. (2011) The variant call format and VCFtools. Bioinformatics 27(15):2156-8.
Danecek P, Bonfield JK, Liddle J, Marshall J, Ohan V, Pollard MO, et al. (2021) Twelve years of SAMtools and BCFtools. Gigascience 10: giab008. https://doi.org/10.1093/gigascience/giab008
Dávila PM, Figueroa D, Müller E (2002) Freshwater input into the coastal ocean and its relation with the salinity distribution off austral Chile (35–55 S). Cont Shelf Res 22(3):521-34. https://doi.org/10.1016/S0278-4343(01)00072-3
Delgado-Gaytán MF, Gómez-Jiménez S, Gámez-Alejo LA, Rosas-Rodríguez JA, Figueroa-Soto CG, Valenzuela-Soto EM (2020) Effect of salinity on the synthesis and concentration of glycine betaine in osmoregulatory tissues from juvenile shrimps Litopenaeus vannamei. Comp Biochem Physiol A Mol Integr Physiol 240:110628. https://doi.org/10.1016/j.cbpa.2019.110628
Delmore KE, Lugo Ramos JS, Van Doren BM, Lundberg M, Bensch S, Irwin DE, et al. (2018) Comparative analysis examining patterns of genomic differentiation across multiple episodes of population divergence in birds. Evol Lett 2(2):76-87. https://doi.org/10.1002/evl3.46
De Jode A, Le Moan A, Johannesson K, Faria R, Stankowski S, Westram AM, et al. (2023) Ten years of demographic modelling of divergence and speciation in the sea. Evol Appl 16(2):542-59. https://doi.org/10.1111/eva.13428
Faria R, Johannesson K, Stankowski S (2021) Speciation in marine environments: Diving under the surface. J Evol Biol 34(1):4-15. https://doi.org/10.1111/jeb.13756
Feder JL, Egan SP, Nosil P (2012) The genomics of speciation-with-gene-flow. Trends Genet 28(7):342-50. https://doi.org/10.1016/j.tig.2012.03.009
Fernandez Iriarte PJ, González-Wevar CA, Segovia NI, Rosenfeld S, Hüne M, Fainburg L, et al. (2020) Quaternary ice sheets and sea level regression drove divergence in a marine gastropod along Eastern and Western coasts of South America. Sci Rep 10(1):844. https://doi.org/10.1038/s41598-020-57543-4
Fraser CI, Velásquez M, Nelson WA, Macaya EC, Hay CH (2020) The Biogeographic Importance of Buoyancy in Macroalgae: A Case Study of the Southern Bull-Kelp Genus Durvillaea (Phaeophyceae), Including Descriptions of Two New Species¹. J Phycol 56(1):23-36. https://doi.org/10.1111/jpy.12939
Freitas R, Pires A, Velez C, Almeida Â, Wrona FJ, Soares AM, et al. (2015) The effects of salinity changes on the Polychaete Diopatra neapolitana: impacts on regenerative capacity and biochemical markers. Aquat Toxicol 163:167-76. https://doi.org/10.1016/j.aquatox.2015.04.006
Garreaud R, Lopez P, Minvielle M, Rojas M (2013) Large-scale control on the Patagonian climate. J Clim 26(1):215-30. https://doi.org/10.1175/JCLI-D-12-00001.1
Gavrilets S, Vose A (2005) Dynamic patterns of adaptive radiation. Proc Natl Acad Sci 102(50):18040-5. https://doi.org/10.1073/pnas.0506330102
Giles EC, González VL, Leiva C, Carimán P, Suescun V, Lemer S, Guillemin ML, Ortiz-Barrientos D, Saenz-Agudelo P. (in review) Comparative genomics points to ecological drivers of genomic divergence among intertidal limpets.
Gompert Z, Lucas LK, Nice CC, Buerkle CA (2013) Genome divergence and the genetic architecture of barriers to gene flow between Lycaeides idas and L. melissa. Evolution 67(9):2498-514. https://doi.org/10.1111/evo.12021
González A, Beltrán J, Hiriart-Bertrand L, Flores V, de Reviers B, Correa JA, et al. (2012) Identification of cryptic species in the Lessonia nigrescens complex (Phaeophyceae, Laminariales). J Phycol 48(5):1153-65. https://doi.org/10.1111/j.1529-8817.2012.01200.x
Gonzalez-Wevar CA, de Aranzamendi MC, Segovia NI, Rosenfeld S, Maturana CS, Molina CR, et al. (2023) Genetic footprints of Quaternary glacial cycles over the patterns of population diversity and structure in three Nacella (Patellogastropoda: Nacellidae) species across the Magellan province in southern South America. Front Mar Sci 10:1154755. https://doi.org/10.3389/fmars.2023.1154755
González-Wevar CA, Gérard K, Rosenfeld S, Saucède T, Naretto J, Díaz A, et al. (2019) Cryptic speciation in Southern Ocean Aequiyoldia eightsii (Jay, 1839): mio-pliocene trans-drake passage separation and diversification. Prog Oceanogr 174:44-54. https://doi.org/10.1016/j.pocean.2018.09.004
González‐Wevar CA, Hüne M, Segovia NI, Nakano T, Spencer HG, Chown SL, et al. (2017) Following the Antarctic Circumpolar Current: patterns and processes in the biogeography of the limpet Nacella (Mollusca: Patellogastropoda) across the Southern Ocean. J Biogeogr 44(4):861-74.
González‐Wevar CA, Nakano T, Cañete JI, Poulin E (2011) Concerted genetic, morphological and ecological diversification in Nacella limpets in the Magellanic Province. Mol Ecol 20(9):1936-51. https://doi.org/10.1111/j.1365-294X.2011.05065.x
Guillemin ML, Contreras-Porcia L, Ramírez ME, Macaya EC, Contador CB, Woods H, et al. (2016) The bladed Bangiales (Rhodophyta) of the South Eastern Pacific: molecular species delimitation reveals extensive diversity. Mol Phylogenet Evol 94:814-26. https://doi.org/10.1016/j.ympev.2015.09.027
Han G dong, Zhang S, Marshall DJ, Ke C huan, Dong Y wei (2013) Metabolic energy sensors (AMPK and SIRT1), protein carbonylation and cardiac failure as biomarkers of thermal stress in an intertidal limpet: linking energetic allocation with environmental temperature during aerial emersion. J Exp Biol 216(17):3273-82. https://doi.org/10.1242/jeb.084269
Hormazabal S, Shaffer G, Leth O (2004) Coastal transition zone off Chile. J Geophys Res Oceans 109(C1). https://doi.org/10.1029/2003JC001956
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. Mol Ecol 25(18):4488-507. https://doi.org/10.1111/mec.13792
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. Mol Ecol 27(23):4839-55. https://doi.org/10.1111/mec.14858
Kinlan BP, Gaines SD (2003) Propagule dispersal in marine and terrestrial environments: a community perspective. Ecology 84(8):2007-20. https://doi.org/10.1890/01-0622
Korunes KL, Samuk K (2021) pixy: Unbiased estimation of nucleotide diversity and divergence in the presence of missing data. Mol Ecol Resour 21(4):1359-68. https://doi.org/10.1111/1755-0998.13326
Lara C, Saldías GS, Cazelles B, Rivadeneira MM, Haye PA, Broitman BR (2019) Coastal biophysical processes and the biogeography of rocky intertidal species along the south-eastern Pacific. J Biogeogr 46(2):420-31. https://doi.org/10.1111/jbi.13492
Lawrence M, Huber W, Pagès H, Aboyoun P, Carlson M, Gentleman R, et al. (2013) Software for computing and annotating genomic ranges. PLoS Comput Biol 9(8):e1003118. https://doi.org/10.1371/journal.pcbi.1003118
Lesser MP (2006) Oxidative stress in marine environments: Biochemistry and Physiological Ecology. Annu Rev Physiol 68(1):253-78. https://doi.org/10.1146/annurev.physiol.68.040104.110001
Li H (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2. https://doi.org/10.48550/arXiv.1303.3997
Li H, Durbin R (2009) Fast and accurate short read alignment with Burrows–Wheeler transform. bioinformatics. 2009;25(14):1754-60. https://doi.org/10.1093/bioinformatics/btp324
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. (2009) The sequence alignment/map format and SAMtools. Bioinformatics 25(16):2078-9. https://doi.org/10.1093/bioinformatics/btp352
Montecinos A, Broitman BR, Faugeron S, Haye PA, Tellier F, Guillemin ML (2012) Species replacement along a linear coastal habitat: phylogeography and speciation in the red alga Mazzaella laminarioides along the south east pacific. BMC Evol Biol 12(1):1-17. https://doi.org/10.1186/1471-2148-12-97
Montecinos AE, Huanel OR, Ramírez ME, Guillemin ML (2021) Molecular data reveal the presence of three Plocamium Lamouroux species with complex patterns of distribution in Southern Chile. Cryptogam Algol 42(1):1-19. https://doi.org/10.5252/cryptogamie-algologie2021v42a1
Nardelli AE, Visch W, Wright JT, Hurd CL (2023) Concise review of genus Lessonia Bory. J Appl Phycol 35(4):1485-98. https://doi.org/10.5252/cryptogamie-algologie2021v42a1
Navarrete SA, Wieters EA, Broitman BR, Castilla JC (2005) Scales of benthic–pelagic coupling and the intensity of species interactions: from recruitment limitation to top-down control. Proc Natl Acad Sci 102(50):18046-51. https://doi.org/10.1073/pnas.0509119102
Nosil P, Funk DJ, Ortiz-Barrientos D (2009) Divergent selection and heterogeneous genomic divergence. Mol Ecol 18(3):375-402. https://doi.org/10.1111/j.1365-294X.2008.03946.x
Ortiz-Barrientos D, Engelstädter J, Rieseberg LH (2016) Recombination rate evolution and the origin of species. Trends Ecol Evol 31(3):226-36. https://doi.org/10.1016/j.tree.2015.12.01
Espinosa EP, Koller A, Allam B (2016) Proteomic characterization of mucosal secretions in the eastern oyster, Crassostrea virginica. J Proteomics 132:63-76. https://doi.org/10.1016/j.jprot.2015.11.018
Park K, Kwak IS (2020) Cadmium-induced developmental alteration and upregulation of serine-type endopeptidase transcripts in wild freshwater populations of Chironomus plumosus. Ecotoxicol Environ Saf 192:110240. https://doi.org/10.1016/j.ecoenv.2020.110240
Paul BK, Rashid H (2017) Salinity intrusion and impacts. Clim Hazards Coast Bangladesh. 2017;153-82. https://doi.org/10.1016/B978-0-12-805276-1.00005-3
Paital B, Chainy GBN (2010) Antioxidant defenses and oxidative stress parameters in tissues of mud crab (Scylla serrata) with reference to changing salinity. Comp Biochem Physiol Part C Toxicol Pharmacol 151(1):142-51. https://doi.org/10.1016/j.cbpc.2009.09.007
Paital B, Chainy GBN (2012) Effects of salinity on O2 consumption, ROS generation and oxidative stress status of gill mitochondria of the mud crab Scylla serrata. Comp Biochem Physiol Part C Toxicol Pharmacol 155(2):228-37. https://doi.org/10.1016/j.cbpc.2011.08.009
Paital B, Chainy GBN (2013) Modulation of expression of SOD isoenzymes in mud crab (Scylla serrata): effects of inhibitors, salinity and season. J Enzyme Inhib Med Chem 28(1):195-204. https://doi.org/10.3109/14756366.2011.645239
Paital B, Chainy GBN (2014) Effects of temperature on complexes I and II mediated respiration, ROS generation and oxidative stress status in isolated gill mitochondria of the mud crab Scylla serrata. J Therm Biol 41:104-11. https://doi.org/10.1016/j.jtherbio.2014.02.013
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. Mol Ecol 32(14):3812-25. https://doi.org/10.1111/mec.16978
Poelstra JW, Vijay N, Bossu CM, Lantz H, Ryll B, Müller I, et al. (2014) The genomic landscape underlying phenotypic integrity in the face of gene flow in crows. Science 344(6190):1410-4. https://doi.org/10.1126/science.1253226
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 of a commercially important tunicate. J Biogeogr 48(10):2604-15. 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-2. https://doi.org/10.1093/bioinformatics/btq033
Rivera-Ingraham GA, Nommick A, Blondeau-Bidet E, Ladurner P, Lignot JH (2016) Salinity stress from the perspective of the energy-redox axis: lessons from a marine intertidal flatworm. Redox Biol 10:53-64. https://doi.org/10.1016/j.redox.2016.09.012
Rivera-Ingraham GA, Lignot JH (2017) Osmoregulation, bioenergetics and oxidative stress in coastal marine invertebrates: raising the questions for future research. J Exp Biol 220(10):1749-60. https://doi.org/10.1242/jeb.135624
Rodríguez-Romero A, Gaitán-Espitía JD, Opitz T, Lardies MA (2022) Heterogeneous environmental seascape across a biogeographic break influences the thermal physiology and tolerances to ocean acidification in an ecosystem engineer. Divers Distrib 28(8):1542-53. https://doi.org/10.1111/ddi.13478
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). Mar Ecol Prog Ser 434:121-31. https://doi.org/10.3354/meps09184
Schächinger PM, Schrödl M, Wilson NG, Moles J (2022) Crossing the polar front—Antarctic species discovery in the nudibranch genus Tritoniella (Gastropoda). Org Divers Evol 22(2):431-56. https://doi.org/10.1007/s13127-022-00541-3
Strub PT, James C, Montecino V, Rutllant JA, Blanco JL (2019) Ocean circulation along the southern Chile transition region (38–46 S): Mean, seasonal and interannual variability, with a focus on 2014–2016. Prog Oceanogr 172:159-98. https://doi.org/10.1016/j.pocean.2019.01.004
Sokolov EP, Sokolova IM (2019) Compatible osmolytes modulate mitochondrial function in a marine osmoconformer Crassostrea gigas (Thunberg, 1793). Mitochondrion 45:29-37. https://doi.org/10.1016/j.mito.2018.02.002
Sokolova I (2018) Mitochondrial adaptations to variable environments and their role in animals’ stress tolerance. Integr Comp Biol 58(3):519-31. https://doi.org/10.1093/icb/icy017
Sorensen JS, Dearing MD (2006) Efflux Transporters as a Novel Herbivore Countermechanism to Plant Chemical Defenses. J Chem Ecol 32(6):1181-96. https://doi.org/10.1007/s10886-006-9079-y
Sun Y, Sun J, Yang Y, Lan Y, Ip JCH, Wong WC, et al. (2021) Genomic signatures supporting the symbiosis and formation of chitinous tube in the deep-sea tubeworm Paraescarpia echinospica. Mol Biol Evol 38(10):4116-34. https://doi.org/10.1093/molbev/msab203
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? Mol Phylogenet Evol 53(3):679-93. https://doi.org/10.1016/j.ympev.2009.07.030
Thiel M, Castilla JC, Fernández Bergia ME, Navarrete S (2007) The Humboldt current system of northern and central Chile.
Van Belleghem SM, Cole JM, Montejo-Kovacevich G, Bacquet CN, McMillan WO, Papa R, et al. (2021) Selection and isolation define a heterogeneous divergence landscape between hybridizing Heliconius butterflies. Evolution 75(9):2251-68. https://doi.org/10.1111/evo.14272
Venier P, Varotto L, Rosani U, Millino C, Celegato B, Bernante F, et al. (2011) Insights into the innate immunity of the Mediterranean mussel Mytilus galloprovincialis. BMC Genomics 12(1):69. https://doi.org/10.1186/1471-2164-12-69
Wang S, Rohwer S, de Zwaan DR, Toews DP, Lovette IJ, Mackenzie J, et al. (2020) Selection on a small genomic region underpins differentiation in multiple color traits between two warbler species. Evol Lett 4(6):502-15. https://doi.org/10.1002/evl3.198
Wang S, Coop G (2022) A complex evolutionary history of genetic barriers to gene flow in hybridizing warblers. bioRxiv 2022-11. https://doi.org/10.1101/2022.11.14.516535
Weir BS, Cockerham CC (1984) Estimating F-statistics for the analysis of population structure. Evolution 1358-70. https://doi.org/10.2307/2408641
Wu CI (2001) The genic view of the process of speciation. J Evol Biol 14(6):851-65. https://doi.org/10.1046/j.1420-9101.2001.00335.x
Xue Q, Itoh N, Schey KL, Cooper RK, La Peyre JF (2009) Evidence indicating the existence of a novel family of serine protease inhibitors that may be involved in marine invertebrate immunity. Fish Shellfish Immunol 27(2):250-9. https://doi.org/10.1016/j.fsi.2009.05.006
Zhou W, Zhang N, Huang K, Lin H, Tu J, Zheng C, et al. (2024) Divergent selection in low recombination regions shapes the genomic islands in two incipient shorebird species. Mol Biol Evol 41(2):msae006. https://doi.org/10.1093/molbev/msae006

Tables 1 to 2 are available in the Supplementary Files section

There is no duality of interest

Table1.xlsx
Table 1. Number of windows with outlier values calculated across the 9 chromosomes of the three species, using the reference genome of S. scurra. The codes TE, CN, CH, and MA represent the Temblador, Concepción, Chiloé, and Magallanes populations, respectively. Outliers were classified as follows: above the 99th percentile of the mean F_ST values, above the 95^th percentile of the mean D_XY values, and below the 95^th percentile of the mean D_XY values.
Table2.xlsx
Table 2. Summary of the number of genomic regions identified from scans of 20 kb windows per model (DGF=Divergence with Gene Flow, DA=Divergence in Allopatry) and species. Genomic regions containing genes are identified, including the number of genes present in these regions, as well as the number of transcripts and the number of transcripts with associated ontology.
TableS1.xlsx
Table S1
TableS2.xlsx
Table S2
TableS3.xlsx
Table S3
FigureS1rev.pdf
Figure S1
FiguraS2rev.pdf
Figure S2
FigureS3rev.pdf
Figure S3

Download PDF

Review #2 received at journal
28 Aug, 2024
Review #1 received at journal
18 Aug, 2024
Reviewer #3 agreed at journal
13 Aug, 2024
Reviewer #2 agreed at journal
13 Aug, 2024
Reviewer #1 agreed at journal
04 Aug, 2024
Reviewers invited by journal
04 Aug, 2024
First submitted to journal
03 Jul, 2024
Editor assigned by journal
03 Jul, 2024

You are reading this latest preprint version

Exploring evolutionary mechanisms of genomic divergence in Scurria limpets

Status:

Version 1

Abstract

Figures

Introduction

Methods

Results

Discussion

Declarations

References

Tables

Additional Declarations

Supplementary Files

Status:

Version 1